Changeset 31154 for trunk/psphot/src/psphotStackImageLoop.c

Timestamp:

Apr 4, 2011, 1:12:39 PM (15 years ago)

Author:

eugene

Message:

use the smoothed image for footprint culling; use a more stringent culling for saturated stars; more tweaks to get sat stars to be fitted; various updates to psphotStack; unify psphotImageLoop varients; be a bit careful about number of stars in a PSF clump

Location:

trunk/psphot

Files:

• . (modified) (2 props)
• src/psphotStackImageLoop.c (modified) (11 diffs)

trunk/psphot

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM

trunk/psphot/src/psphotStackImageLoop.c

@@ -1 +1 @@
 # include "psphotStandAlone.h"
-#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
 
 # define ESCAPE(MESSAGE) {                              \
@@ -8 +7 @@
     }
 
-bool GetAstrometryFPA (pmConfig *config, pmFPAview *view);
-bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry);
-bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry);
+// XXX this implementation is not smart about multi-level astrometry headers
+bool UpdateHeadersForFPA (pmConfig *config, pmFPAview *view);
+bool UpdateHeadersForChip (pmConfig *config, pmFPAview *view);
 
 bool psphotStackImageLoop (pmConfig *config) {
@@ -18 +17 @@
     pmCell *cell;
     pmReadout *readout;
+
+    psMemDump("startloop");
 
     pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
@@ -32 +33 @@
     // XXX for now, just load the full set of images up front
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
-
-    bool bilevelAstrometry = GetAstrometryFPA(config, view);
 
     // for psphot, we force data to be read at the chip level
@@ -51 +50 @@
                 if (! readout->data_exists) { continue; }
 
+                psMemDump("load");
+
                 // PSF matching
                 if (!psphotStackMatchPSFs (config, view)) {
@@ -57 +58 @@
                     return false;
                 }
+                psMemDump("stackmatch");
 
                 // XXX for now, we assume there is only a single chip in the PHU:
@@ -64 +66 @@
                     return false;
                 }
-
+                psMemDump("psphot");
             }
             // drop all versions of the internal files
@@ -78 +80 @@
         }
 
-        GetAstrometryChip(config, view, bilevelAstrometry);
+        UpdateHeadersForChip(config, view);
 
         // save output which is saved at the chip level
         if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot.");
     }
+    psMemDump("doneloop");
 
-    SetAstrometryFPA(config, view, bilevelAstrometry);
-    
+    UpdateHeadersForFPA(config, view);
+
     // save output which is saved at the fpa level
     if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot.");
@@ -93 +96 @@
 
     psFree (view);
+
+    psMemDump("doneoutput");
     return true;
 }
 
-/* 
-   the easiest way to implement this is to assume we can pre-load the full set of images up front.
-   with 5 filters and 6000^2 (image, mask, var = 10 byte per pixel), we need 1.8GB, which is not too bad.
-*/
-
-# define UPDATE_HEADER 1
-
-bool GetAstrometryFPA (pmConfig *config, pmFPAview *view) {
-
-    bool foundAstrometry = false;
-    bool bilevelAstrometry = false;
+bool UpdateHeadersForChip (pmConfig *config, pmFPAview *view) {
 
     int num = psphotFileruleCount(config, "PSPHOT.INPUT");
@@ -122 +117 @@
         psAssert (input, "missing input file");
 
-        // find the FPA phu
-        pmHDU *phu = pmFPAviewThisPHU(view, input->fpa);
-        if (!phu) {
-            psWarning("Unable to read bilevel mosaic astrometry for input FPA entry %d", i);
-            continue;
+        // just copy the input headers to the output headers, then update version info
+        pmChip *inChip = pmFPAviewThisChip(view, input->fpa); ///< Chip in the input
+        pmHDU *inHDU = pmFPAviewThisHDU (view, input->fpa);
+
+        pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
+        pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa);
+        if (!outHDU) {
+            pmFPAAddSourceFromView(output->fpa, view, output->format);
+            outHDU = pmFPAviewThisHDU (view, output->fpa);
+            psAssert (outHDU, "failed to make HDU");
         }
-
-        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
-        if (!ctype) {
-            psWarning("Error in WCS keywords for input FPA entry %d", i);
-            continue;
+        if (!outHDU->header) {
+            outHDU->header = psMetadataCopy(NULL, inHDU->header);
         }
-
-        if (!foundAstrometry) {
-            bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
-        } else {
-            if (bilevelAstrometry != !strcmp (&ctype[4], "-DIS")) {
-                psAbort("astrometry type mis-match");
-            }
-        }
-
-        if (bilevelAstrometry) {
-            // update the output structures
-            if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) {
-                psWarning("Unable to read bilevel mosaic astrometry for input FPA.");
-            }
-        }
+        outChip->toFPA = psMemIncrRefCounter(inChip->toFPA);
+        outChip->fromFPA = psMemIncrRefCounter(inChip->fromFPA);
     }
-    return bilevelAstrometry;
+    return true;
 }
 
-bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) {
+bool UpdateHeadersForFPA (pmConfig *config, pmFPAview *view) {
 
     int num = psphotFileruleCount(config, "PSPHOT.INPUT");
@@ -169 +153 @@
         psAssert (input, "missing input file");
 
-        // Need to update the output for astrometry.  Read WCS data from the corresponding
-        // header and save in the output fpa
-        pmHDU *inHDU = pmFPAviewThisHDU (view, input->fpa);
-        pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
+        output->fpa->toTPA = psMemIncrRefCounter(input->fpa->toTPA);
+        output->fpa->fromTPA = psMemIncrRefCounter(input->fpa->fromTPA);
+        output->fpa->toSky = psMemIncrRefCounter(input->fpa->toSky);
 
-        pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa);
-        if (!outHDU) {
-            pmFPAAddSourceFromView(output->fpa, view, output->format);
-            outHDU = pmFPAviewThisHDU (view, output->fpa);
-            psAssert (outHDU, "failed to make HDU");
-        }
-        if (!outHDU->header) {
-          outHDU->header = psMetadataAlloc();
-        }
+        pmConceptsCopyFPA(output->fpa, input->fpa, true, true);
 
-        if (bilevelAstrometry) {
-            if (!pmAstromReadBilevelChip (outChip, inHDU->header)) {
-                psWarning("Unable to read bilevel chip astrometry for input FPA.");
-                continue;
-            }
-            if (!pmAstromWriteBilevelChip(outHDU->header, outChip, WCS_NONLIN_TOL)) {
-                psWarning("Unable to generate WCS header.");
-                continue;
-            }
-        } else {
-            // we use a default FPA pixel scale of 1.0
-            if (!pmAstromReadWCS (output->fpa, outChip, inHDU->header, 1.0)) {
-                psWarning("Unable to read WCS astrometry for input FPA.");
-                continue;
-            }
-            if (!pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) {
-                psWarning("Unable to generate WCS header.");
-                continue;
-            }
-        }
+        // XXX TEST
+        // pmFPASetFileStatus(output->fpa, true);
+        // pmFPASetDataStatus(output->fpa, true);
+        // pmChip *chip = output->fpa->chips->data[0];
+        // pmCell *cell = chip->cells->data[0];
+        // pmReadout *readout = cell->readouts->data[0];
     }
-
     return true;
 }
 
-bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) {
+/* 
+   the easiest way to implement this is to assume we can pre-load the full set of images up front.
+   with 5 filters and 6000^2 (image, mask, var = 10 byte per pixel), we need 1.8GB, which is not too bad.
+*/
 
-    if (!bilevelAstrometry) return true;
-
-    int num = psphotFileruleCount(config, "PSPHOT.INPUT");
-
-    // loop over the available readouts
-    for (int i = 0; i < num; i++) {
-
-        // find the currently selected readout
-        pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest
-        psAssert (output, "missing file?");
-
-        pmHDU *PHU = pmFPAviewThisPHU(view, output->fpa);
-        if (!PHU) {
-            pmFPAAddSourceFromView(output->fpa, view, output->format);
-            PHU = pmFPAviewThisPHU (view, output->fpa);
-            psAssert (PHU, "failed to make PHU");
-        }
-        if (!PHU->header) {
-          PHU->header = psMetadataAlloc();
-        }
-
-        if (!pmAstromWriteBilevelMosaic(PHU->header, output->fpa, WCS_NONLIN_TOL)) {
-            psWarning("Unable to generate WCS header.");
-        }
-    }
-
-    return true;
-}

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM