Index: /branches/eam_branches/ipp-20111122/psphot/src/psphotDefineFiles.c =================================================================== --- /branches/eam_branches/ipp-20111122/psphot/src/psphotDefineFiles.c (revision 33093) +++ /branches/eam_branches/ipp-20111122/psphot/src/psphotDefineFiles.c (revision 33094) @@ -85,4 +85,26 @@ output->save = true; } +// // optionally save the smoothed variance model (small FITS image) +// if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL")) { +// int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); +// int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); +// pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.VARMDL"); +// if (!output) { +// psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.VARMDL"); +// return false; +// } +// output->save = true; +// } +// // optionally save the smoothed variance model's standard deviation (small FITS image) +// if (psMetadataLookupBool(NULL, recipe, "SAVE.VARMDL.STDEV")) { +// int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); +// int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); +// pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.VARMDL.STDEV"); +// if (!output) { +// psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.VARMDL.STDEV"); +// return false; +// } +// output->save = true; +// } // optionally save the PSF Model if (psMetadataLookupBool(NULL, recipe, "SAVE.PSF")) { Index: /branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c (revision 33093) +++ /branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c (revision 33094) @@ -12,5 +12,5 @@ static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal); -bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, int index, psMetadata *recipe, pmReadout *readout, psArray *sources); +bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources); pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe); bool psphotFreeModelVariance (pmReadout *readout, psArray *sources); @@ -74,5 +74,5 @@ if (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) { // generate the model variance image & source pointers - if (!psphotGenerateModelVariance (config, view, i, recipe, readout, sources)) { + if (!psphotGenerateModelVariance (config, view, file, i, recipe, readout, sources)) { psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); return false; @@ -554,5 +554,14 @@ } -bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) { +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 + bool useVarianceImage + ); + +bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) { bool status = false; @@ -573,21 +582,28 @@ psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32); - // find the pmFPAfile for the readout (background model is saved on PSPHOT.BACKMDL regardless of 'filerule') - pmFPAfile *backModelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest - assert (backModelFile); - - // find the background model readout from the correct location - pmReadout *backModel = READOUT_OR_INTERNAL(view, backModelFile); - psAssert (backModel, "this must exist"); - // find the binning information - psImageBinning *backBinning = psMetadataLookupPtr(&status, backModel->analysis, "PSPHOT.BACKGROUND.BINNING"); + psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config); assert (backBinning); + + psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model + psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model + + if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate background model"); + psFree (varModel); + psFree (varModelStdev); + return false; + } // linear interpolation to full-scale - if (!psImageUnbin (modelVar, backModel->image, backBinning)) { - psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning"); - return false; - } + if (!psImageUnbin (modelVar, varModel, backBinning)) { + psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning"); + psFree (varModel); + psFree (varModelStdev); + return false; + } + + psFree (varModel); + psFree (varModelStdev); // XXX for a test: @@ -609,5 +625,6 @@ psAssert (!source->modelVar, "programming error : modelVar should be NULL here"); psAssert (source->modelFlux, "programming error : modelFlux should not be NULL here"); - source->modelVar = psImageSubset(readout->image, region); + psAssert (source->modelFlux->data.F32, "programming error : modelFlux should not be NULL here"); + source->modelVar = psImageSubset(modelVar, region); // add the source model to the model variance image @@ -615,11 +632,7 @@ } - // XXX add the readnoise to the image? This is a bit problematic because I don't have the - // per-cell information. Alternatively, I could use the median of the variance image or - // something equivalent to a background model to define the base level of the variance - // XXX for a test: psphotSaveImage (NULL, modelVar, "model.var.fits"); - psphotSaveImage (NULL, readout->image, "image.var.fits"); + psphotSaveImage (NULL, readout->variance, "image.var.fits"); // we save the model variance for future reference @@ -638,4 +651,6 @@ assert (modelVar); + psMetadataRemoveKey (readout->analysis, "PSPHOT.MODEL.VAR"); + // clear modelVar pointers for all of the source models for (int i = 0; i < sources->n; i++) { @@ -645,5 +660,4 @@ psFree (source->modelVar); } - psFree (modelVar); return true; Index: /branches/eam_branches/ipp-20111122/psphot/src/psphotModelBackground.c =================================================================== --- /branches/eam_branches/ipp-20111122/psphot/src/psphotModelBackground.c (revision 33093) +++ /branches/eam_branches/ipp-20111122/psphot/src/psphotModelBackground.c (revision 33094) @@ -37,10 +37,11 @@ // 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 +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 + bool useVarianceImage ) { @@ -49,5 +50,6 @@ bool status = true; - psImage *image = readout->image, *mask = readout->mask; // Image and mask for readout + psImage *image = useVarianceImage ? readout->variance : readout->image; + psImage *mask = readout->mask; // Image and mask for readout // select the appropriate recipe information @@ -143,5 +145,7 @@ // 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); + if (!useVarianceImage) { + psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); + } psVector *dQ = psVectorAllocEmpty (100, PS_TYPE_F32); @@ -330,25 +334,30 @@ psVectorStats (statsDQ, dQ, NULL, NULL, 0); - 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); - psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); - psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian); - - // measure image and background stats and save for later output - psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | - PS_STAT_SAMPLE_STDEV | - PS_STAT_MIN | - PS_STAT_MAX); - psImageStats (statsBck, model, NULL, 0); - 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_DEV", PS_META_REPLACE, "sky stdev", ValueStdev); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); - 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); + if (!useVarianceImage) { + 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); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); + psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian); + + // measure image and background stats and save for later output + psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | + PS_STAT_SAMPLE_STDEV | + PS_STAT_MIN | + PS_STAT_MAX); + psImageStats (statsBck, model, NULL, 0); + 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_DEV", PS_META_REPLACE, "sky stdev", ValueStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ", PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian); + 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); + psFree(statsBck); + } else { + psLogMsg ("psphot", PS_LOG_INFO, "variance data : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian); + } psFree(statsDQ); @@ -356,5 +365,4 @@ psFree(stats); - psFree(statsBck); psFree(statsDefaults); psFree(binning); @@ -375,5 +383,5 @@ psImage *modelStdev = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Standard deviation - if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config)) { + if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config, false)) { psFree(model); psFree(modelStdev); @@ -400,5 +408,5 @@ pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index); - if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config)) { + if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate background model"); return false;