Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 36366) +++ trunk/psphot/src/Makefile.am (revision 36375) @@ -25,5 +25,5 @@ libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotForced psphotMinimal psphotMakePSF psphotStack psphotModelTest +bin_PROGRAMS = psphot psphotForced psphotFullForce psphotMinimal psphotMakePSF psphotStack psphotModelTest psmakecff # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy @@ -36,4 +36,8 @@ psphotForced_LDADD = libpsphot.la +psphotFullForce_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +psphotFullForce_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +psphotFullForce_LDADD = libpsphot.la + psphotMinimal_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) psphotMinimal_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) @@ -51,4 +55,9 @@ psphotModelTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) psphotModelTest_LDADD = libpsphot.la + +psmakecff_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +psmakecff_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +psmakecff_LDADD = libpsphot.la + # psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -77,4 +86,13 @@ psphotForced.c \ psphotForcedArguments.c \ + psphotParseCamera.c \ + psphotImageLoop.c \ + psphotMosaicChip.c \ + psphotCleanup.c + +# generalized forced photometry (including Kron, Petro, and Models) of specified positions given a specified psf +psphotFullForce_SOURCES = \ + psphotFullForce.c \ + psphotFullForceArguments.c \ psphotParseCamera.c \ psphotImageLoop.c \ @@ -129,4 +147,8 @@ psphotCleanup.c +# a program that takes a cmf file and makes a cff file (input file for psphotFullForce +psmakecff_SOURCES = \ + psmakecff.c + # psphotTest_SOURCES = \ # psphotTest.c @@ -155,4 +177,6 @@ psphotReadoutMinimal.c \ psphotForcedReadout.c \ + psphotFullForceReadout.c \ + psphotFullForce.SourceStats.c \ psphotMakePSFReadout.c \ psphotModelBackground.c \ @@ -181,4 +205,5 @@ psphotSourceFits.c \ psphotRadiusChecks.c \ + psphotChooseAnalysisOptions.c \ psphotOutput.c \ psphotFakeSources.c \ @@ -198,4 +223,7 @@ psphotRadialPlot.c \ psphotKronIterate.c \ + psphotKronFlux.c \ + psphotPetroFlux.c \ + psphotGalaxyShape.c \ psphotRadialProfileWings.c \ psphotDeblendSatstars.c \ Index: trunk/psphot/src/psmakecff.c =================================================================== --- trunk/psphot/src/psmakecff.c (revision 36375) +++ trunk/psphot/src/psmakecff.c (revision 36375) @@ -0,0 +1,171 @@ +# ifdef HAVE_CONFIG_H +# include +# endif + +#include +#include +#include +#include "psphot.h" + +static pmConfig* psmakecffArguments(int, char**); +static bool psmakecffParseCamera(pmConfig *); +static bool psmakecffImageLoop(pmConfig*); + +int main (int argc, char **argv) { + + psMemInit(); + psTimerStart ("complete"); + pmErrorRegister(); // register psModule's error codes/messages + psphotInit(); + + // load command-line arguments, options, and system config data + pmConfig *config = psmakecffArguments (argc, argv); + assert(config); + +// psphotVersionPrint(); + + // load input data (config and images (signal, noise, mask) + if (!psmakecffParseCamera (config)) { + psErrorStackPrint(stderr, "Error setting up the camera\n"); + exit (1); + } + + // call psphot for each readout + if (!psmakecffImageLoop (config)) { + psErrorStackPrint(stderr, "Error in the psphot image loop\n"); + exit(1); + } + +// psLogMsg ("psphot", PS_LOG_WARN, "complete psphot run: %f sec\n", psTimerMark ("complete")); + +// psErrorCode exit_status = psphotGetExitStatus(); +// psphotCleanup (config); + exit (0); +} + +// all functions which return to this level must raise one of the top-level error codes if they +// exit with an error. these error codes are used to specify the program exit status + +void usage() { + fprintf(stderr, "usage: psmakecff -sources \n"); + exit (1); +} + +static pmConfig* psmakecffArguments(int argc, char **argv) { + + pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE); + int N; + + if (config == NULL) { + psErrorStackPrint(stderr, "Can't read site configuration"); + exit(PS_EXIT_CONFIG_ERROR); + } + if ((N = psArgumentGet (argc, argv, "-sources"))) { + pmConfigFileSetsMD(config->arguments, &argc, argv, "SOURCES", "-sources", "-sourceslist"); + } else { + usage(); + } + if (argc < 2) { + usage(); + } + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", PS_DATA_STRING, "", argv[1]); + + return config; +} +static bool psmakecffParseCamera(pmConfig *config) { + bool status = false; + + pmFPAfile *sources = pmFPAfileDefineFromArgs(&status, config, "PSPHOT.INPUT.CMF", "SOURCES"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "failed to build FPA from sources file"); + return false; + } + + pmFPAfile *output = pmFPAfileDefineOutputFromFile (config, sources, "PSPHOT.OUTPUT.CFF"); + if (!output) { + psError(PS_ERR_UNKNOWN, false, "failed to build output FPA"); + return false; + } + output->save = true; + + return true; +} + +# define ESCAPE(MESSAGE) { \ + psError(PSPHOT_ERR_DATA, false, MESSAGE); \ + psFree (view); \ + return false; \ + } + +bool psmakecffImageLoop (pmConfig *config) { + + bool status; + pmChip *chip; + pmCell *cell; + pmReadout *readout; + + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT.CMF"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); + return false; + } + + pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSPHOT.OUTPUT.CFF"); + if (!output) { + psError(PSPHOT_ERR_PROG, false, "Can't find output data!"); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + pmHDU *lastHDU = NULL; // Last HDU updated + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + +// psImageMaskType maskTest = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + + // for psphot, we force data to be read at the chip level + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psLogMsg ("psmakecmf", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (! chip->process || ! chip->file_exists) { continue; } + +#ifdef notdef + pmFPAfileActivate (config->files, false, NULL); + pmFPAfileActivate (config->files, true, "PSPHOT.LOAD"); +#endif + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psmakecff."); + + // there is now only a single chip (multiple readouts?). loop over it and process + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psLogMsg ("psmakecmf", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + psLogMsg ("psmakecmf", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (! readout->data_exists) { continue; } + + pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); + if (hdu && hdu != lastHDU) { + // psphotVersionHeaderFull(hdu->header); + lastHDU = hdu; + } + } + + } + // Defer output and closing of files until we've (possibly) done the NFrames analysis below +// pmFPAfileActivate (config->files, false, NULL); + + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for Chip in psmakecff."); + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for FPA in psphot."); + + // fail if we failed to handle an error + if (psErrorCodeLast() != PS_ERR_NONE) psAbort ("failed to handle an error!"); + + psFree (view); + return true; +} Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 36366) +++ trunk/psphot/src/psphot.h (revision 36375) @@ -18,4 +18,5 @@ PSPHOT_SINGLE, PSPHOT_FORCED, + PSPHOT_FULL_FORCE, PSPHOT_MAKE_PSF, PSPHOT_MODEL_TEST, @@ -81,4 +82,7 @@ bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool setWindow); +bool psphotFullForceSourceStats (pmConfig *config, const pmFPAview *view, const char *filerule, bool setWindow); +bool psphotFullForceSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool setWindow); + bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); @@ -121,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); @@ -341,4 +349,7 @@ bool psphotForcedReadout(pmConfig *config, const pmFPAview *view, const char *filerule); +pmConfig *psphotFullForceArguments(int argc, char **argv); +bool psphotFullForceReadout(pmConfig *config, const pmFPAview *view, const char *filerule); + pmConfig *psphotMinimalArguments(int argc, char **argv); bool psphotReadoutMinimal(pmConfig *config, const pmFPAview *view, const char *filerule); @@ -414,4 +425,13 @@ } psphotStackOptions; +typedef struct { + float fRmajorMin; + float fRmajorMax; + float fRmajorDel; + float fRminorMin; + float fRminorMax; + float fRminorDel; +} psphotGalaxyShapeOptions; + /*** psphotStackMatchPSF prototypes ***/ bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view); @@ -497,4 +517,21 @@ bool psphotKronIterate_Threaded (psThreadJob *job); +bool psphotKronFlux (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotKronFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources); +bool psphotKronFlux_Threaded (psThreadJob *job); +bool psphotKronFluxSource (pmSource *source, psImageMaskType maskVal); + +bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources); +bool psphotPetroFlux_Threaded (psThreadJob *job); +bool psphotPetroFluxSource (pmSource *source, psImageMaskType maskVal); + +bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf); +bool psphotGalaxyShape_Threaded (psThreadJob *job); +bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize); +bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults); +psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc(); + bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotRadialProfileWingsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources); Index: trunk/psphot/src/psphotArguments.c =================================================================== --- trunk/psphot/src/psphotArguments.c (revision 36366) +++ trunk/psphot/src/psphotArguments.c (revision 36375) @@ -217,4 +217,5 @@ pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist"); pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC", "-src", "-srclist"); + // XXX allow this format? pmConfigFileSetsMD (config->arguments, &argc, argv, "FORCE", "-force", "-forcelist"); pmConfigFileSetsMD (config->arguments, &argc, argv, "EXPNUM", "-expnum", "-expnumlist"); Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 36366) +++ trunk/psphot/src/psphotBlendFit.c (revision 36375) @@ -274,8 +274,9 @@ # if (PS_TRACE_ON) - int TEST_ON = false; # define TEST_X 653 # define TEST_Y 466 - if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { +# define TESTING 1 + int TEST_ON = false; + if (TESTING && (fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { fprintf (stderr, "test object\n"); psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); Index: trunk/psphot/src/psphotChooseAnalysisOptions.c =================================================================== --- trunk/psphot/src/psphotChooseAnalysisOptions.c (revision 36375) +++ trunk/psphot/src/psphotChooseAnalysisOptions.c (revision 36375) @@ -0,0 +1,556 @@ +# include "psphotInternal.h" + +// 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: +// * global override: +// ** PSPHOT.EXT.FIT.ALL.SOURCES +// ** PSPHOT.EXT.FIT.ALL.THRESH (density limit) + +// * 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) + +// * select by S/N? +// ** EXTENDED_SOURCE_MODELS : SNLIM (per model value) + +// * 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)? +// ?? + +// for which sources do we want to run the extended source analysis (petrosian and/or galaxy fits)? +bool psphotChooseAnalysisOptions (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Choose Analysis Options ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + 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; + } + } + 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 +bool psphotChooseAnalysisOptionsReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) +{ + bool status; + + psTimerStart ("psphot.options"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + 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?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + return true; + } + + // 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); + + // Determine if this readout is above the threshold to ext fit all sources + if (extFitAll) { + float maskFrac = psMetadataLookupF32(&status,readout->analysis,"READOUT.MASK.FRAC"); + if (status) maskFrac = 0.0; + if (sources->n * (1.0 - maskFrac) > extFitAllThresh) { + extFitAll = false; + } + } + + // 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; + + // skip sources without a psf model + if (source->modelPSF == NULL) continue; + + // 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; + } + } + } + + psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.options")); + + 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; + + // skip sources without a psf model + if (source->modelPSF == NULL) 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: trunk/psphot/src/psphotDefineFiles.c =================================================================== --- trunk/psphot/src/psphotDefineFiles.c (revision 36366) +++ trunk/psphot/src/psphotDefineFiles.c (revision 36375) @@ -150,4 +150,11 @@ } + if (psMetadataLookupPtr(NULL, config->arguments, "FORCE")) { + if (!pmFPAfileDefineFromArgs (&status, config, "PSPHOT.INPUT.CFF", "FORCE")) { + psError(PSPHOT_ERR_CONFIG, false, "Failed to find/build PSPHOT.INPUT.CFF"); + return status; + } + } + if (psMetadataLookupPtr(NULL, config->arguments, "SRCTEXT")) { // XXX cannot use pmFPAfileDefineFromArgs: this is explicitly a FITS-based I/O function Index: trunk/psphot/src/psphotEllipticalContour.c =================================================================== --- trunk/psphot/src/psphotEllipticalContour.c (revision 36366) +++ trunk/psphot/src/psphotEllipticalContour.c (revision 36375) @@ -64,4 +64,5 @@ psFree (y); psFree (yErr); + source->mode2 |= PM_SOURCE_MODE2_ECONTOUR_FEW_PTS; return false; } Index: trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 36366) +++ trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 36375) @@ -1,5 +1,6 @@ # include "psphotInternal.h" - -// measure the elliptical radial profile and use this to measure the petrosian parameters for the sources +void psphotRadialProfileShowSkips (); + +// measure the petrosian parameters for the sources // for now, let's store the detections on the readout->analysis for each readout @@ -15,6 +16,9 @@ psAssert (recipe, "missing recipe?"); - // perform full non-linear fits / extended source analysis? - if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) { + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + + // measure petrosians? + if (!doPetrosian && !doAnnuli) { psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); return true; @@ -33,4 +37,15 @@ } +/*** for the moment, this test code : it is not thread safe ***/ +int Nall = 0; +int Nskip1 = 0; +int Nskip2 = 0; +int Nskip3 = 0; +int Nskip4 = 0; +int Nskip5 = 0; +int Nskip6 = 0; + +# define SKIP(VALUE) { VALUE++; continue; } + // aperture-like measurements for extended sources bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { @@ -110,5 +125,5 @@ // set this to 0 to run without threading -# if (1) +# if (0) if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); @@ -164,10 +179,18 @@ psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); +# if (PS_TRACE_ON) + fprintf (stderr, "ext analysis skipped @ 1 : %d\n", Nskip1); + fprintf (stderr, "ext analysis skipped @ 2 : %d\n", Nskip2); + fprintf (stderr, "ext analysis skipped @ 3 : %d\n", Nskip3); + fprintf (stderr, "ext analysis skipped @ 4 : %d\n", Nskip4); + fprintf (stderr, "ext analysis skipped @ 5 : %d\n", Nskip5); + fprintf (stderr, "ext analysis skipped @ 6 : %d\n", Nskip6); +#endif + + psphotRadialProfileShowSkips (); + psphotVisualShowResidualImage (readout, false); - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - if (doPetrosian) { - psphotVisualShowPetrosians (sources); - } + psphotVisualShowPetrosians (sources); return true; @@ -190,11 +213,5 @@ float skynoise = PS_SCALAR_VALUE(job->args->data[4],F32); - // 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 doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doPetroStars = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS"); // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -207,42 +224,15 @@ pmSource *source = sources->data[i]; - // if we have checked the source validity on the basis of the object set, then - // we either skip these tests below or we skip the source completely - if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) continue; - if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource; - - // 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; - - // optionally allow non-extended objects to get petrosians as well - if (!doPetroStars) { - if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; - if (source->type == PM_SOURCE_TYPE_STAR) continue; - } - - // limit selection to some SN limit - // assert (source->peak); // how can a source not have a peak? - // limit selection to some SN limit - bool skipSource = false; - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr); - } else { - skipSource = (sqrt(source->peak->detValue) < SN_LIM); - } - if (skipSource) continue; - - // limit selection by analysis region (this automatically apply - if (source->peak->x < region->x0) continue; - if (source->peak->y < region->y0) continue; - if (source->peak->x > region->x1) continue; - if (source->peak->y > region->y1) continue; - - keepSource: + Nall ++; + + // rules for measuring petrosian parameters for specific objects are set in + // psphotChooseAnalysisOptions.c + if (!(source->tmpFlags & PM_SOURCE_TMPF_PETRO)) SKIP (Nskip1); + + // limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?) + if (source->peak->x < region->x0) SKIP (Nskip2); + if (source->peak->y < region->y0) SKIP (Nskip3); + if (source->peak->x > region->x1) SKIP (Nskip4); + if (source->peak->y > region->y1) SKIP (Nskip5); // replace object in image @@ -259,26 +249,24 @@ 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 || doAnnuli) { - 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; - } - Nannuli ++; - source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; - } + // measure the radial profile + if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) { + // 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); + SKIP (Nskip6); + } + + Nannuli ++; + 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; - } + 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; + } } Index: trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceFits.c (revision 36366) +++ trunk/psphot/src/psphotExtendedSourceFits.c (revision 36375) @@ -65,6 +65,4 @@ psphotVisualShowImage(readout); - // psphotSaveImage (NULL, readout->image, "test.01.fits"); - pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); psAssert (detections, "missing detections?"); @@ -111,4 +109,9 @@ if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) { fitMaxTol = 1.0; + } + + float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolution + if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) { + fitNsigmaConv = 5.0; } @@ -137,8 +140,13 @@ fitOptions->minTol = fitMinTol; fitOptions->maxTol = fitMaxTol; + fitOptions->nsigma = fitNsigmaConv; fitOptions->gainFactorMode = gainFactorMode; fitOptions->chisqConvergence = chisqConvergence; fitOptions->isInteractive = isInteractive; + + // use poissonian errors or local-sky errors + fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON"); + if (!status) fitOptions->poissonErrors = true; // maskVal is used to test for rejected pixels, and must include markVal @@ -285,4 +293,5 @@ psArrayAdd(job->args, 1, cells->data[j]); // sources psArrayAdd(job->args, 1, models); + // Allocate a metadata iterator here because psMetadataIteratorAlloc/Free are not thread safe psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); @@ -386,6 +395,4 @@ psphotSersicModelClassCleanup(); - // psphotSaveImage (NULL, readout->image, "test.02.fits"); - psphotVisualShowResidualImage (readout, false); @@ -428,4 +435,12 @@ // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); + pmModelStatus badModel = PM_MODEL_STATUS_NONE; + badModel |= PM_MODEL_STATUS_BADARGS; + badModel |= PM_MODEL_STATUS_OFFIMAGE; + badModel |= PM_MODEL_STATUS_NAN_CHISQ; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID; + badModel |= PM_MODEL_PCM_FAIL_GUESS; + // choose the sources of interest for (int i = 0; i < sources->n; i++) { @@ -433,23 +448,13 @@ pmSource *source = sources->data[i]; - // skip PSF-like and non-astronomical objects - if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; - if (source->type == PM_SOURCE_TYPE_DEFECT) continue; - if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - - // skip saturated stars modeled with a radial profile - if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; - - // XXX this should use peak? + // rules for measuring petrosian parameters for specific objects are set in + // psphotChooseAnalysisOptions.c + if (!(source->tmpFlags & PM_SOURCE_TMPF_EXT_FIT)) continue; + + // limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?) if (source->peak->x < region->x0) continue; if (source->peak->y < region->y0) continue; if (source->peak->x > region->x1) continue; if (source->peak->y > region->y1) continue; - - - // XXX for a test, just do the obvious trail - // XXX if (source->peak->xf < 1100) continue; - // XXX if (source->peak->xf > 1400) continue; - // XXX if (source->peak->yf > 245) continue; // replace object in image @@ -462,7 +467,4 @@ psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs - // XXX WATCH OUT HERE!! - // fitRadius = 30; - // UPDATE : we have changed the moments calculation. There is now an iteration within // psphotKronMasked to determine moments appropriate for a larger object. The values @@ -475,10 +477,7 @@ } +# ifdef TEST_OBJECT bool testObject = false; - // testObject |= ((fabs(source->peak->xf - 179) < 5) && (fabs(source->peak->yf - 1138) < 5)); - // testObject |= ((fabs(source->peak->xf - 5047) < 5) && (fabs(source->peak->yf - 151) < 5)); - // testObject |= ((fabs(source->peak->xf - 3929) < 5) && (fabs(source->peak->yf - 4109) < 5)); - // testObject |= ((fabs(source->peak->xf - 915) < 5) && (fabs(source->peak->yf - 5998) < 5)); - // testObject |= ((fabs(source->peak->xf - 5406) < 5) && (fabs(source->peak->yf - 326) < 5)); + testObject |= ((fabs(source->peak->xf - 179) < 5) && (fabs(source->peak->yf - 1138) < 5)); if (testObject) { fprintf (stderr, "test object @ %f, %f\n", source->peak->xf, source->peak->yf); @@ -486,4 +485,5 @@ psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 5); } +# endif // loop here over the models chosen for each source (exclude by S/N) @@ -501,15 +501,13 @@ assert (status); - // limit selection to some SN limit - bool skipSource = false; - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - skipSource = (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr); - } else { - skipSource = (sqrt(source->peak->detValue) < FIT_SN_LIM); - } - if (skipSource) { + // limit selection to some SN limit for specific models (this value only applies if > EXTENDED_SOURCE_SN_LIM) + if (isfinite(FIT_SN_LIM)) { + if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) { Nfaint ++; continue; + } } + + source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_RUN; // check on the model type @@ -529,4 +527,5 @@ psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My); Nfail ++; + source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_FAIL; continue; } @@ -534,5 +533,5 @@ source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter); Nconvolve ++; - if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { + if (!(modelFit->flags & badModel)) { NconvolvePass ++; source->mode |= PM_SOURCE_MODE_EXTENDED_FIT; @@ -548,9 +547,10 @@ Nfail ++; doneFits = true; + source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_FAIL; continue; } psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter); Nplain ++; - if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { + if (!(modelFit->flags & badModel)) { NplainPass ++; source->mode |= PM_SOURCE_MODE_EXTENDED_FIT; @@ -576,10 +576,10 @@ fprintf (stderr, "update window : %f %f : %f -> %f\n", source->peak->xf, source->peak->yf, fitRadius, 2*fitRadius); psphotSetWindowTrail (&fitRadius, &windowRadius, readout, source, markVal, fitRadius*2.0); + source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_RETRY; } } } } - // XXX really need to do this in a cleaner way: - if (!modelFit) continue; + psAssert (modelFit, "modelFit not set?"); // test for fit quality / result @@ -590,5 +590,16 @@ } + // we are allowed to fit both stars and non-stars here -- if we have fitted + // something which we think is a star, we should use that model to subtract the + // object from the image. + if (source->type == PM_SOURCE_TYPE_STAR) { + // ensure the modelPSF is cached + pmSourceCacheModel (source, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + continue; + } + // evaluate the relative quality of the models, choose one + // the PSF model might be the best fit : allow it to succeed float minChisq = NAN; int minModel = -1; @@ -596,9 +607,9 @@ pmModel *model = source->modelFits->data[i]; + // skip the really bad fits if (!(model->flags & PM_MODEL_STATUS_FITTED)) continue; - - if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue; + if (model->flags & badModel) continue; + // if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue; - if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue; if ((minModel < 0) || (model->chisq < minChisq)) { @@ -618,4 +629,5 @@ pmSourceCacheModel (source, maskVal); +# if (PS_TRACE_ON) pmModel *model = source->modelFits->data[0]; int flags = 0xffffffff; @@ -623,8 +635,9 @@ flags = model->flags; } - fprintf (stderr, "failed to fit extended source model to object %d @ %f, %f (%x)\n", source->id, source->moments->Mx, source->moments->My, flags); +#endif pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_NONE; continue; } @@ -643,9 +656,10 @@ if (source->modelEXT->isPCM) { // fprintf (stderr, "subtract PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My); - pmPCMCacheModel (source, maskVal, psfSize); + pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma); } else { // fprintf (stderr, "subtract non-PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My); pmSourceCacheModel (source, maskVal); } + source->modelEXT->flags |= PM_MODEL_BEST_FIT; // subtract the best fit from the object, leave local sky @@ -655,8 +669,10 @@ psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); +# ifdef TEST_OBJECT if (testObject) { psTraceSetLevel ("psModules.objects.pmPCM_MinimizeChisq", 0); psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 0); } +# endif } Index: trunk/psphot/src/psphotFullForce.SourceStats.c =================================================================== --- trunk/psphot/src/psphotFullForce.SourceStats.c (revision 36375) +++ trunk/psphot/src/psphotFullForce.SourceStats.c (revision 36375) @@ -0,0 +1,198 @@ +# include "psphotInternal.h" + +// construct pixels and measure moments for sources. +// NOTE: this function differs from psphotSourceStats in the following ways: +// 1) detections->allSources are used instead of newSources +// 2) the sources are used directly (peaks are not assigned to sources here) + +// NOTE: this function is meant to be used by psphotFullForce, and the mode of the loaded +// sources should have (MODE_EXTERNAL) set. If so, the Mx,My values will not be recalculated + +bool psphotFullForceSourceStats (pmConfig *config, const pmFPAview *view, const char *filerule, bool setWindow) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Source Stats ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + // int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + // if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + // if (i == chisqNum) continue; // skip chisq image + if (!psphotFullForceSourceStatsReadout (config, view, filerule, i, recipe, setWindow)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +bool psphotFullForceSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool setWindow) { + + bool status = false; + + psTimerStart ("psphot.stats"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + if (psTraceGetLevel("psphot") > 5) { + static int pass = 0; + char name[64]; + sprintf (name, "srstats.v%d.fits", pass); + psphotSaveImage(NULL, readout->image, name); + pass ++; + } + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } + + // determine properties (sky, moments) of initial sources + float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + psAssert (status, "missing SKY_OUTER_RADIUS in recipe?"); + + // XXX this seems like an arbitrary number... + OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius + + float INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); psAssert (status, "missing SKY_INNER_RADIUS"); + float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN"); + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + psAssert (maskVal, "missing MASK.PSPHOT"); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + psAssert (markVal, "missing MARK.PSPHOT"); + + psArray *sources = detections->allSources; + + // generate the array of sources, define the associated pixel + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + psAssert (source->peak, "how can we have a peak-less source?"); + psAssert (source->moments, "how can we have a moments-less source?"); + + // allocate image, weight, mask arrays for each peak (square of radius OUTER) + pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); + } + + if (setWindow) { + // calculate the best gaussian window (keep centroid positions) + // NOTE: if sources->mode,mode2 has MODE_EXTERNAL,MODE2_MATCHED, Mx,My are not changed + if (!psphotSetMomentsWindow(recipe, readout->analysis, sources, maskVal)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!"); + psFree(detections->newSources); + return false; + } + } + + // if we have measured the window, we will be saving the modified version of these recipe values on readout->analysis + float SIGMA = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_GAUSS_SIGMA"); + if (!status) { + SIGMA = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA"); + } + float RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS"); + if (!status) { + RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); + } + float MIN_KRON_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_MIN_KRON"); + if (!status) { + MIN_KRON_RADIUS = 0.75*SIGMA; + } + + // threaded measurement of the source magnitudes + int Nfail = 0; + int Nmoments = 0; + int Nfaint = 0; + + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS"); + + psArrayAdd(job->args, 1, cells->data[j]); // sources + + PS_ARRAY_ADD_SCALAR(job->args, INNER, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, MIN_SN, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, RADIUS, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, SIGMA, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32); + + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS"); + psFree(detections->newSources); + return false; + } + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS"); + psFree(detections->newSources); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[8]; + Nmoments += scalar->data.S32; + scalar = job->args->data[9]; + Nfail += scalar->data.S32; + scalar = job->args->data[10]; + Nfaint += scalar->data.S32; + } + psFree(job); + } + } + + psFree (cellGroups); + + psLogMsg ("psphot", PS_LOG_WARN, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats")); + + psphotVisualShowMoments (sources); + + return true; +} Index: trunk/psphot/src/psphotFullForce.c =================================================================== --- trunk/psphot/src/psphotFullForce.c (revision 36375) +++ trunk/psphot/src/psphotFullForce.c (revision 36375) @@ -0,0 +1,35 @@ +# include "psphotStandAlone.h" + +int main (int argc, char **argv) { + + psTimerStart ("complete"); + pmErrorRegister(); // register psModule's error codes/messages + psphotInit(); + + // load command-line arguments, options, and system config data + pmConfig *config = psphotFullForceArguments (argc, argv); + assert(config); + + psphotVersionPrint(); + + // load input data (config and images (signal, noise, mask) + if (!psphotParseCamera (config)) { + psErrorStackPrint(stderr, "Error setting up the camera\n"); + exit (psphotGetExitStatus()); + } + + // call psphot for each readout + if (!psphotImageLoop (config, PSPHOT_FULL_FORCE)) { + psErrorStackPrint(stderr, "Error in the psphot image loop\n"); + exit (psphotGetExitStatus()); + } + + psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete")); + + psErrorCode exit_status = psphotGetExitStatus(); + psphotCleanup (config); + exit (exit_status); +} + +// all functions which return to this level must raise one of the top-level error codes if they +// exit with an error. these error codes are used to specify the program exit status Index: trunk/psphot/src/psphotFullForceArguments.c =================================================================== --- trunk/psphot/src/psphotFullForceArguments.c (revision 36375) +++ trunk/psphot/src/psphotFullForceArguments.c (revision 36375) @@ -0,0 +1,202 @@ +# include "psphotStandAlone.h" + +static void usage(const char *program, psMetadata *arg, pmConfig *config, int exitCode); +static void writeHelpInfo(const char* program, pmConfig* config, FILE* ofile); + +pmConfig *psphotFullForceArguments(int argc, char **argv) { + + int N; + bool status, status1, status2, status3; + + // load config data from default locations + pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE); + if (config == NULL) { + psErrorStackPrint(stderr, "Can't read site configuration"); + exit(PS_EXIT_CONFIG_ERROR); + } + + PS_ARGUMENTS_GENERIC (psphot, config, argc, argv); + + // save the following additional recipe values based on command-line options + // these options override the PSPHOT recipe values loaded from recipe files + psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE); + + // Number of threads is handled + PS_ARGUMENTS_THREADS( psphot, config, argc, argv ); + + if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(argv[0], config, stdout); + if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(argv[0], config, stdout); + + // visual : interactive display mode + if ((N = psArgumentGet (argc, argv, "-visual"))) { + psArgumentRemove (N, &argc, argv); + pmVisualSetVisual(true); + } + + // break : used from recipe throughout psphotReadout + if ((N = psArgumentGet (argc, argv, "-break"))) { + if (argc<=N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + exit(PS_EXIT_CONFIG_ERROR); + } + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF + if ((N = psArgumentGet (argc, argv, "-region"))) { + if (argc<=N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + exit(PS_EXIT_CONFIG_ERROR); + } + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // chip selection is used to limit chips to be processed + if ((N = psArgumentGet (argc, argv, "-chip"))) { + if (argc<=N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + exit(PS_EXIT_CONFIG_ERROR); + } + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // if these command-line options are supplied, load the file name lists into config->arguments + // override any configuration-specified source for these files + // + pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist"); + pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist"); + + pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist"); + + status1 = pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC", "-src", "-srclist"); + status2 = pmConfigFileSetsMD (config->arguments, &argc, argv, "SRCTEXT", "-srctext", "-srctextlist"); + status3 = pmConfigFileSetsMD (config->arguments, &argc, argv, "FORCE", "-force", "-forcelist"); + + if (!(status1 || status2 || status3)) { + // XXX require -force version? + psError(PSPHOT_ERR_ARGUMENTS, true, "No source list is supplied (use one of -force, -forcelist, -src, -srctext, -srclist, or -srctextlist)"); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + if (argc == 1) { + psError(PSPHOT_ERR_ARGUMENTS, true, "Too few arguments: %d", argc); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + // the input file is a required argument; if not found, we will exit + status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list"); + if (!status) { + psError(PSPHOT_ERR_ARGUMENTS, false, "pmConfigFileSetsMD failed to parse arguments"); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + if (argc != 2) { + psError(PSPHOT_ERR_ARGUMENTS, true, "Expected to see one more argument; saw %d", argc - 1); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + // output position is fixed + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]); + + psTrace("psphot", 1, "Done with psphotFullForceArguments...\n"); + return (config); +} + +static void writeHelpInfo(const char* program, pmConfig* config, FILE* ofile) +{ + fprintf(ofile, + "Usage: one of the following\n" + "%s -file fname1[,fname2,...] -mask maskfile1[,maskfile2,...]\n" + " -variance varfile1[,varfile2,...] OutFileBaseName\n" + "\n" + "%s -list FileNameList [-masklist MaskFileNameList] \n" + " -variancelist VarFileNameList OutFileBaseName\n" + "\n" + "%s -help\n" + "\n" + "%s -version\n" + "\n" + "where:\n" + " FileNameList is a text file containing filenames, one per line\n" + " MaskFileNameList is a text file of mask filenames, one per line\n" + " VarFileNameList is a text file of variance filenames, one per line\n" + " OutFileBaseName is the 'root name' for output files\n" + "also required:\n" + " -src SrcFile1[,SrcFile2,...] or -srclist SrcFileNameList\n" + " specify sources to be measured (required)\n" + "\n" + "additional options:\n" + " -psf PsfFile1[,PsfFile2,...] or -psflist PsfFileNameList\n" + " specify PSF rather than letting %s estimate it\n" + " -chip nn[,nn,...]\n" + " select detector chips to process; default is all.\n" + " Indices correspond to zero-based offset in the FPA metadata table.\n" + " -region RegionString\n" + " specify analysis region. String is of form '[x0:x1,y0:y1]'\n" + " To use this option you must define a default in psphot.config\n" + " -visual\n" + " turns on interactive display mode\n" + " -dumpconfig CfgFileName\n" + " causes config info to be dumped to the named file.\n" + " -break NOTHING|BACKMDL|PEAKS|MOMENTS|PSFMODEL|ENSEMBLE|PASS1\n" + " choose a point at which to exit processing early\n" + " -threads n\n" + " set number of parallel threads of execution\n" + " -F OldFileRule ReplacementFileRule\n" + " change file naming rule; e.g. '-F PSPHOT.OUTPUT PSPHOT.OUT.CMF.MEF'\n" + " -D name stringval\n" + " set a string-valued config parameter\n" + " -Di name intval\n" + " set an integer-valued config parameter\n" + " -Df name fval\n" + " set a float-valued config parameter\n" + " -Db name boolval\n" + " set a boolean-valued config parameter\n" + " -v, -vv, -vvv\n" + " set increasing levels of verbosity\n" + " -logfmt FormatString\n" + " set format string used for log messages\n" + " -trace Fac Lvl\n" + " set tracing for facility Fac to integer Lvl, e.g. '-trace err 10'\n" + " -trace-levels\n" + " print current trace levels\n", + program,program,program,program,program); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(PS_EXIT_SUCCESS); +} + +static void usage(const char *program, // Name of the program + psMetadata *arguments, // Command-line arguments + pmConfig *config, // Configuration + int exitCode + ) +{ + fprintf(stderr, + "Usage: one of the following\n" + "%s -file fname1[,fname2,...] -mask maskfile1[,maskfile2,...]\n" + " -variance varfile1[,varfile2,...] OutFileBaseName\n" + "\n" + "%s -list FileNameList [-masklist MaskFileNameList] \n" + " -variancelist VarFileNameList OutFileBaseName\n" + "also required:\n" + " -src SrcFile1[,SrcFile2,...] or -srclist SrcFileNameList\n" + " specify sources to be measured (required)\n" + "\n" + "Try '%s -help' for more options and explanation\n", + program,program,program); + if (exitCode != PS_EXIT_SUCCESS) + psErrorStackPrint(stderr, "Error reading arguments\n"); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(exitCode); +} + Index: trunk/psphot/src/psphotFullForceImageLoop.c =================================================================== --- trunk/psphot/src/psphotFullForceImageLoop.c (revision 36375) +++ trunk/psphot/src/psphotFullForceImageLoop.c (revision 36375) @@ -0,0 +1,153 @@ +# include "psphotStandAlone.h" + +# define ESCAPE(MESSAGE) { \ + psError(PSPHOT_ERR_DATA, false, MESSAGE); \ + psFree (view); \ + return false; \ + } + +bool psphotFullForceImageLoop (pmConfig *config) { + + bool status; + pmChip *chip; + pmCell *cell; + pmReadout *readout; + + pmFPAfile *load = psMetadataLookupPtr (&status, config->files, "PSPHOT.LOAD"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); + return false; + } + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + pmHDU *lastHDU = NULL; // Last HDU updated + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + psImageMaskType maskTest = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + + // for psphot, we force data to be read at the chip level + while ((chip = pmFPAviewNextChip (view, load->fpa, 1)) != NULL) { + psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (! chip->process || ! chip->file_exists) { continue; } + + // load just the input image data (image, mask, weight) + pmFPAfileActivate (config->files, false, NULL); + pmFPAfileActivate (config->files, true, "PSPHOT.LOAD"); + pmFPAfileActivate (config->files, true, "PSPHOT.MASK"); + pmFPAfileActivate (config->files, true, "PSPHOT.VARIANCE"); + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphot."); + + // mosaic the cells of a chip into a single contiguous (trimmed) chip + if (!psphotMosaicChip(config, view, "PSPHOT.INPUT", "PSPHOT.LOAD")) ESCAPE ("Unable to mosaic chip."); + + // Read WCS if easy. + // XXX Since we're mosaicking cells, we ignore the case where the WCS is defined for a cell. + { + pmChip *inChip = pmFPAviewThisChip(view, input->fpa); // Mosaicked chip + pmHDU *hduLow = pmHDUGetLowest(input->fpa, inChip, NULL); + if (hduLow && !pmAstromReadWCS(input->fpa, inChip, hduLow->header, 1.0)) { + psWarning("Unable to read WCS astrometry from header."); + psErrorClear(); + pmHDU *hduHigh = pmHDUGetHighest(input->fpa, inChip, NULL); + if (hduHigh && hduHigh != hduLow && + !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) { + psWarning("Unable to read WCS astrometry from primary header."); + psErrorClear(); + } + } + } + + // try to load other supporting data (PSF, SRC, etc). + // do not re-load the following three files + pmFPAfileActivate (config->files, true, NULL); + pmFPAfileActivate (config->files, false, "PSPHOT.LOAD"); + pmFPAfileActivate (config->files, false, "PSPHOT.MASK"); + pmFPAfileActivate (config->files, false, "PSPHOT.VARIANCE"); + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphot."); + + // re-activate files so they will be closed and freed below + pmFPAfileActivate (config->files, true, NULL); + + // there is now only a single chip (multiple readouts?). loop over it and process + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (! readout->data_exists) { continue; } + + // Update the header + pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); + if (hdu && hdu != lastHDU) { + psphotVersionHeaderFull(hdu->header); + lastHDU = hdu; + } + + // if an external mask is supplied, ensure that NAN pixels are also masked + if (readout->mask) { + psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels + if (!pmReadoutMaskInvalid(readout, maskTest, maskSat)) { + psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); + psFree(view); + return false; + } + } + + // run the actual photometry analysis on this chip/cell/readout + if (!psphotFullForceReadout (config, view, "PSPHOT.INPUT")) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + } + + // drop all versions of the internal files + status = true; + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); + psFree (view); + return false; + } + } + + // save output which is saved at the chip level + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); + } + // save output which is saved at the fpa level + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot."); + + // fail if we failed to handle an error + if (psErrorCodeLast() != PS_ERR_NONE) psAbort ("failed to handle an error!"); + + psFree (view); + return true; +} + +// I/O files related to psphot: +// PSPHOT.INPUT : input image file(s) +// PSPHOT.RESID : residual image +// PSPHOT.OUTPUT : output object tables (object) + +// PSPHOT.BACKSUB : background subtracted image +// PSPHOT.BACKGND : background model (full-scale image?) +// PSPHOT.BACKMDL : background model (binned image?) +// PSPHOT.PSF : sample PSF images + +// PSPHOT.MASK +// PSPHOT.VARIANCE +// Index: trunk/psphot/src/psphotFullForceReadout.c =================================================================== --- trunk/psphot/src/psphotFullForceReadout.c (revision 36375) +++ trunk/psphot/src/psphotFullForceReadout.c (revision 36375) @@ -0,0 +1,110 @@ +# include "psphotInternal.h" + +bool psphotFullForceReadout(pmConfig *config, const pmFPAview *view, const char *filerule) { + + // measure the total elapsed time in psphotReadout. + psTimerStart ("psphotReadout"); + + // allow objects to be fit with ugly models (central holes, extreme asymmetry, etc) + pmModelClassSetLimits(PM_MODEL_LIMITS_LAX); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + + // set the photcode for this image + if (!psphotAddPhotcode (config, view, filerule)) { + psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); + return false; + } + + // optional break-point for processing + char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); + PS_ASSERT_PTR_NON_NULL (breakPt, false); + + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndVariance (config, view, filerule); + if (!strcasecmp (breakPt, "NOTHING")) { + return psphotReadoutCleanup (config, view, filerule); + } + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, filerule)) { + return psphotReadoutCleanup (config, view, filerule); + } + if (!psphotSubtractBackground (config, view, filerule)) { + return psphotReadoutCleanup (config, view, filerule); + } + if (!strcasecmp (breakPt, "BACKMDL")) { + return psphotReadoutCleanup (config, view, filerule); + } + + if (!psphotLoadPSF (config, view, filerule)) { + // this only happens if we had a programming error in psphotLoadPSF + psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); + return psphotReadoutCleanup (config, view, filerule); + } + + // include externally-supplied sources (loaded onto detections->newSources) + psphotLoadExtSources (config, view, filerule); + + // merge the newly selected sources into the existing list (detections->allSources) + // NOTE: merge OLD and NEW + psphotMergeSources (config, view, filerule); + + // construct sources and measure moments and other basic stats (saved on detections->allSources) + // all sources use the auto-scaled window appropriate to a PSF, except for the saturated + // stars : these use a larger window (3x the basic window) + if (!psphotFullForceSourceStats (config, view, filerule, true)) { // pass 1 + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); + return psphotReadoutCleanup (config, view, filerule); + } + + // generate a psf model for any readouts which need one + // psphotFullForcePSF (config, view, filerule); + // XXX I'm not sure we need a different algorithm here or not : we are supplying the + // sources; we should mark with a flag bit the ones we actually want to use as PSF + // stars (this means we need to supply this info in the load). + psphotChoosePSF (config, view, filerule, false); + + // Construct an initial model for each object, set the radius to fitRadius, set circular + // fit mask. NOTE: only applied to sources without guess models + // XXX this currently only generates a PSF model + psphotGuessModels (config, view, filerule); + + // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) + // XXX this will fit non-PSF sources if they are supplied + psphotFitSourcesLinear (config, view, filerule, false, false); + + // measure moments + // XXX does this set the correct gaussian window (how to do multiple windows?) + // we are calling this above... + // psphotSourceStats (config, view, filerule, false); + + // measure kron fluxes + psphotKronFlux (config, view, filerule); + + // measure petro fluxes + psphotPetroFlux (config, view, filerule); + + // measure radial apertures + psphotRadialApertures (config, view, filerule, 0); + + // measure galaxy shapes + psphotGalaxyShape (config, view, filerule); + + // calculate source magnitudes (psf mag and ap mag) + psphotMagnitudes(config, view, filerule); + + // replace background in residual image + psphotSkyReplace (config, view, filerule); + + // drop the references to the image pixels held by each source + psphotSourceFreePixels (config, view, filerule); + + // create the exported-metadata and free local data + return psphotReadoutCleanup (config, view, filerule); +} Index: trunk/psphot/src/psphotGalaxyShape.c =================================================================== --- trunk/psphot/src/psphotGalaxyShape.c (revision 36375) +++ trunk/psphot/src/psphotGalaxyShape.c (revision 36375) @@ -0,0 +1,400 @@ +# include "psphotInternal.h" + +bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Galaxy Shapes ---"); + + // 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? + // XXX for this to be true for psphotFullForce?? + if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping galaxy shape measurements\n"); + return true; + } + + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + + // 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?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->newSources ? detections->newSources : detections->allSources; + psAssert (sources, "missing sources?"); + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + // psAssert (psf, "missing psf?"); + + if (!psphotGalaxyShapeReadout (config, recipe, view, filerule, readout, sources, psf)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf) { + + bool status = false; + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping galaxy shapes"); + return true; + } + + psTimerStart ("psphot.galaxy"); + + psphotInitRadiusEXT (recipe, readout); + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // bit-mask to mark pixels not used in analysis + psphotGalaxyShapeOptions *opt = psphotGalaxyShapeOptionsAlloc(); + opt->fRmajorMin = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MIN"); if (!status) opt->fRmajorMin = 0.5; + opt->fRmajorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MAX"); if (!status) opt->fRmajorMax = 2.0; + opt->fRmajorDel = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_DEL"); if (!status) opt->fRmajorDel = 0.1; + opt->fRminorMin = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_MIN"); if (!status) opt->fRminorMin = 0.5; + opt->fRminorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_MAX"); if (!status) opt->fRminorMax = 2.0; + opt->fRminorDel = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_DEL"); if (!status) opt->fRminorDel = 0.1; + + // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) + // NOTE: this is only used if we are NOT smoothing with a 1D Gaussian + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + assert (status); + + float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolutio + if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) { + fitNsigmaConv = 5.0; + } + + // Define source fitting parameters for extended source fits + // we are not doing LMM fitting, so most options are irrelevant + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->mode = PM_SOURCE_FIT_EXT_AND_SKY; + fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + fitOptions->nsigma = fitNsigmaConv; + + // Poisson or Constant weight for chisq tests? + fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON"); + if (!status) fitOptions->poissonErrors = true; + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortByFlux); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + return true; + } + + // threaded measurement of the source magnitudes + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_GALAXY_SHAPES"); + + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, fitOptions); + psArrayAdd(job->args, 1, opt); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32); + +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } +# else + if (!psphotGalaxyShape_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + // psFree(AnalysisRegion); + return false; + } + psFree(job); +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } + psFree(job); + } + } + psFree (cellGroups); + psFree (fitOptions); + psFree (opt); + + psLogMsg ("psphot.galaxy", PS_LOG_WARN, "measure galaxy shapes : %f sec for %ld objects\n", psTimerMark ("psphot.galaxy"), sources->n); + return true; +} + +bool psphotGalaxyShape_Threaded (psThreadJob *job) { + + pmReadout *readout = job->args->data[0]; + psArray *sources = job->args->data[1]; + pmSourceFitOptions *fitOptions = job->args->data[2]; + psphotGalaxyShapeOptions *opt = job->args->data[3]; + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); + int psfSize = PS_SCALAR_VALUE(job->args->data[6],S32); + + float fitRadius; + float windowRadius; + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + if (!source->peak) continue; // XXX how can we have a peak-less source? + + // check status of this source's moments + if (!source->moments) continue; + if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + + // psphotSetRadiusMomentsExact sets the radius based on Mrf + if (!isfinite(source->moments->Mrf)) continue; + + // modelFits is allocated if a galaxy fit is requested + if (!source->modelFits) continue; + + psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs + + // skip saturated stars modeled with a radial profile + // XXX worry about this at some point.. + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + // this function populates moments->Mrf,GalaxyShape,GalaxyShapeErr + // do the following for a set of shapes (Ex,Ey) + psphotGalaxyShapeGrid (source, fitOptions, opt, maskVal, psfSize); + + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + + return true; +} + +bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize) { + + // we use the 0th model to define the initial guess shape + pmModel *model = source->modelFits->data[0]; + if (!model) return false; + + pmModelType modelType = model->type; + + // we are using fitOptions->mode : be sure this makes sense + pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); + if (!pcm) return false; + + // we are fitting only PM_PAR_I0; the shape elements are generated from a grid + psF32 *PAR = pcm->modelConv->params->data.F32; + + // double check that the guess is carried along... + + // I have some source guess (e0, e1, e2) + psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType); + + if (!source->galaxyFits) { + source->galaxyFits = pmSourceGalaxyFitsAlloc(); + } + + float fRmajorBest = NAN; + float fRminorBest = NAN; + float chisqBest = NAN; + for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; fRmajor += opt->fRmajorDel) { + for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) { + + psEllipseAxes testAxes = guessAxes; + testAxes.major = guessAxes.major * fRmajor; + testAxes.minor = guessAxes.minor * fRminor; + + pmPSF_AxesToModel (PAR, testAxes, modelType); + + psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, true); + + int i = source->galaxyFits->chisq->n - 1; + float thisChisq = source->galaxyFits->chisq->data.F32[i]; + if (isfinite(thisChisq) && (!isfinite(chisqBest) || thisChisq < chisqBest)) { + chisqBest = thisChisq; + fRmajorBest = fRmajor; + fRminorBest = fRminor; + } + } + } + + if (isfinite(chisqBest)) { + // now save the best fitting model as the source's extended model + psEllipseAxes testAxes = guessAxes; + testAxes.major = guessAxes.major * fRmajorBest; + testAxes.minor = guessAxes.minor * fRminorBest; + + pmPSF_AxesToModel (PAR, testAxes, modelType); + + psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false); + + psFree (source->modelEXT); + + source->modelEXT = psMemIncrRefCounter (pcm->modelConv); + source->type = PM_SOURCE_TYPE_EXTENDED; + source->mode |= PM_SOURCE_MODE_EXTMODEL; + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; + + // adjust the window so the subtraction covers the faint wings + // psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal); + + // cache the model flux + pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma); + } + + psFree (pcm); + return true; +} + +// fit the given model to the source and find chisq & normalization +// XXX is this a single-component model? sersic with a supplied index, Reff, axis ratio, and theta? +bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults) { + + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + + // generate the modelFlux + // need to reset here each time since we assign below to calculate the flux + // XXX note that this does not add sky to model + pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize); + pcm->modelConv->params->data.F32[PM_PAR_I0] = 1.0; + + int nPix = 0; + float YY = 0.0; + float YM = 0.0; + float MM = 0.0; + bool usePoisson = pcm->poissonErrors; + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + // skip masked points + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + continue; + } + // skip zero-variance points + if (source->variance->data.F32[iy][ix] == 0) { + continue; + } + + // skip nan value points + if (!isfinite(source->pixels->data.F32[iy][ix])) { + continue; + } + + float fy = source->pixels->data.F32[iy][ix]; + float fm = source->modelFlux->data.F32[iy][ix]; + float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0; + + YY += PS_SQR(fy) * wt; + YM += fm * fy * wt; + MM += PS_SQR(fm) * wt; + nPix ++; + } + } + + float Io = YM / MM; + float dIo = sqrt (1.0 / MM); + float Chisq = (YY - 2 * Io * YM + Io * Io * MM) / (float) nPix; + // NOTE : if !poissonErrors, Chisq is not really the chisq, but is scaled by the flux. + + pcm->modelConv->params->data.F32[PM_PAR_I0] = Io; + float flux = pcm->modelConv->modelFlux (pcm->modelConv->params); + float dflux = flux * (dIo / Io); + + + if (saveResults) { + psVectorAppend (source->galaxyFits->Flux, flux); + psVectorAppend (source->galaxyFits->dFlux, dflux); + psVectorAppend (source->galaxyFits->chisq, Chisq); + source->galaxyFits->nPix = nPix; + } + + return true; +} + +/**** support functions ****/ + +void psphotGalaxyShapeOptionsFree (psphotGalaxyShapeOptions *opt) { + return; +} + +psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc() +{ + psphotGalaxyShapeOptions *opt = (psphotGalaxyShapeOptions *) psAlloc(sizeof(psphotGalaxyShapeOptions)); + psMemSetDeallocator(opt, (psFreeFunc) psphotGalaxyShapeOptionsFree); + + opt->fRmajorMin = 0.5; + opt->fRmajorMax = 2.0; + opt->fRmajorDel = 0.1; + opt->fRminorMin = 0.5; + opt->fRminorMax = 2.0; + opt->fRminorDel = 0.1; + + return opt; +} + Index: trunk/psphot/src/psphotImageLoop.c =================================================================== --- trunk/psphot/src/psphotImageLoop.c (revision 36366) +++ trunk/psphot/src/psphotImageLoop.c (revision 36375) @@ -132,4 +132,11 @@ } break; + case PSPHOT_FULL_FORCE: + if (!psphotFullForceReadout (config, view, "PSPHOT.INPUT")) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + break; case PSPHOT_MAKE_PSF: if (!psphotMakePSFReadout (config, view, "PSPHOT.INPUT")) { Index: trunk/psphot/src/psphotKronFlux.c =================================================================== --- trunk/psphot/src/psphotKronFlux.c (revision 36375) +++ trunk/psphot/src/psphotKronFlux.c (revision 36375) @@ -0,0 +1,250 @@ +# include "psphotInternal.h" + +bool psphotKronFlux (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Kron Fluxes ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + + // 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?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->newSources ? detections->newSources : detections->allSources; + psAssert (sources, "missing sources?"); + + if (!psphotKronFluxReadout (config, recipe, view, filerule, readout, sources)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); + return false; + } + } + + return true; +} + +bool psphotKronFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources) { + + bool status = false; + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping kron fluxes"); + return true; + } + + psTimerStart ("psphot.kron"); + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortByFlux); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + return true; + } + + // threaded measurement of the source magnitudes + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_FLUX"); + + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, cells->data[j]); // sources + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } +# else + if (!psphotKronFlux_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + // psFree(AnalysisRegion); + return false; + } + psFree(job); +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } + psFree(job); + } + } + psFree (cellGroups); + + psLogMsg ("psphot.kron", PS_LOG_WARN, "measure kron fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); + return true; +} + +bool psphotKronFlux_Threaded (psThreadJob *job) { + + pmReadout *readout = job->args->data[0]; + psArray *sources = job->args->data[1]; + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA); + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + if (!source->peak) continue; // XXX how can we have a peak-less source? + + // check status of this source's moments + if (!source->moments) continue; + if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + + // XXX where are we storing the supplied kron radius? + if (!isfinite(source->moments->Mrf)) continue; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + // is this right?? + float windowRadius = 2.5*source->moments->Mrf ; + + // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS + pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); + psAssert (source->pixels, "WTF?"); + + // this function populates moments->Mrf,KronFlux,KronFluxErr + psphotKronFluxSource (source, maskVal); + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + + return true; +} + +// measure just the flux in the given aperture (this is probably a single common function -- can I use one of the pmSourcePhotometry functions?) +bool psphotKronFluxSource (pmSource *source, psImageMaskType maskVal) { + + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + + // the peak position is less accurate but less subject to extreme deviations + // XXX do this or just use the peak position? + float dX = source->moments->Mx - source->peak->xf; + float dY = source->moments->My - source->peak->yf; + float dR = hypot(dX, dY); + float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf; + float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf; + + // center of mass in subimage. Note: the calculation below uses pixel index, so we correct + // xCM, yCM from pixel coords to pixel index here. + psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage + psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage + + // Calculate the Kron magnitude (make this block optional?) + // XXX set the aperture here + float radKron = 2.5*source->moments->Mrf; + float radKron2 = radKron*radKron; + + int nKronPix = 0; + float Sum = 0.0; + float Var = 0.0; + + // set vPix to the source pixels (it may have been set to the + // smoothed image above) + psF32 **vPix = source->pixels->data.F32; + psF32 **vWgt = source->variance->data.F32; + psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 yDiff = row - yCM; + if (fabs(yDiff) > radKron) continue; + + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + // check mask and value for this pixel + if (vMsk && (vMsk[row][col] & maskVal)) continue; + if (isnan(vPix[row][col])) continue; + + psF32 xDiff = col - xCM; + if (fabs(xDiff) > radKron) continue; + + // radKron is just a function of (xDiff, yDiff) + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > radKron2) continue; + + float pDiff = vPix[row][col]; + psF32 wDiff = vWgt[row][col]; + + Sum += pDiff; + Var += wDiff; + nKronPix ++; + } + } + + // return the flux and error to parameters? + source->moments->KronFlux = Sum; + source->moments->KronFluxErr = sqrt(Var); + + return true; +} Index: trunk/psphot/src/psphotKronIterate.c =================================================================== --- trunk/psphot/src/psphotKronIterate.c (revision 36366) +++ trunk/psphot/src/psphotKronIterate.c (revision 36375) @@ -166,5 +166,7 @@ } +# if (PS_TRACE_ON) fprintf (stderr, "--- starting KRON ---\n"); +#endif // We measure the Kron Radius on a smoothed copy of the readout image Index: trunk/psphot/src/psphotLoadPSF.c =================================================================== --- trunk/psphot/src/psphotLoadPSF.c (revision 36366) +++ trunk/psphot/src/psphotLoadPSF.c (revision 36375) @@ -1,3 +1,20 @@ # include "psphotInternal.h" + +// PSPHOT.PSF.LOAD vs input file -- see note at top +bool psphotLoadPSF (pmConfig *config, const pmFPAview *view, const char *filerule) { + + int num = psphotFileruleCount(config, filerule); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + + // Generate the mask and weight images, including the user-defined analysis region of interest + if (!psphotLoadPSFReadout (config, view, filerule, "PSPHOT.PSF.LOAD", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to load PSF model for PSPHOT.PSF.LOAD entry %d", i); + return false; + } + } + return true; +} // NOTE : pmPSF_IO.c functions must load the psf model onto the chip->analysis metadata because @@ -58,19 +75,2 @@ return true; } - -// PSPHOT.PSF.LOAD vs input file -- see note at top -bool psphotLoadPSF (pmConfig *config, const pmFPAview *view, const char *filerule) { - - int num = psphotFileruleCount(config, filerule); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - - // Generate the mask and weight images, including the user-defined analysis region of interest - if (!psphotLoadPSFReadout (config, view, filerule, "PSPHOT.PSF.LOAD", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to load PSF model for PSPHOT.PSF.LOAD entry %d", i); - return false; - } - } - return true; -} Index: trunk/psphot/src/psphotMergeSources.c =================================================================== --- trunk/psphot/src/psphotMergeSources.c (revision 36366) +++ trunk/psphot/src/psphotMergeSources.c (revision 36375) @@ -39,12 +39,4 @@ psAssert (newSources, "missing sources?"); - // XXX TEST: - if (detections->allSources) { - psphotMaskCosmicRayFootprintCheck(detections->allSources); - } - if (detections->newSources) { - psphotMaskCosmicRayFootprintCheck(detections->newSources); - } - if (!detections->allSources) { detections->allSources = psArrayAllocEmpty(newSources->n); @@ -65,4 +57,7 @@ // Merge the externally supplied sources with the existing sources. Mark them as having mode // PM_SOURCE_MODE_EXTERNAL. + +// XXX this function needs to be updated slightly for psphotFullForce: +// * load the additional parameters to guide the new concepts // XXX This function needs to be updated to loop over set of input files. At the moment, we @@ -73,4 +68,5 @@ bool status; pmDetections *extCMF = NULL; + pmDetections *extCFF = NULL; psArray *extSourcesTXT = NULL; int index = 0; @@ -102,5 +98,5 @@ { pmReadout *readoutCMF = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT.CMF"); - if (!readoutCMF) goto loadTXT; + if (!readoutCMF) goto loadCFF; extCMF = psMetadataLookupPtr (NULL, readoutCMF->analysis, "PSPHOT.DETECTIONS"); @@ -125,4 +121,37 @@ } +loadCFF: + // load data from input CFF file: + { + pmReadout *readoutCFF = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT.CFF"); + if (!readoutCFF) goto loadTXT; + + extCFF = psMetadataLookupPtr (NULL, readoutCFF->analysis, "PSPHOT.DETECTIONS"); + if (extCFF) { + psF32 exptime = psMetadataLookupF32(NULL, readout->parent->concepts, "CELL.EXPOSURE"); + for (int i = 0; i < extCFF->allSources->n; i++) { + pmSource *source = extCFF->allSources->data[i]; + + // setting this bit not only tracks the inputs, it makes pmSourceMoments + // keep the Mx,My values for the centroid. + source->mode |= PM_SOURCE_MODE_EXTERNAL; + + // source->peak->detValue,rawFlux,smoothFlux all set to input flux value which is scaled + // to 1 second exposure time. Scale to this image's exposure. + source->peak->rawFlux *= exptime; + source->peak->smoothFlux *= exptime; + source->peak->detValue *= exptime; + // source->peak->xf,yf, moments->Mx,My all set to input position + + // drop the loaded source modelPSF + psFree (source->modelPSF); + source->modelPSF = NULL; + source->imageID = index; + + psArrayAdd (detections->newSources, 100, source); + } + } + } + loadTXT: @@ -157,5 +186,5 @@ psFree (detections); - if (!extCMF && !extSourcesTXT) { + if (!(extCMF || extCFF || extSourcesTXT)) { psLogMsg ("psphot", 3, "no external sources for this readout"); return true; @@ -163,8 +192,9 @@ int nCMF = extCMF ? extCMF->allSources->n : 0; - int nTXT = extSourcesTXT ? extSourcesTXT->n : 0; - - psLogMsg ("psphot", 3, "%d external sources (%d cmf, %d text) merged to yield %ld total sources", - nCMF + nTXT, nCMF, nTXT, sources->n); + int nCFF = extCFF ? extCFF->allSources->n : 0; + int nTXT = extSourcesTXT ? extSourcesTXT->n : 0; + + psLogMsg ("psphot", 3, "%d external sources (%d cmf, %d cff, %d text) merged to yield %ld total sources", + nCMF + nCFF + nTXT, nCMF, nCFF, nTXT, sources->n); return true; } Index: trunk/psphot/src/psphotModelTestReadout.c =================================================================== --- trunk/psphot/src/psphotModelTestReadout.c (revision 36366) +++ trunk/psphot/src/psphotModelTestReadout.c (revision 36375) @@ -125,5 +125,5 @@ // get the source moments - status = pmSourceMoments (source, radius, 0.25*radius, 0.0, MIN_KRON_RADIUS, maskVal); + status = pmSourceMoments (source, radius, 0.25*radius, 0.0, MIN_KRON_RADIUS, maskVal); if (!status) psAbort("psSourceMoments error"); Index: trunk/psphot/src/psphotOldCode.c =================================================================== --- trunk/psphot/src/psphotOldCode.c (revision 36375) +++ trunk/psphot/src/psphotOldCode.c (revision 36375) @@ -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; +} + Index: trunk/psphot/src/psphotPetroFlux.c =================================================================== --- trunk/psphot/src/psphotPetroFlux.c (revision 36375) +++ trunk/psphot/src/psphotPetroFlux.c (revision 36375) @@ -0,0 +1,259 @@ +# include "psphotInternal.h" + +# define PETROSIAN_RADII 2.0 + +bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Petro Fluxes ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + psMetadataAddBool (recipe, PS_LIST_TAIL, "EXTENDED_SOURCE_ANALYSIS", PS_META_REPLACE, "we measured this, save to disk", true); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + + // 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?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->newSources ? detections->newSources : detections->allSources; + psAssert (sources, "missing sources?"); + + if (!psphotPetroFluxReadout (config, recipe, view, filerule, readout, sources)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources) { + + bool status = false; + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping petrosian fluxes"); + return true; + } + + psTimerStart ("psphot.petro"); + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortByFlux); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + return true; + } + + // threaded measurement of the source magnitudes + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_PETRO_FLUX"); + + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, cells->data[j]); // sources + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } +# else + if (!psphotPetroFlux_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + // psFree(AnalysisRegion); + return false; + } + psFree(job); +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } + psFree(job); + } + } + psFree (cellGroups); + + psLogMsg ("psphot.petro", PS_LOG_WARN, "measure petro fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.petro"), sources->n); + return true; +} + +bool psphotPetroFlux_Threaded (psThreadJob *job) { + + pmReadout *readout = job->args->data[0]; + psArray *sources = job->args->data[1]; + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA); + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + if (!source->peak) continue; // XXX how can we have a peak-less source? + if (!source->extpars) continue; // if this is not set, we did not have a valid petRadius + + // check status of this source's moments + if (!source->moments) continue; + if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + + // XXX where are we storing the supplied petro radius? + if (!isfinite(source->moments->Mrf)) continue; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + // is the right?? XXX need the correct location for Rpet + float windowRadius = 2.0*source->moments->Mrf ; + + // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS + pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); + psAssert (source->pixels, "WTF?"); + + // this function populates moments->Mrf,PetroFlux,PetroFluxErr + psphotPetroFluxSource (source, maskVal); + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + + return true; +} + +// measure just the flux in the given aperture (this is probably a single common function -- can I use one of the pmSourcePhotometry functions?) +bool psphotPetroFluxSource (pmSource *source, psImageMaskType maskVal) { + + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->variance, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->extpars, false); + + // the peak position is less accurate but less subject to extreme deviations + float dX = source->moments->Mx - source->peak->xf; + float dY = source->moments->My - source->peak->yf; + float dR = hypot(dX, dY); + float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf; + float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf; + + // center of mass in subimage. Note: the calculation below uses pixel index, so we correct + // xCM, yCM from pixel coords to pixel index here. + psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage + psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage + + // Calculate the Petro magnitude (make this block optional?) + // XXX set the aperture here + float radPetro = PETROSIAN_RADII*source->extpars->petrosianRadius; + float radPetro2 = radPetro*radPetro; + + int nPetroPix = 0; + float Sum = 0.0; + float Var = 0.0; + + // set vPix to the source pixels (it may have been set to the + // smoothed image above) + psF32 **vPix = source->pixels->data.F32; + psF32 **vWgt = source->variance->data.F32; + psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 yDiff = row - yCM; + if (fabs(yDiff) > radPetro) continue; + + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + // check mask and value for this pixel + if (vMsk && (vMsk[row][col] & maskVal)) continue; + if (isnan(vPix[row][col])) continue; + + psF32 xDiff = col - xCM; + if (fabs(xDiff) > radPetro) continue; + + // radPetro is just a function of (xDiff, yDiff) + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > radPetro2) continue; + + float pDiff = vPix[row][col]; + psF32 wDiff = vWgt[row][col]; + + Sum += pDiff; + Var += wDiff; + nPetroPix ++; + } + } + + // return the flux and error to parameters? + source->extpars->petrosianFlux = Sum; + source->extpars->petrosianFluxErr = sqrt(Var); + source->extpars->petrosianFill = nPetroPix / (M_PI * radPetro2); + + // fprintf (stderr, "petro flux: %f +/- %f\n", Sum, sqrt(Var)); + + return true; +} Index: trunk/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- trunk/psphot/src/psphotPetrosianRadialBins.c (revision 36366) +++ trunk/psphot/src/psphotPetrosianRadialBins.c (revision 36375) @@ -186,4 +186,5 @@ psFree(values); psFree(stats); + source->mode2 |= PM_SOURCE_MODE2_RADBIN_NAN_CENTER; return true; } Index: trunk/psphot/src/psphotPetrosianStats.c =================================================================== --- trunk/psphot/src/psphotPetrosianStats.c (revision 36366) +++ trunk/psphot/src/psphotPetrosianStats.c (revision 36375) @@ -6,4 +6,5 @@ // generate the Petrosian radius and flux from the mean surface brightness (r_i) +float InterpolateValuesQuadratic (float *Xin, float *Yin, float X); float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1); @@ -20,4 +21,5 @@ if (!profile->binSB) { psLogMsg ("psphot", PS_LOG_DETAIL, "no petrosian profile, skipping source %f, %f", source->peak->xf, source->peak->yf); + source->mode2 |= PM_SOURCE_MODE2_PETRO_NO_PROFILE; return true; } @@ -53,5 +55,5 @@ float dFsum2 = 0.0; - float nSigma = 3.0; + float nSigma = 2.0; int lowestSignificantRadius = 0; float lowestSignificantRatio = 1.0; @@ -117,7 +119,19 @@ petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[nOut]); - } else { - petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); - petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]); + source->mode2 |= PM_SOURCE_MODE2_PETRO_RATIO_ZEROBIN; + } else { + // petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + if (nOut > 1) { + petRadius = InterpolateValuesQuadratic (&petRatio->data.F32[nOut-2], &refRadius->data.F32[nOut-2], PETROSIAN_RATIO); + } else { + petRadius = InterpolateValuesQuadratic (&petRatio->data.F32[nOut-3], &refRadius->data.F32[nOut-3], PETROSIAN_RATIO); + } +# if (PS_TRACE_ON) + float petRadiusLinear = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + if (fabs(petRadius - petRadiusLinear) > fabs(refRadius->data.F32[nOut] - refRadius->data.F32[nOut-1])) { + fprintf (stderr, "big difference : %f vs %f\n", petRadius, petRadiusLinear); + } +#endif + petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]); } above = false; @@ -148,4 +162,5 @@ fprintf (stderr, "nan pet radius\n"); } + source->mode2 |= PM_SOURCE_MODE2_PETRO_INSIG_RATIO; } @@ -180,4 +195,5 @@ bool found50 = false; bool found90 = false; + // XXX use bisection to do this faster: for (int i = 0; !(found50 && found90) && i < refRadius->n; i++) { @@ -235,4 +251,23 @@ } +// Lagrange's form of the interpolating polynomial... +float InterpolateValuesQuadratic (float *Xin, float *Yin, float X) { + + float dx01 = Xin[0] - Xin[1]; + float dx02 = Xin[0] - Xin[2]; + float dx12 = Xin[1] - Xin[2]; + + float dx0 = X - Xin[0]; + float dx1 = X - Xin[1]; + float dx2 = X - Xin[2]; + + float y0 = Yin[0]*dx1*dx2/(dx01*dx02); + float y1 = Yin[1]*dx0*dx2/(dx01*dx12); // need - sign + float y2 = Yin[2]*dx0*dx1/(dx02*dx12); + + float Y = y0 - y1 + y2; + return Y; +} + float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) { float dydx = (Y1 - Y0) / (X1 - X0); Index: trunk/psphot/src/psphotRadialBins.c =================================================================== --- trunk/psphot/src/psphotRadialBins.c (revision 36366) +++ trunk/psphot/src/psphotRadialBins.c (revision 36375) @@ -44,12 +44,8 @@ psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - if (!radMin || !radMin->n) { - psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMin missing or empty)"); - return false; - } - if (!radMax || !radMax->n) { - psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMax missing or empty)"); - return false; - } + psAssert (radMin, "RADIAL.ANNULAR.BINS.LOWER is missing from recipe"); + psAssert (radMin->n, "RADIAL.ANNULAR.BINS.LOWER is empty in recipe"); + psAssert (radMax, "RADIAL.ANNULAR.BINS.UPPER is missing from recipe"); + psAssert (radMax->n, "RADIAL.ANNULAR.BINS.UPPER is empty in recipe"); psVector *binSB = psVectorAllocEmpty(radMin->n, PS_TYPE_F32); // surface brightness of radial bin @@ -161,4 +157,5 @@ psFree(values); psFree(stats); + source->mode2 |= PM_SOURCE_MODE2_RADBIN_NAN_CENTER; return true; } Index: trunk/psphot/src/psphotRadialProfile.c =================================================================== --- trunk/psphot/src/psphotRadialProfile.c (revision 36366) +++ trunk/psphot/src/psphotRadialProfile.c (revision 36375) @@ -1,3 +1,11 @@ # include "psphotInternal.h" + +static int Nskip1 = 0; +static int Nskip2 = 0; +static int Nskip3 = 0; +static int Nskip4 = 0; +static int Nskip5 = 0; + +# define SKIP(VALUE) { VALUE++; return false; } bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) { @@ -23,5 +31,5 @@ if (!psphotRadialProfilesByAngles (source, Nsec, Rmax)) { psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian"); - return false; + SKIP (Nskip1); } // allocate: extpars->radFlux->radii,fluxes,theta @@ -32,5 +40,5 @@ if (!psphotRadiiFromProfiles (source, fluxMin, fluxMax)) { psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles"); - return false; + SKIP (Nskip2); } // allocate : extpars->radFlux->isophotalRadii (use profile->radii,fluxes) @@ -40,5 +48,5 @@ if (!psphotEllipticalContour (source)) { // psLogMsg ("psphot", 3, "failed to measure elliptical contour"); - return false; + SKIP (Nskip3); } // use extpars->radFlux->isophotalRadii,theta (result in extpars->axes) @@ -48,5 +56,5 @@ if (!psphotEllipticalProfile (source, RAW_RADIUS)) { psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); - return false; + SKIP (Nskip4); } // allocate extpars->ellipticalFlux->radiusElliptical,fluxElliptical (use axes to scale raw pixels) @@ -55,5 +63,5 @@ if (!psphotRadialBins (recipe, source, Rmax, skynoise)) { psError (PS_ERR_UNKNOWN, false, "failed to generate radial bins"); - return false; + SKIP (Nskip5); } // allocate extpars->radProfile->binSB, binSBstdv, binSum, binFill, radialBins, area (small lengths) @@ -62,2 +70,12 @@ return true; } + +void psphotRadialProfileShowSkips () { +# if (PS_TRACE_ON) + fprintf (stderr, "radial profile skipped @ 1 : %d\n", Nskip1); + fprintf (stderr, "radial profile skipped @ 2 : %d\n", Nskip2); + fprintf (stderr, "radial profile skipped @ 3 : %d\n", Nskip3); + fprintf (stderr, "radial profile skipped @ 4 : %d\n", Nskip4); + fprintf (stderr, "radial profile skipped @ 5 : %d\n", Nskip5); +#endif +} Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 36366) +++ trunk/psphot/src/psphotReadout.c (revision 36375) @@ -316,4 +316,11 @@ psphotSourceSize (config, view, filerule, false); // pass 2 (detections->allSources) + // XXX currently we are doing both the analysis of the size and the assessment of "fit ext" + // in source size. this overloads the bit MODE_EXT_LIMIT to mean "fit ext" not just + // "bigger than a PSF" + + // decide which source(s) are to be fitted with the extended source analysis code. + psphotChooseAnalysisOptions (config, view, filerule); + psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources) psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources) Index: trunk/psphot/src/psphotRoughClass.c =================================================================== --- trunk/psphot/src/psphotRoughClass.c (revision 36366) +++ trunk/psphot/src/psphotRoughClass.c (revision 36375) @@ -87,5 +87,5 @@ if (NX > NY) { NXuse = ScaleForClump; - NYuse = (int) (ScaleForClump * (NY / NX) + 0.5); + NYuse = (int) (ScaleForClump * (NX / NY) + 0.5); } else { NYuse = ScaleForClump; Index: trunk/psphot/src/psphotSetThreads.c =================================================================== --- trunk/psphot/src/psphotSetThreads.c (revision 36366) +++ trunk/psphot/src/psphotSetThreads.c (revision 36375) @@ -47,4 +47,19 @@ psFree(task); + task = psThreadTaskAlloc("PSPHOT_KRON_FLUX", 4); + task->function = &psphotKronFlux_Threaded; + psThreadTaskAdd(task); + psFree(task); + + task = psThreadTaskAlloc("PSPHOT_PETRO_FLUX", 4); + task->function = &psphotPetroFlux_Threaded; + psThreadTaskAdd(task); + psFree(task); + + task = psThreadTaskAlloc("PSPHOT_GALAXY_SHAPES", 7); + task->function = &psphotGalaxyShape_Threaded; + psThreadTaskAdd(task); + psFree(task); + task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 10); task->function = &psphotBlendFit_Threaded; Index: trunk/psphot/src/psphotSourceFits.c =================================================================== --- trunk/psphot/src/psphotSourceFits.c (revision 36366) +++ trunk/psphot/src/psphotSourceFits.c (revision 36375) @@ -19,4 +19,7 @@ static int NfitIterPCM = 0; static int NfitPixPCM = 0; + +bool psphotPCMfitCheckSize (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize); +bool psphotPCMfitRetry (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, float psfSize); bool psphotFitInit (int nThreads) { @@ -579,12 +582,10 @@ maskVal |= markVal; - // allocate the model + // allocate the model (this can only fail on a config error) pmModel *model = pmModelAlloc(modelType); - if (!model) { - return NULL; - } - - float t1, t2, t4, t5; - t1 = t2 = t4 = t5 = 0.0; + psAssert (model, "invalid extended model name"); + + float t1, t2, t3, t4, t5; + t1 = t2 = t3 = t4 = t5 = 0.0; if (TIMING) { psTimerStart ("psphotFitPCM"); } @@ -619,5 +620,5 @@ if (!psphotSersicModelGuessPCM (pcm, source, maskVal, psfSize)) { psFree (pcm); - model->flags |= PM_MODEL_STATUS_BADARGS; + model->flags |= PM_MODEL_SERSIC_PCM_FAIL_GUESS; return model; } @@ -628,5 +629,5 @@ if (!psphotFitSersicShapeAndIndexGridAuto (pcm, readout, source, &options, maskVal, markVal, psfSize)) { psFree (pcm); - model->flags |= PM_MODEL_STATUS_BADARGS; + model->flags |= PM_MODEL_SERSIC_PCM_FAIL_GRID; psError(PS_ERR_UNKNOWN, true, "Failed to find a index & shape"); psErrorClear (); // clear the polynomial error @@ -637,10 +638,10 @@ if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) { psFree (pcm); - model->flags |= PM_MODEL_STATUS_BADARGS; + model->flags |= PM_MODEL_PCM_FAIL_GUESS; return model; } } - if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); } + if (TIMING) { t3 = psTimerMark ("psphotFitPCM"); } // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); @@ -648,9 +649,18 @@ NfitIterPCM += pcm->modelConv->nIter; NfitPixPCM += pcm->modelConv->nDOF; + if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); } + + // XXX we might make this more efficient by setting NITER to be fairly small. if we hit the iteration + // limit, then we could do a small grid search on the size and try again from the best fit + + if (options.isInteractive) psphotPCMfitCheckSize (pcm, source, maskVal, psfSize); + // if (pcm->modelConv->nIter == fitOptions->nIter) { + // psphotPCMfitRetry (pcm, source, &options, maskVal, markVal, psfSize); + // } if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); } if (TIMING) { int nPixBig = source->pixels->numCols * source->pixels->numRows; - fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5); + fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t3: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t3, t4, t5); } if (EXTRA_VERBOSE && !TIMING) { @@ -778,10 +788,14 @@ } -// float indexGuessInv[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083}; -float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125}; -# define N_INDEX_GUESS_INV 7 - -// float reffGuess[] = {3.0, 10.0, 20.0, 30.0, 40.0}; -float reffGuess[] = {0.5, 0.75, 1.0, 1.4, 2.0}; +// 0.5 / n for (1.0, 1.25, 1.66, 2.0, 3.33, 4.0) +// float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125}; + +// 0.5 / n for (0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0) +float indexGuessInv[] = {1.00, 0.50, 0.333, 0.25, 0.166, 0.125, 0.10, 0.0833}; +float indexGuessR1q[] = {1.06, 1.19, 1.335, 1.48, 1.840, 2.290, 2.84, 3.5300}; +# define N_INDEX_GUESS_INV 8 + +// we are going to guess in fractions about the R1-based guess +float reffGuess[] = {0.8, 0.9, 1.0, 1.12, 1.25}; # define N_REFF_GUESS 5 @@ -791,5 +805,5 @@ // we get a reasonable guess from: - // * Reff = Kron R1 + // * Reff = Kron R1 / Q(index) -- Q comes from Graham & Driver // * Rmajor / Rminor & Theta from moments // * Io from total Kron flux @@ -813,18 +827,4 @@ psEllipseAxes momentAxes = psEllipseMomentsToAxes (moments, 20.0); - if (0) { - psEllipseAxes guessAxes; - guessAxes.major = source->moments->Mrf; - guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major; - guessAxes.theta = momentAxes.theta; - - if (!isfinite(guessAxes.major)) return false; - if (!isfinite(guessAxes.minor)) return false; - if (!isfinite(guessAxes.theta)) return false; - - // convert the major,minor,theta to shape parameters for an Reff-like model - pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); - } - // set the model position if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { @@ -846,22 +846,22 @@ // loop over index guesses and find the best fit for (int j = 0; j < N_REFF_GUESS; j++) { - - psEllipseAxes guessAxes; - guessAxes.major = reffGuess[j] * source->moments->Mrf; - guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major); - guessAxes.theta = momentAxes.theta; - - if (!isfinite(guessAxes.major)) return false; - if (!isfinite(guessAxes.minor)) return false; - if (!isfinite(guessAxes.theta)) return false; - - // convert the major,minor,theta to shape parameters for an Reff-like model - pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); - for (int i = 0; i < N_INDEX_GUESS_INV; i++) { PAR[PM_PAR_7] = indexGuessInv[i]; + psEllipseAxes guessAxes; + guessAxes.major = reffGuess[j] * source->moments->Mrf / indexGuessR1q[i]; + guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major); + guessAxes.theta = momentAxes.theta; + + if (!isfinite(guessAxes.major)) return false; + if (!isfinite(guessAxes.minor)) return false; + if (!isfinite(guessAxes.theta)) return false; + + // convert the major,minor,theta to shape parameters for an Reff-like model + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); + // generated the modelFlux - pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize); + // XXX note that this does not add sky to model + pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize); float YY = 0.0; @@ -901,8 +901,10 @@ iMin = Io; sMin = indexGuessInv[i]; - rMin = reffGuess[j]; + rMin = reffGuess[j] / indexGuessR1q[i]; } - // fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin); - // fprintf (stderr, "\n"); + if (EXTRA_VERBOSE) { + fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin); + fprintf (stderr, "\n"); + } } } @@ -1007,5 +1009,20 @@ psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32); - PAR[PM_PAR_7] = indexGuessInv[nStart]; + float Sm = NAN, Sp = NAN, So = NAN; + if (nStart == 0) { + Sm = indexGuessInv[nStart]; + So = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]); + Sp = indexGuessInv[nStart + 1]; + } else if (nStart == N_INDEX_GUESS_INV - 1) { + Sp = indexGuessInv[nStart]; + So = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]); + Sm = indexGuessInv[nStart - 1]; + } else { + Sm = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]); + So = indexGuessInv[nStart]; + Sp = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]); + } + + PAR[PM_PAR_7] = Sm; pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); @@ -1013,5 +1030,5 @@ psVectorAppend (chi2, model->chisqNorm); - PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : 0.1; + PAR[PM_PAR_7] = So; pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); @@ -1019,23 +1036,9 @@ psVectorAppend (chi2, model->chisqNorm); - PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : 0.55; + PAR[PM_PAR_7] = Sp; pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]); psVectorAppend (Sidx, 100*PAR[PM_PAR_7]); psVectorAppend (chi2, model->chisqNorm); - - if (chi2->data.F32[1] < chi2->data.F32[2]) { - if (nStart == N_INDEX_GUESS_INV - 1) { - PAR[PM_PAR_7] = 0.11; - } else { - PAR[PM_PAR_7] = indexGuessInv[nStart + 1]; - } - } else { - if (nStart == 0) { - PAR[PM_PAR_7] = 0.52; - } else { - PAR[PM_PAR_7] = indexGuessInv[nStart - 1]; - } - } psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2); @@ -1051,7 +1054,11 @@ float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0; - // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev) - Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1); + // constrain Smin to be in a valid range: allow the fitted range to go a bit beyond the 3 trial points, but no further + float Smx = Sm - 0.25*(So - Sm); + float Spx = Sp + 0.25*(Sp - So); + Smin = PS_MAX(PS_MIN(Smin, Smx), Spx); PAR[PM_PAR_7] = Smin; + + // XXX I could set the error on PAR_7 here if I knew how to roughly convert these chisq values to true chisq values // return to the original fitting mode (fitOptions) @@ -1246,3 +1253,162 @@ } - +// # define N_REFF_CHECK 11 +// float drefCheck[] = {-0.02, -0.04, -0.06, 0.0, 0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25}; + +// we have an initial fit, check to see if the current size is besst +bool psphotPCMfitCheckSize (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize) { + + // PAR is already at my current best guess + psF32 *PAR = pcm->modelConv->params->data.F32; + + // store best guess as a shape + psEllipseAxes centerAxes; + pmModelParamsToAxes (¢erAxes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); + + float xMin = NAN; + float iMin = NAN; + float rMin = NAN; + + // loop over Reff, keeping the ARatio and Theta constant + for (int j = -4; j <= 4; j++) { + + float dref = j * 0.01; + + psEllipseAxes guessAxes; + guessAxes.major = centerAxes.major + dref; + guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major; + guessAxes.theta = centerAxes.theta; + + if (!isfinite(guessAxes.major)) return false; + if (!isfinite(guessAxes.minor)) return false; + if (!isfinite(guessAxes.theta)) return false; + + // convert the major,minor,theta to shape parameters for an Reff-like model + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); + + // generated the modelFlux + // XXX note that this does not add sky to model + pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize); + + float YY = 0.0; + float YM = 0.0; + float MM = 0.0; + bool usePoisson = false; + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + // skip masked points + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + continue; + } + // skip zero-variance points + if (source->variance->data.F32[iy][ix] == 0) { + continue; + } + // skip nan value points + if (!isfinite(source->pixels->data.F32[iy][ix])) { + continue; + } + + float fy = source->pixels->data.F32[iy][ix]; + float fm = source->modelFlux->data.F32[iy][ix]; + float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0; + + YY += PS_SQR(fy) * wt; + YM += fm * fy * wt; + MM += PS_SQR(fm) * wt; + } + } + + float Io = YM / MM; + float Chisq = YY - 2 * Io * YM + Io * Io * MM; + if (isnan(xMin) || (Chisq < xMin)) { + xMin = Chisq; + iMin = Io; + rMin = dref; + } + // fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin); + } + + psEllipseAxes guessAxes; + guessAxes.major = centerAxes.major + rMin; + guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major; + guessAxes.theta = centerAxes.theta; + + if (!isfinite(guessAxes.major)) return false; + if (!isfinite(guessAxes.minor)) return false; + if (!isfinite(guessAxes.theta)) return false; + + // convert the major,minor,theta to shape parameters for an Reff-like model + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); + PAR[PM_PAR_I0] = iMin; + + return true; +} + +// we have an initial fit, check to see if the current size is besst +bool psphotPCMfitRetry (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, float psfSize) { + + // PAR is already at my current best guess + psF32 *PAR = pcm->modelConv->params->data.F32; + + // store best guess as a shape + psEllipseAxes centerAxes; + pmModelParamsToAxes (¢erAxes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); + + // retry with axes smaller by 1 pixel + psEllipseAxes guessAxes; + guessAxes.major = centerAxes.major - 0.08; + guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major; + guessAxes.theta = centerAxes.theta; + + if (!isfinite(guessAxes.major)) return false; + if (!isfinite(guessAxes.minor)) return false; + if (!isfinite(guessAxes.theta)) return false; + + // convert the major,minor,theta to shape parameters for an Reff-like model + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true); + + // generated the modelFlux + // XXX note that this does not add sky to model + pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize); + + float YY = 0.0; + float YM = 0.0; + float MM = 0.0; + bool usePoisson = false; + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + // skip masked points + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + continue; + } + // skip zero-variance points + if (source->variance->data.F32[iy][ix] == 0) { + continue; + } + // skip nan value points + if (!isfinite(source->pixels->data.F32[iy][ix])) { + continue; + } + + float fy = source->pixels->data.F32[iy][ix]; + float fm = source->modelFlux->data.F32[iy][ix]; + float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0; + + YY += PS_SQR(fy) * wt; + YM += fm * fy * wt; + MM += PS_SQR(fm) * wt; + } + } + + float Io = YM / MM; + PAR[PM_PAR_I0] = Io; + + pmSourceFitPCM (pcm, source, fitOptions, maskVal, markVal, psfSize); // NOTE : 1687 allocs in here + + return true; +} + + Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 36366) +++ trunk/psphot/src/psphotSourceSize.c (revision 36375) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -# include - + +// this structure is only used internally to simplify the function parameters typedef struct { psImageMaskType maskVal; @@ -13,7 +13,4 @@ bool altDiffExt; float altDiffExtThresh; - bool extFitAll; - bool extFitAllReadout; - float extFitAllThresh; float soft; int grow; @@ -33,5 +30,4 @@ bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options); bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal, int maxWindowCR); -bool psphotMaskCosmicRayFootprintCheck (psArray *sources); int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh); float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction); @@ -41,8 +37,5 @@ // we need to call this function after sources have been fitted to the PSF model and -// subtracted. To determine the CR-nature, this function examines the 9 pixels in the 3x3 -// square containing the peak and compares the observed flux to the model. To determine -// the EXT-nature, this function measures the amount of positive or negative total -// deviation from the psf model at the r = FWHM/2 position +// subtracted. // for now, let's store the detections on the readout->analysis for each readout @@ -121,21 +114,16 @@ assert (status); - // XXX recipe name is not great + // XXX recipe name is not great (NOTE : not used!) options.nSigmaMoments = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.MOMENTS"); assert (status); - // Optional extended source measurement algorithm to improve diff image trails + // Optional algorithm to define if a source is extended (used by DIFF analysis) options.altDiffExt = psMetadataLookupBool(&status, recipe, "PSPHOT.EXT.DIFF.ALTERNATE"); assert (status); + // Threshold for this alternate method options.altDiffExtThresh = psMetadataLookupF32(&status, recipe, "PSPHOT.EXT.DIFF.ALTERNATE.THRESH"); assert (status); - // Option to enable fitting of all objects with extended model. - options.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"); - assert (status); - + // location of a single test source options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST"); @@ -422,19 +410,5 @@ psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nskip"); - // Determine if this readout is above the threshold to ext fit all sources - options->extFitAllReadout = false; - if (options->extFitAll) { - float maskFrac = psMetadataLookupF32(&status,readout->analysis,"READOUT.MASK.FRAC"); - if (status) { - maskFrac = 0.0; - } - if (sources->n * (1.0 - maskFrac) > options->extFitAllThresh) { - options->extFitAllReadout = false; - } - else { - options->extFitAllReadout = true; - } - } - + if (!psphotSourceClassRegion (NULL, &psfClump, sources, recipe, psf, options)) { psLogMsg ("psphot", 4, "Failed to determine source classification for full image\n"); @@ -444,4 +418,5 @@ return true; + // NOTE : this section is deactivated (EAM : I think we were getting poor boundary effects?) int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS"); for (int i = 0; i < nRegions; i ++) { @@ -582,10 +557,4 @@ // * SAT stars should not be faint, but defects may? - // 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; - } // Defects may not always match CRs from peak curvature analysis // Defects may also be marked as SATSTAR -- XXX deactivate this flag? @@ -595,5 +564,6 @@ // saturated star (too many saturated pixels or peak above saturation limit). These - // may also be saturated galaxies, or just large saturated regions. + // may also be saturated galaxies, or just large saturated regions. They are never + // marked as 'extended' if (source->mode & PM_SOURCE_MODE_SATSTAR) { psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g SAT\n", @@ -604,5 +574,6 @@ } - // any sources missing a large fraction should just be treated as PSFs + // any sources missing a large fraction should just be treated as PSFs. They are never + // marked as 'extended' if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) { psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g PSF\t%g %g\n", @@ -634,4 +605,5 @@ psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g EXT\n", source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid); + source->type = PM_SOURCE_TYPE_EXTENDED; source->mode |= PM_SOURCE_MODE_EXT_LIMIT; source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; @@ -639,4 +611,5 @@ continue; } + // Alternate extended source limit calculation if (options->altDiffExt) { @@ -652,4 +625,5 @@ source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid, momentRatioVeres,options->altDiffExtThresh); + source->type = PM_SOURCE_TYPE_EXTENDED; source->mode |= PM_SOURCE_MODE_EXT_LIMIT; source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; @@ -968,69 +942,4 @@ } -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; -} - -/**** ------ old versions of cosmic ray masking ----- ****/ - -bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); - -// 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; -} - # define VERBOSE 0 int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh) { @@ -1176,278 +1085,4 @@ } -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; -} - -// 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; -} - float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction) { Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 36366) +++ trunk/psphot/src/psphotSourceStats.c (revision 36375) @@ -58,12 +58,4 @@ pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); psAssert (detections, "missing detections?"); - - // XXX TEST: - if (detections->allSources) { - psphotMaskCosmicRayFootprintCheck(detections->allSources); - } - if (detections->newSources) { - psphotMaskCosmicRayFootprintCheck(detections->newSources); - } // determine the number of allowed threads @@ -246,14 +238,4 @@ psphotVisualShowMoments (sources); - // clear the mark bits - // psImageMaskPixels (readout->mask, "AND", PS_NOT_IMAGE_MASK(markVal)); - - if (detections->allSources) { - psphotMaskCosmicRayFootprintCheck(detections->allSources); - } - if (detections->newSources) { - psphotMaskCosmicRayFootprintCheck(detections->newSources); - } - return true; } @@ -572,4 +554,6 @@ // measure basic source moments (no S/N clipping on input pixels) + // sources with (mode & MODE_EXTERNAL) or (mode2 & MODE2_MATCHED) use the + // supplied Mx,My value for the centroid (not recalculated) status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, 0.0, maskVal); } Index: trunk/psphot/src/psphotStackReadout.c =================================================================== --- trunk/psphot/src/psphotStackReadout.c (revision 36366) +++ trunk/psphot/src/psphotStackReadout.c (revision 36375) @@ -124,19 +124,4 @@ return psphotReadoutCleanup (config, view, STACK_RAW); } - -// XXX TEST for background: - if (!psphotModelBackground (config, view, STACK_RAW)) { - return psphotReadoutCleanup (config, view, STACK_RAW); - } - if (!psphotSubtractBackground (config, view, STACK_RAW)) { - return psphotReadoutCleanup (config, view, STACK_RAW); - } - if (!psphotModelBackground (config, view, STACK_RAW)) { - return psphotReadoutCleanup (config, view, STACK_RAW); - } - if (!psphotSubtractBackground (config, view, STACK_RAW)) { - return psphotReadoutCleanup (config, view, STACK_RAW); - } -// XXX TEST END #ifdef MAKE_CHISQ_IMAGE @@ -319,5 +304,5 @@ psphotStackObjectsUnifyPosition (objects); - psphotStackObjectsSelectForAnalysis (config, view, STACK_RAW, objects); + // psphotStackObjectsSelectForAnalysis (config, view, STACK_RAW, objects); // final linear fit. NOTE: if splitLinearFit is true above, this pass will only fit @@ -346,4 +331,7 @@ // measure kron fluxes for the matched sources only psphotKronIterate(config, view, STACK_RAW, 3); + + // decide which source(s) are to be fitted with the extended source analysis code. + psphotChooseAnalysisOptionsByObject (config, view, STACK_RAW, objects); // measure elliptical apertures, petrosians (objects sorted by S/N) Index: trunk/psphot/test/tap_psphot_galaxygrid.pro =================================================================== --- trunk/psphot/test/tap_psphot_galaxygrid.pro (revision 36366) +++ trunk/psphot/test/tap_psphot_galaxygrid.pro (revision 36375) @@ -30,8 +30,14 @@ $RefOptions = $RefOptions -nx 3000 -ny 3000 +if (not($?PSFMODEL)) + $PSFMODEL = PS1_V1 +end + +macro reset.options # options for the simulated images (using the refimage for the stars) $FakeOptions = $BaseOptions $FakeOptions = $FakeOptions -exptime 30.0 -$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS +# $FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS +$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_$PSFMODEL $FakeOptions = $FakeOptions -nx 3000 -ny 3000 @@ -69,4 +75,7 @@ $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300 $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300 +end + +if (not($?FakeConfig)) reset.options list fwhm @@ -75,4 +84,48 @@ 1.2 1.5 +end + +if (not($?CONVOLVE_NSIGMA)) set CONVOLVE_NSIGMA = 5.0 + +# generate fake images and run psphot on them +macro mkexp.devexp.single + if ($0 != 5) + echo "USAGE: mkexp.devexp.single (basename) (type) (Rmajor) (fwhm)" + break + end + + $basename = $1 + $type = $2 + $Rmajor = $3 + $fwhm = $4 + + $Aratio = 1.0 + + $FakeConfig = -camera SIMTEST + $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN + $FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref + $FakeConfig = $FakeConfig -Db STARS.FAKE F ; # only use stars from catdir.ref + $FakeConfig = $FakeConfig -Db MATCH.DENSITY F + $FakeConfig = $FakeConfig -Db PSF.CONVOLVE T + $FakeConfig = $FakeConfig -Db GALAXY.FAKE T ; # generate a "realistic" distribution of galaxies + $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0 + $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5 + $FakeConfig = $FakeConfig -Db GALAXY.GRID T ; # generate a grid of galaxies (constant mag) + $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0 + $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180 + $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300 + $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300 + $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0 + $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0 + $FakeConfig = $FakeConfig -Df CONVOLVE.NSIGMA $CONVOLVE_NSIGMA + $BaseConfig = $FakeConfig + + $FakeConfig = $BaseConfig + $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor + $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor + $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN $Aratio + $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio + + mkexp $basename $fwhm $type end @@ -272,4 +325,54 @@ label -x sequence -y "M_out| - M_in|" resize 700 320 + + # check on magnitude + set dI = Iot_s - Iin_s + lim -n 4$1 n -1.0 1.0; clear; box; plot n dI + label -x sequence -y "I_out| - I_in|" + resize 700 320 +end + +macro grid.plot.stars + if ($0 != 2) + echo "USAGE: grid.plot.stars (version)" + break + end + # things to examine: theta, Rmajor, AR + + # go.grid.check.devexp + + # check on position + set dX = Xot_s - int(Xin_s) - 0.5 + set dY = Yot_s - int(Yin_s) - 0.5 + set ARin = rin_s / Rin_s + set dR = Rin_s - Rot_s + + create n 0 dR[] + set dRf = dR / Rin_s + lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf + label -x sequence -y "1 - R_out| / R_in|" + resize 700 320 + + # check on A.Ratio + set ARot = rot_s / Rot_s + # lim ARin ARot; clear; box; plot ARin ARot + set fAR = ARot / ARin + lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR + label -x sequence -y "AR_out| / AR_in|" + resize 700 320 + + # check on magnitude + set dM = Mot_s - Min_s + lim -n 3$1 n -0.5 0.5; clear; box; plot n dM + label -x sequence -y "M_out| - M_in|" + resize 700 320 + + lim -n 4$1 n -1.5 1.5; clear; box; plot n dX + label -x sequence -y "X_out| - X_in|" + resize 700 320 + + lim -n 5$1 n -1.5 1.5; clear; box; plot n dY + label -x sequence -y "Y_out| - Y_in|" + resize 700 320 end @@ -487,4 +590,6 @@ end +if (not($?NSIGMA_CONV)) set NSIGMA_CONV = 5.0 + # create a realistic distribution of fake stars, GAUSS PSF macro fitexp @@ -507,4 +612,12 @@ $psphotConfig = $psphotConfig -Db PSPHOT:SAVE.RESID T $psphotConfig = $psphotConfig -D PSPHOT:EXTENDED_SOURCE_MODELS_SELECTION $fitModel + $psphotConfig = $psphotConfig -D PSPHOT:PSF_MODEL PS_MODEL_$PSFMODEL + + $psphotConfig = $psphotConfig -Db PSPHOT:PSF.RESIDUALS F + $psphotConfig = $psphotConfig -Db PSPHOT:POISSON.ERRORS.PHOT.LMM F + $psphotConfig = $psphotConfig -Db PSPHOT:EXTENDED_SOURCE_FITS_POISSON F + $psphotConfig = $psphotConfig -Di PSPHOT:EXT_FIT_ITER 15 + $psphotConfig = $psphotConfig -Df PSPHOT:PSF_FIT_RADIUS_SCALE 3.75 + $psphotConfig = $psphotConfig -Df PSPHOT:EXT_FIT_NSIGMA_CONV $NSIGMA_CONV # ppImage / psphot on the output @@ -516,7 +629,18 @@ end +macro cmf.load.reset + $fields = X Y M T R r MT I + foreach field $fields + foreach set in ot + delete -q $field\$set\_s + end + end + + delete -q min_S Min_S +end + macro cmf.load.concat if ($0 != 4) - echo "USAGE: cmf.load.concat (dat) (cmf) (type)" + echo "USAGE: cmf.load.concat (dat) (cmf) (inType)" break end @@ -525,4 +649,20 @@ read Xin_all 1 Yin_all 2 Fin_all 3 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10 + $TYPE_S = 83 + $TYPE_D = 68 + $TYPE_E = 69 + if ("$3" == "SERSIC") + $InType = $TYPE_S + end + if ("$3" == "DEV") + $InType = $TYPE_D + set IndexIn_all = 0.125 + zero(Xin_all) + end + if ("$3" == "EXP") + $InType = $TYPE_E + set IndexIn_all = 0.5 + zero(Xin_all) + end + + # select only the galaxies subset Xin = Xin_all if (Type == 1) subset Yin = Yin_all if (Type == 1) @@ -539,26 +679,15 @@ subset IndexIn = IndexIn_all if (Type == 1) - $TYPE_S = 83 - $TYPE_D = 68 - $TYPE_E = 69 - if ("$3" == "SERSIC") - $InType = $TYPE_S - end - if ("$3" == "DEV") - $InType = $TYPE_D - end - if ("$3" == "EXP") - $InType = $TYPE_E - end - data $2 break -auto off + output -err /dev/null read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE EXT_PAR_07 $reread = not($STATUS) + output -err stderr break -auto on if ($reread) read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE - set EXT_PAR_07 = (MODEL_TYPE:9 == 68)*4 + (MODEL_TYPE:9 == 69) + set EXT_PAR_07 = (MODEL_TYPE:9 == 68)*0.125 + (MODEL_TYPE:9 == 69)*0.5 end @@ -589,11 +718,7 @@ reindex rin_m = RminIn using index2 - if ("$3" == "SERSIC") - reindex Iot_m = EXT_PAR_07 using index1 - reindex Iin_m = IndexIn using index2 - $fields = X Y M T R r MT I - else - $fields = X Y M T R r MT - end + reindex Iot_m = EXT_PAR_07 using index1 + reindex Iin_m = IndexIn using index2 + $fields = X Y M T R r MT I foreach field $fields @@ -605,4 +730,57 @@ concat min min_S concat Min Min_S +end + +macro cmf.load.stars.concat + if ($0 != 3) + echo "USAGE: cmf.load.concat (dat) (cmf)" + break + end + + data $1 + read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10 + + subset Xin = Xin_all if (Type == 0) + subset Yin = Yin_all if (Type == 0) + subset Min = Min_all if (Type == 0) + subset RmajIn = RmajIn_all if (Type == 0) + subset RminIn = RminIn_all if (Type == 0) + subset IndexIn = IndexIn_all if (Type == 0) + subset Tin_rad = ThetaIn_all if (Type == 0) + set Tin = Tin_rad * 180 / 3.14159265 + + data $2 + + break -auto off + read -fits Chip.psf X_PSF Y_PSF PSF_INST_MAG PSF_MAJOR PSF_MINOR PSF_THETA + set PSF_THETA_ALT = PSF_THETA * (PSF_THETA >= 0.0) + (PSF_THETA + 3.14159265) * (PSF_THETA < 0.0) + set PSF_THETA = PSF_THETA_ALT * 180 / 3.14159265 + + match2d X_PSF Y_PSF Xin Yin 1.0 -index1 index1 -index2 index2 + + reindex Xot_m = X_PSF using index1 + reindex Yot_m = Y_PSF using index1 + + reindex Xin_m = Xin using index2 + reindex Yin_m = Yin using index2 + + reindex Mot_m = PSF_INST_MAG using index1 + reindex Tot_m = PSF_THETA using index1 + + reindex Min_m = Min using index2 + reindex Tin_m = Tin using index2 + + reindex Rot_m = PSF_MAJOR using index1 + reindex rot_m = PSF_MINOR using index1 + + reindex Rin_m = RmajIn using index2 + reindex rin_m = RminIn using index2 + + $fields = X Y M T R r + foreach field $fields + foreach set in ot + concat $field\$set\_m $field\$set\_s + end + end end @@ -1078,6 +1256,6 @@ macro load.normdata - if ($0 != 2) - echo "USAGE: load.normdata (file)" + if ($0 != 3) + echo "USAGE: load.normdata (file) [clear/noclear]" break end @@ -1090,5 +1268,8 @@ set dm = mg - Mg set n = ramp(dm) - lim n dm; clear; box; plot n dm + if ("$2" == "clear") + lim n dm; clear; box; + end + plot n dm end @@ -1097,2 +1278,168 @@ exit 0 end + +# note that t1 = original confi +# t2 = no residuals +# t3 = const weight PSF fit +# t4 = const weight galaxy fits + +macro run.radius.loop + local radius Nrun + + if (1) + $Nrun = 0 + foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0 + mkexp.devexp.single tests.20131120/testrad.ps1v1.exp.$Nrun EXP $radius 1.0 + fitexp tests.20131120/testrad.ps1v1.exp.$Nrun tests.20131120/testrad.ps1v1.exp.$Nrun.t4b EXP_CONV + $Nrun ++ + end + end + + $Nrun = 0 + foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0 + cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4b.cmf EXP + echo -no-return $radius "" + check.fit tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4b.cmf EXP + $Nrun ++ + end + + grid.plots.devexp a +end + +macro check.radius.loop + if ($0 != 2) + echo "USAGE: check.radius.loop (ext)" + break + end + + # cmf.load.reset + + local radius Nrun + + $Nrun = 0 + foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0 + echo -no-return $radius "" + check.fit tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.$1.cmf EXP + $Nrun ++ + end + + grid.plots.devexp a +end + +macro check.fit + if ($0 != 4) + echo "USAGE: check.fit (dat) (cmf) (type)" + break + end + + cmf.load.concat $1 $2 $3 + set dM_m = Mot_m - Min_m + vstat -q dM_m + $mag_off = $MEDIAN + set dR_m = Rot_m - Rin_m + vstat -q dR_m + $rad_off = $MEDIAN + echo $mag_off $rad_off +end + +macro make.circles + + delete Xo Yo Ro + + for ix 0 16 + for iy $ix 16 + $r2 = $ix^2 + $iy^2 + if ($r2 > 15^2) continue + concat $ix Xo + concat $iy Yo + concat $r2 Ro + end + end + + sort Ro Xo Yo + $Rold = -1 + $N = -1 + for i 0 Ro[] + if (Ro[$i] != $Rold) + $N ++ + $Rold = Ro[$i] + end + if ((Xo[$i] == 0) && (Yo[$i] == 0)) + fprintf "// center is 0,0" + continue + end + if (Xo[$i] == 0) + fprintf "ADD_AXIS (%3d, %2d) // r^2 = %3d" $N Yo[$i] Ro[$i] + continue + end + if (Xo[$i] == Yo[$i]) + fprintf "ADD_DIAG (%3d, %2d) // r^2 = %3d" $N Xo[$i] Ro[$i] + continue + end + + fprintf "ADD_RAND (%3d, %2d, %2d) // r^2 = %3d" $N Xo[$i] Yo[$i] Ro[$i] + end +end + +# EXP : tests.20131120/test.nsig +# DEV : tests.20131120/test.dev + +# run.convolve.loop tests.20131120/test.nsig tests.20131120/test.nsig EXP EXP_CONV +# run.convolve.loop tests.20131120/test.dev tests.20131120/test.dev DEV DEV_CONV +# run.convolve.loop tests.20131120/test.nsig tests.20131120/test.serexp EXP SER_CONV +# run.convolve.loop tests.20131120/test.dev tests.20131120/test.serdev DEV SER_CONV + +macro run.convolve.loop + if ($0 != 5) + echo "run.convolve.loop (inName) (fitName) (inType) (fitType)" + break + end + + local radius + $radius = 3.0 + cmf.load.reset + + # foreach Cin 3 5 7 9 11 + foreach Cin 11 + $CONVOLVE_NSIGMA = $Cin + sprintf nameIn %s.%02d $1 $Cin + mkexp.devexp.single $nameIn $3 $radius 1.0 + foreach Cot 3 5 7 9 11 + $NSIGMA_CONV = $Cot + sprintf nameOt %s.%02d.%02d $2 $Cin $Cot + fitexp $nameIn $nameOt $4 + end + end + + # foreach Cin 3 5 7 9 11 + foreach Cin 11 + sprintf nameIn %s.%02d $1 $Cin + foreach Cot 3 5 7 9 11 + sprintf nameOt %s.%02d.%02d $2 $Cin $Cot + check.fit $nameIn.dat $nameOt.cmf + end + end + grid.plots.devexp a +end + +macro check.convolve.loop + if ($0 != 4) + echo "run.convolve.loop (inName) (fitName) (inType)" + break + end + + local radius + $radius = 3.0 + cmf.load.reset + + # foreach Cin 3 5 7 9 11 + foreach Cin 11 + sprintf nameIn %s.%02d $1 $Cin + foreach Cot 3 5 7 9 11 + sprintf nameOt %s.%02d.%02d $2 $Cin $Cot + echo -no-return $Cin $Cot " " + check.fit $nameIn.dat $nameOt.cmf $3 + end + end + grid.plots.devexp a +end