Changeset 26076 for trunk/psModules

Timestamp:

Nov 9, 2009, 2:53:12 PM (17 years ago)

Author:

Paul Price

Message:

Merging branches/pap (unconvolved stacks, reworked combinePixels function into clearer roles, only throw out most variant pixel on each iteration, throw suspect pixels out as first rejection step) stack development. I think I've got everything, but not entirely sure, since I've already merged this branch once before (for dual convolution).

Location:

trunk/psModules/src

Files:

• camera (modified) (1 prop)
• camera/pmFPAMaskWeight.c (modified) (1 diff)
• camera/pmFPAMaskWeight.h (modified) (1 diff)
• camera/pmReadoutFake.c (modified) (1 prop)
• imcombine/pmStack.c (modified) (23 diffs, 1 prop)
• imcombine/pmStack.h (modified) (1 diff, 1 prop)
• objects (modified) (1 prop)
• objects/pmSourceMatch.c (modified) (18 diffs)
• objects/pmSourceMatch.h (modified) (5 diffs)

trunk/psModules/src/camera/pmFPAMaskWeight.c

@@ -444 +444 @@
     }
 
-    return true;
+    return psMetadataAddF32(readout->analysis, PS_LIST_TAIL, PM_READOUT_ANALYSIS_RENORM, 0,
+                            "Renormalisation of variance", PS_SQR(correction));
 }
 

trunk/psModules/src/camera/pmFPAMaskWeight.h

@@ -15 +15 @@
 /// @addtogroup Camera Camera Layout
 /// @{
+
+#define PM_READOUT_ANALYSIS_RENORM "READOUT.RENORM" // Name on analysis metadata for renormalisation
 
 /// Set a temporary readout mask using CELL.SATURATION and CELL.BAD

trunk/psModules/src/imcombine/pmStack.c

@@ -33 +33 @@
 #define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
 
+
 //#define TESTING                         // Enable test output
-//#define TEST_X 1085                     // x coordinate to examine
-//#define TEST_Y 3371                     // y coordinate to examine
+//#define TEST_X 3122-1                     // x coordinate to examine
+//#define TEST_Y 1028-1                     // y coordinate to examine
 
 
@@ -42 +43 @@
 typedef struct {
     psVector *pixels;                   // Pixel values
-    psVector *masks;                    // Pixel masks
     psVector *variances;                // Pixel variances
     psVector *weights;                  // Pixel weightings
     psVector *sources;                  // Pixel sources (which image did they come from?)
     psVector *limits;                   // Rejection limits
+    psVector *suspects;                 // Pixel is suspect?
     psVector *sort;                     // Buffer for sorting (to get a robust estimator of the standard dev)
 } combineBuffer;
@@ -53 +54 @@
 {
     psFree(buffer->pixels);
-    psFree(buffer->masks);
     psFree(buffer->variances);
     psFree(buffer->weights);
     psFree(buffer->sources);
     psFree(buffer->limits);
+    psFree(buffer->suspects);
     psFree(buffer->sort);
     return;
@@ -69 +70 @@
 
     buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32);
-    buffer->masks = psVectorAlloc(numImages, PS_TYPE_VECTOR_MASK);
     buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16);
     buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32);
+    buffer->suspects = psVectorAlloc(numImages, PS_TYPE_U8);
     buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32);
     return buffer;
@@ -143 +144 @@
                                    float *stdev, // Standard deviation value, to return
                                    const psVector *values, // Values to combine
