Changeset 26674

Timestamp:

Jan 22, 2010, 5:32:32 PM (16 years ago)

Author:

Paul Price

Message:

Get running with scaling of kernel parameters.

Location:

branches/eam_branches/20091201/ppStack/src

Files:

• ppStackMatch.c (modified) (1 diff)
• ppStackOptions.c (modified) (3 diffs)
• ppStackOptions.h (modified) (2 diffs)
• ppStackPhotometry.c (modified) (2 diffs)
• ppStackReadout.c (modified) (2 diffs)
• ppStackSetup.c (modified) (1 diff)
• ppStackSources.c (modified) (3 diffs)

branches/eam_branches/20091201/ppStack/src/ppStackMatch.c

@@ -317 +317 @@
             float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor"
 
-            bool scale = psMetadataLookupBool(NULL, recipe, "SCALE");        // Scale kernel parameters?
-            float scaleRef = psMetadataLookupF32(NULL, recipe, "SCALE.REF"); // Reference for scaling
-            float scaleMin = psMetadataLookupF32(NULL, recipe, "SCALE.MIN"); // Minimum for scaling
-            float scaleMax = psMetadataLookupF32(NULL, recipe, "SCALE.MAX"); // Maximum for scaling
+            bool scale = psMetadataLookupBool(NULL, ppsub, "SCALE");        // Scale kernel parameters?
+            float scaleRef = psMetadataLookupF32(NULL, ppsub, "SCALE.REF"); // Reference for scaling
+            float scaleMin = psMetadataLookupF32(NULL, ppsub, "SCALE.MIN"); // Minimum for scaling
+            float scaleMax = psMetadataLookupF32(NULL, ppsub, "SCALE.MAX"); // Maximum for scaling
             if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) {
                 psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                        "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in recipe.",
+                        "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.",
                         scaleRef, scaleMin, scaleMax);
                 return false;

branches/eam_branches/20091201/ppStack/src/ppStackOptions.c

@@ -24 +24 @@
     psFree(options->sourceLists);
     psFree(options->norm);
+    psFree(options->sources);
     psFree(options->cells);
     psFree(options->kernels);
@@ -44 +45 @@
     options->convolve = true;
     options->matchZPs = true;
+    options->photometry = false;
     options->stats = NULL;
     options->statsFile = NULL;
@@ -60 +62 @@
     options->sourceLists = NULL;
     options->norm = NULL;
+    options->sources = NULL;
     options->cells = NULL;
     options->kernels = NULL;

branches/eam_branches/20091201/ppStack/src/ppStackOptions.h

@@ -10 +10 @@
     bool convolve;                      // Convolve images?
     bool matchZPs;                      // Adjust relative fluxes based on transparency analysis?
+    bool photometry;                    // Perform photometry?
     psMetadata *stats;                  // Statistics for output
     FILE *statsFile;                    // File to which to write statistics
@@ -24 +25 @@
     psArray *sourceLists;               // Individual lists of sources for matching
     psVector *norm;                     // Normalisation for each image
+    psArray *sources;                   // Matched sources
     // Convolve
     psArray *cells;                     // Cells for convolved images --- a handle for reading again

branches/eam_branches/20091201/ppStack/src/ppStackPhotometry.c

@@ -19 +19 @@
     psAssert(recipe, "We've thrown an error on this before.");
 
-    bool mdok;                          // Status of MD lookup
-    if (!psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
+    if (!options->photometry) {
         // Nothing to do
         return true;
@@ -60 +59 @@
                            "Bits to use for mark", markValue);
 
-    pmCell *sourcesCell = pmFPAfileThisCell(config->files, photView, "PPSTACK.OUTPUT");
-    psArray *inSources = psMetadataLookupPtr(NULL, sourcesCell->analysis, "PSPHOT.SOURCES");
+    psArray *inSources = options->sources;
     if (!inSources) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find input sources");

branches/eam_branches/20091201/ppStack/src/ppStackReadout.c

@@ -126 +126 @@
     int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
 
-    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
+    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
@@ -219 +219 @@
     bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
 
-    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
+    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits

branches/eam_branches/20091201/ppStack/src/ppStackSetup.c

@@ -22 +22 @@
     psAssert(recipe, "We've thrown an error on this before.");
 
-    options->matchZPs = psMetadataLookupBool(NULL, recipe, "MATCH.ZERO.POINTS"); // Adjust zero points based on tranparency analysis?
+    options->matchZPs = psMetadataLookupBool(NULL, recipe, "ZP"); // Adjust zero points?
+
+    options->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); // Perform photometry?
 
     options->convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images?

