Index: /branches/eam_branches/ipp-20100621/ippconfig/gpc1/ppMerge.config =================================================================== --- /branches/eam_branches/ipp-20100621/ippconfig/gpc1/ppMerge.config (revision 28701) +++ /branches/eam_branches/ipp-20100621/ippconfig/gpc1/ppMerge.config (revision 28702) @@ -60,4 +60,15 @@ # END +# DARK.ORDINATES METADATA +# TEMP METADATA +# ORDER S32 1 +# RULE STR CHIP.TEMP +# END +# TEMP2 METADATA +# ORDER S32 1 +# RULE STR CHIP.TEMP * CHIP.TEMP +# END +# END + # Ordinates for fitting dark current as function of darktime and fpa temp: # Counts = C0 + (CT0 + CT1*Temp)*exptime Index: /branches/eam_branches/ipp-20100621/ippconfig/recipes/psphot.config =================================================================== --- /branches/eam_branches/ipp-20100621/ippconfig/recipes/psphot.config (revision 28701) +++ /branches/eam_branches/ipp-20100621/ippconfig/recipes/psphot.config (revision 28702) @@ -93,6 +93,6 @@ PSF_FIT_MAX_TOL F32 2.00 # Fit tolerance for PSF -PSF_FIT_MIN_VALID_FLUX F32 1e-8 # minimum allow flux for fitted source -PSF_FIT_MAX_VALID_FLUX F32 1e+8 # maximum allow flux for fitted source +PSF_FIT_MIN_VALID_FLUX F32 -100000000.0 # minimum allow flux for fitted source +PSF_FIT_MAX_VALID_FLUX F32 +100000000.0 # maximum allow flux for fitted source # the following is used to require a minimum quality of fit before Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 28701) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 28702) @@ -291,4 +291,6 @@ // convolve model and dmodel arrays with PSF + // XXX speed this up by saving the FFTed psf (for each source, obviously) + psImageConvolveDirect (pcm->modelConvFlux, pcm->modelFlux, pcm->psf); for (int n = 0; n < pcm->dmodelsFlux->n; n++) { Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c (revision 28701) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c (revision 28702) @@ -141,12 +141,16 @@ } -pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, float psfSize) { - - // make sure we savep a cached copy of the psf flux +pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) { + + // make sure we save a cached copy of the psf flux pmSourceCachePSF (source, maskVal); // convert the cached cached psf model for this source to a psKernel psKernel *psf = pmPCMkernelFromPSF (source, psfSize); - if (!psf) return NULL; + if (!psf) { + // NOTE: this only happens if the source is too close to an edge + model->flags |= PM_MODEL_STATUS_BADARGS; + return NULL; + } # if (USE_DELTA_PSF) @@ -154,11 +158,4 @@ psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0; # endif - - // allocate the model - pmModel *modelConv = pmModelAlloc(modelType); - if (!modelConv) { - psFree (psf); - return NULL; - } // count the number of unmasked pixels: @@ -183,10 +180,10 @@ } - psVector *params = modelConv->params; + psVector *params = model->params; // create the minimization constraints psMinConstraint *constraint = psMinConstraintAlloc(); constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK); - constraint->checkLimits = modelConv->modelLimits; + constraint->checkLimits = model->modelLimits; // set parameter mask based on fitting mode @@ -239,17 +236,19 @@ } + if (nPix < nParams + 1) { + psTrace ("psModules.objects", 4, "insufficient valid pixels\n"); + psFree (psf); + psFree (constraint); + model->flags |= PM_MODEL_STATUS_BADARGS; + return NULL; + } + // generate PCM data storage structure pmPCMdata *pcm = pmPCMdataAlloc (params, constraint->paramMask, source); - pcm->modelConv = modelConv; pcm->psf = psf; + pcm->modelConv = psMemIncrRefCounter(model); pcm->constraint = constraint; - if (nPix < nParams + 1) { - psTrace ("psModules.objects", 4, "insufficient valid pixels\n"); - pcm->modelConv->flags |= PM_MODEL_STATUS_BADARGS; - abort (); - // XXX This should not be an abort!! - } pcm->nPix = nPix; pcm->nDOF = nPix - nParams - 1; Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h (revision 28701) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h (revision 28702) @@ -61,5 +61,5 @@ pmSource *source); -pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, float psfSize); +pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize); psImage *pmPCMdataSaveImage (pmPCMdata *pcm); Index: /branches/eam_branches/ipp-20100621/psphot/src/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/Makefile.am (revision 28701) +++ /branches/eam_branches/ipp-20100621/psphot/src/Makefile.am (revision 28702) @@ -163,5 +163,4 @@ psphotOutput.c \ psphotFakeSources.c \ - psphotModelWithPSF.c \ psphotExtendedSourceAnalysis.c \ psphotExtendedSourceAnalysisByObject.c \ Index: /branches/eam_branches/ipp-20100621/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphot.h (revision 28701) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphot.h (revision 28702) @@ -426,5 +426,5 @@ bool psphotStackObjectsUnifyPosition (psArray *objects); -bool psphotFitSersicIndexPCM (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal); +bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinear.c (revision 28701) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotFitSourcesLinear.c (revision 28702) @@ -102,11 +102,11 @@ } - float MIN_VALID_FLUX = psMetadataLookupBool(&status, recipe, "PSF_FIT_MIN_VALID_FLUX"); + float MIN_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MIN_VALID_FLUX"); if (!status) { MIN_VALID_FLUX = 1e-8; } - float MAX_VALID_FLUX = psMetadataLookupBool(&status, recipe, "PSF_FIT_MAX_VALID_FLUX"); - if (!status) { - MAX_VALID_FLUX = 1e-8; + float MAX_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MAX_VALID_FLUX"); + if (!status) { + MAX_VALID_FLUX = 1e+8; } Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c (revision 28701) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c (revision 28702) @@ -479,4 +479,12 @@ pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + // allocate the model + pmModel *model = pmModelAlloc(modelType); + if (!model) { + return NULL; + } + + pmSourceFitOptions options = *fitOptions; + if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) { psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); @@ -489,7 +497,9 @@ // at this stage, skip Gaussian windowing, and do not clip pixels by S/N // this uses the footprint to judge both radius and aperture? - if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false; - - pmSourceFitOptions options = *fitOptions; + if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) { + // XXX set some mask bit/ + model->flags |= PM_MODEL_STATUS_BADARGS; + return model; + } NfitPCM ++; @@ -500,10 +510,10 @@ // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); - pmPCMdata *pcm = pmPCMinit (source, &options, modelType, maskVal, psfSize); + pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize); if (!pcm) { psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); - return NULL; - } - // XXX check for nDOF too small + model->flags |= PM_MODEL_STATUS_BADARGS; // XXX this is probably already set in pmPCMinit + return model; + } // use the source moments, etc to guess basic model parameters @@ -512,5 +522,5 @@ // for sersic models, use a grid search to choose an index, then float the params there if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { - psphotFitSersicIndexPCM (pcm, source, fitOptions, maskVal, markVal); + psphotFitSersicIndexPCM (pcm, source, fitOptions, maskVal, markVal, psfSize); } @@ -520,19 +530,19 @@ options.mode = PM_SOURCE_FIT_EXT; } - pmSourceFitPCM (source, PCM, &options, maskVal); + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); - return (PCM); + psFree (pcm); + + return model; } // note that these should be 1/2n of the standard sersic index float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083}; -define N_INDEX_GUESS 6 +# define N_INDEX_GUESS 6 // A sersic model is very sensitive to the index. attempt to find the index first by grid search in just the index // for a sersic model, attempt to fit just the index and normalization with a modest number of iterations -bool psphotFitSersicIndex (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { - - pmModel *model = pcm->modelConv; +bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) { assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); @@ -544,12 +554,14 @@ options.nIter = 3; - float xMin, chiSquare[N_INDEX_GUESS]; - int iMin; + int iMin = -1; + float xMin = NAN; + float chiSquare[N_INDEX_GUESS]; for (int i = 0; i < N_INDEX_GUESS; i++) { model->params->data.F32[PM_PAR_7] = indexGuess[i]; - pmSourceModelGuessPCM (pcm, source, maskVal, markVal); - - pmSourceFitPCM (pcm, source, &options, maskVal); + + model->modelGuess(model, source); + pmSourceFitModel (source, model, &options, maskVal); + chiSquare[i] = model->chisq; if (i == 0) { @@ -573,5 +585,7 @@ // A sersic model is very sensitive to the index. attempt to find the index first by grid search in just the index // for a sersic model, attempt to fit just the index and normalization with a modest number of iterations -bool psphotFitSersicIndexPCM (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) { +bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); @@ -583,13 +597,17 @@ options.nIter = 3; - float xMin, chiSquare[N_INDEX_GUESS]; - int iMin; - - // XXX we probably cannot be calling model->modelGuess() : this does not include the psf sigma + int iMin = -1; + float xMin = NAN; + float chiSquare[N_INDEX_GUESS]; for (int i = 0; i < N_INDEX_GUESS; i++) { model->params->data.F32[PM_PAR_7] = indexGuess[i]; + model->modelGuess(model, source); - pmSourceFitPCM (source, model, &options, maskVal); + pmSourceFitModel (source, model, &options, maskVal); + + // pmSourceModelGuessPCM(pcm, source, maskVal, markVal); + // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + chiSquare[i] = model->chisq; if (i == 0) { @@ -604,8 +622,11 @@ } + assert (iMin >= 0); + model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it model->params->data.F32[PM_PAR_7] = indexGuess[iMin]; - model->modelGuess(model, source); - - return true; -} + + pmSourceModelGuessPCM(pcm, source, maskVal, markVal); + + return true; +} Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c (revision 28701) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c (revision 28702) @@ -170,4 +170,5 @@ psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32); + psImageMaskType markVal = options->markVal; psImageMaskType maskVal = options->maskVal | options->markVal; @@ -295,4 +296,5 @@ pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + psImageMaskType markVal = options->markVal; psImageMaskType maskVal = options->maskVal | options->markVal; @@ -364,7 +366,7 @@ float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag); - fprintf (stderr, "%f %f : Mxx: %f, Myy: %f, dx: %f, dy: %f, psfMag: %f, apMag: %f, dMag: %f, errMag: %f, nSigmaMag: %f\n", + fprintf (stderr, "%f %f : Mxx: %f, Myy: %f, dx: %f, dy: %f, psfMag: %f, apMag: %f, dMag: %f, errMag: %f, nSigmaMag: %f, nSigmaMxx: %f, nSigmaMyy: %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, source->peak->xf - source->moments->Mx, source->peak->yf - source->moments->My, - source->psfMag, apMag, dMag, source->errMag, nSigmaMAG); + source->psfMag, apMag, dMag, source->errMag, nSigmaMAG, nSigmaMXX, nSigmaMYY); // partially-masked sources are more likely to be mis-measured PSFs