-                                   const psVector *masks, // Mask to apply
                                    psVector *sortBuffer // Buffer for sorting
                                    )
 {
     assert(values);
-    assert(!masks || values->n == masks->n);
     assert(values->type.type == PS_TYPE_F32);
-    assert(!masks || masks->type.type == PS_TYPE_VECTOR_MASK);
     assert(sortBuffer && sortBuffer->nalloc >= values->n && sortBuffer->type.type == PS_TYPE_F32);
 
-    // Need to filter out clipped values
-    int num = 0;            // Number of valid values
-    for (int i = 0; i < values->n; i++) {
-        if (!masks || !masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-            sortBuffer->data.F32[num++] = values->data.F32[i];
-        }
-    }
-    sortBuffer->n = num;
-    if (!psVectorSortInPlace(sortBuffer)) {
+    int num = values->n;                // Number of values
+    sortBuffer = psVectorSortIndex(sortBuffer, values);
+    if (!sortBuffer) {
+        *median = NAN;
+        *stdev = NAN;
         return false;
     }
@@ -167 +161 @@
     if (num == 3) {
         // Attempt to measure standard deviation with only three values (and one of those possibly corrupted)
-        *median = sortBuffer->data.F32[1];
+        *median = values->data.F32[sortBuffer->data.S32[1]];
         if (stdev) {
-            float diff1 = sortBuffer->data.F32[0] - *median;
-            float diff2 = sortBuffer->data.F32[2] - *median;
+            float diff1 = values->data.F32[sortBuffer->data.S32[0]] - *median;
+            float diff2 = values->data.F32[sortBuffer->data.S32[2]] - *median;
             // This factor of sqrt(2) might not be exact, but it's about right
             *stdev = M_SQRT2 * PS_MIN(fabsf(diff1), fabsf(diff2));
         }
     } else {
-        *median = num % 2 ? sortBuffer->data.F32[num / 2] :
-            (sortBuffer->data.F32[num / 2 - 1] + sortBuffer->data.F32[num / 2]) / 2.0;
+        *median = num % 2 ? values->data.F32[sortBuffer->data.S32[num / 2]] :
+            (values->data.F32[sortBuffer->data.S32[num / 2 - 1]] +
+             values->data.F32[sortBuffer->data.S32[num / 2]]) / 2.0;
         if (stdev) {
             if (num <= NUM_DIRECT_STDEV) {
@@ -182 +177 @@
                 double sum = 0.0;
                 for (int i = 0; i < num; i++) {
-                    sum += PS_SQR(sortBuffer->data.F32[i] - *median);
+                    sum += PS_SQR(values->data.F32[sortBuffer->data.S32[i]] - *median);
                 }
                 *stdev = sqrt(sum / (double)(num - 1));
             } else {
                 // Standard deviation from the interquartile range
-                *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] -
-                                 sortBuffer->data.F32[(int)(0.25 * num)]);
+                *stdev = 0.74 * (values->data.F32[sortBuffer->data.S32[(int)(0.75 * num)]] -
+                                 values->data.F32[sortBuffer->data.S32[(int)(0.25 * num)]]);
             }
         }
@@ -195 +190 @@
     return true;
 }
-
-#if 0
-// Return the weighted median for the pixels
-// This does not appear to produce as clean images as the weighted Olympic mean
-static float combinationWeightedMedian(const psVector *values, // Values to combine
-                                       const psVector *weights, // Weights to combine
-                                       const psVector *masks, // Mask to apply
-                                       psVector *sortBuffer // Buffer for sorting
-                                       )
-{
-    double sumWeight = 0.0;             // Sum of weights
-    for (int j = 0; j < values->n; j++) {
-        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
-            continue;
-        }
-        sumWeight += weights->data.F32[j];
-    }
-
-    sortBuffer = psVectorSortIndex(sortBuffer, values);
-    double target = sumWeight / 2.0;    // Target weight
-
-    int dominant = -1;                  // Index of dominant value, if any
-    double cumulativeWeight = 0.0;      // Sum of weights
-    for (int j = 0; j < values->n; j++) {
-        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
-            continue;
-        }
-        int index = sortBuffer->data.S32[j];  // Index of value of interest
-        float weight = weights->data.F32[index]; // Weight for value of interest
-        if (weight >= target) {
-            // Get the weighted median of the rest
-            dominant = index;
-            sumWeight -= weight;
-            target = sumWeight / 2.0;
-            continue;
-        }
-        cumulativeWeight += weight;
-        if (cumulativeWeight >= target) {
-            float median = values->data.F32[index]; // Weighted median median
-            if (dominant != -1) {
-                // In the case that a single value contains a disproportionate weight compared to the rest,
-                // we use a weighted mean between that dominant value and the weighted median of the rest.
-                return (values->data.F32[dominant] * weights->data.F32[dominant] + median * sumWeight) /
-                    (weights->data.F32[dominant] + sumWeight);
-            }
-            return median;
-        }
-    }
-
-    return NAN;
-}
-#endif
 
 // Return the weighted Olympic mean for the pixels
 static float combinationWeightedOlympic(const psVector *values, // Values to combine
                                         const psVector *weights, // Weights to combine
-                                        const psVector *masks, // Mask to apply
                                         float frac, // Fraction to discard
                                         psVector *sortBuffer // Buffer for sorting
                                         )
 {
-    int numGood = 0;                    // Number of good values
-    for (int i = 0; i < values->n; i++) {
-        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-            continue;
-        }
-        numGood++;
-    }
+    int numGood = values->n;            // Number of good values
 
     int numBad = frac * numGood + 0.5;  // Number of bad values
