Index: trunk/psModules/src/imcombine/pmReadoutCombine.c
===================================================================
--- trunk/psModules/src/imcombine/pmReadoutCombine.c	(revision 17137)
+++ trunk/psModules/src/imcombine/pmReadoutCombine.c	(revision 17228)
@@ -18,4 +18,5 @@
 
 //#define SHOW_BUSY 1                   // Show that the function is busy
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -94,5 +95,24 @@
     }
 
-    psStats *stats = psStatsAlloc(params->combine); // The statistics to use in the combination
+    psStatsOptions combineStdev = 0; // Statistics option for weights
+    switch (params->combine) {
+      case PS_STAT_SAMPLE_MEAN:
+      case PS_STAT_SAMPLE_MEDIAN:
+        combineStdev = PS_STAT_SAMPLE_STDEV;
+        break;
+      case PS_STAT_ROBUST_MEDIAN:
+        combineStdev = PS_STAT_ROBUST_STDEV;
+        break;
+      case PS_STAT_FITTED_MEAN:
+        combineStdev = PS_STAT_FITTED_STDEV;
+        break;
+      case PS_STAT_CLIPPED_MEAN:
+        combineStdev = PS_STAT_CLIPPED_STDEV;
+        break;
+      default:
+        psAbort("Should never get here --- checked params->combine before.\n");
+    }
+
+    psStats *stats = psStatsAlloc(params->combine | combineStdev); // The statistics to use in the combination
     if (params->combine == PS_STAT_CLIPPED_MEAN) {
         stats->clipSigma = params->rej;
@@ -107,4 +127,8 @@
         }
     }
+    if (params->weights && first) {
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
+                         "Using input weights to combine images", "");
+    }
 
     int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
@@ -120,31 +144,17 @@
     psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
 
-    psStatsOptions combineStdev = 0; // Statistics option for weights
-    if (params->weights) {
-        if (first) {
-            psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
-                             "Using input weights to combine images", "");
-        }
-
-        // Get the correct statistics option for weights
-        switch (params->combine) {
-        case PS_STAT_SAMPLE_MEAN:
-        case PS_STAT_SAMPLE_MEDIAN:
-            combineStdev = PS_STAT_SAMPLE_STDEV;
-            break;
-        case PS_STAT_ROBUST_MEDIAN:
-            combineStdev = PS_STAT_ROBUST_STDEV;
-            break;
-        case PS_STAT_FITTED_MEAN:
-            combineStdev = PS_STAT_FITTED_STDEV;
-            break;
-        case PS_STAT_CLIPPED_MEAN:
-            combineStdev = PS_STAT_CLIPPED_STDEV;
-            break;
-        default:
-            psAbort("Should never get here --- checked params->combine before.\n");
-        }
-        stats->options |= combineStdev;
-    }
+    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
@@ -273,4 +283,6 @@
                 outputMask[yOut][xOut] = params->blank;
                 outputImage[yOut][xOut] = NAN;
+                counts->data.U16[yOut][xOut] = 0;
+                sigma->data.F32[yOut][xOut] = NAN;
                 continue;
             }
@@ -288,4 +300,5 @@
                         maskData[indexData[k]] = 1;
                         numMasked++;
+                        numValid--;
                     }
                 }
@@ -293,10 +306,12 @@
                 for (int k = pixels->n - 1, numMasked = 0; numMasked < numHigh && k >= 0; k--) {
                     // Don't count the ones that are already masked
-                    if (! maskData[indexData[k]]) {
+                    if (!maskData[indexData[k]]) {
                         maskData[indexData[k]] = 1;
                         numMasked++;
+                        numValid--;
                     }
                 }
             }
+            counts->data.U16[yOut][xOut] = numValid;
 
             // XXXXX this step probably is very expensive : convert errors to variance everywhere?
@@ -314,4 +329,5 @@
                     outputWeight[yOut][xOut] = NAN;
                 }
+                sigma->data.F32[yOut][xOut] = NAN;
             } else {
                 outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
@@ -323,4 +339,5 @@
                     // also, the weighted mean is not obviously the correct thing here
                 }
+                sigma->data.F32[yOut][xOut] = psStatsGetValue(stats, combineStdev);
             }
         }
@@ -338,4 +355,6 @@
     psFree(stats);
     psFree(invScale);
+    psFree(counts);
+    psFree(sigma);
 
     // Update the "concepts"
Index: trunk/psModules/src/imcombine/pmReadoutCombine.h
===================================================================
--- trunk/psModules/src/imcombine/pmReadoutCombine.h	(revision 17137)
+++ trunk/psModules/src/imcombine/pmReadoutCombine.h	(revision 17228)
@@ -5,6 +5,6 @@
  * @author Paul Price, IfA
  *
- * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-06-02 03:51:03 $
+ * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2008-03-29 03:10:17 $
  * Copyright 2004-2006 Institute for Astronomy, University of Hawaii
  */
@@ -20,6 +20,5 @@
 /// These values define how the combination is performed, and should not vary by detector, so that it can be
 /// re-used for multiple combinations.
-typedef struct
-{
+typedef struct {
     psStatsOptions combine;             ///< Statistic to use when performing the combination
     psMaskType maskVal;                 ///< Mask value
@@ -31,6 +30,5 @@
     float rej;                          ///< Rejection threshould for clipping (for CLIPPED_MEAN only)
     bool weights;                       ///< Use the supplied weights (instead of calculated stdev)?
-}
-pmCombineParams;
+} pmCombineParams;
 
 // Allocator for pmCombineParams