Changeset 25975 for branches/pap/psModules/src/imcombine/pmStack.c

Timestamp:

Oct 30, 2009, 5:46:42 PM (17 years ago)

Author:

Paul Price

Message:

Reworking mammoth combinePixels() function into discrete parts.
Adding case to test 3 pixels (occurs often when there are 4 input images) which the general case doesn't treat well, often rejecting all inputs.
Not yet tested!

File:

• branches/pap/psModules/src/imcombine/pmStack.c (modified) (14 diffs)

branches/pap/psModules/src/imcombine/pmStack.c

@@ -30 +30 @@
 #define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
 #define PIXEL_MAP_BUFFER 2              // Number of entries to add to pixel map at a time
-#define ADD_VARIANCE                    // Allow additional variance (besides variance factor)?
+//#define ADD_VARIANCE                    // Allow additional variance (besides variance factor)?
 #define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
 
+
 #define TESTING                         // Enable test output
-#define TEST_X 1500-1                     // x coordinate to examine
-#define TEST_Y 4000-1                     // y coordinate to examine
+#define TEST_X 4177-1                     // x coordinate to examine
+#define TEST_Y 2964-1                     // y coordinate to examine
 
 
@@ -291 +292 @@
 // 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
-                                  int source // Source image index
-                                  )
+static inline void combineMarkInspect(const psArray *inputs, // Stack data
+                                      int x, int y, // Pixel
+                                      int source // Source image index
+                                      )
 {
     pmStackData *data = inputs->data[source]; // Stack data of interest
@@ -308 +309 @@
 // 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
-                                 )
+static inline void combineMarkReject(const psArray *inputs, // Stack data
+                                     int x, int y, // Pixel
+                                     int source // Source image index
+                                     )
 {
     pmStackData *data = inputs->data[source]; // Stack data of interest
@@ -322 +323 @@
 }
 
