Index: trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- trunk/psModules/src/objects/pmPSF_IO.c (revision 30031) +++ trunk/psModules/src/objects/pmPSF_IO.c (revision 31451) @@ -63,4 +63,9 @@ bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header); +bool pmPSFmodelRead_ApTrend (pmPSF *psf, pmFPAfile *file); +bool pmPSFmodelWrite_ApTrend (pmFPAfile *file, pmPSF *psf); + +bool pmPSFmodelRead_GrowthCurve (pmPSF *psf, pmFPAfile *file); +bool pmPSFmodelWrite_GrowthCurve (pmFPAfile *file, pmPSF *psf); bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file) @@ -197,4 +202,80 @@ 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,9 +484,9 @@ 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); @@ -428,45 +509,18 @@ psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); + float PSF_APERTURE = psMetadataLookupF32(&status, roAnalysis, "PSF_APERTURE"); + if (status) { + psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_APERTURE", PS_DATA_F32, "aperture for psf objects", PSF_APERTURE); + } + float PSF_FIT_RADIUS = psMetadataLookupF32(&status, roAnalysis, "PSF_FIT_RADIUS"); + if (status) { + psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_DATA_F32, "aperture for psf objects", PSF_FIT_RADIUS); + } + // build a FITS table of the PSF parameters psArray *psfTable = psArrayAllocEmpty (100); 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); } @@ -547,20 +601,17 @@ } + if (!pmPSFmodelWrite_ApTrend(file, psf)) { + psError(psErrorCodeLast(), false, "Unable to write PSF ApTrend"); + return false; + } + + if (!pmPSFmodelWrite_GrowthCurve(file, psf)) { + psError(psErrorCodeLast(), false, "Unable to write PSF Growth Curve"); + 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; @@ -651,4 +702,6 @@ } + + // if this file needs to have a PHU written out, write one bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config) @@ -934,4 +987,17 @@ psf->skyBias = psMetadataLookupF32 (&status, header, "SKY_BIAS"); + if (roAnalysis) { + float PSF_APERTURE = psMetadataLookupF32(&status, header, "PSF_APERTURE"); + if (status) { + psMetadataAddF32 (roAnalysis, PS_LIST_TAIL, "PSF_APERTURE", PS_DATA_F32, "aperture for psf objects", PSF_APERTURE); + } + float PSF_FIT_RADIUS = psMetadataLookupF32(&status, header, "PSF_FIT_RADIUS"); + if (status) { + psMetadataAddF32 (roAnalysis, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_DATA_F32, "aperture for psf objects", PSF_FIT_RADIUS); + } + } else { + psWarning ("unable to read PSF_APERTURE or PSF_FIT_RADIUS"); + } + // read the raw table data psArray *table = psFitsReadTable (file->fits); @@ -945,7 +1011,6 @@ 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]; @@ -955,26 +1020,5 @@ } - 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); @@ -1063,4 +1107,14 @@ } + if (!pmPSFmodelRead_ApTrend (psf, file)) { + psError(psErrorCodeLast(), false, "Unable to read PSF ApTrend data."); + return false; + } + + if (!pmPSFmodelRead_GrowthCurve(psf, file)) { + psError(psErrorCodeLast(), false, "Unable to read PSF Growth Curve"); + return false; + } + psMetadataAdd (chipAnalysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf); psFree (psf); @@ -1070,4 +1124,217 @@ 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; +} + +// write aperture trend to a FITS table +bool pmPSFmodelWrite_GrowthCurve (pmFPAfile *file, pmPSF *psf) { + + pmGrowthCurve *growth = psf->growth; + if (growth == NULL) { + psWarning ("no PSF Growth Curve to write out, skipping"); + return true; + } + + // we need to write a header for the table, + psMetadata *header = psMetadataAlloc(); + + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_MIN_RAD", 0, "", growth->radius->data.F32[0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_MAX_RAD", 0, "", growth->maxRadius); + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_REF_RAD", 0, "", growth->refRadius); + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_AP_LOSS", 0, "", growth->apLoss); + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_AP_REF", 0, "", growth->apRef); + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_FIT_MAG", 0, "", growth->fitMag); + + // build a FITS table of the ApTrend (only 1) + psArray *table = psArrayAllocEmpty (100); + for (int i = 0; i < growth->apMag->n; i++) { + psMetadata *row = psMetadataAlloc (); + psMetadataAddF32 (row, PS_LIST_TAIL, "RADIUS", 0, "", growth->radius->data.F32[i]); + psMetadataAddF32 (row, PS_LIST_TAIL, "AP_MAG", 0, "", growth->apMag->data.F32[i]); + psArrayAdd (table, 100, row); + psFree (row); + } + + // 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 Growth Curve data %s\n", "GROWTH_CURVE"); + if (!psFitsWriteTable(file->fits, header, table, "GROWTH_CURVE")) { + psError(psErrorCodeLast(), false, "Error writing psf table data %s\n", "GROWTH_CURVE"); + psFree(table); + psFree(header); + return false; + } + + psFree (table); + psFree (header); + return true; +} + +// read aperture trend to a FITS table +bool pmPSFmodelRead_GrowthCurve (pmPSF *psf, pmFPAfile *file) { + + bool status; + + // move fits pointer to AP_TREND section + // advance to the table data extension + if (!psFitsMoveExtNameClean (file->fits, "GROWTH_CURVE")) { + psWarning ("no Growth Curve data in PSF file, skipping"); + return true; + } + + psMetadata *header = psFitsReadHeader (NULL, file->fits); + if (!header) { + psError(psErrorCodeLast(), false, "Unable to read GROWTH_CURVE header."); + return false; + } + + // read the raw table data + psArray *table = psFitsReadTable (file->fits); + if (!table) { + psError(psErrorCodeLast(), false, "Unable to read GROWTH_CURVE table."); + psFree(header); + return false; + } + + float minRadius = psMetadataLookupF32 (&status, header, "GROWTH_MIN_RAD"); if (!status) return false; + float maxRadius = psMetadataLookupF32 (&status, header, "GROWTH_MAX_RAD"); if (!status) return false; + float refRadius = psMetadataLookupF32 (&status, header, "GROWTH_REF_RAD"); if (!status) return false; + + psf->growth = pmGrowthCurveAlloc(minRadius, maxRadius, refRadius); + + psf->growth->apLoss = psMetadataLookupF32 (&status, header, "GROWTH_AP_LOSS"); if (!status) return false; + psf->growth->apRef = psMetadataLookupF32 (&status, header, "GROWTH_AP_REF"); if (!status) return false; + psf->growth->fitMag = psMetadataLookupF32 (&status, header, "GROWTH_FIT_MAG"); if (!status) return false; + + // fill in the matching psf->params entries + for (int i = 0; i < table->n; i++) { + psMetadata *row = table->data[i]; + psf->growth->apMag->data.F32[i] = psMetadataLookupF32 (&status, row, "AP_MAG"); if (!status) return false; + } + + psFree (header); + psFree (table); return true; }