Index: trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- trunk/psModules/src/objects/pmSourceMoments.c (revision 31153) +++ trunk/psModules/src/objects/pmSourceMoments.c (revision 31451) @@ -66,5 +66,5 @@ // if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid) -bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal) +bool pmSourceMoments(pmSource *source, float radius, float sigma, float minSN, float minKronRadius, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -74,5 +74,5 @@ // this function assumes the sky has been well-subtracted for the image - psF32 sky = 0.0; + float sky = 0.0; if (source->moments == NULL) { @@ -80,11 +80,11 @@ } - psF32 Sum = 0.0; - psF32 Var = 0.0; - psF32 SumCore = 0.0; - psF32 VarCore = 0.0; - psF32 R2 = PS_SQR(radius); - psF32 minSN2 = PS_SQR(minSN); - psF32 rsigma2 = 0.5 / PS_SQR(sigma); + float Sum = 0.0; + float Var = 0.0; + float SumCore = 0.0; + float VarCore = 0.0; + float R2 = PS_SQR(radius); + float minSN2 = PS_SQR(minSN); + float rsigma2 = 0.5 / PS_SQR(sigma); // a note about coordinates: coordinates of objects throughout psphot refer to the primary @@ -109,35 +109,48 @@ // 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; - psF32 YY = 0.0; - psF32 XXX = 0.0; - psF32 XXY = 0.0; - psF32 XYY = 0.0; - psF32 YYY = 0.0; - psF32 XXXX = 0.0; - psF32 XXXY = 0.0; - psF32 XXYY = 0.0; - psF32 XYYY = 0.0; - psF32 YYYY = 0.0; + float RFW = 0.0; + float RHW = 0.0; + + float RF = 0.0; + float RH = 0.0; + float RS = 0.0; + float XX = 0.0; + float XY = 0.0; + float YY = 0.0; + float XXX = 0.0; + float XXY = 0.0; + float XYY = 0.0; + float YYY = 0.0; + float XXXX = 0.0; + float XXXY = 0.0; + float XXYY = 0.0; + float XYYY = 0.0; + float YYYY = 0.0; Sum = 0.0; // the second pass may include slightly different pixels, re-determine Sum + + // 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; + float Xo = source->moments->Mx; + float Yo = source->moments->My; // 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 = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage - psF32 yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage + float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage + float 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; + float yDiff = row - yCM; if (fabs(yDiff) > radius) continue; - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->variance->data.F32[row]; + float *vPix = source->pixels->data.F32[row]; + float *vWgt = source->variance->data.F32[row]; + psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; + // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { @@ -151,14 +164,14 @@ if (isnan(*vPix)) continue; - psF32 xDiff = col - xCM; + float xDiff = col - xCM; if (fabs(xDiff) > radius) continue; // radius is just a function of (xDiff, yDiff) - psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); if (r2 > R2) continue; - psF32 fDiff = *vPix - sky; - psF32 pDiff = fDiff; - psF32 wDiff = *vWgt; + float fDiff = *vPix - sky; + float pDiff = fDiff; + float wDiff = *vWgt; // skip pixels below specified significance level. this is allowed, but should be @@ -171,7 +184,7 @@ // weighting over weights the sky for faint sources if (sigma > 0.0) { - psF32 z = r2 * rsigma2; + float z = r2 * rsigma2; assert (z >= 0.0); - psF32 weight = exp(-z); + float weight = exp(-z); wDiff *= weight; @@ -182,30 +195,36 @@ // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?) - psF32 r = sqrt(r2); - psF32 rf = r * fDiff; - psF32 rh = sqrt(r) * fDiff; - psF32 rs = fDiff; - - psF32 x = xDiff * pDiff; - psF32 y = yDiff * pDiff; - - psF32 xx = xDiff * x; - psF32 xy = xDiff * y; - psF32 yy = yDiff * y; - - psF32 xxx = xDiff * xx / r; - psF32 xxy = xDiff * xy / r; - psF32 xyy = xDiff * yy / r; - psF32 yyy = yDiff * yy / r; - - psF32 xxxx = xDiff * xxx / r2; - psF32 xxxy = xDiff * xxy / r2; - psF32 xxyy = xDiff * xyy / r2; - psF32 xyyy = xDiff * yyy / r2; - psF32 yyyy = yDiff * yyy / r2; + float r = sqrt(r2); + float rf = r * fDiff; + float rh = sqrt(r) * fDiff; + float rs = fDiff; + + float rfw = r * pDiff; + float rhw = sqrt(r) * pDiff; + + float x = xDiff * pDiff; + float y = yDiff * pDiff; + + float xx = xDiff * x; + float xy = xDiff * y; + float yy = yDiff * y; + + float xxx = xDiff * xx / r; + float xxy = xDiff * xy / r; + float xyy = xDiff * yy / r; + float yyy = yDiff * yy / r; + + float xxxx = xDiff * xxx / r2; + float xxxy = xDiff * xxy / r2; + float xxyy = xDiff * xyy / r2; + float xyyy = xDiff * yyy / r2; + float yyyy = yDiff * yyy / r2; RF += rf; RH += rh; RS += rs; + + RFW += rfw; + RHW += rhw; XX += xx; @@ -244,8 +263,12 @@ source->moments->Myyyy = YYYY/Sum; - // Calculate the Kron magnitude (make this block optional?) - float radKinner = 1.0*source->moments->Mrf; - float radKron = 2.5*source->moments->Mrf; - float radKouter = 4.0*source->moments->Mrf; + // if Mrf (first radial moment) is very small, we are getting into low-significance + // territory. saturate at minKronRadius. conversely, if Mrf is > radius, we are clearly + // making an error. saturate at radius. + float kronRefRadius = MIN(radius, MAX(minKronRadius, source->moments->Mrf)); + + float radKinner = 1.0*kronRefRadius; + float radKron = 2.5*kronRefRadius; + float radKouter = 4.0*kronRefRadius; int nKronPix = 0; @@ -259,10 +282,12 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { - psF32 yDiff = row - yCM; + float yDiff = row - yCM; if (fabs(yDiff) > radKouter) continue; - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->variance->data.F32[row]; + float *vPix = source->pixels->data.F32[row]; + float *vWgt = source->variance->data.F32[row]; + psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; + // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { @@ -276,12 +301,12 @@ if (isnan(*vPix)) continue; - psF32 xDiff = col - xCM; + float xDiff = col - xCM; if (fabs(xDiff) > radKouter) continue; // radKron is just a function of (xDiff, yDiff) - psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); - - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; + float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + + float pDiff = *vPix - sky; + float wDiff = *vWgt; // skip pixels below specified significance level. this is allowed, but should be @@ -290,5 +315,5 @@ if (PS_SQR(pDiff) < minSN2*wDiff) continue; - psF32 r = sqrt(r2); + float r = sqrt(r2); if (r < radKron) { Sum += pDiff; @@ -335,5 +360,5 @@ } -bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal, float xGuess, float yGuess) { +bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { // First Pass: calculate the first moments (these are subtracted from the coordinates below) @@ -342,15 +367,15 @@ // .. etc - psF32 sky = 0.0; - - psF32 peakPixel = -PS_MAX_F32; + float sky = 0.0; + + float peakPixel = -PS_MAX_F32; psS32 numPixels = 0; - psF32 Sum = 0.0; - psF32 Var = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; - psF32 R2 = PS_SQR(radius); - psF32 minSN2 = PS_SQR(minSN); - psF32 rsigma2 = 0.5 / PS_SQR(sigma); + float Sum = 0.0; + float Var = 0.0; + float X1 = 0.0; + float Y1 = 0.0; + float R2 = PS_SQR(radius); + float minSN2 = PS_SQR(minSN); + float rsigma2 = 0.5 / PS_SQR(sigma); float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage @@ -358,5 +383,5 @@ // we are guaranteed to have a valid pixel and variance at this location (right? right?) - // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); + // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); @@ -370,10 +395,12 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { - psF32 yDiff = row + 0.5 - yPeak; + float yDiff = row + 0.5 - yPeak; if (fabs(yDiff) > radius) continue; - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->variance->data.F32[row]; + float *vPix = source->pixels->data.F32[row]; + float *vWgt = source->variance->data.F32[row]; + psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; + // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { @@ -387,13 +414,13 @@ if (isnan(*vPix)) continue; - psF32 xDiff = col + 0.5 - xPeak; + float xDiff = col + 0.5 - xPeak; if (fabs(xDiff) > radius) continue; // radius is just a function of (xDiff, yDiff) - psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); if (r2 > R2) continue; - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; + float pDiff = *vPix - sky; + float wDiff = *vWgt; // skip pixels below specified significance level. for a PSFs, this @@ -408,7 +435,7 @@ // weighting over weights the sky for faint sources if (sigma > 0.0) { - psF32 z = r2*rsigma2; + float z = r2*rsigma2; assert (z >= 0.0); - psF32 weight = exp(-z); + float weight = exp(-z); wDiff *= weight; @@ -419,6 +446,6 @@ Sum += pDiff; - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; + float xWght = xDiff * pDiff; + float yWght = yDiff * pDiff; X1 += xWght; @@ -477,2 +504,42 @@ return true; } + +float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) { + + psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32); + + for (int i = 0; i < sources->n; i++) { + pmSource *src = sources->data[i]; // Source of interest + if (!src || !src->moments) { + continue; + } + + if (src->mode & PM_SOURCE_MODE_BLEND) { + continue; + } + + if (!src->moments->nPixels) continue; + + if (src->moments->SN < PSF_SN_LIM) continue; + + // XXX put in Mxx,Myy cut based on clump location + + psVectorAppend(radii, src->moments->Mrf); + } + + // find the peak in this image + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + if (!psVectorStats (stats, radii, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n"); + psFree(stats); + return NAN; + } + + float minRadius = stats->sampleMedian; + + psFree(radii); + psFree(stats); + return minRadius; +} +