+// General test for multiple pixels
+// Returns false if we need to re-run without suspect pixels
+static bool combineTestGeneral(int num,      // Number of good pixels
+                               bool suspect, // Are there suspect pixels in the list?
+                               psArray *inputs,       // Original inputs (for flagging)
+                               combineBuffer *buffer, // Buffer with vectors
+                               int x, int y, // Coordinates of interest; frame of output image
+                               int numIter, // Number of rejection iterations
+                               float rej, // Number of standard deviations at which to reject
+                               float sys,    // Relative systematic error
+                               float olympic,// Fraction of values to discard (Olympic weighted mean)
+                               bool useVariance, // Use variance for rejection when combining?
+                               bool safe,    // Combine safely?
+                               bool allowSuspect // Allow suspect values?
+                               )
+{
+    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
+    psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest
+    psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest
+    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
+    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
+    psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
+    psVector *sort = buffer->sort;      // Sort buffer
+
+
+    pixelMasks->n = num;
+    psVectorInit(pixelMasks, 0);
+
+    // Set up rejection limits
+    if (useVariance) {
+        // Convert to rejection limits --- saves doing it later.
+        // Using squared rejection limit because it's cheaper than sqrts
+        float rej2 = PS_SQR(rej); // Rejection level squared
+        double sumWeights = 0.0;
+        for (int i = 0; i < num; i++) {
+            sumWeights += pixelWeights->data.F32[i];
+        }
+        for (int i = 0; i < num; i++) {
+            // Systematic error contributes to the rejection level
+            float sysVar = PS_SQR(sys * pixelData->data.F32[i]);
+#ifdef TESTING
+            // Correct variance for comparison against weighted mean including itself
+            float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights;
+            if (x == TEST_X && y == TEST_Y) {
+                fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i],
+                        pixelVariances->data.F32[i], sysVar, compare);
+            }
+#endif
+            pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar);
+        }
+    }
+
+    int numClipped = INT_MAX;     // Number of pixels clipped per iteration
+    int totalClipped = 0;         // Total number of pixels clipped
+    for (int i = 0; i < numIter && numClipped > 0 && num - totalClipped > 2; i++) {
+        numClipped = 0;
+        float median = NAN;       // Middle of distribution
+        float limit = NAN;        // Rejection limit
+        if (!useVariance) {
+            float stdev;  // Median and stdev of the combination, for rejection
+            if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev,
+                                        pixelData, pixelMasks, sort)) {
+                if (i == 0 && suspect) {
+                    // Something's not right --- repeat
+                    return false;
+                }
+                for (int j = 0; j < num; j++) {
+                    combineMarkReject(inputs, x, y, pixelSources->data.U16[j]);
+                }
+                return true;
+            }
+            limit = rej * stdev;
+#ifdef TESTING
+            if (x == TEST_X && y == TEST_Y) {
+                fprintf(stderr, "Flag without variance; limit: %f\n", limit);
+            }
+#endif
+        } else {
+#ifdef TESTING
+            if (x == TEST_X && y == TEST_Y) {
+                fprintf(stderr, "Flag with variance...\n");
+            }
+#endif
+            median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, olympic, sort);
+        }
+
+#ifdef TESTING
+        if (x == TEST_X && y == TEST_Y) {
+            fprintf(stderr, "Median: %f\n", median);
+        }
+#endif
+
+        // Mask a pixel for inspection
+#define MASK_PIXEL_FOR_INSPECTION()                                     \
+        if (i == 0 && suspect) {                                        \
+            /* Something's inconsistent, so want repeat, throwing out all suspect pixels */ \
+            return false; \
+        }                                                               \
+        pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;            \
+        combineMarkInspect(inputs, x, y, pixelSources->data.U16[j]);    \
+        numClipped++;                                                   \
+        totalClipped++;
+        // End
+
+        for (int j = 0; j < num; j++) {
+            if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+                continue;
+            }
+            float diff = pixelData->data.F32[j] - median; // Difference from expected
+#ifdef TESTING
+            if (x == TEST_X && y == TEST_Y) {
+                fprintf(stderr, "Testing input %d: %f\n", j, diff);
+            }
+#endif
+            if (useVariance) {
+                // Comparing squares --- cheaper than lots of sqrts
+                // pixelVariances includes the rejection limit, from above
+                if (PS_SQR(diff) > pixelLimits->data.F32[j]) {
+                    MASK_PIXEL_FOR_INSPECTION();
+#ifdef TESTING
+                    if (x == TEST_X && y == TEST_Y) {
+                        fprintf(stderr, "Flagging input %d based on variance: %f > %f\n",
+                                j, diff, sqrtf(pixelLimits->data.F32[j]));
+                    }
+#endif
+                }
+            } else if (fabsf(diff) > limit) {
+                MASK_PIXEL_FOR_INSPECTION();
+#ifdef TESTING
+                if (x == TEST_X && y == TEST_Y) {
+                    fprintf(stderr, "Flagging input %d based on distribution: %f > %f\n",
+                            j, diff, limit);
+                }
+#endif
+            }
+        }
+    }
+
+    return true;
+}
+
+
+// Extract vectors for simple combination/rejection operations
+static void combineExtract(int *num,                        // Number of good pixels
+                           bool *suspect,                   // Any suspect pixels?
+                           combineBuffer *buffer, // Buffer with vectors
+                           psImage *image, // Combined image, for output
+                           psImage *mask, // Combined mask, for output
+                           psImage *variance, // Combined variance map, for output
+                           const psArray *inputs, // Stack data
+                           const psVector *weights, // Global (single value) weights for data, or NULL
+                           const psVector *addVariance, // Additional variance for data
+                           const psVector *reject, // Indices of pixels to reject, or NULL
+                           int x, int y, // Coordinates of interest; frame of output image
+                           psImageMaskType maskVal, // Value to mask
+                           psImageMaskType maskSuspect, // Value to suspect
+                           bool allowSuspect        // Allow suspect values?
+                           )
+{
+    // Rudimentary error checking
+    assert(buffer);
+    assert(image);
+    assert(mask);
+    assert(inputs);
+
+    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
+    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
+    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
+    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
+    psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
+
+    // Extract the pixel and mask data
+    int numGood = 0;                    // Number of good pixels
+    for (int i = 0, j = 0; i < inputs->n; i++) {
+        // Check if this pixel has been rejected.  Assumes that the rejection pixel list is sorted --- it
+        // should be because of how pixelMapGenerate works
+        if (reject && reject->data.U16[j] == i) {
+            j++;
+            continue;
+        }
+
+        pmStackData *data = inputs->data[i]; // Stack data of interest
+        if (!data) {
+            continue;
+        }
+
+        int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
+        psImage *mask = data->readout->mask; // Mask of interest
+        if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskVal) {
+            continue;
+        }
+        if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect) {
+            if (!allowSuspect) {
+                combineMarkReject(inputs, x, y, i);
+                continue;
+            }
+            if (suspect) {
+                *suspect = true;
+            }
+        }
+
+        psImage *image = data->readout->image; // Image of interest
+        psImage *variance = data->readout->variance; // Variance map of interest
+        pixelData->data.F32[numGood] = image->data.F32[yIn][xIn];
+        if (variance) {
+            pixelVariances->data.F32[numGood] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i];
+        }
+        pixelWeights->data.F32[numGood] = data->weight;
+        pixelSources->data.U16[numGood] = i;
+        numGood++;
+    }
+    pixelData->n = numGood;
+    if (variance) {
+        pixelVariances->n = numGood;
+    }
+    pixelWeights->n = numGood;
+    pixelSources->n = numGood;
+    pixelLimits->n = numGood;
+    *num = numGood;
+
+#ifdef TESTING
+    if (x == TEST_X && y == TEST_Y) {
+        for (int i = 0; i < numGood; i++) {
+            fprintf(stderr, "Input %d (%" PRIu16 "): %f %f (%f) %f\n",
+                    i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
+                    addVariance->data.F32[i], pixelWeights->data.F32[i]);
+        }
+    }
+#endif
+
+    return;
+}
+
+
+// Combine pixels
+static void combinePixels(psImage *image, // Combined image, for output
+                          psImage *mask, // Combined mask, for output
+                          psImage *variance, // Combined variance map, for output
+                          int num,      // Number of good pixels
+                          combineBuffer *buffer, // Buffer with vectors
+                          int x, int y, // Coordinates of interest; frame of output image
+                          psImageMaskType bad, // Value for bad pixels
+                          bool safe             // Safe combination?
+                          )
+{
+    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
+    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
+    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
+
+    // Default option is that the pixel is bad
+    float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
+    psImageMaskType maskValue = bad;    // Value for combined mask
+    switch (num) {
+      case 0: {
+          // Nothing to combine: it's bad
+#ifdef TESTING
+          if (x == TEST_X && y == TEST_Y) {
+              fprintf(stderr, "No inputs to combine, pixel is bad.\n");
+          }
+#endif
+          break;
+      }
+      case 1: {
+          // Accept the single pixel unless we have to be safe
+          if (!safe) {
+#ifdef TESTING
+              if (x == TEST_X && y == TEST_Y) {
+                  fprintf(stderr, "Single input to combine, safety off.\n");
+              }
+#endif
+              imageValue = pixelData->data.F32[0];
+              if (variance) {
+                  varianceValue = pixelVariances->data.F32[0];
+              }
+              maskValue = 0;
+          }
+#ifdef TESTING
+          else if (x == TEST_X && y == TEST_Y) {
+              fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n");
+          }
+#endif
+          break;
+      }
+      case 2: {
+          // Automatically accept the mean of the pixels only if we're not playing safe
+          if (!safe) {
+              float mean, var;   // Mean and variance from combination
+              if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
+                  imageValue = mean;
+                  varianceValue = var;
+                  maskValue = 0;
+#ifdef TESTING
+                  if (x == TEST_X && y == TEST_Y) {
+                      fprintf(stderr, "Two inputs to combine using unsafe --> %f %f\n", mean, var);
+                  }
+#endif
+              }
+          }
+#ifdef TESTING
+          else {
+              if (x == TEST_X && y == TEST_Y) {
+                  fprintf(stderr, "Two inputs to combine, safety on, pixel is bad\n");
+              }
+          }
+#endif
+          break;
+      }
+      default: {
+          // Can combine without too much worrying
+          float mean, var;           // Mean and variance of the combination
+          if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
+              break;
+          }
+          imageValue = mean;
+          varianceValue = var;
+          maskValue = 0;
+#ifdef TESTING
+          if (x == TEST_X && y == TEST_Y) {
+              fprintf(stderr, "Combined inputs: %f %f\n", mean, var);
+          }
+#endif
+          break;
+      }
+    }
+
+    image->data.F32[y][x] = imageValue;
+    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskValue;
+    if (variance) {
+        variance->data.F32[y][x] = varianceValue;
+    }
+
+    return;
+}
+
+
+// Test pixels to be combined
+// Returns false to repeat without suspect pixels
+static bool combineTest(int num,      // Number of good pixels
+                        bool suspect, // Does the stack contain suspect pixels?
+                        psArray *inputs,       // Original inputs (for flagging)
+                        combineBuffer *buffer, // Buffer with vectors
+                        int x, int y, // Coordinates of interest; frame of output image
+                        int numIter, // Number of rejection iterations
+                        float rej, // Number of standard deviations at which to reject
+                        float sys,    // Relative systematic error
+                        float olympic,// Fraction of values to discard (Olympic weighted mean)
+                        bool useVariance, // Use variance for rejection when combining?
+                        bool safe,    // Combine safely?
+                        bool allowSuspect    // Allow suspect pixels in the combination?
+                        )
+{
+    if (numIter <= 0) {
+        return true;
+    }
+
+    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
+    psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest
+    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
+
+    switch (num) {
+      case 0:
+        break;
+      case 1:
+          if (safe) {
+              combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
+          }
+          break;
+      case 2: {
+          if (useVariance) {
+              // Use variance to check that the two are consistent
+              float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference
+              float var1 = pixelVariances->data.F32[0]; // Variance of first
+              float var2 = pixelVariances->data.F32[1]; // Variance of second
+              // Systematic error contributes to the rejection level
+              var1 += PS_SQR(sys * pixelData->data.F32[0]);
+              var2 += PS_SQR(sys * pixelData->data.F32[1]);
+
+              float sigma2 = var1 + var2; // Combined variance
+              if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
+                  // Not consistent: don't believe either!
+                  if (allowSuspect && suspect) {
+                      combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
+                      combineMarkReject(inputs, x, y, pixelSources->data.U16[1]);
+                  } else {
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
+                  }
+#ifdef TESTING
+                  if (x == TEST_X && y == TEST_Y) {
+                      fprintf(stderr, "Flagged both pixels (%f > %f x %f\n)",
+                              diff, rej, sqrtf(sigma2));
+                  }
+#endif
+              }
+          } else if (safe) {
+              // Can't test them, and we want to be safe, so reject
+              combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
+              combineMarkReject(inputs, x, y, pixelSources->data.U16[1]);
+          }
+          break;
+      }
+      case 3: {
+          // Want to be a bit careful on the rejection than for a larger number of inputs
+          if (!useVariance) {
+              return combineTestGeneral(num, suspect, inputs, buffer, x, y, numIter, rej, sys,
+                                        olympic, useVariance, safe, allowSuspect);
+          }
+
+          // Differences between pixel values
+          float diff01 = pixelData->data.F32[0] - pixelData->data.F32[1];
+          float diff12 = pixelData->data.F32[1] - pixelData->data.F32[2];
+          float diff20 = pixelData->data.F32[2] - pixelData->data.F32[0];
+          // Variance for each pixel
+          float var0 = pixelVariances->data.F32[0] + PS_SQR(sys * pixelData->data.F32[0]);
+          float var1 = pixelVariances->data.F32[1] + PS_SQR(sys * pixelData->data.F32[1]);
+          float var2 = pixelVariances->data.F32[2] + PS_SQR(sys * pixelData->data.F32[2]);
+          float rej2 = PS_SQR(rej); // Rejection level squared
+          // Errors in pixel differences
+          float err01 = var0 + var1;
+          float err12 = var1 + var2;
+          float err20 = var2 + var0;
+
+          int badPairs = 0;         // Number of bad pairs
+          bool bad01 = false, bad12 = false, bad20 = false; // Pair is bad?
+          if (PS_SQR(diff01) > rej2 * err01) {
+              bad01 = true;
+              badPairs++;
+          }
+          if (PS_SQR(diff12) > rej2 * err12) {
+              bad12 = true;
+              badPairs++;
+          }
+          if (PS_SQR(diff20) > rej2 * err20) {
+              bad20 = true;
+              badPairs++;
+          }
+
+          if (badPairs > 0 && allowSuspect && suspect) {
+              return false;
+          }
+
+          switch (badPairs) {
+            case 0:
+              // Nothing to worry about!
+              break;
+            case 1:
+              // Can't tell which image is bad, so be sure to get it
+              if (bad01) {
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
+                  break;
+              }
+              if (bad12) {
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
+                  break;
+              }
+              if (bad20) {
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+                  break;
+              }
+              psAbort("Should never get here");
+            case 2:
+              if (bad01 && bad12) {
+                  // 2 and 0 are good
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
+                  break;
+              }
+              if (bad12 && bad20) {
+                  // 0 and 1 are good
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
+                  break;
+              }
+              if (bad20 && bad01) {
+                  // 1 and 2 are good
+                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+                  break;
+              }
+              psAbort("Should never get here");
+            case 3:
+              // Everything's bad
+              combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+              combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
+              combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
+              break;
+          }
+          break;
+      }
+      default: {
+          return combineTestGeneral(num, suspect, inputs, buffer, x, y, numIter, rej, sys,
+                                    olympic, useVariance, safe, allowSuspect);
+      }
+    }
+
+    return true;
+}
+
+
+
+#if 0
 // Given a stack of images, combine with optional rejection.
 // Pixels in the stack that are rejected are marked for subsequent inspection
