Index: /branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c
===================================================================
--- /branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c	(revision 41923)
+++ /branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c	(revision 41924)
@@ -400,4 +400,6 @@
     //     variance_combination = variance_individual / N
     // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N).
+    // NOTE: in 2012.07.13, the variance calculation was changed without justification to the variance
+    // appropriate to a weighted mean, while the pixel mean was kept as the average unweighted by pixel variance
 
     float sumValueWeight = 0.0;         // Sum of the value multiplied by the weight
@@ -1686,5 +1688,16 @@
 	    // }
 	    // variance on the mean (stdev / sqrt(N))^2
-	    float varValue = varSum / (Npt*Npt);
+
+	    // Note: since we are calculating the average of a subset of a sorted
+	    // list of values, the denominator should not be the number of measurements
+	    // in the calculation above (Npt): in the extreme case of a median, we would
+	    // have a single value (Npt = 1), but the variance of a median is only ~1.4 x
+	    // the variance of the average / sqrt(N).  We should use nGood (the total number
+	    // of values in the sorted list), but the variance should be scaled by a factor
+	    // which depends on the fraction of values included.  
+
+	    // this coefficient varies between 1.4 (for pure median) and 1.05 for 68% range.
+
+	    float varValue = varSum / (float) (nGood*nGood);
 
 	    combined->image->data.F32[y][x] = mean;