Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 21082) +++ trunk/psphot/src/psphot.h (revision 21108) @@ -213,5 +213,5 @@ bool psphotKapaClose (); bool psphotImageBackgroundCellHistogram (psVector *values, float mean, float sigma, int ix, int iy); -bool psphotDiagnosticPlots (pmConfig *config, char *name, ...); +bool psphotDiagnosticPlots (const pmConfig *config, const char *name, ...); bool psphotMakeFluxScale (psImage *image, psMetadata *recipe, pmPSF *psf); @@ -222,3 +222,11 @@ int psphotSupplementStars (psArray *stars, psArray *sources, psImageBinning *binning, int ix, int iy); + +// Create a background model for a readout +// Essentially identical to psphotModelBackground, but with a more convenient API for outsiders. +psImage *psphotBackgroundModel(pmReadout *ro, // Readout for which to generate background model + const pmConfig *config // Configuration + ); + + #endif Index: trunk/psphot/src/psphotDiagnosticPlots.c =================================================================== --- trunk/psphot/src/psphotDiagnosticPlots.c (revision 21082) +++ trunk/psphot/src/psphotDiagnosticPlots.c (revision 21108) @@ -135,5 +135,5 @@ } -bool psphotDiagnosticPlots (pmConfig *config, char *name, ...) { +bool psphotDiagnosticPlots (const pmConfig *config, const char *name, ...) { va_list argPtr; Index: trunk/psphot/src/psphotModelBackground.c =================================================================== --- trunk/psphot/src/psphotModelBackground.c (revision 21082) +++ trunk/psphot/src/psphotModelBackground.c (revision 21108) @@ -1,20 +1,43 @@ # include "psphotInternal.h" static int npass = 0; - -// 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) +static char *defaultStatsName = "FITTED_MEAN"; + +// 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 + const pmConfig *config // Configuration + ) { bool status = true; - static char *defaultStatsName = "FITTED_MEAN"; - + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // 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, "BACKGROUND.XBIN"); + binning->nYbin = psMetadataLookupS32(&status, recipe, "BACKGROUND.YBIN"); + + psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkip(binning, image); + + return binning; +} + +// Generate the background model +static bool backgroundModel(psImage *model, // Model image + psImage *modelStdev, // Model stdev image + psImage *image, // Image for which to generate a background model + psImage *mask, // Mask for image + psImageBinning *binning, // Binning parameters + const pmConfig *config // Configuration + ) +{ psTimerStart ("psphot.background"); - // find the currently selected readout - pmFPAfile *file = psMetadataLookupPtr (&status, config->files, filename); - pmFPA *inFPA = file->fpa; - pmReadout *readout = pmFPAviewThisReadout (view, inFPA); - psImage *image = readout->image; - psImage *mask = readout->mask; + bool status = true; // select the appropriate recipe information @@ -116,22 +139,11 @@ *statsDefaults = *stats; - // 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, "BACKGROUND.XBIN"); - binning->nYbin = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); - - psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER); - psImageBinningSetSkip(binning, image); + // we save the binning structure for use in psphotMagnitudes status = psMetadataAddPtr(recipe, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); PS_ASSERT (status, false); - // we save the binning structure for use in psphotMagnitudes - pmReadout *model = pmFPAGenerateReadout (config, view, "PSPHOT.BACKMDL", inFPA, binning); - pmReadout *modelStdev = pmFPAGenerateReadout (config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning); - - psF32 **modelData = model->image->data.F32; - psF32 **modelStdevData = modelStdev->image->data.F32; + + psF32 **modelData = model->data.F32; + psF32 **modelStdevData = modelStdev->data.F32; // XXXX we can thread this here by running blocks in parallel @@ -140,6 +152,6 @@ psRegion ruffRegion = {0,0,0,0}; psRegion fineRegion = {0,0,0,0}; - for (int iy = 0; iy < model->image->numRows; iy++) { - for (int ix = 0; ix < model->image->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + for (int ix = 0; ix < model->numCols; ix++) { // convert the ruff grid cell to the equivalent fine grid cell @@ -229,5 +241,5 @@ char name[256]; sprintf (name, "backraw.%02d.fits", npass); - psphotSaveImage (NULL, model->image, name); + psphotSaveImage (NULL, model, name); } @@ -237,6 +249,6 @@ double Value = 0; // sum of good pixel's value double ValueStdev = 0; // sum of good pixel's standard deviations - for (int iy = 0; iy < model->image->numRows; iy++) { - for (int ix = 0; ix < model->image->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + for (int ix = 0; ix < model->numCols; ix++) { if (!isnan(modelData[iy][ix])) { Value += modelData[iy][ix]; @@ -250,9 +262,9 @@ for (int jy = iy - 1; jy <= iy + 1; jy++) { if (jy < 0) continue; - if (jy >= model->image->numRows) continue; + if (jy >= model->numRows) continue; for (int jx = ix - 1; jx <= ix + 1; jx++) { if (!jx && !jy) continue; if (jx < 0) continue; - if (jx >= model->image->numCols) continue; + if (jx >= model->numCols) continue; value += modelData[jy][jx]; count += 1.0; @@ -275,6 +287,6 @@ // patch over remaining bad regions (use global average) - for (int iy = 0; iy < model->image->numRows; iy++) { - for (int ix = 0; ix < model->image->numCols; ix++) { + for (int iy = 0; iy < model->numRows; iy++) { + for (int ix = 0; ix < model->numCols; ix++) { if (!isnan(modelData[iy][ix])) continue; modelData[iy][ix] = Value; @@ -294,5 +306,5 @@ PS_STAT_MIN | PS_STAT_MAX); - psImageStats (statsBck, model->image, NULL, 0); + psImageStats (statsBck, model, NULL, 0); psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean", statsBck->sampleMean); @@ -304,7 +316,7 @@ PS_META_REPLACE, "sky model minimum value", statsBck->min); psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NX", - PS_META_REPLACE, "sky model size (x)", model->image->numCols); + PS_META_REPLACE, "sky model size (x)", model->numCols); psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NY", - PS_META_REPLACE, "sky model size (y)", model->image->numRows); + 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); @@ -316,4 +328,54 @@ psFree(rng); + return model; +} + + + +psImage *psphotBackgroundModel(pmReadout *ro, const pmConfig *config) +{ + PM_ASSERT_READOUT_NON_NULL(ro, NULL); + PM_ASSERT_READOUT_IMAGE(ro, NULL); + PS_ASSERT_PTR_NON_NULL(config, NULL); + + psImageBinning *binning = backgroundBinning(ro->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->image, ro->mask, binning, config)) { + psFree(model); + psFree(modelStdev); + psError(PS_ERR_UNKNOWN, false, "Unable to generate background model"); + 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) +{ + bool status = true; + + // find the currently selected readout + pmFPAfile *file = psMetadataLookupPtr (&status, config->files, filename); + pmFPA *inFPA = file->fpa; + pmReadout *readout = pmFPAviewThisReadout (view, inFPA); + psImage *image = readout->image; + psImage *mask = readout->mask; + + psImageBinning *binning = backgroundBinning(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, image, mask, binning, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate background model"); + return false; + } + npass ++; return true;