Changeset 27708

Timestamp:

Apr 16, 2010, 12:52:37 PM (16 years ago)

Author:

Paul Price

Message:

Branching to test fix of variance/covariance calculation in ppSub.

The calculations of the variance and covariance are usually straight-forward --- the variance gets scaled by the square of any scaling factor applied to the image, and the covariance may be calculated from a convolution kernel. However, when two images are combined, it's not so straight-forward. The standard rules for the propagation of uncertainty tell us to sum the (full) covariance matrices, but it's not immediately clear how to do this in our model where we've separated the diagonal and off-diagonal elements. I had been summing the variance maps and averaging the covariance pseudo-matrices, but this appears to be robbing us of nearly 1 mag of depth. I've reworked this calculation to be a bit smarter: I normalise the (mean of the) variance maps and the covariance pseudo-matrices to unity, calculate from the normalisations the desired mean variance, sum the normalised variance maps with the appropriate scaling, and average (with weights) the covariance pseudo-matrices. This results in deeper photometry on the diff, but there are still a few things to consider before we put this into production: the effect on Magic and the effect on warp-stack and stack-stack diffs (should probably review stack variance/covariance in the light of this fix).

Location:

branches/pap

Files:

• . (copied) (copied from trunk )
• ppSub/src/ppSubReadoutSubtract.c (modified) (1 diff)

branches/pap/ppSub/src/ppSubReadoutSubtract.c

@@ -48 +48 @@
     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;