Index: trunk/ppImage/src/ppImage.h =================================================================== --- trunk/ppImage/src/ppImage.h (revision 20439) +++ trunk/ppImage/src/ppImage.h (revision 20444) @@ -35,8 +35,8 @@ bool doFringe; // Fringe subtraction bool doPhotom; // Source identification and photometry + bool doBG; // Background subtraction bool doAstromChip; // per-chip Astrometry bool doAstromMosaic; // full-mosaic Astrometry bool doStats; // call ppStats on the image - bool replaceMasked; // fill in masked values with background model // output files requested @@ -134,5 +134,4 @@ bool ppImageRebinChip (pmConfig *config, pmFPAview *view, ppImageOptions *options, char *outName); -bool ppImageRebinReadout (pmReadout *output, pmReadout *input, pmFPAfile *outFile, ppImageOptions *options); bool ppImagePhotom (pmConfig *config, pmFPAview *view); @@ -140,5 +139,10 @@ bool ppImageAddstar (pmConfig *config); -bool ppImageReplaceBackground (pmConfig *config, pmFPAview *view, ppImageOptions *options); +// Subtract background from the chip-mosaicked image +bool ppImageSubtractBackground( + pmConfig *config, // Configuration + const pmFPAview *view, // View to chip of interest + const ppImageOptions *options // Processing options + ); bool ppImageMosaicChip (pmConfig *config, const ppImageOptions *options, const pmFPAview *view, Index: trunk/ppImage/src/ppImageLoop.c =================================================================== --- trunk/ppImage/src/ppImageLoop.c (revision 20439) +++ trunk/ppImage/src/ppImageLoop.c (revision 20444) @@ -129,7 +129,7 @@ // replace the masked pixels with the background level - if (options->replaceMasked) { - if (!ppImageReplaceBackground(config, view, options)) { - ESCAPE("Unable to replace masked pixels with background level"); + if (options->doBG) { + if (!ppImageSubtractBackground(config, view, options)) { + ESCAPE("Unable to subtract background"); } } Index: trunk/ppImage/src/ppImageOptions.c =================================================================== --- trunk/ppImage/src/ppImageOptions.c (revision 20439) +++ trunk/ppImage/src/ppImageOptions.c (revision 20444) @@ -247,4 +247,5 @@ options->doAstromChip = psMetadataLookupBool(NULL, recipe, "ASTROM.CHIP"); options->doAstromMosaic = psMetadataLookupBool(NULL, recipe, "ASTROM.MOSAIC"); + options->doBG = psMetadataLookupBool(NULL, recipe, "BACKGROUND"); // even if not requested explicitly, if any of these are set, build an internal mask and weight: @@ -253,5 +254,5 @@ options->doMaskBuild = true; options->doWeightBuild = true; - } else if (options->doMask) { + } else if (options->doMask || options->doBG) { options->doMaskBuild = true; } @@ -267,5 +268,4 @@ options->fringeKeep = psMetadataLookupF32(NULL, recipe, "FRINGE.KEEP"); - options->replaceMasked = psMetadataLookupBool(NULL, recipe, "REPLACE.MASKED"); return options; Index: trunk/ppImage/src/ppImageParseCamera.c =================================================================== --- trunk/ppImage/src/ppImageParseCamera.c (revision 20439) +++ trunk/ppImage/src/ppImageParseCamera.c (revision 20444) @@ -74,25 +74,25 @@ return NULL; } - // XXX have ppImageDefineFile return the pmFPAfile? - pmFPAfile *mask = psMetadataLookupPtr(&status, config->files, "PPIMAGE.MASK"); - psAssert (mask, "mask not defined? not possible!"); - - // 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; - } - } + // XXX have ppImageDefineFile return the pmFPAfile? + pmFPAfile *mask = psMetadataLookupPtr(&status, config->files, "PPIMAGE.MASK"); + psAssert (mask, "mask not defined? not possible!"); + + // 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; + } + } } if (options->doShutter) { @@ -114,6 +114,6 @@ int nThreads = psMetadataLookupS32 (&status, config->arguments, "NTHREADS"); if (nThreads > 0) { - int nScanRows = psMetadataLookupS32 (&status, recipe, "SCAN.ROWS"); - pmDetrendSetThreadTasks (nScanRows); + int nScanRows = psMetadataLookupS32 (&status, recipe, "SCAN.ROWS"); + pmDetrendSetThreadTasks (nScanRows); } @@ -274,5 +274,5 @@ // For photometry, we operate on the chip-mosaicked image // we create a copy of the mosaicked image for psphot so we can write out a clean image - if (options->doPhotom) { + if (options->doPhotom || options->doBG) { // this file is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by Index: trunk/ppImage/src/ppImageReplaceBackground.c =================================================================== --- trunk/ppImage/src/ppImageReplaceBackground.c (revision 20439) +++ trunk/ppImage/src/ppImageReplaceBackground.c (revision 20444) @@ -3,108 +3,117 @@ #endif -# include "ppImage.h" - -enum {MODE_NONE, MODE_MODEL, MODE_VALUE}; +#include "ppImage.h" // In this function, we perform the psphot analysis routine for the chip-mosaicked images -bool ppImageReplaceBackground (pmConfig *config, pmFPAview *view, ppImageOptions *options) { +bool ppImageSubtractBackground(pmConfig *config, const pmFPAview *view, const ppImageOptions *options) +{ + psAssert(config, "Need configuration"); + psAssert(view, "Need view to chip"); + psAssert(options, "Need options"); - bool status; - pmCell *cell; - pmReadout *readout; + if (!options->doBG) { + return true; + } - // find the reference source image - pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PPIMAGE.CHIP"); + bool status; // Status of MD lookup + pmFPAfile *input = psMetadataLookupPtr(&status, config->files, "PPIMAGE.CHIP"); // File to correct if (!status) { - psError(PSPHOT_ERR_CONFIG, false, "PSPHOT.INPUT I/O file is not defined"); + psError(PS_ERR_UNEXPECTED_NULL, true, "PPIMAGE.CHIP file is not defined"); return false; } - // select the appropriate recipe information - psMetadata *ppImageRecipe = psMetadataLookupPtr (&status, config->recipes, "PPIMAGE"); - - // select the appropriate recipe information (needed by psphotModelBackground) - psMetadata *psphotRecipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - - char *replaceSkyModeWord = psMetadataLookupStr (&status, ppImageRecipe, "REPLACE.MODE"); - float replaceSkyValue = 0.0; - - int replaceSkyMode = MODE_NONE; - if (!strcasecmp (replaceSkyModeWord, "MODEL")) { - replaceSkyMode = MODE_MODEL; - } - if (!strcasecmp (replaceSkyModeWord, "VALUE")) { - replaceSkyMode = MODE_VALUE; - replaceSkyValue = psMetadataLookupF32 (&status, ppImageRecipe, "REPLACE.VALUE"); - } - psAssert (replaceSkyMode != MODE_NONE, "replace sky mode not defined"); - - // XXX Should this be options->maskValue or & ~options->satMask? the latter will leave saturated pixels high - // psMaskType maskVal = options->maskValue & ~options->satMask; - psMaskType maskVal = options->maskValue; - - pmFPAfile *modelFile = NULL; - if (replaceSkyMode == MODE_MODEL) { - // find the model data, if it exists yet (if not, we build it below) - modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL"); - - // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) - psMetadataAddU8 (psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskVal); + // The background model file may be defined, though the model hasn't been built + // If so, we will build it below. + pmFPAfile *modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL"); // Background model + if (!status) { + psError(PS_ERR_UNEXPECTED_NULL, true, "PSPHOT.BACKMDL file is not defined"); + return false; } - // iterate over the cells and readout for this chip - while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { - psLogMsg ("ppImagePhotom", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); - if (! cell->process || ! cell->file_exists) { continue; } + psMetadata *ppImageRecipe = psMetadataLookupPtr(NULL, config->recipes, RECIPE_NAME); + psAssert(ppImageRecipe, "Need PPIMAGE recipe"); + psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE); + psAssert(psphotRecipe, "Need PSPHOT recipe"); - // process each of the readouts - while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { - if (! readout->data_exists) { continue; } - if (! readout->mask) { continue; } + // XXX Should this be options->maskValue or options->maskValue & ~options->satMask? + // The latter will leave saturated pixels high + psMaskType maskVal = options->maskValue; - // replace masked pixels with values from model (unbinning not needed) - pmReadout *modelReadout = NULL; - psImageBinning *binning = NULL; - psImage *model = NULL; - if (replaceSkyMode == MODE_MODEL) { - // we are using this pmFPAfile as an I/O file: select readout or create - if (modelFile) { - modelReadout = pmFPAviewThisReadout (view, modelFile->fpa); - } - if (!modelReadout) { - psphotModelBackground (config, view, "PPIMAGE.CHIP"); - modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL"); - assert (modelFile); - if (modelFile->mode == PM_FPA_MODE_INTERNAL) { - modelReadout = modelFile->readout; - } else { - modelReadout = pmFPAviewThisReadout (view, modelFile->fpa); - } - assert (modelReadout); - } - binning = psMetadataLookupPtr(&status, psphotRecipe, "PSPHOT.BACKGROUND.BINNING"); - model = modelReadout->image; - } + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskVal); - psImage *mask = readout->mask; - psImage *image = readout->image; + // Since we are working on a chip-mosaicked image, there should only be a single cell and readout + pmChip *chip = pmFPAviewThisChip(view, input->fpa); // Chip of interest + if (!chip) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find chip"); + return false; + } + if (chip->cells->n == 0) { + psWarning("Chip has no cells"); + return true; + } + if (chip->cells->n > 1) { + psWarning("Chip has %ld cells; only the first will be processed", chip->cells->n); + } + pmCell *cell = chip->cells->data[0]; // Cell of interest + if (!cell || !cell->process || !cell->file_exists) { + // Nothing to process + return true; + } + if (cell->readouts->n == 0) { + psWarning("Cell has no readouts"); + return true; + } + if (cell->readouts->n > 1) { + psWarning("Cell has %ld readouts; only the first will be processed", cell->readouts->n); + } + pmReadout *ro = cell->readouts->data[0]; // Readout of interest + if (!ro || !ro->data_exists) { + // Nothing to process + return true; + } + psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (!(mask->data.U8[iy][ix] && maskVal)) continue; - if (replaceSkyMode == MODE_MODEL) { - image->data.F32[iy][ix] = psImageUnbinPixel_V2(ix, iy, model, binning); - } else { - image->data.F32[iy][ix] = replaceSkyValue; - } - } - } - } + + pmFPAview roView = *view; // View to readout + roView.cell = roView.readout = 0; + pmReadout *modelRO = pmFPAviewThisReadout(&roView, modelFile->fpa); // Background model + if (!modelRO) { + // Create the background model + if (!psphotModelBackground(config, &roView, "PPIMAGE.CHIP")) { + psError(PS_ERR_UNKNOWN, false, "Unable to model background"); + return false; + } + modelRO = (modelFile->mode == PM_FPA_MODE_INTERNAL) ? modelFile->readout : + pmFPAviewThisReadout(&roView, modelFile->fpa); + if (!modelRO) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find background model"); + return false; + } + } + psImageBinning *binning = psMetadataLookupPtr(&status, psphotRecipe, + "PSPHOT.BACKGROUND.BINNING"); // Binning for model + psImage *modelImage = modelRO->image; // Background model + + // Do the background subtraction + int numCols = image->numCols, numRows = image->numRows; // Size of image + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (mask && mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) { + image->data.F32[y][x] = 0.0; + } else { + float value = psImageUnbinPixel_V2(x, y, modelImage, binning); // Background value + if (!isfinite(value)) { + image->data.F32[y][x] = NAN; + mask->data.PS_TYPE_MASK_DATA[y][x] |= options->badMask; + } else { + image->data.F32[y][x] -= value; + } + } + } } - if (replaceSkyMode == MODE_MODEL) { - pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); - pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); - } + pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL"); + pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV"); return true;