Index: /branches/eam_branches/psModules.20211015/src/imcombine/pmStack.c
===================================================================
--- /branches/eam_branches/psModules.20211015/src/imcombine/pmStack.c	(revision 41867)
+++ /branches/eam_branches/psModules.20211015/src/imcombine/pmStack.c	(revision 41868)
@@ -1510,4 +1510,256 @@
 }
 
+# define MIN_GOOD_PERCENTILE 2
+bool pmStackCombineByPercentile(
+    pmReadout *combined,
+    psArray *stackData,
+    psF64 minRange,
+    psF64 maxRange,
+    psImageMaskType badMaskBits, 	// treat these bits as 'bad'
+    psImageMaskType suspectMaskBits,	// treat these bits as 'suspect'
+    psImageMaskType blankMaskBits       // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?)
+) {
+    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+    int xSize, ySize;                   // Size of the output image
+
+    // we need to copy the readouts to their own array for validation
+    psArray *stackReadouts = psArrayAlloc(stackData->n);
+    for (int i = 0; i < stackData->n; i++) {
+        pmStackData *data = stackData->data[i]; // Stack data for this input
+        pmReadout *ro = data->readout;  // Readout of interest
+        stackReadouts->data[i] = psMemIncrRefCounter(ro);
+    }
+
+    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, stackReadouts)) {
+	psError(psErrorCodeLast(), false, "Input stack is not valid.");
+	psFree(stackReadouts);
+	return false;
+    }
+    psFree(stackReadouts);
+
+    psVector *pixelData = psVectorAlloc(stackData->n, PS_TYPE_F32);
+    psVector *pixelMask = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK);
+
+    for (int y = minInputRows; y < maxInputRows; y++) {
+	for (int x = minInputCols; x < maxInputCols; x++) {
+
+	    int nGood = 0; 
+	    for (int i = 0; i < stackData->n; i++) {
+
+		pmStackData *data = stackData->data[i]; // Stack data for this input
+		pmReadout *ro = data->readout;  // Readout of interest
+		psImage *image = ro->image;
+		psImage *mask = ro->mask;
+
+		int xIn = x - data->readout->col0;
+		int yIn = y - data->readout->row0; // Coordinates on input readout
+
+		if (!isfinite(image->data.F32[yIn][xIn])) continue;
+		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) continue;
+
+		// 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];
+		nGood ++;
+	    }
+	    pixelData->n = nGood;
+	    psVectorSortInPlace (pixelData);
+
+	    if (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] = 1;
+		continue;
+	    }
+
+# 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 
+	    int Npt = Ne - Ns;
+
+	    float sum = 0.0;
+	    for (int n = Ns; n < Ne; n++) {
+		sum += pixelData->data.F32[n];
+	    }
+	    float mean = sum / (float) Npt;
+
+	    float varSum = 0.0;
+	    for (int n = Ns; n < Ne; n++) {
+		varSum += SQ(pixelData->data.F32[n] - mean);
+	    }
+	    // 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;
+	    combined->variance->data.F32[y][x] = varValue;
+	    combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+	}
+    }
+  
+    psFree(pixelData);
+    psFree(pixelMask);
+
+    return (true);
+}
+
+# define MIN_GOOD_PERCENTILE 2
+bool pmStackCombineByPercentile(
+    pmReadout *combined,
+    psArray *stackData,
+    psF64 minRange,
+    psF64 maxRange,
+    psImageMaskType badMaskBits, 	// treat these bits as 'bad'
+    psImageMaskType suspectMaskBits,	// treat these bits as 'suspect'
+    psImageMaskType blankMaskBits       // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?)
+) {
+    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+    int xSize, ySize;                   // Size of the output image
+
+    // we need to copy the readouts to their own array for validation
+    psArray *stackReadouts = psArrayAlloc(stackData->n);
+    for (int i = 0; i < stackData->n; i++) {
+        pmStackData *data = stackData->data[i]; // Stack data for this input
+        pmReadout *ro = data->readout;  // Readout of interest
+        stackReadouts->data[i] = psMemIncrRefCounter(ro);
+    }
+
+    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, stackReadouts)) {
+	psError(psErrorCodeLast(), false, "Input stack is not valid.");
+	psFree(stackReadouts);
+	return false;
+    }
+    psFree(stackReadouts);
+
+    psVector *pixelData = psVectorAlloc(stackData->n, PS_TYPE_F32);
+    psVector *pixelMask = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK);
+
+    for (int y = minInputRows; y < maxInputRows; y++) {
+	for (int x = minInputCols; x < maxInputCols; x++) {
+
+	    int nGood = 0; 
+	    for (int i = 0; i < stackData->n; i++) {
+
+		pmStackData *data = stackData->data[i]; // Stack data for this input
+		pmReadout *ro = data->readout;  // Readout of interest
+		psImage *image = ro->image;
+		psImage *mask = ro->mask;
+
+		int xIn = x - data->readout->col0;
+		int yIn = y - data->readout->row0; // Coordinates on input readout
+
+		if (!isfinite(image->data.F32[yIn][xIn])) continue;
+		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];
+		nGood ++;
+	    }
+	    pixelData->n = nGood;
+	    psVectorSortInPlace (pixelData);
+
+	    if (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] = 1;
+		continue;
+	    }
+
+# 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 
+	    int Npt = Ne - Ns;
+
+	    float sum = 0.0;
+	    for (int n = Ns; n < Ne; n++) {
+		sum += pixelData->data.F32[n];
+	    }
+	    float mean = sum / (float) Npt;
+
+	    float varSum = 0.0;
+	    for (int n = Ns; n < Ne; n++) {
+		varSum += SQ(pixelData->data.F32[n] - mean);
+	    }
+	    // 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;
+	    combined->variance->data.F32[y][x] = varValue;
+	    combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+	}
+    }
+  
+    psFree(pixelData);
+    psFree(pixelMask);
+
+    return (true);
+}
+
 /// Stack input images
 bool pmStackCombine(