Changeset 19624 for trunk/ppSub/src/ppSubReadout.c

Timestamp:

Sep 19, 2008, 5:07:47 PM (18 years ago)

Author:

Paul Price

Message:

Adding code to update the variance factor for the output. Note that
this can't be exact, so we cheat a bit: the output weight map is
sigma^{2 = sigma1}2 + sigma2^{2, so we calculate vf = (vf1}2 * sigma1^{2 +
vf2}2 * sigma2^{2) / (sigma1}2 + sigma2^{2).
Not perfectly legit, but should be good enough for government work.}

File:

• trunk/ppSub/src/ppSubReadout.c (modified) (1 diff)

trunk/ppSub/src/ppSubReadout.c

@@ -208 +208 @@
                      "Subtraction kernel", kernels->description);
 
+    {
+        // Update variance factors
+        // It's not possible to do this perfectly, since we're combining different images:
+        // vf_sum sigma_sum^2 = vf_1 * sigma_1^2 + vf_2 * sigma_2^2
+        // It's not possible to write vf_sum as a function of vf_1 and vf_2 with no reference to the sigmas.
+        // Instead, we're going to cheat.
+        bool mdok;                      // Status of MD lookup
+        pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, inConv->analysis,
+                                                            PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
+        if (!mdok) {
+            kernels = psMetadataLookupPtr(&mdok, refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
+        }
+        if (!mdok) {
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find subtraction kernels.");
+            psFree(inConv);
+            psFree(refConv);
+            psFree(outRO);
+            return false;
+        }
+        float vfIn = psMetadataLookupF32(NULL, inRO->parent->concepts,
+                                         "CELL.VARFACTOR"); // Variance factor for input
+        if (!isfinite(vfIn)) {
+            vfIn = 1.0;
+        }
+        float vfRef = psMetadataLookupF32(NULL, refRO->parent->concepts,
+                                          "CELL.VARFACTOR"); // Variance factor for ref
+        if (!isfinite(vfRef)) {
+            vfRef = 1.0;
+        }
+        if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+            vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
+        }
+        if (kernels->mode == PM_SUBTRACTION_MODE_2) {
+            vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
+        } else if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+            vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, true);
+        }
+
+        psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
+        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+        if (!psImageBackground(stats, NULL, inConv->weight, inConv->mask, maskVal | maskBad, NULL)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved input");
+            psFree(inConv);
+            psFree(refConv);
+            psFree(outRO);
+            return false;
+        }
+        float inMeanVar = stats->robustMedian; // Mean variance of input
+        if (!psImageBackground(stats, NULL, refConv->weight, refConv->mask, maskVal | maskBad, NULL)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved reference");
+            psFree(inConv);
+            psFree(refConv);
+            psFree(outRO);
+            return false;
+        }
+        float refMeanVar = stats->robustMedian; // Mean variance of reference
+        psFree(rng);
+        psFree(stats);
+
+        float vf = (vfIn * inMeanVar + vfRef * refMeanVar) / (inMeanVar + refMeanVar); // Resulting var factor
+        psMetadataItem *vfItem = psMetadataLookup(outRO->parent->concepts, "CELL.VARFACTOR");
+        vfItem->data.F32 = vf;
+    }
 
     // Statistics on the matching