Index: /branches/eam_branches/20090715/psphot/src/Makefile.am =================================================================== --- /branches/eam_branches/20090715/psphot/src/Makefile.am (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/Makefile.am (revision 25178) @@ -17,13 +17,14 @@ # Force recompilation of psphotVersion.c, since it gets the version information psphotVersion.c: psphotVersionDefinitions.h -psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE - -$(RM) psphotVersionDefinitions.h - $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h -FORCE: ; +# psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE +# -$(RM) psphotVersionDefinitions.h +# $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h +# FORCE: ; libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy +# bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy +bin_PROGRAMS = psphotPetrosianStudy psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) Index: /branches/eam_branches/20090715/psphot/src/pmPetrosian.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/pmPetrosian.c (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/pmPetrosian.c (revision 25178) @@ -24,5 +24,6 @@ psFree(petrosian->fluxElliptical); - psFree(petrosian->binnedFlux); + psFree(petrosian->binSB); + psFree(petrosian->binSBstdev); psFree(petrosian->radialBins); psFree(petrosian->area); @@ -44,5 +45,6 @@ petrosian->radialBins = NULL; petrosian->area = NULL; - petrosian->binnedFlux = NULL; + petrosian->binSB = NULL; + petrosian->binSBstdev = NULL; petrosian->petrosianRadius = NAN; @@ -64,5 +66,6 @@ psFree(petrosian->fluxElliptical); - psFree(petrosian->binnedFlux); + psFree(petrosian->binSB); + psFree(petrosian->binSBstdev); psFree(petrosian->radialBins); psFree(petrosian->area); @@ -78,5 +81,6 @@ petrosian->radialBins = NULL; petrosian->area = NULL; - petrosian->binnedFlux = NULL; + petrosian->binSB = NULL; + petrosian->binSBstdev = NULL; return true; Index: /branches/eam_branches/20090715/psphot/src/pmPetrosian.h =================================================================== --- /branches/eam_branches/20090715/psphot/src/pmPetrosian.h (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/pmPetrosian.h (revision 25178) @@ -20,5 +20,6 @@ psVector *fluxElliptical; // flux for the above radial coordinates - psVector *binnedFlux; // mean surface brightness within radial bins + psVector *binSB; // mean surface brightness within radial bins + psVector *binSBstdev; // scatter of mean surface brightness within radial bins psVector *radialBins; // radii corresponding to above binnedBlux psVector *area; // differential area of the non-overlapping radial bins @@ -52,6 +53,7 @@ bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, psVector *refRadius, psVector *meanSB, - psVector *petRatio, psVector *fluxSum, - float petRadius, float petFlux); + psVector *petRatio, psVector *petRatioErr, psVector *fluxSum, + float petRadius, float ratioForRadius, + float petFlux, float radiusForFlux); bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec); Index: /branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c (revision 25178) @@ -73,13 +73,19 @@ psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, psphotEllipticalContourFunc); + /// XXX rationalize? if epsilon > 1, flip major and minor axes (rotate by 90 degrees) + if (params->data.F32[PAR_EPSILON] < 1.0) { + petrosian->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; + petrosian->axes.minor = params->data.F32[PAR_RMIN]; + petrosian->axes.theta = params->data.F32[PAR_PHI]; + } else { + petrosian->axes.major = params->data.F32[PAR_RMIN]; + petrosian->axes.minor = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; + petrosian->axes.theta = params->data.F32[PAR_PHI] + 0.5*M_PI; + } + fprintf (stderr, "# fitted values:\n"); fprintf (stderr, "Po: %f\n", params->data.F32[PAR_PHI]*PS_DEG_RAD); fprintf (stderr, "Ep: %f\n", params->data.F32[PAR_EPSILON]); fprintf (stderr, "Rm: %f\n", params->data.F32[PAR_RMIN]); - - /// XXX rationalize? if epsilon > 1, flip major and minor axes (rotate by 90 degrees) - petrosian->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; - petrosian->axes.minor = params->data.F32[PAR_RMIN]; - petrosian->axes.theta = params->data.F32[PAR_PHI]; // show the results Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c (revision 25178) @@ -61,5 +61,5 @@ psphotPetrosianFreeVectors(petrosian); - fprintf (stderr, "petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f\n", + fprintf (stdout, "\n petrosian radius: %f\n flux: %f\n axis ratio: %f\n angle: %f\n", petrosian->petrosianRadius, petrosian->petrosianFlux, petrosian->axes.minor/petrosian->axes.major, PS_DEG_RAD*petrosian->axes.theta); Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (revision 25178) @@ -12,7 +12,4 @@ // finely spaced than r_{i+1} = r_i * \beta / \alpha. for the integration, we need to // track the non-overlapping radius values. - -# define PETROSIAN_ALPHA 0.8 -# define PETROSIAN_BETA 1.25 bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax) { @@ -33,10 +30,12 @@ psVector *radBet = psVectorAllocEmpty(nMax, PS_TYPE_F32); - psVector *fluxBin = psVectorAllocEmpty(nMax, PS_TYPE_F32); - psVector *radBin = psVectorAllocEmpty(nMax, PS_TYPE_F32); - psVector *area = psVectorAllocEmpty(nMax, PS_TYPE_F32); + psVector *binSB = psVectorAllocEmpty(nMax, PS_TYPE_F32); // surface brightness of radial bin + psVector *binSBstdev = psVectorAllocEmpty(nMax, PS_TYPE_F32); // surface brightness of radial bin + psVector *binRad = psVectorAllocEmpty(nMax, PS_TYPE_F32); // mean radius of radial bin + psVector *binArea = psVectorAllocEmpty(nMax, PS_TYPE_F32); // area of radial bin (contiguous, non-overlapping) - psVectorInit (fluxBin, 0.0); - psVectorInit (radBin, 0.0); + psVectorInit (binSB, 0.0); + psVectorInit (binSBstdev, 0.0); + psVectorInit (binRad, 0.0); // generate radial bin bounds @@ -56,11 +55,25 @@ radBet->data.F32[2] = 2.0; - int nPts = 3; - for (int i = 3; i < radiusMax; i++) { - radMin->data.F32[nPts] = (i - 1); - radMax->data.F32[nPts] = i; - nPts++; +# define PETROSIAN_ALPHA 0.8 +# define PETROSIAN_BETA 1.25 +# define POWER_LAW_SPACING true + + // power-law spacing with overlapping boundaries at the geometric mid-points + float rBeta = sqrt(PETROSIAN_BETA); + for (int i = 3; radBet->data.F32[i-1] < radiusMax; i++) { + if (POWER_LAW_SPACING) { + radMin->data.F32[i] = radMax->data.F32[i-1]; + radMax->data.F32[i] = radMin->data.F32[i] * PETROSIAN_BETA; + radAlp->data.F32[i] = radMin->data.F32[i] / rBeta; + radBet->data.F32[i] = radMax->data.F32[i] * rBeta; + } else { + radMin->data.F32[i] = radMax->data.F32[i-1]; + radMax->data.F32[i] = radMin->data.F32[i] + 1; + float rMid = 0.5*(radMin->data.F32[i] + radMax->data.F32[i]); + radAlp->data.F32[i] = rMid * PETROSIAN_ALPHA; + radBet->data.F32[i] = rMid * PETROSIAN_BETA; + } + radMin->n = radMax->n = radAlp->n = radBet->n = i + 1; } - radMin->n = radMax->n = radAlp->n = radBet->n = nPts; // generate radial area-weighted mean radius & non-overlapping areas @@ -78,10 +91,11 @@ // XXX calculate area-weighted radius rather than asserting? - radBin->data.F32[i] = rBin; - area->data.F32[i] = M_PI * (rMax2 - rMin2); + binRad->data.F32[i] = rBin; + binArea->data.F32[i] = M_PI * (rMax2 - rMin2); - if (i > 2) { - radAlp->data.F32[i] = rBin*PETROSIAN_ALPHA; - radBet->data.F32[i] = rBin*PETROSIAN_BETA; + if (0) { + fprintf (stderr, "%3d %5.1f %5.1f : %5.1f : %5.1f %5.1f\n", + i, radAlp->data.F32[i], radMin->data.F32[i], binRad->data.F32[i], + radMax->data.F32[i], radBet->data.F32[i]); } } @@ -89,5 +103,5 @@ // storage vector for stats psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32); - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // integrate flux, radius for each of these bins. since flux is sorted by radius, @@ -107,7 +121,15 @@ // calculate the value for the nOut bin psVectorStats (stats, values, NULL, NULL, 0); - fluxBin->data.F32[nOut] = stats->sampleMedian; + // binSB->data.F32[nOut] = stats->sampleMedian; + binSB->data.F32[nOut] = stats->robustMedian; + binSBstdev->data.F32[nOut] = stats->robustStdev / sqrt(values->n); + + if (1) { + fprintf (stderr, "%3d %5.1f %5.1f : %5.1f %5.2f\n", + nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]); + } + nOut ++; - if (nOut >= nMax) break; + if (nOut >= radAlp->n) break; Rmin = radAlp->data.F32[nOut]; Rmax = radBet->data.F32[nOut]; @@ -122,13 +144,14 @@ psVectorAppend (values, flux->data.F32[i]); } - fluxBin->n = radBin->n = area->n = nOut; + binSB->n = binSBstdev->n = binRad->n = binArea->n = nOut; // XXX I think this misses the last radial bin -- do we care? // save the vectors - petrosian->radialBins = radBin; - petrosian->area = area; - petrosian->binnedFlux = fluxBin; + petrosian->radialBins = binRad; + petrosian->area = binArea; + petrosian->binSB = binSB; + petrosian->binSBstdev = binSBstdev; - psphotPetrosianVisualProfileRadii (radius, flux, radBin, fluxBin, 0.0); + psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, 0.0); psFree(radMin); Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c (revision 25178) @@ -1,79 +1,118 @@ # include "psphotInternal.h" -# define PETROSIAN_RATIO 0.05 +# define PETROSIAN_RATIO 0.2 +# define PETROSIAN_RADII 2.0 // generate the Petrosian radius and flux from the mean surface brightness (r_i) +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); + bool psphotPetrosianStats (pmPetrosian *petrosian) { - float petRadius, petFlux; + float petRadius, petFlux; - psVector *fluxBin = petrosian->binnedFlux; - psVector *radBin = petrosian->radialBins; - psVector *area = petrosian->area; + psVector *binSB = petrosian->binSB; + psVector *binSBstdev = petrosian->binSBstdev; + psVector *binRad = petrosian->radialBins; + psVector *area = petrosian->area; - psVector *fluxSum = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32); - psVector *refRadius= psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32); - psVector *petRatio = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32); - psVector *meanSB = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32); - psVector *areaSum = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32); + psVector *fluxSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *refRadius = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *petRatio = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *petRatioErr = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *meanSB = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *areaSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); - bool above = true; - float Fsum = 0.0; - float Asum = 0.0; - float Area = area->data.F32[0]; + bool anyPetro = false; + bool manyPetro = false; + bool above = true; + float Fsum = 0.0; + float Asum = 0.0; - int nOut = 0; - for (int i = 0; i < fluxBin->n; i++) { - // for nan bins, we keep the area for use with the next valid bin - if (!isfinite(fluxBin->data.F32[i])) { - Area += area->data.F32[i]; - continue; - } - Fsum += fluxBin->data.F32[i] * Area; - Asum += Area; - if (i+1 < fluxBin->n) { - Area = area->data.F32[i+1]; + int nOut = 0; + for (int i = 0; i < binSB->n; i++) { + // skip nan bins (do not contribute to flux or area) + if (!isfinite(binSB->data.F32[i])) continue; + + float Area = area->data.F32[i]; + Fsum += binSB->data.F32[i] * Area; + Asum += Area; + + float fluxInner = 0.5 * Area * binSB->data.F32[i]; + float areaInner = 0.5 * Area; + if (nOut > 0) { + fluxInner += fluxSum->data.F32[nOut-1]; + areaInner += areaSum->data.F32[nOut-1]; + } + + psVectorAppend(meanSB, (fluxInner / areaInner)); + psVectorAppend(petRatio, binSB->data.F32[i] / meanSB->data.F32[nOut]); + psVectorAppend(petRatioErr, binSBstdev->data.F32[i] / meanSB->data.F32[nOut]); + psVectorAppend(fluxSum, Fsum); + psVectorAppend(areaSum, Asum); + psVectorAppend(refRadius, binRad->data.F32[i]); + + if (1) { + fprintf (stderr, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n", + i, refRadius->data.F32[nOut], binSB->data.F32[i], meanSB->data.F32[nOut], petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], fluxSum->data.F32[nOut], areaSum->data.F32[nOut]); + } + + // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux + // we will keep and report the last (largest radius) value + if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO)) { + // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec) + if (i == 0) { + // assume Fmax @ R = 0.0 + petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + } else { + petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + } + above = false; + if (anyPetro) manyPetro = true; + anyPetro = true; + } + + // reset on transitions up, but do not re-calculate rad_90, flux_90 + if (!above && (petRatio->data.F32[nOut] >= PETROSIAN_RATIO)) { + above = true; + } + nOut ++; } - psVectorAppend(petRatio, Asum * fluxBin->data.F32[i] / Fsum); - psVectorAppend(fluxSum, Fsum); - psVectorAppend(meanSB, Fsum / Asum); - psVectorAppend(areaSum, Asum); - psVectorAppend(refRadius, radBin->data.F32[i]); + if (!anyPetro) { + // set default radius: + } - // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux - // we will keep and report the last (largest radius) value - if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO)) { - // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec) - if (i == 0) { - // assume Fmax @ R = 0.0 - petRadius = radBin->data.F32[nOut] * (PETROSIAN_RATIO - 1.0) / (petRatio->data.F32[nOut] - 1.0); - petFlux = fluxSum->data.F32[nOut] * (PETROSIAN_RATIO - 1.0) / (petRatio->data.F32[nOut] - 1.0); - } else { - petRadius = radBin->data.F32[nOut-1] + (radBin->data.F32[nOut] - radBin->data.F32[nOut-1]) * (PETROSIAN_RATIO - petRatio->data.F32[nOut-1]) / (petRatio->data.F32[nOut] - petRatio->data.F32[nOut-1]); - petFlux = fluxSum->data.F32[nOut-1] + (fluxSum->data.F32[nOut] - fluxSum->data.F32[nOut-1]) * (PETROSIAN_RATIO - petRatio->data.F32[nOut-1]) / (petRatio->data.F32[nOut] - petRatio->data.F32[nOut-1]); - } - above = false; + // now measure the flux within PETROSIAN_RADII * petRadius + float apRadius = PETROSIAN_RADII * petRadius; + for (int i = 0; i < refRadius->n; i++) { + // XXX use bisection to do this faster: + if (refRadius->data.F32[i] > apRadius) { + if (i == 0) { + psAbort ("does this case make any sense?"); + } else { + petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius); + break; + } + } } - - // reset on transitions up, but do not re-calculate rad_90, flux_90 - if (!above && (petRatio->data.F32[nOut] >= PETROSIAN_RATIO)) { - above = true; - } - nOut ++; - } - // save petRadius, petFlux - petrosian->petrosianRadius = petRadius; - petrosian->petrosianFlux = petFlux; + // save petRadius, petFlux + petrosian->petrosianRadius = petRadius; + petrosian->petrosianFlux = petFlux; - psphotPetrosianVisualStats (radBin, fluxBin, refRadius, meanSB, petRatio, fluxSum, petRadius, petFlux); + psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius); - psFree(fluxSum); - psFree(petRatio); - psFree(meanSB); - psFree(areaSum); + psFree(fluxSum); + psFree(petRatio); + psFree(meanSB); + psFree(areaSum); - return true; + return true; } + +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) { + float Y = Y0 + (Y1 - Y0) * (X - X0) / (X1 - X0); + return Y; +} + Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c (revision 25178) @@ -133,5 +133,5 @@ source = pmSourceAlloc(); source->peak = psphotLocalPeak(readout, Xo, Yo); - pmSourceDefinePixels (source, readout, Xo, Yo, 64); + pmSourceDefinePixels (source, readout, Xo, Yo, 128); } Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c (revision 25178) @@ -161,6 +161,8 @@ bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, psVector *refRadius, psVector *meanSB, - psVector *petRatio, psVector *fluxSum, - float petRadius, float petFlux) + psVector *petRatio, psVector *petRatioErr, + psVector *fluxSum, + float petRadius, float ratioForRadius, + float petFlux, float radiusForFlux) { Graphdata graphdata; @@ -229,5 +231,7 @@ graphdata.color = KapaColorByName ("black"); - pmVisualLimitsFromVectors (&graphdata, radBin, petRatio); + graphdata.ymax = +1.05; + graphdata.ymin = -0.05; + pmVisualLimitsFromVectors (&graphdata, radBin, NULL); KapaSetLimits (kapa, &graphdata); @@ -240,9 +244,12 @@ graphdata.ptype = 0; graphdata.size = 1.0; + graphdata.etype = 0x01; KapaPrepPlot (kapa, refRadius->n, &graphdata); KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x"); KapaPlotVector (kapa, refRadius->n, petRatio->data.F32, "y"); - - float ratio_90 = 0.1; + KapaPlotVector (kapa, refRadius->n, petRatioErr->data.F32, "dym"); + KapaPlotVector (kapa, refRadius->n, petRatioErr->data.F32, "dyp"); + graphdata.etype = 0; + graphdata.color = KapaColorByName ("red"); graphdata.style = 2; @@ -251,5 +258,5 @@ KapaPrepPlot (kapa, 1, &graphdata); KapaPlotVector (kapa, 1, &petRadius, "x"); - KapaPlotVector (kapa, 1, &ratio_90, "y"); + KapaPlotVector (kapa, 1, &ratioForRadius, "y"); // *** section 3: radius vs integrated flux @@ -283,5 +290,5 @@ graphdata.size = 2.0; KapaPrepPlot (kapa, 1, &graphdata); - KapaPlotVector (kapa, 1, &petRadius, "x"); + KapaPlotVector (kapa, 1, &radiusForFlux, "x"); KapaPlotVector (kapa, 1, &petFlux, "y"); Index: /branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c (revision 25178) @@ -4,5 +4,8 @@ // separations -// XXX choices for rMax and Nsec? +float psphotMeanSectorValue (psImage *image, int x, int y, float dL, float dW, float theta); +psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta); +psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW); +psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords); bool psphotRadialProfilesByAngles (pmPetrosian *petrosian, pmSource *source, int Nsec, float Rmax) { @@ -28,7 +31,10 @@ psVector *flux = psVectorAllocEmpty(Rmax, PS_TYPE_F32); - // start at Xo,Yo and find the x,y locations for r_i, theta where r_i increments by 1 pixel - // XXX at large radii ( > 10-15) use stats in a patch rather than sub-pixel interpolation - for (float r = 0; r < Rmax; r += 1.0) { + // Start at Xo,Yo and find the x,y locations for r_i, theta where r_i initially + // increments by 1 pixel. At large radii (r*dtheta > 2) use stats in a box rather than + // sub-pixel interpolation + + int dR = 1.0; + for (float r = 0; r < Rmax; r += dR) { float Xo = source->peak->xf + 0.5; @@ -38,18 +44,27 @@ float x = r * cos (theta) + Xo; float y = r * sin (theta) + Yo; - - if (x < 0) continue; - if (y < 0) continue; - if (x >= source->pixels->parent->numCols) continue; - if (y >= source->pixels->parent->numRows) continue; - - // value is NAN if we run off the image - float value = psImageInterpolatePixelBilinear(x, y, source->pixels); - - // keep the nan values so all vectors are matched - // if (isnan(value)) continue; - + dR = 2*(int)(0.5*r*sin(dtheta)) + 1; + + if (x < 0) goto badvalue; + if (y < 0) goto badvalue; + if (x >= source->pixels->parent->numCols) goto badvalue; + if (y >= source->pixels->parent->numRows) goto badvalue; + + float value = NAN; + if (dR < 2) { + // value is NAN if we run off the image + value = psImageInterpolatePixelBilinear(x, y, source->pixels); + } else { + value = psphotMeanSectorValue(source->pixels, x, y, dR, dR, theta); + } + + // keep the all values (even NAN) so all vectors are matched in length psVectorAppend (radius, r); psVectorAppend (flux, value); + continue; + + badvalue: + psVectorAppend (radius, r); + psVectorAppend (flux, NAN); } @@ -103,2 +118,112 @@ return true; } + +float psphotMeanSectorValue (psImage *image, int x, int y, float dL, float dW, float theta) { + + psVector *values = psphotBoxValues (image, x, y, dL, dW, theta); + + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + psVectorStats (stats, values, NULL, NULL, 0); + + float value = stats->sampleMedian; + + psFree (stats); + psFree (values); + + return value; +} + +psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta) { + + // extract pixels from a series of lines (from -0.5*dW to +0.5*dW) of length dL, + // centered on x0, y0 in parent pixel coordinates (not pixel indicies) + + float xs = x0 - image->col0 - 0.5*dL*cos(theta); + float ys = y0 - image->row0 - 0.5*dL*sin(theta); + + float xe = xs + 0.5*dL*cos(theta); + float ye = ys + 0.5*dL*sin(theta); + + psVector *values = psphotLineValues (image, xs, ys, xe, ye, (int) dW); + return values; +} + +/** + * identify the quadrant and draw the correct line + */ +psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW) { + + int FlipDirect, FlipCoords; + int X1, Y1, X2, Y2, dX, dY; + + /* rather than draw the line from float positions, we find the closest + integer end-points and draw the line between those pixels */ + + X1 = ROUND(x1); + Y1 = ROUND(y1); + X2 = ROUND(x2); + Y2 = ROUND(y2); + + dX = X2 - X1; + dY = Y2 - Y1; + + FlipCoords = (abs(dX) < abs(dY)); + FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2); + + psVector *values = NULL; + if (!FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X1, Y1, X2, Y2, dW, FALSE); + if ( FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X2, Y2, X1, Y1, dW, FALSE); + if (!FlipDirect && FlipCoords) values = psphotLineValuesBresen (image, Y1, X1, Y2, X2, dW, TRUE); + if ( FlipDirect && FlipCoords) values = psphotLineValuesBresen (image, Y2, X2, Y1, X1, dW, TRUE); + + return values; +} + +/** + * use the Bresenham line drawing technique + * integer-only Bresenham line-draw version which is fast + */ +psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) { + + int X, Y, dX, dY; + int e, e2; + + psVector *values = psVectorAllocEmpty(100, PS_TYPE_F32); + + dX = X2 - X1; + dY = Y2 - Y1; + + Y = Y1; + e = 0; + for (X = X1; X <= X2; X++) { + if (X > 0) { + if (swapcoords) { + if (X >= image->numRows) continue; + for (int y = Y - dW; y <= Y + dW; y++) { + if (y < 0) continue; + if (y >= image->numCols) continue; + psVectorAppend(values, image->data.F32[X][y]); + } + } else { + if (X >= image->numCols) continue; + for (int y = Y - dW; y <= Y + dW; y++) { + if (y < 0) continue; + if (y >= image->numRows) continue; + psVectorAppend(values, image->data.F32[y][X]); + } + } + } + e += dY; + e2 = 2 * e; + if (e2 > dX) { + Y++; + e -= dX; + } + if (e2 < -dX) { + Y--; + e += dX; + } + } + return values; +} + Index: /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (revision 25177) +++ /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (revision 25178) @@ -31,7 +31,7 @@ // examine data in the two ranges Fm - Fo and Fo - Fp to define the bin size // XXX reconsider the fractional isophote value - float Fm = fluxMin + 0.25*fluxRange; - float Fp = fluxMin + 0.75*fluxRange; - float Fo = fluxMin + 0.50*fluxRange; + float Fm = fluxMin + 0.10*fluxRange; + float Fo = fluxMin + 0.25*fluxRange; + float Fp = fluxMin + 0.50*fluxRange; int Rbin = 1;