Index: trunk/psphot/src/psphotKronIterate.c =================================================================== --- trunk/psphot/src/psphotKronIterate.c (revision 34545) +++ trunk/psphot/src/psphotKronIterate.c (revision 34559) @@ -1,10 +1,9 @@ # include "psphotInternal.h" -bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow); -bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels, bool measureRadius); - - -bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) -{ +bool psphotKronRadiusMeasure (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal, psImage *smoothedPixels); +bool psphotKronFluxMeasure (pmSource *source, psImageMaskType maskVal); + + +bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) { bool status = true; @@ -170,19 +169,5 @@ } - // generate the window image: multiply the flux by this to downweight neighbors - // XXX: we don't need this image if we aren't going to apply the window - psImage *kronWindow = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32); - psImageInit (kronWindow, 1.0); - - // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // set a window function for each source based on the moments - // (this skips really bad sources (no peak, no moments, DEFECT) - psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW); - } - + // pass 3 we only measure fluxes for matched sources so we don't need smooth pixels or multiple iterations if (pass == 3) { KRON_SMOOTH = false; @@ -240,5 +225,4 @@ psArrayAdd(job->args, 1, readout); - psArrayAdd(job->args, 1, kronWindow); psArrayAdd(job->args, 1, cells->data[j]); // sources psArrayAdd(job->args, 1, smoothedImage); @@ -248,6 +232,4 @@ PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32); PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32); PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA, PS_TYPE_F32); PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32); @@ -286,5 +268,4 @@ } psFree (cellGroups); - psFree (kronWindow); if (KRON_SMOOTH) { for (int i = 0; i < sources->n; i++) { @@ -303,17 +284,14 @@ pmReadout *readout = job->args->data[0]; - psImage *kronWindow = job->args->data[1]; - psArray *sources = job->args->data[2]; - psImage *smoothedImage = job->args->data[3]; - psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); - psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); - float RADIUS = PS_SCALAR_VALUE(job->args->data[6],F32); - float MIN_KRON_RADIUS = PS_SCALAR_VALUE(job->args->data[7],F32); - int KRON_ITERATIONS = PS_SCALAR_VALUE(job->args->data[8],S32); - bool KRON_APPLY_WEIGHT = PS_SCALAR_VALUE(job->args->data[9],S32); - bool KRON_APPLY_WINDOW = PS_SCALAR_VALUE(job->args->data[10],S32); - float KRON_SMOOTH_SIGMA = PS_SCALAR_VALUE(job->args->data[11],F32); - float KRON_SB_MIN_DIVISOR = PS_SCALAR_VALUE(job->args->data[12],F32); - int pass = PS_SCALAR_VALUE(job->args->data[13],S32); + psArray *sources = job->args->data[1]; + psImage *smoothedImage = job->args->data[2]; + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); + float RADIUS = PS_SCALAR_VALUE(job->args->data[5],F32); + float MIN_KRON_RADIUS = PS_SCALAR_VALUE(job->args->data[6],F32); + int KRON_ITERATIONS = PS_SCALAR_VALUE(job->args->data[7],S32); + float KRON_SMOOTH_SIGMA = PS_SCALAR_VALUE(job->args->data[8],F32); + float KRON_SB_MIN_DIVISOR = PS_SCALAR_VALUE(job->args->data[9],F32); + int pass = PS_SCALAR_VALUE(job->args->data[10],S32); bool measureRadius = true; @@ -327,18 +305,19 @@ pmSource *source = sources->data[i]; if (!source->peak) continue; // XXX how can we have a peak-less source? - - if (pass == 3) { - // in pass 3 we only measure the flux for matched sources - if (!(source->mode2 & PM_SOURCE_MODE2_MATCHED)) { - continue; - } - } + if (!source->moments) continue; + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + // in pass 3 we only measure the flux for matched sources + if ((pass == 3) && !(source->mode2 & PM_SOURCE_MODE2_MATCHED)) continue; # if (0) # define TEST_X 653 # define TEST_Y 466 - if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { - fprintf (stderr, "test object\n"); - } + if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { + fprintf (stderr, "test object\n"); + } # undef TEST_X # undef TEST_Y @@ -347,5 +326,5 @@ if (measureRadius) { // check status of this source's moments - if (!source->moments) continue; + // XXX: I don't think that we have to apply these restrictions since we dropped the window function if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue; if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; @@ -354,10 +333,7 @@ if (!isfinite(source->moments->Mrf)) { // Once we save a bad Mrf measurement we give up on this source - // checking here allows us to avoid adding and subtracting the model + // XXX: is this the right thing to do? continue; } - - // skip saturated stars modeled with a radial profile - if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; // replace object in image @@ -402,54 +378,40 @@ float windowRadius = PS_MAX(RADIUS, maxWindow); -#ifdef REVERT_ON_BAD_MEASURMENT - // save previous measurements. We might revert back to them if this round fails - float MrfPrior = source->moments->Mrf; - float KronFluxPrior = source->moments->KronFlux; - float KronFluxErrPrior = source->moments->KronFluxErr; -#endif - // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); - psAssert (source->pixels, "WTF?"); + psAssert (source->pixels, "redefine pixels failed?"); if (extend && smoothedPixels) { psFree(source->tmpPtr); smoothedPixels = psImageSubset(smoothedImage, source->region); - psAssert (smoothedPixels, "WTF?"); + psAssert (smoothedPixels, "redefine smoothed pixels failed?"); source->tmpPtr = (psPtr) smoothedPixels ; } - - // clear the window function for this source based on the moments - // Note: this function also applies cuts on the source and returns false if it - // does not meet the requirements for measuring the Kron Radius or Magnitude. - // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false - if (!measureRadius || psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) { - - // this function populates moments->Mrf,KronFlux,KronFluxErr - psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, - KRON_APPLY_WEIGHT, smoothedPixels, measureRadius); - - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // set a window function for each source based on the moments - psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW); - } - -#ifdef REVERT_ON_BAD_MEASUREMENT - // on pass 2 if we get an invalid measurement on a pass 1 source where we had a good one previously - // in pass 1 keep that measurement - bool reverted = false; - if (pass > 1 && !isfinite(source->moments->Mrf) && (source->mode2 & PM_SOURCE_MODE2_PASS1_SRC)) { - source->moments->Mrf = MrfPrior; // This is finite otherwise we wouldn't have gotten here - source->moments->KronFlux = KronFluxPrior; - source->moments->KronFluxErr = KronFluxErrPrior; - reverted = true; - } -#endif + bool measureFlux = true; + if (measureRadius) { + measureFlux = psphotKronRadiusMeasure (source, windowRadius, MIN_KRON_RADIUS, maskVal, smoothedPixels); + } + + if (measureFlux) { + // Make sure the sources images are large enough for the measured Kron Radius + bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, + 2.5 * source->moments->Mrf + 2); + + psAssert (source->pixels, "redefine pixels failed?"); + if (extend && smoothedPixels) { + psFree(source->tmpPtr); + smoothedPixels = psImageSubset(smoothedImage, source->region); + psAssert (smoothedPixels, "redefine smoothed pixels failed?"); + source->tmpPtr = (psPtr) smoothedPixels ; + } + + // this function populates moments->Mrf,KronFlux,KronFluxErr, KronFinner, and KronFouter + psphotKronFluxMeasure (source, maskVal); + } + + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + #ifdef DUMP_KRS -#ifndef REVERT_ON_BAD_MEASUREMENT - bool reverted = false; -#endif - fprintf(dumpFile, "%7d %1d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, reverted, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID); + fprintf(dumpFile, "%7d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID); #endif @@ -467,6 +429,6 @@ } -bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, - bool applyWeight, psImage *smoothedPixels, bool measureRadius) { +bool psphotKronRadiusMeasure (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal, + psImage *smoothedPixels) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -476,5 +438,4 @@ psF32 R2 = PS_SQR(radius); - float rsigma2 = applyWeight ? 0.5 / R2 : 0; // a note about coordinates: coordinates of objects throughout psphot refer to the primary @@ -489,5 +450,4 @@ // Xn = SUM (x - xc)^n * (z - sky) - psF32 RF = 0.0; psF32 RS = 0.0; @@ -505,7 +465,4 @@ psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage - int Xwo = source->pixels->col0; - int Ywo = source->pixels->row0; - psF32 **vPix; if (smoothedPixels) { @@ -514,10 +471,8 @@ vPix = source->pixels->data.F32; } - psF32 **vWin = kronWindow->data.F32; - psF32 **vWgt = source->variance->data.F32; psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; - if (measureRadius) { + { for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -551,18 +506,10 @@ if (r2 > R2) continue; - // flux * window - float z = r2 * rsigma2; - assert (z >= 0.0); - - // weight by window image and wide Gaussian - float weight1 = vWin[row+Ywo][col+Xwo]*exp(-z); - float weight2 = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z); - - float fDiff1 = vPix[row][col]*weight1; - float fDiff2 = vPix[yFlip][xFlip]*weight2; + float fDiff1 = vPix[row][col]; + float fDiff2 = vPix[yFlip][xFlip]; float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2)); - // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?) + // Kron Flux uses the 1st radial moment psF32 rf = pDiff * sqrt(r2); psF32 rs = 0.5 * (fDiff1 + fDiff2); @@ -589,16 +536,71 @@ source->moments->Mrf = Mrf; } - - // Calculate the Kron magnitude (make this block optional?) - float radKron = 2.5 * source->moments->Mrf; - float radKron2 = radKron*radKron; - - int nKronPix = 0; + return true; +} + +bool psphotKronFluxMeasure(pmSource *source, psImageMaskType maskVal) { + + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + + // a note about coordinates: coordinates of objects throughout psphot refer to the primary + // image coordinates. the source->pixels image has an offset relative to its parent of + // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in + // this subimage. we subtract off the peak coordinates, adjusted to this subimage, to have + // minimal round-off error in the sums. since these values are subtracted just to minimize + // the dynamic range and are added back below, the exact value does not matter. these are + // (int) so they can be used in the image index below. + + // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates + // Xn = SUM (x - xc)^n * (z - sky) + + // the peak position is less accurate but less subject to extreme deviations + float dX = source->moments->Mx - source->peak->xf; + float dY = source->moments->My - source->peak->yf; + float dR = hypot(dX, dY); + float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf; + float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf; + + // center of mass in subimage. Note: the calculation below uses pixel index, so we correct + // xCM, yCM from pixel coords to pixel index here. + psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage + psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage + + psF32 **vPix = source->pixels->data.F32; + psF32 **vWgt = source->variance->data.F32; + + psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + + source->moments->KronFlux = NAN; + source->moments->KronFluxErr = NAN; + source->moments->KronFinner = NAN; + source->moments->KronFouter = NAN; + + // Calculate the Kron fluxes + float radKinner = 1.0 * source->moments->Mrf; + float radKron = 2.5 * source->moments->Mrf; + float radKouter = 4.0 * source->moments->Mrf; + + float limitRadius = MIN (radKouter, source->windowRadius); + if (radKouter > source->windowRadius) { + // This happens but the measurement isn't important enough to allocate the extra pixels + // psWarning ("outer kron radius: %f is larger than windowRadius: %f for %d\n", + // radKouter, source->windowRadius, source->id); + limitRadius = MIN (radKron, source->windowRadius); + } + if (radKron > source->windowRadius) { + // caller should have prevented this from happening + psWarning ("kron radius: %f is larger than windowRadius: %f for %d\n", + radKron, source->windowRadius, source->id); + return false; + } + float Sum = 0.0; float Var = 0.0; - float Win = 0.0; - - // set vPix to the source pixels (it may have been set to the - // smoothed image above) + float SumInner = 0.0; + float SumOuter = 0.0; + + // set vPix to the source pixels (it may have been set to the smoothed image above) vPix = source->pixels->data.F32; @@ -606,5 +608,5 @@ psF32 yDiff = row - yCM; - if (fabs(yDiff) > radKron) continue; + if (fabs(yDiff) > limitRadius) continue; for (psS32 col = 0; col < source->pixels->numCols ; col++) { @@ -614,20 +616,22 @@ psF32 xDiff = col - xCM; - if (fabs(xDiff) > radKron) continue; - - // radKron is just a function of (xDiff, yDiff) + if (fabs(xDiff) > limitRadius) continue; + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); - if (r2 > radKron2) continue; - - float weight1 = vWin[row+Ywo][col+Xwo]; - float fDiff1 = vPix[row][col]*weight1; - - float pDiff = fDiff1; - psF32 wDiff = vWgt[row][col] * weight1; - - Sum += pDiff; - Var += wDiff; - Win += weight1; - nKronPix ++; + psF32 r = sqrt(r2); + + float pDiff = vPix[row][col]; + psF32 wDiff = vWgt[row][col]; + + if (r > radKinner && r < radKron) { + SumInner += pDiff; + } + if (r < radKron) { + Sum += pDiff; + Var += wDiff; + } + if (r > radKron && r < radKouter) { + SumOuter += pDiff; + } } } @@ -635,89 +639,11 @@ source->moments->KronFlux = Sum; source->moments->KronFluxErr = sqrt(Var); + source->moments->KronFinner = SumInner; + + // only save radKouter if the radius is inside the integration radius limit + if (radKouter <= limitRadius) { + source->moments->KronFouter = SumOuter; + } return true; } - -bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) { - - if (!source) return false; - if (!source->peak) return false; // XXX how can we have a peak-less source? - if (!source->moments) return false; - if (source->type == PM_SOURCE_TYPE_DEFECT) return false; - if (source->type == PM_SOURCE_TYPE_SATURATED) return false; - if (source->mode2 & PM_SOURCE_MODE2_MATCHED) return true; - - if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) return false; - if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false; - psAssert(kronWindow, "need a window"); - - // XXX: If we are not applying the window then we don't need to check for valid Mrf here. - // We should give the this module a chance to measure a good value. - // However experiments show that it hardly ever succeeds in getting a better value - if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false; - - if (!applyWindow) { - return true; - } - - // we have a source with moments Mx, My, Mxx, Mxy, Myy. we just need to define a Gaussian that has - // these values (and peak of 1.0) - - int Nx = kronWindow->numCols; - int Ny = kronWindow->numRows; - - float Xo = source->moments->Mx; - float Yo = source->moments->My; - - psEllipseMoments moments; - moments.x2 = source->moments->Mxx; - moments.y2 = source->moments->Myy; - moments.xy = source->moments->Mxy; - - psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0); - if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) { - // Shall we log a proper warning? This happens often with matched sources (forced photometry) - // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta); - return false; - } - - // Why this factor of 0.5 ? - float scale = 0.5 * source->moments->Mrf / axes.major; - axes.major *= scale; - axes.minor *= scale; - - psEllipseShape shape = psEllipseAxesToShape(axes); - if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) { - // Shall we log a proper warning? This happens often with matched sources (forced photometry) - // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy); - return false; - } - - float Smajor = axes.major; - - int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1); - int maxX = PS_MIN(PS_MAX(Xo + 5*Smajor, 0), Nx - 1); - int minY = PS_MIN(PS_MAX(Yo - 5*Smajor, 0), Ny - 1); - int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1); - - float rMxx = 0.5 / PS_SQR(shape.sx); - float rMyy = 0.5 / PS_SQR(shape.sy); - float Sxy = -1. * shape.sxy; // factor of -1 is included to match the previous window function - // implementation. XXX: Is this correct? - - for (int iy = minY; iy < maxY; iy++) { - for (int ix = minX; ix < maxX; ix++) { - - float dX = (ix + 0.5 - Xo); - float dY = (iy + 0.5 - Yo); - - float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY; - - float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z)); - - kronWindow->data.F32[iy][ix] *= f; - } - } - - return true; -} Index: trunk/psphot/src/psphotSetThreads.c =================================================================== --- trunk/psphot/src/psphotSetThreads.c (revision 34545) +++ trunk/psphot/src/psphotSetThreads.c (revision 34559) @@ -30,5 +30,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 14); + task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 11); task->function = &psphotKronIterate_Threaded; psThreadTaskAdd(task);