Index: trunk/psModules/src/imcombine/pmImageCombine.c
===================================================================
--- trunk/psModules/src/imcombine/pmImageCombine.c	(revision 10551)
+++ trunk/psModules/src/imcombine/pmImageCombine.c	(revision 11115)
@@ -8,6 +8,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-12-08 11:39:30 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-01-16 23:51:51 $
  *
  *  XXX: pmRejectPixels() has a known bug with the pmImageTransform() call.
@@ -26,17 +26,163 @@
 #include "pslib.h"
 
-//
-// The following macros define how big the initial pixel list will be, and
-// how much it should be incremented when realloc'ed.
-//
-#define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH 100
-#define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC 100
-#define PS_COMBINE_IMAGE_MAX_QUESTIONABLE_PIXELS 1000
-/******************************************************************************
-pmCombineImages(combine, questionablePixels, images, errors, masks, maskVal,
-                pixels, numIter, sigmaClip, stats)
- 
-XXX: Allocate a dummy psStats structure so that we don't destroy away its data.
- *****************************************************************************/
+#define PIXEL_LIST_BUFFER 100           // Size of the pixel list buffer
+
+// Data structure for use as a buffer in combining pixels
+typedef struct
+{
+    psVector *pixels;                   // Pixel values
+    psVector *masks;                    // Pixel masks
+    psVector *errors;                   // Pixel errors
+    psStats *stats;                     // Statistics to use with combination
+}
+combineBuffer;
+
+void combineBufferFree(combineBuffer *buffer)
+{
+    psFree(buffer->pixels);
+    psFree(buffer->masks);
+    psFree(buffer->errors);
+    psFree(buffer->stats);
+}
+
+combineBuffer *combineBufferAlloc(long numImages // Number of images that will be combined
+                                 )
+{
+    combineBuffer *buffer = psAlloc(sizeof(combineBuffer));
+    psMemSetDeallocator(buffer, (psFreeFunc)combineBufferFree);
+
+    buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32);
+    buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK);
+    buffer->errors = psVectorAlloc(numImages, PS_TYPE_F32);
+    buffer->stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+
+    return buffer;
+}
+
+
+static bool combinePixels(psImage *combine, // Combined image, for output
+                          psArray *questionablePixels, // Array of rejection masks
+                          int x, int y, // Position in the images
+                          const psArray *images, // Array of input images
+                          const psArray *errors, // Array of input error images
+                          const psArray *masks, // Array of input masks
+                          psU32 maskVal, // Mask value
+                          psS32 numIter, // Number of rejection iterations
+                          psF32 sigmaClip, // Number of standard deviations at which to reject
+                          combineBuffer *buffer // Buffer for combination; to avoid multiple allocations
+                         )
+{
+    assert(combine);
+    assert(x >= 0 && x < combine->numCols);
+    assert(y >= 0 && y < combine->numRows);
+    assert(images);
+    int numImages = images->n;          // Number of images to combine
+    if (masks) {
+        assert(masks->n == numImages);
+    }
+    if (errors) {
+        assert(errors->n == numImages);
+    }
+    assert(numIter >= 0);
+    assert(sigmaClip > 0);
+    assert(!questionablePixels || (questionablePixels && questionablePixels->n == numImages));
+
+    if (buffer) {
+        psMemIncrRefCounter(buffer);
+    } else {
+        buffer = combineBufferAlloc(numImages);
+    }
+
+    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
+    psVector *pixelErrors = buffer->errors; // Errors for the pixel of interest
+    psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest
+    psStats *stats = buffer->stats;     // Statistics for combination
+
+    //
+    // Loop through each image, extract the pixel/mask/error data into psVectors.
+    //
+    if (!masks) {
+        psVectorInit(pixelMasks, 0);
+    }
+    if (!errors) {
+        pixelErrors = NULL;
+    }
+    for (int i = 0; i < numImages; i++) {
+        // Set the pixel data
+        psImage *image = images->data[i]; // Image of interest
+        pixelData->data.F32[i] = image->data.F32[y][x];
+        // Set the pixel mask data, if necessary
+        if (masks) {
+            psImage *mask = masks->data[i]; // Mask of interest
+            pixelMasks->data.U8[i] = mask->data.U8[y][x];
+        }        // Set the pixel error data, if necessary
+        if (errors) {
+            psImage *error = errors->data[i]; // Error image of interest
+            pixelErrors->data.F32[i] = error->data.F32[y][x];
+        }
+    }
+
+    //
+    // Iterate on the pixels, rejecting outliers
+    //
+    for (int iter = 0; iter < numIter; iter++) {
+        // Combine all the pixels, using the specified stat.
+        if (!psVectorStats(stats, pixelData, pixelErrors, pixelMasks, maskVal)) {
+            combine->data.F32[y][x] = NAN;
+            psFree(buffer);
+            return false;
+        }
+        float combinedPixel = stats->sampleMean; // Value of the combination
+
+        if (iter == 0) {
+            // Save the value produced with no rejection, since it may be useful later
+            // (if the rejection turns out to be unnecessary)
+            combine->data.F32[y][x] = combinedPixel;
+        }
+
+        //
+        // Reject all pixels that lie more that sigmaClip standard deviations from
+        // the combined pixel value.
+        //
+        int numRejects = 0;     // Number of rejections
+        float stdev = stats->sampleStdev;
+        for (int i = 0; i < numImages; i++) {
+            if (!(pixelMasks->data.U8[i] & maskVal) &&
+                    fabs(pixelData->data.F32[i] - combinedPixel) > sigmaClip * stdev) {
+                // Reject pixel as questionable
+                numRejects++;
+                pixelMasks->data.U8[i] = maskVal;
+                if (questionablePixels) {
+                    // Mark the pixel as questionable
+                    psPixels *qp = questionablePixels->data[i]; // Questionable pixels for this image
+                    int qpNum = qp->n; // Number of QPs in the image of interest
+                    if (qpNum >= qp->nalloc) {
+                        // Grow dynamically, if required
+                        qp = psPixelsRealloc(qp, qp->nalloc + PIXEL_LIST_BUFFER);
+                        questionablePixels->data[i] = qp;
+                    }
+                    qp->data[qpNum].x = x;
+                    qp->data[qpNum].y = y;
+                    qp->n++;
+                }
+            }
+        }
+
+        //
+        // If the number of rejected pixels is zero, then there's no point continuing the loop.
+        //
+        if (numRejects == 0) {
+            break;
+        }
+
+        //
+        // XXX: Is it possible to have all pixels rejected?  If so, we should exit the loop.
+        //
+    }
+
+    psFree(buffer);
+    return true;
+}
+
 psImage *pmCombineImages(
     psImage *combine,                   ///< Combined image (output)
@@ -48,69 +194,53 @@
     const psPixels *pixels,             ///< Pixels to combine
     psS32 numIter,                      ///< Number of rejection iterations
-    psF32 sigmaClip,                    ///< Number of standard deviations at which to reject
-    psStats *stats)               ///< Statistics to use in the combination
-{
-
-    PS_ASSERT_PTR_NON_NULL(images, combine);
-    psU32 numImages = images->n;
+    psF32 sigmaClip                     ///< Number of standard deviations at which to reject
+)
+{
     psTrace("psModules.imcombine", 3, "Calling pmCombineImages(%ld)\n", images->n);
 
-    if (errors != NULL) {
-        if (images->n != errors->n) {
-            psError(PS_ERR_UNKNOWN, true, "images and errors args must have same length (%ld != %ld)\n",
-                    images->n, errors->n);
-            return(combine);
-        }
-    }
-    if (masks != NULL) {
-        if (images->n != masks->n) {
-            psError(PS_ERR_UNKNOWN, true, "images and masks args must have same length (%ld != %ld)\n",
-                    images->n, masks->n);
-            return(combine);
-        }
-    }
-
-    psImage *tmpImg = (psImage *) images->data[0];
-    psU32 numRows = tmpImg->numRows;
-    psU32 numCols = tmpImg->numCols;
-
-    //
-    // Check that all images have the appropriate size and type.
-    //
-    for (psS32 i = 0 ; i < numImages ; i++) {
-        psImage *tmpDataImg = (psImage *) images->data[i];
-        PS_ASSERT_IMAGE_NON_NULL(tmpDataImg, combine);
-        PS_ASSERT_IMAGE_TYPE(tmpDataImg, PS_TYPE_F32, combine);
-        if ((tmpDataImg->numRows != numRows) || (tmpDataImg->numCols != numCols)) {
-            psError(PS_ERR_UNKNOWN, true, "image %d has size (%d, %d); should be (%d, %d).\n",
-                    i, tmpDataImg->numRows, tmpDataImg->numCols, numRows, numCols);
-        }
-
-        if (errors != NULL) {
-            psImage *tmpErrorImg = (psImage *) errors->data[i];
-            PS_ASSERT_IMAGE_NON_NULL(tmpErrorImg, combine);
-            PS_ASSERT_IMAGE_TYPE(tmpErrorImg, PS_TYPE_F32, combine);
-            PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpErrorImg, NULL);
-        }
-
-        if (masks != NULL) {
-            psImage *tmpMaskImg = (psImage *) masks->data[i];
-            PS_ASSERT_IMAGE_NON_NULL(tmpMaskImg, combine);
-            PS_ASSERT_IMAGE_TYPE(tmpMaskImg, PS_TYPE_U8, combine);
-            PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpMaskImg, NULL);
-        }
-    }
-    PS_ASSERT_PTR_NON_NULL(stats, combine);
-    psStatsOptions statistic = psStatsSingleOption(stats->options);
-    if (statistic == 0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %x\n", statistic);
-        return combine;
-    }
+    PS_ASSERT_ARRAY_NON_NULL(images, NULL);
+    PS_ASSERT_INT_POSITIVE(images->n, NULL);
+    long numImages = images->n;          // Number of images
+    int numCols = ((psImage*)images->data[0])->numCols; // Size in x
+    int numRows = ((psImage*)images->data[0])->numRows; // Size in y
+
+    if (combine) {
+        PS_ASSERT_IMAGE_NON_NULL(combine, NULL);
+        PS_ASSERT_IMAGE_SIZE(combine, numCols, numRows, NULL);
+    }
+    if (questionablePixels && !*questionablePixels) {
+        PS_ASSERT_ARRAY_NON_NULL(*questionablePixels, NULL);
+        PS_ASSERT_ARRAY_SIZE(*questionablePixels, numImages, 0);
+    }
+    for (int i = 1; i < images->n; i++) {
+        psImage *image = images->data[i]; // Image of interest
+        PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
+        PS_ASSERT_IMAGE_SIZE(image, numCols, numRows, NULL);
+    }
+    if (errors) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(images, errors, NULL);
+        for (int i = 0; i < images->n; i++) {
+            psImage *error = errors->data[i];
+            PS_ASSERT_IMAGE_SIZE(error, numCols, numRows, NULL);
+            PS_ASSERT_IMAGE_TYPE(error, PS_TYPE_F32, NULL);
+        }
+    }
+    if (masks) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(images, masks, NULL);
+        for (int i = 0; i < images->n; i++) {
+            psImage *mask  = masks->data[i];
+            PS_ASSERT_IMAGE_SIZE(mask, numCols, numRows, NULL);
+            PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL);
+        }
+    }
+    PS_ASSERT_INT_POSITIVE(numIter, NULL);
+    PS_ASSERT_FLOAT_LARGER_THAN(sigmaClip, 0.0, NULL);
 
     // Allocate and initialize the combined image, if necessary.
