Changeset 23239 for branches/cnb_branches/cnb_branch_20090215/ppStack/src/ppStackSources.c

Timestamp:

Mar 9, 2009, 3:41:26 PM (17 years ago)

Author:

beaumont

Message:

sync with trunk

Location:

branches/cnb_branches/cnb_branch_20090215

Files:

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

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

@@ -81 +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
@@ -110 +117 @@
         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);
@@ -153 +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");
@@ -164 +172 @@
 #endif
 
-    psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
-                                           transThresh, starRej, starSys); // Transparencies for each image
+    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");
@@ -184 +193 @@
     // 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;
+        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);
+        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);