Index: /branches/pap/ppSub/src/ppSubReadoutSubtract.c
===================================================================
--- /branches/pap/ppSub/src/ppSubReadoutSubtract.c	(revision 27708)
+++ /branches/pap/ppSub/src/ppSubReadoutSubtract.c	(revision 27709)
@@ -50,11 +50,10 @@
     outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);
 
-
-    // Combining the variance is tricky.  The problem is, we're representing the variance by a variance map
-    // and a covariance pseudo-matrix, and we cannot simply add them and get the correct result, because they
-    // don't combine nicely.  It's like we want dD = aA + bB, but the upper- and lower-case variables can't be
-    // merged.  Our solution is to combine the variance maps and covariance pseudo-matrices all normalised to
-    // unity, and then scale the result by the desired variance normalisation.  It's not perfectly correct,
-    // but it's gotta be good enough.
+    // Combining the variance is not straight-forward.  The problem is, we're representing the variance by a
+    // variance map and a covariance pseudo-matrix, and we cannot simply add them and get the correct result,
+    // because they don't combine nicely.  It's like we want dD = aA + bB, but the upper- and lower-case
+    // variables can't be merged.  Our solution is to pull the normalisation of the covariance pseudo-matrices
+    // into the variance map, sum the scaled variance maps and (weighted) average the covariance
+    // pseudo-matrices.  It's not perfectly correct, but it's gotta be good enough.
 
     // Measure the variance scales
@@ -89,16 +88,18 @@
     psFree(covarWeights);
 
+    // Restore the scales, in case someone cares to look at the noise in the convolved images
+    psBinaryOp(minuend->covariance->image, minuend->covariance->image, "*",
+               psScalarAlloc(minuendSum, PS_TYPE_F32));
+    psBinaryOp(subtrahend->covariance->image, subtrahend->covariance->image, "*",
+               psScalarAlloc(subtrahendSum, PS_TYPE_F32));
+
     // Combine the variance maps
     int numCols = outRO->image->numCols, numRows = outRO->image->numRows; // Size of image
-    outRO->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-    // The factor of 0.5 in the target variance is because there's an addition here *and* an addition below,
-    // and we only want the first to count (the second should add to 2 in the background).
-    float target = (minuendVar * minuendSum + subtrahendVar * subtrahendSum); // Target mean variance
-    fprintf(stderr, "%f %f %f %f --> %f %f\n", minuendVar, minuendSum, subtrahendVar, subtrahendSum, target, psImageCovarianceFactor(outRO->covariance));
+    psImage *outVar = outRO->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Output variance map
+    psImage *minuendVariance = minuend->variance, *subtrahendVariance = subtrahend->variance; // Variance maps
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
-            outRO->variance->data.F32[y][x] = target *
-                (minuend->variance->data.F32[y][x] / minuendVar +
-                 subtrahend->variance->data.F32[y][x] / subtrahendVar);
+            outVar->data.F32[y][x] = minuendVariance->data.F32[y][x] * minuendSum +
+                subtrahendVariance->data.F32[y][x] * subtrahendSum;
         }
     }