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

Timestamp:

May 18, 2010, 4:11:53 PM (16 years ago)

Author:

eugene

Message:

merge changes from branches/eam_branches/psphot,psModules.20100506 (finish basic psphotStack)

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src/psphotStackImageLoop.c (modified) (7 diffs)

trunk/psphot/src/psphotStackImageLoop.c

@@ -1 +1 @@
 # include "psphotStandAlone.h"
-
-# define ESCAPE(MESSAGE) { \
-  psError(PSPHOT_ERR_DATA, false, MESSAGE); \
-  psFree (view); \
-  return false; \
-}
+#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
+
+# define ESCAPE(MESSAGE) {                              \
+        psError(PSPHOT_ERR_DATA, false, MESSAGE);       \
+        psFree (view);                                  \
+        return false;                                   \
+    }
+
+bool GetAstrometryFPA (pmConfig *config, pmFPAview *view);
+bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry);
+bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry);
 
 bool psphotStackImageLoop (pmConfig *config) {
@@ -14 +19 @@
     pmReadout *readout;
 
-    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
-    if (!status) {
+    pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
+    pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
+    pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+
+    if (!input) {
         psError(PSPHOT_ERR_PROG, false, "Can't find input data!");
         return false;
@@ -24 +32 @@
     // 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
@@ -41 +51 @@
                 if (! readout->data_exists) { continue; }
 
-# if (0)                
-                // uncomment to generate matched psfs
+                // PSF matching
                 if (!psphotStackMatchPSFs (config, view)) {
                     psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
@@ -48 +57 @@
                     return false;
                 }
-# endif
 
                 // XXX for now, we assume there is only a single chip in the PHU:
@@ -69 +77 @@
             }
         }
+
+        GetAstrometryChip(config, view, bilevelAstrometry);
+
         // save output which is saved at the chip level
         if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot.");
     }
+
+    SetAstrometryFPA(config, view, bilevelAstrometry);
+    
     // save output which is saved at the fpa level
     if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot.");
@@ -87 +101 @@
 */
 
+# define UPDATE_HEADER 0
+
+bool GetAstrometryFPA (pmConfig *config, pmFPAview *view) {
+
+    bool status = false;
+
+    bool foundAstrometry = false;
+    bool bilevelAstrometry = false;
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // 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?");
+
+        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
+        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
+        pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+        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;
+        }
+
+        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
+        if (!ctype) {
+            psWarning("Error in WCS keywords for input FPA entry %d", i);
+            continue;
+        }
+
+        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.");
+            }
+        }
+    }
+    return bilevelAstrometry;
+}
+
+bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) {
+
+    bool status = false;
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // 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?");
+
+        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
+        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
+        pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+        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
+
+# if (UPDATE_HEADER)
+        pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa);
+        if (!outHDU) {
+            pmFPAAddSourceFromView(output->fpa, "name", view, output->format);
+            outHDU = pmFPAviewThisHDU (view, output->fpa);
+            psAssert (outHDU, "failed to make HDU");
+        }
+# endif
+
+        if (bilevelAstrometry) {
+            if (!pmAstromReadBilevelChip (outChip, inHDU->header)) {
+                psWarning("Unable to read bilevel chip astrometry for input FPA.");
+                continue;
+            }
+# if (UPDATE_HEADER)
+            if (!pmAstromWriteBilevelChip(outHDU->header, outChip, WCS_NONLIN_TOL)) {
+                psWarning("Unable to generate WCS header.");
+                continue;
+            }
+# endif
+        } 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 (UPDATE_HEADER)
+            if (UPDATE_HEADER && !pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) {
+                psWarning("Unable to generate WCS header.");
+                continue;
+            }
+# endif
+        }
+    }
+
+    return true;
+}
+
+bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) {
+
+    bool status = false;
+
+    if (!bilevelAstrometry) return true;
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // 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?");
+
+# if (UPDATE_HEADER)
+        pmHDU *PHU = pmFPAviewThisPHU(view, output->fpa);
+        if (!PHU) {
+            pmFPAAddSourceFromView(output->fpa, "name", view, output->format);
+            PHU = pmFPAviewThisPHU (view, output->fpa);
+            psAssert (PHU, "failed to make PHU");
+        }
+
+        if (!pmAstromWriteBilevelMosaic(PHU->header, output->fpa, WCS_NONLIN_TOL)) {
+            psWarning("Unable to generate WCS header.");
+        }
+# endif
+    }
+
+    return true;
+}