Index: /branches/eam_branches/psModules.20211028/src/imcombine/pmStack.c
===================================================================
--- /branches/eam_branches/psModules.20211028/src/imcombine/pmStack.c	(revision 41871)
+++ /branches/eam_branches/psModules.20211028/src/imcombine/pmStack.c	(revision 41872)
@@ -1514,8 +1514,26 @@
 
 # define MIN_GOOD_PERCENTILE 2
+# define SUSPECT_FRACTION 0.65
+
+// Comparison and swap functions for sorting values directly
+#define SORT_VV_COMPARE(A,B) (pixelData->data.F32[A] < pixelData->data.F32[B])
+#define SORT_VV_SWAP(TYPE,A,B) {					\
+    if (A != B) {							\
+      psF32 tempVal = pixelData->data.F32[A];				\
+      pixelData->data.F32[A] = pixelData->data.F32[B];			\
+      pixelData->data.F32[B] = tempVal;					\
+      psF32 tempVar = pixelVariances->data.F32[A];			\
+      pixelVariances->data.F32[A] = pixelVariances->data.F32[B];		\
+      pixelVariances->data.F32[B] = tempVar;				\
+    }									\
+}
+
+#define SORT_VALUE_AND_VAR(A,B) { PSSORT(A->n, SORT_VV_COMPARE, SORT_VV_SWAP, F32); }
+  
+
 bool pmStackCombineByPercentile(
     pmReadout *combined,
     psArray *stackData,
-    psF64 minRange,
+    psF64 rejectFraction,
     psF64 maxRange,
     psImageMaskType badMaskBits, 	// treat these bits as 'bad'
@@ -1541,6 +1559,14 @@
     psFree(stackReadouts);
 
-    psVector *pixelData = psVectorAlloc(stackData->n, PS_TYPE_F32);
-    psVector *pixelMask = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK);
+    // make sure the output readout matches the inputs, and set to blank by default
+    pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, blankMaskBits);
+
+    psVector *pixelData      = psVectorAlloc(stackData->n, PS_TYPE_F32);
+    psVector *pixelVariances = psVectorAlloc(stackData->n, PS_TYPE_F32);
+    psVector *pixelMask      = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK);
+
+    int nGoodBits[16]; // accumulate the good pixel bits here for fuzzy logic
+    psAssert (sizeof(psImageMaskType) == 2, "psImageMaskType is not the expected size");
+    memset (nGoodBits, 0, 16*sizeof(int));
 
     for (int y = minInputRows; y < maxInputRows; y++) {
@@ -1552,82 +1578,101 @@
 		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;
+
+		psImage *image    = ro->image;
+		psImage *variance = ro->variance;
+		psImage *mask     = ro->mask;
 
 		int xIn = x - data->readout->col0;
 		int yIn = y - data->readout->row0; // Coordinates on input readout
 
+		// skip obviously bad input data
 		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)
+		if (fabs(image->data.F32[yIn][xIn]) > 1e5) continue; // XXX this cut is a bit arbitrary..
+		if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits) continue;
+
 		// 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;
+		psImageMaskType value = 0x0001;
+		for (int nbit = 0; nbit < 16; nbit ++) {
+		  if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & value) {
+		    nGoodBits[nbit] ++;
 		  }
+		  value <<= 1;
 		}
-# endif
+
+		// pixelWeights->data.F32[nGood] = data->weight;
+		// pixelExps->data.F32[nGood] = data->exp;
+		// pixelSources->data.U16[nGood] = i;
 
 		pixelData->data.F32[nGood] = image->data.F32[yIn][xIn];
+		if (variance) {
+		  // pixelVariances->data.F32[numGood] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i];
+		  pixelVariances->data.F32[nGood] = variance->data.F32[yIn][xIn];
+		}
+
 		nGood ++;
 	    }
 	    pixelData->n = nGood;
-	    psVectorSortInPlace (pixelData);
-
-	    if (nGood < MIN_GOOD_PERCENTILE) {
+
+
+	    // rather than define a min and max value,
+	    // what we really want is a symmetric selection about the middle,
+	    // or the output is biased.  If N % 2 = 1, then 
+	    
+	    // we are going to exclude rejectFraction*nGood measurements.  But the
+	    // rejection needs to be symmetric
+
+# define AT_LEAST 0
+# if (AT_LEAST)
+	    int Ns = MIN(MAX(1, 0.5*rejectFraction * nGood), nGood);
+	    int Ne = nGood - Ns;
+	    int Npt = Ne - Ns;
+# else
+	    int Ns = MIN(MAX(0, 0.5*rejectFraction * nGood), nGood);
+	    int Ne = nGood - Ns;
+	    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] = 1;
+		combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = blankMaskBits; // probably not needed since it was set above.
 		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;
+	    // sort both pixelData and pixelVariance
+	    SORT_VALUE_AND_VAR (pixelData, pixelVariances);
+
+	    // 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
+	    psImageMaskType value = 0x0001;
+	    psImageMaskType outputMask = 0x0000;
+	    for (int nbit = 0; nbit < 16; nbit ++) {
+	      if (nGoodBits[nbit] > SUSPECT_FRACTION*nGood) {
+		outputMask |= 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;
+	    float varSum = 0.0;
 	    for (int n = Ns; n < Ne; n++) {
 		sum += pixelData->data.F32[n];
+		varSum += pixelVariances->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);
-	    }
+	    // alternative: calculate the stdev of the pixel values
+	    // 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)
+	    float varValue = varSum / (Npt*Npt);
 
 	    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;
+	    combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = outputMask;
 	}
     }