Index: /branches/eam_branches/ipp-20110710/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphot.h (revision 31989) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphot.h (revision 31990) @@ -469,3 +469,9 @@ psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc); +bool psphotSersicModelClassGuessPCM (pmPCMdata *pcm, pmSource *source); +void psphotSersicModelClassInit (); +void psphotSersicModelClassCleanup (); + +bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal); + #endif Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c (revision 31989) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c (revision 31990) @@ -45,4 +45,6 @@ int Nfail = 0; + psphotSersicModelClassInit(); + psTimerStart ("psphot.extended"); @@ -103,6 +105,6 @@ if (item->type != PS_DATA_METADATA) { + // XXX we could cull the bad entries or build a validated model folder psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name); - // XXX we could cull the bad entries or build a validated model folder } @@ -247,4 +249,6 @@ psFree(AnalysisRegion); + psphotSersicModelClassCleanup(); + psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); psLogMsg ("psphot", PS_LOG_INFO, " %d convolved models (%d passed)\n", Nconvolve, NconvolvePass); @@ -265,6 +269,5 @@ int Nplain = 0; int NplainPass = 0; - bool savePics = false; - float radius; + float fitRadius, windowRadius; psScalar *scalar = NULL; pmMoments psfMoments; @@ -279,6 +282,4 @@ psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); - // pthread_t tid = pthread_self(); // Thread identifier - // Define source fitting parameters for extended source fits pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); @@ -309,10 +310,4 @@ if (source->peak->y > region->y1) continue; - // if model is NULL, we don't have a starting guess - // XXX use the parameters guessed from moments - // if (source->modelEXT == NULL) continue; - - // fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid); - // replace object in image if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { @@ -321,19 +316,17 @@ Next ++; - // set the radius based on the footprint (also sets the mask pixels) - if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) { + // set the radius based on the first radial moment (also sets the mask pixels) + if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) { fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - // XXX raise an error of some kind + // XXX raise an error of some kind or set a flag bit continue; } - // 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 + // 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. Save the + // psf-based moments for output psfMoments = *source->moments; - if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) { + if (!pmSourceMoments (source, windowRadius, 0.0, 0.0, 0.0, maskVal & PS_NOT_IMAGE_MASK(markVal))) { // subtract the best fit from the object, leave local sky fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf); @@ -344,5 +337,6 @@ } - // save the modelFlux here in case we need to subtract it (for failure) + // save the PSF-based modelFlux here in case we need to subtract it (for failure). If + // it does not exist, attempt to generate it. Give up if we cannot even do that. psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux); if (!modelFluxStart) { @@ -358,8 +352,4 @@ } - if (savePics) { - psphotSaveImage (NULL, readout->image, "image.xp.fits"); - } - // array to store the pointers to the model flux images while the models are being fitted psArray *modelFluxes = psArrayAllocEmpty (models->list->n); @@ -368,4 +358,8 @@ if (source->modelFits == NULL) { source->modelFits = psArrayAllocEmpty (models->list->n); + } + // extFitPars are the non-parametric data measured for this model fit (moments, kron, etc) + if (source->extFitPars == NULL) { + source->extFitPars = psArrayAllocEmpty (models->list->n); } @@ -397,10 +391,4 @@ bool convolved = psMetadataLookupBool (&status, model, "PSF_CONVOLVED_VALUE"); assert (status); - - // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6); - // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6); - - // XXX this does not make sense in a threaded context - psTimerStart ("psphot.extended.fit"); // fit the model as convolved or not @@ -439,11 +427,23 @@ } - psLogMsg ("psphot", PS_LOG_INFO, "model fit: %7.5f sec, %5d npix, %2d iter, %s type\n", psTimerMark ("psphot.extended.fit"), modelFit->nPix, modelFit->nIter, pmModelClassGetName (modelFit->type)); - - // save each of the model flux images and store the best + // save each of the model flux images for now, and (below) store the best psArrayAdd (modelFluxes, 4, source->modelFlux); + pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc(); + extFitPars->Mxx = source->moments->Mxx; + extFitPars->Mxy = source->moments->Mxy; + extFitPars->Myy = source->moments->Myy; + extFitPars->Mrf = source->moments->Mrf; + extFitPars->Mrh = source->moments->Mrh; + extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN; + + // save kron mag, but assign apMag & psfMag on output (not yet calculated) + extFitPars->krMag = -2.5*log10(source->moments->KronFlux); + + psArrayAdd (source->extFitPars, 4, extFitPars); + psFree (extFitPars); + // test for fit quality / result - modelFit->fitRadius = radius; + modelFit->fitRadius = fitRadius; psArrayAdd (source->modelFits, 4, modelFit); @@ -470,6 +470,9 @@ } + // clear the circular mask + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + if (minModel == -1) { - // re-subtract the object, leave local sky + // no valid extended fit; re-subtract the object, leave local sky psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->Mx, source->moments->My); @@ -483,14 +486,4 @@ psFree (modelFluxes); - if (savePics) { - psphotSaveImage (NULL, readout->image, "image.xp.fits"); - char key[10]; - fprintf (stdout, "continue? "); - if (!fgets (key, 8, stdin)) { - psWarning("Couldn't read anything."); - } else if (key[0] == 'n') { - savePics = false; - } - } continue; } @@ -516,19 +509,6 @@ psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq); psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - - if (savePics) { - psphotSaveImage (NULL, readout->image, "image.xm.fits"); - char key[10]; - fprintf (stdout, "continue? "); - if (!fgets (key, 8, stdin)) { - psWarning("Couldn't read anything."); - } else if (key[0] == 'n') { - savePics = false; - } - } } psFree (fitOptions); - - // fprintf (stderr, "xfit : tried %ld objects\n", sources->n); // change the value of a scalar on the array (wrap this and put it in psArray.h) Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c (revision 31989) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c (revision 31990) @@ -142,4 +142,29 @@ } +# define MIN_WINDOW 5.0 +# define SCALE1 5.0 +# define SCALE2 12.5 + +// call this function whenever you (re)-define the EXT model +bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal) { + + psAssert (source, "source not defined??"); + psAssert (source->moments, "moments not defined??"); + + *fitRadius = SCALE1 * source->moments->Mrf; + *fitRadius = PS_MIN (PS_MAX(*fitRadius, MIN_WINDOW), EXT_FIT_MAX_RADIUS); + + *windowRadius = SCALE2 * source->moments->Mrf; + *windowRadius = PS_MIN (PS_MAX(*windowRadius, 2.5*MIN_WINDOW), 2.5*EXT_FIT_MAX_RADIUS); + + // redefine the pixels if needed + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, *windowRadius); + + // set the mask to flag the excluded pixels + psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal); + + return true; +} + // call this function whenever you (re)-define the EXT model bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor) { Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c (revision 31989) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c (revision 31990) @@ -505,7 +505,10 @@ } -void psphotWriteGuess(pmModel *model, float dKronMag, float PSFratio, float ASPratio); +# define TIMING false bool psphotSersicModelClassGuessPCM (pmPCMdata *pcm, pmSource *source) { + + float t1, t2, t4, t5; + if (TIMING) { psTimerStart ("SersicGuess"); } psF32 *PAR = pcm->modelConv->params->data.F32; @@ -549,4 +552,6 @@ // find the containing model class: + + if (TIMING) { t1 = psTimerMark ("SersicGuess"); } psphotSersicModelClass *class = NULL; @@ -561,4 +566,6 @@ } + if (TIMING) { t2 = psTimerMark ("SersicGuess"); } + psAssert (class, "PSFratio and ASPratio must be in range"); @@ -610,6 +617,14 @@ PAR[PM_PAR_I0] = 1.0; + if (TIMING) { t4 = psTimerMark ("SersicGuess"); } + // set the normalization by linear fit between model and data psphotSersicModelNorm (pcm, source); + + if (TIMING) { t5 = psTimerMark ("SersicGuess"); } + + if (TIMING) { + fprintf (stderr, "guess, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", t1, t2, t4, t5); + } // float flux = pcm->modelConv->modelFlux(pcm->modelConv->params); @@ -617,6 +632,4 @@ // float Io = pow(10.0, -0.4*(totalMag - normMag)); // PAR[PM_PAR_I0] = Io; - - psphotWriteGuess(pcm->modelConv, dKronMag, PSFratio, ASPratio); return true; Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c (revision 31989) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c (revision 31990) @@ -211,5 +211,5 @@ bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { - float radius; + float fitRadius, windowRadius; bool okEXT, okDBL; pmModel *ONE = NULL; @@ -225,6 +225,17 @@ if (source->type == PM_SOURCE_TYPE_SATURATED) return false; +# define TEST_X -540.0 +# define TEST_Y 540.0 + + if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { + fprintf (stderr, "test galaxy\n"); + } + +# undef TEST_X +# undef TEST_Y + // set the radius based on the footprint (also sets the mask pixels) - if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false; + if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) return false; + // fprintf (stderr, "rad: %6.1f %6.1f | %5.2f %5.2f %5.2f ", source->peak->xf, source->peak->yf, source->moments->Mrf, fitRadius, windowRadius); // XXX note that this changes the source moments that are published... @@ -234,9 +245,15 @@ // this uses the footprint to judge both radius and aperture? // XXX save the psf-based moments for output + // XXX do not limit to pixels in the marked circle psfMoments = *source->moments; - if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) { + if (!pmSourceMoments (source, windowRadius, 0.0, 0.5, 0.0, maskVal & PS_NOT_IMAGE_MASK(markVal))) { *source->moments = psfMoments; + fprintf (stderr, "skip (bad moments)\n"); return false; } + // fprintf (stderr, "r1: %f\n", source->moments->Mrf); + + // XXX note that we do not re-try to measure the moments if the resulting Mrf is large compared + // to the initial value psTrace ("psphot", 5, "trying blob...\n"); @@ -263,7 +280,7 @@ ONE = DBL->data[0]; if (ONE) { - if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); + psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit"); chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison - ONE->fitRadius = radius; + ONE->fitRadius = fitRadius; } @@ -271,6 +288,6 @@ ONE = DBL->data[1]; if (ONE) { - if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - ONE->fitRadius = radius; + psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit"); + ONE->fitRadius = fitRadius; } } @@ -284,4 +301,5 @@ okEXT = psphotEvalEXT (tmpSrc, EXT); chiEXT = EXT ? EXT->chisqNorm : NAN; + EXT->fitRadius = fitRadius; } @@ -476,11 +494,5 @@ // for sersic models, use a grid search to choose an index, then float the params there if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { - // 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); } @@ -499,4 +511,6 @@ } +# define TIMING 0 + pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { @@ -517,4 +531,14 @@ return NULL; } + +# define TEST_X -540.0 +# define TEST_Y 540.0 + + if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { + psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5); + } + + float t1, t2, t4, t5; + if (TIMING) { psTimerStart ("psphotFitPCM"); } pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize); @@ -524,27 +548,24 @@ return model; } - - // use the source moments, etc to guess basic model parameters - 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 (TIMING) { t1 = psTimerMark ("psphotFitPCM"); } + + // get the guess for sersic models if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { - // 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)) { + // use the source moments, etc to guess basic model parameters + if (!psphotSersicModelClassGuessPCM (pcm, source)) { psFree (pcm); model->flags |= PM_MODEL_STATUS_BADARGS; return model; } - } - - if (!psphotSetRadiusModel (model, readout, source, markVal, true)) { - psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0); - } + } else { + // use the source moments, etc to guess basic model parameters + if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) { + psFree (pcm); + model->flags |= PM_MODEL_STATUS_BADARGS; + return model; + } + } + + if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); } if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { @@ -555,8 +576,17 @@ // update the pcm elements if we have changed the circumstance (options.mode or source->pixels) pmPCMupdate(pcm, source, &options, model); + if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); } // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); - // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix); + if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); } + + if (TIMING) { + fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, t1, t2, t4, t5); + } + + if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { + psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0); + } // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); @@ -566,7 +596,12 @@ } +# undef TEST_X +# undef TEST_Y + // 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 +// float indexGuess[] = {0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0}; +// float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083}; +float indexGuess[] = {1.0, 2.0, 3.0, 4.0}; +# define N_INDEX_GUESS 4 // A sersic model is very sensitive to the index. attempt to find the index first by grid search in just the index @@ -586,6 +621,13 @@ float chiSquare[N_INDEX_GUESS]; +# define TEST_X -540.0 +# define TEST_Y 540.0 + + if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { + psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5); + } + for (int i = 0; i < N_INDEX_GUESS; i++) { - model->params->data.F32[PM_PAR_7] = indexGuess[i]; + model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i]; if (!model->modelGuess(model, source)) { @@ -594,12 +636,6 @@ } - // 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); + fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix); chiSquare[i] = model->chisqNorm; @@ -616,12 +652,11 @@ assert (iMin >= 0); + if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { + psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 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->params->data.F32[PM_PAR_7] = 0.5/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; @@ -648,4 +683,8 @@ float xMin = NAN; float chiSquare[N_INDEX_GUESS]; + + if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { + psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5); + } for (int i = 0; i < N_INDEX_GUESS; i++) { @@ -657,11 +696,6 @@ } - // 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); - } - # if (0) + // this block is to test the relative speed of straight and PCM fits pmSourceFitModel (source, model, &options, maskVal); # else @@ -669,4 +703,5 @@ pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); # endif + fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix); // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix); @@ -684,4 +719,8 @@ assert (iMin >= 0); + if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { + psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 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];