Changeset 29473

Timestamp:

Oct 17, 2010, 12:31:48 PM (16 years ago)

Author:

eugene

Message:

the measured centroid, with a Gaussian window, is biased by the error in the guess position; do an initial iteration on the centroid using an un-windowed measurement, then refine with the gaussian window

Location:

branches/eam_branches/ipp-20100823/psModules/src/objects

Files:

• pmSource.h (modified) (1 diff)
• pmSourceMoments.c (modified) (7 diffs)

branches/eam_branches/ipp-20100823/psModules/src/objects/pmSource.h

@@ -249 +249 @@
   psF32 sigma, 
   psF32 minSN, 
-  psImageMaskType maskVal);
+  psImageMaskType maskVal, 
+  float xGuess, float yGuess);
 
 pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source);

branches/eam_branches/ipp-20100823/psModules/src/objects/pmSourceMoments.c

@@ -64 +64 @@
 void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
 
-bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal);
-
 // if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
 
@@ -96 +94 @@
     // (int) so they can be used in the image index below.
 
-    if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal)) {
+    // 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;
     }
@@ -314 +318 @@
 }
 
-bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal) { 
+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)
@@ -333 +337 @@
     psF32 rsigma2 = 0.5 / PS_SQR(sigma);
 
-    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
+    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?)
@@ -344 +346 @@
     // 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 - yPeak;
+        psF32 yDiff = row + 0.5 - yPeak;
         if (fabs(yDiff) > radius) continue;
 
@@ -363 +370 @@
             if (isnan(*vPix)) continue;
 
-            psF32 xDiff = col - xPeak;
+            psF32 xDiff = col + 0.5 - xPeak;
             if (fabs(xDiff) > radius) continue;
 
@@ -439 +446 @@
         source->moments->My = source->peak->yf;
     } else {
-        source->moments->Mx = Mx + xOff + 0.5;
-        source->moments->My = My + yOff + 0.5;
+        source->moments->Mx = Mx + xGuess;
+        source->moments->My = My + yGuess;
     }