Index: trunk/psModules/src/imcombine/pmStack.c
===================================================================
--- trunk/psModules/src/imcombine/pmStack.c	(revision 16479)
+++ trunk/psModules/src/imcombine/pmStack.c	(revision 16604)
@@ -8,6 +8,4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-02-14 23:33:09 $
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  *
@@ -23,4 +21,5 @@
 #include "pmHDU.h"
 #include "pmFPA.h"
+#include "pmReadoutStack.h"
 #include "pmConceptsAverage.h"
 
@@ -42,4 +41,8 @@
     psVector *weights;                  // Pixel weights
     psVector *sort;                     // Buffer for sorting (to get a robust estimator of the standard dev)
+#if 0
+    int x0, y0;                         // Offset from original image to combination region
+    int nx, ny;                         // Number of pixels to combine
+#endif
 } combineBuffer;
 
@@ -64,4 +67,11 @@
     buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32);
+
+#if 0
+    buffer->x0 = 0;
+    buffer->y0 = 0;
+    buffer->nx = 0;
+    buffer->ny = 0;
+#endif
 
     return buffer;
@@ -136,7 +146,20 @@
         return false;
     }
-    *median = num % 2 ? (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0 :
-        sortBuffer->data.F32[num / 2];
-    *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]);
+    *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
+    } else {
+        *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]);
+    }
+#endif
 
     return true;
@@ -186,5 +209,5 @@
 
 // Given a stack of images, combine with optional rejection.
-// Pixels in the stack that are rejected are marked for subsequent
+// Pixels in the stack that are rejected are marked for subsequent inspection
 static bool combinePixels(psImage *image, // Combined image, for output
                           psImage *mask, // Combined mask, for output
@@ -193,5 +216,5 @@
                           const psVector *weights, // Global (single value) weights for data, or NULL
                           const psVector *reject, // Indices of pixels to reject, or NULL
-                          int x, int y, // Coordinates of interest
+                          int x, int y, // Coordinates of interest; frame of output image
                           psMaskType maskVal, // Value to mask
                           psMaskType bad, // Value to give bad pixels
@@ -223,12 +246,13 @@
             continue;
         }
+        int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
         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[y][x];
+        pixelData->data.F32[i] = image->data.F32[yIn][xIn];
         if (weight) {
-            pixelWeights->data.F32[i] = weight->data.F32[y][x];
-        }
-        pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[y][x];
+            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++;
@@ -563,4 +587,24 @@
     }
 
+    // Get the sizes
+    psArray *stack = psArrayAlloc(input->n);
+    for (int i = 0; i < stack->n; i++) {
+        pmStackData *data = input->data[i]; // Stack data
+        stack->data[i] = psMemIncrRefCounter(data->readout);
+    }
+    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+    int xSize, ySize;                   // Size of the output image
+    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
+                                stack)) {
+        psError(PS_ERR_UNKNOWN, false, "Input stack is not valid.");
+        psFree(stack);
+        return false;
+    }
+    psFree(stack);
+    pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, bad);
+    psTrace("psModules.imcombine", 1, "Combining [%d:%d,%d:%d] (%dx%d)\n",
+            minInputCols, maxInputCols, minInputRows, maxInputRows, xSize, ySize);
+
+
     // Buffer for combination
     combineBuffer *buffer = combineBufferAlloc(num, numIter == 0 ? PS_STAT_SAMPLE_MEAN :
@@ -573,5 +617,5 @@
         psImage *combinedWeight = combined->weight; // Combined mask
 
-        psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source
+        psArray *pixelMap = pixelMapGenerate(input, maxInputCols, maxInputRows); // Map of pixels to source
         psPixels *pixels = NULL;            // Total list of pixels, with no duplicates
         for (int i = 0; i < num; i++) {
@@ -584,5 +628,9 @@
         }
         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, x, y); // Inspect these images closely
             combinePixels(combinedImage, combinedMask, combinedWeight, input, weights, reject, x, y,
@@ -622,6 +670,6 @@
         }
 
-        for (int y = 0; y < numRows; y++) {
-            for (int x = 0; x < numCols; x++) {
+        for (int y = minInputRows; y < maxInputRows; y++) {
+            for (int x = minInputCols; x < maxInputCols; x++) {
                 combinePixels(combinedImage, combinedMask, combinedWeights, input, weights, NULL, x, y,
                               maskVal, bad, numIter, rej, buffer);