Changeset 35535 for branches/eam_branches/ipp-20130419/pswarp/src/pswarpUpdateMetadata.c

Timestamp:

May 8, 2013, 4:41:42 PM (13 years ago)

Author:

eugene

Message:

warping of background model appears to work

File:

• branches/eam_branches/ipp-20130419/pswarp/src/pswarpUpdateMetadata.c (modified) (4 diffs)

branches/eam_branches/ipp-20130419/pswarp/src/pswarpUpdateMetadata.c

@@ -10 +10 @@
 
 #include "pswarp.h"
+#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
 
-// once the output fpa elements have been built, loop over the fpa and generate stats
-// for each readout
-bool pswarpUpdateMetadata (pmFPA *output, pmFPA *input, pmFPA *astrom, pmConfig *config)  {
+// once the output fpa elements have been built, loop over the fpa and generate metadata
+bool pswarpUpdateMetadata (pmFPA *output, pmFPA *skycell, pmFPA *input, pmFPA *astrom, pmConfig *config, bool fullImage)  {
 
     pmFPAview *view = pmFPAviewAlloc(0);
     
+    // check output astrometry mode
+    bool bilevelAstrometry = psMetadataLookupBool (NULL, skycell->analysis, "ASTROMETRY.BILEVEL");
+
     pmChip *chip;
     while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) {
@@ -39 +42 @@
                 if (!readout->data_exists) continue;
     
-                // Set covariance matrix for output
-# if (0)
-                psList *covariances = psMetadataLookupPtr(&mdok, output->analysis, PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices
-                psAssert(covariances, "Should be there");
-                psArray *covars = psListToArray(covariances); // Array of covariance matrices
-                psKernel *covar = psImageCovarianceAverage(covars);
-                psFree(covars);
-                psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
+                if (fullImage) {
+                    // Set covariance matrix for output
+                    bool status = false;
+                    psList *covariances = psMetadataLookupPtr(&status, readout->analysis, PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices
+                    psAssert(covariances, "Should be there");
+                    psArray *covars = psListToArray(covariances); // Array of covariance matrices
+                    psKernel *covar = psImageCovarianceAverage(covars);
+                    psFree(covars);
+                    psMetadataRemoveKey(readout->analysis, PSWARP_ANALYSIS_COVARIANCES);
 
-                // Correct covariance matrix scale for the mean (square root of the) Jacobian
-                double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
-                int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
-                jacobian /= goodPixels;
-                output->covariance = psImageCovarianceScale(covar, jacobian);
-                psFree(covar);
+                    // Correct covariance matrix scale for the mean (square root of the) Jacobian
+                    double jacobian = psMetadataLookupF64(NULL, readout->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
+                    int goodPixels = psMetadataLookupS32(NULL, readout->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
+                    jacobian /= goodPixels;
+                    readout->covariance = psImageCovarianceScale(covar, jacobian);
+                    psFree(covar);
 
-                if (output->variance) {
-                  psImageCovarianceTransfer(output->variance, output->covariance);
+                    if (readout->variance) {
+                        psImageCovarianceTransfer(readout->variance, readout->covariance);
+                    }
                 }
-# endif
 
-                {
-                  // Add MD5 information for readout
-                  psString headerName = NULL; ///< Header name for MD5
-                  psVector *md5 = psImageMD5(readout->image); ///< md5 hash
-                  psString md5string = psMD5toString(md5); ///< String
-
-                  pmHDU *hdu = pmHDUFromReadout(readout);
-                  psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
-                  psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, "Image MD5", md5string);
-                  psFree(md5);
-                  psFree(md5string);
-                  psFree(headerName);
-                }
+                // Add MD5 information for readout
+                psString headerName = NULL; ///< Header name for MD5
+                psVector *md5 = psImageMD5(readout->image); ///< md5 hash
+                psString md5string = psMD5toString(md5); ///< String
+                
+                pmHDU *hdu = pmHDUFromReadout(readout);
+                psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
+                psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, "Image MD5", md5string);
+                psFree(md5);
+                psFree(md5string);
+                psFree(headerName);
             }
 
-            psList *inputCells = psMetadataLookupPtr (NULL, cell->analysis, "INPUT.CELLS");
-            if (!pmConceptsAverageCells(cell, inputCells, NULL, NULL, false)) {
-                psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
-                psFree(view);
-                return false;
+            if (fullImage) {
+                psList *inputCells = psMetadataLookupPtr (NULL, cell->analysis, "INPUT.CELLS");
+                if (!pmConceptsAverageCells(cell, inputCells, NULL, NULL, false)) {
+                    psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
+                    psFree(view);
+                    return false;
+                }
             }
 
@@ -85 +89 @@
             psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); ///< Trim section
             trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
+        }
 
+        // update astrometry headers
+        pmHDU *hdu = pmHDUFromChip(chip); ///< HDU for source
+        if (!hdu || !hdu->header) {
+            psError(PM_ERR_PROG, false, "Unable to find header for output.");
+            psFree(view);
+            return false;
+        }
+
+        if (bilevelAstrometry) {
+            if (!pmAstromWriteBilevelChip(hdu->header, chip, WCS_NONLIN_TOL)) {
+                psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for skycell.");
+                psFree(view);
+                return false;
+            }
+        } else {
+            // we use a default FPA pixel scale of 1.0
+            if (!pmAstromWriteWCS(hdu->header, output, chip, WCS_NONLIN_TOL)) {
+                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell.");
+                psFree(view);
+                return false;
+            }
         }
     }
@@ -93 +119 @@
         return false;
     }
-
+    
     // Update ZP from the astrometry
-    {
+    if (fullImage) {
         psMetadataItem *item = psMetadataLookup(output->concepts, "FPA.ZP");
         item->data.F32 = psMetadataLookupF32(NULL, astrom->concepts, "FPA.ZP");
     }
 
-    // need to update the chip and fpa level astrometry appropriately, see
-    // code in psastro
-    
-# if (0)
-    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
-        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
-        return false;
+    // apply the bilevel astrometry elements to the target
+    if (bilevelAstrometry) {
+        pmHDU *phu = pmFPAviewThisPHU(view, output); ///< Astrometry PHU
+        if (!phu->header) {
+            phu->header = psMetadataAlloc ();
+        }
+        if (!pmAstromWriteBilevelMosaic(phu->header, output, WCS_NONLIN_TOL)) {
+            psError(psErrorCodeLast(), false, "Unable to write bilevel mosaic astrometry for skycell.");
+            psFree(view);
+            return false;
+        }
     }
-# endif
 
     return true;