Index: trunk/psphot/src/psphotMakeResiduals.c
===================================================================
--- trunk/psphot/src/psphotMakeResiduals.c	(revision 20307)
+++ trunk/psphot/src/psphotMakeResiduals.c	(revision 20451)
@@ -125,4 +125,5 @@
     }
     pmResiduals *resid = pmResidualsAlloc (xSize, ySize, xBin, yBin);
+    psImageInit (resid->mask, 0);
 
     // x(resid) = (x(image) - Xo)*xBin + xCenter
@@ -195,4 +196,5 @@
 
             // mark input pixels which are more than N sigma from the median
+	    int nKeep = 0;
             for (int i = 0; i < fluxes->n; i++) {
                 float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
@@ -204,4 +206,5 @@
                     fmasks->data.U8[i] = clippedMask;
                 }
+		if (!fmasks->data.U8[i]) nKeep++;
             }
 
@@ -216,7 +219,9 @@
                 //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
 
-                if (resid->Ro->data.F32[oy][ox] < pixelSN*fluxStats->sampleStdev) {
+                if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*fluxStats->sampleStdev/sqrt(nKeep)) {
                   resid->mask->data.U8[oy][ox] = 1;
                 }
+
+		// fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n", ox, oy, resid->Ro->data.F32[oy][ox], fluxStats->sampleStdev, fluxStats->sampleStdev/sqrt(nKeep), nKeep, resid->mask->data.U8[oy][ox]);
 
             } else {
@@ -255,5 +260,7 @@
 
                 float dRo = sqrt(A->data.F32[0][0]);
-                if (resid->Ro->data.F32[oy][ox] < pixelSN*dRo) {
+		// fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n", ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.U8[oy][ox]);
+
+                if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) {
                   resid->mask->data.U8[oy][ox] = 1;
                 }