Index: trunk/psModules/src/imcombine/pmStack.c
===================================================================
--- trunk/psModules/src/imcombine/pmStack.c	(revision 16608)
+++ trunk/psModules/src/imcombine/pmStack.c	(revision 16624)
@@ -8,6 +8,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-02-22 20:05:17 $
+ *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-02-23 02:33:59 $
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  *
@@ -34,4 +34,6 @@
 #define MASK_BAD 0x01                   // Mask value for pixels that have formerly been masked as bad
 #define MASK_SUSPECT 0x02               // Mask value for pixels that are suspected bad (by pmStackCombine)
+
+#define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
 
 
@@ -42,4 +44,5 @@
     psVector *masks;                    // Pixel masks
     psVector *weights;                  // Pixel weights
+    psVector *sources;                  // Pixel sources (which image did they come from?)
     psVector *sort;                     // Buffer for sorting (to get a robust estimator of the standard dev)
 } combineBuffer;
@@ -50,4 +53,5 @@
     psFree(buffer->masks);
     psFree(buffer->weights);
+    psFree(buffer->sources);
     psFree(buffer->sort);
     return;
@@ -64,4 +68,5 @@
     buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK);
     buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
+    buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16);
     buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32);
     return buffer;
@@ -110,13 +115,13 @@
 }
 
-// Generate a standard deviation for the combination
+// Return the median and standard deviation for the pixels
 // Not using psVectorStats because it has additional allocations which slow things down
-static bool combinationStdev(float *median, // Median value, to return
-                             float *stdev, // Standard deviation value, to return
-                             const psVector *values, // Values to combine
-                             const psVector *masks, // Mask to apply
-                             psMaskType maskVal, // Value to mask
-                             psVector *sortBuffer // Buffer for sorting
-                             )
+static bool combinationMedianStdev(float *median, // Median value, to return
+                                   float *stdev, // Standard deviation value, to return
+                                   const psVector *values, // Values to combine
+                                   const psVector *masks, // Mask to apply
+                                   psMaskType maskVal, // Value to mask
+                                   psVector *sortBuffer // Buffer for sorting
+                                   )
 {
     assert(values && masks);
@@ -126,7 +131,8 @@
     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->data.PS_TYPE_MASK_DATA[i] & maskVal)) {
+        if (!(masks->data.PS_TYPE_MASK_DATA[i])) {
             sortBuffer->data.F32[num++] = values->data.F32[i];
         }
@@ -136,20 +142,32 @@
         return false;
     }
-    *median = num % 2 ? sortBuffer->data.F32[num / 2] :
-        (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0 ;
-#if 0
-    if (num < NUM_DIRECT_STDEV) {
-#endif
-        // If there are not many values, the direct standard deviation is better
-        double sum = 0.0;
-        for (int i = 0; i < num; i++) {
-            sum += PS_SQR(sortBuffer->data.F32[i] - *median);
-        }
-        *stdev = sqrt(sum / (float)(num - 1));
-#if 0
+
+    if (num == 3) {
+        // Attempt to measure standard deviation with only three values (and one of those possibly corrupted)
+        *median = sortBuffer->data.F32[1];
+        if (stdev) {
+            float diff1 = sortBuffer->data.F32[0] - *median;
+            float diff2 = sortBuffer->data.F32[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 {
-        *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]);
-    }
-#endif
+        *median = num % 2 ? sortBuffer->data.F32[num / 2] :
+            (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0;
+        if (stdev) {
+            if (num <= NUM_DIRECT_STDEV) {
+                // If there are not many values, the direct standard deviation is better
+                double sum = 0.0;
+                for (int i = 0; i < num; i++) {
+                    sum += PS_SQR(sortBuffer->data.F32[i] - *median);
+                }
+                *stdev = sqrt(sum / (float)(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)]);
+            }
+        }
+    }
 
     return true;
@@ -198,7 +216,37 @@
 }
 
+// Common path for a bad combined pixel
+static inline bool combineBad(psImage *image, // Combined image
+                              psImage *mask, // Combined mask
+                              psImage *weight, // Combined variance map
+                              int x, int y, // Coordinates of interest
+                              psMaskType bad // Value for bad pixels
+                              )
+{
+    image->data.F32[y][x] = NAN;
+    if (weight) {
+        weight->data.F32[y][x] = NAN;
+    }
+    mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
+    return false;
+}
+
+// Mark a pixel for inspection
+static inline void combineInspect(const psArray *inputs, // Stack data
+                                  int x, int y, // Pixel
+                                  int source // Source image index
+                                  )
+{
+    pmStackData *data = inputs->data[source]; // Stack data of interest
+    if (!data) {
+        return;
+    }
+    data->pixels = psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y);
+    return;
+}
+
 // Given a stack of images, combine with optional rejection.
 // Pixels in the stack that are rejected are marked for subsequent inspection