branches/eam_branches/20091201/ppStack/src/ppStackSources.c

@@ -61 +61 @@
     PS_ASSERT_PTR_NON_NULL(config, false);
 
-    if (!options->matchZPs) {
-        int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
-        options->norm = psVectorAlloc(num, PS_TYPE_F32);
-        psVectorInit (options->norm, 0.0);
-
-        // XXX do I need to set this?
-        // options->sumExposure = sumExpTime;
-
-        return true;
+    if (!options->matchZPs && !options->photometry) {
+        int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+        options->norm = psVectorAlloc(num, PS_TYPE_F32);
+        psVectorInit (options->norm, 0.0);
+
+        // XXX do I need to set this?
+        // options->sumExposure = sumExpTime;
+
+        return true;
     }
 
@@ -191 +191 @@
 #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 false;
-    }
+    if (options->matchZPs) {
+        psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
+                                               iter2, starRej2, starSys2, starLimit,
+                                               transIter, transRej, transThresh); // Transparencies per image
+        if (!trans) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
+            return false;
+        }
+        for (int i = 0; i < trans->n; i++) {
+            if (!isfinite(trans->data.F32[i])) {
+                inputMask->data.U8[i] |= PPSTACK_MASK_CAL;
+            }
+        }
+        // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
+        // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
+        // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
+        // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
+        // We don't need to know the magnitude zero point for the filter, since it cancels out
+        if (options->matchZPs) {
+            options->norm = psVectorAlloc(num, PS_TYPE_F32);
+            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];
+                options->norm->data.F32[i] = magCorr;
+                psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n",
+                         i, magCorr);
+
+                for (int j = 0; j < sources->n; j++) {
+                    pmSource *source = sources->data[j]; // Source of interest
+                    if (!source) {
+                        continue;
+                    }
+                    source->psfMag += magCorr;
+                }
+            }
+        }
+
+#ifdef TESTING
+        dumpMatches("source_mags.dat", num, matches, zp, trans);
+#endif
+        psFree(trans);
+
+#ifdef TESTING
+        // Double check: all transparencies should be zero
+        {
+            psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
+            if (!matches) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
+                psFree(zp);
+                return false;
+            }
+            psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
+                                                   transThresh, starRej, starSys);
+            for (int i = 0; i < num; i++) {
+                fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
+            }
+            psFree(trans);
+        }
+#endif
+    }
+
 
 #if 0
@@ -216 +273 @@
 #endif
 
-#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(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 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);
+    if (options->photometry) {
+        // Save best matches for future photometry
+        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
+        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 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);
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
+    }
 
     psFree(matches);
-
-    // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
-    // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
-    // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
-    // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
-    // We don't need to know the magnitude zero point for the filter, since it cancels out
-    options->norm = psVectorAlloc(num, PS_TYPE_F32);
-    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];
-        options->norm->data.F32[i] = magCorr;
-        psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr);
-
-        for (int j = 0; j < sources->n; j++) {
-            pmSource *source = sources->data[j]; // Source of interest
-            if (!source) {
-                continue;
-            }
-            source->psfMag += magCorr;
-        }
-    }
-    psFree(trans);
-
-#ifdef TESTING
-    // Double check: all transparencies should be zero
-    {
-        psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
-        if (!matches) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
-            psFree(zp);
-            return false;
-        }
-        psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
-                                               transThresh, starRej, starSys);
-        for (int i = 0; i < num; i++) {
-            fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
-        }
-        psFree(trans);
-    }
-#endif
 
     return true;