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

Timestamp:

Apr 9, 2013, 1:33:36 PM (13 years ago)

Author:

watersc1

Message:

ppStack changes to fix source matching bug.

File:

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

trunk/ppStack/src/ppStackMatch.c

@@ -8 +8 @@
 // #define TESTING                         // Enable debugging output
 // #define TESTING_REUSE_CONV           // Enable debugging for re-using convolved outputs from previous run
+
+// #define TESTING_CZW 
 
 #ifdef TESTING_REUSE_CONV
@@ -264 +266 @@
             // Add the background into the target image
             psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
-            psBinaryOp(fake->image, fake->image, "+", bgImage);
+            psBinaryOp(fake->image, fake->image, "+", bgImage);
             psFree(bgImage);
 #endif
 
-#ifdef TESTING
+#ifdef TESTING_CZW
             {
                 pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
-                psStringAppend(&name, "fake_%03d.fits", index);
+                psStringAppend(&name, "/local/tmp/fake_%03d.fits", index);
                 pmStackVisualPlotTestImage(fake->image, name);
                 psFits *fits = psFitsOpen(name, "w");
@@ -282 +284 @@
                 pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
-                psStringAppend(&name, "real_%03d.fits", index);
+                psStringAppend(&name, "/local/tmp/real_%03d.fits", index);
                 pmStackVisualPlotTestImage(readout->image, name);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
                 psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
+                psFitsClose(fits);
+            }
+            {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
+                psString name = NULL;
+                psStringAppend(&name, "/local/tmp/real_var_%03d.fits", index);
+                pmStackVisualPlotTestImage(readout->variance, name);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, hdu->header, readout->variance, 0, NULL);
                 psFitsClose(fits);
             }
@@ -355 +367 @@
 
 
-#ifdef TESTING
+#ifdef TESTING_CZW
             {
                 pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
-                psStringAppend(&name, "conv_%03d.fits", index);
+                psStringAppend(&name, "/local/tmp/conv_%03d.fits", index);
                 pmStackVisualPlotTestImage(conv->image, name);
                 psFits *fits = psFitsOpen(name, "w");
@@ -369 +381 @@
                 pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
-                psStringAppend(&name, "diff_%03d.fits", index);
+                psStringAppend(&name, "/local/tmp/conv_var_%03d.fits", index);
+                pmStackVisualPlotTestImage(conv->variance, name);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, hdu->header, conv->variance, 0, NULL);
+                psFitsClose(fits);
+            }
+            {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
+                psString name = NULL;
+                psStringAppend(&name, "/local/tmp/diff_%03d.fits", index);
                 pmStackVisualPlotTestImage(fake->image, name);
                 psFits *fits = psFitsOpen(name, "w");
@@ -510 +532 @@
     }
 
-    if (!stackRenormaliseReadout(config, readout)) {
-        psFree(rng);
-        psFree(bg);
-        return false;
-    }
-
+#ifdef TESTING_CZW
+    // CZW begin check of weight before renormalize step
     // Measure the variance level for the weighting
     if (psMetadataLookupBool(NULL, recipe, "WEIGHTS")) {
@@ -529 +547 @@
         options->weightings->data.F32[index] = 1.0;
     }
-    psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n",
-             index, options->weightings->data.F32[index]);
+    psLogMsg("ppStack", PS_LOG_INFO, "Prenorm Weighting for image %d is %g covar %g median %g\n",
+             index, options->weightings->data.F32[index],psImageCovarianceFactor(readout->covariance),
+             psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN));
+    // CZW end
+#endif
+    if (!stackRenormaliseReadout(config, readout)) {
+        psFree(rng);
+        psFree(bg);
+        return false;
+    }
+    
+    // Measure the variance level for the weighting
+    if (psMetadataLookupBool(NULL, recipe, "WEIGHTS")) {
+        if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
+            psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image.");
+            psFree(rng);
+            psFree(bg);
+            return false;
+        }
+        options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
+                                                      psImageCovarianceFactor(readout->covariance));
+    } else {
+        options->weightings->data.F32[index] = 1.0;
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %g covar %g median %g\n",
+             index, options->weightings->data.F32[index],psImageCovarianceFactor(readout->covariance),
+             psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN));
 
     psFree(rng);