-static bool combinePixels(psImage *image, // Combined image, for output
+static void combinePixels(psImage *image, // Combined image, for output
                           psImage *mask, // Combined mask, for output
                           psImage *weight, // Combined variance map, for output
@@ -211,4 +259,6 @@
                           int numIter, // Number of rejection iterations
                           float rej, // Number of standard deviations at which to reject
+                          bool useVariance, // Use variance for rejection when combining?
+                          bool safe,    // Combine safely?
                           combineBuffer *buffer // Buffer for combination; to avoid multiple allocations
                          )
@@ -221,103 +271,152 @@
     assert(rej > 0);
     assert(buffer);
-
-    int num = inputs->n;          // Number of images to combine
+    assert((useVariance && weight) || !useVariance);
 
     psVector *pixelData = buffer->pixels; // Values for the pixel of interest
     psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest
     psVector *pixelWeights = buffer->weights; // Weights for the pixel of interest
+    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
     psVector *sort = buffer->sort;      // Sort buffer
 
     // Extract the pixel and mask data
-    int numBad = 0;                     // Number of bad pixels
-    for (int i = 0; i < num; i++) {
+    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_MASK_DATA[yIn][xIn] & maskVal) {
+            continue;
+        }
+
         psImage *image = data->readout->image; // Image of interest
         psImage *weight = data->readout->weight; // Weight map of interest
-        psImage *mask = data->readout->mask; // Mask of interest
-        pixelData->data.F32[i] = image->data.F32[yIn][xIn];
+        pixelData->data.F32[num] = image->data.F32[yIn][xIn];
         if (weight) {
-            pixelWeights->data.F32[i] = weight->data.F32[yIn][xIn];
-        }
-        pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[yIn][xIn];
-        if (pixelMasks->data.PS_TYPE_MASK_DATA[i] & maskVal) {
-            numBad++;
-        }
-    }
-    if (numBad == num) {
-        // Catch the case where everything is masked
-        image->data.F32[y][x] = NAN;
-        if (weight) {
-            weight->data.F32[y][x] = NAN;
-        }
-        mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
-        return false;
-    }
-
-    if (reject) {
-        for (int i = 0; i < reject->n; i++) {
-            pixelMasks->data.PS_TYPE_MASK_DATA[reject->data.U16[i]] |= maskVal;
-        }
-    }
-
-    // Record the value derived with no clipping, because pixels rejected using the harsh clipping applied in
-    // the first pass might later be accepted.
-    float mean, variance;               // Mean and variance of the combination
-    if (!combinationMean(&mean, pixelData, weights, pixelMasks, maskVal) ||
-        (weight && !combinationVariance(&variance, pixelWeights, weights, pixelMasks, maskVal))) {
-        image->data.F32[y][x] = NAN;
-            mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
-        if (weight) {
-            weight->data.F32[y][x] = NAN;
-        }
-        return false;
-    }
-
-    image->data.F32[y][x] = mean;
+            pixelWeights->data.F32[num] = weight->data.F32[yIn][xIn];
+        }
+        pixelSources->data.U16[num] = i;
+        num++;
+    }
+    pixelData->n = num;
     if (weight) {
-        weight->data.F32[y][x] = variance;
-    }
-    mask->data.PS_TYPE_MASK_DATA[y][x] = 0;
-
-    // The clipping that follows is solely to identify suspect pixels.
-    // These suspect pixels will be inspected in more detail by other functions.
-    long numClipped = LONG_MAX;         // Number of pixels clipped
-    psMaskType ignore = maskVal | bad;  // Ignore values with this mask value
-    for (int i = 0; i < numIter && numClipped > 0; i++) {
-        float median, stdev;                    // Median and stdev of the combination, for rejection
-        // Only get the mean if it's not the first iteration (because we got the mean earlier)
-        if (!combinationStdev(&median, &stdev, pixelData, pixelMasks, maskVal, sort)) {
-            image->data.F32[y][x] = NAN;
-            mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
-            weight->data.F32[y][x] = NAN;
-            return false;
-        }
-
-        float limit = rej * stdev; // Rejection limit
-        numClipped = 0;
-        for (int j = 0; j < num; j++) {
-            if (pixelMasks->data.PS_TYPE_MASK_DATA[j] & ignore) {
-                continue;
-            }
-            float diff = pixelData->data.F32[j] - median;
-            if (fabsf(diff) > limit) {
-                // Add the pixel as one to inspect
-                pmStackData *data = inputs->data[j]; // Stack data of interest
-                if (!data) {
-                    continue;
-                }
-                data->pixels = psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y);
-                // Mask it so it's not considered in other iterations within this function
-                pixelMasks->data.PS_TYPE_MASK_DATA[j] |= bad;
-                numClipped++;
-            }
-        }
-    }
-
-    return true;
+        pixelWeights->n = num;
+    }
+    pixelSources->n = num;
+
+    // The sensible thing to do varies according to how many good pixels there are.
+    // Default option is that the pixel is bad
+    float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
+    psMaskType maskValue = bad;         // Value for combined mask
+    switch (num) {
+      case 0:
+        // Nothing to combine: it's bad
+        break;
+      case 1: {
+          // Accept the single pixel unless we have to be safe
+          if (!safe) {
+              imageValue = pixelData->data.F32[0];
+              if (weight) {
+                  varianceValue = pixelWeights->data.F32[0];
+              }
+              maskValue = 0;
+          }
+          break;
+      }
+      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) {
+              float mean, variance;   // Mean and variance from combination
+              if (combinationMean(&mean, pixelData, weights, pixelMasks, maskVal) &&
+                  (!weight || combinationVariance(&variance, pixelWeights, weights, pixelMasks, maskVal))) {
+                  imageValue = mean;
+                  varianceValue = variance;
+                  maskValue = 0;
+              }
+          }
+          if (useVariance) {
+              // Use variance to check that the two are consistent
+              float diff = pixelData->data.F32[0] - pixelData->data.F32[1];
+              float sigma = sqrtf(pixelWeights->data.F32[0] + pixelWeights->data.F32[1]);
+              if (fabs(diff) > rej * sigma) {
+                  // Not consistent: mark both for inspection
+                  combineInspect(inputs, x, y, pixelSources->data.U16[0]);
+                  combineInspect(inputs, x, y, pixelSources->data.U16[1]);
+              }
+          }
+          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.
+          float mean, variance;           // Mean and variance of the combination
+          if (!combinationMean(&mean, pixelData, weights, pixelMasks, maskVal) ||
+              (weight && !combinationVariance(&variance, pixelWeights, weights, pixelMasks, maskVal))) {
+              break;
+          }
+          imageValue = mean;
+          varianceValue = variance;
+          maskValue = 0;
+
+          // Prepare for clipping iteration
+          if (numIter > 0) {
+              pixelMasks->n = num;
+              psVectorInit(pixelMasks, 0);
+              if (useVariance) {
+                  // Convert to rejection limits
+                  for (int i = 0; i < num; i++) {
+                      pixelWeights->data.F32[i] = rej * sqrtf(pixelWeights->data.F32[i]);
+                  }
+              }
+          }
+
+          // 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, stdev;    // Median and stdev of the combination, for rejection
+
+              if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev, pixelData, pixelMasks,
+                                          maskVal, sort)) {
+                  psWarning("Bad median/stdev at %d,%d", x, y);
+                  break;
+              }
+
+              float limit;                // Rejection limit, when rejecting based on data properties
+              if (!useVariance) {
+                  limit = rej * stdev;
+              }
+
+              for (int j = 0; j < num; j++) {
+                  if (pixelMasks->data.PS_TYPE_MASK_DATA[j]) {
+                      continue;
+                  }
+                  float diff = fabsf(pixelData->data.F32[j] - median); // Difference from expected
+                  if (diff > (useVariance ? pixelWeights->data.F32[i] : limit)) {
+                      pixelMasks->data.PS_TYPE_MASK_DATA[j] = 0xff;
+                      combineInspect(inputs, x, y, pixelSources->data.U16[j]);
+                      numClipped++;
+                      totalClipped++;
+                  }
+              }
+          }
+          break;
+      }
+    }
+
+    return;
 }
 
