Index: /branches/eam_branch_20080719/psModules/src/imcombine/pmReadoutCombine.c
===================================================================
--- /branches/eam_branch_20080719/psModules/src/imcombine/pmReadoutCombine.c	(revision 18810)
+++ /branches/eam_branch_20080719/psModules/src/imcombine/pmReadoutCombine.c	(revision 18811)
@@ -42,4 +42,122 @@
 }
 
+// check the input parameters and set up the output images
+bool pmReadoutCombinePrepare(pmReadout *output, const psArray *inputs, const pmCombineParams *params)
+{
+    // Check inputs
+    PS_ASSERT_PTR_NON_NULL(output, false);
+    PS_ASSERT_ARRAY_NON_NULL(inputs, false);
+    PS_ASSERT_PTR_NON_NULL(params, false);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracLow, 0.0, 1.0, false);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracHigh, 0.0, 1.0, false);
+
+    // valid combintion statistic?
+    bool valid = false;
+    valid |= (params->combine == PS_STAT_SAMPLE_MEAN);
+    valid |= (params->combine == PS_STAT_SAMPLE_MEDIAN);
+    valid |= (params->combine == PS_STAT_ROBUST_MEDIAN);
+    valid |= (params->combine == PS_STAT_FITTED_MEAN);
+    valid |= (params->combine == PS_STAT_CLIPPED_MEAN);
+    if (!valid) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Combination method is not SAMPLE_MEAN, SAMPLE_MEDIAN, "
+                "ROBUST_MEDIAN, FITTED_MEAN or CLIPPED_MEAN.\n");
+        return false;
+    }
+
+    // weights exist if weights desired?
+    for (int i = 0; i < inputs->n; i++) {
+        pmReadout *readout = inputs->data[i]; // Readout of interest
+        if (params->weights && !readout->weight) {
+            psError(PS_ERR_UNEXPECTED_NULL, true,
+                    "Rejection based on weights requested, but no weights supplied for image %d.\n", i);
+            return false;
+        }
+    }
+
+    pmHDU *hdu = pmHDUFromReadout(output); // Output HDU
+    if (!hdu) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find HDU for readout.\n");
+        return false;
+    }
+
+    //  set the output header metadata
+    psString comment = NULL;        // Comment to add to header
+    psStringAppend(&comment, "Combining using statistic: %x", params->combine);
+    if (!hdu->header) {
+	hdu->header = psMetadataAlloc();
+    }
+    psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+    psFree(comment);
+
+    // note the clipping parameters, if used
+    if (params->combine == PS_STAT_CLIPPED_MEAN) {
+	psString comment = NULL;    // Comment to add to header
+	psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", params->iter, params->rej);
+	psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+	psFree(comment);
+    }
+
+    // note the use of weights
+    if (params->weights) {
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
+                         "Using input weights to combine images", "");
+    }
+
+    // note the rejection fraction
+    float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep
+    if (keepFrac != 1.0) {
+        psString comment = NULL;        // Comment to add to header
+        psStringAppend(&comment, "Min/max rejection: %f high, %f low, keep %d",
+                       params->fracHigh, params->fracLow, params->nKeep);
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+        psFree(comment);
+    }
+
+    // note the mask value actually used
+    psMaskType maskVal = params->maskVal; // The mask value
+    if (maskVal) {
+        psString comment = NULL;        // Comment to add to header
+        psStringAppend(&comment, "Mask for combination: %x", maskVal);
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+        psFree(comment);
+    }
+
+    // determine the output image size based on the input images
+    int row0, col0, numCols, numRows;
+    if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, inputs)) {
+        psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
+        return false;
+    }
+
+    // generate the required output images based on the specified sizes
+    pmReadoutStackDefineOutput(output, col0, row0, numCols, numRows, true, params->weights, params->blank);
+    psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
+
+    // these calls allocate and save the requested images on the output analysis metadata
+    psImage *counts = pmReadoutSetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT, numCols, numRows, PS_TYPE_U16, 0);
+    if (!counts) {
+        return false;
+    }
+    psImage *sigma = pmReadoutSetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA, numCols, numRows, PS_TYPE_F32, NAN);
+    if (!sigma) {
+        return false;
+    }
+
+    // Update the "concepts"
+    psList *inputCells = psListAlloc(NULL); // List of cells
+    for (long i = 0; i < inputs->n; i++) {
+        pmReadout *readout = inputs->data[i]; // Readout of interest
+        psListAdd(inputCells, PS_LIST_TAIL, readout->parent);
+    }
+    bool success = pmConceptsAverageCells(output->parent, inputCells, NULL, NULL, true);
+    psFree(inputCells);
+
+    // set these even though the values are not yet set
+    output->data_exists = true;
+    output->parent->data_exists = true;
+    output->parent->parent->data_exists = true;
+
+    return success;
+}
 
 // XXX: Maybe add support for S16 and S32 types.  Currently, only F32 supported.
@@ -61,21 +179,4 @@
     PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracLow, 0.0, 1.0, false);
     PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracHigh, 0.0, 1.0, false);
-    if (params->combine != PS_STAT_SAMPLE_MEAN && params->combine != PS_STAT_SAMPLE_MEDIAN &&
-            params->combine != PS_STAT_ROBUST_MEDIAN && params->combine != PS_STAT_FITTED_MEAN &&
-            params->combine != PS_STAT_CLIPPED_MEAN) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Combination method is not SAMPLE_MEAN, SAMPLE_MEDIAN, "
-                "ROBUST_MEDIAN, FITTED_MEAN or CLIPPED_MEAN.\n");
-        return false;
-    }
-    for (int i = 0; i < inputs->n; i++) {
-        pmReadout *readout = inputs->data[i]; // Readout of interest
-        if (params->weights && !readout->weight) {
-            psError(PS_ERR_UNEXPECTED_NULL, true,
-                    "Rejection based on weights requested, but no weights supplied for image %d.\n", i);
-            return false;
-        }
-    }
-
-    bool first = !output->image;        // First pass through?
 
     pmHDU *hdu = pmHDUFromReadout(output); // Output HDU
