Index: trunk/psModules/src/objects/pmSourceMoments.c
===================================================================
--- trunk/psModules/src/objects/pmSourceMoments.c	(revision 21363)
+++ trunk/psModules/src/objects/pmSourceMoments.c	(revision 24577)
@@ -57,5 +57,5 @@
 # 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,4 +84,6 @@
     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,4 +99,7 @@
 
     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,5 +118,5 @@
 
             psF32 xDiff = col - xPeak;
-            psF32 yDiff = row - yPeak;
+	    if (fabs(xDiff) > radius) continue;
 
             // radius is just a function of (xDiff, yDiff)
@@ -121,15 +126,33 @@
             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,4 +211,7 @@
     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,5 +229,5 @@
 
             psF32 xDiff = col - xCM;
-            psF32 yDiff = row - yCM;
+	    if (fabs(xDiff) > radius) continue;
 
             // radius is just a function of (xDiff, yDiff)
@@ -215,6 +241,21 @@
             // 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;