Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 5771) +++ /trunk/psphot/Makefile (revision 5772) @@ -34,5 +34,4 @@ $(SRC)/psphotOutput.$(ARCH).o \ $(SRC)/psphotMarkPSF.$(ARCH).o \ -$(SRC)/psphotSubtractPSF.$(ARCH).o \ $(SRC)/psphotSortBySN.$(ARCH).o \ $(SRC)/psphotDefinePixels.$(ARCH).o\ @@ -41,4 +40,9 @@ $(SRC)/psphotBasicDeblend.$(ARCH).o \ $(SRC)/psphotModelGroupInit.$(ARCH).o \ +$(SRC)/psphotRadiusChecks.$(ARCH).o \ +$(SRC)/psphotReplaceUnfit.$(ARCH).o \ +$(SRC)/psphotEvalPSF.$(ARCH).o \ +$(SRC)/psphotEvalFLT.$(ARCH).o \ +$(SRC)/psphotFullFit.$(ARCH).o \ $(SRC)/pmSourceContour.$(ARCH).o \ $(SRC)/psLine.$(ARCH).o \ Index: /trunk/psphot/src/models/pmModel_SGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_SGAUSS.c (revision 5771) +++ /trunk/psphot/src/models/pmModel_SGAUSS.c (revision 5772) @@ -81,5 +81,5 @@ params_min[0][0].data.F32[6] = -5.0; params_min[0][0].data.F32[7] = 0.5; - params_min[0][0].data.F32[8] = 0.001; + params_min[0][0].data.F32[8] = 0.05; params_max[0][0].data.F32[0] = 1e5; Index: /trunk/psphot/src/pmPeaksSigmaLimit.c =================================================================== --- /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 5771) +++ /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 5772) @@ -24,6 +24,6 @@ NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT"); - // threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; - threshold = NSIGMA*sky->sampleStdev; + threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; + // threshold = NSIGMA*sky->sampleStdev; psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold); Index: /trunk/psphot/src/psModulesUtils.c =================================================================== --- /trunk/psphot/src/psModulesUtils.c (revision 5771) +++ /trunk/psphot/src/psModulesUtils.c (revision 5772) @@ -101,2 +101,17 @@ return true; } + +pmModel *pmModelCopy (pmModel *model) { + + pmModel *new = pmModelAlloc (model->type); + + new->chisq = model->chisq; + new->nIter = model->nIter; + + for (int i = 0; i < new->params->n; i++) { + new->params->data.F32[i] = model->params->data.F32[i]; + new->dparams->data.F32[i] = model->dparams->data.F32[i]; + } + + return (new); +} Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 5771) +++ /trunk/psphot/src/psphot.c (revision 5772) @@ -43,8 +43,17 @@ psfClump = pmSourcePSFClump (sources, config); - // group into STAR, COSMIC, GALAXY, SATURATED + // group into STAR, COSMIC, GALAXY, SATURATED, etc. pmSourceRoughClass (sources, config, psfClump); + + // optional dump of all rough source data + char *output = psMetadataLookupPtr (&status, config, "MOMENTS_OUTPUT_FILE"); + if (status && (output != NULL) && (output[0])) + { + pmMomentsWriteText (sources, output); + psFree (output); + } if (!strcasecmp (breakPt, "CLASS")) exit (0); + // mark blended peaks PS_SOURCE_BLEND psphotBasicDeblend (sources, config, sky); if (!strcasecmp (breakPt, "DEBLEND")) exit (0); @@ -57,24 +66,22 @@ if (!strcasecmp (breakPt, "PSFFIT")) exit (0); + // XXX add this as conditional output + psphotSamplePSFs (psf, imdata->image); + int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE"); if (!status) FITMODE = 2; switch (FITMODE) { - case -2: + case 0: psphotEnsemblePSF (imdata, config, sources, psf, sky); - psphotReapplyPSF (imdata, config, sources, psf, sky); - break; - - case -1: - psphotEnsemblePSF (imdata, config, sources, psf, sky); - break; - - case 0: - psphotFixedPSF (imdata, config, sources, psf, sky); + sources = psArraySort (sources, psphotSortBySN); + psphotFullFit (imdata, config, sources, psf, sky); +// psphotReplaceUnfit (sources); break; case 1: - // test PSF on all except SATURATE and DEFECT (artifacts) + // fit extended objects with galaxy models psphotApplyPSF (imdata, config, sources, psf, sky); + psphotFitGalaxies (imdata, config, sources, sky); break; @@ -82,5 +89,4 @@ // fit extended objects with galaxy models psphotApplyPSF (imdata, config, sources, psf, sky); - psphotFitGalaxies (imdata, config, sources, sky); break; Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 5771) +++ /trunk/psphot/src/psphot.h (revision 5772) @@ -74,2 +74,13 @@ psMetadata *psphotTestArguments (int *argc, char **argv); bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky); + +bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotInitLimitsPSF (psMetadata *config); +bool psphotEvalPSF (pmSource *source); +bool psphotEvalFLT (pmSource *source); +bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type); +bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source); +bool psphotInitRadiusFLT (psMetadata *config, psStats *sky, pmModelType type); +bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source); +bool psphotSamplePSFs (pmPSF *psf, psImage *image); +bool psphotReplaceUnfit (psArray *sources); Index: /trunk/psphot/src/psphotApplyPSF.c =================================================================== --- /trunk/psphot/src/psphotApplyPSF.c (revision 5771) +++ /trunk/psphot/src/psphotApplyPSF.c (revision 5772) @@ -2,34 +2,16 @@ // fit psf model to all objects -// PSFSTAR objects will be refitted -// run this function to a specific flux limit? - bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { - bool status; float x; float y; - int Nfit = 0; - int Nsub = 0; + int Nfit = 0; + int Nsub = 0; int Niter = 0; psTimerStart ("psphot"); - // we may set this differently here from the value used to mark likely saturated stars - float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); - float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); - - float PSF_MIN_SN = psMetadataLookupF32 (&status, config, "PSF_MIN_SN"); - float PSF_MAX_CHI = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI"); - float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); - float PSF_FIT_PADDING = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING"); - float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); - - // set the object surface-brightness limit for fitted pixels - float FLUX_LIMIT = PSF_FIT_NSIGMA * sky->sampleStdev; - psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); - - // this function specifies the radius at this the model hits the given flux - pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type); + psphotInitLimitsPSF (config); + psphotInitRadiusPSF (config, sky, psf->type); for (int i = 0; i < sources->n; i++) { @@ -38,5 +20,5 @@ // skip non-astronomical objects (very likely defects) - // XXX EAM : should we try these anyway? + if (source->mode & PM_SOURCE_BLEND) continue; if (source->type == PM_SOURCE_DEFECT) continue; if (source->type == PM_SOURCE_SATURATED) continue; @@ -47,52 +29,32 @@ // set PSF parameters for this model pmModel *model = pmModelFromPSF (modelFLT, psf); + psFree (modelFLT); + + source->modelPSF = model; + + // sets the model radius (via source->model) and adjusts pixels as needed + psphotCheckRadiusPSF (imdata, source); + x = model->params->data.F32[2]; y = model->params->data.F32[3]; - psFree (modelFLT); - // set the fit radius based on the object flux limit and the model - // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median) - model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING; - if (isnan(model->radius)) { - psAbort ("apply_psf_model", "error in radius"); - } - - // check if we need to redefine the pixels - // XXX EAM : a better test would examine the source pixels - if (model->radius > OUTER_RADIUS) { - // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER) - psphotDefinePixels (source, imdata, x, y, model->radius); - } - - // if (i > 66) psTraceSetLevel (".psLib.dataManip.psMinimizeLMChi2", 5); - // fit PSF model (set/unset the pixel mask) psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitModel (source, model, true); + pmSourceFitModel (source, model, true); psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); - if (!status || (model->params->data.F32[1] < 0)) { - psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter); - source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? - psFree (model); - continue; - } - // XXX EAM : this was an attempt to correct for fit-sky biases - // pmModelSkyOffset (model, x, y, model->radius); - - // model succeeded : tentatively keep it - source->modelPSF = model; Niter += model[0].nIter; Nfit ++; - // is it a good model? - psphotMarkPSF (source, PSF_SHAPE_NSIGMA, PSF_MIN_SN, PSF_MAX_CHI, SATURATION); - if (psphotSubtractPSF (source)) { - Nsub ++; + // check if model fit is acceptable + if (psphotEvalPSF (source)) { + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + Nsub ++; } } psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); - psLogMsg ("psphot", 4, "subtracted %d PSF objects\n", Nsub); + psLogMsg ("psphot", 4, "subtracted %d PSF models\n", Nsub); return (true); } Index: /trunk/psphot/src/psphotBasicDeblend.c =================================================================== --- /trunk/psphot/src/psphotBasicDeblend.c (revision 5771) +++ /trunk/psphot/src/psphotBasicDeblend.c (revision 5772) @@ -14,11 +14,11 @@ psTimerStart ("psphot"); - // this should be added to the PM_SOURCE flags: - int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1; + float FRACTION = psMetadataLookupF32 (&status, config, "DEBLEND_PEAK_FRACTION"); + if (!status) FRACTION = 0.25; - float FRACTION = psMetadataLookupF32 (&status, config, "DEBLEND_PEAK_FRACTION"); float NSIGMA = psMetadataLookupF32 (&status, config, "DEBLEND_SKY_NSIGMA"); + if (!status) NSIGMA = 5.0; + float minThreshold = NSIGMA*sky->sampleStdev; - fprintf (stderr, "min threshold: %f\n", minThreshold); // we need sources spatially-sorted to find overlaps @@ -43,5 +43,5 @@ source = sources->data[N]; - if (source->type == PM_SOURCE_BLEND) continue; + if (source->mode & PM_SOURCE_BLEND) continue; overlap->n = 0; @@ -119,5 +119,5 @@ ); - testSource->type = PM_SOURCE_BLEND; + testSource->mode |= PM_SOURCE_BLEND; Nblend ++; j = xv->n; Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 5771) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 5772) @@ -13,5 +13,6 @@ // array to store candidate PSF stars int NSTARS = psMetadataLookupS32 (&status, config, "PSF_MAX_NSTARS"); - if (!status) NSTARS = sources->n; + if (!status) NSTARS = PS_MIN (sources->n, 200); + stars = psArrayAlloc (sources->n); stars->n = 0; @@ -20,8 +21,7 @@ for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) { pmSource *source = sources->data[i]; - if (source->type != PM_SOURCE_PSFSTAR) continue; - psArrayAdd (stars, 200, source); + if (source->mode & PM_SOURCE_PSFSTAR) psArrayAdd (stars, 200, source); } - psTrace (".psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n); + psLogMsg ("psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n); // get the fixed PSF fit radius @@ -66,4 +66,6 @@ } } + // XXX EAM : need to unflag the PSF stars which are not used to build the PSF + // XXX EAM : each pmPSFtry needs to have its own mask array // keep only the selected model: Index: /trunk/psphot/src/psphotEnsemblePSF.c =================================================================== --- /trunk/psphot/src/psphotEnsemblePSF.c (revision 5771) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 5772) @@ -12,7 +12,4 @@ // source analysis is done in spatial order sources = psArraySort (sources, psphotSortByY); - - // this should be added to the PM_SOURCE flags: - int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1; float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); @@ -35,5 +32,5 @@ // skip non-astronomical objects (very likely defects) // XXX EAM : should we try these anyway? - if (inSource->type == PM_SOURCE_BLEND) continue; + if (inSource->mode & PM_SOURCE_BLEND) continue; if (inSource->type == PM_SOURCE_DEFECT) continue; if (inSource->type == PM_SOURCE_SATURATED) continue; @@ -73,7 +70,10 @@ psImage *mask = otSource->mask; - // XXX EAM : set model to unit peak, zero sky, maybe use peak (x,y) + // XXX EAM : use these lines to fit to the peak // model->params->data.F32[2] = inSource->peak->x; // model->params->data.F32[3] = inSource->peak->y; + // XXX EAM : better option: improve the peak with 2D poly fit 3x3 + + // set model to unit peak, zero sky (we assume sky is constant) model->params->data.F32[0] = 0.0; model->params->data.F32[1] = 1.0; Index: /trunk/psphot/src/psphotEvalFLT.c =================================================================== --- /trunk/psphot/src/psphotEvalFLT.c (revision 5772) +++ /trunk/psphot/src/psphotEvalFLT.c (revision 5772) @@ -0,0 +1,51 @@ +# include "psphot.h" + +bool psphotEvalFLT (pmSource *source) +{ + int keep; + + pmModel *model = source->modelFLT; + + // do we actually have a valid FLT model? + if (model == NULL) { + source->mode &= ~PM_SOURCE_FITTED; + return false; + } + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_UNTRIED: + source->mode &= ~PM_SOURCE_FITTED; + return false; + case PM_MODEL_BADARGS: + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + default: + psLogMsg ("psphot", 5, "GAL fail fit\n"); + source->mode |= PM_SOURCE_FAIL; + return false; + } + + // unless we prove otherwise, this object is a galaxy. + // XXX EAM : do we need to save the old type? + source->type = PM_SOURCE_GALAXY; + + // the following source->mode information pertains to modelPSF: + source->mode |= PM_SOURCE_FLTMODEL; + + // if the object has a fitted peak below 0, the fit did not converge cleanly + if (model->params->data.F32[1] < 0) { + source->mode |= PM_SOURCE_FAIL; + return false; + } + + keep = pmModelFitStatus (model); + if (keep) return true; + + // poor-quality fit; only keep if nothing else works... + psLogMsg ("psphot", 5, "GAL poor fit\n"); + source->mode |= PM_SOURCE_POOR; + return false; +} Index: /trunk/psphot/src/psphotEvalPSF.c =================================================================== --- /trunk/psphot/src/psphotEvalPSF.c (revision 5772) +++ /trunk/psphot/src/psphotEvalPSF.c (revision 5772) @@ -0,0 +1,150 @@ +# include "psphot.h" + +// given a pmSource which has been fitted using modelPSF, evaluate the +// resulting fit: did the fit succeed? is this object PSF-like? is this object +// extended? return status is TRUE for a valid PSF, FALSE otherwise. + +// identify objects consistent with PSF shape/magnitude distribution +// we expect dparams[4],dparams[5] to have a scatter of: +// sigma_x / (S/N) * sqrt(2) +// 1 / (params[4],params[5])*(S/N) + +// sigma_x : 1 / SX +// dsx : 1 / (SX * SN) +// dsx_o = hypot (1/(SX*SN), MIN_DSX) + +// any objects which is consistent with the PSF should have +// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag) +// this includes a minimum buffer (DS) for the brighter objects + +// saturated stars should fall outside (larger), but have peaks above SATURATION +// galaxies should be larger, cosmic rays smaller +// we also reject objects with S/N too low or ChiSquare to high + +// floor for DS value +// XXX EAM : add to configuration? +# define MIN_DS 0.01 + +static float SATURATION; +static float PSF_MIN_SN; +static float PSF_MAX_CHI; +static float PSF_SHAPE_NSIGMA; + +bool psphotInitLimitsPSF (psMetadata *config) { + + bool status; + + // we may set this differently here from the value used to mark likely saturated stars + SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); + PSF_MIN_SN = psMetadataLookupF32 (&status, config, "PSF_MIN_SN"); + PSF_MAX_CHI = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI"); + PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); + + return true; +} + +// examine the model->status, fit parameters, etc and decide if the model succeeded +// set the source->type and source->mode appropriately +bool psphotEvalPSF (pmSource *source) { + + int keep; + float dSX, dSY, SX, SY, SN; + float nSx, nSy, Chi; + + pmModel *model = source->modelPSF; + + // do we actually have a valid PSF model? + if (model == NULL) { + source->mode &= ~PM_SOURCE_FITTED; + return false; + } + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_UNTRIED: + source->mode &= ~PM_SOURCE_FITTED; + return false; + case PM_MODEL_BADARGS: + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + default: + source->mode |= PM_SOURCE_FAIL; + return false; + } + + // unless we prove otherwise, this object is a star. + source->type = PM_SOURCE_STAR; + + // the following source->mode information pertains to modelPSF: + source->mode |= PM_SOURCE_PSFMODEL; + + // if the object has fitted peak above saturation, label as SATSTAR + // this is a valid PSF object, but ignore the other quality tests + // remember: fit does not use saturated pixels (masked) + // XXX no extended object can saturate and stay extended... + if (model->params->data.F32[1] >= SATURATION) { + if (source->mode & PM_SOURCE_PSFSTAR) { + psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); + } + source->mode |= PM_SOURCE_SATSTAR; + return true; + } + + // if the object has a fitted peak below 0, the fit did not converge cleanly + if (model->params->data.F32[1] < 0) { + source->mode |= PM_SOURCE_FAIL; + return false; + } + + // if the source was predicted to be a SATSTAR, but it fitted below saturation, + // make a note to the user + if (source->mode & PM_SOURCE_SATSTAR) { + psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n"); + source->mode &= ~PM_SOURCE_SATSTAR; + } + + SN = model->params->data.F32[1]/model->dparams->data.F32[1]; + SX = model->params->data.F32[4]; + SY = model->params->data.F32[5]; + dSX = model->dparams->data.F32[4]; + dSY = model->dparams->data.F32[5]; + Chi = model->chisq / model->nDOF; + + // swing of sigma_x,y in sigmas + nSx = dSX / hypot (MIN_DS, 1 / (SX * SN)); + nSy = dSY / hypot (MIN_DS, 1 / (SY * SN)); + + // assign PM_SOURCE_GOODSTAR to bright objects within PSF region of dparams[] + keep = TRUE; + keep &= (fabs(nSx) < PSF_SHAPE_NSIGMA); + keep &= (fabs(nSy) < PSF_SHAPE_NSIGMA); + keep &= (SN > PSF_MIN_SN); + keep &= (Chi < PSF_MAX_CHI); + if (keep) return true; + + // this source is not a star, unflag as PSFSTAR + // XXX : if this object was used to build the PSF, this flag should + // be set even if the object is not a star... + if (source->mode & PM_SOURCE_PSFSTAR) { + source->mode &= ~PM_SOURCE_PSFSTAR; + psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality\n"); + } + + // object appears to be small, suspected defect + if ((nSx <= -PSF_SHAPE_NSIGMA) || (nSy <= -PSF_SHAPE_NSIGMA)) { + source->type = PM_SOURCE_DEFECT; + return false; + } + + // object appears to be large, suspected galaxy + if ((nSx >= PSF_SHAPE_NSIGMA) || (nSy >= PSF_SHAPE_NSIGMA)) { + source->type = PM_SOURCE_GALAXY; + return false; + } + + // poor-quality fit; only keep if nothing else works... + source->mode |= PM_SOURCE_POOR; + return false; +} Index: /trunk/psphot/src/psphotFitGalaxies.c =================================================================== --- /trunk/psphot/src/psphotFitGalaxies.c (revision 5771) +++ /trunk/psphot/src/psphotFitGalaxies.c (revision 5772) @@ -3,53 +3,34 @@ bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats) { - bool status, goodfit; + bool status; float x; float y; - int Nfit = 0; - int Nfail = 0; + int Nfit = 0; + int Nsub = 0; int Niter = 0; - float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); + psTimerStart ("psphot"); float GAL_MIN_SN = psMetadataLookupF32 (&status, config, "GAL_MIN_SN"); - float GAL_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "GAL_FIT_NSIGMA"); - float GAL_FIT_PADDING = psMetadataLookupF32 (&status, config, "GAL_FIT_PADDING"); float GAL_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS"); + char *modelName = psMetadataLookupPtr (&status, config, "GAL_MODEL"); + pmModelType modelType = pmModelSetType (modelName); - float FLUX_LIMIT = GAL_FIT_NSIGMA * skyStats->sampleStdev; - - char *modelName = psMetadataLookupPtr (&status, config, "GAL_MODEL"); - pmModelType modelType = pmModelSetType (modelName); - pmModelRadius modelRadius = pmModelRadius_GetFunction (modelType); - - psTimerStart ("psphot"); + psphotInitRadiusFLT (config, skyStats, modelType); for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - // sources which should not be fitted - // skip all valid stars - if (source->type == PM_SOURCE_PSFSTAR) continue; - if (source->type == PM_SOURCE_SATSTAR) continue; - if (source->type == PM_SOURCE_GOODSTAR) continue; + // only fit suspected extended sources + if (source->type != PM_SOURCE_GALAXY) continue; + if (source->mode & PM_SOURCE_BLEND) continue; - // skip all likely defects - if (source->type == PM_SOURCE_DEFECT) continue; - if (source->type == PM_SOURCE_SATURATED) continue; - - // skip poorly fitted stars - if (source->type == PM_SOURCE_FAINTSTAR) continue; - if (source->type == PM_SOURCE_POOR_FIT_PSF) continue; - - // XXX EAM when do we pick these up again? - if (source->moments->SN < GAL_MIN_SN) { - source->type = PM_SOURCE_FAINT_GALAXY; // better choice? - continue; - } + // skip possible extended sources below fit threshold + if (source->moments->SN < GAL_MIN_SN) continue; // recalculate the source moments using the larger galaxy moments radius status = pmSourceMoments (source, GAL_MOMENTS_RAD); if (!status) { - source->type = PM_SOURCE_DROP_GALAXY; // better choice? + source->mode |= PM_SOURCE_FAIL; // better choice? continue; } @@ -57,53 +38,28 @@ // use the source moments, etc to guess basic model parameters pmModel *model = pmSourceModelGuess (source, modelType); + source->modelFLT = model; x = model->params->data.F32[2]; y = model->params->data.F32[3]; - // set the fit radius based on the object flux limit and the model - // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median) - model->radius = modelRadius (model->params, FLUX_LIMIT) + GAL_FIT_PADDING; - if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius"); - - // check if we need to redefine the pixels - // XXX EAM : a better test would examine the source pixels - if (model->radius > OUTER_RADIUS) { - // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER) - psphotDefinePixels (source, imdata, x, y, model->radius); - } + // sets the model radius (via source->model) and adjusts pixels as needed + psphotCheckRadiusFLT (imdata, source); // fit FLT (not PSF) model (set/unset the pixel mask) psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitModel (source, model, false); + pmSourceFitModel (source, model, false); psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); - if (!status) { - // if the fit fails, we need to change the classification - psLogMsg ("psphot", 5, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); - source->type = PM_SOURCE_FAIL_FIT_GAL; // better choice? - source->modelFLT = model; - Nfail ++; - continue; - } - // check if model fit is acceptable - goodfit = pmModelFitStatus (model); - if (!goodfit) { - // if the fit fails, we need to change the classification - psLogMsg ("psphot", 5, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); - source->type = PM_SOURCE_POOR_FIT_GAL; // better choice? - source->modelFLT = model; - Nfail ++; - continue; - } - - // accept the model - source->type = PM_SOURCE_GALAXY; - source->modelFLT = model; Niter += model[0].nIter; Nfit++; - // subtract object, leave local sky - pmSourceSubModel (source->pixels, source->mask, model, false, false); + // check if model fit is acceptable + if (pmModelFitStatus (model)) { + pmSourceSubModel (source->pixels, source->mask, model, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + Nsub ++; + } } - psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter); + psLogMsg ("psphot", 3, "fit GAL models: %f sec for %d galaxies (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); + psLogMsg ("psphot", 4, "subtracted %d GAL models\n", Nsub); return (true); } Index: /trunk/psphot/src/psphotFixedPSF.c =================================================================== --- /trunk/psphot/src/psphotFixedPSF.c (revision 5771) +++ /trunk/psphot/src/psphotFixedPSF.c (revision 5772) @@ -66,5 +66,5 @@ if (!status || (model->params->data.F32[1] < 0)) { psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (peak %f\n", x, y, model->params->data.F32[1]); - source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? + source->mode |= PM_SOURCE_FAIL; // better choice? psFree (model); continue; Index: /trunk/psphot/src/psphotFullFit.c =================================================================== --- /trunk/psphot/src/psphotFullFit.c (revision 5772) +++ /trunk/psphot/src/psphotFullFit.c (revision 5772) @@ -0,0 +1,128 @@ +# include "psphot.h" + +bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) +{ + bool skip, status; + float x, y; + int Nfit = 0; + int Nsub = 0; + int Niter = 0; + + psTimerStart ("psphot"); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + // galaxy model parameters + float GAL_MIN_SN = psMetadataLookupF32 (&status, config, "GAL_MIN_SN"); + float GAL_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS"); + + // 'galaxy' model descriptions + char *modelNameFLT = psMetadataLookupPtr (&status, config, "GAL_MODEL"); + pmModelType modelTypeFLT = pmModelSetType (modelNameFLT); + + psphotInitLimitsPSF (config); + psphotInitRadiusPSF (config, sky, psf->type); + psphotInitRadiusFLT (config, sky, modelTypeFLT); + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + // XXX EAM : should we try these anyway? + if (source->mode & PM_SOURCE_BLEND) continue; + if (source->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; + + // save the PSF model from the Ensemble fit + pmModel *save = pmModelCopy (source->modelPSF); + + // attempt to fit the PSF model + pmModel *modelPSF = source->modelPSF; + + // replace object in image + pmSourceAddModel (source->pixels, source->mask, modelPSF, false, false); + + psphotCheckRadiusPSF (imdata, source); + + x = modelPSF->params->data.F32[2]; + y = modelPSF->params->data.F32[3]; + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, modelPSF->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitModel (source, modelPSF, true); + psImageKeepCircle (source->mask, x, y, modelPSF->radius, "AND", ~PSPHOT_MASK_MARKED); + + // track model evaluations + Niter += modelPSF[0].nIter; + Nfit ++; + + // does the PSF model succeed? + if (psphotEvalPSF (source)) { + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + psFree (save); + Nsub ++; + continue; + } + + // skip the source if we don't think it is extended + skip = (source->type != PM_SOURCE_GALAXY); + skip |= (source->moments->SN < GAL_MIN_SN); + if (skip) { + source->modelPSF = save; + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + psFree (modelPSF); + Nsub ++; + continue; + } + + // recalculate the source moments using the larger galaxy moments radius + if (!pmSourceMoments (source, GAL_MOMENTS_RAD)) { + source->mode |= PM_SOURCE_FAIL; // better choice? + continue; + } + + // use the source moments, etc to guess basic model parameters + source->modelFLT = pmSourceModelGuess (source, modelTypeFLT); + pmModel *modelFLT = source->modelFLT; + + psphotCheckRadiusFLT (imdata, source); + + x = modelFLT->params->data.F32[2]; + y = modelFLT->params->data.F32[3]; + + // fit FLT (not PSF) model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, modelFLT->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitModel (source, modelFLT, false); + psImageKeepCircle (source->mask, x, y, modelFLT->radius, "AND", ~PSPHOT_MASK_MARKED); + + // track model evaluations + Niter += modelFLT[0].nIter; + Nfit++; + + // does the FLT model succeed? + if (psphotEvalFLT (source)) { + pmSourceSubModel (source->pixels, source->mask, source->modelFLT, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + psFree (save); + Nsub ++; + continue; + } + + // subtract object, leave local sky + psFree (source->modelPSF); + source->modelPSF = save; + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + source->mode |= PM_SOURCE_SUBTRACTED; + } + + psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); + psLogMsg ("psphot", 4, "subtracted %d PSF objects\n", Nsub); + return (true); +} + + +// XXX EAM : can we distinguish the outcomes enough to re-add poor fits, etc? Index: /trunk/psphot/src/psphotImageBackground.c =================================================================== --- /trunk/psphot/src/psphotImageBackground.c (revision 5771) +++ /trunk/psphot/src/psphotImageBackground.c (revision 5772) @@ -64,9 +64,9 @@ } } + sky->sampleMean = 0.0; psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot")); psphotSaveImage (imdata->header, imdata->image, "backsub.fits"); - return (skyModel); } Index: /trunk/psphot/src/psphotMagnitudes.c =================================================================== --- /trunk/psphot/src/psphotMagnitudes.c (revision 5771) +++ /trunk/psphot/src/psphotMagnitudes.c (revision 5772) @@ -49,7 +49,5 @@ // use PSF model with stars - case PM_SOURCE_PSFSTAR: - case PM_SOURCE_SATSTAR: - case PM_SOURCE_GOODSTAR: + case PM_SOURCE_STAR: model = source->modelPSF; break; @@ -61,10 +59,4 @@ // skip defects and poorly fitted stars or galaxies - case PM_SOURCE_DEFECT: - case PM_SOURCE_SATURATED: - case PM_SOURCE_FAINTSTAR: - case PM_SOURCE_POOR_FIT_PSF: - case PM_SOURCE_POOR_FIT_GAL: - case PM_SOURCE_FAIL_FIT_GAL: default: model = NULL; Index: /trunk/psphot/src/psphotMarkPSF.c =================================================================== --- /trunk/psphot/src/psphotMarkPSF.c (revision 5771) +++ /trunk/psphot/src/psphotMarkPSF.c (revision 5772) @@ -1,3 +1,7 @@ # include "psphot.h" + +// given a pmSource which has been fitted using modelPSF, evaluate the +// resulting fit: did the fit succeed? is this object PSF-like? is this object +// extended? return status is TRUE for a valid PSF, FALSE otherwise. // identify objects consistent with PSF shape/magnitude distribution @@ -14,9 +18,7 @@ // this includes a minimum buffer (DS) for the brighter objects -// saturated stars should fall outside (but are already IDed) -// galaxies should be larger, cosmic rays smaller, but need to test? +// saturated stars should fall outside (larger), but have peaks above SATURATE +// galaxies should be larger, cosmic rays smaller // we also reject objects with S/N too low or ChiSquare to high - -// any object which meets the criterion is marked as PM_SOURCE_BRIGHTSTAR // floor for DS value @@ -33,11 +35,16 @@ // do we actually have a valid PSF model? - if (model == NULL) return (false); + if (model == NULL) { + source->mode &= ~PM_SOURCE_FITTED; + return false; + } // did the model fit fail for one or another reason? switch (model->status) { + case PM_MODEL_UNTRIED: + source->mode &= ~PM_SOURCE_FITTED; + return false; case PM_MODEL_BADARGS: - case PM_MODEL_UNTRIED: - source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? + source->mode |= PM_SOURCE_FAIL; return false; case PM_MODEL_SUCCESS: @@ -46,5 +53,5 @@ case PM_MODEL_OFFIMAGE: psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", model->params->data.F32[2], model->params->data.F32[3], model->nIter); - source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? + source->mode |= PM_SOURCE_FAIL; return false; default: @@ -52,26 +59,31 @@ } + // unless we prove otherwise, this object is a star. + source->type = PM_SOURCE_STAR; + // if the object has fitted peak above saturation, label as SATSTAR // this is a valid PSF object, but ignore the other quality tests // remember: fit does not use saturated pixels (masked) + // XXX no extended object can saturate and stay extended... if (model->params->data.F32[1] >= SATURATE) { - if (source->type == PM_SOURCE_PSFSTAR) { + if (source->mode & PM_SOURCE_PSFSTAR) { psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); } - source->type = PM_SOURCE_SATSTAR; - return (true); + source->mode &= ~PM_SOURCE_PSFSTAR; + source->mode |= PM_SOURCE_SATSTAR; + return true; } // if the object has a fitted peak below 0, the fit did not converge cleanly if (model->params->data.F32[1] < 0) { - source->type = PM_SOURCE_FAIL_FIT_PSF; - return (false); + source->mode |= PM_SOURCE_FAIL; + return false; } // if the source was predicted to be a SATSTAR, but it fitted below saturation, // make a note to the user - if (source->type == PM_SOURCE_SATSTAR) { - psLogMsg ("psphot", 5, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n"); - source->type = PM_SOURCE_GOODSTAR; + if (source->mode & PM_SOURCE_SATSTAR) { + psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n"); + source->mode &= ~PM_SOURCE_SATSTAR; } @@ -93,11 +105,11 @@ keep &= (SN > minSN); keep &= (Chi < maxChi); - if (keep) { - if (source->type == PM_SOURCE_PSFSTAR) return (true); - source->type = PM_SOURCE_GOODSTAR; - return (true); - } - - if (source->type == PM_SOURCE_PSFSTAR) { + if (keep) return true; + + // this source is not a star, unflag as PSFSTAR + // XXX : if this object was used to build the PSF, this flag should + // be set even if the object is not a star... + if (source->mode & PM_SOURCE_PSFSTAR) { + source->mode &= ~PM_SOURCE_PSFSTAR; psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality\n"); } @@ -106,5 +118,5 @@ if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) { source->type = PM_SOURCE_DEFECT; - return (false); + return false; } @@ -112,15 +124,9 @@ if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) { source->type = PM_SOURCE_GALAXY; - return (false); + return false; } - // object appears to be extremely faint: what is this? - if (SN < minSN) { - source->type = PM_SOURCE_FAINTSTAR; - return (false); - } - - // these are pooly fitted, probable stars near other stars? - source->type = PM_SOURCE_POOR_FIT_PSF; - return (false); + // poor-quality fit; only keep if nothing else works... + source->mode |= PM_SOURCE_POOR; + return false; } Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 5771) +++ /trunk/psphot/src/psphotOutput.c (revision 5772) @@ -389,11 +389,9 @@ pmSource *source = (pmSource *) sources->data[i]; - // valid source types for this function - if (source->type == PM_SOURCE_SATSTAR) goto valid; - if (source->type == PM_SOURCE_PSFSTAR) goto valid; - if (source->type == PM_SOURCE_GOODSTAR) goto valid; - continue; - - valid: + // write out sources fitted as PSFs + if (source->type != PM_SOURCE_STAR) continue; + if (source->mode & PM_SOURCE_FAIL) continue; + if (source->mode & PM_SOURCE_POOR) continue; + model = pmSourceMagnitudes (source, psf, RADIUS); if (model == NULL) continue; @@ -452,8 +450,8 @@ model = (pmModel *) source->modelFLT; if (model == NULL) continue; - if (source->type == PM_SOURCE_GALAXY) goto valid; - continue; - - valid: + if (source->type != PM_SOURCE_GALAXY) continue; + if (source->mode & PM_SOURCE_FAIL) continue; + if (source->mode & PM_SOURCE_POOR) continue; + params = model->params; dparams = model->dparams; @@ -513,9 +511,17 @@ // skip these sources (in PSF or FLT) - if (source->type == PM_SOURCE_GALAXY) continue; - if (source->type == PM_SOURCE_PSFSTAR) continue; - if (source->type == PM_SOURCE_SATSTAR) continue; - if (source->type == PM_SOURCE_GOODSTAR) continue; - + if (source->type == PM_SOURCE_STAR) { + if (source->mode & PM_SOURCE_FAIL) goto isNULL; + if (source->mode & PM_SOURCE_POOR) goto isNULL; + continue; + } + if (source->type == PM_SOURCE_GALAXY) { + if (source->mode & PM_SOURCE_FAIL) goto isNULL; + if (source->mode & PM_SOURCE_POOR) goto isNULL; + continue; + } + + isNULL: + if (source->moments == NULL) { fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %8.1f %7.1f %7.1f %7.1f %4d %2d\n", @@ -579,26 +585,13 @@ return (8); - case PM_SOURCE_SATSTAR: - return (10); - - case PM_SOURCE_PSFSTAR: - case PM_SOURCE_GOODSTAR: + case PM_SOURCE_STAR: + if (source->mode & PM_SOURCE_SATSTAR) return (10); + if (source->mode & PM_SOURCE_POOR) return (7); + if (source->mode & PM_SOURCE_FAIL) return (4); return (1); - case PM_SOURCE_POOR_FIT_PSF: - return (7); - - case PM_SOURCE_FAIL_FIT_PSF: - case PM_SOURCE_FAINTSTAR: - return (4); - case PM_SOURCE_GALAXY: - case PM_SOURCE_FAINT_GALAXY: - case PM_SOURCE_DROP_GALAXY: - case PM_SOURCE_FAIL_FIT_GAL: - case PM_SOURCE_POOR_FIT_GAL: return (2); - case PM_SOURCE_OTHER: default: return (0); Index: /trunk/psphot/src/psphotRadiusChecks.c =================================================================== --- /trunk/psphot/src/psphotRadiusChecks.c (revision 5772) +++ /trunk/psphot/src/psphotRadiusChecks.c (revision 5772) @@ -0,0 +1,78 @@ +# include "psphot.h" + +static float PSF_OUTER_RADIUS; +static float PSF_FIT_PADDING; +static float PSF_FLUX_LIMIT; +static pmModelRadius modelRadiusPSF; + +bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type) { + + bool status; + + float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); + PSF_OUTER_RADIUS = psMetadataLookupF32 (&status, config, "PSF_OUTER_RADIUS"); + PSF_FIT_PADDING = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING"); + PSF_FLUX_LIMIT = PSF_FIT_NSIGMA * sky->sampleStdev; + + // this function specifies the radius at this the model hits the given flux + modelRadiusPSF = pmModelRadius_GetFunction (type); + return true; +} + +bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source) { + + pmModel *model = source->modelPSF; + + // set the fit radius based on the object flux limit and the model + // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median) + model->radius = modelRadiusPSF (model->params, PSF_FLUX_LIMIT) + PSF_FIT_PADDING; + if (isnan(model->radius)) psAbort ("apply_psf_model", "error in radius"); + + // check if we need to redefine the pixels + // XXX EAM : a better test would examine the source pixels + if (model->radius > PSF_OUTER_RADIUS) { + // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER) + psphotDefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius); + return true; + } + + return false; +} + + +static float GAL_OUTER_RADIUS; +static float GAL_FIT_PADDING; +static float GAL_FLUX_LIMIT; +static pmModelRadius modelRadiusGAL; + +bool psphotInitRadiusFLT (psMetadata *config, psStats *sky, pmModelType type) { + + bool status; + + float GAL_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "GAL_FIT_NSIGMA"); + GAL_OUTER_RADIUS = psMetadataLookupF32 (&status, config, "GAL_OUTER_RADIUS"); + GAL_FIT_PADDING = psMetadataLookupF32 (&status, config, "GAL_FIT_PADDING"); + GAL_FLUX_LIMIT = GAL_FIT_NSIGMA * sky->sampleStdev; + + // this function specifies the radius at this the model hits the given flux + modelRadiusGAL = pmModelRadius_GetFunction (type); + return true; +} + +bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source) { + + pmModel *model = source->modelFLT; + + // set the fit radius based on the object flux limit and the model + model->radius = modelRadiusGAL (model->params, GAL_FLUX_LIMIT) + GAL_FIT_PADDING; + if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius"); + + // check if we need to redefine the pixels + if (model->radius > GAL_OUTER_RADIUS) { + // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER) + // do I need to save the old pixel and mask data? + psphotDefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius); + return true; + } + return false; +} Index: /trunk/psphot/src/psphotReapplyPSF.c =================================================================== --- /trunk/psphot/src/psphotReapplyPSF.c (revision 5771) +++ /trunk/psphot/src/psphotReapplyPSF.c (revision 5772) @@ -14,7 +14,4 @@ sources = psArraySort (sources, psphotSortBySN); - // this should be added to the PM_SOURCE flags: - int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1; - // we may set this differently here from the value used to mark likely saturated stars float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); @@ -29,5 +26,5 @@ // skip non-astronomical objects (very likely defects) // XXX EAM : should we try these anyway? - if (source->type == PM_SOURCE_BLEND) continue; + if (source->mode |= PM_SOURCE_BLEND) continue; if (source->type == PM_SOURCE_DEFECT) continue; if (source->type == PM_SOURCE_SATURATED) continue; @@ -70,17 +67,2 @@ return (true); } - -pmModel *pmModelCopy (pmModel *model) { - - pmModel *new = pmModelAlloc (model->type); - - new->chisq = model->chisq; - new->nIter = model->nIter; - - for (int i = 0; i < new->params->n; i++) { - new->params->data.F32[i] = model->params->data.F32[i]; - new->dparams->data.F32[i] = model->dparams->data.F32[i]; - } - - return (new); -}