Index: trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- trunk/psModules/src/objects/pmSourceMoments.c (revision 24577) +++ trunk/psModules/src/objects/pmSourceMoments.c (revision 25754) @@ -36,6 +36,4 @@ #include "pmSource.h" -bool pmSourceMoments_Old(pmSource *source, psF32 radius); - /****************************************************************************** pmSourceMoments(source, radius): this function takes a subImage defined in the @@ -50,8 +48,7 @@ pmSource->mask (optional) -XXX: The peak calculations are done in image coords, not subImage coords. - -this version clips input pixels on S/N -XXX EAM : this version returns false for several reasons +The peak calculations are done in image coords, not subImage coords. + +this version optionally clips input pixels on S/N *****************************************************************************/ # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) @@ -64,5 +61,5 @@ PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); - // XXX supply sky in a different way? + // use sky from moments if defined, 0.0 otherwise psF32 sky = 0.0; if (source->moments == NULL) { @@ -91,10 +88,18 @@ // 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 - - psF32 xOff = source->peak->x; - psF32 yOff = source->peak->y; - psF32 xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage - psF32 yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage + // 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. + + int xOff = source->peak->x; + int yOff = source->peak->y; + int xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage + int yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage + + // 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]); + // 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"); + // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel"); for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -126,21 +131,20 @@ psF32 wDiff = *vWgt; - // skip pixels below specified significance level + // 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; + // Apply a Gaussian window function. Be careful with the window function. S/N + // weighting over weights the sky for faint sources if (sigma > 0.0) { - // apply a pseudo-gaussian weight - - // XXX a lot of extra flops; can we do pre-calculate? - psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; - assert (z >= 0.0); - psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0)); - psF32 weight = 1.0 / t; - - // fprintf (stderr, "%6.1f %6.1f %8.1f %8.1f %5.3f ", xDiff, yDiff, pDiff, wDiff, weight); - - wDiff *= weight; - pDiff *= weight; + // XXX a lot of extra flops; can we pre-calculate? + psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + assert (z >= 0.0); + psF32 weight = exp(-z); + + wDiff *= weight; + pDiff *= weight; } @@ -154,9 +158,7 @@ Y1 += yWght; - // fprintf (stderr, " : %6.1f %6.1f %8.1f %8.1f\n", X1, Y1, Sum, Var); - - peakPixel = PS_MAX (*vPix, peakPixel); - numPixels++; - } + peakPixel = PS_MAX (*vPix, peakPixel); + numPixels++; + } } @@ -164,6 +166,6 @@ // XXX EAM - the limit is a bit arbitrary. make it user defined? if ((numPixels < 0.75*R2) || (Sum <= 0)) { - psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum); - return (false); + psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum); + return (false); } @@ -172,6 +174,6 @@ float My = Y1/Sum; if ((fabs(Mx) > radius) || (fabs(My) > radius)) { - psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y); - return (false); + psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y); + return (false); } @@ -179,6 +181,9 @@ // add back offset of peak in primary image - source->moments->Mx = Mx + xOff; - source->moments->My = My + yOff; + // also offset from pixel index to pixel coordinate + // (the calculation above uses pixel index instead of coordinate) + // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords + source->moments->Mx = Mx + xOff + 0.5; + source->moments->My = My + yOff + 0.5; source->moments->Sum = Sum; @@ -205,5 +210,6 @@ Sum = 0.0; // the second pass may include slightly different pixels, re-determine Sum - // center of mass in subimage + // center of mass in subimage. Note: the calculation below uses pixel index, so we do not + // correct xCM, yCM to pixel coordinate here. psF32 xCM = Mx + xPeak; // coord of peak in subimage psF32 yCM = My + yPeak; // coord of peak in subimage @@ -211,87 +217,85 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { - psF32 yDiff = row - yCM; + psF32 yDiff = row - yCM; if (fabs(yDiff) > radius) 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) { - vMsk++; - continue; - } - vMsk++; - } - if (isnan(*vPix)) continue; - - psF32 xDiff = col - xCM; + 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) { + 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); - psF32 r = sqrt(r2); - if (r > radius) continue; - - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; - - // XXX EAM : should this limit be user-defined? - - if (PS_SQR(pDiff) < minSN2*wDiff) continue; - if (pDiff < 0) continue; - + // radius is just a function of (xDiff, yDiff) + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + psF32 r = sqrt(r2); + if (r > radius) continue; + + 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 (pDiff < 0) continue; + + // Apply a Gaussian window function. Be careful with the window function. S/N + // weighting over weights the sky for faint sources if (sigma > 0.0) { - // apply a pseudo-gaussian weight - - // XXX a lot of extra flops; can we do pre-calculate? - psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; - assert (z >= 0.0); - psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0)); - psF32 weight = 1.0 / t; - - // fprintf (stderr, "%6.1f %6.1f %8.1f %8.1f %5.3f ", xDiff, yDiff, pDiff, wDiff, weight); - - wDiff *= weight; - pDiff *= weight; + // XXX a lot of extra flops; can we do pre-calculate? + psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + assert (z >= 0.0); + psF32 weight = exp(-z); + + wDiff *= weight; + pDiff *= weight; } - Sum += pDiff; - - 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; - - XX += xx; - XY += xy; - YY += yy; - - XXX += xxx; - XXY += xxy; - XYY += xyy; - YYY += yyy; - - XXXX += xxxx; - XXXY += xxxy; - XXYY += xxyy; - XYYY += xyyy; - YYYY += yyyy; - } + Sum += pDiff; + + 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; + + XX += xx; + XY += xy; + YY += yy; + + XXX += xxx; + XXY += xxy; + XYY += xyy; + YYY += yyy; + + XXXX += xxxx; + XXXY += xxxy; + XXYY += xxyy; + XYYY += xyyy; + YYYY += yyyy; + } } @@ -312,151 +316,18 @@ if (source->moments->Mxx < 0) { - fprintf (stderr, "error: neg second moment??\n"); + fprintf (stderr, "error: neg second moment??\n"); } if (source->moments->Myy < 0) { - fprintf (stderr, "error: neg second moment??\n"); + 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", - source->moments->Mxx, source->moments->Mxy, source->moments->Myy, - source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, - source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy); + source->moments->Mxx, source->moments->Mxy, source->moments->Myy, + source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, + source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy); psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f My: %f Sum: %f Mxx: %f Mxy: %f Myy: %f sky: %f Npix: %d\n", - source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx, source->moments->My, Sum, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, sky, numPixels); - - if (source->moments->Mxx < 0) { - fprintf (stderr, "error: neg second moment??\n"); - } - if (source->moments->Myy < 0) { - fprintf (stderr, "error: neg second moment??\n"); - } - - // XXX TEST: - // pmSourceMoments_Old (source, radius); + source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx, source->moments->My, Sum, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, sky, numPixels); + return(true); } - - -bool pmSourceMoments_Old(pmSource *source, psF32 radius) -{ - 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 sky = 0.0; - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } else { - sky = source->moments->Sky; - } - - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // X2 = SUM (x - xc)^2 * (z - sky) - // XY = SUM (x - xc)*(y - yc)*(z - sky) - psF32 peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; - psF32 Sum = 0.0; - psF32 Var = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; - psF32 X2 = 0.0; - psF32 Y2 = 0.0; - psF32 XY = 0.0; - psF32 x = 0; - psF32 y = 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 - - psF32 xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage - psF32 yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - - 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) { - vMsk++; - continue; - } - vMsk++; - } - if (isnan(*vPix)) continue; - - psF32 xDiff = col - xPeak; - psF32 yDiff = row - yPeak; - - // radius is just a function of (xDiff, yDiff) - if (!VALID_RADIUS(xDiff, yDiff, R2)) continue; - - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; - - // XXX EAM : should this limit be user-defined? - if (PS_SQR(pDiff) < wDiff) continue; - - Var += wDiff; - Sum += pDiff; - - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; - - X1 += xWght; - Y1 += yWght; - - X2 += xDiff * xWght; - XY += xDiff * yWght; - Y2 += yDiff * yWght; - - peakPixel = PS_MAX (*vPix, peakPixel); - numPixels++; - } - } - - // if we have less than (1/4) of the possible pixels, force a retry - // XXX EAM - the limit is a bit arbitrary. make it user defined? - if ((numPixels < 0.75*R2) || (Sum <= 0)) { - psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum); - return (false); - } - - psTrace ("psModules.objects", 4, "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", - sky, Sum, X1, Y1, X2, Y2, XY, numPixels); - - // - // first moment X = X1/Sum + xc - // second moment X = sqrt (X2/Sum - (X1/Sum)^2) - // Sxy = XY / Sum - // - x = X1/Sum; - y = Y1/Sum; - if ((fabs(x) > radius) || (fabs(y) > radius)) { - psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", - source->peak->x, source->peak->y); - return (false); - } - -# if (PS_TRACE_ON) - float Sxx = PS_MAX(X2/Sum - PS_SQR(x), 0); - float Sxy = XY / Sum; - float Syy = PS_MAX(Y2/Sum - PS_SQR(y), 0); - - psTrace ("psModules.objects", 3, - "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", - sky, Sum, x, y, Sxx, Sxy, Syy); -# endif - - return(true); -} -