@@ -333 +835 @@
                           int x, int y, // Coordinates of interest; frame of output image
                           psImageMaskType maskVal, // Value to mask
+                          psImageMaskType suspect, // Value to suspect
                           psImageMaskType bad, // Value to give bad pixels
                           int numIter, // Number of rejection iterations
@@ -341 +844 @@
                           bool safe,    // Combine safely?
                           bool rejection, // Reject values marked for inspection from combination?
-                          combineBuffer *buffer // Buffer for combination; to avoid multiple allocations
+                          combineBuffer *buffer, // Buffer for combination; to avoid multiple allocations
+                          bool allowSuspect    // Allow suspect pixels in the combination?
                          )
 {
@@ -361 +865 @@
     psVector *sort = buffer->sort;      // Sort buffer
 
-    // Extract the pixel and mask data
-    int num = 0;                        // Number of good images
-    for (int i = 0, j = 0; i < inputs->n; i++) {
-        // Check if this pixel has been rejected.  Assumes that the rejection pixel list is sorted --- it
-        // should be because of how pixelMapGenerate works
-        if (reject && reject->data.U16[j] == i) {
-            j++;
-            continue;
-        }
-
-        pmStackData *data = inputs->data[i]; // Stack data of interest
-        if (!data) {
-            continue;
-        }
-
-        int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
-        psImage *mask = data->readout->mask; // Mask of interest
-        if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskVal) {
-            continue;
-        }
-
-        psImage *image = data->readout->image; // Image of interest
-        psImage *variance = data->readout->variance; // Variance map of interest
-        pixelData->data.F32[num] = image->data.F32[yIn][xIn];
-        if (variance) {
-            pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i];
-        }
-        pixelWeights->data.F32[num] = data->weight;
-        pixelSources->data.U16[num] = i;
-        num++;
-    }
-    pixelData->n = num;
-    if (variance) {
-        pixelVariances->n = num;
-    }
-    pixelWeights->n = num;
-    pixelSources->n = num;
-    pixelLimits->n = num;
-
-#ifdef TESTING
-    if (x == TEST_X && y == TEST_Y) {
-        for (int i = 0; i < num; i++) {
-            fprintf(stderr, "Input %d (%" PRIu16 "): %f %f %f\n",
-                    i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
-                    pixelWeights->data.F32[i]);
-        }
-    }
-    int numRejected = 0;                // Number of rejected inputs
-#endif
 
     // The sensible thing to do varies according to how many good pixels there are.
