Index: /branches/eam_branch_20080324/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- /branches/eam_branch_20080324/psphot/src/psphotExtendedSourceAnalysis.c (revision 17276) +++ /branches/eam_branch_20080324/psphot/src/psphotExtendedSourceAnalysis.c (revision 17276) @@ -0,0 +1,131 @@ +# include "psphot.h" + +// aperture-like measurements for extended sources +bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + int Next = 0; + int Nconvolve = 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)) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate radial profile"); + return false; + } + } + + // XXX probably should not give up if we fail on any of the measurements unless it is due to + // a programming error. + + // Isophotal Mags + if (doIsophotal) { + if (!psphotIsophotal (source, recipe, maskVal)) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure in Isophotal analysis"); + return false; + } + psTrace ("psphot", 5, "measure 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)) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure in Petrosian analysis"); + return false; + } + 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)) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure in Kron analysis"); + return false; + } + 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: /branches/eam_branch_20080324/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /branches/eam_branch_20080324/psphot/src/psphotExtendedSourceFits.c (revision 17276) +++ /branches/eam_branch_20080324/psphot/src/psphotExtendedSourceFits.c (revision 17276) @@ -0,0 +1,170 @@ +# include "psphot.h" + +// non-linear model fitting for extended sources +bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + int Next = 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_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; + + // check on the model type + char *modelName = psMetadataLookupStr (&status, model, "MODEL"); + int modelType = psModelClassGetType (modelName); + if (modelType < 0) { + psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, 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, modelType); + } + bool convolved = !strcasecmp (convolvedWord, "true"); + psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved); + } + + // 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"); + + // 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; + + // limit selection to some SN limit + // XXX move this inside model loop + assert (source->peak); // how can a source not have a peak? + if (source->peak->SN < SN_LIM) 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 ++; + + // 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; + + // check the SNlim and skip model if source is too faint + float SNlim = psMetadataLookupF32 (&status, model, "SNLIM_VALUE"); + assert (status); + + if (X < SNlim) { + continue; + } + + // check on the model type + char *modelName = psMetadataLookupStr (&status, model, "MODEL"); + assert (status); + + // check on the PSF-Convolution status + bool convolved = psMetadataLookupBool (&status, model, "PSF_CONVOLVED_VALUE"); + assert (status); + + // fit the model as convolved or not + // XXX use the values above, rework this function + if (convolved) { + pmModel *EXT = psphotPSFConvModel (source, recipe, maskVal); + } else { + pmModel *EXT = psphotFitExt (readout, source, maskVal); + } + if (!EXT) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure to fit extended model"); + return false; + } + + // XXX do something with the resulting fit + + + + } + + // evaluate the relative quality of the models, choose one + + // cache only the chosen model + pmSourceCacheModel (source, maskVal); // XXX put this in the source model function? + psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + Nconvolve ++; + + // 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 convolved models\n", Nconvolve); + 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; +}