Index: branches/eam_branches/ipp-20120601/psphot/src/psphot.h =================================================================== --- branches/eam_branches/ipp-20120601/psphot/src/psphot.h (revision 34051) +++ branches/eam_branches/ipp-20120601/psphot/src/psphot.h (revision 34053) @@ -61,4 +61,7 @@ bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index); @@ -99,10 +102,5 @@ bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final); - -# if (HAVE_MODEL_VAR) bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode); -# else -bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final); -# endif bool psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize); Index: branches/eam_branches/ipp-20120601/psphot/src/psphotEfficiency.c =================================================================== --- branches/eam_branches/ipp-20120601/psphot/src/psphotEfficiency.c (revision 34051) +++ branches/eam_branches/ipp-20120601/psphot/src/psphotEfficiency.c (revision 34053) @@ -420,9 +420,5 @@ // psphotFitSourcesLinearReadout subtracts the model fits -# if (HAVE_MODEL_VAR) if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) { -# else - if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) { -# endif psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); psFree(fakeSources); Index: branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c (revision 34051) +++ branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c (revision 34053) @@ -12,9 +12,7 @@ static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal); -# if (HAVE_MODEL_VAR) 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); -# endif // for now, let's store the detections on the readout->analysis for each readout @@ -30,5 +28,4 @@ assert (recipe); -# if (HAVE_MODEL_VAR) pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe); if (!fitVarMode) { @@ -40,5 +37,4 @@ // do a single pass. pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode; -# endif int num = psphotFileruleCount(config, filerule); @@ -68,15 +64,9 @@ psAssert (psf, "missing psf?"); -# if (HAVE_MODEL_VAR) - if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) -# else - if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) -# endif - { + if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) { psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); return false; } -# if (HAVE_MODEL_VAR) // the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance // after we have determined the initial set of fits, then we can generate the variance image and @@ -104,5 +94,4 @@ } } -# endif psphotVisualShowResidualImage (readout, (num > 0)); @@ -113,5 +102,4 @@ } -# if (HAVE_MODEL_VAR) // look up the fit variance mode from the recipe; older recipes do not have the value // 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and @@ -136,4 +124,7 @@ return PM_SOURCE_PHOTFIT_IMAGE_VAR; } + if (!strcasecmp(fitVarModeString, "SKY") || !strcasecmp(fitVarModeString, "MODEL_SKY")) { + return PM_SOURCE_PHOTFIT_MODEL_SKY; + } if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) { return PM_SOURCE_PHOTFIT_MODEL_VAR; @@ -142,12 +133,6 @@ return PM_SOURCE_PHOTFIT_NONE; } -# endif - -# if (HAVE_MODEL_VAR) -bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) -# else -bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) -# endif -{ + +bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) { bool status; float x; @@ -185,11 +170,4 @@ if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); -# if (!HAVE_MODEL_VAR) - bool CONSTANT_PHOTOMETRIC_WEIGHTS = - psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); - if (!status) { - psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); - } -# endif int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER"); if (!status) { @@ -319,5 +297,4 @@ psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder); -# if (HAVE_MODEL_VAR) // if fitVarMode is MODEL_VAR, then we need to generate the model image variance // XXX we have two possibilities here: @@ -327,5 +304,4 @@ // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?) -# endif // fill out the sparse matrix elements and border elements (B) @@ -336,14 +312,9 @@ // diagonal elements of the sparse matrix (auto-cross-product) -# if (HAVE_MODEL_VAR) f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal); -# else - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); -# endif psSparseMatrixElement (sparse, i, i, f); -# if (HAVE_MODEL_VAR) - // if we have used CONSTANT errors, then we need to calculate the value of the parameter error - if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) { + // if we have used CONSTANT errors, then we need to re-calculate the value of the parameter error + if (fitVarMode == PM_SOURCE_PHOTFIT_CONST) { float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal); errors->data.F32[i] = 1.0 / sqrt(var); @@ -351,20 +322,7 @@ errors->data.F32[i] = 1.0 / sqrt(f); } -# else - // the formal error depends on the weighting scheme - if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal); - errors->data.F32[i] = 1.0 / sqrt(var); - } else { - errors->data.F32[i] = 1.0 / sqrt(f); - } -# endif // find the image x model value -# if (HAVE_MODEL_VAR) f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal); -# else - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); -# endif psSparseVectorElement (sparse, i, f); @@ -372,22 +330,11 @@ switch (SKY_FIT_ORDER) { case 1: -# if (HAVE_MODEL_VAR) f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal); psSparseBorderElementB (border, i, 1, f); f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal); psSparseBorderElementB (border, i, 2, f); -# else - f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); - psSparseBorderElementB (border, i, 1, f); - f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); - psSparseBorderElementB (border, i, 2, f); -# endif case 0: -# if (HAVE_MODEL_VAR) f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal); -# else - f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); -# endif psSparseBorderElementB (border, i, 0, f); break; @@ -409,9 +356,5 @@ // got an overlap; calculate cross-product and add to output array -# if (HAVE_MODEL_VAR) f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal); -# else - f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); -# endif psSparseMatrixElement (sparse, j, i, f); } @@ -451,9 +394,5 @@ // set the sky, sky_x, sky_y components of border matrix -# if (HAVE_MODEL_VAR) SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal); -# else - SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal); -# endif psSparseConstraint constraint; @@ -614,5 +553,4 @@ } -# if (HAVE_MODEL_VAR) bool psphotModelBackgroundReadout(psImage *model, // Model image psImage *modelStdev, // Model stdev image @@ -650,4 +588,5 @@ psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model + // generate an image of the mean variance image in DN if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate background model"); @@ -668,6 +607,13 @@ psFree (varModelStdev); - // XXX for a test: - psphotSaveImage (NULL, modelVar, "model.bck.fits"); + float gain = 1.0; // accept 1.0 as a default since it is not critical to the analysis + pmCell *cell = readout->parent; // The parent cell + if (cell) { + gain = psMetadataLookupF32(&status, cell->concepts, "CELL.GAIN"); // Cell gain + if (!status) { + gain = 1.0; // set note above + } + } + if (gain > 2.0) { /* warn? */ } // insert all of the source models @@ -690,10 +636,9 @@ // add the source model to the model variance image + // XXX note that this should be added with gain applied + // var_DN = flux_DN / gain [e/DN] + // to do this requires an API upgrade... pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal); } - - // XXX for a test: - psphotSaveImage (NULL, modelVar, "model.var.fits"); - psphotSaveImage (NULL, readout->variance, "image.var.fits"); // we save the model variance for future reference @@ -724,3 +669,2 @@ return true; } -# endif Index: branches/eam_branches/ipp-20120601/psphot/src/psphotMaskReadout.c =================================================================== --- branches/eam_branches/ipp-20120601/psphot/src/psphotMaskReadout.c (revision 34051) +++ branches/eam_branches/ipp-20120601/psphot/src/psphotMaskReadout.c (revision 34053) @@ -94,5 +94,5 @@ // test output of files at this stage - if (psTraceGetLevel("psphot") >= 5) { + if (psTraceGetLevel("psphot.imsave") >= 5) { psphotSaveImage (NULL, readout->image, "image.fits"); psphotSaveImage (NULL, readout->mask, "mask.fits"); @@ -105,2 +105,101 @@ return true; } + +// XXX this function and support below was created to test the theory that the faint-end +// bias results from the Poisson variation of the background pixels. This is NOT the +// case. Using the code below maintains the faint-end bias. +bool psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule) { + + bool status = false; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + + // Generate the mask and weight images, including the user-defined analysis region of interest + if (!psphotUpdateVarianceReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +// determine the mean variance image (equivalent to the background model, but for the variance image) +// set the variance image to the MAX(input, mean) +bool psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { + + bool status; + + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + // find the currently selected readout + pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); + psAssert (readout, "missing readout?"); + + pmSourceFitVarMode varMode = psphotGetFitVarMode (recipe); + if (varMode == PM_SOURCE_PHOTFIT_NONE) { + psError (PSPHOT_ERR_CONFIG, false, "need valid LINEAR_FIT_VARIANCE_MODE"); + return false; + } + + // make this an option via the recipe + if (varMode != PM_SOURCE_PHOTFIT_MODEL_SKY) return true; + + // create a model variance image (full-scale image to take result of psImageUnbin below) + psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32); + + // find the binning information + 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, varModel, backBinning)) { + psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning"); + psFree (varModel); + psFree (varModelStdev); + return false; + } + + // XXX save these? + psFree (varModel); + psFree (varModelStdev); + + psImage *im = readout->image; + psImage *wt = readout->variance; + for (int j = 0; j < im->numRows; j++) { + for (int i = 0; i < im->numCols; i++) { + if (!isfinite(im->data.F32[j][i])) continue; + if (!isfinite(wt->data.F32[j][i])) continue; + // XXX for a test, make variance constant wt->data.F32[j][i] = PS_MAX(wt->data.F32[j][i], modelVar->data.F32[j][i]); + wt->data.F32[j][i] = modelVar->data.F32[j][i]; + } + } + + // test output of files at this stage + if (psTraceGetLevel("psphot.imsave") >= 5) { + psphotSaveImage (NULL, readout->image, "image.varsky.fits"); + psphotSaveImage (NULL, readout->mask, "mask.varsky.fits"); + psphotSaveImage (NULL, readout->variance, "variance.varsky.fits"); + } + + psFree (modelVar); + + return true; +}