Index: trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- trunk/psModules/src/objects/pmPSF_IO.c (revision 13212) +++ trunk/psModules/src/objects/pmPSF_IO.c (revision 13414) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-05-04 00:53:15 $ + * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-05-18 00:17:12 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -49,9 +49,9 @@ if (view->chip == -1) { - bool exists = pmFPACheckDataStatusForPSFmodel (fpa); + bool exists = pmFPACheckDataStatusForPSFmodel (fpa); return exists; } if (view->chip >= fpa->chips->n) { - psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n); + psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n); return false; } @@ -63,5 +63,5 @@ } if (view->cell >= chip->cells->n) { - psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n); + psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n); return false; } @@ -74,11 +74,10 @@ if (view->readout >= cell->readouts->n) { - psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouds->n == %ld", view->readout, cell->readouts->n); + psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouds->n == %ld", view->readout, cell->readouts->n); return false; } pmReadout *readout = cell->readouts->data[view->readout]; - bool exists = pmReadoutCheckDataStatusForPSFmodel (readout); - return exists; + return pmReadoutCheckDataStatusForPSFmodel (readout); } @@ -86,7 +85,7 @@ for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (!chip) continue; - if (pmChipCheckDataStatusForPSFmodel (chip)) return true; + pmChip *chip = fpa->chips->data[i]; + if (!chip) continue; + if (pmChipCheckDataStatusForPSFmodel (chip)) return true; } return false; @@ -96,7 +95,7 @@ for (int i = 0; i < chip->cells->n; i++) { - pmCell *cell = chip->cells->data[i]; - if (!cell) continue; - if (pmCellCheckDataStatusForPSFmodel (cell)) return true; + pmCell *cell = chip->cells->data[i]; + if (!cell) continue; + if (pmCellCheckDataStatusForPSFmodel (cell)) return true; } return false; @@ -106,7 +105,7 @@ for (int i = 0; i < cell->readouts->n; i++) { - pmReadout *readout = cell->readouts->data[i]; - if (!readout) continue; - if (pmReadoutCheckDataStatusForPSFmodel (readout)) return true; + pmReadout *readout = cell->readouts->data[i]; + if (!readout) continue; + if (pmReadoutCheckDataStatusForPSFmodel (readout)) return true; } return false; @@ -118,7 +117,6 @@ // select the psf of interest - pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); - if (!psf) return false; - return true; + pmPSF *psf = psMetadataLookupPtr(&status, readout->analysis, "PSPHOT.PSF"); + return psf ? true : false; } @@ -129,5 +127,5 @@ if (view->chip == -1) { - if (!pmFPAWritePSFmodel (fpa, view, file, config)) { + if (!pmFPAWritePSFmodel(fpa, view, file, config)) { psError(PS_ERR_IO, false, "Failed to write PSF for fpa"); return false; @@ -168,6 +166,6 @@ if (!pmReadoutWritePSFmodel (readout, view, file, config)) { - psError(PS_ERR_IO, false, "Failed to write PSF for readout"); - return false; + psError(PS_ERR_IO, false, "Failed to write PSF for readout"); + return false; } return true; @@ -179,5 +177,4 @@ for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; if (!pmChipWritePSFmodel (chip, view, file, config)) { @@ -192,7 +189,5 @@ 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)) { @@ -207,7 +202,5 @@ 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)) { @@ -219,5 +212,5 @@ } -// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, +// 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) @@ -228,6 +221,6 @@ pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); if (!psf) { - psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout"); - return false; + psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout"); + return false; } @@ -238,12 +231,12 @@ psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); if (!menu) { - psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); - return false; + 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; + psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config"); + return false; } char *tableName = pmFPAfileNameFromRule (rule, file, view); @@ -251,5 +244,5 @@ // write the PSF model parameters in a FITS table - // we need to write a header for the table, + // we need to write a header for the table, psMetadata *header = psMetadataAlloc(); @@ -264,10 +257,10 @@ // save the dimensions of each parameter for (int i = 0; i < nPar; i++) { - char name[9]; + char name[9]; psPolynomial2D *poly = psf->params_NEW->data[i]; if (poly == NULL) continue; - snprintf (name, 9, "PAR%02d_NX", i); + snprintf (name, 9, "PAR%02d_NX", i); psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX); - snprintf (name, 9, "PAR%02d_NY", i); + snprintf (name, 9, "PAR%02d_NY", i); psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY); } @@ -287,34 +280,34 @@ 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); - } - } - } - + 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?) + // 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; - } + 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); @@ -325,15 +318,15 @@ 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; + 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) + // write the residual images (3D) header = psMetadataAlloc (); if (psf->residuals == NULL) { - // set some header keywords to make it clear there are no residuals? + // set some header keywords to make it clear there are no residuals? psFitsWriteBlank (file->fits, header, imageName); - psFree (imageName); + psFree (imageName); psFree (header); return true; @@ -348,17 +341,17 @@ // 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); + // 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); - } - + // this call creates an extension with NAXIS3 = 1 + psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, imageName); + } + psFree (imageName); psFree (header); @@ -373,21 +366,20 @@ 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); - } + 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) { @@ -396,35 +388,31 @@ if (view->chip == -1) { - pmFPAReadPSFmodel (fpa, view, file, config); - return true; + return pmFPAReadPSFmodel(fpa, view, file, config); } if (view->chip >= fpa->chips->n) { - return false; + psAbort("Programming error: view does not apply to FPA."); } pmChip *chip = fpa->chips->data[view->chip]; if (view->cell == -1) { - pmChipReadPSFmodel (chip, view, file, config); - return true; + return pmChipReadPSFmodel(chip, view, file, config); } if (view->cell >= chip->cells->n) { - return false; + psAbort("Programming error: view does not apply to FPA."); } pmCell *cell = chip->cells->data[view->cell]; if (view->readout == -1) { - pmCellReadPSFmodel (cell, view, file, config); - return true; + return pmCellReadPSFmodel(cell, view, file, config); } if (view->readout >= cell->readouts->n) { - return false; + psAbort("Programming error: view does not apply to FPA."); } pmReadout *readout = cell->readouts->data[view->readout]; - pmReadoutReadPSFmodel (readout, view, file, config); - return true; + return pmReadoutReadPSFmodel(readout, view, file, config); } @@ -432,11 +420,10 @@ bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { - + bool success = true; // Was everything successful? for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - pmChipReadPSFmodel (chip, view, file, config); - } - return true; + success &= pmChipReadPSFmodel(chip, view, file, config); + } + return success; } @@ -444,11 +431,10 @@ bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { - + bool success = true; // Was everything successful? for (int i = 0; i < chip->cells->n; i++) { - pmCell *cell = chip->cells->data[i]; - pmCellReadPSFmodel (cell, view, file, config); - } - return true; + success &= pmCellReadPSFmodel (cell, view, file, config); + } + return success; } @@ -456,18 +442,17 @@ bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { - + bool success = true; // Was everything successful? 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, + success &= pmReadoutReadPSFmodel(readout, view, file, config); + } + return success; +} + +// 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; +bool pmReadoutReadPSFmodel(pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +{ + bool status; char *rule = NULL; psMetadata *header = NULL; @@ -476,12 +461,12 @@ psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); if (!menu) { - psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); - return false; + 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; + psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config"); + return false; } char *tableName = pmFPAfileNameFromRule (rule, file, view); @@ -489,6 +474,6 @@ 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; + psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config"); + return false; } char *imageName = pmFPAfileNameFromRule (rule, file, view); @@ -498,5 +483,5 @@ // 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); + psAbort("cannot find data extension %s in %s", tableName, file->filename); } @@ -521,10 +506,10 @@ // load the dimensions of each parameter for (int i = 0; i < nPar; i++) { - char name[9]; - snprintf (name, 9, "PAR%02d_NX", i); + char name[9]; + snprintf (name, 9, "PAR%02d_NX", i); int nXorder = psMetadataLookupS32 (&status, header, name); - snprintf (name, 9, "PAR%02d_NY", i); + snprintf (name, 9, "PAR%02d_NY", i); int nYorder = psMetadataLookupS32 (&status, header, name); - psf->params_NEW->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder); + psf->params_NEW->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder); } @@ -543,15 +528,15 @@ // 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 + 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"); + 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); @@ -562,5 +547,5 @@ // 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); + psAbort("cannot find data extension %s in %s", imageName, file->filename); } @@ -569,26 +554,26 @@ 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 - } + 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 + } }