Index: trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceFits.c (revision 31673) +++ trunk/psphot/src/psphotExtendedSourceFits.c (revision 32348) @@ -6,4 +6,7 @@ bool status = true; + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Extended Source Fits ---"); + // select the appropriate recipe information psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); @@ -12,5 +15,5 @@ // perform full extended source non-linear fits? if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) { - psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); + psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source fits\n"); return true; } @@ -45,4 +48,6 @@ int Nfail = 0; + psphotSersicModelClassInit(); + psTimerStart ("psphot.extended"); @@ -103,6 +108,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,5 +252,9 @@ psFree(AnalysisRegion); - psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); + psphotSersicModelClassCleanup(); + + psphotVisualShowResidualImage (readout, false); + + psLogMsg ("psphot", PS_LOG_WARN, "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); psLogMsg ("psphot", PS_LOG_INFO, " %d plain models (%d passed)\n", Nplain, NplainPass); @@ -265,8 +274,6 @@ int Nplain = 0; int NplainPass = 0; - bool savePics = false; - float radius; + float fitRadius, windowRadius; psScalar *scalar = NULL; - pmMoments psfMoments; // arguments: readout, sources, models, region, psfSize, maskVal, markVal @@ -275,9 +282,7 @@ psMetadata *models = job->args->data[2]; psRegion *region = job->args->data[3]; - int psfSize = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); + int psfSize = PS_SCALAR_VALUE(job->args->data[4],S32); psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); 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 @@ -299,5 +304,5 @@ // skip PSF-like and non-astronomical objects - if (source->type == PM_SOURCE_TYPE_STAR) continue; + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; @@ -309,10 +314,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,47 +320,11 @@ Next ++; - // set the radius based on the footprint (also sets the mask pixels) - if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) { - 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 - 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 - psfMoments = *source->moments; - if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) { - // subtract the best fit from the object, leave local sky - fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf); - *source->moments = psfMoments; - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - // XXX raise an error flag of some kind - continue; - } - - // save the modelFlux here in case we need to subtract it (for failure) - psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux); - if (!modelFluxStart) { - pmSourceCacheModel (source, maskVal); - modelFluxStart = psMemIncrRefCounter (source->modelFlux); - if (!modelFluxStart) { - fprintf (stderr, "skipping (3) %f, %f\n", source->peak->xf, source->peak->yf); - *source->moments = psfMoments; - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - // XXX raise an error of some kind? - continue; - } - } - - 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); + // set the fit radius based on the first radial moment (also sets the mask pixels) + psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); + + // UPDATE : we have changed the moments calculation. There is now an iteration within + // psphotKronMasked to determine moments appropriate for a larger object. The values + // Mrf, KronFlux, and KronFluxErr are calculated for the iterated radius. The other + // values are left at the psf-based values. // allocate the array to store the model fits @@ -380,10 +343,17 @@ // check the SNlim and skip model if source is too faint - float SNlim = psMetadataLookupF32 (&status, model, "SNLIM_VALUE"); + float FIT_SN_LIM = psMetadataLookupF32 (&status, model, "SNLIM_VALUE"); assert (status); // limit selection to some SN limit // assert (source->peak); // how can a source not have a peak? - if (sqrt(source->peak->detValue) < SNlim) { + // limit selection to some SN limit + bool skipSource = false; + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + skipSource = (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr); + } else { + skipSource = (sqrt(source->peak->detValue) < FIT_SN_LIM); + } + if (skipSource) { Nfaint ++; continue; @@ -397,7 +367,4 @@ bool convolved = psMetadataLookupBool (&status, model, "PSF_CONVOLVED_VALUE"); assert (status); - - // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6); - // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6); // fit the model as convolved or not @@ -410,5 +377,6 @@ continue; } - psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq); + psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", + source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter); Nconvolve ++; if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { @@ -425,6 +393,6 @@ continue; } - pmSourceCacheModel (source, maskVal); - psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq); + psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", + source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter); Nplain ++; if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { @@ -434,9 +402,21 @@ } - // save each of the model flux images and store the best - psArrayAdd (modelFluxes, 4, source->modelFlux); + // XXX deprecate? + // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc(); + // XXX extFitPars->Mxx = source->moments->Mxx; + // XXX extFitPars->Mxy = source->moments->Mxy; + // XXX extFitPars->Myy = source->moments->Myy; + // XXX extFitPars->Mrf = source->moments->Mrf; + // XXX extFitPars->Mrh = source->moments->Mrh; + // XXX 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) + // XXX 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); @@ -463,27 +443,16 @@ } + // 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); - // replace original model, subtract it - psFree (source->modelFlux); - source->modelFlux = modelFluxStart; - - *source->moments = psfMoments; + // ensure the modelEXT is cached + pmSourceCacheModel (source, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - 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; } @@ -493,35 +462,21 @@ source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]); - // save the modelFlux for the best model - psFree (source->modelFlux); - source->modelFlux = psMemIncrRefCounter (modelFluxes->data[minModel]); - - // replace the original moments - *source->moments = psfMoments; + // adjust the window so the subtraction covers the faint wings + psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal); + + // cache the model flux + if (source->modelEXT->isPCM) { + pmPCMCacheModel (source, maskVal, psfSize); + } else { + pmSourceCacheModel (source, maskVal); + } // subtract the best fit from the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - // the initial model flux is no longer needed - psFree (modelFluxStart); - psFree (modelFluxes); - 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)