Index: trunk/psModules/src/imcombine/pmReadoutCombine.c
===================================================================
--- trunk/psModules/src/imcombine/pmReadoutCombine.c	(revision 9730)
+++ trunk/psModules/src/imcombine/pmReadoutCombine.c	(revision 10248)
@@ -32,4 +32,5 @@
     params->iter = 1;
     params->rej = INFINITY;
+    params->weights = false;
 
     return params;
@@ -77,5 +78,4 @@
     int xSize = 0, ySize = 0;           // The size of the output image
 
-    bool haveWeights = false;           // Do we have weight images?
     bool valid = false;                 // Do we have a single valid input?
     for (long i = 0; i < inputs->n; i++) {
@@ -97,16 +97,10 @@
         }
 
-        if (readout->weight) {
-            if (valid && !haveWeights) {
-                psLogMsg(__func__, PS_LOG_WARN, "Readout %ld has a weight map, but others don't --- "
-                         "weights ignored.\n", i);
-            } else {
-                haveWeights = true;
-            }
-        } else if (haveWeights) {
-            psLogMsg(__func__, PS_LOG_WARN, "Readout %ld doesn't have a weight map, but others do --- "
-                     "weights ignored.\n", i);
-            haveWeights = false;
-        }
+        if (params->weights && !readout->weight) {
+            psError(PS_ERR_UNEXPECTED_NULL, true,
+                    "Rejection based on weights requested, but no weights supplied for image %ld.\n", i);
+            return false;
+        }
+
         valid = true;
 
@@ -160,5 +154,5 @@
 
     psStatsOptions combineStdev = 0; // Statistics option for weights
-    if (haveWeights) {
+    if (params->weights) {
 
         if (!output->weight) {
@@ -207,6 +201,7 @@
 
     psVector *weights = NULL;           // Stack of weights
+    psVector *errors = NULL;            // Stack of errors (sqrt of variance/weights), for psVectorStats
     psF32 *weightsData = NULL;          // Dereference weights
-    if (haveWeights) {
+    if (params->weights) {
         weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights
         weightsData = weights->data.F32;
@@ -272,4 +267,8 @@
                 #endif
                 pixelsData[r] = image->data.F32[yIn][xIn];
+                if (params->weights) {
+                    weightsData[r] = readout->weight->data.F32[yIn][xIn];
+                }
+
                 if (zero) {
                     pixelsData[r] -= zero->data.F32[r];
@@ -277,4 +276,7 @@
                 if (scale) {
                     pixelsData[r] *= invScale->data.F32[r];
+                    if (params->weights) {
+                        weightsData[r] *= invScale->data.F32[r] * invScale->data.F32[r];
+                    }
                 }
 
@@ -319,8 +321,12 @@
             }
 
+            if (params->weights) {
+                errors = (psVector*)psUnaryOp(errors, weights, "sqrt");
+            }
+
             // Combination
-            psVectorStats(stats, pixels, weights, mask, 1);
+            psVectorStats(stats, pixels, errors, mask, 1);
             outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
-            if (haveWeights) {
+            if (params->weights) {
                 float stdev = psStatsGetValue(stats, combineStdev);
                 outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance
@@ -337,4 +343,5 @@
     psFree(mask);
     psFree(weights);
+    psFree(errors);
     psFree(stats);
     psFree(invScale);