Index: /branches/eam_branches/ipp-20100621/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphot.h (revision 28781) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphot.h (revision 28782) @@ -180,11 +180,11 @@ // functions to set the correct source pixels -bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type); +bool psphotInitRadiusPSF (psMetadata *recipe, pmReadout *readout); +bool psphotInitRadiusEXT (psMetadata *recipe, pmReadout *readout); bool psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal); bool psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal, float dR); -bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type); -bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal); -float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal); +bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor); +bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep); bool psphotDumpMoments (psMetadata *recipe, psArray *sources); @@ -202,5 +202,5 @@ // functions to support the source fitting process bool psphotInitLimitsPSF (psMetadata *recipe, pmReadout *readout); -bool psphotInitLimitsEXT (psMetadata *recipe); +bool psphotInitLimitsEXT (psMetadata *recipe, pmReadout *readout); bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); @@ -426,5 +426,7 @@ bool psphotStackObjectsUnifyPosition (psArray *objects); -bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); +bool psphotFitSersicIndex (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); + +bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmReadout *readout, 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/psphotAddNoise.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotAddNoise.c (revision 28781) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotAddNoise.c (revision 28782) @@ -34,7 +34,4 @@ bool status = false; - psEllipseShape oldshape; - psEllipseShape newshape; - psEllipseAxes axes; // find the currently selected readout @@ -78,4 +75,6 @@ } + psphotSaveImage (NULL, readout->variance, "test.var0.fits"); + // loop over all source for (int i = 0; i < sources->n; i++) { @@ -86,44 +85,5 @@ if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; - // select appropriate model - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) continue; // model must be defined - - if (add) { - psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); - } else { - psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); - } - - psF32 *PAR = model->params->data.F32; - - // save original values - float oldI0 = PAR[PM_PAR_I0]; - oldshape.sx = PAR[PM_PAR_SXX]; - oldshape.sy = PAR[PM_PAR_SYY]; - oldshape.sxy = PAR[PM_PAR_SXY]; - - // XXX can this be done more intelligently? - if (oldI0 == 0.0) continue; - if (!isfinite(oldI0)) continue; - - // increase size and height of source - axes = psEllipseShapeToAxes (oldshape, 20.0); - axes.major *= SIZE; - axes.minor *= SIZE; - newshape = psEllipseAxesToShape (axes); - PAR[PM_PAR_I0] = FACTOR*oldI0; - PAR[PM_PAR_SXX] = newshape.sx; - PAR[PM_PAR_SYY] = newshape.sy; - PAR[PM_PAR_SXY] = newshape.sxy; - - // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected - pmSourceOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, add, maskVal, 0, 0); - - // restore original values - PAR[PM_PAR_I0] = oldI0; - PAR[PM_PAR_SXX] = oldshape.sx; - PAR[PM_PAR_SYY] = oldshape.sy; - PAR[PM_PAR_SXY] = oldshape.sxy; + pmSourceNoiseOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, FACTOR, SIZE, add, maskVal, 0, 0); } if (add) { @@ -132,4 +92,6 @@ psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); } + + psphotSaveImage (NULL, readout->variance, "test.var1.fits"); return true; } Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotBlendFit.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotBlendFit.c (revision 28781) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotBlendFit.c (revision 28782) @@ -90,6 +90,6 @@ psphotInitLimitsPSF (recipe, readout); - psphotInitLimitsEXT (recipe); - psphotInitRadiusPSF (recipe, readout->analysis, psf->type); + psphotInitLimitsEXT (recipe, readout); + psphotInitRadiusPSF (recipe, readout); // starts the timer, sets up the array of fitSets Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotExtendedSourceFits.c (revision 28781) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotExtendedSourceFits.c (revision 28782) @@ -39,4 +39,5 @@ int NplainPass = 0; bool savePics = false; + float radius; // find the currently selected readout @@ -164,4 +165,19 @@ Next ++; + // set the radius based on the footprint (also sets the mask pixels) + if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false; + + // XXX note that this changes the source moments that are published... + // recalculate the source moments using the larger extended-source moments radius + // at this stage, skip Gaussian windowing, and do not clip pixels by S/N + // this uses the footprint to judge both radius and aperture? + // XXX save the psf-based moments for output + if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) { + // subtract the best fit from the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + // XXX raise an error of some kind + continue; + } + // save the modelFlux here in case we need to subtract it (for failure) psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux); @@ -242,4 +258,5 @@ // test for fit quality / result + modelFit->fitRadius = radius; psArrayAdd (source->modelFits, 4, modelFit); Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotGuessModels.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotGuessModels.c (revision 28781) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotGuessModels.c (revision 28782) @@ -80,5 +80,5 @@ // setup the PSF fit radius details - psphotInitRadiusPSF (recipe, readout->analysis, psf->type); + psphotInitRadiusPSF (recipe, readout); // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotRadiusChecks.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotRadiusChecks.c (revision 28781) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotRadiusChecks.c (revision 28782) @@ -8,5 +8,5 @@ // and a per-object radius is calculated) -bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) { +bool psphotInitRadiusPSF(psMetadata *recipe, pmReadout *readout) { bool status = true; @@ -15,10 +15,10 @@ PSF_FIT_PADDING = psMetadataLookupF32(&status, recipe, "PSF_FIT_PADDING"); - PSF_FIT_RADIUS = psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS"); + PSF_FIT_RADIUS = psMetadataLookupF32(&status, readout->analysis, "PSF_FIT_RADIUS"); if (!status) { PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); } - PSF_APERTURE = psMetadataLookupF32(&status, analysis, "PSF_APERTURE"); + PSF_APERTURE = psMetadataLookupF32(&status, readout->analysis, "PSF_APERTURE"); if (!status) { PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); @@ -28,5 +28,5 @@ if (PSF_FIT_RADIUS == 0.0) { - float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA"); + float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA"); if (!status) { gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); @@ -37,5 +37,5 @@ if (PSF_APERTURE == 0.0) { - float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA"); + float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA"); if (!status) { gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); @@ -122,9 +122,10 @@ } +static float EXT_FIT_SKY_SIG; static float EXT_FIT_NSIGMA; static float EXT_FIT_PADDING; static float EXT_FIT_MAX_RADIUS; -bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type) { +bool psphotInitRadiusEXT (psMetadata *recipe, pmReadout *readout) { bool status; @@ -134,9 +135,16 @@ EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS"); + float skyMean = psMetadataLookupF32 (&status, readout->analysis, "SKY_MEAN"); + float skyStdev = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV"); + + fprintf (stderr, "sky: %f +/- %f\n", skyMean, skyStdev); + + EXT_FIT_SKY_SIG = skyStdev; + return true; } // call this function whenever you (re)-define the EXT model -float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal) { +bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor) { psAssert (source, "source not defined??"); @@ -146,11 +154,11 @@ // set the radius based on the footprint: - if (!peak->footprint) goto escape; + if (!peak->footprint) return false; pmFootprint *footprint = peak->footprint; - if (!footprint->spans) goto escape; - if (footprint->spans->n < 1) goto escape; + if (!footprint->spans) return false; + if (footprint->spans->n < 1) return false; // find the max radius - float radius = 0.0; + float rawRadius = 0.0; for (int j = 0; j < footprint->spans->n; j++) { pmSpan *span = footprint->spans->data[j]; @@ -160,30 +168,22 @@ float dX1 = span->x1 - peak->xf; - radius = PS_MAX (radius, hypot(dY, dX0)); - radius = PS_MAX (radius, hypot(dY, dX1)); - } - - radius += EXT_FIT_PADDING; - if (isnan(radius)) psAbort("error in radius"); - - radius = PS_MIN (radius, EXT_FIT_MAX_RADIUS); + rawRadius = PS_MAX (rawRadius, hypot(dY, dX0)); + rawRadius = PS_MAX (rawRadius, hypot(dY, dX1)); + } + if (isnan(rawRadius)) return false; + rawRadius = PS_MIN (factor*rawRadius + EXT_FIT_PADDING, EXT_FIT_MAX_RADIUS); // redefine the pixels if needed - pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, radius); - - // set the mask to flag the excluded pixels - psImageKeepCircle (source->maskObj, peak->xf, peak->yf, radius, "OR", markVal); - return radius; - -escape: - return NAN; - // bool result = psphotCheckRadiusEXT (readout, source, model, markVal); - // return result; + pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, rawRadius); + + // set the mask to flag the excluded pixels + psImageKeepCircle (source->maskObj, peak->xf, peak->yf, rawRadius, "OR", markVal); + + *radius = rawRadius; + return true; } // alternative EXT radius based on model guess (for use without footprints) -bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal) { - - psAbort ("do not use this function"); +bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) { psF32 *PAR = model->params->data.F32; @@ -193,14 +193,17 @@ // set the fit radius based on the object flux limit and the model - float rawRadius = model->modelRadius (model->params, EXT_FIT_NSIGMA*moments->dSky); - - model->fitRadius = rawRadius + EXT_FIT_PADDING; - if (isnan(model->fitRadius)) psAbort("error in radius"); + float flux = deep ? EXT_FIT_NSIGMA*EXT_FIT_SKY_SIG : 0.1 * model->params->data.F32[PM_PAR_I0]; + + float rawRadius = model->modelRadius (model->params, flux); + if (isnan(rawRadius)) return false; + + rawRadius = PS_MIN (rawRadius + EXT_FIT_PADDING, EXT_FIT_MAX_RADIUS); + model->fitRadius = rawRadius; // redefine the pixels if needed - bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius); + pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius); // set the mask to flag the excluded pixels psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal); - return status; -} + return true; +} Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c (revision 28781) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c (revision 28782) @@ -1,4 +1,3 @@ # include "psphotInternal.h" -bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal); // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful @@ -19,5 +18,5 @@ psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n", - NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits")); + NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits")); return true; } @@ -205,7 +204,8 @@ } +// save a local, static copy of the EXT model type so we don't have to lookup for each object static pmModelType modelTypeEXT; -bool psphotInitLimitsEXT (psMetadata *recipe) { +bool psphotInitLimitsEXT (psMetadata *recipe, pmReadout *readout) { bool status; @@ -214,5 +214,6 @@ char *modelNameEXT = psMetadataLookupStr (&status, recipe, "EXT_MODEL"); modelTypeEXT = pmModelClassGetType (modelNameEXT); - psphotInitRadiusEXT (recipe, modelTypeEXT); + + psphotInitRadiusEXT (recipe, readout); return true; @@ -221,4 +222,5 @@ bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { + float radius; bool okEXT, okDBL; float chiEXT, chiDBL; @@ -230,4 +232,5 @@ // skip the source if we don't think it is extended + // XXX are these robust, or do we have better info from the source-size analysis?? if (source->type == PM_SOURCE_TYPE_UNKNOWN) return false; if (source->type == PM_SOURCE_TYPE_DEFECT) return false; @@ -235,10 +238,12 @@ // set the radius based on the footprint (also sets the mask pixels) - float radius = psphotSetRadiusEXT (readout, source, markVal); + if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false; // XXX note that this changes the source moments that are published... + // XXX all published moments should use the same measurement // recalculate the source moments using the larger extended-source moments radius // at this stage, skip Gaussian windowing, and do not clip pixels by S/N // this uses the footprint to judge both radius and aperture? + // XXX save the psf-based moments for output if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false; @@ -250,42 +255,47 @@ // this temporary source is used as a place-holder by the psphotEval functions below tmpSrc = pmSourceAlloc (); - - // XXX need to handle failures better here - EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal); - if (!EXT) goto escape; - if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape; - - okEXT = psphotEvalEXT (tmpSrc, EXT); - chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN; - - // DBL will always be defined, but DBL->data[n] might not - DBL = psphotFitDBL (readout, source, fitOptions, maskVal, markVal); - if (!DBL) goto escape; - if (!DBL->n) goto escape; - - okDBL = psphotEvalDBL (tmpSrc, DBL->data[0]); - okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]); - // XXX should I keep / save the flags set in the eval functions? + { + // DBL will always be defined, but DBL->data[n] might not + DBL = psphotFitDBL (readout, source, fitOptions, maskVal, markVal); + if (!DBL) goto escape; + if (!DBL->n) goto escape; + + okDBL = psphotEvalDBL (tmpSrc, DBL->data[0]); + okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]); + // XXX should I keep / save the flags set in the eval functions? + + // correct first model chisqs for flux trend + chiDBL = NAN; + ONE = DBL->data[0]; + if (ONE) { + if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); + chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); + ONE->chisqNorm = ONE->chisq / chiTrend; + chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison + ONE->fitRadius = radius; + } + + // correct second model chisqs for flux trend + ONE = DBL->data[1]; + if (ONE) { + if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); + chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); + ONE->chisqNorm = ONE->chisq / chiTrend; + ONE->fitRadius = radius; + } + } + + { + // XXX need to handle failures better here + EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal); + if (!EXT) goto escape; + if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape; + + okEXT = psphotEvalEXT (tmpSrc, EXT); + chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN; + } // clear the circular mask psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // correct first model chisqs for flux trend - chiDBL = NAN; - ONE = DBL->data[0]; - if (ONE) { - if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); - ONE->chisqNorm = ONE->chisq / chiTrend; - chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison - } - - // correct second model chisqs for flux trend - ONE = DBL->data[1]; - if (ONE) { - if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); - ONE->chisqNorm = ONE->chisq / chiTrend; - } psFree (tmpSrc); @@ -314,6 +324,6 @@ // save new model + // XXX save the correct radius... source->modelEXT = EXT; - source->modelEXT->fitRadius = radius; source->type = PM_SOURCE_TYPE_EXTENDED; source->mode |= PM_SOURCE_MODE_EXTMODEL; @@ -343,10 +353,8 @@ source->modelPSF = psMemIncrRefCounter (DBL->data[0]); source->mode |= PM_SOURCE_MODE_PAIR; - source->modelPSF->fitRadius = radius; // copy most data from the primary source (modelEXT, blends stay NULL) pmSource *newSrc = pmSourceCopyData (source); newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]); - newSrc->modelPSF->fitRadius = radius; // build cached models and subtract @@ -452,9 +460,7 @@ maskVal |= markVal; - // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); - // use the source moments, etc to guess basic model parameters - pmModel *EXT = pmSourceModelGuess (source, modelType); - if (!EXT) { + pmModel *model = pmSourceModelGuess (source, modelType); + if (!model) { psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); return NULL; @@ -463,5 +469,11 @@ // for sersic models, use a grid search to choose an index, then float the params there if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { - psphotFitSersicIndex (source, EXT, fitOptions, maskVal); + // for the test fits, use a somewhat smaller radius + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5); + psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal); + } + + if (!psphotSetRadiusModel (model, readout, source, markVal, true)) { + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0); } @@ -471,11 +483,25 @@ options.mode = PM_SOURCE_FIT_EXT; } - pmSourceFitModel (source, EXT, &options, maskVal); + + // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); + pmSourceFitModel (source, model, &options, maskVal); + fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix); // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); - return (EXT); + return (model); } pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + 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); + } + + pmSourceFitOptions options = *fitOptions; + + NfitPCM ++; + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; // allocate the model @@ -484,29 +510,4 @@ 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); - } - - float radius = psphotSetRadiusEXT (readout, source, markVal); - - // XXX note that this changes the source moments that are published... - // recalculate the source moments using the larger extended-source moments radius - // 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)) { - // XXX set some mask bit/ - model->flags |= PM_MODEL_STATUS_BADARGS; - return model; - } - - NfitPCM ++; - - // maskVal is used to test for rejected pixels, and must include markVal - maskVal |= markVal; - - // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize); @@ -518,9 +519,22 @@ // use the source moments, etc to guess basic model parameters - pmSourceModelGuessPCM (pcm, source, maskVal, markVal); + if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) { + psFree (pcm); + model->flags |= PM_MODEL_STATUS_BADARGS; + return model; + } // 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, psfSize); + // for the test fits, use a somewhat smaller radius + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5); + if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) { + model->flags |= PM_MODEL_STATUS_BADARGS; + return model; + } + } + + if (!psphotSetRadiusModel (model, readout, source, markVal, true)) { + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0); } @@ -530,5 +544,8 @@ options.mode = PM_SOURCE_FIT_EXT; } + + // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix); // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); @@ -544,5 +561,5 @@ // 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 (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) { +bool psphotFitSersicIndex (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); @@ -552,5 +569,5 @@ // fit EXT (not PSF) model (set/unset the pixel mask) options.mode = PM_SOURCE_FIT_NO_INDEX; - options.nIter = 3; + options.nIter = 4; int iMin = -1; @@ -561,6 +578,80 @@ model->params->data.F32[PM_PAR_7] = indexGuess[i]; - model->modelGuess(model, source); + if (!model->modelGuess(model, source)) { + model->flags |= PM_MODEL_STATUS_BADARGS; + return false; + } + + // each time we change the model guess, we need to adjust the radius + // XXX this did not work : we do not need such a large radius -- just uses moments-based radius + if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) { + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5); + } + pmSourceFitModel (source, model, &options, maskVal); + fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix); + + chiSquare[i] = model->chisqNorm; + if (i == 0) { + xMin = chiSquare[i]; + iMin = i; + } else { + if (chiSquare[i] < xMin) { + xMin = chiSquare[i]; + iMin = i; + } + } + } + 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); + + // each time we change the model guess, we need to adjust the radius + // if (!psphotSetRadiusModel (model, readout, source, markVal, true)) { + // psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal); + // } + + return true; +} + +// 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 (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + + // fit EXT (not PSF) model (set/unset the pixel mask) + options.mode = PM_SOURCE_FIT_NO_INDEX; + options.nIter = 4; + + 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]; + + if (!model->modelGuess(model, source)) { + model->flags |= PM_MODEL_STATUS_BADARGS; + return false; + } + + // each time we change the model guess, we need to adjust the radius + // XXX this did not work : we do not need such a large radius -- just uses moments-based radius + if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) { + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5); + } + + pmSourceFitModel (source, model, &options, maskVal); + fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix); + + // pmSourceModelGuessPCM(pcm, source, maskVal, markVal); + // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); chiSquare[i] = model->chisq; @@ -575,51 +666,4 @@ } } - - 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; -} - -// 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 (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { - - pmModel *model = pcm->modelConv; - - assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); - - pmSourceFitOptions options = *fitOptions; - - // fit EXT (not PSF) model (set/unset the pixel mask) - options.mode = PM_SOURCE_FIT_NO_INDEX; - options.nIter = 3; - - 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); - pmSourceFitModel (source, model, &options, maskVal); - - // pmSourceModelGuessPCM(pcm, source, maskVal, markVal); - // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); - - chiSquare[i] = model->chisq; - if (i == 0) { - xMin = chiSquare[i]; - iMin = i; - } else { - if (chiSquare[i] < xMin) { - xMin = chiSquare[i]; - iMin = i; - } - } - } - assert (iMin >= 0);