@@ -85,13 +186,8 @@
     }
 
-    if (first) {
-        psString comment = NULL;        // Comment to add to header
-        psStringAppend(&comment, "Combining using statistic: %x", params->combine);
-        if (!hdu->header) {
-            hdu->header = psMetadataAlloc();
-        }
-        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
-        psFree(comment);
-    }
+    pthread_t id = pthread_self();
+    char name[64];
+    sprintf (name, "%x", (unsigned int) id);
+    psTimerStart (name);
 
     psStatsOptions combineStdev = 0; // Statistics option for weights
@@ -118,43 +214,23 @@
         stats->clipSigma = params->rej;
         stats->clipIter = params->iter;
-
-        if (first) {
-            psString comment = NULL;    // Comment to add to header
-            psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma",
-                           params->iter, params->rej);
-            psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
-            psFree(comment);
-        }
-    }
-    if (params->weights && first) {
-        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
-                         "Using input weights to combine images", "");
-    }
+    }
+
+    psImage *counts = pmReadoutGetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT);
+    if (!counts) {
+        return false;
+    }
+    psImage *sigma = pmReadoutGetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA);
+    if (!sigma) {
+        return false;
+    }
+
+    stats->options |= combineStdev;
 
     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,
-                                inputs)) {
+    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, inputs)) {
         psError(PS_ERR_UNKNOWN, false, "No valid input readouts.");
         return false;
     }
-
-    pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true, params->weights,
-                        params->blank);
-    psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
-
-    psImage *counts = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT, xSize, ySize,
-                                             PS_TYPE_U16, 0);
-    if (!counts) {
-        return false;
-    }
-    psImage *sigma = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA, xSize, ySize,
-                                            PS_TYPE_F32, NAN);
-    if (!sigma) {
-        psFree(counts);
-        return false;
-    }
-
-    stats->options |= combineStdev;
 
     // We loop through each pixel in the output image.  We loop through each input readout.  We determine if
@@ -180,19 +256,5 @@
 
     float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep
-    if (keepFrac != 1.0 && first) {
-        psString comment = NULL;        // Comment to add to header
-        psStringAppend(&comment, "Min/max rejection: %f high, %f low, keep %d",
-                       params->fracHigh, params->fracLow, params->nKeep);
-        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
-        psFree(comment);
-    }
-
     psMaskType maskVal = params->maskVal; // The mask value
-    if (maskVal && first) {
-        psString comment = NULL;        // Comment to add to header
-        psStringAppend(&comment, "Mask for combination: %x", maskVal);
-        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
-        psFree(comment);
-    }
 
     #ifndef PS_NO_TRACE
@@ -226,5 +288,5 @@
         int yOut = i - output->row0; // y position on output readout
         #ifdef SHOW_BUSY
-
+	
         if (psTraceGetLevel("psModules.imcombine") > 9) {
             printf("Processing row %d\r", i);
@@ -242,11 +304,4 @@
                 int xIn = j - readout->col0; // x position on input readout
                 psImage *image = readout->image; // The readout image
-
-                #if 0 // This should have been taken care of already:
-                // Check bounds
-                if (xIn < 0 || xIn >= image->numCols || yIn < 0 || yIn >= image->numRows) {
-                    continue;
-                }
-                #endif
 
                 pixelsData[r] = image->data.F32[yIn][xIn];
@@ -348,4 +403,5 @@
     }
     #endif
+
     psFree(index);
     psFree(pixels);
@@ -355,21 +411,7 @@
     psFree(stats);
     psFree(invScale);
-    psFree(counts);
-    psFree(sigma);
-
-    // Update the "concepts"
-    psList *inputCells = psListAlloc(NULL); // List of cells
-    for (long i = 0; i < inputs->n; i++) {
-        pmReadout *readout = inputs->data[i]; // Readout of interest
-        psListAdd(inputCells, PS_LIST_TAIL, readout->parent);
-    }
-    bool success = pmConceptsAverageCells(output->parent, inputCells, NULL, NULL, true);
-    psFree(inputCells);
-
-    output->data_exists = true;
-    output->parent->data_exists = true;
-    output->parent->parent->data_exists = true;
-
-    return success;
+
+    // fprintf (stderr, "done with combine %x : %f sec\n", (unsigned int) id, psTimerMark (name));
+    return true;
 }
 
Index: /branches/eam_branch_20080719/psModules/src/imcombine/pmReadoutCombine.h
===================================================================
--- /branches/eam_branch_20080719/psModules/src/imcombine/pmReadoutCombine.h	(revision 18810)
+++ /branches/eam_branch_20080719/psModules/src/imcombine/pmReadoutCombine.h	(revision 18811)
@@ -5,6 +5,6 @@
  * @author Paul Price, IfA
  *
- * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-03-29 03:10:17 $
+ * @version $Revision: 1.13.20.1 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2008-07-31 04:01:40 $
  * Copyright 2004-2006 Institute for Astronomy, University of Hawaii
  */
@@ -36,4 +36,7 @@
                                      );
 
+// check the input parameters and set up the output images
+bool pmReadoutCombinePrepare(pmReadout *output, const psArray *inputs, const pmCombineParams *params);
+
 /// Combine multiple readouts, applying zero and scale, with optional minmax clipping
 bool pmReadoutCombine(pmReadout *output,///< Output readout; altered and returned