Changeset 23352 for branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackSources.c

Timestamp:

Mar 17, 2009, 12:08:50 PM (17 years ago)

Author:

beaumont

Message:

merged with head

Location:

branches/cnb_branches/cnb_branch_20090301

Files:

• . (modified) (1 prop)
• ppStack (modified) (1 prop)
• ppStack/src/ppStackSources.c (modified) (6 diffs)

branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackSources.c

@@ -13 +13 @@
 #define FAKE_ROWS 4913
 
-float ppStackSourcesTransparency(const psArray *sourceLists, const pmFPAview *view, const pmConfig *config)
+#ifdef TESTING
+// Dump matches to a file
+static void dumpMatches(const char *filename, // File to which to dump
+                        int num,        // Number of inputs
+                        psArray *matches, // Star matches
+                        psVector *zp,   // Zero points
+                        psVector *trans // Transparencies
+                        )
+{
+    FILE *outMatches = fopen(filename, "w"); // Output matches
+    psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star
+    psVector *magErr = psVectorAlloc(num, PS_TYPE_F32); // Errors for each star
+    for (int i = 0; i < matches->n; i++) {
+        pmSourceMatch *match = matches->data[i]; // Match of interest
+        psVectorInit(mag, NAN);
+        psVectorInit(magErr, NAN);
+        for (int j = 0; j < match->num; j++) {
+            if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+                continue;
+            }
+            int index = match->image->data.U32[j]; // Image index
+            mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index];
+            if (trans) {
+                mag->data.F32[index] -= trans->data.F32[index];
+            }
+            magErr->data.F32[index] = match->magErr->data.F32[j];
+        }
+        for (int j = 0; j < num; j++) {
+            fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]);
+        }
+        fprintf(outMatches, "\n");
+    }
+    psFree(mag);
+    psFree(magErr);
+    fclose(outMatches);
+    return;
+}
+#endif
+
+
+float ppStackSourcesTransparency(const psArray *sourceLists, psVector *inputMask,
+                                 const pmFPAview *view, const pmConfig *config)
 {
     PS_ASSERT_ARRAY_NON_NULL(sourceLists, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(inputMask, NAN);
+    PS_ASSERT_VECTOR_TYPE(inputMask, PS_TYPE_U8, NAN);
+    PS_ASSERT_VECTOR_SIZE(inputMask, sourceLists->n, NAN);
     PS_ASSERT_PTR_NON_NULL(view, NAN);
     PS_ASSERT_PTR_NON_NULL(config, NAN);
 
-#ifdef TESTING
+#if defined(TESTING) && 0
     {
         // Deliberately induce a major transparency difference
@@ -37 +81 @@
 
     float radius = psMetadataLookupF32(NULL, recipe, "ZP.RADIUS"); // Radius (pixels) for matching sources
-    int iter = psMetadataLookupS32(NULL, recipe, "ZP.ITER"); // Maximum iterations
+    int iter1 = psMetadataLookupS32(NULL, recipe, "ZP.ITER.1"); // Maximum iterations for pass 1
+    int iter2 = psMetadataLookupS32(NULL, recipe, "ZP.ITER.2"); // Maximum iterations for pass 2
     float tol = psMetadataLookupF32(NULL, recipe, "ZP.TOL"); // Tolerance for zero point iterations
     int transIter = psMetadataLookupS32(NULL, recipe, "ZP.TRANS.ITER"); // Iterations for transparency
     float transRej = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.REJ");// Rejection threshold for transparency
     float transThresh = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.THRESH"); // Threshold for transparency
-    float starRej = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ"); // Rejection threshold for stars
+
+    float starRej1 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ.1"); // Rejection threshold for stars
+    float starSys1 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS.1"); // Estimated systematic error
+    float starRej2 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ.2"); // Rejection threshold for stars
+    float starSys2 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS.2"); // Estimated systematic error
+
     float starLimit = psMetadataLookupF32(NULL, recipe, "ZP.STAR.LIMIT"); // Limit on star rejection fraction
-    float starSys = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS"); // Estimated systematic error
+
+    float fracMatch = psMetadataLookupF32(NULL, recipe, "ZP.MATCH"); // Fraction of images to match for star
 
     psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms
@@ -63 +114 @@
         pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
 
-#ifdef TESTING
+#if defined(TESTING) && 0
         pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
         pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
-        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sourceLists->data[i], NULL, NULL, psf, 5, 0, false, true);
+        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sourceLists->data[i],
+                                 NULL, NULL, psf, 5, 0, false, true);
         psString name = NULL;
         psStringAppend(&name, "start_%03d.fits", i);
+        pmStackVisualPlotTestImage(fake->image, name);
         psFits *fits = psFitsOpen(name, "w");
         psFree(name);
@@ -108 +161 @@
     }
 
