Index: trunk/ppImage/src/Makefile.am =================================================================== --- trunk/ppImage/src/Makefile.am (revision 13924) +++ trunk/ppImage/src/Makefile.am (revision 13970) @@ -26,4 +26,5 @@ ppImageAstrom.c \ ppImageAddstar.c \ + ppImageStats.c \ ppImageDefineFile.c \ ppImageVersion.c @@ -48,4 +49,5 @@ ppImageAstrom.c \ ppImageAddstar.c \ + ppImageStats.c \ ppImageDefineFile.c \ ppImageVersion.c Index: trunk/ppImage/src/ppImage.h =================================================================== --- trunk/ppImage/src/ppImage.h (revision 13924) +++ trunk/ppImage/src/ppImage.h (revision 13970) @@ -65,3 +65,9 @@ bool ppImageDefineFile (pmConfig *config, pmFPA *input, char *filerule, char *argname, pmFPAfileType fileType, pmDetrendType detrendType); +// calculate stats, including MD5 +bool ppImageStats (pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options); + +// write stats to output file +bool ppImageStatsOutput (pmConfig *config, pmFPA *fpa, const ppImageOptions *options); + #endif Index: trunk/ppImage/src/ppImageDetrendFringe.c =================================================================== --- trunk/ppImage/src/ppImageDetrendFringe.c (revision 13924) +++ trunk/ppImage/src/ppImageDetrendFringe.c (revision 13970) @@ -10,5 +10,5 @@ #include "ppImageDetrendFringe.h" -bool ppImageDetrendFringeMeasure(pmReadout *readout, pmCell *fringe, const ppImageOptions *options) +bool ppImageDetrendFringeMeasure(pmReadout *readout, pmCell *fringe, const bool isResidual, const ppImageOptions *options) { PS_ASSERT_PTR_NON_NULL(readout, false); @@ -49,6 +49,13 @@ psBinaryOp(measurements->df, measurements->df, "*", psScalarAlloc(expTime, PS_TYPE_F32)); + char *scienceFringes = NULL; + if (isResidual) { + scienceFringes = psStringCopy ("FRINGE.RESIDUALS"); + } else { + scienceFringes = psStringCopy ("FRINGE.MEASUREMENTS"); + } + // Science fringe measurements - pmFringeStats *previous = psMetadataLookupPtr(NULL, readout->parent->analysis, "FRINGE.MEASUREMENTS"); + pmFringeStats *previous = psMetadataLookupPtr(NULL, readout->parent->analysis, scienceFringes); if (previous) { // Multiple readouts: concatenate @@ -73,8 +80,8 @@ } - psMetadataAdd(readout->parent->analysis, PS_LIST_TAIL, "FRINGE.MEASUREMENTS", + psMetadataAdd(readout->parent->analysis, PS_LIST_TAIL, scienceFringes, PS_DATA_UNKNOWN | PS_META_REPLACE, "Fringe measurements", measurements); psFree(measurements); - + psFree(scienceFringes); psFree(references); @@ -84,6 +91,6 @@ // Pull the fringes out of the cell analysis FRINGE.MEASUREMENTS for a chip -static psArray *getFringes(const pmChip *chip // Chip of interest - ) +// XXX need some error checks +static psArray *getFringes(const pmChip *chip, const char *source) { psArray *cells = chip->cells; // Component cells @@ -91,6 +98,5 @@ for (int i = 0; i < cells->n; i++) { pmCell *cell = cells->data[i]; // Cell of interest - fringes->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, cell->analysis, - "FRINGE.MEASUREMENTS")); + fringes->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, cell->analysis, source)); } @@ -101,5 +107,5 @@ // Solve the fringe system: we have science fringe measurements for each cell, and an array of reference // fringe measurements for each cell. Need to concatenate these together first, and then solve. -bool ppImageDetrendFringeSolve(pmChip *scienceChip, const pmChip *refChip, const ppImageOptions *options) +bool ppImageDetrendFringeSolve(pmChip *scienceChip, const pmChip *refChip, const bool isResidual, const ppImageOptions *options) { PS_ASSERT_PTR_NON_NULL(scienceChip, NULL); @@ -107,5 +113,11 @@ PS_ASSERT_PTR_NON_NULL(options, NULL); - psArray *science = getFringes(scienceChip); // Fringe measurements on science chip + psArray *science = NULL; + if (isResidual) { + science = getFringes(scienceChip, "FRINGE.RESIDUALS"); // Fringe residuals on science chip + } else { + science = getFringes(scienceChip, "FRINGE.MEASUREMENTS"); // Fringe measurements on science chip + } + pmFringeStats *scienceCat = pmFringeStatsConcatenate(science, NULL, NULL); // Science fringes psFree(science); @@ -113,5 +125,5 @@ // Need to transform the array of cells each with an array of fringes --> array of fringes for the chip as // a whole - psArray *references = getFringes(refChip); // Fringe measurements on reference chip + psArray *references = getFringes(refChip, "FRINGE.MEASUREMENTS"); // Fringe measurements on reference chip int numRefs = ((psArray*)references->data[0])->n; // Number of reference fringes psArray *referencesCat = psArrayAlloc(numRefs); // Reference fringes @@ -132,6 +144,36 @@ options->fringeIter, options->fringeKeep); - psMetadataAdd(scienceChip->analysis, PS_LIST_TAIL, "FRINGE.SOLUTION", PS_DATA_UNKNOWN, - "Fringe solution", solution); + if (isResidual) { + psMetadataAdd(scienceChip->analysis, PS_LIST_TAIL, "FRINGE.RESIDUAL.SOLUTION", PS_DATA_UNKNOWN, "Fringe solution", solution); + } else { + psMetadataAdd(scienceChip->analysis, PS_LIST_TAIL, "FRINGE.SOLUTION", PS_DATA_UNKNOWN, "Fringe solution", solution); + } + +# if (0) + // write the fringe amplitude or residual amplitude to the header + // XXX this is measured per cell, but we only have headers per chip + pmHDU *hdu = pmHDUFromCell(science);// HDU of interest + for (int i = 0; i < solution->nFringeFrames; i++) { + // write metadata header value + psString keyword = NULL; + if (isResidual) { + psStringAppend (&keyword, "FRES_%02dV", i); + } else { + psStringAppend (&keyword, "FRNG_%02dV", i); + } + psMetadataAddF32(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "Fringe Amplitude", solution->coeff->data.F32[i + 1]); + psFree (keyword); + + keyword = NULL; + if (isResidual) { + psStringAppend (&keyword, "FRES_%02dE", i); + } else { + psStringAppend (&keyword, "FRNG_%02dE", i); + } + psMetadataAddF32(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "Fringe Amplitude error", solution->coeffErr->data.F32[i + 1]); + psFree (keyword); + } +# endif + psFree(solution); psFree(scienceCat); @@ -142,5 +184,5 @@ -bool ppImageDetrendFringeGenerate(pmCell *science, pmCell *fringes) +bool ppImageDetrendFringeGenerate(pmCell *science, pmCell *fringes, const ppImageOptions *options) { PS_ASSERT_PTR_NON_NULL(science, false); @@ -191,14 +233,18 @@ psFree(md5string); +# if (0) // write metadata header value + // XXX this is measured per cell, but we only have headers per chip psString keyword = NULL; - psStringAppend (&keyword, "FRNG_%02d", i); + psStringAppend (&keyword, "FRNG_%02dV", i); psMetadataAddF32(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "Fringe Amplitude", solution->coeff->data.F32[i + 1]); psFree (keyword); keyword = NULL; - psStringAppend (&keyword, "FRNG_%02dD", i); + psStringAppend (&keyword, "FRNG_%02dE", i); psMetadataAddF32(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "Fringe Amplitude error", solution->coeffErr->data.F32[i + 1]); psFree (keyword); +# endif + } if (expTime != 1.0) { @@ -219,4 +265,7 @@ return false; } + + // measure residual fringe amplitude. results go to FRINGE.RESIDUALS + ppImageDetrendFringeMeasure (readout, fringes, true, options); } psFree(sumFringe); @@ -233,2 +282,62 @@ } +bool ppImageDetrendFringeApply (pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options) { + + pmCell *cell = NULL; + + assert (options->doFringe); // do not call if not needed + assert (inputView->chip != -1); + assert (inputView->cell == -1); + assert (inputView->readout == -1); + + pmFPAview *view = pmFPAviewAlloc(0); // View for local processing + *view = *inputView; + + // select the reference chip + pmChip *fringe = pmFPAfileThisChip(config->files, view, "PPIMAGE.FRINGE"); + if (!fringe) { + psError(PS_ERR_UNKNOWN, false, "missing fringe reference data.\n"); + psFree (view); + return false; + } + + // Solve the fringe system + if (!ppImageDetrendFringeSolve(chip, fringe, false, options)) { + psError(PS_ERR_UNKNOWN, false, "failed to solve the fringe system.\n"); + psFree (view); + return false; + } + + // Go back over the cells to apply the fringe correction + view->cell = view->readout = -1; + while ((cell = pmFPAviewNextCell(view, chip->parent, 1)) != NULL) { + if (!cell->process || !cell->file_exists) { + continue; + } + + // Apply the fringe correction + psTrace("ppImage", 3, "Applying fringe correction...\n"); + pmCell *fringeCell = pmFPAfileThisCell(config->files, view, "PPIMAGE.FRINGE"); + if (!fringeCell) { + psError(PS_ERR_UNKNOWN, false, "missing fringe reference data.\n"); + psFree (view); + return false; + } + + if (!ppImageDetrendFringeGenerate(cell, fringeCell, options)) { + psError(PS_ERR_UNKNOWN, false, "failed to apply fringe image.\n"); + psFree (view); + return false; + } + } + + // Solve the residual fringe system + if (!ppImageDetrendFringeSolve(chip, fringe, true, options)) { + psError(PS_ERR_UNKNOWN, false, "failed to solve the residual fringe system.\n"); + psFree (view); + return false; + } + + psFree (view); + return true; +} Index: trunk/ppImage/src/ppImageDetrendFringe.h =================================================================== --- trunk/ppImage/src/ppImageDetrendFringe.h (revision 13924) +++ trunk/ppImage/src/ppImageDetrendFringe.h (revision 13970) @@ -9,4 +9,5 @@ bool ppImageDetrendFringeMeasure(pmReadout *readout, // Readout to measure pmCell *fringe, // Fringe cell (each readout is a different component) + const bool isResidual, const ppImageOptions *options // Options ); @@ -15,4 +16,5 @@ bool ppImageDetrendFringeSolve(pmChip *scienceChip, // Chip with science const pmChip *refChip, // Chip with reference fringes + const bool isResidual, const ppImageOptions *options // Options ); @@ -20,7 +22,14 @@ // Generate fringe frame bool ppImageDetrendFringeGenerate(pmCell *science, // Science cell - pmCell *fringes // Fringe cell, one readout per fringe component + pmCell *fringes, // Fringe cell, one readout per fringe component + const ppImageOptions *options // Options ); +bool ppImageDetrendFringeApply (pmConfig *config, // config + pmChip *chip, // science chip + const pmFPAview *inputView, // current view + const ppImageOptions *options // options +); + #endif Index: trunk/ppImage/src/ppImageDetrendReadout.c =================================================================== --- trunk/ppImage/src/ppImageDetrendReadout.c (revision 13924) +++ trunk/ppImage/src/ppImageDetrendReadout.c (revision 13970) @@ -88,5 +88,5 @@ if (options->doFringe) { pmCell *fringe = pmFPAfileThisCell(config->files, detview, "PPIMAGE.FRINGE"); - if (!ppImageDetrendFringeMeasure(input, fringe, options)) { + if (!ppImageDetrendFringeMeasure(input, fringe, false, options)) { return false; } Index: trunk/ppImage/src/ppImageLoop.c =================================================================== --- trunk/ppImage/src/ppImageLoop.c (revision 13924) +++ trunk/ppImage/src/ppImageLoop.c (revision 13970) @@ -3,5 +3,4 @@ #endif -#include #include "ppImage.h" #include "ppImageDetrendFringe.h" @@ -10,24 +9,8 @@ bool ppImageLoop (pmConfig *config, ppImageOptions *options) { - bool mdok; // Status of MD lookup bool status; pmChip *chip; pmCell *cell; pmReadout *readout; - psMetadata *stats = NULL; - - const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); // Filename for statistics - FILE *statsFile = NULL; // File stream for statistics - if (mdok && statsName && strlen(statsName) > 0) { - psString resolved = pmConfigConvertFilename(statsName, config, true); // Resolved filename - statsFile = fopen(resolved, "w"); - if (!statsFile) { - psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); - psFree(resolved); - return false; - } - stats = psMetadataAlloc(); - psFree(resolved); - } pmFPAfile *input = psMetadataLookupPtr(&status, config->files, "PPIMAGE.INPUT"); @@ -96,75 +79,18 @@ } - // Solve the fringe system + // Apply the fringe correction if (options->doFringe) { - pmChip *fringe = pmFPAfileThisChip(config->files, view, "PPIMAGE.FRINGE"); - if (!ppImageDetrendFringeSolve(chip, fringe, options)) { + if (!ppImageDetrendFringeApply (config, chip, view, options)) { psFree (view); return false; } - } - - // Go back over the cells to apply the fringe and do statistics - view->cell = view->readout = -1; - while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) { - if (!cell->process || !cell->file_exists) { - continue; - } - - // Apply the fringe correction - if (options->doFringe) { - psTrace("ppImage", 3, "Applying fringe correction...\n"); - pmCell *fringeCell = pmFPAfileThisCell(config->files, view, "PPIMAGE.FRINGE"); - if (!ppImageDetrendFringeGenerate(cell, fringeCell)) { - psFree (view); - return false; - } - } - - // Perform statistics on the detrended cell - if (stats) { - bool mdok; // Status of MD lookup - pmFPAfile *output = psMetadataLookupPtr(&mdok, config->files, "PPIMAGE.OUTPUT"); - if (!mdok || !output) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PPIMAGE.OUTPUT."); - psFree (view); - return false; - } - - if (!ppStats(stats, output->fpa, view, - options->satMask | options->badMask | options->maskValue, - config)) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate stats for image."); - psFree(stats); - psFree(view); - return false; - } - - if (options->doFringe && !ppStatsFringe(stats, chip, "FRINGE.SOLUTION")) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to extract fringe solution for image."); - psFree(stats); - psFree(view); - return false; - } - } - - // Add MD5 information for cell - pmHDU *hdu = pmHDUFromCell(cell); // HDU that owns the cell - while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) { - const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); - const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); - - psString headerName = NULL; // Header name for MD5 - psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout); - - psVector *md5 = psImageMD5(readout->image); // md5 hash - psString md5string = psMD5toString(md5); // String - psFree(md5); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, - "Image MD5", md5string); - psFree(md5string); - psFree(headerName); - } - } + } + + // measure various statistics for this image + if (!ppImageStats (config, chip, view, options)) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to measures stats for image"); + psFree (view); + return false; + } if (!ppImageMosaicChip(config, options, view, "PPIMAGE.CHIP", "PPIMAGE.OUTPUT")) { @@ -214,14 +140,8 @@ // Write out summary statistics - if (stats) { - char *statsMDC = psMetadataConfigFormat(stats); - if (!statsMDC || strlen(statsMDC) == 0) { - psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n"); - } else { - fprintf(statsFile, "%s", statsMDC); - } - psFree(statsMDC); - fclose(statsFile); - psFree(stats); + if (!ppImageStatsOutput (config, input->fpa, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write statistics file.\n"); + psFree (view); + return false; } Index: trunk/ppImage/src/ppImageOptions.c =================================================================== --- trunk/ppImage/src/ppImageOptions.c (revision 13924) +++ trunk/ppImage/src/ppImageOptions.c (revision 13970) @@ -39,4 +39,6 @@ options->doAstromChip = false; // Astrometry options->doAstromMosaic = false; // Astrometry + + options->doStats = false; // Measure and save image statistics options->BaseFITS = false; // create binned image (scale 1) @@ -187,4 +189,10 @@ options->doShutter = psMetadataLookupBool(NULL, recipe, "SHUTTER"); + options->doStats = false; + char *statsName = psMetadataLookupStr(&status, config->arguments, "STATS"); // Filename for statistics + if (statsName) { + options->doStats = true; + } + // binned image options options->xBin1 = psMetadataLookupS32(&status, recipe, "BIN1.XBIN"); Index: trunk/ppImage/src/ppImageOptions.h =================================================================== --- trunk/ppImage/src/ppImageOptions.h (revision 13924) +++ trunk/ppImage/src/ppImageOptions.h (revision 13970) @@ -29,4 +29,6 @@ void *nonLinearSource; + bool doStats; + bool BaseFITS; bool ChipFITS; Index: trunk/ppImage/src/ppImageStats.c =================================================================== --- trunk/ppImage/src/ppImageStats.c (revision 13970) +++ trunk/ppImage/src/ppImageStats.c (revision 13970) @@ -0,0 +1,147 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppImage.h" +#include "ppStats.h" + +// calculate stats, including MD5 +bool ppImageStats (pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options) { + + bool mdok; // Status of MD lookup + pmCell *cell = NULL; + pmReadout *readout = NULL; + + pmFPAview *view = pmFPAviewAlloc(0); // View for local processing + *view = *inputView; + + // select or create the fpa-level analysis stats metadata: + psMetadata *stats = psMetadataLookupPtr (&mdok, chip->parent->analysis, "PPIMAGE.STATS"); + psMemIncrRefCounter (stats); + if (!stats) { + stats = psMetadataAlloc (); + if (!psMetadataAdd (chip->parent->analysis, PS_LIST_TAIL, "PPIMAGE.STATS", PS_DATA_METADATA, "stats container", stats)) { + psError(PS_ERR_UNKNOWN, false, "Unable to push stats on fpa.analysis."); + psFree (view); + psFree (stats); + return false; + } + } + + view->cell = view->readout = -1; + while ((cell = pmFPAviewNextCell(view, chip->parent, 1)) != NULL) { + if (!cell->process || !cell->file_exists) { + continue; + } + + // Perform statistics on the detrended cell + if (options->doStats) { + + pmFPAfile *output = psMetadataLookupPtr(&mdok, config->files, "PPIMAGE.OUTPUT"); + if (!mdok || !output) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PPIMAGE.OUTPUT."); + psFree (view); + psFree(stats); + return false; + } + + if (!ppStats(stats, output->fpa, view, + options->satMask | options->badMask | options->maskValue, + config)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate stats for image."); + psFree(view); + psFree(stats); + return false; + } + + // extract the fringe amplitude from the analysis + if (options->doFringe && !ppStatsFringe(stats, chip, "FRINGE", "FRINGE.SOLUTION")) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to extract fringe solution for image."); + psFree(view); + psFree(stats); + return false; + } + + // extract the fringe residual amplitude from the analysis + if (options->doFringe && !ppStatsFringe(stats, chip, "FRINGE_RESID", "FRINGE.RESIDUAL.SOLUTION")) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to extract fringe solution for image."); + psFree(view); + psFree(stats); + return false; + } + } + + // Add MD5 information for cell + pmHDU *hdu = pmHDUFromCell(cell); // HDU that owns the cell + while ((readout = pmFPAviewNextReadout(view, chip->parent, 1)) != NULL) { + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); + + psString headerName = NULL; // Header name for MD5 + psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout); + + psVector *md5 = psImageMD5(readout->image); // md5 hash + psString md5string = psMD5toString(md5); // String + psFree(md5); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, + "Image MD5", md5string); + psFree(md5string); + psFree(headerName); + } + } + psFree(view); + psFree(stats); + return true; +} + +// write stats to output file +bool ppImageStatsOutput (pmConfig *config, pmFPA *fpa, const ppImageOptions *options) { + + bool mdok; + + // measure statistics, or ignore? + if (!options->doStats) return true; + + // get the output stats filename + const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); // Filename for statistics + if (!statsName && !strlen(statsName)) { + psError (PS_ERR_UNEXPECTED_NULL, false, "missing STATS entry in arguments list."); + return false; + } + + // convert to a real UNIX filename + psString resolved = pmConfigConvertFilename(statsName, config, true); // Resolved filename + FILE *statsFile = fopen (resolved, "w"); + if (!statsFile) { + psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); + psFree(resolved); + return false; + } + psFree(resolved); + + // select the fpa-level analysis stats metadata: + psMetadata *stats = psMetadataLookupPtr (&mdok, fpa->analysis, "PPIMAGE.STATS"); + if (!stats) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PPIMAGE.STATS entry in fpa analysis.\n"); + return false; + } + + // convert the stats MDC to a string + char *statsMDC = psMetadataConfigFormat(stats); + if (!statsMDC || strlen(statsMDC) == 0) { + psError(PS_ERR_IO, false, "Unable to serialize stats metadata.\n"); + return false; + } + + // write the stats MDC to a file + // XXX why does this not call psMetadataConfigPrint? + fprintf(statsFile, "%s", statsMDC); + + psFree(statsMDC); + fclose(statsFile); + return true; +}