Index: branches/pap/ppSub/src/ppSubReadoutSubtract.c
===================================================================
--- trunk/ppSub/src/ppSubReadoutSubtract.c	(revision 26982)
+++ branches/pap/ppSub/src/ppSubReadoutSubtract.c	(revision 27708)
@@ -48,14 +48,60 @@
     pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
     outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
-    if (minuend->variance && subtrahend->variance) {
-        outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance);
-    }
     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.
+
+    // Measure the variance scales
+    psStats *varStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance images
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);           // Random number generator
+    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config) |
+        pmConfigMaskGet("BLANK", config); // Bits to mask in inputs
+    psImageBackground(varStats, NULL, minuend->variance, minuend->mask, maskVal, rng);
+    float minuendVar = varStats->robustMedian; // Mean variance for minuend
+    psImageBackground(varStats, NULL, subtrahend->variance, subtrahend->mask, maskVal, rng);
+    float subtrahendVar = varStats->robustMedian; // Mean variance for subtrahend
+    psFree(varStats);
+    psFree(rng);
+
+    // Measure and apply the covariance scales
+    float minuendSum = psImageCovarianceSum(minuend->covariance); // Sum of covariance for minuend
+    float subtrahendSum = psImageCovarianceSum(subtrahend->covariance); // Sum of covariance for subtrahend
+    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 scaled covariances
     psArray *covars = psArrayAlloc(2);  // Covariance pseudo-matrices
+    psVector *covarWeights = psVectorAlloc(2, PS_TYPE_F32); // Weights for covariances
     covars->data[0] = psMemIncrRefCounter(minuend->covariance);
     covars->data[1] = psMemIncrRefCounter(subtrahend->covariance);
-    outRO->covariance = psImageCovarianceAverage(covars);
+    covarWeights->data.F32[0] = minuendSum * minuendVar;
+    covarWeights->data.F32[1] = subtrahendSum * subtrahendVar;
+    outRO->covariance = psImageCovarianceAverageWeighted(covars, covarWeights);
     psFree(covars);
+    psFree(covarWeights);
+
+    // 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));
+    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);
+        }
+    }
+    psImageCovarianceTransfer(outRO->variance, outRO->covariance);
 
     outRO->data_exists = true;