Changeset 35383

Timestamp:

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

Author:

watersc1

Message:

ppStack changes to fix source matching bug.

Location:

trunk

Files:

• ppStack/src/ppStackMatch.c (modified) (7 diffs)
• ppStack/src/ppStackSources.c (modified) (11 diffs)
• ppStack/src/ppStackTarget.c (modified) (2 diffs)
• psModules/src/objects/pmSourceMatch.c (modified) (3 diffs)
• psModules/src/objects/pmSourceMatch.h (modified) (1 diff)

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);

trunk/ppStack/src/ppStackSources.c

@@ -11 +11 @@
 #endif
 
-#ifdef TESTING
+// #define TESTING_CZW
+
+#ifdef TESTING_CZW
 // Dump matches to a file
 static void dumpMatches(const char *filename, // File to which to dump
@@ -28 +30 @@
         psVectorInit(magErr, NAN);
         for (int j = 0; j < match->num; j++) {
+          if (match->mask) {
             if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
-                continue;
-            }
+              continue;
+            }
+          }
             int index = match->image->data.U32[j]; // Image index
             mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index];
@@ -39 +43 @@
         }
         for (int j = 0; j < num; j++) {
-            fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]);
+          fprintf(outMatches, "%f %f %f %f %f %f\t",
+                  match->mag->data.F32[j],
+                  zp->data.F32[j],
+                  trans ? trans->data.F32[j] : NAN,
+                  magErr->data.F32[j],
+                  match->x->data.F32[j],
+                  match->y->data.F32[j]);
         }
         fprintf(outMatches, "\n");
@@ -271 +281 @@
         return false;
     }
-
-#ifdef TESTING
+    options->sources = psMemIncrRefCounter(matches);
+    psArray *matchedSources = psArrayAlloc(0);
+
+#ifdef TESTING_CZW
     dumpMatches("source_match.dat", num, matches, zp, NULL);
 #endif
 
     if (options->matchZPs) {
-        psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
+
+      psVector *trans = pmSourceMatchRelphot(matches, matchedSources, zp, tol, iter1, starRej1, starSys1,
                                                iter2, starRej2, starSys2, starLimit,
                                                transIter, transRej, transThresh); // Transparencies per image
@@ -324 +337 @@
 
         // m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i
-
+#ifdef TESTING_CZW
+    dumpMatches("source_match2.dat", num, matches, zp, trans);
+#endif
         if (options->matchZPs) {
             options->norm = psVectorAlloc(num, PS_TYPE_F32);
             for (int i = 0; i < num; i++) {
-                if (inputMask->data.U8[i] || !isfinite(trans->data.F32[i])) {
-                    continue;
-                }
+              if (!isfinite(trans->data.F32[i])) {
+                continue;
+              }
                 psArray *sources = sourceLists->data[i]; // Sources of interest
                 float magCorr = zp->data.F32[i] - trans->data.F32[i] - 2.5*log10(options->sumExposure) - airmassTerm * airmassTarget;
@@ -341 +356 @@
                 options->norm->data.F32[i] = magCorr;
                 psLogMsg("ppStack", PS_LOG_INFO,
-                         "Applying scale correction to image %d: %f mag (%f)\n",
-                         i, magCorr, trans->data.F32[i]);
+                         "Applying scale correction to image %d: %f mag (%f) (with inputMask value %d)\n",
+                         i, magCorr, trans->data.F32[i],inputMask->data.U8[i]);
 
                 for (int j = 0; j < sources->n; j++) {
                     pmSource *source = sources->data[j]; // Source of interest
+                    source->psfMag += magCorr;
+                    source->apMag += magCorr;
+
                     if (!source) {
                         continue;
                     }
                     // XXX need to apply to apMag as well
-                    source->psfMag += magCorr;
-                    source->apMag += magCorr;
+                    psTrace("ppStack",5,"Source corrections: %d %d %f %f %f %f %f %f %f %f\n",i,j,
+                            source->peak->xf,source->peak->yf,
+                            source->psfMag,source->psfMagErr,
+                            airmassTerm,airmassTarget,2.5 * log10(options->sumExposure),
+                            source->psfMag + magCorr);
                 }
             }
         }
