Changeset 41945

Timestamp:

Dec 1, 2021, 9:27:09 AM (5 years ago)

Author:

eugene

Message:

add coarse outlier rejection

File:

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

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

@@ -1538 +1538 @@
 #define SORT_VALUES(NVALUES) { PSSORT(NVALUES, SORT_VV_COMPARE, SORT_VV_SWAP, F32); }
 
+#define ESCAPE                                                          \
+  combined->image->data.F32[y][x] = NAN;                                \
+  combined->variance->data.F32[y][x] = NAN;                             \
+  combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = blankMaskBits;   \
+  if (expmaps) {                                                        \
+    expmaps->image->data.F32[y][x] = 0.0;                               \
+    expmaps->variance->data.F32[y][x] = 0.0;                            \
+    expmaps->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;              \
+  }                                                                     \
+  continue;
+
 bool pmStackCombineByPercentile(
     pmReadout *combined,                // output stacked readout
@@ -1631 +1642 @@
                 nGood ++;
             }
+# define MIN_GOOD_PERCENTILE 3
+
+            if (nGood < MIN_GOOD_PERCENTILE) ESCAPE;
+            
             pixelData->n = nGood;
+
+            // sort pixelData, pixelVariance, expTime (if it exists)
+            SORT_VALUES (pixelData->n);
+
+            // we now have a sorted vector of values.  We can make a very coarse outlier
+            // cut based on the median and interquartile range.  This will let us reject
+            // some extreme outliers that bias the signal otherwise.
+
+            int Nlo = 0;     
+            int Nhi = nGood; 
+
+            if (nGood >= 5) {
+              int midPoint = nGood / 2;
+              float rawMedian = (nGood % 2) ? pixelData->data.F32[midPoint] : 0.5*(pixelData->data.F32[midPoint] + pixelData->data.F32[midPoint-1]);
+              
+              // XXX measure interquartile range
+              int P25 = 0.25*nGood;
+              int P75 = 0.75*nGood;
+              float rawSigma = 0.74*(pixelData->data.F32[P75] - pixelData->data.F32[P25]);
+              float minThresh = rawMedian - 7.0*rawSigma;
+              float maxThresh = rawMedian + 7.0*rawSigma;
+
+              // find the entries which are in the range
+              // these should be safe since minThresh and maxThresh are guaranteed to contain the median
+              while (pixelData->data.F32[Nlo    ] < minThresh) { Nlo ++; }
+              while (pixelData->data.F32[Nhi - 1] > maxThresh) { Nhi --; }
+            }
+
+            // if we did not clip above (either nGood < 5 or no rejection), Nlo will be 0, Nhi will be nGood
+            // In either case, Nlo is the offset of the first unclipped point.
+            int nGoodClip = Nhi - Nlo;
 
             // rather than define a min and max value,
@@ -1645 +1691 @@
 # define MIN_GOOD_PERCENTILE 3
 # if (AT_LEAST)
-            int Ns = MIN(MAX(1, 0.5*rejectFraction * nGood), nGood);
-            int Ne = nGood - Ns;
+            int Ns = MIN(MAX(1, 0.5*rejectFraction * nGoodClip), nGoodClip) + Nlo;
+            int Ne = nGood - Ns + Nlo;
             int Npt = Ne - Ns;
 # else
-            int Ns = MIN(MAX(0, 0.5*rejectFraction * nGood), nGood);
-            int Ne = nGood - Ns;
+            int Ns = MIN(MAX(0, 0.5*rejectFraction * nGood), nGood) + Nlo;
+            int Ne = nGood - Ns + Nlo;
             int Npt = Ne - Ns;
 # endif
-            if ((Npt < 1) || (nGood < MIN_GOOD_PERCENTILE)) {
-                combined->image->data.F32[y][x] = NAN;
-                combined->variance->data.F32[y][x] = NAN;
-                combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = blankMaskBits; // probably not needed since it was set above.
-                if (expmaps) {
-                    expmaps->image->data.F32[y][x] = 0.0; // XXX should this value be NAN or 0?
-                    expmaps->variance->data.F32[y][x] = 0.0;
-                    expmaps->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
-                }
-                continue;
-            }
-
-            // sort pixelData, pixelVariance, expTime (if it exists)
-            SORT_VALUES (pixelData->n);
-
+            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
@@ -1705 +1739 @@
             // this coefficient varies between 1.4 (for pure median) and 1.05 for 68% range.
 
-            float varValue = varSum / (float) (nGood*nGood);
+            float varValue = varSum / (float) (nGoodClip*nGoodClip);
 
             combined->image->data.F32[y][x] = mean;
@@ -1719 +1753 @@
             if (expTime) { 
                 float sum = 0.0;
-                for (int n = 0; n < nGood; n++) {
+                for (int n = Nlo; n < Nhi; n++) {
                     sum += expTime->data.F32[n];
                 }
                 expmaps->image->data.F32[y][x] = sum;
             }
-            if (expmaps) { expmaps->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = nGood; }
+            if (expmaps) { expmaps->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = nGoodClip; }
         }
     }