Changeset 31339 for branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF_IO.c

Timestamp:

Apr 21, 2011, 4:21:45 PM (15 years ago)

Author:

eugene

Message:

add ApTrend I/O code

File:

• branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF_IO.c (modified) (10 diffs)

branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF_IO.c

@@ -63 +63 @@
 
 bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header);
+bool pmPSFmodelRead_ApTrend (pmPSF *psf, pmFPAfile *file);
+bool pmPSFmodelWrite_ApTrend (pmFPAfile *file, pmPSF *psf);
 
 bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file)
@@ -197 +199 @@
         psError(psErrorCodeLast(), false, "Failed to write PSF for chip");
         return false;
+    }
+    return true;
+}
+
+// XXX we save the model term identifiers (item) as S32, but they probably should be more flexible
+bool pmTrend2DtoTable (psArray *table, pmTrend2D *trend, char *label, int item) {
+
+    if (trend == NULL) return true; 
+
+    if (trend->mode == PM_TREND_MAP) {
+        // write the image components into a table: this is needed because they may each be a different size
+        psImageMap *map = trend->map;
+        for (int ix = 0; ix < map->map->numCols; ix++) {
+            for (int iy = 0; iy < map->map->numRows; iy++) {
+                psMetadata *row = psMetadataAlloc ();
+                psMetadataAddS32 (row, PS_LIST_TAIL, label,        0, "", item);
+                psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
+                psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
+                psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
+
+                psArrayAdd (table, 100, row);
+                psFree (row);
+            }
+        }
+    } else {
+        // write the polynomial components into a table
+        psPolynomial2D *poly = trend->poly;
+        for (int ix = 0; ix <= poly->nX; ix++) {
+            for (int iy = 0; iy <= poly->nY; iy++) {
+                psMetadata *row = psMetadataAlloc ();
+                psMetadataAddS32 (row, PS_LIST_TAIL, label,        0, "", item);
+                psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
+                psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
+                psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->coeffMask[ix][iy]);
+
+                psArrayAdd (table, 100, row);
+                psFree (row);
+            }
+        }
+    }
+    return true;
+}
+
+// extra trend2D elements from a row
+bool pmTrend2DfromTableRow (pmTrend2D *trend, psMetadata *row) {
+
+    bool status = false;
+
+    int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
+    int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
+
+    if (trend->mode == PM_TREND_MAP) {
+        psImageMap *map = trend->map;
+        assert (map);
+        assert (map->map);
+        assert (map->error);
+        assert (xPow >= 0);
+        assert (yPow >= 0);
+        assert (xPow < map->map->numCols);
+        assert (yPow < map->map->numRows);
+        map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
+        map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
+    } else {
+        psPolynomial2D *poly = trend->poly;
+        assert (poly);
+        assert (xPow >= 0);
+        assert (yPow >= 0);
+        assert (xPow <= poly->nX);
+        assert (yPow <= poly->nY);
+        poly->coeff[xPow][yPow]     = psMetadataLookupF32 (&status, row, "VALUE");
+        poly->coeffErr[xPow][yPow]  = psMetadataLookupF32 (&status, row, "ERROR");
+        poly->coeffMask[xPow][yPow] = psMetadataLookupU8  (&status, row, "MASK");
     }
     return true;
@@ -403 +481 @@
 
             if (trend->mode == PM_TREND_MAP) {
-              nX = trend->map->map->numCols;
-              nY = trend->map->map->numRows;
+                nX = trend->map->map->numCols;
+                nY = trend->map->map->numRows;
             } else {
-              nX = trend->poly->nX;
-              nY = trend->poly->nY;
+                nX = trend->poly->nX;
+                nY = trend->poly->nY;
             }
             snprintf (name, 9, "PAR%02d_NX", i);
@@ -441 +519 @@
         for (int i = 0; i < nPar; i++) {
             pmTrend2D *trend = psf->params->data[i];
-            if (trend == NULL) continue; // skip unset parameters (eg, XPOS)
-
-            if (trend->mode == PM_TREND_MAP) {
-                // write the image components into a table: this is needed because they may each be a different size
-                psImageMap *map = trend->map;
-                for (int ix = 0; ix < map->map->numCols; ix++) {
-                    for (int iy = 0; iy < map->map->numRows; iy++) {
-                        psMetadata *row = psMetadataAlloc ();
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
-                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
-                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
-                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
-
-                        psArrayAdd (psfTable, 100, row);
-                        psFree (row);
-                    }
-                }
-            } else {
-                // write the polynomial components into a table
-                psPolynomial2D *poly = trend->poly;
-                for (int ix = 0; ix <= poly->nX; ix++) {
-                    for (int iy = 0; iy <= poly->nY; iy++) {
-                        psMetadata *row = psMetadataAlloc ();
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
-                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
-                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
-                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->coeffMask[ix][iy]);
-
-                        psArrayAdd (psfTable, 100, row);
-                        psFree (row);
-                    }
-                }
-            }
+            pmTrend2DtoTable (psfTable, trend, "MODEL_TERM", i);
         }
 
