Index: /branches/eam_branches/ipp-20130904/psphot/src/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/Makefile.am (revision 36342) +++ /branches/eam_branches/ipp-20130904/psphot/src/Makefile.am (revision 36343) @@ -205,4 +205,5 @@ psphotSourceFits.c \ psphotRadiusChecks.c \ + psphotChooseAnalysisOptions.c \ psphotOutput.c \ psphotFakeSources.c \ Index: /branches/eam_branches/ipp-20130904/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphot.h (revision 36342) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphot.h (revision 36343) @@ -125,4 +125,8 @@ bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add); bool psphotAddOrSubNoise_Threaded (psThreadJob *job); + +bool psphotChooseAnalysisOptions (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotChooseAnalysisOptionsReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotChooseAnalysisOptionsByObject(pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects); bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule); Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotChooseAnalysisOptions.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotChooseAnalysisOptions.c (revision 36342) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotChooseAnalysisOptions.c (revision 36343) @@ -3,4 +3,7 @@ // choose which sources will be processed for the petrosian radii and/or galaxy/trail model fits // a source for which we want galaxy models or petrosians gets + +// currently we also skip the following items (regardless of recipe values): +// TYPE_DEFECT, TYPE_SATURATED, MODE_DEFECT, MODE_SATSTAR, SATSTAR_PROFILE // RULES: @@ -11,5 +14,5 @@ // * select by Star/Galaxy? // ** PSPHOT.EXT.NSIGMA.LIMIT (sets EXT_LIMIT bit -- star or galaxy) -// ** PSPHOT.FIT.EXT_LIMIT (bool : ignore or not EXT_LIMIT) +// -- PSPHOT.FIT.EXT_LIMIT (bool : ignore or not EXT_LIMIT) // * select by S/N? @@ -17,20 +20,14 @@ // * select by Flux? +// * needs to depend on the filter +// * define a metadata block +// -- EXT.ANALYSIS.MAG.LIMITS (metadata block per filter) + // * select by galactic latitude (global override)? +// * if (|b| > LIM) do not limit the number of sources by density +// -- GLAT_MAX + // * de-select by density (global override)? - -// Fit all sources? -// * PSPHOT.EXT.FIT.ALL.SOURCES -// * PSPHOT.EXT.FIT.ALL.THRESH (density limit) - -// galactic lat limit? -// - -// density limit? -// magnitude limit? -// S/N limit? - -// DIFF MODE : -// * +// ?? // for which sources do we want to run the extended source analysis (petrosian and/or galaxy fits)? @@ -55,5 +52,5 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotSourceSizeReadout (config, view, filerule, i, recipe)) { + if (!psphotChooseAnalysisOptionsReadout (config, view, filerule, i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for %s entry %d", filerule, i); return false; @@ -62,4 +59,6 @@ return true; } + +bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos); // this function use an internal flag to mark sources which have already been measured @@ -77,4 +76,7 @@ psAssert (readout, "missing readout?"); + pmChip *chip = pmFPAviewThisChip(view, file->fpa); + psAssert (chip, "missing chip?"); + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); psAssert (detections, "missing detections?"); @@ -88,88 +90,152 @@ } + // we do the petrosian analysis for the same sources as the extended source fits IFF + // this recipe value is turned on (otherwise we only make a selection based on the mag limits) + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + // Option to enable fitting of all objects with extended model. - options.extFitAll = psMetadataLookupBool(&status, recipe, "PSPHOT.EXT.FIT.ALL.SOURCES"); + bool extFitAll = psMetadataLookupBool(&status, recipe, "PSPHOT.EXT.FIT.ALL.SOURCES"); assert (status); // Fitting everything is fine, but if the source density is high, we probably shouldn't. - options.extFitAllThresh = psMetadataLookupF32(&status, recipe, "PSPHOT.EXT.FIT.ALL.THRESH"); + float extFitAllThresh = psMetadataLookupF32(&status, recipe, "PSPHOT.EXT.FIT.ALL.THRESH"); assert (status); // Determine if this readout is above the threshold to ext fit all sources - options->extFitAllReadout = false; - if (options->extFitAll) { + if (extFitAll) { float maskFrac = psMetadataLookupF32(&status,readout->analysis,"READOUT.MASK.FRAC"); - if (status) { - maskFrac = 0.0; + if (status) maskFrac = 0.0; + if (sources->n * (1.0 - maskFrac) > extFitAllThresh) { + extFitAll = false; } - if (sources->n * (1.0 - maskFrac) > options->extFitAllThresh) { - options->extFitAllReadout = false; - } - else { - options->extFitAllReadout = true; - } - } - + } + + // use EXT_LIMIT bit to select objects? + bool useEXT_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.NSIGMA.LIMIT.USE"); + assert (status); + + float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); + assert (status); + + // use GAL_LIMIT to select / skip objects? + psSphereRot *toGal = NULL; + float GAL_LIMIT = 0.0; + bool useGAL_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.USE"); + assert (status); + if (useGAL_LIMIT) { + GAL_LIMIT = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT"); + toGal = psSphereRotICRSToGalactic(); + } + + // use MAG.LIMITS + psMetadata *magLimits = psMetadataLookupPtr (&status, recipe, "EXT.ANALYSIS.MAG.LIMITS"); + if (!status || !magLimits) { + psLogMsg ("psphot", PS_LOG_WARN, "EXT.ANALYSIS.MAG.LIMITS not found in the recipe, will use other criteria."); + magLimits = NULL; + } + + float extFitFluxLim = NAN; + + if (magLimits) { + float extFitMagLimDefault = NAN; + float extFitMagLim = NAN; + + // match to the given filter + psString filterID = psMetadataLookupStr(&status, file->fpa->concepts, "FPA.FILTERID"); + psAssert (filterID, "missing FPA.FILTERID?"); + + psMetadataIterator *iter = psMetadataIteratorAlloc(magLimits, PS_LIST_HEAD, NULL); + psMetadataItem *item = NULL; + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + if (item->type != PS_DATA_METADATA) { + psAbort ("Invalid type for EXT.ANALYSIS.MAG.LIMITS: %s, not a metadata folder", item->name); + } + psString thisFilter = psMetadataLookupStr (&status, item->data.md, "FILTER.ID"); + psAssert(thisFilter, "missing FILTER.ID"); + + // find a matching filter or default to 'any' + if (!strcasecmp (thisFilter, "any")) { + psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT"); + psAssert(extFitMagLimStr, "missing MAG.LIMIT"); + extFitMagLimDefault = atof (extFitMagLimStr); + } + + // find a matching filter or default to 'any' + if (!strcasecmp (thisFilter, filterID)) { + psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT"); + psAssert(extFitMagLimStr, "missing MAG.LIMIT"); + extFitMagLim = atof (extFitMagLimStr); + break; + } + } + if (!isfinite (extFitMagLim)) extFitMagLim = extFitMagLimDefault; + + // now I need to convert the mag limits into instrumental flux limits + // I need to get a zero point and exposure time for this image + + // select the exposure time + float exptime = psMetadataLookupF32(&status, file->fpa->concepts, "CELL.EXPOSURE"); + if (!status) { + exptime = psMetadataLookupF32(&status, file->fpa->concepts, "FPA.EXPOSURE"); + psAssert (status, "missing CELL.EXPOSURE and FPA.EXPOSURE?"); + } + + // select the exposure time + float zeropt = psMetadataLookupF32(&status, file->fpa->concepts, "FPA.ZP"); + psAssert (status, "missing FPA.ZP?"); + + extFitFluxLim = exptime * pow (10.0, 0.4*(zeropt - extFitMagLim)); + } + + pmSourceTmpF clearBits = ~(PM_SOURCE_TMPF_EXT_FIT | PM_SOURCE_TMPF_PETRO); for (psS32 i = 0 ; i < sources->n ; i++) { pmSource *source = (pmSource *) sources->data[i]; + + // clear the 2 relevant bits + source->tmpFlags &= clearBits; + + // 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; // skip saturated stars modeled with a radial profile if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; - // If the recipe requests we do extended source fits to everything, set - // the EXT_LIMIT flag - if (options->extFitAllReadout) { - psTrace("psphotSourceClassRegion.EXTALT",10,"In extFitAll: %d %d\n",options->extFitAll,options->extFitAllReadout); - source->mode |= PM_SOURCE_MODE_EXT_LIMIT; - } - - // saturated star (too many saturated pixels or peak above saturation limit). These - // may also be saturated galaxies, or just large saturated regions. - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; - continue; - } - - // any sources missing a large fraction should just be treated as PSFs - if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) { - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; - continue; - } - - bool isCR = (source->moments->SN > 7.0) && (Mminor < options->sizeLimitCR); - if (isCR) { - psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g CR\n", - source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid); - source->mode |= PM_SOURCE_MODE_DEFECT; - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE; - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; - Ncr ++; - continue; - } - - source->mode |= PM_SOURCE_MODE_EXT_LIMIT; - - // Alternate extended source limit calculation - if (options->altDiffExt) { - // ratio of major to minor axes - // MRV = Major / Minor - // MRV = (0.5 * (Mxx + Myy) + 0.5 * sqrt( (Mxx + Myy)^2 + 4 Mxy^2)) / - // (0.5 * (Mxx + Myy) - 0.5 * sqrt( (Mxx + Myy)^2 + 4 Mxy^2)) - // MRV = (2 * 0.5 * (Mxx + Myy) - Minor) / Minor - float momentRatioVeres = (Mxx + Myy - Mminor) / Mminor; - bool isAltEXT = (momentRatioVeres > options->altDiffExtThresh); - if (isAltEXT) { - psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g ALTEXT\t%g %g\n", - source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid, - momentRatioVeres,options->altDiffExtThresh); - source->mode |= PM_SOURCE_MODE_EXT_LIMIT; - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; - Next ++; - continue; - } - } - - // Everything else should just be treated as a PSF + // Do the fits if the recipe requests we do extended source fits to everything + if (extFitAll) { + source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT; + if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO; + continue; + } + + if (useEXT_LIMIT) { + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) { + // not extended so skip + continue; + } + } + + if (useGAL_LIMIT) { + psSphere ptGal, ptSky; + GetGalacticCoords (&ptGal, &ptSky, toGal, chip, source->peak->xf, source->peak->yf); + if (fabs(ptGal.d) < GAL_LIMIT) continue; + // include an exception for low density skycells below the limit? + } + + // for petro and extFit, we will either use the mag limits or the S/N + if (isfinite(extFitFluxLim)) { + if (source->moments->KronFlux > extFitFluxLim) { + source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT; + if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO; + } + } else { + if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) { + source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT; + if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO; + } + } } @@ -178,2 +244,307 @@ return true; } + +// this function use an internal flag to mark sources which have already been measured +bool psphotChooseAnalysisOptionsByObject(pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) +{ + bool status; + + psTimerStart ("psphot.options"); + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Choose Analysis Options (By Object) ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // S/N lim to perform radial aperture analysis + float SN_LIM_RADIAL = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM"); + + // we do the petrosian analysis for the same sources as the extended source fits IFF + // this recipe value is turned on (otherwise we only make a selection based on the mag limits) + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + + // Option to enable fitting of all objects with extended model. + bool extFitAll = psMetadataLookupBool(&status, recipe, "PSPHOT.EXT.FIT.ALL.SOURCES"); + assert (status); + + // Fitting everything is fine, but if the source density is high, we probably shouldn't. + float extFitAllThresh = psMetadataLookupF32(&status, recipe, "PSPHOT.EXT.FIT.ALL.THRESH"); + assert (status); + + // use EXT_LIMIT bit to select objects? + bool useEXT_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.NSIGMA.LIMIT.USE"); + assert (status); + + float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); + assert (status); + + // use GAL_LIMIT to select / skip objects? + psSphereRot *toGal = NULL; + float GAL_LIMIT = 0.0; + bool useGAL_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.USE"); + assert (status); + if (useGAL_LIMIT) { + GAL_LIMIT = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT"); + toGal = psSphereRotICRSToGalactic(); + } + + // use MAG.LIMITS + psMetadata *magLimits = psMetadataLookupPtr (&status, recipe, "EXT.ANALYSIS.MAG.LIMITS"); + if (!status || !magLimits) { + psLogMsg ("psphot", PS_LOG_WARN, "EXT.ANALYSIS.MAG.LIMITS not found in the recipe, will use other criteria."); + magLimits = NULL; + } + + int num = psphotFileruleCount(config, filerule); + + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // Find the filter for each of the inputs in fpa concepts + psArray *inputFilters = psArrayAlloc(num); + psArray *chips = psArrayAlloc(num); + psVector *zeropt = psVectorAlloc (inputFilters->n, PS_TYPE_F32); + psVector *exptime = psVectorAlloc (inputFilters->n, PS_TYPE_F32); + + // we will not use the chisq image to set the fitting limits + if (chisqNum >= 0) { + inputFilters->data[chisqNum] = psStringCopy("chisq"); + chips->data[chisqNum] = NULL; + zeropt->data.F32[chisqNum] = NAN; + exptime->data.F32[chisqNum] = NAN; + } + + // get the needed metadata for the non-chisq images + for (int i = 0 ; i < num; i++) { + if (i == chisqNum) continue; + + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); + psAssert (file, "missing file?"); + + pmChip *chip = pmFPAviewThisChip(view, file->fpa); + psAssert (chip, "missing chip?"); + chips->data[i] = psMemIncrRefCounter (chip); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + float maskFrac = psMetadataLookupF32(&status,readout->analysis,"READOUT.MASK.FRAC"); + if (!status) maskFrac = 0.0; + if (sources->n * (1.0 - maskFrac) > extFitAllThresh) { + extFitAll = false; + } + + // select the filterID for this image + char *filterID = psMetadataLookupStr(&status, file->fpa->concepts, "FPA.FILTERID"); + psAssert (status, "missing FPA.FILTERID?"); + inputFilters->data[i] = psStringCopy (filterID); + + // select the exposure time for this image + exptime->data.F32[i] = psMetadataLookupF32(&status, file->fpa->concepts, "CELL.EXPOSURE"); + if (!status) { + exptime->data.F32[i] = psMetadataLookupF32(&status, file->fpa->concepts, "FPA.EXPOSURE"); + psAssert (status, "missing CELL.EXPOSURE and FPA.EXPOSURE?"); + } + + // select the zero point for this image + zeropt->data.F32[i] = psMetadataLookupF32(&status, file->fpa->concepts, "FPA.ZP"); + psAssert (status, "missing FPA.ZP?"); + } + + // find extFitFluxLim->data.F32[i] for i == image number + psVector *extFitFluxLim = NULL; + if (magLimits) { + extFitFluxLim = psVectorAlloc (inputFilters->n, PS_TYPE_F32); + psVectorInit (extFitFluxLim, NAN); + + float extFitMagLimDefault = NAN; + + // match mag limits (flux limits) to the filters + psMetadataIterator *iter = psMetadataIteratorAlloc(magLimits, PS_LIST_HEAD, NULL); + psMetadataItem *item = NULL; + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + if (item->type != PS_DATA_METADATA) { + psAbort ("Invalid type for EXT.ANALYSIS.MAG.LIMITS: %s, not a metadata folder", item->name); + } + psString thisFilter = psMetadataLookupStr (&status, item->data.md, "FILTER.ID"); + psAssert(thisFilter, "missing FILTER.ID"); + + // save the default value to assign to unset filters + if (!strcasecmp (thisFilter, "any")) { + psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT"); + psAssert(magString, "missing MAG.LIMIT"); + extFitMagLimDefault = atof (magString); + continue; + } + + // not every entry in the metadata block needs to match to an image in our list + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; + if (!strcasecmp (thisFilter, inputFilters->data[i])) { + psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT"); + psAssert(magString, "missing MAG.LIMIT"); + float magvalue = atof (magString); + extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - magvalue)); + break; + } + } + } + + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; + if (isfinite(extFitFluxLim->data.F32[i])) continue; + extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - extFitMagLimDefault)); + } + } + + pmSourceTmpF clearBits = ~(PM_SOURCE_TMPF_EXT_FIT | PM_SOURCE_TMPF_PETRO); + + for (int i = 0; i < objects->n; i++) { + pmPhotObj *object = objects->data[i]; + if (!object) continue; + if (!object->sources) continue; + + // we check each source for an object and keep the object if any source is valid + + bool doObjectRadial = false; + bool doObjectExtFit = false; + for (int j = 0; !doObjectExtFit && !doObjectRadial && (j < object->sources->n); j++) { + + pmSource *source = object->sources->data[j]; + if (!source) continue; + if (!source->peak) continue; + if (source->imageID < 0) continue; // skip sources which come from other images? + if (source->imageID >= num) continue; // skip sources which come from other images? + + // 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; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + // XXX should I fit all even if one of the detections matches the above? + + // check on radial aperture analysis (fewer options for now) + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + bool doSourceRadial = (source->moments->KronFlux > SN_LIM_RADIAL * source->moments->KronFluxErr); + doObjectRadial = doObjectRadial | doSourceRadial; + } else { + bool doSourceRadial = (sqrt(source->peak->detValue) > SN_LIM_RADIAL); + doObjectRadial = doObjectRadial | doSourceRadial; + } + + // Do the fits if the recipe requests we do extended source fits to everything + if (extFitAll) { + doObjectExtFit = true; + continue; + } + + if (useEXT_LIMIT) { + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) { + continue; // not extended so skip + } + } + + int imageID = source->imageID; + + if (useGAL_LIMIT) { + psSphere ptGal, ptSky; + GetGalacticCoords (&ptGal, &ptSky, toGal, chips->data[imageID], source->peak->xf, source->peak->yf); + if (fabs(ptGal.d) < GAL_LIMIT) continue; + // include an exception for low density skycells below the limit? + } + + float fluxLim = extFitFluxLim ? extFitFluxLim->data.F32[imageID] : NAN; + // for petro and extFit, we will either use the mag limits or the S/N + if (isfinite(fluxLim)) { + if (source->moments->KronFlux > extFitFluxLim->data.F32[imageID]) { + doObjectExtFit = true; + } + } else { + if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) { + doObjectExtFit = true; + } + } + } + + for (int j = 0; j < object->sources->n; j++) { + + pmSource *source = object->sources->data[j]; + if (!source) continue; + if (!source->peak) continue; + + // clear the 2 relevant bits + source->tmpFlags &= clearBits; + + // 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; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + // Do the fits if the recipe requests we do extended source fits to everything + if (doObjectExtFit) { + source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT; + if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO; + } + + // Do the fits if the recipe requests we do extended source fits to everything + if (doObjectRadial) { + source->tmpFlags |= PM_SOURCE_TMPF_RADIAL_KEEP; + source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_SKIP; + } else { + source->tmpFlags |= PM_SOURCE_TMPF_RADIAL_SKIP; + source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_KEEP; + } + } + } + + psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld objects: %f sec\n", objects->n, psTimerMark ("psphot.options")); + + return true; +} + +bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos) { + + pmFPA *fpa = chip->parent; + + if (!chip->toFPA) goto escape; + if (!fpa->toTPA) goto escape; + if (!fpa->toSky) goto escape; + + psPlane ptCH, ptFP, ptTP; + + // calculate the astrometry for the coordinate of interest + ptCH.x = xPos; + ptCH.y = yPos; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); + psDeproject (ptSky, &ptTP, fpa->toSky); + psSphereRotApply (ptGal, toGal, ptSky); + + return true; + +escape: + // no astrometry calibration, give up + ptSky->r = NAN; + ptSky->d = NAN; + + ptGal->r = NAN; + ptGal->d = NAN; + + return false; +} Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotOldCode.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotOldCode.c (revision 36343) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotOldCode.c (revision 36343) @@ -0,0 +1,350 @@ + +/*************** + + In the process of developing psphot, I have written a bunch of code that has been used + for some tests or for a period, but which we no longer think is appropriate. I am + trying to clean up the kruft of old psphot bits and will but various such functions in + this file if they are worth keeping + +****************/ + +bool psphotMaskCosmicRayFootprintCheck (psArray *sources); + + +// maybe move this into psModules, pmFootprints? +bool psphotMaskCosmicRayFootprintCheck (psArray *sources) { +#ifdef CHECK_FOOTPRINTS + // This gets really expensive for complex images + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + pmPeak *peak = source->peak; + pmFootprint *footprint = peak->footprint; + if (!footprint) continue; + for (int j = 0; j < footprint->spans->n; j++) { + pmSpan *sp = footprint->spans->data[j]; + psAssert (sp, "missing span"); + } + } +#endif + return true; +} + +// mark the isophotal boundary +bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { + + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + pmPeak *peak = source->peak; + psAssert (peak, "NULL peak"); + + psImage *mask = source->maskView; + psImage *pixels = source->pixels; + psImage *variance = source->variance; + + // XXX This should be a recipe variable +# define SN_LIMIT 5.0 + + int xo = peak->x - pixels->col0; + int yo = peak->y - pixels->row0; + + // mark the pixels in this row to the left, then the right + for (int ix = xo; ix >= 0; ix--) { + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } + } + for (int ix = xo + 1; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } + } + + // for each of the neighboring rows, mark the high pixels if they have a marked neighbor + // first go up: + for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) { + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } + } + // next go down: + for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) { + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } + } + return true; +} + +// This attempt to mask the cosmic rays used the isophotal boundary +bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { + + // replace the source flux + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + // flag this as a CR + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + pmPeak *peak = source->peak; + psAssert (peak, "NULL peak"); + + // grab the matching footprint + pmFootprint *footprint = peak->footprint; + if (!footprint) { + psTrace("psphot.czw",2,"Using isophot CR mask code."); + + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + return true; + } + + if (!footprint->spans) { + psTrace("psphot.czw",2,"Using isophot CR mask code."); + + // if we have no footprint, use the old code to mask by isophot + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + return true; + } + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + // mask all of the pixels covered by the spans of the footprint + for (int j = 1; j < footprint->spans->n; j++) { + pmSpan *span1 = footprint->spans->data[j]; + + int iy = span1->y; + int xs = span1->x0; + int xe = span1->x1; + + for (int ix = xs; ix < xe; ix++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } + } + return true; +} + +// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total +// deviation in sigmas. This is measured on the residual image - should we ignore negative +// deviations? NOTE: This function was an early attempt to classify extended objects, and is +// no longer used by psphot. +float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, + psImageMaskType maskVal, const pmModel *model, float Ro) +{ + psF32 *PAR = model->params->data.F32; // Model parameters + float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters + + // We treat the contour as an ellipse: + // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY + // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0; + // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro + // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so: + // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2] + + // min/max value of x is where B^2 - 4AC = 0; solve this for x + float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0); + if (Q < 0.0) { + // ellipse is imaginary + return NAN; + } + + int radius = sqrtf(Q) + 0.5; // Radius of ellipse + int nPts = 0; // Number of points in ellipse + float nSigma = 0.0; // + + for (int x = -radius; x <= radius; x++) { + // Polynomial coefficients + // XXX Should we be using the centre of the pixel as x or x+0.5? + float A = PS_SQR (1.0 / syy); + float B = x * sxy; + float C = PS_SQR (x / sxx) - Ro; + float T = PS_SQR(B) - 4*A*C; + if (T < 0.0) { + continue; + } + + // y position in source frame + float yP = (-B + sqrt (T)) / (2.0 * A); + float yM = (-B - sqrt (T)) / (2.0 * A); + + // Get the closest pixel positions (image frame) + 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 || xPix >= image->numCols) { + continue; + } + + if (yPixM >= 0 && yPixM < image->numRows && + !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]); + nSigma += dSigma; + nPts++; + } + + if (yPixM == yPixP) { + continue; + } + + if (yPixP >= 0 && yPixP < image->numRows && + !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]); + nSigma += dSigma; + nPts++; + } + } + nSigma /= nPts; + return nSigma; +} + +// this was an old attempt to identify cosmic rays based on the peak curvature +bool psphotSourcePeakCurvature (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) { + + // classify the sources based on the CR test (place this in a function?) + // XXX use an internal flag to mark sources which have already been measured + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // skip source if it was already measured + if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { + psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); + continue; + } + + // source must have been subtracted + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { + source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; + psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); + continue; + } + + psF32 **resid = source->pixels->data.F32; + psF32 **variance = source->variance->data.F32; + psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + + // Integer position of peak + 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. These are mostly caught as DEFECT + if (xPeak < 1 || xPeak > source->pixels->numCols - 2 || + yPeak < 1 || yPeak > source->pixels->numRows - 2) { + psTrace("psphot", 7, "Not calculating crNsigma due to edge\n"); + continue; + } + + // Skip sources with masked pixels. These are mostly caught as DEFECT + bool keep = true; + for (int iy = -1; (iy <= +1) && keep; iy++) { + for (int ix = -1; (ix <= +1) && keep; ix++) { + if (mask[yPeak+iy][xPeak+ix] & options->maskVal) { + keep = false; + } + } + } + if (!keep) { + psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); + continue; + } + + // Compare the central pixel with those on either side, for the four possible lines through it. + + // Soften variances (add systematic error) + float softening = options->soft * PS_SQR(source->peak->rawFlux); // Softening for variances + + // Across the middle: y = 0 + float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1]; + float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1]; + float nX = cX / sqrtf(dcX + softening); + + // Up the centre: x = 0 + float cY = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak+0] - resid[yPeak+1][xPeak+0]; + float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0]; + float nY = cY / sqrtf(dcY + softening); + + // Diagonal: x = y + float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; + float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1]; + float nL = cL / sqrtf(dcL + softening); + + // Diagonal: x = - y + float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; + float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1]; + float nR = cR / sqrtf(dcR + softening); + + // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) + // Ndof = 4 ? (four measurements, no free parameters) + // XXX this value is going to be biased low because of systematic errors. + // we need to calibrate it somehow + // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq); + + // not strictly accurate: overcounts the chisq contribution from the center pixel (by + // factor of 4); also biases a bit low if any pixels are masked + // XXX I am not sure I want to keep this value... + source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR); + + float fCR = 0.0; + int nCR = 0; + if (nX > 0.0) { + fCR += nX; + nCR ++; + } + if (nY > 0.0) { + fCR += nY; + nCR ++; + } + if (nL > 0.0) { + fCR += nL; + nCR ++; + } + if (nR > 0.0) { + fCR += nR; + nCR ++; + } + source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + + if (!isfinite(source->crNsigma)) { + continue; + } + + // this source is thought to be a cosmic ray. flag the detection and mask the pixels + if (source->crNsigma > options->nSigmaCR) { + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask); + // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask); + } + } + + // now that we have masked pixels associated with CRs, we can grow the mask + if (options->grow > 0) { + bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask + psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow); + psImageConvolveSetThreads(oldThreads); + if (!newMask) { + psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask"); + return false; + } + psFree(readout->mask); + readout->mask = newMask; + } + return true; +} +