Index: trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- trunk/psModules/src/objects/pmPSF_IO.c (revision 12949) +++ trunk/psModules/src/objects/pmPSF_IO.c (revision 13034) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-04-21 19:47:14 $ + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-04-26 01:20:29 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -44,4 +44,488 @@ #include "pmSourceIO.h" +bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + if (!pmFPAWritePSFmodel (fpa, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write PSF for fpa"); + return false; + } + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + if (!pmChipWritePSFmodel (chip, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write PSF for chip"); + return false; + } + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + if (!pmCellWritePSFmodel (cell, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write PSF for cell"); + return false; + } + return true; + } + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + if (!pmReadoutWritePSFmodel (readout, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write PSF for readout"); + return false; + } + return true; +} + +// read in all chip-level PSFmodel files for this FPA +bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->chips->data[i]; + if (!pmChipWritePSFmodel (chip, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write PSF for %dth chip", i); + return false; + } + } + return true; +} + +// read in all cell-level PSFmodel files for this chip +bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + + for (int i = 0; i < chip->cells->n; i++) { + + pmCell *cell = chip->cells->data[i]; + if (!pmCellWritePSFmodel (cell, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write PSF for %dth cell", i); + return false; + } + } + return true; +} + +// read in all readout-level PSFmodel files for this cell +bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + + for (int i = 0; i < cell->readouts->n; i++) { + + pmReadout *readout = cell->readouts->data[i]; + if (!pmReadoutWritePSFmodel (readout, view, file, config)) { + psError(PS_ERR_IO, false, "Failed to write PSF for %dth readout", i); + return false; + } + } + return true; +} + +// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, +// and header + image for the PSF residual images +bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + bool status; + + // select the psf of interest + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + if (!psf) { + psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout"); + return false; + } + + // lookup the EXTNAME values used for table data and residual image segments + char *rule = NULL; + + // Menu of EXTNAME rules + psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); + if (!menu) { + psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); + return false; + } + // EXTNAME for table data + rule = psMetadataLookupStr(&status, menu, "PSF.TABLE"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config"); + return false; + } + char *tableName = pmFPAfileNameFromRule (rule, file, view); + + // write the PSF model parameters in a FITS table + + // we need to write a header for the table, + psMetadata *header = psMetadataAlloc(); + + char *modelName = pmModelGetType (psf->type); + psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName); + + psMetadataAddBool (header, PS_LIST_TAIL, "POISSON", 0, "Use Poisson errors in fits?", psf->poissonErrors); + + int nPar = pmModelParameterCount (psf->type) ; + psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); + + // save the dimensions of each parameter + for (int i = 0; i < nPar; i++) { + char name[9]; + psPolynomial2D *poly = psf->params_NEW->data[i]; + if (poly == NULL) continue; + snprintf (name, 9, "PAR%02d_NX", i); + psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX); + snprintf (name, 9, "PAR%02d_NY", i); + psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY); + } + + // other required information describing the PSF + psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid); + psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid); + psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ", 0, "chi-square for fit", psf->chisq); + psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", 0, "number of stars used to measure PSF", psf->nPSFstars); + + // XXX can we drop this now? + psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); + + // build a FITS table of the PSF parameters + psArray *psfTable = psArrayAllocEmpty (100); + for (int i = 0; i < nPar; i++) { + psPolynomial2D *poly = psf->params_NEW->data[i]; + if (poly == NULL) continue; // skip unset parameters (eg, XPOS) + 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->mask[ix][iy]); + + psArrayAdd (psfTable, 100, row); + psFree (row); + } + } + } + + // write an empty FITS segment if we have no PSF information + if (psfTable->n == 0) { + // XXX this is probably an error (if we have a PSF, how do we have no data?) + psFitsWriteBlank (file->fits, header, tableName); + } else { + psTrace ("pmFPAfile", 5, "writing psf data %s\n", tableName); + if (!psFitsWriteTable (file->fits, header, psfTable, tableName)) { + psError(PS_ERR_IO, false, "writing psf table data %s\n", tableName); + psFree (tableName); + psFree (psfTable); + psFree (header); + return false; + } + } + psFree (tableName); + psFree (psfTable); + psFree (header); + + // EXTNAME for residual images + rule = psMetadataLookupStr(&status, menu, "PSF.RESID"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config"); + return false; + } + char *imageName = pmFPAfileNameFromRule (rule, file, view); + + // write the residual images (3D) + header = psMetadataAlloc (); + if (psf->residuals == NULL) { + // set some header keywords to make it clear there are no residuals? + psFitsWriteBlank (file->fits, header, imageName); + psFree (imageName); + psFree (header); + return true; + } + + psMetadataAddS32 (header, PS_LIST_TAIL, "XBIN", 0, "", psf->residuals->xBin); + psMetadataAddS32 (header, PS_LIST_TAIL, "YBIN", 0, "", psf->residuals->yBin); + psMetadataAddS32 (header, PS_LIST_TAIL, "XCENTER", 0, "", psf->residuals->xCenter); + psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter); + + // write the residuals as three planes of the image + // this call creates an extension with NAXIS3 = 3 + if (psf->residuals->Rx) { + // this call creates an extension with NAXIS3 = 3 + psArray *images = psArrayAllocEmpty (3); + psArrayAdd (images, 1, psf->residuals->Ro); + psArrayAdd (images, 1, psf->residuals->Rx); + psArrayAdd (images, 1, psf->residuals->Ry); + + psFitsWriteImageCube (file->fits, header, images, imageName); + + // psFitsWriteImage(file->fits, header, psf->residuals->Ro, 3, imageName); + // psFitsUpdateImage(file->fits, psf->residuals->Rx, 0, 0, 1); + // psFitsUpdateImage(file->fits, psf->residuals->Ry, 0, 0, 2); + } else { + // this call creates an extension with NAXIS3 = 1 + psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, imageName); + } + + psFree (imageName); + psFree (header); + return true; + + // XXX save the growth curve + // XXX save the ApResid fit + +# if (0) + // build a FITS table of the fit to the Aperture Residuals + psArray *apresTable = psArrayAllocEmpty (100); + psPolynomial4D *poly = psf->ApTrend; + for (int ix = 0; ix < poly->nX; ix++) { + for (int iy = 0; iy < poly->nY; iy++) { + + row = psMetadataAlloc (); + 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->mask[ix][iy]); + + psArrayAdd (psfTable, 100, row); + psFree (row); + } + } +# endif +} + +// XXX add in error handling +bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + pmFPAReadPSFmodel (fpa, view, file, config); + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipReadPSFmodel (chip, view, file, config); + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellReadPSFmodel (cell, view, file, config); + return true; + } + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + pmReadoutReadPSFmodel (readout, view, file, config); + return true; +} + +// read in all chip-level PSFmodel files for this FPA +bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->chips->data[i]; + pmChipReadPSFmodel (chip, view, file, config); + } + return true; +} + +// read in all cell-level PSFmodel files for this chip +bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + + for (int i = 0; i < chip->cells->n; i++) { + + pmCell *cell = chip->cells->data[i]; + pmCellReadPSFmodel (cell, view, file, config); + } + return true; +} + +// read in all readout-level PSFmodel files for this cell +bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + + for (int i = 0; i < cell->readouts->n; i++) { + + pmReadout *readout = cell->readouts->data[i]; + pmReadoutReadPSFmodel (readout, view, file, config); + } + return true; +} + +// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, +// and header + image for the PSF residual images +bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + bool status; + char *rule = NULL; + psMetadata *header = NULL; + + // Menu of EXTNAME rules + psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); + if (!menu) { + psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); + return false; + } + // EXTNAME for table data + rule = psMetadataLookupStr(&status, menu, "PSF.TABLE"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config"); + return false; + } + char *tableName = pmFPAfileNameFromRule (rule, file, view); + // EXTNAME for residual images + rule = psMetadataLookupStr(&status, menu, "PSF.RESID"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config"); + return false; + } + char *imageName = pmFPAfileNameFromRule (rule, file, view); + + // move fits pointer to table and read header + // advance to the table data extension + // since we have read the IMAGE header, the TABLE header should exist + if (!psFitsMoveExtName (file->fits, tableName)) { + psAbort("cannot find data extension %s in %s", tableName, file->filename); + } + + // load the PSF model table header + header = psFitsReadHeader (NULL, file->fits); + if (!header) psAbort("cannot read table header"); + + // load the PSF model parameters from the FITS table + char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME"); + pmModelType type = pmModelSetType (modelName); + + bool poissonErrors = psMetadataLookupBool (&status, header, "POISSON"); + + // we determine the PSF parameter polynomials from the MD-defined polynomials + pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL); + + // check the number of expected parameters + int nPar = psMetadataLookupS32 (&status, header, "PSF_NPAR"); + if (nPar != pmModelParameterCount (psf->type)) + psAbort("mismatch model par count"); + + // load the dimensions of each parameter + for (int i = 0; i < nPar; i++) { + char name[9]; + snprintf (name, 9, "PAR%02d_NX", i); + int nXorder = psMetadataLookupS32 (&status, header, name); + snprintf (name, 9, "PAR%02d_NY", i); + int nYorder = psMetadataLookupS32 (&status, header, name); + psf->params_NEW->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder); + } + + // other required information describing the PSF + psf->ApResid = psMetadataLookupF32 (&status, header, "AP_RESID"); + psf->dApResid = psMetadataLookupF32 (&status, header, "AP_ERROR"); + psf->chisq = psMetadataLookupF32 (&status, header, "CHISQ"); + psf->nPSFstars = psMetadataLookupS32 (&status, header, "NSTARS"); + + // XXX can we drop this now? + psf->skyBias = psMetadataLookupF32 (&status, header, "SKY_BIAS"); + + // read the raw table data + psArray *table = psFitsReadTable (file->fits); + + // fill in the matching psf->params entries + 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"); + // XXX sanity check here + + psPolynomial2D *poly = psf->params_NEW->data[iPar]; + + poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); + poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR"); + poly->mask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); + } + psFree (header); + psFree (table); + + // move fits pointer to residual image and read header + // advance to the table data extension + // since we have read the IMAGE header, the TABLE header should exist + if (!psFitsMoveExtName (file->fits, imageName)) { + psAbort("cannot find data extension %s in %s", imageName, file->filename); + } + + header = psFitsReadHeader (NULL, file->fits); + int Naxis = psMetadataLookupS32 (&status, header, "NAXIS"); + if (Naxis != 0) { + + int Nx = psMetadataLookupS32 (&status, header, "NAXIS1"); + int Ny = psMetadataLookupS32 (&status, header, "NAXIS2"); + int Nz = psMetadataLookupS32 (&status, header, "NAXIS3"); + + int xBin = psMetadataLookupS32 (&status, header, "XBIN"); + int yBin = psMetadataLookupS32 (&status, header, "YBIN"); + + int xSize = Nx / xBin; + int ySize = Ny / yBin; + + psf->residuals = pmResidualsAlloc (xSize, ySize, xBin, yBin); + + psf->residuals->xCenter = psMetadataLookupS32 (&status, header, "XCENTER"); + psf->residuals->yCenter = psMetadataLookupS32 (&status, header, "YCENTER"); + + psRegion fullImage = {0, 0, 0, 0}; + psFitsReadImageBuffer(psf->residuals->Ro, file->fits, fullImage, 0); // Desired pixels + if (Nz > 1) { + assert (Nz == 3); + psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1); // Desired pixels + psFitsReadImageBuffer(psf->residuals->Ry, file->fits, fullImage, 2); // Desired pixels + } + } + + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf); + psFree (psf); + + psFree (tableName); + psFree (imageName); + + return true; +} + +/************ old support functions, deprecate? **************/ + psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf) { @@ -124,81 +608,6 @@ } -bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) -{ - - pmFPA *fpa = file->fpa; - - if (view->chip == -1) { - pmFPAWritePSFmodel (fpa, view, file, config); - return true; - } - - if (view->chip >= fpa->chips->n) { - return false; - } - pmChip *chip = fpa->chips->data[view->chip]; - - if (view->cell == -1) { - pmChipWritePSFmodel (chip, view, file, config); - return true; - } - - if (view->cell >= chip->cells->n) { - return false; - } - pmCell *cell = chip->cells->data[view->cell]; - - if (view->readout == -1) { - pmCellWritePSFmodel (cell, view, file, config); - return true; - } - - if (view->readout >= cell->readouts->n) { - return false; - } - pmReadout *readout = cell->readouts->data[view->readout]; - - pmReadoutWritePSFmodel (readout, view, file, config); - return true; -} - -// read in all chip-level PSFmodel files for this FPA -bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) -{ - - for (int i = 0; i < fpa->chips->n; i++) { - - pmChip *chip = fpa->chips->data[i]; - pmChipWritePSFmodel (chip, view, file, config); - } - return true; -} - -// read in all cell-level PSFmodel files for this chip -bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) -{ - - for (int i = 0; i < chip->cells->n; i++) { - - pmCell *cell = chip->cells->data[i]; - pmCellWritePSFmodel (cell, view, file, config); - } - return true; -} - -// read in all readout-level PSFmodel files for this cell -bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) -{ - - for (int i = 0; i < cell->readouts->n; i++) { - - pmReadout *readout = cell->readouts->data[i]; - pmReadoutWritePSFmodel (readout, view, file, config); - } - return true; -} - // read in all readout-level Objects files for this cell -bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +bool pmReadoutWritePSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { bool status; @@ -228,83 +637,6 @@ } - - -bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) -{ - - pmFPA *fpa = file->fpa; - - if (view->chip == -1) { - pmFPAReadPSFmodel (fpa, view, file, config); - return true; - } - - if (view->chip >= fpa->chips->n) { - return false; - } - pmChip *chip = fpa->chips->data[view->chip]; - - if (view->cell == -1) { - pmChipReadPSFmodel (chip, view, file, config); - return true; - } - - if (view->cell >= chip->cells->n) { - return false; - } - pmCell *cell = chip->cells->data[view->cell]; - - if (view->readout == -1) { - pmCellReadPSFmodel (cell, view, file, config); - return true; - } - - if (view->readout >= cell->readouts->n) { - return false; - } - pmReadout *readout = cell->readouts->data[view->readout]; - - pmReadoutReadPSFmodel (readout, view, file, config); - return true; -} - -// read in all chip-level PSFmodel files for this FPA -bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) -{ - - for (int i = 0; i < fpa->chips->n; i++) { - - pmChip *chip = fpa->chips->data[i]; - pmChipReadPSFmodel (chip, view, file, config); - } - return true; -} - -// read in all cell-level PSFmodel files for this chip -bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) -{ - - for (int i = 0; i < chip->cells->n; i++) { - - pmCell *cell = chip->cells->data[i]; - pmCellReadPSFmodel (cell, view, file, config); - } - return true; -} - -// read in all readout-level PSFmodel files for this cell -bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) -{ - - for (int i = 0; i < cell->readouts->n; i++) { - - pmReadout *readout = cell->readouts->data[i]; - pmReadoutReadPSFmodel (readout, view, file, config); - } - return true; -} - // read in all readout-level Objects files for this cell -bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +bool pmReadoutReadPSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) {