Index: branches/eam_branches/ipp-20130419/pswarp/src/pswarpMakePSF.c =================================================================== --- branches/eam_branches/ipp-20130419/pswarp/src/pswarpMakePSF.c (revision 35511) +++ branches/eam_branches/ipp-20130419/pswarp/src/pswarpMakePSF.c (revision 35511) @@ -0,0 +1,96 @@ +/** @file pswarpLoop.c + * + * @brief mail processing loop for pswarp + * @ingroup pswarp + * + * @author IfA + * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-13 21:54:32 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "pswarp.h" +#include +#include "pswarpFileNames.h" // Lists of file rules used at different stages + +#define WCS_NONLIN_TOL 0.001 // Non-linear tolerance for header WCS +#define TESTING 0 // Testing output? + +// We need a new PSF model for the warped frame. It would be good to generate this analytically, but +// that's going to be tricky. We have a list of sources, so we use those to redetermine the PSF model. +bool pswarpMakePSF (pmConfig *config) { + + // load the recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE); + if (!recipe) { + psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE); + return false; + } + + if (psMetadataLookupBool(&mdok, recipe, "PSF")) { + pswarpFileActivation(config, photFiles, true); + if (!pswarpIOChecksBefore(config)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto DONE; + } + + // supply the readout and fpa of interest to psphot + pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); + pmFPACopy(photFile->fpa, outFPA); + + pmFPAview *view = pmFPAviewAlloc(0); ///< View into skycell + view->chip = view->cell = view->readout = 0; + + // grab the sources of interest from the storage location (pmFPAfile PSPHOT.INPUT.CMF) + psArray *sources = psphotLoadPSFSources (config, view); + if (!sources) { + psError(psErrorCodeLast(), false, "No sources supplied to measure PSF"); + goto DONE; + } + + pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT); + + // measure the PSF using these sources + if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) { + // This is likely a data quality issue + // XXX Split into multiple cases using error codes? + psErrorStackPrint(stderr, "Unable to determine PSF"); + psWarning("Unable to determine PSF --- suspect bad data quality."); + if (stats && psMetadataLookupS32(NULL, stats, "QUALITY") == 0) { + psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, + "Unable to determine PSF", psErrorCodeLast()); + } + psErrorClear(); + psphotFilesActivate(config, false); + } + + // Ensure seeing is carried over + pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa); // Chip with seeing + psMetadataItem *item = psMetadataLookup(outChip->concepts, "CHIP.SEEING"); // Concept with seeing + item->data.F32 = psMetadataLookupF32(NULL, photChip->concepts, "CHIP.SEEING"); + +// XXX EAM : put this in a visualization function +#if (TESTING) + { + #define PSF_SIZE 20 ///< Half-size of PSF + #define PSF_FLUX 10000 ///< Central flux for PSF + pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa); + pmPSF *psf = psMetadataLookupPtr(NULL, photChip->analysis, "PSPHOT.PSF"); + psImage *image = psImageAlloc(2 * PSF_SIZE + 1, 2 * PSF_SIZE + 1, PS_TYPE_F32); + psImageInit(image, 0); + pmModel *model = pmModelFromPSFforXY(psf, PSF_SIZE, PSF_SIZE, PSF_FLUX); + pmModelAdd(image, NULL, model, PM_MODEL_OP_FULL, 0); + psFree(model); + psFits *fits = psFitsOpen("psf.fits", "w"); + psFitsWriteImage(fits, NULL, image, 0, NULL); + psFitsClose(fits); + psFree(image); + } +#endif + + psFree(view); + } + return true; +} +# endif + Index: branches/eam_branches/ipp-20130419/pswarp/src/pswarpUpdateStatistics.c =================================================================== --- branches/eam_branches/ipp-20130419/pswarp/src/pswarpUpdateStatistics.c (revision 35511) +++ branches/eam_branches/ipp-20130419/pswarp/src/pswarpUpdateStatistics.c (revision 35511) @@ -0,0 +1,179 @@ +/** @file pswarpLoop.c + * + * @brief mail processing loop for pswarp + * @ingroup pswarp + * + * @author IfA + * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-13 21:54:32 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "pswarp.h" +#include + +// once the output fpa elements have been built, loop over the fpa and generate stats +// for each readout +bool pswarpUpdateStatistics (pmFPA *output, psMetadata *fpaStats, pmFPA *input, pmFPA *astrom, pmConfig *config) { + + // load the recipe + bool status = false; + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE); + if (!recipe) { + psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE); + return false; + } + + // output mask bits + psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT"); + psAssert (status, "MASK.OUTPUT was not defined"); + + int Nreadout = 0; // count the number of readouts with data + + bool doMaskStats = psMetadataLookupBool(&status, recipe, "MASK.STATS"); + + pmFPAview *view = pmFPAviewAlloc(0); + + pmChip *chip; + while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) { + psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + + // use this in output metadata info (MD5 sums) + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); + + psMetadata *chipStats = psMetadataAlloc (); + psMetadataAdd (fpaStats, PS_LIST_TAIL, chipName, PS_DATA_METADATA, "chip stats folder", chipStats); + + pmCell *cell; + while ((cell = pmFPAviewNextCell (view, output, 1)) != NULL) { + psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + + // use this in output metadata info (MD5 sums) + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); + + psMetadata *cellStats = psMetadataAlloc (); + psMetadataAdd (chipStats, PS_LIST_TAIL, cellName, PS_DATA_METADATA, "cell stats folder", cellStats); + + // process each of the readouts + pmReadout *readout; + while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) { + // skip empty output readouts + if (!readout->data_exists) { + Nreadout ++; + continue; + } + + psString readoutName = NULL; + psStringAppend (&readoutName, "READOUT.%02d", view->readout); + + psMetadata *readoutStats = psMetadataAlloc (); + psMetadataAdd (cellStats, PS_LIST_TAIL, readoutName, PS_DATA_METADATA, "readout stats folder", readoutStats); + + if (!pswarpPixelsLit(readout, readoutStats, config)) { + psError(psErrorCodeLast(), false, "Unable to calculate pixel regions."); + psFree(view); + return false; + } + + if (doMaskStats && !pswarpMaskStats(readout, readoutStats, config)) { + psError(psErrorCodeLast(), false, "Unable to calculate mask stats."); + psFree(view); + return false; + } + + // Set covariance matrix for output + // XXX is this the right place for the covariance calculation? +# if (0) + psList *covariances = psMetadataLookupPtr(&mdok, output->analysis, PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices + psAssert(covariances, "Should be there"); + psArray *covars = psListToArray(covariances); // Array of covariance matrices + psKernel *covar = psImageCovarianceAverage(covars); + psFree(covars); + psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES); + + // Correct covariance matrix scale for the mean (square root of the) Jacobian + double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian + int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX); // Good pixels + jacobian /= goodPixels; + output->covariance = psImageCovarianceScale(covar, jacobian); + psFree(covar); + + if (output->variance) { + psImageCovarianceTransfer(output->variance, output->covariance); + } +# endif + + { + // Add MD5 information for readout + psString headerName = NULL; ///< Header name for MD5 + psVector *md5 = psImageMD5(readout->image); ///< md5 hash + psString md5string = psMD5toString(md5); ///< String + + pmHDU *hdu = pmHDUFromReadout(readout); + psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, "Image MD5", md5string); + psFree(md5); + psFree(md5string); + psFree(headerName); + } + psFree (readoutStats); + psFree (readoutName); + } + + psList *inputCells = psMetadataLookupPtr (NULL, cell->analysis, "INPUT.CELLS"); + if (!pmConceptsAverageCells(cell, inputCells, NULL, NULL, false)) { + psError(psErrorCodeLast(), false, "Unable to average cell concepts."); + psFree(view); + return false; + } + + // XXX Is this too ad-hoc? + psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); ///< Trim section + trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels + + psFree (cellStats); + } + psFree (chipStats); + } + + // Perform statistics on the output image (ppStatsFPA loops down the hierarchy) + if (fpaStats) { + pmFPAviewReset (view); + if (!ppStatsFPA(fpaStats, output, view, maskValue, config)) { + psWarning("Unable to perform statistics on warped image."); + } + } + psFree(view); + + if (!psMetadataCopy(output->concepts, input->concepts)) { + psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output."); + return false; + } + + if (Nreadout == 0) { + psWarning("No overlap between input and skycell."); + if (fpaStats) { + psMetadataAddS32(fpaStats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, "No overlap between input and skycell", PSWARP_ERR_NO_OVERLAP); + } + } + + // Update ZP from the astrometry + { + psMetadataItem *item = psMetadataLookup(output->concepts, "FPA.ZP"); + item->data.F32 = psMetadataLookupF32(NULL, astrom->concepts, "FPA.ZP"); + } + + // need to update the chip and fpa level astrometry appropriately, see + // code in psastro + +# if (0) + if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) { + psError(psErrorCodeLast(), false, "Unable to generate WCS header."); + return false; + } +# endif + + return true; +}