Index: trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- trunk/psModules/src/objects/pmSourceMoments.c (revision 27531) +++ trunk/psModules/src/objects/pmSourceMoments.c (revision 29004) @@ -22,16 +22,23 @@ #include #include + #include "pmHDU.h" #include "pmFPA.h" -#include "pmFPAMaskWeight.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" #include "pmMoments.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.h" +#include "pmModelFuncs.h" #include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" #include "pmSource.h" @@ -54,4 +61,7 @@ # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) +static bool beVerbose = false; +void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; } + bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal) { @@ -61,11 +71,15 @@ PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); - // use sky from moments if defined, 0.0 otherwise + // use sky from moments if defined, 0.0 otherwise + + // XXX this value comes from the sky model at the source center, and tends to over-estimate + // the sky in the vicinity of bright sources. we are better off assuming the model worked + // well: psF32 sky = 0.0; - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } else { - sky = source->moments->Sky; - } + // XXX if (source->moments == NULL) { + // XXX source->moments = pmMomentsAlloc(); + // XXX } else { + // XXX sky = source->moments->Sky; + // XXX } // First Pass: calculate the first moments (these are subtracted from the coordinates below) @@ -131,9 +145,11 @@ psF32 wDiff = *vWgt; - // skip pixels below specified significance level. this is allowed, but should be - // avoided -- the over-weights the wings of bright stars compared to those of faint - // stars. - if (PS_SQR(pDiff) < minSN2*wDiff) continue; - // if (pDiff < 0) continue; // XXX : MWV says I should include < 0.0 valued points... + // skip pixels below specified significance level. for a PSFs, this + // over-weights the wings of bright stars compared to those of faint stars. + // for the estimator used for extended source analysis (where the window + // function is allowed to be arbitrarily large), we need to clip to avoid + // negative second moments. + if (PS_SQR(pDiff) < minSN2*wDiff) continue; // + if ((minSN > 0.0) && (pDiff < 0)) continue; // // Apply a Gaussian window function. Be careful with the window function. S/N @@ -197,4 +213,7 @@ // Xn = SUM (x - xc)^n * (z - sky) + psF32 RF = 0.0; + psF32 RH = 0.0; + psF32 RS = 0.0; psF32 XX = 0.0; psF32 XY = 0.0; @@ -244,5 +263,6 @@ if (r > radius) continue; - psF32 pDiff = *vPix - sky; + psF32 fDiff = *vPix - sky; + psF32 pDiff = fDiff; psF32 wDiff = *vWgt; @@ -257,5 +277,7 @@ if (sigma > 0.0) { // XXX a lot of extra flops; can we do pre-calculate? - psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + // XXX we were re-calculating r2 (maybe the compiler caught this?) + // psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + psF32 z = r2 * rsigma2; assert (z >= 0.0); psF32 weight = exp(-z); @@ -266,4 +288,17 @@ Sum += pDiff; + +# if (1) +# if (0) + if (fDiff < 0) continue; +# endif + psF32 rf = r * fDiff; + psF32 rh = sqrt(r) * fDiff; + psF32 rs = fDiff; +# else + psF32 rf = r * pDiff; + psF32 rh = sqrt(r) * pDiff; + psF32 rs = pDiff; +# endif psF32 x = xDiff * pDiff; @@ -284,4 +319,8 @@ psF32 xyyy = xDiff * yyy / r2; psF32 yyyy = yDiff * yyy / r2; + + RF += rf; + RH += rh; + RS += rs; XX += xx; @@ -302,4 +341,7 @@ } + source->moments->Mrf = RF/RS; + source->moments->Mrh = RH/RS; + source->moments->Mxx = XX/Sum; source->moments->Mxy = XY/Sum; @@ -317,12 +359,71 @@ source->moments->Myyyy = YYYY/Sum; - // if (source->moments->Mxx < 0) { - // fprintf (stderr, "error: neg second moment??\n"); - // } - // if (source->moments->Myy < 0) { - // fprintf (stderr, "error: neg second moment??\n"); - // } - - psTrace ("psModules.objects", 4, "Mxx: %f Mxy: %f Myy: %f Mxxx: %f Mxxy: %f Mxyy: %f Myyy: %f Mxxxx: %f Mxxxy: %f Mxxyy: %f Mxyyy: %f Mxyyy: %f\n", + // Calculate the Kron magnitude (make this block optional?) + // float radKron = 2.5*source->moments->Mrf; + float radKinner = 1.0*source->moments->Mrf; + float radKron = 2.5*source->moments->Mrf; + float radKouter = 4.0*source->moments->Mrf; + + int nKronPix = 0; + Sum = Var = 0.0; + float SumInner = 0.0; + float SumOuter = 0.0; + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 yDiff = row - yCM; + if (fabs(yDiff) > radKron) continue; + + psF32 *vPix = source->pixels->data.F32[row]; + psF32 *vWgt = source->variance->data.F32[row]; + psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; + + for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { + if (vMsk) { + if (*vMsk & maskVal) { + vMsk++; + continue; + } + vMsk++; + } + if (isnan(*vPix)) continue; + + psF32 xDiff = col - xCM; + if (fabs(xDiff) > radKron) continue; + + // radKron is just a function of (xDiff, yDiff) + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + psF32 r = sqrt(r2); + + psF32 pDiff = *vPix - sky; + psF32 wDiff = *vWgt; + + // skip pixels below specified significance level. this is allowed, but should be + // avoided -- the over-weights the wings of bright stars compared to those of faint + // stars. + if (PS_SQR(pDiff) < minSN2*wDiff) continue; + + if (r < radKron) { + Sum += pDiff; + Var += wDiff; + nKronPix ++; + // if (beVerbose) fprintf (stderr, "mome: %d %d %f %f %f\n", col, row, sky, *vPix, Sum); + } + + if ((r > radKinner) && (r < radKron)) { + SumInner += pDiff; + } + if ((r > radKron) && (r < radKouter)) { + SumOuter += pDiff; + } + } + } + source->moments->KronFlux = Sum; + source->moments->KronFinner = SumInner; + source->moments->KronFouter = SumOuter; + source->moments->KronFluxErr = sqrt(Var); + + psTrace ("psModules.objects", 4, "Mrf: %f KronFlux: %f Mxx: %f Mxy: %f Myy: %f Mxxx: %f Mxxy: %f Mxyy: %f Myyy: %f Mxxxx: %f Mxxxy: %f Mxxyy: %f Mxyyy: %f Mxyyy: %f\n", + source->moments->Mrf, source->moments->KronFlux, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy,