Index: /branches/eam_branches/ipp-20110710/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphot.h (revision 32156) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphot.h (revision 32157) @@ -19,4 +19,5 @@ PSPHOT_FORCED, PSPHOT_MAKE_PSF, + PSPHOT_MODEL_TEST, } psphotImageLoopMode; @@ -328,4 +329,7 @@ bool psphotMakePSFReadout(pmConfig *config, const pmFPAview *view, const char *filerule); +pmConfig *psphotModelTestArguments(int argc, char **argv); +bool psphotModelTestReadout(pmConfig *config, const pmFPAview *view, const char *filerule); + int psphotFileruleCount(const pmConfig *config, const char *filerule); bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num); @@ -476,4 +480,5 @@ bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal); +bool psphotSetRadiusMomentsExact (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 32156) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c (revision 32157) @@ -316,11 +316,6 @@ Next ++; - // 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 or set a flag bit - continue; - } + // set the fit radius based on the first radial moment (also sets the mask pixels) + psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // Recalculate the source moments using the larger extended-source moments radius. At @@ -332,26 +327,10 @@ fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf); *source->moments = psfMoments; + + // probably need to bump up the size for subtraction pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); // XXX raise an error flag of some kind continue; } - - // 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) { - 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; - } - } - - // array to store the pointers to the model flux images while the models are being fitted - psArray *modelFluxes = psArrayAllocEmpty (models->list->n); // allocate the array to store the model fits @@ -417,5 +396,4 @@ continue; } - pmSourceCacheModel (source, maskVal); 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); @@ -426,7 +404,4 @@ } } - - // save each of the model flux images for now, and (below) store the best - psArrayAdd (modelFluxes, 4, source->modelFlux); pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc(); @@ -477,12 +452,9 @@ 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; + // ensure the modelEXT is cached + pmSourceCacheModel (source, maskVal); *source->moments = psfMoments; pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - - psFree (modelFluxes); continue; @@ -493,7 +465,13 @@ source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]); - // save the modelFlux for the best model - psFree (source->modelFlux); - source->modelFlux = psMemIncrRefCounter (modelFluxes->data[minModel]); + // 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); + } // replace the original moments @@ -502,8 +480,4 @@ // 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); Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotKronMasked.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotKronMasked.c (revision 32156) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotKronMasked.c (revision 32157) @@ -96,4 +96,14 @@ if (!source->moments) continue; + float windowRadius = RADIUS; + if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) { + windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf); + } + + // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS + pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); + + pmMoments psfMoments = *source->moments; + // replace object in image if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { @@ -101,5 +111,5 @@ } - psphotKronMag (source, RADIUS, MIN_KRON_RADIUS, maskVal); + psphotKronMag (source, windowRadius, MIN_KRON_RADIUS, maskVal); // re-subtract the object, leave local sky @@ -136,4 +146,5 @@ } } + // psphotSaveImage (NULL, kronMask, "kronmask.fits"); psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c (revision 32156) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c (revision 32157) @@ -144,7 +144,8 @@ # define MIN_WINDOW 5.0 # define SCALE1 5.0 -# define SCALE2 12.5 +# define SCALE2 12.0 // call this function whenever you (re)-define the EXT model +// XXX this function does not shrink the window bool psphotSetRadiusMoments (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal) { @@ -160,4 +161,38 @@ // 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; +} +# undef SCALE1 +# undef SCALE2 +# undef MIN_WINDOW + +# define MIN_WINDOW 5.0 +# define SCALE1 5.0 +# define PAD_WINDOW 3.0 + +// call this function whenever you (re)-define the EXT model +// XXX alternate function to set exactly the desired window size +bool psphotSetRadiusMomentsExact (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal) { + + psRegion newRegion; + + 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 = *fitRadius + PAD_WINDOW; + + // check to see if new region is completely contained within old region + newRegion = psRegionForSquare (source->peak->xf, source->peak->yf, *windowRadius); + newRegion = psRegionForImage (readout->image, newRegion); + + // redefine the pixels to match + pmSourceRedefinePixelsByRegion (source, readout, newRegion); // set the mask to flag the excluded pixels Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c (revision 32156) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c (revision 32157) @@ -137,5 +137,7 @@ } - // construct sources and measure basic stats (saved on detections->newSources) + // construct sources and measure moments and other basic stats (saved on detections->newSources) + // all sources use the auto-scaled window appropriate to a PSF, except for the saturated + // stars : these use a larger window (3x the basic window) if (!psphotSourceStats (config, view, filerule, true)) { // pass 1 psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); @@ -200,5 +202,7 @@ psphotDumpChisqs (config, view, filerule); - // XXX re-measure the kron mags with models subtracted + // re-measure the kron mags with models subtracted. this pass uses a circular window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments) + + // but this is chosen above to be appropriate for the PSF objects (not galaxies) psphotKronMasked(config, view, filerule); @@ -214,4 +218,7 @@ psphotDumpChisqs (config, view, filerule); + // XXX note above the extended sources have their moments remeasured with a larger window, + // but the result is replaced with the original PSF-scaled window + // replace all sources psphotReplaceAllSources (config, view, filerule); // pass 1 (detections->allSources) @@ -240,5 +247,5 @@ psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources) - // define new sources based on only the new peaks + // define new sources based on only the new peaks & measure moments // NOTE: new sources are saved on detections->newSources psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources) @@ -313,4 +320,5 @@ // XXX re-measure the kron mags with models subtracted + // Note that this uses the PSF_MOMENTS_RADIUS as a window psphotKronMasked(config, view, filerule); Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c (revision 32156) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceFits.c (revision 32157) @@ -511,5 +511,5 @@ } -# define TIMING 0 +# define TIMING 1 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { @@ -583,5 +583,6 @@ 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); + int nPixBig = source->pixels->numCols * source->pixels->numRows; + fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5); }