@@ -549,5 +648,5 @@
 /// Stack input images
 bool pmStackCombine(pmReadout *combined, psArray *input, psMaskType maskVal, psMaskType bad,
-                    int numIter, float rej)
+                    int numIter, float rej, bool useVariance, bool safe)
 {
     PS_ASSERT_PTR_NON_NULL(combined, false);
@@ -570,4 +669,8 @@
         PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false);
     }
+    if (useVariance && !haveWeights) {
+        psWarning("Unable to use variance in rejection if no variance maps supplied --- option turned off");
+        useVariance = false;
+    }
 
     // Pull out the weightings
@@ -629,5 +732,5 @@
             psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely
             combinePixels(combinedImage, combinedMask, combinedWeight, input, weights, reject, x, y,
-                          maskVal, bad, numIter, rej, buffer);
+                          maskVal, bad, numIter, rej, useVariance, safe, buffer);
         }
         psTrace("psModules.imcombine", 5, "Additional %ld pixels fixed.\n", pixels->n);
@@ -667,5 +770,5 @@
             for (int x = minInputCols; x < maxInputCols; x++) {
                 combinePixels(combinedImage, combinedMask, combinedWeights, input, weights, NULL, x, y,
-                              maskVal, bad, numIter, rej, buffer);
+                              maskVal, bad, numIter, rej, useVariance, safe, buffer);
             }
         }