@@ -556 +598 @@
     }
 
+    if (!pmPSFmodelWrite_ApTrend(file, psf)) {
+        psError(psErrorCodeLast(), false, "Unable to write PSF ApTrend");
+        return false;
+    }
+
     // write a representation of the psf model
     {
         psMetadata *header = psMetadataAlloc ();
-
-        if (0) {
-            // set some header keywords to make it clear there are no residuals?
-            if (!psFitsWriteBlank (file->fits, header, residName)) {
-                psError(psErrorCodeLast(), false, "Unable to write blank PSF residuals.");
-                psFree(residName);
-                psFree(header);
-                return false;
-            }
-            psFree (residName);
-            psFree (header);
-            return true;
-        }
 
         int DX = 65;
@@ -660 +694 @@
 }
 
+
+
 // if this file needs to have a PHU written out, write one
 bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config)
@@ -967 +1003 @@
     for (int i = 0; i < table->n; i++) {
         psMetadata *row = table->data[i];
+
         int iPar = psMetadataLookupS32 (&status, row, "MODEL_TERM");
-        int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
-        int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
 
         pmTrend2D *trend = psf->params->data[iPar];
@@ -977 +1012 @@
         }
 
-        if (trend->mode == PM_TREND_MAP) {
-            psImageMap *map = trend->map;
-            assert (map);
-            assert (map->map);
-            assert (map->error);
-            assert (xPow >= 0);
-            assert (yPow >= 0);
-            assert (xPow < map->map->numCols);
-            assert (yPow < map->map->numRows);
-            map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
-            map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
-        } else {
-            psPolynomial2D *poly = trend->poly;
-            assert (poly);
-            assert (xPow >= 0);
-            assert (yPow >= 0);
-            assert (xPow <= poly->nX);
-            assert (yPow <= poly->nY);
-            poly->coeff[xPow][yPow]     = psMetadataLookupF32 (&status, row, "VALUE");
-            poly->coeffErr[xPow][yPow]  = psMetadataLookupF32 (&status, row, "ERROR");
-            poly->coeffMask[xPow][yPow] = psMetadataLookupU8  (&status, row, "MASK");
-        }
+        pmTrend2DfromTableRow(trend, row);
     }
     psFree (header);
@@ -1085 +1099 @@
     }
 