-    psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
+    psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
     if (!matches) {
         psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
@@ -114 +167 @@
         return NAN;
     }
-    psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
-                                           transThresh, starRej, starSys); // Transparencies for each image
-
-#ifdef TESTING
-    {
-        // Dump the corrected magnitudes
-        FILE *outMatches = fopen("source_match.dat", "w"); // Output matches
-        psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star
-        psVector *magErr = psVectorAlloc(num, PS_TYPE_F32); // Errors for each star
-        for (int i = 0; i < matches->n; i++) {
-            pmSourceMatch *match = matches->data[i]; // Match of interest
-            psVectorInit(mag, NAN);
-            psVectorInit(magErr, NAN);
-            for (int j = 0; j < match->num; j++) {
-                if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
-                    continue;
-                }
-                int index = match->image->data.U32[j]; // Image index
-                mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index] - trans->data.F32[index];
-                magErr->data.F32[index] = match->magErr->data.F32[j];
-            }
-            for (int j = 0; j < num; j++) {
-                fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]);
-            }
-            fprintf(outMatches, "\n");
-        }
-        psFree(mag);
-        psFree(magErr);
-        fclose(outMatches);
-    }
-#endif
+
+#ifdef TESTING
+    dumpMatches("source_match.dat", num, matches, zp, NULL);
+#endif
+
+    psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
+                                           iter2, starRej2, starSys2, starLimit,
+                                           transIter, transRej, transThresh); // Transparencies for each image
+    if (!trans) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
+        return NAN;
+    }
+
+#ifdef TESTING
+    dumpMatches("source_mags.dat", num, matches, zp, trans);
+#endif
+
+    for (int i = 0; i < trans->n; i++) {
+        if (!isfinite(trans->data.F32[i])) {
+            inputMask->data.U8[i] = PPSTACK_MASK_CAL;
+        }
+    }
 
     // Save best matches SOMEWHERE for future photometry
     // XXX this is a really poor output location; clean up the pmFPAfiles used in ppStack
     pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT");
-    psArray *sourcesBest = psArrayAllocEmpty (100);
-
-    // XXX something of a hack: require at 2 detections or 1/2 of the max possible
-    int minMatches = PS_MAX (2, 0.5*num);
+    psArray *sourcesBest = 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
     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
-        int nImage = match->image->data.S32[0];
-        int nIndex = match->index->data.S32[0];
-        psArray *sources = sourceLists->data[nImage];
-        pmSource *source = sources->data[nIndex];
-        
-        // stick this sample source on sourcesBest
-        psArrayAdd (sourcesBest, 100, source);
-    }
-    psMetadataAdd (sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, "psphot sources", sourcesBest);
-    psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis\n", sourcesBest->n);
-    psFree (sourcesBest);
+        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
+        int image = match->image->data.S32[0]; // Index of image
+        int index = match->index->data.S32[0]; // 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);
+    }
+    psMetadataAdd(sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE,
+                  "psphot sources", sourcesBest);
+    psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
+    psFree(sourcesBest);
 
     psFree(matches);
-    if (!trans) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
-        return NAN;
-    }
 
     // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
@@ -182 +224 @@
     // We don't need to know the magnitude zero point for the filter, since it cancels out
     for (int i = 0; i < num; i++) {
+        if (!isfinite(trans->data.F32[i])) {
+            continue;
+        }
         psArray *sources = sourceLists->data[i]; // Sources of interest
         float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];