Index: branches/pap/ppImage/notes.txt =================================================================== --- branches/pap/ppImage/notes.txt (revision 23948) +++ branches/pap/ppImage/notes.txt (revision 25027) @@ -1,23 +1,15 @@ -ppImage pmFPAfiles: +20090611 : adjusting things to allow the readnoise to be an image map for each cell - pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PPIMAGE.INPUT", "INPUT"); - pmFPAfile *inputMask = pmFPAfileBindFromArgs(&status, input, config, "PPIMAGE.INPUT.MASK", "INPUT.MASK"); - pmFPAfile *inputWeight = pmFPAfileBindFromArgs(&status, input, config, "PPIMAGE.INPUT.WEIGHT", "INPUT.WEIGHT"); - pmFPAfile *output = pmFPAfileDefineOutput(config, input->fpa, "PPIMAGE.OUTPUT"); - pmFPAfile *outMask = pmFPAfileDefineOutput(config, input->fpa, "PPIMAGE.OUTPUT.MASK"); - pmFPAfile *outWeight = pmFPAfileDefineOutput(config, input->fpa, "PPIMAGE.OUTPUT.WEIGHT"); - pmFPAfile *chipImage = pmFPAfileDefineChipMosaic(config, input->fpa, "PPIMAGE.CHIP"); - pmFPAfile *chipMask = pmFPAfileDefineOutput(config, chipImage->fpa, "PPIMAGE.CHIP.MASK"); - pmFPAfile *chipWeight = pmFPAfileDefineOutput(config, chipImage->fpa, "PPIMAGE.CHIP.WEIGHT"); - pmFPAfile *byFPA1 = pmFPAfileDefineFPAMosaic(config, input->fpa, "PPIMAGE.OUTPUT.FPA1"); - pmFPAfile *byFPA2 = pmFPAfileDefineFPAMosaic(config, input->fpa, "PPIMAGE.OUTPUT.FPA2"); + * pmReadoutSetVaraince : when variance is created, we are suppling a single CELL.READNOISE value -- change this to an optional image (carried through concepts?) - pmFPAfile *psphotInput = pmFPAfileDefineFromFPA (config, chipImage->fpa, 1, 1, "PSPHOT.INPUT"); - pmFPAfile *psphotOutput = psMetadataLookupPtr (&status, config->files, "PSPHOT.OUTPUT"); - pmFPAfile *psastroInput = pmFPAfileDefineInput (config, psphotOutput->fpa, "PSASTRO.INPUT"); + I need to: - pmFPAfile *bin1 = pmFPAfileDefineFromFPA (config, chipImage->fpa, options->xBin1, options->yBin1, "PPIMAGE.BIN1"); - pmFPAfile *bin2 = pmFPAfileDefineFromFPA (config, chipImage->fpa, options->xBin2, options->yBin2, "PPIMAGE.BIN2"); - pmFPAfile *jpg1 = pmFPAfileDefineOutput (config, byFPA1->fpa, "PPIMAGE.JPEG1"); - pmFPAfile *jpg2 = pmFPAfileDefineOutput (config, byFPA2->fpa, "PPIMAGE.JPEG2"); + * define a format for the readnoise map (image map in SPLIT/MEF format?) + * load the readnoise map + + * supply the map for each cell to the pmReadoutSetVariance (pointer to ppImageMap?) + pmReadoutSetVariance -> add psImageMap to API (also pmReadoutGenerateVariance) + + * add recipe information to define the source of the READNOISE info. + Index: branches/pap/ppImage/src/Makefile.am =================================================================== --- branches/pap/ppImage/src/Makefile.am (revision 23948) +++ branches/pap/ppImage/src/Makefile.am (revision 25027) @@ -54,4 +54,5 @@ ppImageDefineFile.c \ ppImageSetMaskBits.c \ + ppImageParityFlip.c \ ppImageCheckCTE.c \ ppImageFileCheck.c \ Index: branches/pap/ppImage/src/ppImage.h =================================================================== --- branches/pap/ppImage/src/ppImage.h (revision 23948) +++ branches/pap/ppImage/src/ppImage.h (revision 25027) @@ -31,4 +31,5 @@ bool doNonLin; // Non-linearity correction bool doOverscan; // Overscan subtraction + bool doNoiseMap; // Bias subtraction bool doBias; // Bias subtraction bool doDark; // Dark subtraction @@ -36,4 +37,5 @@ bool doShutter; // Shutter correction bool doFlat; // Flat-field normalisation + bool doPattern; // Pattern noise subtraction bool doFringe; // Fringe subtraction bool doPhotom; // Source identification and photometry @@ -43,7 +45,8 @@ bool doStats; // call ppStats on the image bool checkCTE; // measure pixel-based variance - - bool doCrosstalkMeasure; // measure crosstalk signal - bool doCrosstalkCorrect; // apply crosstalk correction + bool applyParity; // Apply Cell parities + + bool doCrosstalkMeasure; // measure crosstalk signal + bool doCrosstalkCorrect; // apply crosstalk correction // output files requested @@ -69,4 +72,5 @@ psImageMaskType lowMask; // Mask value to give bad pixels psImageMaskType flatMask; // Mask value to give bad flat pixels + psImageMaskType darkMask; // Mask value to give bad dark pixels psImageMaskType blankMask; // Mask value to give blank pixels @@ -89,4 +93,12 @@ int fringeIter; // Fringe iterations float fringeKeep; // Fringe keep fraction + + // Pattern noise subtraction + int patternOrder; // Polynomial order + int patternIter; // Clipping iterations + float patternRej; // Clipping threshold + float patternThresh; // Ignore threshold + psStatsOptions patternMean; // Statistic for mean + psStatsOptions patternStdev; // Statistic for stdev int remnanceSize; // Size for remnance detection @@ -117,4 +129,7 @@ bool ppImageSetMaskBits (pmConfig *config, ppImageOptions *options); +// apply the cell flips to the input data before analysis +bool ppImageParityFlip (pmConfig *config, const ppImageOptions *options, const pmFPAview *view); + // Loop over the input bool ppImageLoop(pmConfig *config, ppImageOptions *options); Index: branches/pap/ppImage/src/ppImageArguments.c =================================================================== --- branches/pap/ppImage/src/ppImageArguments.c (revision 23948) +++ branches/pap/ppImage/src/ppImageArguments.c (revision 25027) @@ -17,4 +17,5 @@ fprintf(stderr, "\n"); fprintf(stderr, "Input options (single file / file list):\n"); + fprintf(stderr, "\t-noisemap/-noisemaplist: Noise Map image.\n"); fprintf(stderr, "\t-bias/-biaslist: Bias image.\n"); fprintf(stderr, "\t-dark/-darklist: Dark image.\n"); @@ -112,4 +113,5 @@ // if these command-line options are supplied, load the file name lists into config->arguments // override any configuration-specified source for these files + pmConfigFileSetsMD (config->arguments, &argc, argv, "NOISEMAP", "-noisemap", "-noisemaplist"); pmConfigFileSetsMD (config->arguments, &argc, argv, "BIAS", "-bias", "-biaslist"); pmConfigFileSetsMD (config->arguments, &argc, argv, "DARK", "-dark", "-darklist"); Index: branches/pap/ppImage/src/ppImageCheckCTE.c =================================================================== --- branches/pap/ppImage/src/ppImageCheckCTE.c (revision 23948) +++ branches/pap/ppImage/src/ppImageCheckCTE.c (revision 25027) @@ -60,8 +60,6 @@ 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; - } + if (fineRegion.x0 >= image->numCols) continue; + if (fineRegion.y0 >= image->numRows) continue; psImage *subset = psImageSubset (image, fineRegion); @@ -85,8 +83,13 @@ float normVariance = PS_SQR(stats->robustStdev) / cellMedian; // 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; } Index: branches/pap/ppImage/src/ppImageDetrendFree.c =================================================================== --- branches/pap/ppImage/src/ppImageDetrendFree.c (revision 23948) +++ branches/pap/ppImage/src/ppImageDetrendFree.c (revision 25027) @@ -8,4 +8,5 @@ static char *detrendTypes[] = { "PPIMAGE.MASK", + "PPIMAGE.NOISEMAP", "PPIMAGE.BIAS", "PPIMAGE.DARK", Index: branches/pap/ppImage/src/ppImageDetrendFringe.c =================================================================== --- branches/pap/ppImage/src/ppImageDetrendFringe.c (revision 23948) +++ branches/pap/ppImage/src/ppImageDetrendFringe.c (revision 25027) @@ -11,15 +11,9 @@ PS_ASSERT_PTR_NON_NULL(options, false); - // Reference fringe measurements - psArray *references = psMemIncrRefCounter(psMetadataLookupPtr(NULL, fringe->analysis, - "FRINGE.MEASUREMENTS")); + // Reference fringe measurements (stored on the reference cell->analysis) + psArray *references = psMemIncrRefCounter(psMetadataLookupPtr(NULL, fringe->analysis, "FRINGE.MEASUREMENTS")); if (!references) { - references = pmFringesParse(fringe); // Reference fringes - if (!references) { - psError(PS_ERR_IO, false, "Unable to find fringe references.\n"); - return false; - } - psMetadataAdd(fringe->analysis, PS_LIST_TAIL, "FRINGE.MEASUREMENTS", PS_DATA_UNKNOWN, - "Fringe measurements", references); + psError(PS_ERR_IO, false, "Unable to find fringe references.\n"); + return false; } @@ -75,6 +69,5 @@ } - psMetadataAdd(readout->parent->analysis, PS_LIST_TAIL, scienceFringes, - PS_DATA_UNKNOWN | PS_META_REPLACE, "Fringe measurements", measurements); + psMetadataAdd(readout->parent->analysis, PS_LIST_TAIL, scienceFringes, PS_DATA_UNKNOWN | PS_META_REPLACE, "Fringe measurements", measurements); psFree(measurements); psFree(scienceFringes); @@ -92,10 +85,12 @@ psArray *fringes = psArrayAlloc(cells->n); // Fringes, to return for (int i = 0; i < cells->n; i++) { + fringes->data[i] = NULL; + pmCell *cell = cells->data[i]; // Cell of interest // XXX for now, skip the video cells (cell->readouts->n > 1) if (cell->readouts->n > 1) { - psWarning ("Skipping Video Cell for ppImageDetrendFringe.c:getFringes"); - continue; + psWarning ("Skipping Video Cell (%d) for ppImageDetrendFringe.c:getFringes", i); + continue; } @@ -197,4 +192,6 @@ } assert(fringes->readouts->n == solution->nFringeFrames); + + psImageMaskType maskVal = options->flatMask; bool mdok; // Status of MD lookup @@ -262,9 +259,19 @@ } + // subtract fringe and update mask if fringe value is NAN + for (int iy = 0; iy < readout->image->numRows; iy++) { + for (int ix = 0; ix < readout->image->numCols; ix++) { + readout->image->data.F32[iy][ix] -= sumFringe->data.F32[iy][ix]; + if (!isfinite(sumFringe->data.F32[iy][ix]) && readout->mask) { + readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= maskVal; + } + } + } + // XXX: Make generic, so subregions may be subtracted as well - if (!psBinaryOp(readout->image, readout->image, "-", sumFringe)) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to subtract fringe.\n"); - return false; - } + // if (!psBinaryOp(readout->image, readout->image, "-", sumFringe)) { + // psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to subtract fringe.\n"); + // return false; + // } // measure residual fringe amplitude. results go to FRINGE.RESIDUALS Index: branches/pap/ppImage/src/ppImageDetrendReadout.c =================================================================== --- branches/pap/ppImage/src/ppImageDetrendReadout.c (revision 23948) +++ branches/pap/ppImage/src/ppImageDetrendReadout.c (revision 25027) @@ -17,6 +17,6 @@ // Masking on the basis of pixel value needs to be done before anything else, so the values are pristine. if (options->doMaskBuild) { - psImageMaskType satMask = options->doMaskSat ? options->satMask : 0; - psImageMaskType lowMask = options->doMaskLow ? options->lowMask : 0; + psImageMaskType satMask = options->doMaskSat ? options->satMask : 0; + psImageMaskType lowMask = options->doMaskLow ? options->lowMask : 0; pmReadoutGenerateMask(input, satMask, lowMask); } @@ -56,5 +56,5 @@ if (!pmBiasSubtract(input, options->overscan, bias, oldDark, view)) { psError(PS_ERR_UNKNOWN, false, "Unable to subtract bias."); - psFree(detview); + psFree(detview); return false; } @@ -64,11 +64,33 @@ if (options->doVarianceBuild) { // create the target mask and variance images - pmReadoutGenerateVariance(input, true); + psImage *noiseImage = NULL; + if (options->doNoiseMap) { + // XXX convert the noiseMap image to a binned image + pmReadout *noiseMap = NULL; + noiseMap = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.NOISEMAP"); + noiseImage = psImageCopy (NULL, input->image, PS_TYPE_F32); + psImageInit (noiseImage, 0.0); + + // XXX this works, but is not really quite right: the model shoud include the + // offset information, we are not really getting exactly the right mapping from the + // original file. + psImageBinning *binning = psImageBinningAlloc(); + binning->nXruff = noiseMap->image->numCols; + binning->nYruff = noiseMap->image->numRows; + binning->nXfine = input->image->numCols; + binning->nYfine = input->image->numRows; + psImageBinningSetScale(binning, PS_IMAGE_BINNING_LEFT); + + psImageUnbin (noiseImage, noiseMap->image, binning); + psFree (binning); + } + pmReadoutGenerateVariance(input, noiseImage, true); + psFree (noiseImage); } if (options->doDark && dark) { - if (!pmDarkApply(input, dark, options->maskValue)) { + if (!pmDarkApply(input, dark, options->darkMask)) { psError(PS_ERR_UNKNOWN, false, "Unable to subtract dark."); - psFree(detview); + psFree(detview); return false; } @@ -79,5 +101,5 @@ options->remnanceSize, options->remnanceThresh)) { psError(PS_ERR_UNKNOWN, false, "Unable to mask remnance."); - psFree(detview); + psFree(detview); return false; } @@ -88,5 +110,5 @@ pmReadout *shutter = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.SHUTTER"); if (!pmShutterCorrectionApply(input, shutter, pmConfigMaskGet("FLAT", config))) { - psFree(detview); + psFree(detview); return false; } @@ -97,5 +119,15 @@ pmReadout *flat = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.FLAT"); if (!pmFlatField(input, flat, options->flatMask)) { - psFree(detview); + psFree(detview); + return false; + } + } + + // Pattern noise correction + if (options->doPattern) { + if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, + options->patternThresh, options->patternMean, options->patternStdev, + options->maskValue, options->darkMask)) { + psFree(detview); return false; } @@ -119,35 +151,35 @@ psMetadata *normlist = psMetadataLookupMetadata(&mdok, config->arguments, "NORMALIZATION.TABLE"); if (normlist) { - pmFPAfile *inputFile = psMetadataLookupPtr(&mdok, config->files, "PPIMAGE.INPUT"); - - // get the menu of class IDs - psMetadata *menu = psMetadataLookupMetadata(&mdok, inputFile->camera, "CLASSID"); + pmFPAfile *inputFile = psMetadataLookupPtr(&mdok, config->files, "PPIMAGE.INPUT"); + + // get the menu of class IDs + psMetadata *menu = psMetadataLookupMetadata(&mdok, inputFile->camera, "CLASSID"); if (!menu) { psError(PS_ERR_IO, false, "Unable to find CLASSID metadata in camera configuration"); - psFree(detview); - return false; - } - // get the rule for class_id for the desired class - const char *rule = psMetadataLookupStr(&mdok, menu, options->normClass); + psFree(detview); + return false; + } + // get the rule for class_id for the desired class + const char *rule = psMetadataLookupStr(&mdok, menu, options->normClass); if (!rule) { psError(PS_ERR_IO, false, "Unable to find NORM.CLASS value %s in CLASSID in camera configuration", options->normClass); - psFree(detview); - return false; - } - // get the class_id from the rule + psFree(detview); + return false; + } + // get the class_id from the rule char *classID = pmFPAfileNameFromRule(rule, inputFile, view); if (!classID) { psError(PS_ERR_IO, false, "error converting CLASSID rule %s to name\n", rule); - psFree(detview); - return false; - } - - // get normalization from the class_id - float norm = psMetadataLookupF32 (&mdok, normlist, classID); - if (!mdok) { - psError(PS_ERR_IO, false, "failed to find class ID %s in normalization table\n", classID); - psFree(detview); - return false; - } + psFree(detview); + return false; + } + + // get normalization from the class_id + float norm = psMetadataLookupF32 (&mdok, normlist, classID); + if (!mdok) { + psError(PS_ERR_IO, false, "failed to find class ID %s in normalization table\n", classID); + psFree(detview); + return false; + } pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest @@ -157,8 +189,8 @@ psFree(comment); - // apply the normalization + // apply the normalization psBinaryOp(input->image, input->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); - psFree (classID); + psFree (classID); } # endif @@ -167,5 +199,5 @@ pmCell *fringe = pmFPAfileThisCell(config->files, detview, "PPIMAGE.FRINGE"); if (!ppImageDetrendFringeMeasure(input, fringe, false, options)) { - psFree(detview); + psFree(detview); return false; } Index: branches/pap/ppImage/src/ppImageDetrendRecord.c =================================================================== --- branches/pap/ppImage/src/ppImageDetrendRecord.c (revision 23948) +++ branches/pap/ppImage/src/ppImageDetrendRecord.c (revision 25027) @@ -67,12 +67,11 @@ psMetadataAddMetadata(cell->analysis, PS_LIST_TAIL, "DETREND", 0, "Detrend information", detrend); - detrendRecord(options->doMask, detrend, config, view, "PPIMAGE.MASK", "DETREND.MASK", "Mask filename"); - detrendRecord(options->doBias, detrend, config, view, "PPIMAGE.BIAS", "DETREND.BIAS", "Bias filename"); - detrendRecord(options->doDark, detrend, config, view, "PPIMAGE.DARK", "DETREND.DARK", "Dark filename"); - detrendRecord(options->doShutter, detrend, config, view, "PPIMAGE.SHUTTER", "DETREND.SHUTTER", - "Shutter correction filename"); - detrendRecord(options->doFlat, detrend, config, view, "PPIMAGE.FLAT", "DETREND.FLAT", "Flat filename"); - detrendRecord(options->doFringe, detrend, config, view, "PPIMAGE.FRINGE", "DETREND.FRINGE", - "Fringe filename"); + detrendRecord(options->doMask, detrend, config, view, "PPIMAGE.MASK", "DETREND.MASK", "Mask filename"); + detrendRecord(options->doNoiseMap, detrend, config, view, "PPIMAGE.NOISEMAP", "DETREND.NOISEMAP", "Noise Map filename"); + detrendRecord(options->doBias, detrend, config, view, "PPIMAGE.BIAS", "DETREND.BIAS", "Bias filename"); + detrendRecord(options->doDark, detrend, config, view, "PPIMAGE.DARK", "DETREND.DARK", "Dark filename"); + detrendRecord(options->doShutter, detrend, config, view, "PPIMAGE.SHUTTER", "DETREND.SHUTTER", "Shutter correction filename"); + detrendRecord(options->doFlat, detrend, config, view, "PPIMAGE.FLAT", "DETREND.FLAT", "Flat filename"); + detrendRecord(options->doFringe, detrend, config, view, "PPIMAGE.FRINGE", "DETREND.FRINGE", "Fringe filename"); psFree (detrend); Index: branches/pap/ppImage/src/ppImageLoop.c =================================================================== --- branches/pap/ppImage/src/ppImageLoop.c (revision 23948) +++ branches/pap/ppImage/src/ppImageLoop.c (revision 25027) @@ -98,4 +98,12 @@ continue; } + + // perform the detrend analysis + if (!ppImageParityFlip(config, options, view)) { + ESCAPE("Unable to detrend readout"); + } + + // XXX TEST: + // psphotSaveImage (NULL, readout->image, "test.image.fits"); // XXX set the options->*Mask values here (after the mask images have been loaded Index: branches/pap/ppImage/src/ppImageMeasureCrosstalk.c =================================================================== --- branches/pap/ppImage/src/ppImageMeasureCrosstalk.c (revision 23948) +++ branches/pap/ppImage/src/ppImageMeasureCrosstalk.c (revision 25027) @@ -149,5 +149,8 @@ psStatsInit (stats); - psVectorStats (stats, vector, NULL, NULL, 0); + if (!psVectorStats (stats, vector, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } float background = psMetadataLookupF32 (&status, cell->analysis, "XTALK.REF"); Index: branches/pap/ppImage/src/ppImageOptions.c =================================================================== --- branches/pap/ppImage/src/ppImageOptions.c (revision 23948) +++ branches/pap/ppImage/src/ppImageOptions.c (revision 25027) @@ -25,4 +25,5 @@ options->doNonLin = false; // Non-linearity correction options->doOverscan = false; // Overscan subtraction + options->doNoiseMap = false; // Apply Read Noise Map options->doBias = false; // Bias subtraction options->doDark = false; // Dark subtraction @@ -30,4 +31,5 @@ options->doShutter = false; // Shutter correction options->doFlat = false; // Flat-field normalisation + options->doPattern = false; // Pattern noise subtraction options->doFringe = false; // Fringe subtraction options->doPhotom = false; // Source identification and photometry @@ -36,4 +38,5 @@ options->doStats = false; // Measure and save image statistics options->checkCTE = false; // Measure pixel-based variance + options->applyParity = false; // Apply Cell parities // output files requested @@ -58,4 +61,5 @@ options->lowMask = 0x00; // out-of-bounds (low) pixels (supplied to pmReadoutGenerateMask) options->flatMask = 0x00; // Bad flat pixels (supplied to pmFlatField) + options->darkMask = 0x00; // Bad dark pixels (supplied to pmDarkApply) options->blankMask = 0x00; // Blank (no data, cell gap) pixels (supplied to pmChipMosaic, pmFPAMosaic) options->markValue = 0x00; // A safe bit for internal marking @@ -83,4 +87,13 @@ options->fringeIter = 0; // Fringe iterations options->fringeKeep = 1.0; // Fringe keep fraction + + // Pattern noise values + + options->patternOrder = 0; // Polynomial order + options->patternIter = 0; // Clipping iterations + options->patternRej = NAN; // Clipping rejection threshold + options->patternThresh = NAN; // Threshold for ignoring pixels + options->patternMean = PS_STAT_NONE; // Statistic for mean + options->patternStdev = PS_STAT_NONE; // Statistic for standard deviation // Remnance values @@ -215,4 +228,5 @@ options->doCrosstalkCorrect = psMetadataLookupBool(NULL, recipe, "CROSSTALK.CORRECT"); + options->doNoiseMap = psMetadataLookupBool(NULL, recipe, "NOISEMAP"); options->doBias = psMetadataLookupBool(NULL, recipe, "BIAS"); options->doDark = psMetadataLookupBool(NULL, recipe, "DARK"); @@ -221,4 +235,5 @@ options->doFringe = psMetadataLookupBool(NULL, recipe, "FRINGE"); options->doShutter = psMetadataLookupBool(NULL, recipe, "SHUTTER"); + options->doPattern = psMetadataLookupBool(NULL, recipe, "PATTERN"); options->doStats = false; @@ -228,4 +243,6 @@ } + options->applyParity = psMetadataLookupBool(NULL, recipe, "APPLY.CELL.PARITY"); + // binned image options options->xBin1 = psMetadataLookupS32(&status, recipe, "BIN1.XBIN"); @@ -274,5 +291,5 @@ // 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 || + if (options->doNoiseMap || options->doBias || options->doOverscan || options->doDark || options->doShutter || options->doFlat || options->doPhotom) { options->doMaskBuild = true; @@ -292,4 +309,12 @@ options->fringeKeep = psMetadataLookupF32(NULL, recipe, "FRINGE.KEEP"); + // Pattern noise + options->patternOrder = psMetadataLookupS32(NULL, recipe, "PATTERN.ORDER"); + options->patternIter = psMetadataLookupS32(NULL, recipe, "PATTERN.ITER"); + options->patternRej = psMetadataLookupF32(NULL, recipe, "PATTERN.REJ"); + options->patternThresh = psMetadataLookupF32(NULL, recipe, "PATTERN.THRESH"); + options->patternMean = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "PATTERN.MEAN")); + options->patternStdev = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "PATTERN.STDEV")); + // Remnance options options->remnanceSize = psMetadataLookupF32(NULL, recipe, "REMNANCE.SIZE"); Index: branches/pap/ppImage/src/ppImageParityFlip.c =================================================================== --- branches/pap/ppImage/src/ppImageParityFlip.c (revision 25027) +++ branches/pap/ppImage/src/ppImageParityFlip.c (revision 25027) @@ -0,0 +1,127 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "ppImage.h" + +// Update a concept to the assumed value +#define FIX_CONCEPT(SOURCE, NAME, TYPE, VALUE) { \ +psMetadataItem *item = psMetadataLookup(SOURCE, NAME); \ +item->data.TYPE = VALUE; } + +bool ppImageParityFlip (pmConfig *config, const ppImageOptions *options, const pmFPAview *view) { + + bool status; + + if (!options->applyParity) return true; + + // find the currently selected readout + pmReadout *readout = pmFPAfileThisReadout(config->files, view, "PPIMAGE.INPUT"); + + pmCell *cell = readout->parent; + + int xParity = psMetadataLookupS32(&status, cell->concepts, "CELL.XPARITY"); + if (!status || (xParity != -1 && xParity != 1)) { + psWarning("CELL.XPARITY is not set for the input cell; assuming 1.\n"); + FIX_CONCEPT(cell->concepts, "CELL.XPARITY", S32, 1); + xParity = 1; + } + int yParity = psMetadataLookupS32(&status, cell->concepts, "CELL.YPARITY"); + if (!status || (yParity != -1 && yParity != 1)) { + psWarning("CELL.YPARITY is not set for the input cell; assuming 1.\n"); + FIX_CONCEPT(cell->concepts, "CELL.YPARITY", S32, 1); + yParity = 1; + } + + // for this case, nothing to be done: + if ((xParity == 1) && (yParity == 1)) return true; + + psImage *image = readout->image; + psImage *var = readout->variance; + psImage *mask = readout->mask; + + if (var) psAssert (image->numCols == var->numCols, "image and variance sizes are mismatched"); + if (var) psAssert (image->numRows == var->numRows, "image and variance sizes are mismatched"); + if (mask) psAssert (image->numCols == mask->numCols, "image and mask sizes are mismatched"); + if (mask) psAssert (image->numRows == mask->numRows, "image and mask sizes are mismatched"); + + int Nx = image->numCols; + int Ny = image->numRows; + + // we need working buffers for the image, variance, and mask: + psF32 *imrow = (psF32 *) psAlloc (Nx*sizeof(psF32)); + psF32 *wtrow = (psF32 *) psAlloc (Nx*sizeof(psF32)); + psImageMaskType *mkrow = (psImageMaskType *) psAlloc (Nx*sizeof(psImageMaskType)); + + // the three cases (+1,-1), (-1,+1), (-1,-1) should be handled independently + if ((xParity == -1) && (yParity == +1)) { + for (int iy = 0; iy < Ny; iy++) { + for (int ix = 0; ix < Nx; ix++) { + imrow[Nx-1-ix] = image->data.F32[iy][ix]; + } + for (int ix = 0; (var != NULL) && (ix < Nx); ix++) { + wtrow[Nx-1-ix] = var->data.F32[iy][ix]; + } + for (int ix = 0; (mask != NULL) && (ix < Nx); ix++) { + mkrow[Nx-1-ix] = mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]; + } + memcpy (image->data.F32[iy], imrow, Nx*sizeof(psF32)); + if (var) memcpy (var->data.F32[iy], wtrow, Nx*sizeof(psF32)); + if (mask) memcpy (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy], mkrow, Nx*sizeof(psImageMaskType)); + } + } + if ((xParity == +1) && (yParity == -1)) { + for (int iy = 0; iy < Ny/2; iy++) { + memcpy (imrow, image->data.F32[iy], Nx*sizeof(psF32)); + memcpy (image->data.F32[iy], image->data.F32[Ny-1-iy], Nx*sizeof(psF32)); + memcpy (image->data.F32[Ny-1-iy], imrow, Nx*sizeof(psF32)); + + if (var) { + memcpy (wtrow, var->data.F32[iy], Nx*sizeof(psF32)); + memcpy (var->data.F32[iy], var->data.F32[Ny-1-iy], Nx*sizeof(psF32)); + memcpy (var->data.F32[Ny-1-iy], wtrow, Nx*sizeof(psF32)); + } + + if (mask) { + memcpy (mkrow, mask->data.PS_TYPE_IMAGE_MASK_DATA[iy], Nx*sizeof(psImageMaskType)); + memcpy (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy], mask->data.PS_TYPE_IMAGE_MASK_DATA[Ny-1-iy], Nx*sizeof(psImageMaskType)); + memcpy (mask->data.PS_TYPE_IMAGE_MASK_DATA[Ny-1-iy], mkrow, Nx*sizeof(psImageMaskType)); + } + } + } + if ((xParity == -1) && (yParity == -1)) { + for (int iy = 0; iy < Ny/2; iy++) { + for (int ix = 0; ix < Nx; ix++) { + imrow[Nx-1-ix] = image->data.F32[iy][ix]; + } + for (int ix = 0; var && (ix < Nx); ix++) { + wtrow[Nx-1-ix] = var->data.F32[iy][ix]; + } + for (int ix = 0; mask && (ix < Nx); ix++) { + mkrow[Nx-1-ix] = mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]; + } + + for (int ix = 0; ix < Nx; ix++) { + image->data.F32[iy][Nx-1-ix] = image->data.F32[Ny-1-iy][ix]; + } + for (int ix = 0; var && (ix < Nx); ix++) { + var->data.F32[iy][Nx-1-ix] = var->data.F32[Ny-1-iy][ix]; + } + for (int ix = 0; mask && (ix < Nx); ix++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][Nx-1-ix] = mask->data.PS_TYPE_IMAGE_MASK_DATA[Ny-1-iy][ix]; + } + + memcpy (image->data.F32[Ny-1-iy], imrow, Nx*sizeof(psF32)); + if (var) memcpy (var->data.F32[Ny-1-iy], wtrow, Nx*sizeof(psF32)); + if (mask) memcpy (mask->data.PS_TYPE_IMAGE_MASK_DATA[Ny-1-iy], mkrow, Nx*sizeof(psImageMaskType)); + } + } + + // FIX_CONCEPT(cell->concepts, "CELL.XPARITY", S32, 1); + // FIX_CONCEPT(cell->concepts, "CELL.YPARITY", S32, 1); + + psFree (imrow); + psFree (wtrow); + psFree (mkrow); + return true; +} Index: branches/pap/ppImage/src/ppImageParseCamera.c =================================================================== --- branches/pap/ppImage/src/ppImageParseCamera.c (revision 23948) +++ branches/pap/ppImage/src/ppImageParseCamera.c (revision 25027) @@ -25,4 +25,12 @@ // otherwise they revert to the config information // not all input or output images are used in a given recipe + if (options->doNoiseMap) { + if (!ppImageDefineFile(config, input->fpa, "PPIMAGE.NOISEMAP", "NOISEMAP", + PM_FPA_FILE_IMAGE, PM_DETREND_TYPE_NOISEMAP)) { + psError(PS_ERR_IO, false, "Can't find a noise map image source"); + psFree(options); + return NULL; + } + } if (options->doBias) { if (!ppImageDefineFile(config, input->fpa, "PPIMAGE.BIAS", "BIAS", @@ -48,30 +56,4 @@ return NULL; } - -#if 0 - // I think this is now done automatically in the pmFPAfileDefine and pmFPAfileIOChecks -- PAP. - - // XXX have ppImageDefineFile return the pmFPAfile? - pmFPAfile *mask = psMetadataLookupPtr(&status, config->files, "PPIMAGE.MASK"); - psAssert(mask, "Just defined the mask!"); - - // Need to read the names of bit masks from the mask header and set them in the - // recipe. If we are loading this from the detrend db, this action will happen - // when the file is resolved. - if (!mask->detrend) { - // XXX need to load the mask bit names from one of the headers - // this grabs the first available hdu : no guarantee that it will be valid, though - pmHDU *hdu = pmHDUGetFirst(mask->fpa); - if (!hdu) { - psError(PS_ERR_IO, true, "no valid HDU for PPIMAGE.INPUT.MASK"); - return NULL; - } - // XXX is this consistent with the pmConfigMaskReadHeader call above? - if (!pmConfigMaskReadHeader(config, hdu->header)) { - psError(PS_ERR_IO, false, "error in mask bits"); - return NULL; - } - } -#endif } if (options->doShutter) { Index: branches/pap/ppImage/src/ppImageReplaceBackground.c =================================================================== --- branches/pap/ppImage/src/ppImageReplaceBackground.c (revision 23948) +++ branches/pap/ppImage/src/ppImageReplaceBackground.c (revision 25027) @@ -144,11 +144,15 @@ for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - float value = backData[y][x]; - if (!isfinite(value)) { - image->data.F32[y][x] = NAN; - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= options->lowMask; + if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) { + image->data.F32[y][x] = 0.0; } else { - image->data.F32[y][x] -= value; - } + float value = backData[y][x]; + if (!isfinite(value)) { + image->data.F32[y][x] = NAN; + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= options->lowMask; + } else { + image->data.F32[y][x] -= value; + } + } } } Index: branches/pap/ppImage/src/ppImageSetMaskBits.c =================================================================== --- branches/pap/ppImage/src/ppImageSetMaskBits.c (revision 23948) +++ branches/pap/ppImage/src/ppImageSetMaskBits.c (revision 25027) @@ -14,7 +14,11 @@ // at this point we know we have valid values for required entries SAT, BAD, FLAT, BLANK: - // mask for non-linear flat corrections + // mask for bad flat corrections options->flatMask = pmConfigMaskGet("FLAT", config); psAssert (options->flatMask, "flat mask not set"); + + // mask for bad dark corrections + options->darkMask = pmConfigMaskGet("DARK", config); + psAssert (options->darkMask, "dark mask not set"); // mask for non-existent data (default to DETECTOR if not defined)