Index: branches/eam_branches/ipp-dev-20210817/psModules/src/imcombine/pmStack.c
===================================================================
--- branches/eam_branches/ipp-dev-20210817/psModules/src/imcombine/pmStack.c	(revision 41851)
+++ branches/eam_branches/ipp-dev-20210817/psModules/src/imcombine/pmStack.c	(revision 41852)
@@ -756,4 +756,7 @@
 	// *goodMask &= mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn]; // save the mask bits still used
 	// check for set bits and increment counter as appropriate
+	// count the number of times a given mask bit is set in the input pixels.
+	// NOTE: since we have explicitly skipped the pixels with any bad bits, these are only
+	// the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs'
 	{ 
 	    psImageMaskType value = 0x0001;
@@ -1550,4 +1553,5 @@
 		pmReadout *ro = data->readout;  // Readout of interest
 		psImage *image = ro->image;
+		psImage *mask = ro->mask;
 
 		int xIn = x - data->readout->col0;
@@ -1557,4 +1561,22 @@
 		if (fabs(image->data.F32[yIn][xIn]) > 1e5) continue;
 		// XXX need to test the input mask as well here
+		if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits);
+
+		// XXX save the count of suspect bits
+
+# if (0)
+		// count the number of times a given mask bit is set in the input pixels.
+		// NOTE: since we have explicitly skipped the pixels with any bad bits, these are only
+		// the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs'
+		if (1) { 
+		  psImageMaskType value = 0x0001;
+		  for (int nbit = 0; nbit < 16; nbit ++) {
+		    if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & value) {
+		      // XXX nGoodBits[nbit] ++;
+		    }
+		    value <<= 1;
+		  }
+		}
+# endif
 
 		pixelData->data.F32[nGood] = image->data.F32[yIn][xIn];
@@ -1571,4 +1593,18 @@
 	    }
 
+# if (0)
+# define SUSPECT_FRACTION 0.65
+	    // set the mask bits if nGoodBits[i] > f*numGood
+	    if (0) {
+	      psImageMaskType value = 0x0001;
+	      for (int nbit = 0; nbit < 16; nbit ++) {
+		if (nGoodBits[nbit] > SUSPECT_FRACTION*numGood) {
+		  *goodMask |= value;
+		}
+		value <<= 1;
+	      }
+	    }
+# endif
+
 	    int Ns = MIN(MAX(0, minRange * nGood),nGood); // e.g., 0.1 * 50 = 5
 	    int Ne = MIN(MAX(0, maxRange * nGood),nGood); // e.g., 0.9 * 50 = 45 
@@ -1587,4 +1623,7 @@
 	    // variance on the mean (stdev / sqrt(N))^2
 	    float varValue = varSum / (Npt - 1) / Npt;
+
+	    // XXX this is probably not the correct variance to save
+	    // the predicted variance should be the 1 / sum (1/var)
 
 	    combined->image->data.F32[y][x] = mean;