@@ -272 +208 @@
     for (int i = 0, j = 0; i < values->n; i++) {
         int index = sortBuffer->data.S32[i]; // Index of interest
-        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[index]) {
-            continue;
-        }
         j++;
         if (j > high) {
@@ -290 +223 @@
 
 // Mark a pixel for inspection
-static inline void combineInspect(const psArray *inputs, // Stack data
-                                  int x, int y, // Pixel
-                                  int source // Source image index
-                                  )
-{
+// Value in pixel doesn't seem to agree with the stack, so need to look closer
+static inline void combineMarkInspect(const psArray *inputs, // Stack data
+                                      int x, int y, // Pixel
+                                      int source // Source image index
+                                      )
+{
+#ifdef TESTING
+    if (x == TEST_X && y == TEST_Y) {
+        fprintf(stderr, "Marking image %d for inspection\n", source);
+    }
+#endif
     pmStackData *data = inputs->data[source]; // Stack data of interest
     if (!data) {
@@ -304 +243 @@
 }
 
-// Given a stack of images, combine with optional rejection.
-// Pixels in the stack that are rejected are marked for subsequent inspection
-static void combinePixels(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 bad, // Value to give bad pixels
-                          int numIter, // Number of rejection iterations
-                          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)
-                          bool useVariance, // Use variance for rejection when combining?
-                          bool safe,    // Combine safely?
-                          bool rejectInspect, // Reject values marked for inspection from combination?
-                          combineBuffer *buffer // Buffer for combination; to avoid multiple allocations
-                         )
+// 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 combineMarkReject(const psArray *inputs, // Stack data
+                                     int x, int y, // Pixel
+                                     int source // Source image index
+                                     )
+{
+#ifdef TESTING
+    if (x == TEST_X && y == TEST_Y) {
+        fprintf(stderr, "Marking pixel image %d for rejection\n", source);
+    }
+#endif
+    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;
+}
+
+
+// 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
+                           )
 {
     // Rudimentary error checking
+    assert(buffer);
     assert(image);
     assert(mask);
     assert(inputs);
-    assert(numIter >= 0);
-    assert(buffer);
-    assert(addVariance);
-    assert((useVariance && variance) || !useVariance);
 
     psVector *pixelData = buffer->pixels; // Values for the pixel of interest
-    psVector *pixelMasks = buffer->masks; // Masks 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
-    psVector *sort = buffer->sort;      // Sort buffer
+    psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect?
+
+    if (suspect) {
+        *suspect = false;
+    }
 
     // Extract the pixel and mask data
-    int num = 0;                        // Number of good images
+    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
@@ -364 +320 @@
         }
 
+        pixelSuspects->data.U8[numGood] = mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect ?
+            true : false;
+
         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];
+        pixelData->data.F32[numGood] = 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;
+            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 = num;
-    }
-    pixelWeights->n = num;
-    pixelSources->n = num;
-    pixelLimits->n = num;
+        pixelVariances->n = numGood;
+    }
+    pixelWeights->n = numGood;
+    pixelSources->n = numGood;
+    pixelLimits->n = numGood;
+    pixelSuspects->n = numGood;
+    *num = numGood;
 
 #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",
+        for (int i = 0; i < numGood; i++) {
+            fprintf(stderr, "Input %d (%" PRIu16 "): %f %f (%f) %f %d\n",
                     i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
-                    pixelWeights->data.F32[i]);
-        }
-    }
-#endif
-
-    // The sensible thing to do varies according to how many good pixels there are.
+                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[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 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
@@ -427 +407 @@
       }
       case 2: {
-          // Accept the mean of the pixels only if we're going to reject based on the variance, or we're not
-          // playing safe
-          if (useVariance || !safe) {
+          // 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)) {
@@ -437 +416 @@
 #ifdef TESTING
                   if (x == TEST_X && y == TEST_Y) {
-                      fprintf(stderr, "Two inputs to combine using variance/unsafe --> %f %f\n",
-                              mean, var);
+                      fprintf(stderr, "Two inputs to combine using unsafe --> %f %f\n", mean, var);
                   }
 #endif
               }
           }
-          if (useVariance && numIter > 0) {
-              // 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: mark both for inspection
-                  if (rejectInspect) {
-                      imageValue = NAN;
-                      varianceValue = NAN;
-                      maskValue = bad;
-                  } else {
-                      combineInspect(inputs, x, y, pixelSources->data.U16[0]);
-                      combineInspect(inputs, x, y, pixelSources->data.U16[1]);
-                  }
-#ifdef TESTING
-                  if (x == TEST_X && y == TEST_Y) {
-                      fprintf(stderr, "Both pixels marked for inspection (%f > %f x %f\n)",
-                              diff, rej, sqrtf(sigma2));
-                  }
-#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: {
-          // Record the value derived with no clipping, because pixels rejected using the harsh clipping
-          // applied in the first pass might later be accepted.
+          // Can combine without too much worrying
           float mean, var;           // Mean and variance of the combination
           if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
@@ -488 +444 @@
           }
 #endif
-
-          // 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, discard, 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 (!rejectInspect) { \
-        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
-                  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]));
-                          }
-#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);
-                      }
-#endif
-                  }
-              }
-          }
-
-          if (rejectInspect && 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;
       }
