Index: trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- trunk/psModules/src/objects/pmPSF_IO.c (revision 13414) +++ trunk/psModules/src/objects/pmPSF_IO.c (revision 13526) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-05-18 00:17:12 $ + * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-05-26 02:52:27 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -20,4 +20,7 @@ /* INCLUDE FILES */ /*****************************************************************************/ + +#include +#include /* for strn?casecmp */ #include @@ -212,148 +215,234 @@ } -// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, -// and header + image for the PSF residual images +// for each Readout (ie, analysed image), we write out +// - image header : FITS Image NAXIS = 0 +// - psf table (+header) : FITS Table +// - psf resid (+header) : FITS Image +// if needed, we also write out a PHU blank header bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { bool status; + pmHDU *hdu; + pmHDU *phu; + char *headName, *tableName, *residName; + + // write a PHU? (only if input image is MEF) + // write a header? (only if this is the first readout for cell) + // note that the file->header is set to track the last hdu->header written + // write the data? (always?) + + // get the current header + hdu = pmFPAviewThisHDU (view, file->fpa); + + // if file does not yet have a PHU, attempt to write one to disk + // we only need a PHU if chips->n > 1 and file->fileLevel == FPA + // otherwise, the chip header fills the PHU location + // XXX this code could be placed in a 'pmPSF_WritePHU' function and called + // from pmFPAfileIO.c. + if ((file->fileLevel == PM_FPA_LEVEL_FPA) && (file->fpa->chips->n > 1) && !file->phu) { + + // find the FPA phu + phu = pmFPAviewThisPHU (view, file->fpa); + + // if there is no PHU, this is a single header+image (extension-less) file + // if there is a PHU, write it out as a 'blank' + psMetadata *outhead = psMetadataAlloc(); + if (phu) { + psMetadataCopy (outhead, phu->header); + } else { + pmConfigConformHeader (outhead, file->format); + + psMetadata *fileData = psMetadataLookupMetadata(NULL, file->format, "FILE"); // File information + const char *fpaNameHdr = psMetadataLookupStr(NULL, fileData, "FPA.NAME"); + if (fpaNameHdr && strlen(fpaNameHdr) > 0) { + const char *fpaName = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.NAME"); + psMetadataAddStr(outhead, PS_LIST_TAIL, fpaNameHdr, PS_META_REPLACE, "FPA name", fpaName); + } + } + + psMetadataAddBool (outhead, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true); + psFitsWriteBlank (file->fits, outhead, ""); + file->phu = true; + psTrace ("pmFPAfile", 5, "wrote phu %s (type: %d)\n", file->filename, file->type); + psFree (outhead); + } + + // define the EXTNAME values used for image header, table data, and residual image segments + { + // lookup the EXTNAME values used for table data and image header 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 image header + rule = psMetadataLookupStr(&status, menu, "PSF.HEAD"); + if (!rule) { + psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.HEAD in EXTNAME.RULES in camera.config"); + return false; + } + headName = pmFPAfileNameFromRule (rule, file, view); + + // EXTNAME for table data + rule = psMetadataLookupStr(&status, menu, "PSF.TABLE"); + if (!rule) { + psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config"); + psFree (headName); + return false; + } + tableName = pmFPAfileNameFromRule (rule, file, view); + + // EXTNAME for resid data + rule = psMetadataLookupStr(&status, menu, "PSF.RESID"); + if (!rule) { + psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config"); + psFree (headName); + psFree (tableName); + return false; + } + residName = pmFPAfileNameFromRule (rule, file, view); + } + + // write out the IMAGE header segment + // this header block is new, write it to disk + if (hdu->header != file->header) { + // add EXTNAME, EXTHEAD, EXTTYPE to header + psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTABLE", PS_META_REPLACE, "name of table extension", tableName); + psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTRESID", PS_META_REPLACE, "name of resid extension", residName); + psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTYPE", PS_META_REPLACE, "extension type", "IMAGE"); + if (!file->phu) { + // this hdu->header acts as the PHU: set EXTEND to be true + psMetadataAddBool (hdu->header, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true); + } + + psFitsWriteBlank (file->fits, hdu->header, headName); + psTrace ("pmFPAfile", 5, "wrote ext head %s (type: %d)\n", file->filename, file->type); + file->header = hdu->header; + psFree (headName); + } // 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); + psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout"); + psFree (tableName); + psFree (residName); + return false; + } // 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); + { + // 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 (residName); + psFree (psfTable); + psFree (header); + return false; + } + } + psFree (tableName); + psFree (psfTable); + psFree (header); + } // 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); - psFree (images); - } else { - // this call creates an extension with NAXIS3 = 1 - psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, imageName); - } - - psFree (imageName); - psFree (header); + { + psMetadata *header = psMetadataAlloc (); + if (psf->residuals == NULL) { + // set some header keywords to make it clear there are no residuals? + psFitsWriteBlank (file->fits, header, residName); + psFree (residName); + 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, residName); + psFree (images); + } else { + // this call creates an extension with NAXIS3 = 1 + psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, residName); + } + psFree (residName); + psFree (header); + } + return true; @@ -458,4 +547,6 @@ psMetadata *header = NULL; + psTrace ("psModules.objects", 5, "read psf model for %s\n", file->filename); + // Menu of EXTNAME rules psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); @@ -493,4 +584,8 @@ char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME"); pmModelType type = pmModelSetType (modelName); + if (type == -1) { + psError(PS_ERR_UNKNOWN, true, "invalid model name %s in psf file %s", modelName, file->filename); + return false; + } bool poissonErrors = psMetadataLookupBool (&status, header, "POISSON"); @@ -509,6 +604,11 @@ snprintf (name, 9, "PAR%02d_NX", i); int nXorder = psMetadataLookupS32 (&status, header, name); + if (!status) continue; // not all parameters are defined snprintf (name, 9, "PAR%02d_NY", i); int nYorder = psMetadataLookupS32 (&status, header, name); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but no NY", i); + return false; + } psf->params_NEW->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder); } @@ -527,5 +627,5 @@ // fill in the matching psf->params entries - for (int i = 0; i > table->n; i++) { + for (int i = 0; i < table->n; i++) { psMetadata *row = table->data[i]; int iPar = psMetadataLookupS32 (&status, row, "MODEL_TERM"); @@ -535,4 +635,8 @@ psPolynomial2D *poly = psf->params_NEW->data[iPar]; + if (poly == NULL) { + psError(PS_ERR_UNKNOWN, true, "values for parameter %d, but missing NX", iPar); + return false; + } poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); @@ -583,4 +687,5 @@ psFree (tableName); psFree (imageName); + psFree (header); return true;