Index: branches/eam_branches/ipp-20110404/psphot/src/psphotKronMasked.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotKronMasked.c (revision 31328) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotKronMasked.c (revision 31337) @@ -1,3 +1,5 @@ # include "psphotInternal.h" + +bool psphotKronMag (pmSource *source, float radius, psImageMaskType maskVal); bool psphotKronMasked (pmConfig *config, const pmFPAview *view, const char *filerule) @@ -29,5 +31,5 @@ psAssert (detections, "missing detections?"); - psArray *sources = detections->newSources; + psArray *sources = detections->allSources; psAssert (sources, "missing sources?"); @@ -58,4 +60,9 @@ if (!status) { nThreads = 0; + } + + float RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS"); + if (!status) { + RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); } @@ -84,4 +91,17 @@ if (!source->moments) continue; + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + psphotKronMag (source, RADIUS, maskVal); + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + continue; + + // XXX skip this code // generate the pixel masks // int Xo = source->moments->Mx; @@ -117,2 +137,137 @@ return true; } + +bool psphotKronMag (pmSource *source, float radius, 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); + PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); + + psF32 Sum = 0.0; + psF32 Var = 0.0; + psF32 R2 = PS_SQR(radius); + + // 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) + + psF32 RF = 0.0; + psF32 RS = 0.0; + +# if (1) + 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; +# else + float Xo = source->moments->Mx; + float Yo = source->moments->My; +# endif + + // 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 + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 yDiff = row - yCM; + if (fabs(yDiff) > radius) continue; + + psF32 *vPix = source->pixels->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++) { + if (vMsk) { + if (*vMsk & maskVal) { + vMsk++; + continue; + } + vMsk++; + } + if (isnan(*vPix)) continue; + + psF32 xDiff = col - xCM; + if (fabs(xDiff) > radius) continue; + + // radius is just a function of (xDiff, yDiff) + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > R2) continue; + + psF32 pDiff = *vPix; + + Sum += pDiff; + + // Kron Flux uses the 1st radial moment (NOT Gaussian windowed) + psF32 rf = pDiff * sqrt(r2); + psF32 rs = pDiff; + + RF += rf; + RS += rs; + } + } + + // Saturate the 1st radial moment + float Mrf = MIN(radius, MAX(0.25*radius, RF/RS)); + + // Calculate the Kron magnitude (make this block optional?) + float radKron = 2.5*Mrf; + float radKron2 = radKron*radKron; + + int nKronPix = 0; + Sum = Var = 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); + if (r2 > radKron2) continue; + + psF32 pDiff = *vPix; + psF32 wDiff = *vWgt; + + Sum += pDiff; + Var += wDiff; + nKronPix ++; + } + } + // XXX for a test, save the old values here: + source->moments->KronCore = source->moments->KronFlux; + source->moments->KronCoreErr = source->moments->KronFluxErr; + + source->moments->KronFlux = Sum * M_PI * PS_SQR(radKron) / nKronPix; + source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix; + + return true; +} Index: branches/eam_branches/ipp-20110404/psphot/src/psphotLoadSRCTEXT.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotLoadSRCTEXT.c (revision 31328) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotLoadSRCTEXT.c (revision 31337) @@ -52,5 +52,5 @@ // NOTE: most of these values are irrelevant for loaded source positions - source->seq = 0; + source->seq = source->id; PAR[PM_PAR_XPOS] = X; PAR[PM_PAR_YPOS] = Y; Index: branches/eam_branches/ipp-20110404/psphot/src/psphotMergeSources.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotMergeSources.c (revision 31328) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotMergeSources.c (revision 31337) @@ -137,4 +137,9 @@ pmSource *source = extSourcesTXT->data[i]; source->mode |= PM_SOURCE_MODE_EXTERNAL; + + // the supplied peak flux needs to be re-normalized + source->peak->rawFlux = 1.0; + source->peak->smoothFlux = 1.0; + source->peak->detValue = 1.0; // drop the loaded source modelPSF Index: branches/eam_branches/ipp-20110404/psphot/src/psphotReadout.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotReadout.c (revision 31328) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotReadout.c (revision 31337) @@ -146,6 +146,4 @@ } - // psphotKronMasked(config, view, filerule); - // find blended neighbors of very saturated stars (detections->newSources) // if (!psphotDeblendSatstars (config, view, filerule)) { @@ -201,4 +199,7 @@ psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources) psphotDumpChisqs (config, view, filerule); + + // XXX re-measure the kron mags with models subtracted + psphotKronMasked(config, view, filerule); // identify CRs and extended sources (only unmeasured sources are measured) @@ -311,4 +312,7 @@ pass1finish: + // XXX re-measure the kron mags with models subtracted + psphotKronMasked(config, view, filerule); + // measure source size for the remaining sources // NOTE: applies only to NEW (unmeasured) sources