@@ -633 +454 @@
     }
 
+#ifdef TESTING
+    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
+#endif
+
+
     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
+                        float iter, // Number of rejection iterations per input
+                        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?
+                        )
+{
+    if (iter <= 0) {
+        return true;
+    }
+
+    int numIter = PS_MAX(iter * num, 1); // Number of iterations
+
+#ifdef TESTING
+    if (x == TEST_X && y == TEST_Y) {
+        fprintf(stderr, "Testing pixel %d,%d: %d %f %f %f %d %d\n",
+                x, y, numIter, rej, sys, olympic, useVariance, safe);
+    }
+#endif
+
+    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
+    psVector *pixelWeights = buffer->weights; // Is the pixel suspect?
+    psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest
+    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
+    psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect?
+    psVector *pixelLimits = buffer->limits; // Is the pixel suspect?
+
+    // Set up rejection limits
+    float rej2 = PS_SQR(rej); // Rejection level squared
+    if (num > 2 && useVariance) {
+        // Convert rejection limits --- saves doing it later multiple times
+        // Using squared rejection limit because it's cheaper than sqrts
+        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 maskIndex = 0;                  // Index of pixel to mask
+    int totalClipped = 0;               // Total number of pixels clipped
+    for (int i = 0; i < numIter && maskIndex >= 0; i++) {
+        maskIndex = -1;                 // Nothing to reject
+
+        switch (num) {
+          case 0:
+            break;
+          case 1:
+            if (i == 0 && 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) > rej2 * sigma2) {
+                      // Not consistent: don't believe either!
+                      if (i == 0 && 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 (i == 0 && 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;
+          }
+#if 0
+          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]);
+              // Errors in pixel differences
+              float err01 = var0 + var1;
+              float err12 = var1 + var2;
+              float err20 = var2 + var0;
+
+#ifdef TESTING
+              if (x == TEST_X && y == TEST_Y) {
+                  fprintf(stderr, "Diff 0-1: %f %f\n", diff01, err01);
+                  fprintf(stderr, "Diff 1-2: %f %f\n", diff12, err12);
+                  fprintf(stderr, "Diff 2-0: %f %f\n", diff20, err20);
+              }
+#endif
+
+              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;
+          }
+#endif
+          default: {
+              if (useVariance) {
+                  float median = combinationWeightedOlympic(pixelData, pixelWeights,
+                                                            olympic, buffer->sort); // Median for stack
+#ifdef TESTING
+                  if (x == TEST_X && y == TEST_Y) {
+                      fprintf(stderr, "Flag with variance, median = %f\n", median);
+                  }
+#endif
+                  float worst = -INFINITY; // Largest deviation
+                  for (int j = 0; j < num; j++) {
+                      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
+
+                      // Comparing squares --- cheaper than lots of sqrts
+                      // pixelVariances includes the rejection limit, from above
+                      float diff2 = PS_SQR(diff); // Square difference
+                      if (diff2 > pixelLimits->data.F32[j]) {
+                          float dev = diff2 / pixelLimits->data.F32[j]; // Deviation
+                          if (dev > worst) {
+                              worst = dev;
+                              maskIndex = j;
+                          }
+                      }
+                  }
+              } else {
+                  float median = NAN, stdev = NAN;  // Median and stdev of the combination, for rejection
+                  combinationMedianStdev(&median, &stdev, pixelData, buffer->sort);
+                  float limit = rej * stdev; // Rejection limit
+#ifdef TESTING
+                  if (x == TEST_X && y == TEST_Y) {
+                      fprintf(stderr, "Flag without variance; median = %f, stdev = %f, limit = %f\n",
+                              median, stdev, limit);
+                  }
+#endif
+                  float worst = -INFINITY; // Largest deviation
+                  for (int j = 0; j < num; j++) {
+                      float diff = fabsf(pixelData->data.F32[j] - median); // Difference from expected
+
+                      if (diff > limit) {
+                          float dev = diff / limit; // Deviation
+                          if (dev > worst) {
+                              worst = dev;
+                              maskIndex = j;
+                          }
+                      }
+                  }
+              }
+          }
+        }
+
+        // Do the actual rejection of the pixel
+        if (maskIndex >= 0) {
+            if (suspect) {
+#ifdef TESTING
+                if (x == TEST_X && y == TEST_Y) {
+                    fprintf(stderr, "Throwing out all suspect pixels\n");
+                }
+#endif
+                // Throw out all suspect pixels
+                int numGood = 0;        // Number of good pixels
+                for (int j = 0; j < num; j++) {
+                    if (pixelSuspects->data.U8[j]) {
+                        combineMarkReject(inputs, x, y, pixelSources->data.U16[j]);
+                        continue;
+                    }
+                    if (numGood == j) {
+                        numGood++;
+                        continue;
+                    }
+                    pixelData->data.F32[numGood] = pixelData->data.F32[j];
+                    pixelWeights->data.F32[numGood] = pixelWeights->data.F32[j];
+                    pixelSources->data.U16[numGood] = pixelSources->data.U16[j];
+                    pixelLimits->data.F32[numGood] = pixelLimits->data.F32[j];
+                    pixelVariances->data.F32[numGood] = pixelVariances->data.F32[j];
+                    numGood++;
+                }
+                pixelData->n = numGood;
+                pixelWeights->n = numGood;
+                pixelSources->n = numGood;
+                pixelLimits->n = numGood;
+                pixelVariances->n = numGood;
+                totalClipped += num - numGood;
+                num = numGood;
+                suspect = false;
+            } else {
+                // Throw out masked pixel
+#ifdef TESTING
+                if (x == TEST_X && y == TEST_Y) {
+                    fprintf(stderr, "Throwing out input %d\n", maskIndex);
+                }
+#endif
+                combineMarkInspect(inputs, x, y, pixelSources->data.U16[maskIndex]);
+                int numGood = 0;        // Number of good pixels
+                for (int j = 0; j < num; j++) {
+                    if (j == maskIndex) {
+                        continue;
+                    }
+                    if (numGood == j) {
+                        numGood++;
+                        continue;
+                    }
+                    pixelData->data.F32[numGood] = pixelData->data.F32[j];
+                    pixelWeights->data.F32[numGood] = pixelWeights->data.F32[j];
+                    pixelSources->data.U16[numGood] = pixelSources->data.U16[j];
+                    pixelLimits->data.F32[numGood] = pixelLimits->data.F32[j];
+                    pixelVariances->data.F32[numGood] = pixelVariances->data.F32[j];
+                    numGood++;
+                }
+                pixelData->n = numGood;
+                pixelWeights->n = numGood;
+                pixelSources->n = numGood;
+                pixelLimits->n = numGood;
+                pixelVariances->n = numGood;
+                totalClipped++;
+                num--;
+            }
+        }
+    }
+
+    return true;
 }
 
