Index: /branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c
===================================================================
--- /branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c	(revision 41944)
+++ /branches/eam_branches/ipp-20211108/psModules/src/imcombine/pmStack.c	(revision 41945)
@@ -1538,4 +1538,15 @@
 #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,5 +1642,40 @@
 		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,27 +1691,15 @@
 # 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,5 +1739,5 @@
 	    // 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,10 +1753,10 @@
 	    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; }
 	}
     }