Index: branches/eam_branches/ipp-20111122/ppImage/src/Makefile.am =================================================================== --- branches/eam_branches/ipp-20111122/ppImage/src/Makefile.am (revision 32855) +++ branches/eam_branches/ipp-20111122/ppImage/src/Makefile.am (revision 33638) @@ -59,4 +59,5 @@ ppImageParityFlip.c \ ppImageCheckCTE.c \ + ppImageCheckNoise.c \ ppImageFileCheck.c \ ppImageVersion.c \ Index: branches/eam_branches/ipp-20111122/ppImage/src/ppImage.h =================================================================== --- branches/eam_branches/ipp-20111122/ppImage/src/ppImage.h (revision 32855) +++ branches/eam_branches/ipp-20111122/ppImage/src/ppImage.h (revision 33638) @@ -40,4 +40,5 @@ bool doPatternRow; // Row pattern correction bool doPatternCell; // Cell pattern correction + bool doPatternContinuity; // Cell continuity correction bool doFringe; // Fringe subtraction bool doPhotom; // Source identification and photometry @@ -47,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 @@ -109,4 +111,6 @@ psStatsOptions patternCellMean; // Statistic for mean + int patternContinuityEdgeWidth; // Size of box to use for edge matching. + int remnanceSize; // Size for remnance detection float remnanceThresh; // Threshold for remnance detection @@ -169,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: branches/eam_branches/ipp-20111122/ppImage/src/ppImageCheckNoise.c =================================================================== --- branches/eam_branches/ipp-20111122/ppImage/src/ppImageCheckNoise.c (revision 33638) +++ branches/eam_branches/ipp-20111122/ppImage/src/ppImageCheckNoise.c (revision 33638) @@ -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: branches/eam_branches/ipp-20111122/ppImage/src/ppImageDetrendPattern.c =================================================================== --- branches/eam_branches/ipp-20111122/ppImage/src/ppImageDetrendPattern.c (revision 32855) +++ branches/eam_branches/ipp-20111122/ppImage/src/ppImageDetrendPattern.c (revision 33638) @@ -18,5 +18,5 @@ pmCell *cell = NULL; - assert(options->doPatternRow || options->doPatternCell); // do not call if not needed + assert(options->doPatternRow || options->doPatternCell || options->doPatternContinuity); // do not call if not needed assert(inputView->chip != -1); assert(inputView->cell == -1); @@ -40,5 +40,5 @@ if (psMetadataLookupBool(NULL,hdu->header,"PTRN_ROW")) { psLogMsg("ppImage", PS_LOG_INFO, "Not performing row pattern correction as it has already been done."); - goto pattern_cell; + goto pattern_continuity; } @@ -98,8 +98,54 @@ } + pattern_continuity: + + // see the comment for PATTERN.ROW; the same rules apply for PATTERN.CELL + + if (options->doPatternContinuity) { + int numCells = chip->cells->n; // Number of cells + psVector *tweak = psVectorAlloc(numCells, PS_TYPE_U8); // Tweak cell? + pmFPAview *view = pmFPAviewAlloc(0); // View for local processing + *view = *inputView; + + pmHDU *hdu = pmHDUFromChip(chip); + if (psMetadataLookupBool(NULL,hdu->header,"PTRN_CON")) { + psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction as it has already been done."); + goto pattern_cell; + } + + for (int i = 0; i < chip->cells->n; i++) { + view->cell = i; + + pmCell *cell = chip->cells->data[i]; // Cell of interest + + if (cell->readouts->n > 1) { + psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction on video cell."); + continue; + } + + bool doPattern = false; + if (!doPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CONTINUITY.SUBSET")) { + ESCAPE(false, "Unable to determine whether row pattern matching should be applied."); + } + if (doPattern) { + tweak->data.U8[i] = 0xFF; + } + } + + // Tweak the cells + if (!pmPatternContinuity(chip, tweak, options->patternCellBG, options->patternCellMean, + options->maskValue, options->darkMask,options->patternContinuityEdgeWidth)) { + psFree(tweak); + psFree(view); + return false; + } + psFree(tweak); + psFree(view); + + psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CON",PS_META_REPLACE,"PATTERN.CONTINUITY correction applied",true); + } + pattern_cell: - - // see the comment for PATTERN.ROW; the same rules apply for PATTERN.CELL - + if (options->doPatternCell) { int numCells = chip->cells->n; // Number of cells @@ -131,5 +177,4 @@ tweak->data.U8[i] = 0xFF; } - } Index: branches/eam_branches/ipp-20111122/ppImage/src/ppImageDetrendReadout.c =================================================================== --- branches/eam_branches/ipp-20111122/ppImage/src/ppImageDetrendReadout.c (revision 32855) +++ branches/eam_branches/ipp-20111122/ppImage/src/ppImageDetrendReadout.c (revision 33638) @@ -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: branches/eam_branches/ipp-20111122/ppImage/src/ppImageLoop.c =================================================================== --- branches/eam_branches/ipp-20111122/ppImage/src/ppImageLoop.c (revision 32855) +++ branches/eam_branches/ipp-20111122/ppImage/src/ppImageLoop.c (revision 33638) @@ -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)){ @@ -170,5 +174,5 @@ // Apply the pattern correction - if (options->doPatternRow || options->doPatternCell) { + if (options->doPatternRow || options->doPatternCell || options->doPatternContinuity) { if (!ppImageDetrendPatternApply(config,chip,view,options)) { ESCAPE("Unable to apply pattern corrections"); Index: branches/eam_branches/ipp-20111122/ppImage/src/ppImageOptions.c =================================================================== --- branches/eam_branches/ipp-20111122/ppImage/src/ppImageOptions.c (revision 32855) +++ branches/eam_branches/ipp-20111122/ppImage/src/ppImageOptions.c (revision 33638) @@ -34,4 +34,5 @@ options->doPatternRow = false; // Row pattern correction options->doPatternCell = false; // Cell pattern correction + options->doPatternContinuity = false; // Cell continuity correction options->doFringe = false; // Fringe subtraction options->doPhotom = false; // Source identification and photometry @@ -40,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 @@ -253,4 +255,5 @@ options->doPatternRow = psMetadataLookupBool(NULL, recipe, "PATTERN.ROW"); options->doPatternCell = psMetadataLookupBool(NULL, recipe, "PATTERN.CELL"); + options->doPatternContinuity = psMetadataLookupBool(NULL, recipe, "PATTERN.CONTINUITY"); options->doMaskStats = psMetadataLookupBool(NULL, recipe, "MASK.STATS"); @@ -317,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' @@ -408,4 +412,11 @@ } + if (psMetadataLookup(format, "PATTERN.CONTINUITY.WIDTH")) { + options->patternContinuityEdgeWidth = psMetadataLookupS32(NULL, format, "PATTERN.CONTINUITY.WIDTH"); + } + else { + options->patternContinuityEdgeWidth = psMetadataLookupS32(NULL, recipe, "PATTERN.CONTINUITY.WIDTH"); + } + // Remnance options