@@ -639 +796 @@
 // 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
     )
 {
@@ -668 +825 @@
         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
@@ -681 +838 @@
             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.");
@@ -697 +854 @@
             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);
     }
 
@@ -782 +948 @@
 
 /// 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);
+bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect,
+                    psImageMaskType bad, int kernelSize, float iter, 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;
     }
     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
     PS_ASSERT_INT_POSITIVE(bad, false);
-    PS_ASSERT_INT_NONNEGATIVE(numIter, false);
     if (isnan(rej)) {
-        PS_ASSERT_INT_EQUAL(numIter, 0, false);
+        PS_ASSERT_FLOAT_EQUAL(iter, 0, false);
     } else {
+        PS_ASSERT_FLOAT_LARGER_THAN(iter, 0, false);
         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) {
@@ -833 +989 @@
         }
 #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);
+            }
         }
     }
@@ -863 +1025 @@
     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++) {
+#ifdef TESTING
+            if (x == TEST_X && y == TEST_Y) {
+                fprintf(stderr, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n",
+                        x, y, maskVal, bad, iter, rej, sys, olympic, useVariance, safe, rejection);
+            }
+#endif
+            psVector *reject = NULL; // Images to reject for this pixel
+            if (rejection) {
+                reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, x, y);
+#ifdef TESTING
+                if (x == TEST_X && y == TEST_Y) {
+                    fprintf(stderr, "Rejected inputs: ");
+                    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");
+                    }
+                }
+#endif
+            }
+
+            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);
+            combinePixels(combinedImage, combinedMask, combinedVariance, num, buffer, x, y, bad, safe);
+
+            if (iter > 0) {
+                combineTest(num, suspect, input, buffer, x, y, iter, rej, sys, olympic,
+                            useVariance, safe);
+            }
+        }
+    }
+
+    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;

trunk/psModules/src/imcombine/pmStack.h

@@ -44 +44 @@
                     psArray *input,     ///< Input array of pmStackData
                     psImageMaskType maskVal, ///< Mask value of bad pixels
+                    psImageMaskType suspect, ///< Mask value of suspect pixels
                     psImageMaskType bad,     ///< Mask value to give rejected pixels
                     int kernelSize,     ///< Half-size of the convolution kernel
-                    int numIter,        ///< Number of iterations
+                    float iter,         ///< Number of iterations per input
                     float rej,          ///< Rejection limit (standard deviations)
                     float sys,          ///< Relative systematic error
                     float discard,      ///< Fraction of values to discard for Olympic weighted mean