+        matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
 
         if (zpExpNum == numGoodImages) {
@@ -366 +388 @@
 
 
-#ifdef TESTING
-        dumpMatches("source_mags.dat", num, matches, zp, trans);
-#endif
         psFree(trans);
 
-#ifdef TESTING
+#ifdef TESTING_CZW
         // Double check: all transparencies should be zero
         {
             psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
+            psArray *matchedSources = psArrayAlloc(matches->n);
             if (!matches) {
                 psError(PPSTACK_ERR_DATA, false, "Unable to match sources");
@@ -380 +400 @@
                 return false;
             }
-            psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
+            psVector *trans = pmSourceMatchRelphot(matches, matchedSources, zp, tol, iter1, starRej1, starSys1,
                                                    iter2, starRej2, starSys2, starLimit,
                                                    transIter, transRej, transThresh); // Transparencies
@@ -386 +406 @@
                 fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
             }
+
             psFree(trans);
             psFree(matches);
-        }
-#endif
-    }
-
-    psFree(zp);
+            psFree(matchedSources);
+        }
+#endif
+    }
+
+
     psFree(zpExp);
 
 #ifdef ASTROMETRY
+    // CZW: This is off by default. 
     // Position offsets
     {
@@ -436 +459 @@
 
 #if (defined TESTING && defined ASTROMETRY)
+    // CZW: This is off by default.
         // Double check: all offsets should be zero
         {
@@ -458 +482 @@
     if (options->photometry) {
         // Save best matches for future photometry
-        psArray *sourcesBest = options->sources = psArrayAllocEmpty(matches->n);
-
+      psArray *sourcesBest = options->sources = psArrayAllocEmpty(matches->n);
         // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible
         //int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required
         //MEH - if many poor images, will hit oddity with num fraction. use numGoodImages instead
-        int minMatches = PS_MAX(2, fracMatch * numGoodImages);// Minimum number of matches required
-
-        int good0 = 0;
-        for (int i = 0; i < num; i++) {
-            if (inputMask->data.U8[i]) {
-                continue;
-            }
-            good0 = i;
-            break;
+      int minMatches = PS_MAX(2, fracMatch * numGoodImages);// Minimum number of matches required
+      
+      int good0 = 0;
+      for (int i = 0; i < num; i++) {
+        if (inputMask->data.U8[i]) {
+          continue;
         }
-
-        for (int i = 0; i < matches->n; i++) {
-            pmSourceMatch *match = matches->data[i]; // Match of interest
-            if (match->num < minMatches) {
-                continue;
-            }
-
-            // We need to grab a single instance of this source: just take the first available
-            // MEH - if first available is a rejected image, then problem. so take first non-rejected instance above
-            //int image = match->image->data.S32[0]; // Index of image
-            //int index = match->index->data.S32[0]; // Index of source within image
-            int image = match->image->data.S32[good0]; // Index of image
-            int index = match->index->data.S32[good0]; // Index of source within image
-            psArray *sources = sourceLists->data[image]; // Sources for image
-            pmSource *source = sources->data[index]; // Source of interest
-
-            psArrayAdd(sourcesBest, sourcesBest->n, source);
-        }
-        psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
-    }
-
+        good0 = i;
+        break;
+      }
+      
+      for (int i = 0; i < matches->n; i++) {
+        pmSourceMatch *match = matches->data[i]; // Match of interest
+         if (match->num < minMatches) {
+           continue;
+         }
+
+         // We need to grab a single instance of this source: just take the first available
+         // MEH - if first available is a rejected image, then problem. so take first non-rejected instance above
+         //      int image = match->image->data.S32[0]; // Index of image
+         // int index = match->index->data.S32[0]; // Index of source within image
+         int image = match->image->data.S32[good0]; // Index of image
+         int index = match->index->data.S32[good0]; // Index of source within image
+         psArray *sources = sourceLists->data[image]; // Sources for image
+         pmSource *source = sources->data[index]; // Source of interest
+
+         // CZW: Accept the first source for all properties, except for the magnitudes.
+         //      We have average magnitudes from relphot, so we should use those
+         pmSource *photsource = matchedSources->data[i];
+         if ((photsource)&&(isfinite(photsource->psfMag))) {
+           // CZW: photsource is the average relphot mag, which is equal to <m_i + zpIN_j>
+           //      (primed are averaged values for this discussion)
+           //      Therefore, the correct magnitude for psfMag' is:
+           //         psfMag' = (photsource - zpIN_j + magCorr)
+           //      apMag - psfMag should be retained, so:
+           //         apMag - psfMag = apMag' - psfMag'
+           //         apMag' = apMag + psfMag' - psfMag
+           source->apMag  = source->apMag + (photsource->psfMag + options->norm->data.F32[image] - zp->data.F32[image] - source->psfMag);
+           source->psfMag = photsource->psfMag + options->norm->data.F32[image] - zp->data.F32[image];
+           
+           psTrace("ppStack", 5, "Best Source %d/%ld @ %f %f %f %f %f %d",i,sourcesBest->n,source->peak->xf,source->peak->yf,
+                   source->psfMag,
+                   photsource->psfMag + options->norm->data.F32[image] - zp->data.F32[image] ,
+                   source->psfMagErr,-1);
+           psArrayAdd(sourcesBest, sourcesBest->n, source);
+         }
+        }
+      psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
+
+    }
+    psFree(zp);
     psFree(matches);
 

trunk/ppStack/src/ppStackTarget.c

@@ -60 +60 @@
         }
         float xSource, ySource;         // Coordinates of source
-        coordsFromSource(&xSource, &ySource, NULL, source);
+        float Smag;                      // magnitude of source
+        coordsFromSource(&xSource, &ySource, &Smag, source);
 
         coords->data.F64[0] = xSource;
@@ -66 +67 @@
 
         psVector *indices = psTreeAllWithin(tree, coords, exclusion); // Number within exclusion zone
-        psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
-                coords->data.F64[0], coords->data.F64[1], indices->n);
+        psTrace("ppStack", 5, "Filtering: Source at %.0lf %.0lf with mag %f has %ld sources in exclusion zone",
+                coords->data.F64[0], coords->data.F64[1], Smag, indices->n);
         if (indices->n == 1) {
             // Only itself inside the exclusion zone

trunk/psModules/src/objects/pmSourceMatch.c

@@ -630 +630 @@
 
 psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches
+                               psArray *matchedSources, // Array of averaged sources
                                const psVector *zp, // Zero points for each image (including airmass term)
                                float tol, // Relative tolerance for convergence
@@ -648 +649 @@
     PS_ASSERT_VECTOR_TYPE(zp, PS_TYPE_F32, NULL);
     PS_ASSERT_FLOAT_LARGER_THAN(transClip, 0.0, NULL);
-
+    PS_ASSERT_ARRAY_NON_NULL(matchedSources, NULL);
+    
     sys1 *= sys1;
     sys2 *= sys2;
@@ -742 +744 @@
     }
 
+    for (int i = 0; i < matches->n; i++) {
+      pmSource *source = pmSourceAlloc();
+      source->psfMag = stars->data.F32[i];
+      psArrayAdd(matchedSources,matchedSources->n,source);
+    }
+    
     psFree(photo);
     psFree(badImage);

trunk/psModules/src/objects/pmSourceMatch.h

@@ -74 +74 @@
 /// Perform relative photometry to calibrate images
 psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches
+                               psArray *matchedSources, // Array of average sources
                                const psVector *zp, // Zero points for each image (including airmass term)
                                float tol, // Relative tolerance for convergence