Changeset 29546 for trunk/psModules/src/objects/pmSourceMoments.c

Timestamp:

Oct 25, 2010, 2:54:22 PM (16 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-20100823

File:

• trunk/psModules/src/objects/pmSourceMoments.c (modified) (11 diffs)

trunk/psModules/src/objects/pmSourceMoments.c

@@ -64 +64 @@
 void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
 
+// 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)
 {
@@ -71 +73 @@
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
-    // 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:
+    // this function assumes the sky has been well-subtracted for the image
     psF32 sky = 0.0;
+
     if (source->moments == NULL) {
       source->moments = pmMomentsAlloc();
     }
-    // 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)
-    // Sum = SUM (z - sky)
-    // X1  = SUM (x - xc)*(z - sky)
-    // .. etc
-
-    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 R2 = PS_SQR(radius);
     psF32 minSN2 = PS_SQR(minSN);
@@ -109 +94 @@
     // (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++) {
-
-        psF32 yDiff = row - yPeak;
-        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 & maskVal) {
-                    vMsk++;
-                    continue;
-                }
-                vMsk++;
-            }
-            if (isnan(*vPix)) continue;
-
-            psF32 xDiff = col - xPeak;
-            if (fabs(xDiff) > radius) continue;
-
-            // radius is just a function of (xDiff, yDiff)
-            if (!VALID_RADIUS(xDiff, yDiff, R2)) continue;
-
-            psF32 pDiff = *vPix - sky;
-            psF32 wDiff = *vWgt;
-
-            // 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
-            // weighting over weights the sky for faint sources
-            if (sigma > 0.0) {
-                // 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;
-            } 
-
-            Var += wDiff;
-            Sum += pDiff;
-
-            psF32 xWght = xDiff * pDiff;
-            psF32 yWght = yDiff * pDiff;
-
-            X1  += xWght;
-            Y1  += yWght;
-
-            peakPixel = PS_MAX (*vPix, peakPixel);
-            numPixels++;
-        }
-    }
-
-    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
-    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
-
-    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
-    if ((numPixels < minPixels) || (Sum <= 0)) {
-        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
-        return (false);
-    }
-
-    // calculate the first moment.
-    float Mx = X1/Sum;
-    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", 5, "sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", sky, Mx, My, Sum, X1, Y1, numPixels);
-
-    // add back offset of peak in primary image
-    // 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;
-    source->moments->SN  = Sum / sqrt(Var);
-    source->moments->Peak = peakPixel;
-    source->moments->nPixels = numPixels;
+    // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to 
+    // get an unbiased (but probably noisy) centroid
+    if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
+        return false;
+    }
+    // second pass applies the Gaussian window and uses the centroid from the first pass
+    if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) {
+        return false;
+    }
 
     // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
@@ -234 +125 @@
     Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
 
-    // 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
+    // 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
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -263 +154 @@
             // radius is just a function of (xDiff, yDiff)
             psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-            psF32 r  = sqrt(r2);
-            if (r > radius) continue;
+            if (r2 > R2) continue;
 
             psF32 fDiff = *vPix - sky;
@@ -274 +164 @@
             // stars.
             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
-            // if (pDiff < 0) continue;
+            if ((minSN > 0.0) && (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) {
-                // XXX a lot of extra flops; can we do pre-calculate?
-                // XXX we were re-calculating r2 (maybe the compiler caught this?)
-                // psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
-                psF32 z  = r2 * rsigma2;
+                psF32 z = r2 * rsigma2;
                 assert (z >= 0.0);
                 psF32 weight  = exp(-z);
@@ -292 +179 @@
             Sum += pDiff;
 
-# if (1)
-# if (0)
-            if (fDiff < 0) continue;
-# endif
+            // 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;
-# else
-            psF32 rf = r * pDiff;
-            psF32 rh = sqrt(r) * pDiff;
-            psF32 rs = pDiff;
-# endif
 
             psF32 x = xDiff * pDiff;
@@ -363 +243 @@
 
     // 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;
@@ -397 +276 @@
             // radKron is just a function of (xDiff, yDiff)
             psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-            psF32 r  = sqrt(r2);
 
             psF32 pDiff = *vPix - sky;
@@ -407 +285 @@
             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
 
+            psF32 r  = sqrt(r2);
             if (r < radKron) {
                 Sum += pDiff;
@@ -434 +313 @@
 
     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);
+             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, source->moments->nPixels);
 
     return(true);
 }
+
+bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
+
+    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
+    // Sum = SUM (z - sky)
+    // X1  = SUM (x - xc)*(z - sky)
+    // .. etc
+
+    psF32 sky = 0.0;
+
+    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 R2 = PS_SQR(radius);
+    psF32 minSN2 = PS_SQR(minSN);
+    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
+
+    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
+    float yPeak = yGuess - 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");
+
+    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
+    // not depend on the fractional pixel location of the source.  However, the aperture
+    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
+    // position of the expected centroid
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+        psF32 yDiff = row + 0.5 - yPeak;
+        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 & maskVal) {
+                    vMsk++;
+                    continue;
+                }
+                vMsk++;
+            }
+            if (isnan(*vPix)) continue;
+
+            psF32 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);
+            if (r2 > R2) continue;
+
+            psF32 pDiff = *vPix - sky;
+            psF32 wDiff = *vWgt;
+
+            // 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
+            // weighting over weights the sky for faint sources
+            if (sigma > 0.0) {
+                psF32 z  = r2*rsigma2;
+                assert (z >= 0.0);
+                psF32 weight  = exp(-z);
+
+                wDiff *= weight;
+                pDiff *= weight;
+            } 
+
+            Var += wDiff;
+            Sum += pDiff;
+
+            psF32 xWght = xDiff * pDiff;
+            psF32 yWght = yDiff * pDiff;
+
+            X1  += xWght;
+            Y1  += yWght;
+
+            peakPixel = PS_MAX (*vPix, peakPixel);
+            numPixels++;
+        }
+    }
+
+    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
+    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
+
+    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
+    if ((numPixels < minPixels) || (Sum <= 0)) {
+        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
+        return (false);
+    }
+
+    // calculate the first moment.
+    float Mx = X1/Sum;
+    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", 5, "sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", sky, Mx, My, Sum, X1, Y1, numPixels);
+
+    // add back offset of peak in primary image
+    // 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
+
+    // we only update the centroid if the position is not supplied from elsewhere
+    bool skipCentroid = false;
+    skipCentroid |= (source->mode  & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions
+    skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions
+
+    if (skipCentroid) {
+        source->moments->Mx = source->peak->xf;
+        source->moments->My = source->peak->yf;
+    } else {
+        source->moments->Mx = Mx + xGuess;
+        source->moments->My = My + yGuess;
+    }
+
+    source->moments->Sum = Sum;
+    source->moments->SN  = Sum / sqrt(Var);
+    source->moments->Peak = peakPixel;
+    source->moments->nPixels = numPixels;
+
+    return true;
+}