-                    bool entire,        ///< Combine entire image even if rejection lists provided?
                     bool useVariance,   ///< Use variance values for rejection?
                     bool safe,          ///< Play safe with small numbers of input pixels (mask if N <= 2)?

trunk/psModules/src/objects/pmSourceMatch.c

@@ -17 +17 @@
 #define SOURCES_MAX_LEAF 2              // Maximum number of points on a tree leaf
 #define ARRAY_BUFFER 16                 // Buffer for array
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -124 +125 @@
     match->mag = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
     match->magErr = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
+    match->x = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
+    match->y = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32);
     match->image = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32);
     match->index = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32);
@@ -133 +136 @@
 void pmSourceMatchAdd(pmSourceMatch *match, // Match data
                       float mag, float magErr, // Magnitude and error
+                      float x, float y,        // Position
                       int image, // Image index
                       int index // Source index
@@ -141 +145 @@
     match->mag = psVectorExtend(match->mag, match->mag->nalloc, 1);
     match->magErr = psVectorExtend(match->magErr, match->magErr->nalloc, 1);
+    match->x = psVectorExtend(match->x, match->x->nalloc, 1);
+    match->y = psVectorExtend(match->y, match->y->nalloc, 1);
     match->image = psVectorExtend(match->image, match->image->nalloc, 1);
     match->index = psVectorExtend(match->index, match->index->nalloc, 1);
@@ -147 +153 @@
     match->mag->data.F32[num] = mag;
     match->magErr->data.F32[num] = magErr;
+    match->x->data.F32[num] = x;
+    match->y->data.F32[num] = y;
     match->image->data.S32[num] = image;
     match->index->data.S32[num] = index;
@@ -192 +200 @@
                 pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
                 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
-                                 i, indices->data.S32[j]);
+                                 xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
                 matches->data[j] = match;
             }
@@ -212 +220 @@
                 pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
                 pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
-                                 i, indices->data.S32[j]);
+                                 xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
                 matches->data[k] = match;
             }
@@ -238 +246 @@
                     pmSourceMatch *match = matches->data[index]; // Match data
                     pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
-                                     i, indices->data.S32[j]);
+                                     xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
                     numMatch++;
                 } else {
@@ -244 +252 @@
                     pmSourceMatch *match = pmSourceMatchAlloc(); // Match data
                     pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j],
-                                     i, indices->data.S32[j]);
+                                     xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]);
                     xMaster->data.F32[numMaster] = xImage->data.F32[j];
                     yMaster->data.F32[numMaster] = yImage->data.F32[j];
@@ -304 +312 @@
         pmSourceMatch *match = matches->data[i]; // Match of interest
         for (int j = 0; j < match->num && !source; j++) {
-            if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j])) {
+            if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j]) ||
+                !isfinite(match->x->data.F32[j]) || !isfinite(match->y->data.F32[j])) {
                 continue;
             }
@@ -365 +374 @@
         double star = 0.0, starErr = 0.0; // Accumulators for star
         for (int j = 0; j < match->num; j++) {
-            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
                 continue;
             }
@@ -396 +405 @@
         pmSourceMatch *match = matches->data[i]; // Matched stars
         for (int j = 0; j < match->num; j++) {
-            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
                 continue;
             }
@@ -424 +433 @@
         pmSourceMatch *match = matches->data[i]; // Matched stars
         for (int j = 0; j < match->num; j++) {
-            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
                 continue;
             }
@@ -543 +552 @@
         pmSourceMatch *match = matches->data[i]; // Matched stars
         for (int j = 0; j < match->num; j++) {
-            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) {
                 continue;
             }
@@ -564 +573 @@
                 if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) {
                     numRejected++;
-                    match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;
+                    match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_PHOT;
                 }
             }
@@ -627 +636 @@
             return NULL;
         }
-        psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages);
+        psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric\n", numPhoto, numImages);
 
         fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej1, sys1);
-        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100);
+        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100);
 
         chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys1);
@@ -649 +658 @@
             return NULL;
         }
-        psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric", numPhoto, numImages);
+        psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric\n", numPhoto, numImages);
 
         fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej2, sys2);
-        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected", fracRej * 100);
+        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100);
 
         chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys2);
@@ -668 +677 @@
     return trans;
 }
