Changeset 21477 for trunk/ppStack/src/ppStackMatch.c

Timestamp:

Feb 13, 2009, 1:53:56 PM (17 years ago)

Author:

Paul Price

Message:

Using covariance factor to set the correct variance level for the combination. Adding additional recipe parameter to provide the fraction of values to discard when performing a 'olympic weighted mean'.

File:

• trunk/ppStack/src/ppStackMatch.c (modified) (4 diffs)

trunk/ppStack/src/ppStackMatch.c

@@ -170 +170 @@
     assert(kernels && !*kernels);
     assert(config);
+    *weighting = 0.0;
 
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
@@ -260 +261 @@
         pmSubtractionAnalysis(readout->analysis, kernels, region,
                               readout->image->numCols, readout->image->numRows);
+
+        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
+        psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // New covariance matrix
+        psFree(readout->covariance);
+        readout->covariance = covar;
+        psFree(kernel);
+
     } else {
 #endif
@@ -541 +549 @@
         }
         psFree(iter);
-        *chi2 /= num;
+        *chi2 /= psImageCovarianceFactor(readout->covariance) * num;
     }
 
@@ -579 +587 @@
         (void)psBinaryOp(readout->image, readout->image, "-",
                          psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
-        *weighting = 1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
-        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
-                         "Weighting by 1/noise^2 for stack",
-                         1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)));
-    }
+    }
+
+    // Measure the variance level for the weighting
+    if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image.");
+        psFree(output);
+        return false;
+    }
+    *weighting = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
+                        psImageCovarianceFactor(readout->covariance));
+    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
+                     "Weighting by 1/noise^2 for stack", *weighting);
+
     psFree(bg);