Changeset 41957

Timestamp:

Dec 3, 2021, 8:02:04 AM (5 years ago)

Author:

eugene

Message:

calculate alternate variance value if number of measurements > 8

File:

• branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c (modified) (4 diffs)

branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c

@@ -1715 +1715 @@
             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 +1735 @@
             }
             float mean = sum / (float) Npt;
+            float varValue = varSum / (float) (nGoodClip*nGoodClip);
 
             // alternative: calculate the stdev of the pixel values
@@ -1742 +1743 @@
             // }
             // 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 +1768 @@
 
             // 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;