+
+
+// Iterate on the star positions and image shifts
+// Returns the solution chi^2
+static float sourceMatchRelastroIterate(psVector *xShift, psVector *yShift, // Shift for image
+                                        psVector *xStar, psVector *yStar,   // Position for star
+                                        const psArray *matches // Array of matches
+    )
+{
+    psAssert(matches, "Need list of matches");
+
+    int numImages = xShift->n;          // Number of images
+    int numStars = matches->n;          // Number of stars
+
+    psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32,
+             "Need shifts");
+    psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n);
+    psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32,
+             "Need star positions");
+    psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n");
+
+    // Solve the star positions
+    psVectorInit(xStar, NAN);
+    psVectorInit(yStar, NAN);
+    psVector *starMask = psVectorAlloc(numStars, PS_TYPE_U8); // Mask for stars
+    psVectorInit(starMask, 0xFF);
+    int numGoodStars = 0;               // Number of stars with good measurements
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        int numMeasurements = 0;        // Number of unmasked measurements for star
+        double xSum = 0.0, ySum = 0.0;  // Accumulators for star
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
+                continue;
+            }
+            numMeasurements++;
+            int index = match->image->data.U32[j]; // Image index
+
+            xSum += match->x->data.F32[j] - xShift->data.F32[index];
+            ySum += match->y->data.F32[j] - yShift->data.F32[index];
+        }
+        if (numMeasurements > 1) {
+            // It's only a good star (contributing to the chi^2) if there's more than 1 measurement
+            numGoodStars++;
+            xStar->data.F32[i] = xSum / numMeasurements;
+            yStar->data.F32[i] = ySum / numMeasurements;
+            starMask->data.U8[i] = 0;
+        }
+    }
+
+    // Solve for the shifts
+    psVectorInit(xShift, 0.0);
+    psVectorInit(yShift, 0.0);
+    psVector *num = psVectorAlloc(numImages, PS_TYPE_S32);    // Number of stars
+    psVectorInit(num, 0);
+    for (int i = 0; i < numStars; i++) {
+        if (starMask->data.U8[i]) {
+            continue;
+        }
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
+                continue;
+            }
+            int index = match->image->data.U32[j]; // Image index
+
+            xShift->data.F32[index] += match->x->data.F32[j] - xStar->data.F32[i];
+            yShift->data.F32[index] += match->y->data.F32[j] - yStar->data.F32[i];
+            num->data.S32[index]++;
+        }
+    }
+    for (int i = 0; i < numImages; i++) {
+        xShift->data.F32[i] /= num->data.S32[i];
+        yShift->data.F32[i] /= num->data.S32[i];
+        psTrace("psModules.objects", 3, "Shift for image %d: %f,%f\n",
+                i, xShift->data.F32[i], yShift->data.F32[i]);
+    }
+    psFree(num);
+
+    // Once more through to evaluate chi^2
+    float chi2 = 0.0;                   // chi^2 for iteration
+    int dof = 0;                        // Degrees of freedom
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        if (starMask->data.U8[i]) {
+            continue;
+        }
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+                continue;
+            }
+
+            int index = match->image->data.U32[j]; // Image index
+            float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i];
+            float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i];
+
+            chi2 += PS_SQR(dx) + PS_SQR(dy);
+            dof++;
+        }
+    }
+    dof -= numGoodStars + numImages;
+    chi2 /= dof;
+
+    return chi2;
+}
+
+// Reject star measurements
+// Returns the fraction of measurements that were rejected
+static float sourceMatchRelastroReject(const psVector *xShift, const psVector *yShift, // Shifts for each image
+                                       const psVector *xStar, const psVector *yStar, // Positions for each star
+                                       const psArray *matches, // Array of matches
+                                       float chi2,             // chi^2 from fit
+                                       float rej               // Rejection threshold
+                                )
+{
+    psAssert(matches, "Need list of matches");
+
+    int numImages = xShift->n;          // Number of images
+    int numStars = matches->n;          // Number of stars
+
+    psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32,
+             "Need shifts");
+    psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n);
+    psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32,
+             "Need star positions");
+    psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n");
+
+    int numRejected = 0;                // Number rejected
+    int numMeasurements = 0;            // Number of measurements
+
+    float thresh = PS_SQR(rej) * chi2;    // Threshold for rejection
+
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) {
+                continue;
+            }
+            numMeasurements++;
+            int index = match->image->data.U32[j]; // Image index
+
+            float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i];
+            float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i];
+
+            if (PS_SQR(dx) + PS_SQR(dy) > thresh) {
+                numRejected++;
+                match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_ASTRO;
+            }
+        }
+    }
+
+    return (float)numRejected / (float)numMeasurements;
+}
+
+psArray *pmSourceMatchRelastro(const psArray *matches, // Array of matches
+                               int numImages,          // Number of images
+                               float tol, // Relative tolerance for convergence
+                               int iter1, // Number of iterations for pass 1
+                               float rej1, // Limit on rejection between iterations for pass 1
+                               int iter2, // Number of iterations for pass 2
+                               float rej2, // Limit on rejection between iterations for pass 2
+                               float rejLimit // Limit on rejection between iterations
+    )
+{
+    PS_ASSERT_ARRAY_NON_NULL(matches, NULL);
+
+    int numStars = matches->n;          // Number of stars
+    psVector *xShift = psVectorAlloc(numImages, PS_TYPE_F32); // x shift for each image
+    psVector *yShift = psVectorAlloc(numImages, PS_TYPE_F32); // y shift for each image
+    psVectorInit(xShift, 0.0);
+    psVectorInit(yShift, 0.0);
+    psVector *xStar = psVectorAlloc(numStars, PS_TYPE_F32); // x position for each star
+    psVector *yStar = psVectorAlloc(numStars, PS_TYPE_F32); // y position for each star
+
+    float chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches); // chi^2 for solution
+    psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2);
+    float lastChi2 = INFINITY;          // chi^2 on last iteration
+    float fracRej = INFINITY;           // Fraction of measurements rejected
+
+    // In the first passes, the shifts are not well deteremined: use high systematic error and
+    // rejection thresholds
+    for (int i = 0; i < iter1; i++) {
+        fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej1);
+        psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100);
+
+        chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches);
+        psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
+    }
+
+    for (int i = 0; i < iter2 && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) {
+        lastChi2 = chi2;
+
+        fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej2);
+        psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100);
+
+        chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches);
+        psTrace("psModules.objects", 1, "Pass 2: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
+    }
+
+    psFree(xStar);
+    psFree(yStar);
+
+    if (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit) {
+        psWarning("Unable to converge to relphot solution (%f,%f)", (lastChi2 - chi2) / chi2, fracRej);
+    }
+
+    psArray *results = psArrayAlloc(numImages); // Array of results
+    for (int i = 0; i < numImages; i++) {
+        psVector *offset = results->data[i] = psVectorAlloc(2, PS_TYPE_F32); // Offset for image
+        offset->data.F32[0] = xShift->data.F32[i];
+        offset->data.F32[1] = yShift->data.F32[i];
+    }
+    psFree(xShift);
+    psFree(yShift);
+
+    return results;
+}
+

