Changeset 25027 for branches/pap/psModules/src/objects/pmSourceMoments.c

Timestamp:

Aug 7, 2009, 4:08:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects/pmSourceMoments.c (modified) (8 diffs)

branches/pap/psModules/src/objects/pmSourceMoments.c

@@ -57 +57 @@
 # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
 
-bool pmSourceMoments(pmSource *source, psF32 radius)
+bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN)
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -84 +84 @@
     psF32 Y1 = 0.0;
     psF32 R2 = PS_SQR(radius);
+    psF32 minSN2 = PS_SQR(minSN);
+    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -97 +99 @@
 
     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];
@@ -113 +118 @@
 
             psF32 xDiff = col - xPeak;
-            psF32 yDiff = row - yPeak;
+            if (fabs(xDiff) > radius) continue;
 
             // radius is just a function of (xDiff, yDiff)
@@ -121 +126 @@
             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;
+            // skip pixels below specified significance level
+            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
+            if (pDiff < 0) continue;
+
+            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;
+            } 
+
+            Var += wDiff;
+            Sum += pDiff;
+
+            psF32 xWght = xDiff * pDiff;
+            psF32 yWght = yDiff * pDiff;
+
+            X1  += xWght;
+            Y1  += yWght;
+
+            // fprintf (stderr, " : %6.1f %6.1f  %8.1f %8.1f\n", X1, Y1, Sum, Var);
 
             peakPixel = PS_MAX (*vPix, peakPixel);
@@ -188 +211 @@
     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];
         psF32 *vWgt = source->variance->data.F32[row];
@@ -203 +229 @@
 
             psF32 xDiff = col - xCM;
-            psF32 yDiff = row - yCM;
+            if (fabs(xDiff) > radius) continue;
 
             // radius is just a function of (xDiff, yDiff)
@@ -215 +241 @@
             // XXX EAM : should this limit be user-defined?
 
-            if (PS_SQR(pDiff) < wDiff) continue;
+            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
             if (pDiff < 0) continue;
+
+            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;
+            } 
 
             Sum += pDiff;