Index: branches/pap/psModules/src/imcombine/pmStack.c
===================================================================
--- branches/pap/psModules/src/imcombine/pmStack.c	(revision 25960)
+++ branches/pap/psModules/src/imcombine/pmStack.c	(revision 25964)
@@ -290,4 +290,5 @@
 
 // Mark a pixel for inspection
+// Value in pixel doesn't seem to agree with the stack, so need to look closer
 static inline void combineInspect(const psArray *inputs, // Stack data
                                   int x, int y, // Pixel
@@ -301,4 +302,21 @@
     }
     data->inspect = psPixelsAdd(data->inspect, data->inspect->nalloc, x, y);
+    return;
+}
+
+// Mark a pixel for rejection
+// Cannot possibly inspect this pixel and confirm that it's good.
+// e.g., Only a single input
+static inline void combineReject(const psArray *inputs, // Stack data
+                                 int x, int y, // Pixel
+                                 int source // Source image index
+                                 )
+{
+    pmStackData *data = inputs->data[source]; // Stack data of interest
+    if (!data) {
+        psWarning("Can't find input data for source %d", source);
+        return;
+    }
+    data->reject = psPixelsAdd(data->reject, data->reject->nalloc, x, y);
     return;
 }
@@ -319,8 +337,8 @@
                           float rej, // Number of standard deviations at which to reject
                           float sys,    // Relative systematic error
-                          float discard,// Fraction of values to discard (Olympic weighted mean)
+                          float olympic,// Fraction of values to discard (Olympic weighted mean)
                           bool useVariance, // Use variance for rejection when combining?
                           bool safe,    // Combine safely?
-                          bool rejectInspect, // Reject values marked for inspection from combination?
+                          bool rejection, // Reject values marked for inspection from combination?
                           combineBuffer *buffer // Buffer for combination; to avoid multiple allocations
                          )
@@ -419,13 +437,15 @@
               }
               maskValue = 0;
-          }
-#ifdef TESTING
-          else {
+          } else {
+              if (!rejection) {
+                  combineReject(inputs, x, y, pixelSources->data.U16[0]);
+              }
+#ifdef TESTING
               numRejected = 1;
               if (x == TEST_X && y == TEST_Y) {
                   fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n");
               }
+#endif
           }
-#endif
           break;
       }
@@ -459,5 +479,5 @@
               if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
                   // Not consistent: mark both for inspection
-                  if (rejectInspect) {
+                  if (rejection) {
                       imageValue = NAN;
                       varianceValue = NAN;
@@ -550,5 +570,5 @@
                   }
 #endif
-                  median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, discard, sort);
+                  median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, olympic, sort);
               }
 
@@ -563,5 +583,5 @@
 #define MASK_PIXEL_FOR_INSPECTION() \
     pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; \
-    if (!rejectInspect) { \
+    if (!rejection) { \
         combineInspect(inputs, x, y, pixelSources->data.U16[j]); \
     } \
@@ -605,5 +625,5 @@
           }
 
-          if (rejectInspect && totalClipped > 0) {
+          if (rejection && totalClipped > 0) {
               // Get rid of the masked values
               // The alternative to this is to make combinationMeanVariance() accept a mask
@@ -658,8 +678,8 @@
 // Ensure the input array of pmStackData is valid, and get some details out of it
 static bool validateInputData(bool *haveVariances, // Do we have variance maps in the sky images?
-                              bool *haveRejects, // Do we have lists of rejected pixels?
                               int *num,    // Number of inputs
                               int *numCols, int *numRows, // Size of (sky) images
-                              psArray *input // Input array of pmStackData to validate
+                              const psArray *input, // Input array of pmStackData to validate
+                              const pmReadout *output // Output readout
     )
 {
@@ -687,5 +707,5 @@
         PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false);
     }
-    *haveRejects = (data->reject != NULL);
+    bool haveRejects = (data->reject != NULL); // Do we have rejected pixels?
 
     // Make sure the rest correspond with the first
@@ -700,5 +720,5 @@
             return false;
         }
-        if ((*haveRejects && !data->reject) || (data->reject && !*haveRejects)) {
+        if ((haveRejects && !data->reject) || (data->reject && !haveRejects)) {
             psError(PS_ERR_UNEXPECTED_NULL, true,
                     "The rejected pixels are specified in some but not all inputs.");
@@ -716,4 +736,13 @@
             PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false);
         }