trunk/psModules/src/objects/pmSourceMatch.h

@@ -1 +1 @@
 #ifndef PM_SOURCE_MATCH_H
 #define PM_SOURCE_MATCH_H
+
+#include <pslib.h>
+
+/// Mask values for matched sources
+typedef enum {
+    PM_SOURCE_MATCH_MASK_PHOT = 0x01,   // Source was rejected from photometry fit
+    PM_SOURCE_MATCH_MASK_ASTRO = 0x02,     // Source was rejected from astrometry fit
+} pmSourceMatchMask;
 
 /// Matched sources
@@ -12 +20 @@
     psVector *mag;                      // Magnitudes
     psVector *magErr;                   // Magnitude errors
+    psVector *x, *y;                    // Positions
     psVector *mask;                     // Mask for measurements
 } pmSourceMatch;
@@ -26 +35 @@
     PS_ASSERT_VECTOR_NON_NULL((MATCH)->mag, RVAL); \
     PS_ASSERT_VECTOR_NON_NULL((MATCH)->magErr, RVAL); \
+    PS_ASSERT_VECTOR_NON_NULL((MATCH)->x, RVAL); \
+    PS_ASSERT_VECTOR_NON_NULL((MATCH)->y, RVAL); \
     PS_ASSERT_VECTOR_SIZE((MATCH)->image, (MATCH)->num, RVAL); \
     PS_ASSERT_VECTOR_SIZE((MATCH)->index, (MATCH)->num, RVAL); \
     PS_ASSERT_VECTOR_SIZE((MATCH)->mag, (MATCH)->num, RVAL); \
     PS_ASSERT_VECTOR_SIZE((MATCH)->magErr, (MATCH)->num, RVAL); \
+    PS_ASSERT_VECTOR_SIZE((MATCH)->x, (MATCH)->num, RVAL); \
+    PS_ASSERT_VECTOR_SIZE((MATCH)->y, (MATCH)->num, RVAL); \
 }
 
@@ -38 +51 @@
 void pmSourceMatchAdd(pmSourceMatch *match, // Match data
                       float mag, float magErr, // Magnitude and error
+                      float x, float y,        // Position
                       int image, // Image index
                       int index // Source index
@@ -75 +89 @@
     );
 
+/// Perform relative astrometry to calibrate images
+psArray *pmSourceMatchRelastro(const psArray *matches, // Array of matches
+                               int numImages,          // Number of images
+                               float tol, // Relative tolerance for convergence
+                               int iter1, // Number of iterations for pass 1
+                               float rej1, // Limit on rejection between iterations for pass 1
+                               int iter2, // Number of iterations for pass 2
+                               float rej2, // Limit on rejection between iterations for pass 2
+                               float rejLimit // Limit on rejection between iterations
+    );
+
 #endif