Changeset 35383 for trunk/ppStack/src/ppStackSources.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/ppStackSources.c (modified) (11 diffs)

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