Index: /branches/eam_branches/psphot.20100506/src/Makefile.am =================================================================== --- /branches/eam_branches/psphot.20100506/src/Makefile.am (revision 27938) +++ /branches/eam_branches/psphot.20100506/src/Makefile.am (revision 27939) @@ -164,4 +164,5 @@ psphotModelWithPSF.c \ psphotExtendedSourceAnalysis.c \ + psphotExtendedSourceAnalysisByObject.c \ psphotExtendedSourceFits.c \ psphotKernelFromPSF.c \ @@ -192,4 +193,5 @@ psphotRadialBins.c \ psphotRadialApertures.c \ + psphotRadialAperturesByObject.c \ psphotPetrosian.c \ psphotPetrosianRadialBins.c \ Index: /branches/eam_branches/psphot.20100506/src/psphot.h =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphot.h (revision 27938) +++ /branches/eam_branches/psphot.20100506/src/psphot.h (revision 27939) @@ -462,3 +462,6 @@ bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax); +bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule); +bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule); + #endif Index: /branches/eam_branches/psphot.20100506/src/psphotExtendedSourceAnalysis.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotExtendedSourceAnalysis.c (revision 27938) +++ /branches/eam_branches/psphot.20100506/src/psphotExtendedSourceAnalysis.c (revision 27939) @@ -1,3 +1,9 @@ # include "psphotInternal.h" + +// ?? these cannot happen here --> we would need to do this in psphotExtendedSourceAnalysis +// XXX option to choose a consistent position +// XXX option to choose a consistent elliptical contour +// XXX SDSS uses the r-band petrosian radius to measure petrosian fluxes in all bands +// XXX consistent choice of extendedness... // for now, let's store the detections on the readout->analysis for each readout Index: /branches/eam_branches/psphot.20100506/src/psphotExtendedSourceAnalysisByObject.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotExtendedSourceAnalysisByObject.c (revision 27939) +++ /branches/eam_branches/psphot.20100506/src/psphotExtendedSourceAnalysisByObject.c (revision 27939) @@ -0,0 +1,164 @@ +# include "psphotInternal.h" + +// ?? these cannot happen here --> we would need to do this in psphotExtendedSourceAnalysis +// XXX option to choose a consistent position +// XXX option to choose a consistent elliptical contour +// XXX SDSS uses the r-band petrosian radius to measure petrosian fluxes in all bands +// XXX consistent choice of extendedness... + +// aperture-like measurements for extended sources +bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule) { + + bool status; + int Next = 0; + int Npetro = 0; + int Nannuli = 0; + + psTimerStart ("psphot.extended"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // 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; + } + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // XXX require petrosian analysis for non-linear fits? + + // XXX temporary user-supplied systematic sky noise measurement (derive from background model) + float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + + // 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"); + + // number of images used to define sources + int nImages = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // generate look-up arrays for readouts + psArray *readouts = psArrayAlloc(nImages); + for (int i = 0; i < nImages; i++) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + readouts->data[i] = psMemIncrRefCounter(readout); + } + + // source analysis is done in S/N order (brightest first) + // XXX add this in (need to put S/N in pmPhotObj) : objects = psArraySort (objects, pmPhotObjSortBySN); + + // process the objects in order. + for (int i = 0; i < objects->n; i++) { + pmPhotObj *object = objects->data[i]; + if (!object) continue; + if (!object->sources) continue; + + // choose the sources of interest + for (int j = 0; j < object->sources->n; j++) { + + pmSource *source = object->sources->data[j]; + + // 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; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) 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; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Next ++; + + int index = source->imageID; + pmReadout *readout = readouts->data[index]; + + // force source image to be a bit larger... + float radius = source->peak->xf - source->pixels->col0; + radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius); + + // if we request any of these measurements, we require the radial profile + if (doPetrosian || doIsophotal || doAnnuli || doKron) { + if (!psphotRadialProfile (source, recipe, skynoise, 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->Mx, source->moments->My); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + continue; + } + source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; + } + + // Petrosian Mags + if (doPetrosian) { + if (!psphotPetrosian (source, recipe, skynoise, maskVal)) { + psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + } else { + psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + Npetro ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; + } + } + + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + if (source->extpars) { + psFree(source->extpars->radFlux); + psFree(source->extpars->ellipticalFlux); + } + } + } + + psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); + psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); + psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); + + psFree(readouts); + return true; +} + +// sort by X (ascending) +# if (0) +int pmPhotObjSortBySN (const void **a, const void **b) +{ + pmPhotObj *objA = *(pmPhotObj **)a; + pmPhotObj *objB = *(pmPhotObj **)b; + + psF32 fA = objA->x; + psF32 fB = objB->x; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (+1); + if (diff < FLT_EPSILON) return (-1); + return (0); +} +# endif Index: /branches/eam_branches/psphot.20100506/src/psphotFitSourcesLinearStack.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotFitSourcesLinearStack.c (revision 27938) +++ /branches/eam_branches/psphot.20100506/src/psphotFitSourcesLinearStack.c (revision 27939) @@ -28,5 +28,4 @@ // analysis is done in spatial order (to speed up overlap search) - // sort by first element in each source list objects = psArraySort (objects, pmPhotObjSortByX); Index: /branches/eam_branches/psphot.20100506/src/psphotRadialAperturesByObject.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotRadialAperturesByObject.c (revision 27939) +++ /branches/eam_branches/psphot.20100506/src/psphotRadialAperturesByObject.c (revision 27939) @@ -0,0 +1,112 @@ +# include "psphotInternal.h" + +// aperture-like measurements for extended sources +// flux in simple, circular apertures +bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule) { + + bool status; + int Nradial = 0; + + psTimerStart ("psphot.radial"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping radial apertures\n"); + return true; + } + + // number of images used to define sources + int nImages = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // radMax stores the upper bounds of the annuli + // XXX keep the same name here as for the petrosian / elliptical apertures? + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe"); + psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // XXX temporary user-supplied systematic sky noise measurement (derive from background model) + float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + + // S/N limit to perform full non-linear fits + float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM"); + + // source analysis is done in S/N order (brightest first) + // XXX are we getting the objects out of order? does it matter? + // XXX objects = psArraySort (objects, pmPhotObjSortBySN); + + // generate look-up arrays for readouts + psArray *readouts = psArrayAlloc(nImages); + for (int i = 0; i < nImages; i++) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + readouts->data[i] = psMemIncrRefCounter(readout); + } + + // process the objects in order. + for (int i = 0; i < objects->n; i++) { + pmPhotObj *object = objects->data[i]; + if (!object) continue; + if (!object->sources) continue; + + // choose the sources of interest + for (int j = 0; j < object->sources->n; j++) { + + pmSource *source = object->sources->data[j]; + + // skip PSF-like and non-astronomical objects + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) 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; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Nradial ++; + + int index = source->imageID; + pmReadout *readout = readouts->data[index]; + + // force source image to be a bit larger... + float radius = source->peak->xf - source->pixels->col0; + radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); + radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius); + + if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) { + psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + } else { + source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; + } + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + } + + psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial); + + psFree(readouts); + return true; +} Index: /branches/eam_branches/psphot.20100506/src/psphotSourceMatch.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotSourceMatch.c (revision 27938) +++ /branches/eam_branches/psphot.20100506/src/psphotSourceMatch.c (revision 27939) @@ -1,5 +1,6 @@ # include "psphotInternal.h" -bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects); +bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects); +bool psphotMatchSourcesSetIDs (psArray *objects); psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view, const char *filerule) @@ -21,5 +22,9 @@ } - psphotMatchSourcesGenerate (config, view, filerule, objects); + // create sources for images where an object has been detected in the other images + psphotMatchSourcesAddMissing (config, view, filerule, objects); + + // choose a consistent position; set common sequence values + psphotMatchSourcesSetIDs (objects); return objects; @@ -145,5 +150,5 @@ } -bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) { +bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) { bool status = false; @@ -255,2 +260,16 @@ return true; } + +bool psphotMatchSourcesSetIDs (psArray *objects) { + + for (int i = 0; i < objects->n; i++) { + pmPhotObj *obj = objects->data[i]; + + // set the source->seq values + for (int j = 0; j < obj->sources->n; j++) { + pmSource *src = obj->sources->data[j]; + src->seq = i; + } + } + return true; +} Index: /branches/eam_branches/psphot.20100506/src/psphotStackReadout.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotStackReadout.c (revision 27938) +++ /branches/eam_branches/psphot.20100506/src/psphotStackReadout.c (revision 27939) @@ -150,5 +150,4 @@ // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) psphotFitSourcesLinearStack (config, objects, FALSE); - psFree (objects); // identify CRs and extended sources @@ -157,10 +156,11 @@ // measure aperture photometry corrections if (!psphotApResid (config, view, STACK_OUT)) { + psFree (objects); psLogMsg ("psphot", 3, "failed on psphotApResid"); return psphotReadoutCleanup (config, view, STACK_OUT); } - psphotRadialApertures (config, view, STACK_OUT); - psphotExtendedSourceAnalysis (config, view, STACK_OUT); // pass 1 (detections->allSources) + psphotRadialAperturesByObject (config, objects, view, STACK_OUT); + psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_OUT); // pass 1 (detections->allSources) psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources) @@ -186,4 +186,6 @@ psphotStackRemoveChisqFromInputs(config, STACK_RAW); + psFree (objects); + // create the exported-metadata and free local data return psphotReadoutCleanup (config, view, STACK_OUT);