Index: /branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c
===================================================================
--- /branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c	(revision 41956)
+++ /branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c	(revision 41957)
@@ -1715,8 +1715,8 @@
 	    if ((Npt < 1) || (nGood < MIN_GOOD_PERCENTILE)) ESCAPE;
 
-	    // XXX note that this does not respect the clipping above
-	    // set the given (suspect) mask bit if nGoodBits[i] > f*nGood
-	    // in other words if more than 65% of the good inputs had one of
-	    // these bits set, then we should set that bit in the output mask
+	    // Set the given (suspect) mask bit if nGoodBits[i] > f*nGood in other words
+	    // if more than 65% of the good inputs had one of these bits set, then we
+	    // should set that bit in the output mask.  Note that this analysis counts the
+	    // mask bits of pixels rejected by the clipping above.
 	    psImageMaskType value = 0x0001;
 	    psImageMaskType outputMask = 0x0000;
@@ -1735,4 +1735,5 @@
 	    }
 	    float mean = sum / (float) Npt;
+	    float varValue = varSum / (float) (nGoodClip*nGoodClip);
 
 	    // alternative: calculate the stdev of the pixel values
@@ -1742,4 +1743,19 @@
 	    // }
 	    // variance on the mean (stdev / sqrt(N))^2
+
+	    // the reported variance values can be extremely high / wrong.
+	    // if we have enough measurements, let's just use the interquartile range
+	    // of the data to estimate the per-pixel variance.  NOTE: this is not valid
+	    // if the inputs have been significantly smoothed.  In that case we need
+	    // to include the covariance explicitly.  But this algorithm should be used
+	    // without convolution.
+	    // XXX How do we choose the cutoff here?
+	    if (nGoodClip >= 9) {
+	      // Measure interquartile range
+	      int P25 = 0.25*nGoodClip + Nlo;
+	      int P75 = 0.75*nGoodClip + Nlo;
+	      float rawSigma = 0.74*(pixelData->data.F32[P75] - pixelData->data.F32[P25]);
+	      varValue = PS_SQR(rawSigma) / (float) nGoodClip; // sigma_mean = sigma_meas / sqrt(Nmeas) -> var_mean = var_meas / Nmeas
+	    }
 
 	    // Note: since we are calculating the average of a subset of a sorted
@@ -1752,6 +1768,4 @@
 
 	    // this coefficient varies between 1.4 (for pure median) and 1.05 for 68% range.
-
-	    float varValue = varSum / (float) (nGoodClip*nGoodClip);
 
 	    combined->image->data.F32[y][x] = mean;