+    }
+
+    PM_ASSERT_READOUT_NON_NULL(output, false);
+    if (output->image) {
+        PS_ASSERT_IMAGE_NON_NULL(output->image, false);
+        PS_ASSERT_IMAGE_TYPE(output->image, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGE_NON_NULL(output->mask, false);
+        PS_ASSERT_IMAGE_TYPE(output->mask, PS_TYPE_IMAGE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(output->image, output->mask, false);
     }
 
@@ -802,13 +831,11 @@
 /// Stack input images
 bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad,
-                    int kernelSize, int numIter, float rej, float sys, float discard,
-                    bool entire, bool useVariance, bool safe, bool rejectInspect)
-{
-    PS_ASSERT_PTR_NON_NULL(combined, false);
+                    int kernelSize, int numIter, float rej, float sys, float olympic,
+                    bool useVariance, bool safe, bool rejection)
+{
     bool haveVariances;                 // Do we have the variance maps?
-    bool haveRejects;                   // Do we have lists of rejected pixels?
     int num;                            // Number of inputs
     int numCols, numRows;               // Size of (sky) images
-    if (!validateInputData(&haveVariances, &haveRejects, &num, &numCols, &numRows, input)) {
+    if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) {
         return false;
     }
@@ -821,12 +848,4 @@
         PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
     }
-    if (haveRejects) {
-        // This is a subsequent combination, so expect that the image and mask already exist
-        PS_ASSERT_IMAGE_NON_NULL(combined->image, false);
-        PS_ASSERT_IMAGE_TYPE(combined->image, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGE_NON_NULL(combined->mask, false);
-        PS_ASSERT_IMAGE_TYPE(combined->mask, PS_TYPE_IMAGE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false);
-    }
     if (useVariance && !haveVariances) {
         psWarning("Unable to use variance in rejection if no variance maps supplied --- option turned off");
@@ -852,6 +871,12 @@
         }
 #endif
-        if (!haveRejects && !data->inspect) {
-            data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
+        if (!rejection) {
+            // Ensure pixels can be put on the appropriate list
+            if (!data->inspect) {
+                data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
+            }
+            if (!data->reject) {
+                data->reject = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
+            }
         }
     }
@@ -882,91 +907,56 @@
     combineBuffer *buffer = combineBufferAlloc(num);
 
-    if (haveRejects) {
-        psImage *combinedImage = combined->image; // Combined image
-        psImage *combinedMask = combined->mask; // Combined mask
-        psImage *combinedVariance = combined->variance; // Combined variance map
-
-        psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols,
-                                             minInputRows, maxInputRows); // Map of pixels to source
-        psPixels *pixels = NULL;            // Total list of pixels, with no duplicates
+    // Pull the products out, allocate if necessary
+    psImage *combinedImage = combined->image; // Combined image
+    if (!combinedImage) {
+        combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        combinedImage = combined->image;
+    }
+    psImage *combinedMask = combined->mask; // Combined mask
+    if (!combinedMask) {
+        combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+        combinedMask = combined->mask;
+    }
+
+    psImage *combinedVariance = combined->variance; // Combined variance map
+    if (haveVariances && !combinedVariance) {
+        combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        combinedVariance = combined->variance;
+    }
+
+    // Set up rejection list
+    psArray *pixelMap = NULL;           // Map of pixels to source
+    if (rejection) {
+        pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, minInputRows, maxInputRows);
+    }
+
+    // Combine each pixel
+    for (int y = minInputRows; y < maxInputRows; y++) {
+        for (int x = minInputCols; x < maxInputCols; x++) {
+            psVector *reject = NULL; // Images to reject for this pixel
+            if (rejection) {
+                reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, x, y);
+            }
+            combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
+                          addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, olympic,
+                          useVariance, safe, rejection, buffer);
+        }
+    }
+
+    psFree(pixelMap);
+    psFree(weights);
+    psFree(buffer);
+
+#ifndef PS_NO_TRACE
+    if (!rejection && psTraceGetLevel("psModules.imcombine") >= 5) {
         for (int i = 0; i < num; i++) {
-            pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels
-            if (!data) {
+            pmStackData *data = input->data[i]; // Stack data for this input
+            if (!data || !data->inspect) {
                 continue;
             }
-            pixels = psPixelsConcatenate(pixels, data->reject);
-        }
-        pixels = psPixelsDuplicates(pixels, pixels);
-
-        if (entire) {
-            // Combine entire image
-            for (int y = minInputRows; y < maxInputRows; y++) {
-                for (int x = minInputCols; x < maxInputCols; x++) {
-                    psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows,
-                                                     x, y); // Reject these images
-                    combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
-                                  addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
-                                  useVariance, safe, rejectInspect, buffer);
-                }
-            }
-        } else {
-            // Only combine previously rejected pixels
-            for (int i = 0; i < pixels->n; i++) {
-                // Pixel coordinates are in the frame of the original image
-                int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest
-                if (x < minInputCols || x >= maxInputCols || y < minInputRows || y >= maxInputRows) {
-                    continue;
-                }
-                psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows,
-                                                 x, y); // Reject these images
-                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
-                              addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
-                              useVariance, safe, rejectInspect, buffer);
-            }
-        }
-        psFree(pixels);
-        psFree(pixelMap);
-    } else {
-        // Pull the products out, allocate if necessary
-        psImage *combinedImage = combined->image; // Combined image
-        if (!combinedImage) {
-            combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            combinedImage = combined->image;
-        }
-        psImage *combinedMask = combined->mask; // Combined mask
-        if (!combinedMask) {
-            combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-            combinedMask = combined->mask;
-        }
-
-        psImage *combinedVariance = combined->variance; // Combined variance map
-        if (haveVariances && !combinedVariance) {
-            combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            combinedVariance = combined->variance;
-        }
-
-        for (int y = minInputRows; y < maxInputRows; y++) {
-            for (int x = minInputCols; x < maxInputCols; x++) {
-                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
-                              addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard,
-                              useVariance, safe, rejectInspect, buffer);
-            }
-        }
-
-#ifndef PS_NO_TRACE
-        if (psTraceGetLevel("psModules.imcombine") >= 5) {
-            for (int i = 0; i < num; i++) {
-                pmStackData *data = input->data[i]; // Stack data for this input
-                if (!data || !data->inspect) {
-                    continue;
-                }
-                psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n);
-            }
-        }
-#endif
-    }
-
-    psFree(weights);
-    psFree(buffer);
+            psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n);
+        }
+    }
+#endif
 
     return true;