Index: trunk/ppImage/src/Makefile.am =================================================================== --- trunk/ppImage/src/Makefile.am (revision 33243) +++ trunk/ppImage/src/Makefile.am (revision 33590) @@ -59,4 +59,5 @@ ppImageParityFlip.c \ ppImageCheckCTE.c \ + ppImageCheckNoise.c \ ppImageFileCheck.c \ ppImageVersion.c \ Index: trunk/ppImage/src/ppImage.h =================================================================== --- trunk/ppImage/src/ppImage.h (revision 33243) +++ trunk/ppImage/src/ppImage.h (revision 33590) @@ -48,4 +48,5 @@ bool doStats; // call ppStats on the image bool checkCTE; // measure pixel-based variance + bool checkNoise; // measure cell-level variance bool applyParity; // Apply Cell parities @@ -172,4 +173,6 @@ bool ppImageCheckCTE(pmConfig *config, ppImageOptions *options, pmFPAview *view); +bool ppImageCheckNoise(pmConfig *config, ppImageOptions *options, pmFPAview *view); + bool ppImageBurntoolMask(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmReadout *mask); Index: trunk/ppImage/src/ppImageCheckNoise.c =================================================================== --- trunk/ppImage/src/ppImageCheckNoise.c (revision 33590) +++ trunk/ppImage/src/ppImageCheckNoise.c (revision 33590) @@ -0,0 +1,109 @@ +#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 ppImageCheckNoise(pmConfig *config, ppImageOptions *options, pmFPAview *view) +{ + bool status; + + // this step is complete optional. + if (!options->checkNoise) { + 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, "NOISE.XBIN"); + binning->nYbin = psMetadataLookupS32 (&status, recipe, "NOISE.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) continue; + if (fineRegion.y0 >= 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 = stats->robustStdev; + // float normVariance = PS_SQR(stats->sampleStdev) / cellMedian; + // if (!isfinite(normVariance)) fprintf (stderr, "** normVariance is nan **\n"); + + // apply resulting value to the input pixels + for (int jy = fineRegion.y0; jy < fineRegion.y1; jy++) { + if (jy < 0) continue; + if (jy >= image->numRows) continue; + for (int jx = fineRegion.x0; jx < fineRegion.x1; jx++) { + if (jx < 0) continue; + if (jx >= image->numCols) continue; + image->data.F32[jy][jx] = normVariance; + } + } + + psFree (subset); + psFree (submask); + } + } + + psFree (rng); + psFree (binning); + psFree (stats); + psFree (statsDefaults); + + return true; +} Index: trunk/ppImage/src/ppImageDetrendReadout.c =================================================================== --- trunk/ppImage/src/ppImageDetrendReadout.c (revision 33243) +++ trunk/ppImage/src/ppImageDetrendReadout.c (revision 33590) @@ -110,4 +110,6 @@ // offset information, we are not really getting exactly the right mapping from the // original file. + // CZW 2012-03-21: I do not believe this is true anymore. In any case, this seems to be what + // ppImageReplaceBackground does to do sky subtraction. psImageBinning *binning = psImageBinningAlloc(); binning->nXruff = noiseMap->image->numCols; @@ -119,4 +121,39 @@ psImageUnbin (noiseImage, noiseMap->image, binning); psFree (binning); + // Stolen from pmSkySubtract.c + // CZW: Unneeded, as psImageUnbin does the bilinear interpolation? It still looks blocky. +/* psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, RECIPE_NAME); // Recipe */ +/* psAssert(recipe, "Should be there!"); */ + +/* psS32 xBinFactor = psMetadataLookupS32(NULL,recipe,"NOISE.XBIN"); */ +/* psS32 yBinFactor = psMetadataLookupS32(NULL,recipe,"NOISE.YBIN"); */ +/* psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, */ +/* noiseMap->image, */ +/* NULL, NULL, */ +/* 0, 0.0, 0.0, 0, 0, 0.0); */ +/* for (psS32 row = 0; row < input->image->numRows; row++) { */ +/* for (psS32 col = 0; col < input->image->numCols; col++) { */ +/* // We calculate the F32 value of the pixel coordinates in the */ +/* // binned image and then use a pixel interpolation routine to */ +/* // determine the value of the pixel at that location. */ +/* psF32 binRowF64 = ((psF32) row) / ((psF32) yBinFactor); */ +/* psF32 binColF64 = ((psF32) col) / ((psF32) xBinFactor); */ +/* // We add 0.5 to the pixel locations since the pixel */ +/* // interpolation routine defines the location of pixel */ +/* // (i, j) as (i+0.5, j+0.5). */ +/* binRowF64+= 0.5; */ +/* binColF64+= 0.5; */ + +/* double binPixel; */ +/* if (!psImagePixelInterpolate(&binPixel, NULL, NULL, binColF64, binRowF64, interp)) { */ +/* psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); */ +/* psFree(interp); */ +/* psFree(noiseImage); */ +/* return NULL; */ +/* } */ +/* noiseImage->data.F32[row][col] = binPixel; */ +/* } */ +/* } */ +/* psFree(interp); */ } Index: trunk/ppImage/src/ppImageLoop.c =================================================================== --- trunk/ppImage/src/ppImageLoop.c (revision 33243) +++ trunk/ppImage/src/ppImageLoop.c (revision 33590) @@ -126,4 +126,8 @@ } + if (!ppImageCheckNoise (config, options, view)) { + ESCAPE("Unable to generate noisemap"); + } + // optionally degrade a MD image to 3pi exposure times if (!ppImageAddNoise(config, options, view, input->fpa)){ Index: trunk/ppImage/src/ppImageOptions.c =================================================================== --- trunk/ppImage/src/ppImageOptions.c (revision 33243) +++ trunk/ppImage/src/ppImageOptions.c (revision 33590) @@ -41,4 +41,5 @@ options->doStats = false; // Measure and save image statistics options->checkCTE = false; // Measure pixel-based variance + options->checkNoise = false; // Measure cell-level variances. options->applyParity = false; // Apply Cell parities options->doMaskStats = false; // Calculate mask fractions @@ -319,4 +320,5 @@ options->checkCTE = psMetadataLookupBool(NULL, recipe, "CHECK.CTE"); + options->checkNoise = psMetadataLookupBool(NULL, recipe, "CHECK.NOISE"); /* doMaskBuild : there are some cases where we require a mask, so we force doMaskBuild to be set even if the user specified 'FALSE'