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

Timestamp:

Sep 11, 2008, 6:12:42 PM (18 years ago)

Author:

Paul Price

Message:

Correcting stack variances for the imperfect convolution. Changed (again) the way the deviations are calculated; they're now straight chi^{2, which allows us to estimate how the variance factor needs to be changed in order to correct for it. Moving additional values into pmSubtractionKernels so that this can be calculated from a kernel that has been read in (very useful for debugging, also nice to have a record of the quality of the fit).}

File:

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

trunk/ppStack/src/ppStackMatch.c

@@ -15 +15 @@
                      PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
 
-//#define TESTING                         // Enable debugging output
+#define TESTING                         // Enable debugging output
 
 
@@ -48 +48 @@
 #endif
 
-bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels,
+bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2,
                   const psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
 {
@@ -111 +111 @@
                                                                 PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
             float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
-            psMetadataItem *vfItem = psMetadataLookup(output->parent->concepts, "CELL.VARFACTOR");
+            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
             if (!isfinite(vf)) {
                 vf = 1.0;
@@ -341 +341 @@
     assert((*regions)->n == (*kernels)->n);
 
+    // Correct the variance for the chi^2
+    {
+        *chi2 = 0.0;
+        int num = 0;                    // Number of measurements of chi^2
+        psString regex = NULL;          // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
+        psFree(regex);
+        psMetadataItem *item = NULL;// Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            assert(item->type == PS_DATA_UNKNOWN);
+            pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
+            *chi2 += kernels->mean;
+            num++;
+        }
+        psFree(iter);
+
+        float vf = psMetadataLookupF32(NULL, readout->parent->concepts, "CELL.VARFACTOR"); // Variance factor
+        *chi2 /= vf * num;
+    }
+
     // Renormalise the variances if desired
     if (renorm) {