-    if (combine == NULL) {
+    if (!combine) {
         combine = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        if (pixels != NULL) {
-            psImageInit(combine, 0.0);
+        if (pixels) {
+            // Set everything we're not working on to NAN
+            psImageInit(combine, NAN);
         }
     }
@@ -120,140 +250,31 @@
     // struct for each image.
     //
-    if (*questionablePixels == NULL) {
-        *questionablePixels = psArrayAlloc(numImages);
-    } else if ((*questionablePixels)->n != numImages) {
-        *questionablePixels = psArrayRealloc(*questionablePixels, numImages);
-    }
-    for (psS32 im = 0 ; im < numImages ; im++) {
-        psFree((*questionablePixels)->data[im]);
-        ((*questionablePixels)->data[im]) =
-            (psPtr *) psPixelsAlloc(PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH);
-        ((psPixels *) ((*questionablePixels)->data[im]))->n = 0;
-    }
-    //
-    // qpPtr is used to maintain a count of the questionable pixels for each image.
-    //
-    psVector *qpPtr = psVectorAlloc(numImages, PS_TYPE_S32);
-    psVectorInit(qpPtr, 0);
-    //
-    // Allocate the necessary psVectors for the call to psVectorStats().
-    // These vectors will be used whether we are combining a list of pixels,
-    // or every pixel in the input images.
-    //
-    psVector *pixelData = psVectorAlloc(numImages, PS_TYPE_F32);
-    psVector *pixelMask = NULL;
-    if (masks != NULL) {
-        pixelMask = psVectorAlloc(numImages, PS_TYPE_U8);
-        psVectorInit(pixelMask, 0);
-    }
-
-    psVector *pixelErrors = NULL;
-    if (errors != NULL) {
-        pixelErrors = psVectorAlloc(numImages, PS_TYPE_F32);
-        psVectorInit(pixelErrors, 1.0);
-    }
-
-    if (pixels != NULL) {
+    if (questionablePixels) {
+        if (*questionablePixels == NULL) {
+            *questionablePixels = psArrayAlloc(numImages);
+        } else if ((*questionablePixels)->n != numImages) {
+            *questionablePixels = psArrayRealloc(*questionablePixels, numImages);
+        }
+        for (int i = 0; i < numImages; i++) {
+            psFree((*questionablePixels)->data[i]);
+            (*questionablePixels)->data[i] = psPixelsAlloc(PIXEL_LIST_BUFFER);
+        }
+    }
+
+    combineBuffer *buffer = combineBufferAlloc(numImages); // Buffer for combination
+
+    if (pixels) {
         // Only those specified pixels should be combined.
 
-        psStats *stdevStats = psStatsAlloc(PS_STAT_SAMPLE_STDEV);
-
-        for (psS32 p = 0 ; p < pixels->n ; p++) {
-            // Must initialize the mask to 0 for every pixel.
-            psVectorInit(pixelMask, 0);
-            psS32 col = (pixels->data[p]).x;
-            psS32 row = (pixels->data[p]).y;
-
-            //
-            // Loop through each image, extract the pixel/mask/error data
-            // into psVectors.
-            //
-            for (psS32 im = 0 ; im < numImages ; im++) {
-                // Set the pixel data
-                pixelData->data.F32[im] =       ((psImage *) images->data[im])->data.F32[row][col];
-                // Set the pixel mask data, if necessary
-                if (masks != NULL) {
-                    pixelMask->data.U8[im] =   ((psImage *) masks->data[im])->data.F32[row][col];
-                }
-
-                // Set the pixel error data, if necessary
-                if (errors != NULL) {
-                    pixelErrors->data.F32[im] = ((psImage *) errors->data[im])->data.F32[row][col];
-                }
+        for (int p = 0; p < pixels->n; p++) {
+            int x = pixels->data[p].x; // Column of interest
+            int y = pixels->data[p].y; // Row of interest
+
+            if (!combinePixels(combine, questionablePixels ? *questionablePixels : NULL, x, y,
+                               images, errors, masks, maskVal, numIter, sigmaClip, buffer)) {
+                // Bad pixel --- no big deal
+                psErrorClear();
             }
-
-            //
-            // Iterate on the pixels, rejecting outliers
-            //
-            for (psS32 iter = 0 ; iter < numIter ; iter++) {
-                // Combine all the pixels, using the specified stat.
-
-                psVectorStats(stats, pixelData, pixelErrors, pixelMask, maskVal);
-                psF64 combinedPixel = psStatsGetValue(stats, statistic);
-
-                if (iter == 0) {
-                    combine->data.F32[row][col] = (psF32) combinedPixel;
-                }
-
-                //
-                // Reject all pixels that lie more that sigmaClip standard deviations from
-                // the combined pixel value.
-                //
-                psS32 numRejects = 0;
-                for (psS32 im = 0 ; im < numImages ; im++) {
-                    psVectorStats(stdevStats, pixelData, pixelErrors, pixelMask, maskVal);
-                    psF64 stdev = stdevStats->sampleStdev;
-
-                    if (!(pixelMask->data.U8[im] & maskVal)) {
-                        if (fabs(pixelData->data.F32[im]) >
-                                (fabs(sigmaClip * stdev) + fabs(combinedPixel))) {
-                            numRejects++;
-                            pixelMask->data.U8[im] = maskVal;
-                            //
-                            // XXX: These data structures indirections are getting complicated.
-                            //
-                            psS32 ptr = qpPtr->data.S32[im];
-                            psPixels *pixelListPtr = ((psPixels *) ((*questionablePixels)->data[im]));
-
-                            if (ptr >= pixelListPtr->nalloc) {
-                                (*questionablePixels)->data[im] =
-                                    (psPtr *) psPixelsRealloc(((psPixels *) ((*questionablePixels)->data[im])),
-                                                              ((((psPixels *) ((*questionablePixels)->data[im]))->nalloc) +
-                                                               PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC));
-                                // XXX: Can the realloc() fail?  Must we check for NULL?
-                            }
-                            ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].x = col;
-                            ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].y = row;
-                            (qpPtr->data.S32[im])++;
-                            // XXX: this pixel ->n increment is wierd
-                            ((psPixels *) ((*questionablePixels)->data[im]))->n = qpPtr->data.S32[im];
-                        }
-                    }
-                }
-
-                //
-                // If the number of rejected pixels is zero, then there's no point
-                // continuing the loop.
-                //
-                if (numRejects == 0) {
-                    break;
-                }
-                psS32 totalRejects = 0;
-                for (psS32 im = 0 ; im < numImages ; im++) {
-                    if (pixelMask->data.U8[im] & maskVal) {
-                        totalRejects++;
-                    }
-                }
-
-                //
-                // XXX: Is it possible to have all pixels rejected?  If so, we should
-                // exit the loop.
-                //
-                if (totalRejects == numImages) {
-                    break;
-                }
-            }
-        }
-        psFree(stdevStats);
+        }
     } else {
         //
@@ -267,35 +288,19 @@
         // combine image.
         //
-        for (psS32 row = 0 ; row < numRows ; row++) {
-            for (psS32 col = 0 ; col < numCols ; col++) {
-                for (psS32 im = 0 ; im < numImages ; im++) {
-                    // Set the pixel data
-                    pixelData->data.F32[im] =       ((psImage *) images->data[im])->data.F32[row][col];
-
-                    // Set the pixel mask data, if necessary
-                    if (masks != NULL) {
-                        pixelMask->data.U8[im] =   ((psImage *) masks->data[im])->data.F32[row][col];
-                    }
-
-                    // Set the pixel error data, if necessary
-                    if (errors != NULL) {
-                        pixelErrors->data.F32[im] = ((psImage *) errors->data[im])->data.F32[row][col];
-                    }
-
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                if (!combinePixels(combine, questionablePixels ? *questionablePixels : NULL, x, y,
+                                   images, errors, masks, maskVal, numIter, sigmaClip, buffer)) {
+                    // Bad pixel --- no big deal
+                    psErrorClear();
                 }
-                // Combine all the pixels, using the specified stat.
-                psVectorStats(stats, pixelData, pixelErrors, pixelMask, maskVal);
-                combine->data.F32[row][col] = psStatsGetValue(stats, statistic);
             }
         }
     }
 
-    psFree(pixelData);
-    psFree(pixelMask);
-    psFree(pixelErrors);
-    psFree(qpPtr);
+    psFree(buffer);
 
     psTrace("psModules.imcombine", 3, "Exiting pmCombineImages(%ld)\n", images->n);
-    return(combine);
+    return combine;
 }