Index: branches/simmosaic_branches/psphot/src/psphotModelBackground.c =================================================================== --- branches/simmosaic_branches/psphot/src/psphotModelBackground.c (revision 24860) +++ branches/simmosaic_branches/psphot/src/psphotModelBackground.c (revision 27839) @@ -5,5 +5,5 @@ // Determine the binning characteristics for the background model // Sets the ruff image size and skip based on recipe values for the binning -static psImageBinning *backgroundBinning(const psImage *image, // Image for which to generate a bg model +psImageBinning *psphotBackgroundBinning(const psImage *image, // Image for which to generate a bg model const pmConfig *config // Configuration ) @@ -28,12 +28,20 @@ } + // Generate the background model -static bool backgroundModel(psImage *model, // Model image - psImage *modelStdev, // Model stdev image - psMetadata *analysis, // Analysis metadata for outputs - pmReadout *readout, // Readout for which to generate a background model - psImageBinning *binning, // Binning parameters - const pmConfig *config // Configuration - ) +// generate the median in NxN boxes, clipping heavily +// linear interpolation to generate full-scale model +// +// NOTE that the 'analysis' metedata pass in here is used to store the binning information. +// This may be the analysis for this readout, but it may be the analysis for the pmFPAfile +// corresponding to the model. Other information about the background model is saved on the +// readout->analysis +static bool psphotModelBackgroundReadout(psImage *model, // Model image + psImage *modelStdev, // Model stdev image + psMetadata *analysis, // Analysis metadata for outputs + pmReadout *readout, // Readout for which to generate a background model + psImageBinning *binning, // Binning parameters + const pmConfig *config // Configuration + ) { psTimerStart ("psphot.background"); @@ -137,6 +145,5 @@ // we save the binning structure for use in psphotMagnitudes - psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", - PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); + psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); psF32 **modelData = model->data.F32; @@ -144,4 +151,7 @@ // XXXX we can thread this here by running blocks in parallel + + int nFailures = 0; + psImageBackgroundInit(); // measure clipped median for subimages @@ -211,7 +221,10 @@ stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV; if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) { - psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for " - "(%dx%d, (row0,col0) = (%d,%d)", - subset->numRows, subset->numCols, subset->row0, subset->col0); + if ((nFailures < 3) || (nFailures % 100 == 0)) { + psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for " + "(%dx%d, (row0,col0) = (%d,%d)", + subset->numRows, subset->numCols, subset->row0, subset->col0); + } + nFailures ++; modelData[iy][ix] = modelStdevData[iy][ix] = NAN; } else { @@ -232,4 +245,8 @@ psFree (submask); } + } + + if (nFailures) { + psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background for %d of %d subimages", nFailures, (model->numRows*model->numCols)); } @@ -293,6 +310,6 @@ psLogMsg ("psphot", PS_LOG_INFO, "built background image: %f sec\n", psTimerMark ("psphot.background")); - psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value); - psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev); psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f stdev %f", Value, ValueStdev); @@ -303,16 +320,10 @@ PS_STAT_MAX); psImageStats (statsBck, model, NULL, 0); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MN", - PS_META_REPLACE, "sky model mean", statsBck->sampleMean); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_SIG", - PS_META_REPLACE, "sky model stdev", statsBck->sampleStdev); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MAX", - PS_META_REPLACE, "sky model maximum value", statsBck->max); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MIN", - PS_META_REPLACE, "sky model minimum value", statsBck->min); - psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NX", - PS_META_REPLACE, "sky model size (x)", model->numCols); - psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NY", - PS_META_REPLACE, "sky model size (y)", model->numRows); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean", statsBck->sampleMean); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev", statsBck->sampleStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)", model->numCols); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)", model->numRows); psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f", statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev); @@ -327,17 +338,16 @@ } - - -psImage *psphotBackgroundModel(pmReadout *ro, const pmConfig *config) +// generate a background model for a single readout. do not save an associated pmFPAfile for possible output +psImage *psphotModelBackgroundReadoutNoFile(pmReadout *readout, const pmConfig *config) { - PM_ASSERT_READOUT_NON_NULL(ro, NULL); - PM_ASSERT_READOUT_IMAGE(ro, NULL); + PM_ASSERT_READOUT_NON_NULL(readout, NULL); + PM_ASSERT_READOUT_IMAGE(readout, NULL); PS_ASSERT_PTR_NON_NULL(config, NULL); - psImageBinning *binning = backgroundBinning(ro->image, config); // Image binning parameters + psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters psImage *model = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Background model psImage *modelStdev = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Standard deviation - if (!backgroundModel(model, modelStdev, ro->analysis, ro, binning, config)) { + if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config)) { psFree(model); psFree(modelStdev); @@ -345,27 +355,24 @@ return NULL; } - psFree(modelStdev); - return model; } - -// generate the median in NxN boxes, clipping heavily -// linear interpolation to generate full-scale model -bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filename) +// generate a background model for readout number index; save an associated pmFPAfile for possible output +bool psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filename, int index) { - bool status = true; - // find the currently selected readout - pmFPAfile *file = psMetadataLookupPtr (&status, config->files, filename); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + pmFPA *inFPA = file->fpa; pmReadout *readout = pmFPAviewThisReadout(view, inFPA); - - psImageBinning *binning = backgroundBinning(readout->image, config); // Image binning parameters - pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning); - pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning); - - if (!backgroundModel(model->image, modelStdev->image, model->analysis, readout, binning, config)) { + psAssert (readout, "missing readout?"); + + psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters + pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning, index); + pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index); + + if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate background model"); return false; @@ -375,2 +382,20 @@ return true; } + +// XXX supply filename or keep PSPHOT.INPUT fixed? +bool psphotModelBackground (pmConfig *config, const pmFPAview *view) +{ + bool status = false; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotModelBackgroundReadoutFileIndex(config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +}