Index: /trunk/ippconfig/recipes/psphot.config =================================================================== --- /trunk/ippconfig/recipes/psphot.config (revision 33836) +++ /trunk/ippconfig/recipes/psphot.config (revision 33837) @@ -38,4 +38,5 @@ SKY_CLIP_SIGMA F32 2.0 # statistic used to measure background SKY_SIG F32 1.0 # optional sky error for +SKY_SLOPE_MIN F32 3.0 # exit radial profile loop when abs(slope) is less than this # allowed values for SKY_STAT: @@ -188,4 +189,7 @@ KRON_ITERATIONS S32 2 +KRON_APPLY_WEIGHT BOOL TRUE +KRON_APPLY_WINDOW BOOL TRUE +KRON_SMOOTH BOOL FALSE # Extended source fit parameters Index: /trunk/psphot/src/psphotKronIterate.c =================================================================== --- /trunk/psphot/src/psphotKronIterate.c (revision 33836) +++ /trunk/psphot/src/psphotKronIterate.c (revision 33837) @@ -1,15 +1,10 @@ # include "psphotInternal.h" -# ifndef ROUND -# define ROUND(X) ((int) ((X) + 0.5*SIGN(X))) -# endif - -bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert); -bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal); + +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, bool smooth); bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule) { bool status = true; - - // return true; fprintf (stdout, "\n"); @@ -83,4 +78,23 @@ if (!status) { MIN_KRON_RADIUS = 0.25*RADIUS; + } + + int KRON_ITERATIONS = psMetadataLookupS32 (&status, recipe, "KRON_ITERATIONS"); + if (!status) { + KRON_ITERATIONS = 1; + } + + bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT"); + if (!status) { + KRON_APPLY_WEIGHT = true; + } + + bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW"); + if (!status) { + KRON_APPLY_WINDOW = false; + } + bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH"); + if (!status) { + KRON_SMOOTH = false; } @@ -114,5 +128,5 @@ // 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); + psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW); } @@ -140,7 +154,11 @@ PS_ARRAY_ADD_SCALAR(job->args, RADIUS, PS_TYPE_F32); 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, (psS32) KRON_SMOOTH, PS_TYPE_S32); // set this to 0 to run without threading -# if (1) +# if (0) if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); @@ -188,7 +206,10 @@ float RADIUS = PS_SCALAR_VALUE(job->args->data[5],F32); float MIN_KRON_RADIUS = PS_SCALAR_VALUE(job->args->data[6],F32); - - // XXX TEST : set iteration to 1 - for (int j = 0; j < 1; j++) { + int KRON_ITERATIONS = PS_SCALAR_VALUE(job->args->data[7],S32); + bool KRON_APPLY_WEIGHT = PS_SCALAR_VALUE(job->args->data[8],S32); + bool KRON_APPLY_WINDOW = PS_SCALAR_VALUE(job->args->data[9],S32); + bool KRON_SMOOTH = PS_SCALAR_VALUE(job->args->data[10],S32); + + for (int j = 0; j < KRON_ITERATIONS; j++) { for (int i = 0; i < sources->n; i++) { @@ -196,6 +217,8 @@ if (!source->peak) continue; // XXX how can we have a peak-less source? - // allocate space for moments + // check status of this source's moments if (!source->moments) continue; + if (!source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; // replace object in image @@ -219,11 +242,12 @@ // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS // if we lack the skyRadius (eg MATCHED sources), go to the default value - float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS; + float maxWindow; + if (j == 0) { + maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS; + } else { + maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS; + } float windowRadius = PS_MAX(RADIUS, maxWindow); -#ifdef notdef - fprintf(stderr, "Redefining pixels for source: %d %4d %4d new radius: %f\n", - i, source->peak->x, source->peak->y, windowRadius+2); -#endif // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); @@ -231,12 +255,13 @@ // clear the window function for this source based on the moments - psphotKronWindowSetSource (source, kronWindow, (j > 0), false); - - // this function populates moments->Mrf,KronFlux,KronFluxErr - psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal); - 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); + if (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, KRON_SMOOTH); + 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); + } // if we subtracted it above, re-subtract the object, leave local sky @@ -249,5 +274,6 @@ } -bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal) { +bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, +bool applyWeight, bool smooth) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -257,5 +283,5 @@ psF32 R2 = PS_SQR(radius); - float rsigma2 = 0.5 / PS_SQR(radius/2.0); + float rsigma2 = applyWeight ? 0.5 / R2 : 0; // a note about coordinates: coordinates of objects throughout psphot refer to the primary @@ -288,5 +314,13 @@ int Ywo = source->pixels->row0; - psF32 **vPix = source->pixels->data.F32; + psF32 **vPix; + psImage *smoothedImage = NULL;; + if (smooth) { + smoothedImage = psImageCopy(NULL, source->pixels, PS_TYPE_F32); + psImageSmooth(smoothedImage, 1.7, 2); + vPix = smoothedImage->data.F32; + } else { + vPix = source->pixels->data.F32; + } psF32 **vWin = kronWindow->data.F32; psF32 **vWgt = source->variance->data.F32; @@ -361,4 +395,8 @@ float Var = 0.0; float Win = 0.0; + + // set vPix to the source pixels (it may have been set to the + // smoothed image above) + vPix = source->pixels->data.F32; for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -395,9 +433,12 @@ source->moments->KronFlux = Sum; source->moments->KronFluxErr = sqrt(Var); + if (smoothedImage) { + psFree(smoothedImage); + } return true; } -bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert) { +bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) { if (!source) return false; @@ -406,5 +447,13 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) return false; if (source->type == PM_SOURCE_TYPE_SATURATED) return false; + 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"); + + if (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 @@ -417,19 +466,31 @@ float Yo = source->moments->My; - float Mxx = source->moments->Mxx; - float Mxy = source->moments->Mxy; - float Myy = source->moments->Myy; - - float Mmajor = 0.5*(Mxx + Myy) + 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy)); - float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy)); - - // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width - float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor; - - float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2 - float Sxy = Mxy / (scale * Mmajor * Mminor); - float Syy = scale * Mmajor * Mminor / Mxx; // sigma_y^2 - - float Smajor = sqrt(Mmajor); + 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; + } + + float scale = fabs(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 Sxx = PS_SQR(shape.sx); + float Syy = PS_SQR(shape.sy); + float Sxy = -1. * shape.sxy; + float Smajor = axes.major; int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1); @@ -440,5 +501,5 @@ float rMxx = 0.5 / Sxx; float rMyy = 0.5 / Syy; - + for (int iy = minY; iy < maxY; iy++) { for (int ix = minX; ix < maxX; ix++) { @@ -451,7 +512,8 @@ float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z)); - kronWindow->data.F32[iy][ix] *= f; + kronWindow->data.F32[iy][ix] *= f; } } + return true; }