@@ -439 +894 @@
           } else {
               if (!rejection) {
-                  combineReject(inputs, x, y, pixelSources->data.U16[0]);
+                  combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
               }
 #ifdef TESTING
@@ -478 +933 @@
               float sigma2 = var1 + var2; // Combined variance
               if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
-                  // Not consistent: mark both for inspection
+                  // Not consistent: don't believe either!
                   if (rejection) {
                       imageValue = NAN;
                       varianceValue = NAN;
                       maskValue = bad;
+                  } else if (allowSuspect && numSuspect > 0) {
+                      combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
+                      combineMarkReject(inputs, x, y, pixelSources->data.U16[1]);
                   } else {
-                      combineInspect(inputs, x, y, pixelSources->data.U16[0]);
-                      combineInspect(inputs, x, y, pixelSources->data.U16[1]);
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
+                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
                   }
 #ifdef TESTING
@@ -494 +952 @@
                   numRejected = 2;
 #endif
+              }
+          }
+          break;
+      }
+      case 3: {
+          // Can combine without too much worrying, but want to be a bit careful on the rejection
+          float mean, var;           // Mean and variance of the combination
+          if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
+              break;
+          }
+          imageValue = mean;
+          varianceValue = var;
+          maskValue = 0;
+#ifdef TESTING
+          if (x == TEST_X && y == TEST_Y) {
+              fprintf(stderr, "Combined inputs: %f %f\n", mean, var);
+          }
+#endif
+
+          if (numIter > 0) {
+              if (!useVariance) {
+                  combineRejectionGeneral();
+              } else {
+                  // Differences between pixel values
+                  float diff01 = pixelData->data.F32[0] - pixelData->data.F32[1];
+                  float diff12 = pixelData->data.F32[1] - pixelData->data.F32[2];
+                  float diff20 = pixelData->data.F32[2] - pixelData->data.F32[0];
+                  // Variance for each pixel
+                  float var0 = pixelVariances->data.F32[0] + PS_SQR(sys * pixelData->data.F32[0]);
+                  float var1 = pixelVariances->data.F32[1] + PS_SQR(sys * pixelData->data.F32[1]);
+                  float var2 = pixelVariances->data.F32[2] + PS_SQR(sys * pixelData->data.F32[2]);
+                  float rej2 = PS_SQR(rej); // Rejection level squared
+                  // Errors in pixel differences
+                  float err01 = var0 + var1;
+                  float err12 = var1 + var2;
+                  float err20 = var2 + var0;
+
+                  int badPairs = 0;         // Number of bad pairs
+                  bool bad01 = false, bad12 = false, bad20 = false; // Pair is bad?
+                  if (PS_SQR(diff01) > rej2 * err01) {
+                      bad01 = true;
+                      badPairs++;
+                  }
+                  if (PS_SQR(diff12) > rej2 * err12) {
+                      bad12 = true;
+                      badPairs++;
+                  }
+                  if (PS_SQR(diff20) > rej2 * err20) {
+                      bad20 = true;
+                      badPairs++;
+                  }
+
+                  switch (badPairs) {
+                    case 0:
+                      // Nothing to worry about!
+                      break;
+                    case 1:
+                      // Can't tell which image is bad, so be sure to get it
+                      if (bad01) {
+                          combineInspect(inputs, x, y, pixelSources->data.U16[0]);
+                          combineInspect(inputs, x, y, pixelSources->data.U16[1]);
+                          break;
+                      }
+                      if (bad12) {
+                          combineInspect(inputs, x, y, pixelSources->data.U16[1]);
+                          combineInspect(inputs, x, y, pixelSources->data.U16[2]);
+                          break;
+                      }
+                      if (bad20) {
+                          combineInspect(inputs, x, y, pixelSources->data.U16[2]);
+                          combineInspect(inputs, x, y, pixelSources->data.U16[0]);
+                          break;
+                      }
+                      psAbort("Should never get here");
+                    case 2:
+                      if (bad01 && bad12) {
+                          // 2 and 0 are good
+                          combineInspect(inputs, x, y, pixelSources->data.U16[1]);
+                          break;
+                      }
+                      if (bad12 && bad20) {
+                          // 0 and 1 are good
+                          combineInspect(inputs, x, y, pixelSources->data.U16[2]);
+                          break;
+                      }
+                      if (bad20 && bad01) {
+                          // 1 and 2 are good
+                          combineInspect(inputs, x, y, pixelSources->data.U16[0]);
+                          break;
+                      }
+                      psAbort("Should never get here");
+                    case 3:
+                      // Everything's bad
+                      combineInspect(inputs, x, y, pixelSources->data.U16[0]);
+                      combineInspect(inputs, x, y, pixelSources->data.U16[1]);
+                      combineInspect(inputs, x, y, pixelSources->data.U16[2]);
+                      break;
+                  }
               }
           }
