Index: trunk/ppImage/src/Makefile.am =================================================================== --- trunk/ppImage/src/Makefile.am (revision 23366) +++ trunk/ppImage/src/Makefile.am (revision 23411) @@ -38,4 +38,5 @@ ppImageDefineFile.c \ ppImageSetMaskBits.c \ + ppImageCheckCTE.c \ ppImageFileCheck.c \ ppImageVersion.c \ Index: trunk/ppImage/src/ppImage.h =================================================================== --- trunk/ppImage/src/ppImage.h (revision 23366) +++ trunk/ppImage/src/ppImage.h (revision 23411) @@ -40,4 +40,5 @@ bool doAstromMosaic; // full-mosaic Astrometry bool doStats; // call ppStats on the image + bool checkCTE; // measure pixel-based variance // output files requested @@ -133,4 +134,6 @@ bool ppImageFringeFree(pmConfig *config, pmFPAview *view); +bool ppImageCheckCTE(pmConfig *config, ppImageOptions *options, pmFPAview *view); + // Record which detrend file was used for the detrending bool ppImageDetrendRecord( Index: trunk/ppImage/src/ppImageCheckCTE.c =================================================================== --- trunk/ppImage/src/ppImageCheckCTE.c (revision 23411) +++ trunk/ppImage/src/ppImageCheckCTE.c (revision 23411) @@ -0,0 +1,106 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "ppImage.h" + +// XXX I originally coded this to create a new pmFPAfile, but in retrospect it makes more sense +// to treat this function as an operation on the input image + +// XXX make the choice of stats optional +bool ppImageCheckCTE(pmConfig *config, ppImageOptions *options, pmFPAview *view) +{ + bool status; + + // this step is complete optional. + if (!options->checkCTE) { + return true; + } + + // add recipe options supplied on command line + psMetadata *recipe = psMetadataLookupPtr(&status, config->recipes, RECIPE_NAME); + + // find the currently selected readout + pmReadout *inReadout = pmFPAfileThisReadout(config->files, view, "PPIMAGE.INPUT"); + + psImage *image = inReadout->image; + psImage *mask = inReadout->mask; + // psImage *variance = inReadout->variance; + + // I have the fine image size, I know the binning factor, determine the ruff image size + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = image->numCols; + binning->nYfine = image->numRows; + binning->nXbin = psMetadataLookupS32 (&status, recipe, "CTE.XBIN"); + binning->nYbin = psMetadataLookupS32 (&status, recipe, "CTE.YBIN"); + + psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkip(binning, image); + + psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN); + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psImageBackground (stats, NULL, image, mask, 0xffff, rng); + float cellMedian = stats->robustMedian; + psFree (stats); + + // stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV); + // psStats *statsDefaults = psStatsAlloc (PS_STAT_SAMPLE_STDEV); + + stats = psStatsAlloc (PS_STAT_ROBUST_STDEV); + psStats *statsDefaults = psStatsAlloc (PS_STAT_ROBUST_STDEV); + + // measure median and variance for subimages + psRegion ruffRegion = {0,0,0,0}; + psRegion fineRegion = {0,0,0,0}; + for (int iy = 0; iy < binning->nYruff; iy++) { + for (int ix = 0; ix < binning->nXruff; ix++) { + + // convert the ruff grid cell to the equivalent fine grid cell + ruffRegion = psRegionSet (ix, ix + 1, iy, iy + 1); + fineRegion = psImageBinningSetFineRegion (binning, ruffRegion); + fineRegion = psRegionForImage (image, fineRegion); + if (fineRegion.x0 >= image->numCols || fineRegion.x1 >= image->numCols || + fineRegion.y0 >= image->numRows || fineRegion.y1 >= image->numRows) { + continue; + } + + psImage *subset = psImageSubset (image, fineRegion); + if (!subset->numCols || !subset->numRows) { + psFree (subset); + continue; + } + psImage *submask = NULL; + if (mask) { + submask = psImageSubset (mask, fineRegion); + } + + // reset the default values + statsDefaults->tmpData = stats->tmpData; // XXX this is fairly hackish: tmpData is internal storage for stats; the assign drops the reference... + *stats = *statsDefaults; + statsDefaults->tmpData = NULL; + + psImageStats (stats, subset, submask, 0xffff); + + // XXX need to apply the gain as well + float normVariance = PS_SQR(stats->robustStdev) / cellMedian; + // float normVariance = PS_SQR(stats->sampleStdev) / cellMedian; + + // apply resulting value to the input pixels + for (int jy = fineRegion.y0; jy < fineRegion.y1; jy++) { + for (int jx = fineRegion.x0; jx < fineRegion.x1; jx++) { + image->data.F32[jy][jx] = normVariance; + } + } + + psFree (subset); + psFree (submask); + } + } + + psFree (rng); + psFree (binning); + psFree (stats); + psFree (statsDefaults); + + return true; +} Index: trunk/ppImage/src/ppImageCheckCTE.v1.c =================================================================== --- trunk/ppImage/src/ppImageCheckCTE.v1.c (revision 23411) +++ trunk/ppImage/src/ppImageCheckCTE.v1.c (revision 23411) @@ -0,0 +1,101 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "ppImage.h" + +bool ppImageCheckCTE(pmConfig *config, ppImageOptions *options, pmFPAview *view) +{ + bool status; + + // this step is complete optional. + if (!options->checkCTE) { + return true; + } + + // add recipe options supplied on command line + psMetadata *recipe = psMetadataLookupPtr(&status, config->recipes, RECIPE_NAME); + + // find the currently selected readout + pmReadout *inReadout = pmFPAfileThisReadout(config->files, view, "PPIMAGE.INPUT"); + + psImage *image = inReadout->image; + psImage *mask = inReadout->mask; + // psImage *variance = inReadout->variance; + + // I have the fine image size, I know the binning factor, determine the ruff image size + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = image->numCols; + binning->nYfine = image->numRows; + binning->nXbin = psMetadataLookupS32 (&status, recipe, "CTE.XBIN"); + binning->nYbin = psMetadataLookupS32 (&status, recipe, "CTE.YBIN"); + + psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkip(binning, image); + + pmCell *inCell = pmFPAfileThisCell (config->files, view, "PPIMAGE.INPUT"); + pmCell *outCell = pmFPAfileThisCell (config->files, view, "PPIMAGE.CTEMAP"); + if (!pmCellCopyStructure(outCell, inCell, binning->nXbin, binning->nYbin)) { + psError(PS_ERR_UNKNOWN, false, "Unable to copy cell structure."); + return false; + } + + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPIMAGE.CTEMAP"); + psImage *output = outRO->image; + + // Don't care about the bias: get rid of it, if present + psFree(outRO->bias); + outRO->bias = psListAlloc(NULL); + psMetadataItem *biassec = psMetadataLookup(outCell->concepts, "CELL.BIASSEC"); + psFree(biassec->data.V); + biassec->data.V = psListAlloc(NULL); + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + psStats *statsDefaults = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + + // measure median and variance for subimages + psRegion ruffRegion = {0,0,0,0}; + psRegion fineRegion = {0,0,0,0}; + for (int iy = 0; iy < output->numRows; iy++) { + for (int ix = 0; ix < output->numCols; ix++) { + + // convert the ruff grid cell to the equivalent fine grid cell + ruffRegion = psRegionSet (ix, ix + 1, iy, iy + 1); + fineRegion = psImageBinningSetFineRegion (binning, ruffRegion); + fineRegion = psRegionForImage (image, fineRegion); + if (fineRegion.x0 >= image->numCols || fineRegion.x1 >= image->numCols || + fineRegion.y0 >= image->numRows || fineRegion.y1 >= image->numRows) { + continue; + } + + psImage *subset = psImageSubset (image, fineRegion); + if (!subset->numCols || !subset->numRows) { + psFree (subset); + continue; + } + psImage *submask = NULL; + if (mask) { + submask = psImageSubset (mask, fineRegion); + } + + // reset the default values + statsDefaults->tmpData = stats->tmpData; // XXX this is fairly hackish: tmpData is internal storage for stats; the assign drops the reference... + *stats = *statsDefaults; + statsDefaults->tmpData = NULL; + + psImageStats (stats, subset, submask, 0); + + // XXX need to apply the gain as well + output->data.F32[iy][ix] = PS_SQR(stats->sampleStdev) / stats->sampleMedian; + + psFree (subset); + psFree (submask); + } + } + + psFree (binning); + psFree (stats); + psFree (statsDefaults); + + return true; +} Index: trunk/ppImage/src/ppImageLoop.c =================================================================== --- trunk/ppImage/src/ppImageLoop.c (revision 23366) +++ trunk/ppImage/src/ppImageLoop.c (revision 23411) @@ -94,4 +94,9 @@ if (!ppImageDetrendFree (config, view)) { ESCAPE("Unable to free detrend images"); + } + + // optionally measure CTE by examining the variance in a box + if (!ppImageCheckCTE (config, options, view)) { + ESCAPE("Unable to measure CTE"); } } Index: trunk/ppImage/src/ppImageOptions.c =================================================================== --- trunk/ppImage/src/ppImageOptions.c (revision 23366) +++ trunk/ppImage/src/ppImageOptions.c (revision 23411) @@ -33,4 +33,5 @@ options->doAstromMosaic = false; // Astrometry (full-mosaic) options->doStats = false; // Measure and save image statistics + options->checkCTE = false; // Measure pixel-based variance // output files requested @@ -258,4 +259,6 @@ options->doBG = psMetadataLookupBool(NULL, recipe, "BACKGROUND"); + options->checkCTE = psMetadataLookupBool(NULL, recipe, "CHECK.CTE"); + // even if not requested explicitly, if any of these are set, build an internal mask and variance: if (options->doBias || options->doOverscan || options->doDark || options->doShutter || options->doFlat || Index: trunk/ppImage/src/ppImageParseCamera.c =================================================================== --- trunk/ppImage/src/ppImageParseCamera.c (revision 23366) +++ trunk/ppImage/src/ppImageParseCamera.c (revision 23411) @@ -154,4 +154,15 @@ } + if (options->checkCTE && false) { + int DX = psMetadataLookupS32 (&status, recipe, "CTE.XBIN"); + int DY = psMetadataLookupS32 (&status, recipe, "CTE.YBIN"); + pmFPAfile *outCTE = pmFPAfileDefineFromFile (config, input, DX, DY, "PPIMAGE.CTEMAP"); + if (!outCTE) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PPIMAGE.CTEMAP"); + return false; + } + outCTE->save = true; + } + // the following files are output targets pmFPAfile *outImage = pmFPAfileDefineOutput(config, input->fpa, "PPIMAGE.OUTPUT");