Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 17288) +++ trunk/psphot/src/Makefile.am (revision 17396) @@ -18,58 +18,59 @@ libpsphot_la_SOURCES = \ - psphotErrorCodes.c \ - pmFootprint.c \ - psphotVersion.c \ - psphotModelGroupInit.c \ - psphotMaskReadout.c \ - psphotDefineFiles.c \ - psphotReadout.c \ - psphotModelBackground.c \ - psphotSubtractBackground.c \ - psphotFindDetections.c \ - psphotFindPeaks.c \ - psphotSourceStats.c \ - psphotRoughClass.c \ - psphotBasicDeblend.c \ - psphotChoosePSF.c \ - psphotGuessModels.c \ - psphotFitSourcesLinear.c \ - psphotBlendFit.c \ - psphotReplaceUnfit.c \ - psphotApResid.c \ - psphotMagnitudes.c \ - psphotSkyReplace.c \ - psphotEvalPSF.c \ - psphotEvalFLT.c \ - psphotSourceFits.c \ - psphotRadiusChecks.c \ - psphotOutput.c \ - psphotGrowthCurve.c \ - psphotFakeSources.c \ - psphotModelWithPSF.c \ - psphotExtendedSources.c \ - psphotRadialProfile.c \ - psphotPetrosian.c \ - psphotIsophotal.c \ - psphotAnnuli.c \ - psphotKron.c \ - psphotKernelFromPSF.c \ - psphotPSFConvModel.c \ - psphotModelTest.c \ - psphotFitSet.c \ - psphotSourceFreePixels.c \ - psphotSummaryPlots.c \ - psphotMergeSources.c \ - psphotLoadPSF.c \ - psphotReadoutCleanup.c \ - psphotSourcePlots.c \ - psphotRadialPlot.c \ - psphotDeblendSatstars.c \ - psphotMosaicSubimage.c \ - psphotMakeResiduals.c \ - psphotSourceSize.c \ - psphotDiagnosticPlots.c \ - psphotMakeFluxScale.c \ - psphotCheckStarDistribution.c \ + psphotErrorCodes.c \ + pmFootprint.c \ + psphotVersion.c \ + psphotModelGroupInit.c \ + psphotMaskReadout.c \ + psphotDefineFiles.c \ + psphotReadout.c \ + psphotModelBackground.c \ + psphotSubtractBackground.c \ + psphotFindDetections.c \ + psphotFindPeaks.c \ + psphotSourceStats.c \ + psphotRoughClass.c \ + psphotBasicDeblend.c \ + psphotChoosePSF.c \ + psphotGuessModels.c \ + psphotFitSourcesLinear.c \ + psphotBlendFit.c \ + psphotReplaceUnfit.c \ + psphotApResid.c \ + psphotMagnitudes.c \ + psphotSkyReplace.c \ + psphotEvalPSF.c \ + psphotEvalFLT.c \ + psphotSourceFits.c \ + psphotRadiusChecks.c \ + psphotOutput.c \ + psphotGrowthCurve.c \ + psphotFakeSources.c \ + psphotModelWithPSF.c \ + psphotExtendedSourceAnalysis.c \ + psphotExtendedSourceFits.c \ + psphotRadialProfile.c \ + psphotPetrosian.c \ + psphotIsophotal.c \ + psphotAnnuli.c \ + psphotKron.c \ + psphotKernelFromPSF.c \ + psphotPSFConvModel.c \ + psphotModelTest.c \ + psphotFitSet.c \ + psphotSourceFreePixels.c \ + psphotSummaryPlots.c \ + psphotMergeSources.c \ + psphotLoadPSF.c \ + psphotReadoutCleanup.c \ + psphotSourcePlots.c \ + psphotRadialPlot.c \ + psphotDeblendSatstars.c \ + psphotMosaicSubimage.c \ + psphotMakeResiduals.c \ + psphotSourceSize.c \ + psphotDiagnosticPlots.c \ + psphotMakeFluxScale.c \ + psphotCheckStarDistribution.c \ psphotAddNoise.c Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 17288) +++ trunk/psphot/src/psphot.h (revision 17396) @@ -45,5 +45,6 @@ bool psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background); bool psphotSkyReplace (pmConfig *config, const pmFPAview *view); -bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe); // used by psphotFindDetections @@ -87,5 +88,5 @@ bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psMaskType maskVal); bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal); -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, psMaskType maskVal); +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal); psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psMaskType maskVal); @@ -117,5 +118,6 @@ bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal); -bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal); +pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal, int psfSize); + psKernel *psphotKernelFromPSF (pmSource *source, int nPix); Index: trunk/psphot/src/psphotAnnuli.c =================================================================== --- trunk/psphot/src/psphotAnnuli.c (revision 17288) +++ trunk/psphot/src/psphotAnnuli.c (revision 17396) @@ -1,3 +1,3 @@ -# include "psphot.h" +# include "psphotInternal.h" bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal) { @@ -55,4 +55,6 @@ source->extpars->annuli->fluxVar = fluxSquare; + psFree (pixelCount); + return true; } Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 17288) +++ trunk/psphot/src/psphotBlendFit.c (revision 17396) @@ -40,4 +40,5 @@ // skip non-astronomical objects (very likely defects) if (source->mode & PM_SOURCE_MODE_BLEND) continue; + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; @@ -79,14 +80,19 @@ // XXX re-consider conditions under which the source has EXT fit: // I should try EXT if the source size measurement says it is large - if (psphotFitBlend (readout, source, psf, maskVal)) { - psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); - Npsf ++; - continue; - } - if (psphotFitBlob (readout, source, sources, psf, maskVal)) { - psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y); - Next ++; - continue; - } + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + if (psphotFitBlob (readout, source, sources, psf, maskVal)) { + source->type = PM_SOURCE_TYPE_EXTENDED; + psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y); + Next ++; + continue; + } + } else { + if (psphotFitBlend (readout, source, psf, maskVal)) { + source->type = PM_SOURCE_TYPE_STAR; + psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); + Npsf ++; + continue; + } + } psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y); Index: trunk/psphot/src/psphotEvalFLT.c =================================================================== --- trunk/psphot/src/psphotEvalFLT.c (revision 17288) +++ trunk/psphot/src/psphotEvalFLT.c (revision 17396) @@ -8,4 +8,5 @@ if (model == NULL) { source->mode &= ~PM_SOURCE_MODE_FITTED; + psTrace ("psphot", 5, "no model fitted?\n"); return false; } @@ -14,6 +15,8 @@ if (!(model->flags & PM_MODEL_STATUS_FITTED)) { source->mode &= ~PM_SOURCE_MODE_FITTED; + psTrace ("psphot", 5, "no model fitted?\n"); return false; } + // did the model fit fail for one or another reason? if (model->flags & (PM_MODEL_STATUS_BADARGS | @@ -22,4 +25,15 @@ source->mode |= PM_SOURCE_MODE_FAIL; psLogMsg ("psphot", 5, "EXT fail fit\n"); + psTrace ("psphot", 5, "EXT fail fit\n"); + + if (model->flags & PM_MODEL_STATUS_OFFIMAGE) { + psTrace ("psphot", 5, "off image\n"); + } + if (model->flags & PM_MODEL_STATUS_BADARGS) { + psTrace ("psphot", 5, "bad args\n"); + } + if (model->flags & PM_MODEL_STATUS_NONCONVERGE) { + psTrace ("psphot", 5, "non converge\n"); + } return false; } @@ -35,4 +49,5 @@ if (model->params->data.F32[PM_PAR_I0] <= 0.02) { source->mode |= PM_SOURCE_MODE_FAIL; + psTrace ("psphot", 5, "model central intensity ~ zero\n"); return false; } @@ -43,4 +58,6 @@ // poor-quality fit; only keep if nothing else works... psLogMsg ("psphot", 5, "EXT poor fit\n"); + psTrace ("psphot", 5, "EXT poor fit\n"); + source->mode |= PM_SOURCE_MODE_POOR; return false; Index: trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 17396) +++ trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 17396) @@ -0,0 +1,131 @@ +# include "psphotInternal.h" + +// aperture-like measurements for extended sources +bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + int Next = 0; + int Npetro = 0; + int Nisophot = 0; + int Nannuli = 0; + int Nkron = 0; + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); + return true; + } + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + + // S/N limit to perform full non-linear fits + float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); + + // which extended source analyses should we perform? + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortBySN); + + // XXX some init functions for the extended source recipe options? + + // choose the sources of interest + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip PSF-like and non-astronomical objects + if (source->type == PM_SOURCE_TYPE_STAR) continue; + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // limit selection to some SN limit + assert (source->peak); // how can a source not have a peak? + if (source->peak->SN < SN_LIM) continue; + + // limit selection by analysis region + if (source->peak->x < AnalysisRegion.x0) continue; + if (source->peak->y < AnalysisRegion.y0) continue; + if (source->peak->x > AnalysisRegion.x1) continue; + if (source->peak->y > AnalysisRegion.y1) continue; + + // replace object in image + if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Next ++; + + // if we request any of these measurements, we require the radial profile + if (doPetrosian || doIsophotal || doAnnuli || doKron) { + if (!psphotRadialProfile (source, recipe, maskVal)) { + // all measurements below require the radial profile; skip them all + // re-subtract the object, leave local sky + psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + continue; + } + } + + // Isophotal Mags + if (doIsophotal) { + if (!psphotIsophotal (source, recipe, maskVal)) { + psTrace ("psphot", 5, "failed to measure isophotal mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + } else { + psTrace ("psphot", 5, "measured isophotal mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + Nisophot ++; + } + } + + // Petrosian Mags + if (doPetrosian) { + if (!psphotPetrosian (source, recipe, maskVal)) { + psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + } else { + psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + Npetro ++; + } + } + + // Kron Mags + if (doKron) { + if (!psphotKron (source, recipe, maskVal)) { + psTrace ("psphot", 5, "failed to measure kron mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + } else { + psTrace ("psphot", 5, "measure kron mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + Nkron ++; + } + } + + // Radial Annuli + if (doAnnuli) { + if (!psphotAnnuli (source, recipe, maskVal)) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure in Annuli analysis"); + return false; + } + psTrace ("psphot", 5, "measured annuli for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + Nannuli ++; + } + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + } + + psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next); + psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); + psLogMsg ("psphot", PS_LOG_INFO, " %d isophotal\n", Nisophot); + psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); + psLogMsg ("psphot", PS_LOG_INFO, " %d kron\n", Nkron); + return true; +} Index: trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceFits.c (revision 17396) +++ trunk/psphot/src/psphotExtendedSourceFits.c (revision 17396) @@ -0,0 +1,270 @@ +# include "psphotInternal.h" + +// non-linear model fitting for extended sources +bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + int Next = 0; + int Nconvolve = 0; + int NconvolvePass = 0; + int Nplain = 0; + int NplainPass = 0; + bool savePics = false; + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // perform full extended source non-linear fits? + if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); + return true; + } + + // select the collection of desired models + psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS"); + if (!status) { + psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n"); + return true; + } + if (models->list->n == 0) { + psWarning ("extended source model fits requested but no models are specified\n"); + return true; + } + + // validate the model entries + psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); + psMetadataItem *item = NULL; + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + + if (item->type != PS_DATA_METADATA) { + 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 + } + + psMetadata *model = (psMetadata *) item->data.md; + + // check on the model type + char *modelName = psMetadataLookupStr (&status, model, "MODEL"); + int modelType = pmModelClassGetType (modelName); + if (modelType < 0) { + psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName); + } + psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType); + + // check on the SNLIM, set a float value + char *SNword = psMetadataLookupStr (&status, model, "SNLIM"); + if (!status) { + psAbort("SNLIM not defined for extended source model %s\n", item->name); + } + float SNlim = atof (SNword); + psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim); + + // check on the PSF-Convolution status + char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED"); + if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) { + psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name); + } + bool convolved = !strcasecmp (convolvedWord, "true"); + psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved); + } + psFree (iter); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + + // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + assert (status); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortBySN); + + // choose the sources of interest + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip PSF-like and non-astronomical objects + if (source->type == PM_SOURCE_TYPE_STAR) continue; + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // XXX this should use peak? + if (source->peak->x < AnalysisRegion.x0) continue; + if (source->peak->y < AnalysisRegion.y0) continue; + if (source->peak->x > AnalysisRegion.x1) continue; + if (source->peak->y > AnalysisRegion.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; + + // replace object in image + if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Next ++; + + // save the modelFlux here in case we need to subtract it (for failure) + psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux); + + 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 (4); + + // allocate the array to store the model fits + if (source->modelFits == NULL) { + source->modelFits = psArrayAllocEmpty (4); + } + + // loop here over the models chosen for each source (exclude by S/N) + psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); + psMetadataItem *item = NULL; + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + + // XXX this should have been forced above + assert (item->type == PS_DATA_METADATA); + psMetadata *model = (psMetadata *) item->data.md; + + // check the SNlim and skip model if source is too faint + float SNlim = 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 (source->peak->SN < SNlim) continue; + + // check on the model type + pmModelType modelType = psMetadataLookupS32 (&status, model, "MODEL_TYPE"); + assert (status); + + // check on the PSF-Convolution status + 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 + pmModel *modelFit = NULL; + if (convolved) { + modelFit = psphotPSFConvModel (readout, source, modelType, maskVal, psfSize); + if (!modelFit) { + psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->x, source->moments->y); + continue; + } + psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq); + Nconvolve ++; + if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { + NconvolvePass ++; + } + } else { + psFree (source->modelFlux); + source->modelFlux = NULL; + modelFit = psphotFitEXT (readout, source, modelType, maskVal); + if (!modelFit) { + psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->x, source->moments->y); + continue; + } + pmSourceCacheModel (source, maskVal); + psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq); + Nplain ++; + if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { + NplainPass ++; + } + } + + // save each of the model flux images and store the best + psArrayAdd (modelFluxes, 4, source->modelFlux); + + // test for fit quality / result + psArrayAdd (source->modelFits, 4, modelFit); + + psFree (modelFit); + } + psFree (iter); + + // evaluate the relative quality of the models, choose one + float minChisq = NAN; + int minModel = -1; + for (int i = 0; i < source->modelFits->n; i++) { + pmModel *model = source->modelFits->data[i]; + + if (!(model->flags & PM_MODEL_STATUS_FITTED)) continue; + + if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue; + if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue; + if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue; + + if ((minModel < 0) || (model->chisq < minChisq)) { + minChisq = model->chisq; + minModel = i; + } + } + + if (minModel == -1) { + // re-subtract the object, leave local sky + psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->x, source->moments->y); + + // replace original model, subtract it + psFree (source->modelFlux); + source->modelFlux = modelFluxStart; + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + + psFree (modelFluxes); + + if (savePics) { + psphotSaveImage (NULL, readout->image, "image.xp.fits"); + char key[10]; + fprintf (stdout, "continue? "); + fgets (key, 8, stdin); + if (key[0] == 'n') { + savePics = false; + } + } + continue; + } + + // save the best extended model in modelEXT + psFree (source->modelEXT); + source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]); + + // save the modelFlux for the best model + psFree (source->modelFlux); + source->modelFlux = psMemIncrRefCounter (modelFluxes->data[minModel]); + + // subtract the best fit from the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + + // 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->x, source->moments->y, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq); + psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + + if (savePics) { + psphotSaveImage (NULL, readout->image, "image.xm.fits"); + char key[10]; + fprintf (stdout, "continue? "); + fgets (key, 8, stdin); + if (key[0] == 'n') { + savePics = false; + } + } + } + + psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), 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); + return true; +} Index: trunk/psphot/src/psphotFindDetections.c =================================================================== --- trunk/psphot/src/psphotFindDetections.c (revision 17288) +++ trunk/psphot/src/psphotFindDetections.c (revision 17396) @@ -8,4 +8,6 @@ bool status; int pass; + + psTimerStart ("psphot"); // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -67,5 +69,5 @@ psphotCullPeaks(readout->image, readout->weight, recipe, detections->footprints); detections->peaks = pmFootprintArrayToPeaks(detections->footprints); - psLogMsg ("psphot", PS_LOG_INFO, "peaks: %ld, total footprints: %ld\n", detections->peaks->n, detections->footprints->n); + psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot")); return detections; Index: trunk/psphot/src/psphotFindPeaks.c =================================================================== --- trunk/psphot/src/psphotFindPeaks.c (revision 17288) +++ trunk/psphot/src/psphotFindPeaks.c (revision 17396) @@ -11,5 +11,5 @@ // smooth the image and weight map - psTimerStart ("psphot"); + psTimerStart ("peaks"); // XXX if we have been supplied a PSF, we can use that to set the FWHM_X,FWHM_Y values @@ -32,10 +32,10 @@ psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32); psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH); - psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("peaks")); // smooth the weight, applying the mask as we go psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32); psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH); - psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("peaks")); psImage *mask = readout->mask; @@ -61,5 +61,5 @@ } } - psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("peaks")); // optionally save example images under trace @@ -70,5 +70,5 @@ } - psTimerStart ("psphot"); + psTimerStart ("peaks"); // set peak threshold @@ -140,5 +140,5 @@ pmPeaksWriteText (peaks, output); } - psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot")); + psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("peaks")); // If they asked us to return a set of pmFootprints, not a raw pmPeak list, @@ -167,5 +167,5 @@ peaks = footprints; // well, you know what I mean - psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("psphot")); + psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("peaks")); } Index: trunk/psphot/src/psphotIsophotal.c =================================================================== --- trunk/psphot/src/psphotIsophotal.c (revision 17288) +++ trunk/psphot/src/psphotIsophotal.c (revision 17396) @@ -1,3 +1,5 @@ -# include "psphot.h" +# include "psphotInternal.h" + +static float ISOPHOT_FLUX = NAN; bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal) { @@ -14,9 +16,12 @@ // flux at which to measure isophotal parameters - // XXX cache this? - float ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX"); + if (!isfinite(ISOPHOT_FLUX)) { + // XXX ISOPHOTAL_FLUX should be specified in mags, need the zero point to get counts/sec + ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX"); + assert (status); + } // find the first bin below the flux level and the last above the level - // XXX can this be done faster with disection? + // XXX can this be done faster with bisection? // XXX do I need to worry about crazy outliers? // XXX should i be smoothing or fitting the curve? @@ -27,5 +32,17 @@ if ((firstBelow < 0) && (flux->data.F32[i] < ISOPHOT_FLUX)) firstBelow = i; } - + // if we don't go out far enough, we have a problem... + if (lastAbove == flux->n - 1) { + psTrace ("psphot", 5, "did not go out far enough to reach isophotal magnitude"); + // XXX raise a flag ? + return false; + } + if (firstBelow < 0) { + psTrace ("psphot", 5, "did not go out far enough to bound isophotal magnitude: error unmeasured"); + // XXX raise a flag ? + lastAbove = firstBelow; + return false; + } + // need to examine pixels in this vicinity float isophotalFluxFirst = 0; @@ -56,4 +73,9 @@ source->extpars->isophot->radErr = isophotalRadErr; + psTrace ("psphot", 5, "Isophot flux:%f +/- %f @ %f +/- %f for %f, %f\n", + source->extpars->isophot->mag, source->extpars->isophot->magErr, + source->extpars->isophot->rad, source->extpars->isophot->radErr, + source->peak->xf, source->peak->yf); + return true; Index: trunk/psphot/src/psphotKernelFromPSF.c =================================================================== --- trunk/psphot/src/psphotKernelFromPSF.c (revision 17288) +++ trunk/psphot/src/psphotKernelFromPSF.c (revision 17396) @@ -1,20 +1,45 @@ -# include "psphot.h" +# include "psphotInternal.h" psKernel *psphotKernelFromPSF (pmSource *source, int nPix) { - assert (source); - assert (source->psfFlux); // XXX build if needed? + assert (source); + assert (source->psfFlux); // XXX build if needed? - int x0 = source->peak->xf - source->psfFlux->col0; - int y0 = source->peak->yf - source->psfFlux->row0; + int x0 = source->peak->xf - source->psfFlux->col0; + int y0 = source->peak->yf - source->psfFlux->row0; - // need to decide on the size: dynamically? statically? - psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix); + // need to decide on the size: dynamically? statically? + psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix); - for (int j = psf->yMin; j <= psf->yMax; j++) { - for (int i = psf->xMin; i <= psf->xMax; i++) { - psf->kernel[j][i] = source->psfFlux->data.F32[y0 + j][x0 + i]; + // XXX we should just re-construct a PSF at this location + // psModelAdd (psf->image, NULL, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM | PM_MODEL_OP_CENTER); + + // if the realized PSF for this object does not cover the full kernel, give up for now + if (x0 + psf->xMin < 0) goto escape; + if (x0 + psf->xMax >= source->psfFlux->numCols) goto escape; + if (y0 + psf->yMin < 0) goto escape; + if (y0 + psf->yMax >= source->psfFlux->numRows) goto escape; + + double sum = 0.0; + for (int j = psf->yMin; j <= psf->yMax; j++) { + for (int i = psf->xMin; i <= psf->xMax; i++) { + double value = source->psfFlux->data.F32[y0 + j][x0 + i]; + psf->kernel[j][i] = value; + sum += value; + } } - } - return psf; + assert (sum > 0.0); + + // psf must be normalized (integral = 1.0) + for (int i = 0; i < psf->image->numRows; i++) { + for (int j = 0; j < psf->image->numCols; j++) { + psf->image->data.F32[i][j] /= sum; + } + } + + return psf; + +escape: + psFree (psf); + return NULL; } Index: trunk/psphot/src/psphotKron.c =================================================================== --- trunk/psphot/src/psphotKron.c (revision 17288) +++ trunk/psphot/src/psphotKron.c (revision 17396) @@ -1,3 +1,3 @@ -# include "psphot.h" +# include "psphotInternal.h" bool psphotKron (pmSource *source, psMetadata *recipe, psMaskType maskVal) { Index: trunk/psphot/src/psphotModelTest.c =================================================================== --- trunk/psphot/src/psphotModelTest.c (revision 17288) +++ trunk/psphot/src/psphotModelTest.c (revision 17396) @@ -208,6 +208,10 @@ source->modelPSF = pmModelFromPSF (model, psf); source->modelEXT = model; - status = psphotPSFConvModel (source, recipe, maskVal); - model = source->modelConv; + + // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + assert (status); + + model = psphotPSFConvModel (readout, source, modelType, maskVal, psfSize); params = model->params->data.F32; } else { Index: trunk/psphot/src/psphotModelWithPSF.c =================================================================== --- trunk/psphot/src/psphotModelWithPSF.c (revision 17288) +++ trunk/psphot/src/psphotModelWithPSF.c (revision 17396) @@ -1,3 +1,4 @@ -# include "psphot.h" +# include "psphotInternal.h" +# define SAVE_IMAGES 0 bool psphotModelWithPSF_LMM ( @@ -36,8 +37,15 @@ // allocate internal arrays (current vs Guess) - psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32); - psImage *Alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32); - psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32); - psVector *Beta = psVectorAlloc(params->n, PS_TYPE_F32); + psImage *Alpha = NULL; + psVector *Beta = NULL; + + // Alpha & Beta only contain elements to represent the unmasked parameters + if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) { + psAbort ("programming error: no unmasked parameters to be fit\n"); + } + + // allocate internal arrays (current vs Guess) + psImage *alpha = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32); + psVector *beta = psVectorAlloc(Beta->n, PS_TYPE_F32); psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32); @@ -79,6 +87,7 @@ // dump some useful info if trace is defined if (psTraceGetLevel("psphot") >= 6) { - p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)"); - p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)"); + p_psImagePrint(psTraceGetDestination(), Alpha, "Alpha guess (1)"); + p_psVectorPrint(psTraceGetDestination(), Beta, "Beta guess (1)"); + p_psVectorPrint(psTraceGetDestination(), beta, "beta current (1)"); } if (psTraceGetLevel("psphot") >= 5) { @@ -116,5 +125,5 @@ beta = psVectorCopy(beta, Beta, PS_TYPE_F32); params = psVectorCopy(params, Params, PS_TYPE_F32); - lambda *= 0.1; + lambda *= 0.25; // save the new convolved model image @@ -130,7 +139,21 @@ // construct & return the covariance matrix (if requested) if (covar != NULL) { - if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) { + if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) { psTrace ("psphot", 5, "failure to calculate covariance matrix\n"); } + // set covar values which are not masked + psImageInit (covar, 0.0); + for (int j = 0, J = 0; j < params->n; j++) { + if (paramMask && (paramMask->data.U8[j])) { + covar->data.F32[j][j] = 1.0; + continue; + } + for (int k = 0, K = 0; k < params->n; k++) { + if (paramMask && (paramMask->data.U8[k])) continue; + covar->data.F32[j][k] = Alpha->data.F32[J][K]; + K++; + } + J++; + } } @@ -177,5 +200,5 @@ } - // generate the model and derivative images for this parameter set + // 1 *** generate the model and derivative images for this parameter set // storage for model derivatives @@ -223,7 +246,7 @@ for (int n = 0; n < params->n; n++) { - if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } - psImage *dmodel = pcm->dmodels->data[n]; - dmodel->data.F32[i][j] = deriv->data.F32[n]; + if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } + psImage *dmodel = pcm->dmodels->data[n]; + dmodel->data.F32[i][j] = deriv->data.F32[n]; } } @@ -231,17 +254,16 @@ psFree(coord); psFree(deriv); - // convolve model and dmodel arrays with PSF psImageConvolveDirect (pcm->modelConv, pcm->model, psf); for (int n = 0; n < pcm->dmodels->n; n++) { - if (pcm->dmodels->data[n] == NULL) continue; - psImage *dmodel = pcm->dmodels->data[n]; - psImage *dmodelConv = pcm->dmodelsConv->data[n]; - psImageConvolveDirect (dmodelConv, dmodel, psf); + if (pcm->dmodels->data[n] == NULL) continue; + psImage *dmodel = pcm->dmodels->data[n]; + psImage *dmodelConv = pcm->dmodelsConv->data[n]; + psImageConvolveDirect (dmodelConv, dmodel, psf); } // XXX TEST : SAVE IMAGES -# if (1) +# if (SAVE_IMAGES) psphotSaveImage (NULL, psf->image, "psf.fits"); psphotSaveImage (NULL, pcm->model, "model.fits"); @@ -252,4 +274,6 @@ # endif + // 2 *** accumulate alpha & beta + // zero alpha and beta for summing below psImageInit (alpha, 0.0); @@ -259,5 +283,5 @@ for (psS32 i = 0; i < source->pixels->numRows; i++) { for (psS32 j = 0; j < source->pixels->numCols; j++) { - // XXX are we doing the right thing with the mask? + // XXX are we doing the right thing with the mask? // skip masked points if (source->maskObj->data.U8[i][j]) { @@ -279,23 +303,20 @@ chisq += PS_SQR(delta) * yweight; - if (isnan(delta)) - psAbort("nan in delta"); - if (isnan(chisq)) - psAbort("nan in chisq"); - - for (psS32 n1 = 0; n1 < params->n; n1++) { - if ((paramMask != NULL) && (paramMask->data.U8[n1])) { - continue; - } - psImage *dmodel = pcm->dmodelsConv->data[n1]; - float weight = dmodel->data.F32[i][j] * yweight; - for (psS32 n2 = 0; n2 <= n1; n2++) { - if ((paramMask != NULL) && (paramMask->data.U8[n2])) { - continue; - } - dmodel = pcm->dmodelsConv->data[n2]; - alpha->data.F32[n1][n2] += weight * dmodel->data.F32[i][j]; - } - beta->data.F32[n1] += weight * delta; + if (isnan(delta)) psAbort("nan in delta"); + if (isnan(chisq)) psAbort("nan in chisq"); + + // alpha & beta only contain unmasked elements + for (int n1 = 0, N1 = 0; n1 < params->n; n1++) { + if ((paramMask != NULL) && (paramMask->data.U8[n1])) continue; + psImage *dmodel = pcm->dmodelsConv->data[n1]; + float weight = dmodel->data.F32[i][j] * yweight; + for (int n2 = 0, N2 = 0; n2 <= n1; n2++) { + if ((paramMask != NULL) && (paramMask->data.U8[n2])) continue; + dmodel = pcm->dmodelsConv->data[n2]; + alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j]; + N2++; + } + beta->data.F32[N1] += weight * delta; + N1++; } } @@ -303,17 +324,7 @@ // calculate lower-left half of alpha - for (psS32 j = 1; j < params->n; j++) { + for (psS32 j = 1; j < alpha->numCols; j++) { for (psS32 k = 0; k < j; k++) { alpha->data.F32[k][j] = alpha->data.F32[j][k]; - } - } - - // fill in pivots if we apply a mask - if (paramMask != NULL) { - for (psS32 j = 0; j < params->n; j++) { - if (paramMask->data.U8[j]) { - alpha->data.F32[j][j] = 1; - beta->data.F32[j] = 1; - } } } @@ -345,6 +356,7 @@ pcm->dmodels = psArrayAlloc (params->n); for (psS32 n = 0; n < params->n; n++) { - if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } - pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + pcm->dmodels->data[n] = NULL; + if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } + pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); } @@ -353,6 +365,7 @@ pcm->dmodelsConv = psArrayAlloc (params->n); for (psS32 n = 0; n < params->n; n++) { - if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } - pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + pcm->dmodelsConv->data[n] = NULL; + if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } + pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); } @@ -378,8 +391,8 @@ while () { - GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda) - dLinear = dLinear(Beta, beta, lambda); - chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func) - convergence tests... + GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda) + dLinear = dLinear(Beta, beta, lambda); + chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func) + convergence tests... } @@ -388,16 +401,16 @@ ** GuessABP: - f_c = sum_i (kern_i * func (x_i; p_o)) - -df_c/dp_o = d/dp_o [sum_i (kern_i * func (x_i; p_o))] - -df_c/dp_o = sum_i (d/dp_o [kern_i * func (x_i; p_o)]) - -df_c/dp_o = sum_i (kern_i * d/dp_o [func (x_i; p_o)]) - -- generate image arrays for func, dfunc/dp_j (not masked) -- convolve each with psf -- measure delta = f_conv - data -- etc + f_c = sum_i (kern_i * func (x_i; p_o)) + + df_c/dp_o = d/dp_o [sum_i (kern_i * func (x_i; p_o))] + + df_c/dp_o = sum_i (d/dp_o [kern_i * func (x_i; p_o)]) + + df_c/dp_o = sum_i (kern_i * d/dp_o [func (x_i; p_o)]) + + - generate image arrays for func, dfunc/dp_j (not masked) + - convolve each with psf + - measure delta = f_conv - data + - etc */ Index: trunk/psphot/src/psphotPSFConvModel.c =================================================================== --- trunk/psphot/src/psphotPSFConvModel.c (revision 17288) +++ trunk/psphot/src/psphotPSFConvModel.c (revision 17396) @@ -1,52 +1,24 @@ -# include "psphot.h" +# include "psphotInternal.h" +# define USE_DELTA_PSF 0 // save as static values so they may be set externally static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; -// static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0; -// static bool PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true; // input source has both modelPSF and modelEXT. on successful exit, we set the // modelConv to contain the fitted parameters, and the modelFlux to contain the // convolved model image. -bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal) { +pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal, int psfSize) { - bool status; - - int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); - if (!status) { - psfSize = 2; - } - // make sure we save a cached copy of the psf flux pmSourceCachePSF (source, maskVal); // convert the cached cached psf model for this source to a psKernel - // XXX for the moment, hard-wire the kernel to be 5x5 (2 pix radius) - // XXX for the moment, hard-wire the kernel to be 9x9 (4 pix radius) psKernel *psf = psphotKernelFromPSF (source, psfSize); + if (!psf) return NULL; - // psf must be normalized (integral = 1.0) - double sum = 0.0; - for (int i = 0; i < psf->image->numRows; i++) { - for (int j = 0; j < psf->image->numCols; j++) { - sum += psf->image->data.F32[i][j]; - } - } - assert (sum > 0.0); - for (int i = 0; i < psf->image->numRows; i++) { - for (int j = 0; j < psf->image->numCols; j++) { - psf->image->data.F32[i][j] /= sum; - } - } - -# if (0) - // XXX sanity check: convolve with delta function should behave like unconvolved version - for (int i = 0; i < psf->image->numRows; i++) { - for (int j = 0; j < psf->image->numCols; j++) { - psf->image->data.F32[i][j] = 0.0; - } - } - psf->image->data.F32[2][2] = 1.0; +# if (USE_DELTA_PSF) + psImageInit (psf->image, 0.0); + psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0; # endif @@ -54,7 +26,44 @@ // XXX we could modify the parameter values or even the model // here based on the observed seeing (some lookup table...) - pmModel *modelConv = pmModelCopy (source->modelEXT); + + // use the source moments, etc to guess basic model parameters + pmModel *modelConv = pmSourceModelGuess (source, modelType); + if (!modelConv) { + psFree (psf); + return NULL; + } + + // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: + psEllipseShape psfShape; + psfShape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; + psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; + psfShape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; + psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0); + + // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: + psEllipseShape extShape; + extShape.sx = modelConv->params->data.F32[PM_PAR_SXX] / M_SQRT2; + extShape.sxy = modelConv->params->data.F32[PM_PAR_SXY]; + extShape.sy = modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2; + psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0); + + // decrease the initial guess ellipse by psf_minor axis: + psEllipseAxes extAxesMod; + extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor))); + extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor))); + extAxesMod.theta = extAxes.theta; + + psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod); + modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2; + modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy; + modelConv->params->data.F32[PM_PAR_SYY] = extShapeMod.sy * M_SQRT2; + + // increase the initial guess central intensity by 2pi r^2: + modelConv->params->data.F32[PM_PAR_I0] *= (1.0 + PS_SQR(psfAxes.minor) / PS_SQR(extAxesMod.minor)); + psVector *params = modelConv->params; psVector *dparams = modelConv->dparams; + + psphotCheckRadiusEXT (readout, source, modelConv); // create the minimization constraints @@ -91,5 +100,5 @@ psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); - // renormalize output model image + // renormalize output model image (generated by fitting process) float Io = params->data.F32[PM_PAR_I0]; for (int iy = 0; iy < source->modelFlux->numRows; iy++) { @@ -115,17 +124,13 @@ onPic &= (params->data.F32[PM_PAR_YPOS] >= source->pixels->row0); onPic &= (params->data.F32[PM_PAR_YPOS] < source->pixels->row0 + source->pixels->numRows); - if (!onPic) { - modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE; - } + if (!onPic) modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE; - source->mode |= PM_SOURCE_MODE_FITTED; - source->modelConv = modelConv; + source->mode |= PM_SOURCE_MODE_FITTED; // XXX is this needed? + psFree(psf); psFree(myMin); psFree(covar); psFree(constraint); - bool retval = (onPic && fitStatus); - psTrace("psphot", 5, "---- %s(%d) end ----\n", __func__, retval); - return(retval); + return modelConv; } Index: trunk/psphot/src/psphotPetrosian.c =================================================================== --- trunk/psphot/src/psphotPetrosian.c (revision 17288) +++ trunk/psphot/src/psphotPetrosian.c (revision 17396) @@ -1,3 +1,6 @@ -# include "psphot.h" +# include "psphotInternal.h" + +static float PETROSIAN_R0 = NAN; +static float PETROSIAN_RF = NAN; bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) { @@ -14,14 +17,27 @@ // flux at which to measure isophotal parameters - // XXX cache this? - float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0"); - float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO"); + if (!isfinite(PETROSIAN_R0)) { + PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0"); + PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO"); + assert (status); + } // first find flux at R0 - int firstBelow = -1; - int lastAbove = -1; + int firstAbove = -1; + int lastBelow = -1; for (int i = 0; i < radius->n; i++) { - if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i; - if ((firstBelow < 0) && (flux->data.F32[i] < PETROSIAN_R0)) firstBelow = i; + if (radius->data.F32[i] < PETROSIAN_R0) lastBelow = i; + if ((firstAbove < 0) && (radius->data.F32[i] > PETROSIAN_R0)) firstAbove = i; + } + // if we don't go out far enough, we have a problem... + if (lastBelow == radius->n - 1) { + psTrace ("psphot", 5, "did not go out far enough to reach petrosian reference radius..."); + // XXX skip object? raise a flag ? + return false; + } + if (firstAbove < 0) { + psTrace ("psphot", 5, "did not go out far enough to bound petrosian reference radius"); + // XXX raise a flag ? + return false; } @@ -29,5 +45,5 @@ float fluxR0 = 0.0; int fluxRn = 0; - for (int i = PS_MIN(firstBelow, lastAbove); i <= PS_MAX(firstBelow, lastAbove); i++) { + for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) { fluxR0 += flux->data.F32[i]; fluxRn ++; @@ -42,9 +58,20 @@ // XXX do I need to worry about crazy outliers? // XXX should i be smoothing or fitting the curve? - firstBelow = -1; - lastAbove = -1; + int firstBelow = -1; + int lastAbove = -1; for (int i = 0; i < flux->n; i++) { if (flux->data.F32[i] > fluxRP) lastAbove = i; if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i; + } + // if we don't go out far enough, we have a problem... + if (lastAbove == radius->n - 1) { + psTrace ("psphot", 5, "did not go out far enough to reach petrosian radius..."); + // XXX skip object? raise a flag ? + return false; + } + if (firstBelow < 0) { + psTrace ("psphot", 5, "did not go out far enough to bound petrosian radius"); + // XXX raise a flag ? + return false; } @@ -78,4 +105,9 @@ source->extpars->petrosian->radErr = radErr; + psTrace ("psphot", 5, "Petrosian flux:%f +/- %f @ %f +/- %f for %f, %f\n", + source->extpars->petrosian->mag, source->extpars->petrosian->magErr, + source->extpars->petrosian->rad, source->extpars->petrosian->radErr, + source->peak->xf, source->peak->yf); + return true; Index: trunk/psphot/src/psphotRadialProfile.c =================================================================== --- trunk/psphot/src/psphotRadialProfile.c (revision 17288) +++ trunk/psphot/src/psphotRadialProfile.c (revision 17396) @@ -41,6 +41,15 @@ int n = 0; - float Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; - float Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; + + float Xo = 0.0; + float Yo = 0.0; + + if (source->modelEXT) { + Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; + Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; + } else { + Xo = source->peak->xf - source->pixels->col0; + Yo = source->peak->yf - source->pixels->row0; + } for (int iy = 0; iy < source->pixels->numRows; iy++) { for (int ix = 0; ix < source->pixels->numCols; ix++) { Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 17288) +++ trunk/psphot/src/psphotReadout.c (revision 17396) @@ -117,6 +117,5 @@ psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); - // XXX test the CR/EXT measurement - // XXX we need to call this here and after the second pass: option for starting number? + // identify CRs and extended sources psphotSourceSize (readout, sources, recipe, 0); if (!strcasecmp (breakPt, "ENSEMBLE")) { @@ -127,11 +126,13 @@ psphotBlendFit (readout, sources, recipe, psf); - // replace all sources (make this part of psphotFitSourcesLinear?) + // replace all sources psphotReplaceAll (sources, recipe); // linear fit to include all sources psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); + + // if we only do one pass, skip to extended source analysis if (!strcasecmp (breakPt, "PASS1")) { - goto finish; + goto pass1finish; } @@ -145,7 +146,8 @@ psphotAddNoise (readout, sources, recipe); + // find fainter sources (pass 2) detections = psphotFindDetections (detections, readout, recipe); - // remove noise for subtracted objects + // remove noise for subtracted objects (ie, return to normal noise level) psphotSubNoise (readout, sources, recipe); @@ -172,8 +174,12 @@ psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); +pass1finish: + // measure source size for the remaining sources psphotSourceSize (readout, sources, recipe, 0); - psphotExtendedSources (readout, sources, recipe); + psphotExtendedSourceAnalysis (readout, sources, recipe); + + psphotExtendedSourceFits (readout, sources, recipe); finish: Index: trunk/psphot/src/psphotSourceFits.c =================================================================== --- trunk/psphot/src/psphotSourceFits.c (revision 17288) +++ trunk/psphot/src/psphotSourceFits.c (revision 17396) @@ -220,5 +220,5 @@ pmSource *tmpSrc = pmSourceAlloc (); - pmModel *EXT = psphotFitEXT (readout, source, maskVal); + pmModel *EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal); okEXT = psphotEvalEXT (tmpSrc, EXT); chiEXT = EXT->chisq / EXT->nDOF; @@ -267,4 +267,6 @@ // save new model source->modelEXT = EXT; + source->type = PM_SOURCE_TYPE_EXTENDED; + source->mode |= PM_SOURCE_MODE_EXTMODEL; // build cached model and subtract @@ -349,10 +351,10 @@ } -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, psMaskType maskVal) { +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal) { NfitEXT ++; // use the source moments, etc to guess basic model parameters - pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); + pmModel *EXT = pmSourceModelGuess (source, modelType); PS_ASSERT (EXT, NULL); Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 17288) +++ trunk/psphot/src/psphotSourceSize.c (revision 17396) @@ -1,4 +1,6 @@ # include "psphotInternal.h" # include + +float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro); // we need to call this function after sources have been fitted to the PSF model and @@ -14,5 +16,8 @@ psTimerStart ("psphot"); - float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CRNSIGMA.LIMIT"); + float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT"); + assert (status); + + float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT"); assert (status); @@ -24,20 +29,24 @@ if (isfinite(source->crNsigma)) continue; - source->crNsigma = -1.0; - source->extNsigma = -1.0; - // source must have been subtracted - source->crNsigma = -3.0; if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; - psF32 **resid = source->pixels->data.F32; + psF32 **resid = source->pixels->data.F32; psF32 **weight = source->weight->data.F32; - psU8 **mask = source->maskObj->data.U8; + psU8 **mask = source->maskObj->data.U8; + + // check for extendedness: measure the delta flux significance at the 1 sigma contour + source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, source->modelPSF, 1.0); + + // XXX prevent a source from being both CR and EXT? + if (source->extNsigma > EXT_NSIGMA_LIMIT) { + source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + } int xPeak = source->peak->xf - source->pixels->col0 + 0.5; int yPeak = source->peak->yf - source->pixels->row0 + 0.5; - // skip sources which are too close to a boundary - source->crNsigma = -4.0; + // XXX for now, skip sources which are too close to a boundary + // XXX raise a flag? if (xPeak < 1) continue; if (xPeak > source->pixels->numCols - 2) continue; @@ -46,5 +55,5 @@ // XXX for now, just skip any sources with masked pixels - source->crNsigma = -5.0; + // XXX raise a flag? bool keep = true; for (int iy = -1; (iy <= +1) && keep; iy++) { @@ -113,5 +122,4 @@ } source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; - source->extNsigma = (nEXT > 0) ? fabs(fEXT) / nEXT : 0.0; // NOTE: abs needed to make the Nsigma value positive @@ -122,5 +130,5 @@ if (source->crNsigma > CR_NSIGMA_LIMIT) { - source->mode |= PM_SOURCE_MODE_CRLIMIT; + source->mode |= PM_SOURCE_MODE_CR_LIMIT; } } @@ -169,2 +177,65 @@ */ + +// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total +// deviation in sigmas. This is measure on the residual image - should we ignore negative +// deviations? +float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro) { + + psF32 *PAR = model->params->data.F32; + + // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY; + // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0; + // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]; + // y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A]; + + // min/max value of x is where T -> 0 + // solve this for x2: + float Q = Ro * PS_SQR(PAR[PM_PAR_SXX]) / (1.0 - PS_SQR(PAR[PM_PAR_SXX]*PAR[PM_PAR_SYY]*PAR[PM_PAR_SXY]) / 4.0); + if (Q < 0.0) return NAN; // ellipse is imaginary + + int xMax = sqrt(Q); + int xMin = -1.0*xMax; + + int nPts = 0; + float nSigma = 0.0; + + for (int x = xMin; x <= xMax; x++) { + float A = PS_SQR (1.0 / PAR[PM_PAR_SYY]); + float B = x * PAR[PM_PAR_SXY]; + float C = PS_SQR (x / PAR[PM_PAR_SXX]) - Ro; + + float T = PS_SQR(B) - 4*A*C; + if (T < 0.0) continue; + + float yP = (-B + sqrt (T)) / (2.0 * A); + float yM = (-B - sqrt (T)) / (2.0 * A); + + int xPix = x + PAR[PM_PAR_XPOS] - image->col0 + 0.5; + int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5; + int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5; + + if (xPix < 0) continue; + if (xPix >= image->numCols) continue; + + if ((yPixM >= 0) && (yPixM < image->numRows)) { + if (!mask || !mask->data.U8[yPixM][xPix]) { + float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]); + nSigma += dSigma; + nPts ++; + } + } + + if (yPixM == yPixP) continue; + + if ((yPixP >= 0) && (yPixP < image->numRows)) { + if (!mask || !mask->data.U8[yPixP][xPix]) { + float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]); + nSigma += dSigma; + nPts ++; + } + } + } + nSigma /= nPts; + return nSigma; +} Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 17288) +++ trunk/psphot/src/psphotSourceStats.c (revision 17396) @@ -37,11 +37,15 @@ for (int i = 0; i < peaks->n; i++) { + pmPeak *peak = peaks->data[i]; + if (peak->assigned) continue; + // create a new source, add peak pmSource *source = pmSourceAlloc(); - source->peak = (pmPeak *)psMemIncrRefCounter(peaks->data[i]); + source->peak = psMemIncrRefCounter(peak); // allocate image, weight, mask arrays for each peak (square of radius OUTER) pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); if (!strcasecmp (breakPt, "PEAKS")) { + peak->assigned = true; psArrayAdd (sources, 100, source); psFree (source); @@ -49,6 +53,7 @@ } - // skip faint sources + // skip faint sources for moments measurement if (source->peak->SN < MIN_SN) { + peak->assigned = true; psArrayAdd (sources, 100, source); psFree (source); @@ -57,5 +62,5 @@ // measure a local sky value - // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken + // the local sky is now ignored; kept here for reference only status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, mark); if (!status) { @@ -80,4 +85,5 @@ if (status) { // add to the source array + peak->assigned = true; psArrayAdd (sources, 100, source); psFree (source); @@ -93,4 +99,5 @@ if (status) { // add to the source array + peak->assigned = true; psArrayAdd (sources, 100, source); psFree (source);