Changeset 27955

Timestamp:

May 13, 2010, 6:31:45 PM (16 years ago)

Author:

eugene

Message:

update astrometry

Location:

branches/eam_branches/psphot.20100506/src

Files:

• psphotPetrosianStats.c (modified) (1 diff)
• psphotStackImageLoop.c (modified) (5 diffs)

branches/eam_branches/psphot.20100506/src/psphotPetrosianStats.c

@@ -195 +195 @@
     source->extpars->petrosianR90Err    = NAN;
 
-    fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
+    // fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
     psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
 

branches/eam_branches/psphot.20100506/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) {
@@ -16 +21 @@
     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;
 
@@ -28 +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
@@ -71 +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.");
@@ -89 +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;
+}