@@ -515 +1071 @@
 
           // Prepare for clipping iteration
-          if (numIter > 0) {
-              pixelMasks->n = num;
-              psVectorInit(pixelMasks, 0);
-              if (useVariance) {
-                  // Convert to rejection limits --- saves doing it later.
-                  // Using squared rejection limit because it's cheaper than sqrts
-                  float rej2 = PS_SQR(rej); // Rejection level squared
-                  double sumWeights = 0.0;
-                  for (int i = 0; i < num; i++) {
-                      sumWeights += pixelWeights->data.F32[i];
-                  }
-                  for (int i = 0; i < num; i++) {
-                      // Systematic error contributes to the rejection level
-                      float sysVar = PS_SQR(sys * pixelData->data.F32[i]);
-#ifdef TESTING
-                      // Correct variance for comparison against weighted mean including itself
-                      float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights;
-                      if (x == TEST_X && y == TEST_Y) {
-                          fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i],
-                                  pixelVariances->data.F32[i], sysVar, compare);
-                      }
-#endif
-
-                      pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar);
-                  }
-              }
-          }
-
-          // The clipping that follows is solely to identify suspect pixels.
-          // These suspect pixels will be inspected in more detail by other functions.
-          int numClipped = INT_MAX;     // Number of pixels clipped per iteration
-          int totalClipped = 0;         // Total number of pixels clipped
-          for (int i = 0; i < numIter && numClipped > 0 && num - totalClipped > 2; i++) {
-              numClipped = 0;
-              float median = NAN;       // Middle of distribution
-              float limit = NAN;        // Rejection limit
-              if (!useVariance) {
-                  float stdev;  // Median and stdev of the combination, for rejection
-                  if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev,
-                                              pixelData, pixelMasks, sort)) {
-                      psWarning("Bad median/stdev at %d,%d", x, y);
-                      break;
-                  }
-                  limit = rej * stdev;
-#ifdef TESTING
-                  if (x == TEST_X && y == TEST_Y) {
-                      fprintf(stderr, "Rejecting without variance; rejection limit: %f\n", limit);
-                  }
-#endif
-              } else {
-#ifdef TESTING
-                  if (x == TEST_X && y == TEST_Y) {
-                      fprintf(stderr, "Rejecting with variance...\n");
-                  }
-#endif
-                  median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, olympic, sort);
-              }
-
-#ifdef TESTING
-              if (x == TEST_X && y == TEST_Y) {
-                  fprintf(stderr, "Median: %f\n", median);
-              }
-#endif
-
-
-// Mask a pixel for inspection
-#define MASK_PIXEL_FOR_INSPECTION() \
-    pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; \
-    if (!rejection) { \
-        combineInspect(inputs, x, y, pixelSources->data.U16[j]); \
-    } \
-    numClipped++; \
-    totalClipped++;
-
-              for (int j = 0; j < num; j++) {
-                  if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
-                      continue;
-                  }
-                  float diff = pixelData->data.F32[j] - median; // Difference from expected
-#ifdef TESTING
-                  if (x == TEST_X && y == TEST_Y) {
-                      fprintf(stderr, "Testing input %d for rejection: %f\n", j, diff);
-                  }
-#endif
-                  if (useVariance) {
-                      // Comparing squares --- cheaper than lots of sqrts
-                      // pixelVariances includes the rejection limit, from above
-                      if (PS_SQR(diff) > pixelLimits->data.F32[j]) {
-                          MASK_PIXEL_FOR_INSPECTION();
-#ifdef TESTING
-                          if (x == TEST_X && y == TEST_Y) {
-                              fprintf(stderr, "Rejecting input %d based on variance: %f > %f\n",
-                                      j, diff, sqrtf(pixelLimits->data.F32[j]));
-                          }
-                          numRejected++;
-#endif
-                      }
-                  } else if (fabsf(diff) > limit) {
-                      MASK_PIXEL_FOR_INSPECTION();
-#ifdef TESTING
-                      if (x == TEST_X && y == TEST_Y) {
-                          fprintf(stderr, "Rejecting input %d based on distribution: %f > %f\n",
-                                  j, diff, limit);
-                      }
-                      numRejected++;
-#endif
-                  }
-              }
-          }
-
-          if (rejection && totalClipped > 0) {
-              // Get rid of the masked values
-              // The alternative to this is to make combinationMeanVariance() accept a mask
-              int good = 0;            // Index of good value
-              for (int i = 0; i < num; i++) {
-                  if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-                      continue;
-                  }
-                  if (i != good) {
-                      pixelData->data.F32[good] = pixelData->data.F32[i];
-                      pixelVariances->data.F32[good] = pixelVariances->data.F32[i];
-                      pixelWeights->data.F32[good] = pixelWeights->data.F32[i];
-                      pixelData->data.F32[good] = pixelData->data.F32[i];
-                  }
-                  good++;
-              }
-              pixelData->n = good;
-              pixelVariances->n = good;
-              pixelWeights->n = good;
-              if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
-                  imageValue = mean;
-                  varianceValue = var;
-                  maskValue = 0;
-              } else {
-                  imageValue = NAN;
-                  varianceValue = NAN;
-                  maskValue = bad;
-              }
-          }
 
           break;
@@ -659 +1076 @@
     }
 
-    image->data.F32[y][x] = imageValue;
-    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskValue;
-    if (variance) {
-        variance->data.F32[y][x] = varianceValue;
-    }
-
-#ifdef TESTING
-    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
-    if (variance) {
-        variance->data.F32[y][x] = num > 0 ? 1.0 - (float)numRejected / (float)num : 0.0;
-    }
-#endif
-
     return;
 }
-
+#endif
 
 // Ensure the input array of pmStackData is valid, and get some details out of it
@@ -830 +1234 @@
 
 /// Stack input images
-bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad,
-                    int kernelSize, int numIter, float rej, float sys, float olympic,
+bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect,
+                    psImageMaskType bad, int kernelSize, int numIter, float rej, float sys, float olympic,
                     bool useVariance, bool safe, bool rejection)
 {
@@ -946 +1350 @@
                 if (x == TEST_X && y == TEST_Y) {
                     fprintf(stderr, "Rejected inputs: ");
-                    for (int i = 0; i < reject->n; i++) {
-                        fprintf(stderr, "%d ", reject->data.U16[i]);
+                    if (!reject) {
+                        fprintf(stderr, "<none>\n");
+                    } else {
+                        for (int i = 0; i < reject->n; i++) {
+                            fprintf(stderr, "%d ", reject->data.U16[i]);
+                        }
+                        fprintf(stderr, "\n");
                     }
-                    fprintf(stderr, "\n");
                 }
 #endif
             }
-            combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
-                          addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, olympic,
-                          useVariance, safe, rejection, buffer);
+
+            int num;                    // Number of good pixels
+            bool suspect;               // Suspect pixels in stack?
+            combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance,
+                           input, weights, addVariance, reject, x, y, maskVal, maskSuspect, true);
+            if (numIter == 0) {
+                combinePixels(combinedImage, combinedMask, combinedVariance,
+                              num, buffer, x, y, bad, safe);
+            } else {
+                if (combineTest(num, suspect, input, buffer, x, y, numIter, rej, sys, olympic,
+                                useVariance, safe, true)) {
+                    // Need to repeat without suspect pixels
+                    combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance,
+                           input, weights, addVariance, reject, x, y, maskVal, maskSuspect, false);
+                    combineTest(num, suspect, input, buffer, x, y, numIter, rej, sys, olympic,
+                                useVariance, safe, false);
+                }
+            }
         }
     }