+    if (!pmPSFmodelRead_ApTrend (psf, file)) {
+        psError(psErrorCodeLast(), false, "Unable to read PSF ApTrend data.");
+        return false;
+    }
+
     psMetadataAdd (chipAnalysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
     psFree (psf);
@@ -1092 +1111 @@
     psFree (header);
 
+    return true;
+}
+
+// write aperture trend to a FITS table
+bool pmPSFmodelWrite_ApTrend (pmFPAfile *file, pmPSF *psf) {
+
+    pmTrend2D *trend = psf->ApTrend;
+    if (trend == NULL) { 
+        psWarning ("no PSF ApTrend to write out, skipping");
+        return true; 
+    }
+
+    // we need to write a header for the table,
+    psMetadata *header = psMetadataAlloc();
+
+    int nX = 0, nY = 0;
+    if (trend->mode == PM_TREND_MAP) {
+        nX = trend->map->map->numCols;
+        nY = trend->map->map->numRows;
+    } else {
+        nX = trend->poly->nX;
+        nY = trend->poly->nY;
+    }
+    psMetadataAddS32 (header, PS_LIST_TAIL, "TREND_NX", 0, "", nX);
+    psMetadataAddS32 (header, PS_LIST_TAIL, "TREND_NY", 0, "", nY);
+    char *modeName = pmTrend2DModeToString (trend->mode);
+    psMetadataAddStr (header, PS_LIST_TAIL, "TREND_MD", 0, "", modeName);
+    psFree (modeName);
+
+    // build a FITS table of the ApTrend (only 1)
+    psArray *table = psArrayAllocEmpty (100);
+    pmTrend2DtoTable (table, trend, "APTREND", 0);
+
+    // write an empty FITS segment if we have no PSF information
+    if (table->n == 0) {
+        psError(PM_ERR_PROG, true, "No PSF data to write.");
+        psFree(table);
+        psFree(header);
+        return false;
+    } 
+
+    psTrace ("pmFPAfile", 5, "writing psf ApTrend data %s\n", "AP_TREND");
+    if (!psFitsWriteTable(file->fits, header, table, "AP_TREND")) {
+        psError(psErrorCodeLast(), false, "Error writing psf table data %s\n", "AP_TREND");
+        psFree(table);
+        psFree(header);
+        return false;
+    }
+
+    psFree (table);
+    psFree (header);
+    return true;
+}
+
+// read aperture trend to a FITS table
+bool pmPSFmodelRead_ApTrend (pmPSF *psf, pmFPAfile *file) {
+
+    bool status;
+
+    // move fits pointer to AP_TREND section
+    // advance to the table data extension
+    if (!psFitsMoveExtNameClean (file->fits, "AP_TREND")) {
+        psWarning ("no Aperture Trend data in PSF file, skipping");
+        return true;
+    }
+
+    psMetadata *header = psFitsReadHeader (NULL, file->fits);
+    if (!header) {
+        psError(psErrorCodeLast(), false, "Unable to read AP_TREND header.");
+        return false;
+    }
+        
+    // read the raw table data
+    psArray *table = psFitsReadTable (file->fits);
+    if (!table) {
+        psError(psErrorCodeLast(), false, "Unable to read AP_TREND table.");
+        psFree(header);
+        return false;
+    }
+
+    // XXX allow user to set this optionally?
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    psImageBinning *binning = psImageBinningAlloc();
+    binning->nXfine = psf->fieldNx;
+    binning->nYfine = psf->fieldNy;
+    binning->nXruff = psMetadataLookupS32 (&status, header, "TREND_NX");
+    binning->nYruff = psMetadataLookupS32 (&status, header, "TREND_NY");
+    psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
+    char *modeName  = psMetadataLookupStr (&status, header, "TREND_MD");
+    if (!status) {
+        psError(PM_ERR_PROG, true, "inconsistent PSF header: NX & NY defined for AP TREND, but not MD");
+        psFree (header);
+        psFree (stats);
+        psFree (table);
+        return false;
+    }
+    pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName);
+    if (psfTrendMode == PM_TREND_NONE) {
+        psfTrendMode = PM_TREND_POLY_ORD;
+    }
+
+    // measure Trend2D for the current spatial scale
+    pmTrend2D *apTrend = pmTrend2DNoImageAlloc (PM_TREND_MAP, binning, stats);
+
+    // fill in the matching psf->params entries
+    for (int i = 0; i < table->n; i++) {
+        psMetadata *row = table->data[i];
+        pmTrend2DfromTableRow(apTrend, row);
+    }
+    psf->ApTrend = apTrend;
+
+    psFree (binning);
+    psFree (header);
+    psFree (stats);
+    psFree (table);
     return true;
 }