Index: /branches/eam_branches/psphot.stack.20100120/src/psphot.h =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphot.h (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphot.h (revision 26681) @@ -75,4 +75,6 @@ bool psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal); +bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources); + bool psphotReplaceAllSources (psArray *sources, psMetadata *recipe); bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe); @@ -126,4 +128,6 @@ bool psphotSetMaskAndVarianceReadout (pmConfig *config, pmReadout *readout, psMetadata *recipe); void psphotSourceFreePixels (psArray *sources); + +bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *target, psMetadata *recipe, const bool havePSF); // functions to set the correct source pixels Index: /branches/eam_branches/psphot.stack.20100120/src/psphotAddNoise.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotAddNoise.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotAddNoise.c (revision 26681) @@ -1,13 +1,35 @@ # include "psphotInternal.h" -bool psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) { - return psphotAddOrSubNoise (readout, sources, recipe, true); +bool psphotAddNoise (pmConfig *config, const pmFPAview *view) { + return psphotAddOrSubNoise (config, view, true); } -bool psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) { - return psphotAddOrSubNoise (readout, sources, recipe, false); +bool psphotSubNoise (pmConfig *config, const pmFPAview *view) { + return psphotAddOrSubNoise (config, view, false); } -bool psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotAddOrSubNoiseReadout (config, view, "PSPHOT.INPUT", i, recipe, add)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add) { bool status = false; @@ -16,7 +38,15 @@ psEllipseAxes axes; - PS_ASSERT (readout, false); - PS_ASSERT (readout->parent, false); - PS_ASSERT (readout->parent->concepts, false); + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (readout, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + psAssert (readout->parent, "missing cell?"); + psAssert (readout->parent->concepts, "missing concepts?"); + + psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); + psAssert (sources, "missing sources?"); psTimerStart ("psphot.noise"); @@ -24,5 +54,5 @@ // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels - assert (maskVal); + psAssert (maskVal, "missing mask value?"); // increase variance by factor*(object noise): Index: /branches/eam_branches/psphot.stack.20100120/src/psphotBasicDeblend.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotBasicDeblend.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotBasicDeblend.c (revision 26681) @@ -1,3 +1,21 @@ # include "psphotInternal.h" + +// for now, let's store the detections on the readout->analysis for each readout +bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotBasicDeblendReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view, const char *filename, int index) { @@ -21,4 +39,9 @@ psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping basic deblend"); + return true; + } // select the appropriate recipe information @@ -142,20 +165,2 @@ return true; } - -// for now, let's store the detections on the readout->analysis for each readout -bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view) -{ - bool status = true; - - int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); - psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - if (!psphotBasicDeblendReadout (config, view, "PSPHOT.INPUT", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); - return false; - } - } - return true; -} Index: /branches/eam_branches/psphot.stack.20100120/src/psphotBlendFit.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotBlendFit.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotBlendFit.c (revision 26681) @@ -1,6 +1,28 @@ # include "psphotInternal.h" +// for now, let's store the detections on the readout->analysis for each readout +bool psphotBlendFit (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotBlendFitReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + // XXX I don't like this name -bool psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) { +bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { int Nfit = 0; @@ -12,7 +34,21 @@ psTimerStart ("psphot.fit.nonlinear"); - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - assert (recipe); + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (readout, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); // determine the number of allowed threads @@ -47,4 +83,8 @@ // source analysis is done in S/N order (brightest first) sources = psArraySort (sources, pmSourceSortBySN); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + return true; + } // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) Index: /branches/eam_branches/psphot.stack.20100120/src/psphotChoosePSF.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotChoosePSF.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotChoosePSF.c (revision 26681) @@ -1,8 +1,36 @@ # include "psphotInternal.h" +// generate a PSF model for inputs without PSF models already loaded +bool psphotChoosePSF (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + // try PSF models and select best option -bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { +bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; + + // do not generate a PSF if we already were supplied one + if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { + psLogMsg ("psphot", PS_LOG_DETAIL, "psf model supplied for input file %d", i); + return true; + } psTimerStart ("psphot.choose.psf"); @@ -18,7 +46,8 @@ psAssert (sources, "missing sources?"); - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - psAssert (recipe, "missing recipe?"); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping PSF model"); + return true; + } // bit-masks to test for good/bad pixels @@ -322,4 +351,8 @@ return false; } + psFree (psf); // XXX double-check + + // move into psphotChoosePSF + psphotVisualShowPSFModel (readout, psf); return true; @@ -466,20 +499,2 @@ return true; } - -// for now, let's store the detections on the readout->analysis for each readout -bool psphotChoosePSF (pmConfig *config, const pmFPAview *view) -{ - bool status = true; - - int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); - psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, havePSF)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); - return false; - } - } - return true; -} Index: /branches/eam_branches/psphot.stack.20100120/src/psphotDeblendSatstars.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotDeblendSatstars.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotDeblendSatstars.c (revision 26681) @@ -1,3 +1,21 @@ # include "psphotInternal.h" + +// for now, let's store the detections on the readout->analysis for each readout +bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotDeblendSatstarsReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { @@ -20,4 +38,9 @@ psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend"); + return true; + } // select the appropriate recipe information @@ -183,20 +206,2 @@ return true; } - -// for now, let's store the detections on the readout->analysis for each readout -bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view) -{ - bool status = true; - - int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); - psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - if (!psphotDeblendSatstarsReadout (config, view, "PSPHOT.INPUT", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); - return false; - } - } - return true; -} Index: /branches/eam_branches/psphot.stack.20100120/src/psphotFitSourcesLinear.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotFitSourcesLinear.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotFitSourcesLinear.c (revision 26681) @@ -12,5 +12,23 @@ static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal); -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final) +{ + bool status = true; + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotFitSourcesLinearReadout (config, view, "PSPHOT.INPUT", i, final)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { bool status; @@ -19,4 +37,26 @@ float f; // float r; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (readout, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (sources, "missing psf?"); psTimerStart ("psphot.linear"); @@ -248,4 +288,8 @@ psphotVisualPlotChisq (sources); // psphotVisualShowFlags (sources); + + // We have to place this visualization here because the models are not realized until + // psphotGuessModels or fitted until psphotFitSourcesLinear. + psphotVisualShowPSFStars (recipe, psf, sources); return true; Index: /branches/eam_branches/psphot.stack.20100120/src/psphotGuessModels.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotGuessModels.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotGuessModels.c (revision 26681) @@ -7,10 +7,46 @@ // 3) define the threaded function to work with sources for a given cell +// for now, let's store the detections on the readout->analysis for each readout +bool psphotGuessModels (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotGuessModelsReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + // construct an initial PSF model for each object -bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) { +bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { bool status; psTimerStart ("psphot.models"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (readout, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping model guess"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (sources, "missing psf?"); // select the appropriate recipe information Index: /branches/eam_branches/psphot.stack.20100120/src/psphotImageQuality.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotImageQuality.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotImageQuality.c (revision 26681) @@ -4,7 +4,49 @@ // XXX Lots of code duplication here +// for now, let's store the detections on the readout->analysis for each readout +bool psphotImageQuality (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotImageQualityReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + // selecting the 'good' stars (likely to be psf stars), measure the M_cn, M_sn terms for n = 2,3,4 -bool psphotImageQuality(psMetadata *recipe, psArray *sources) +bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (readout, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // if we have a PSF, skip this analysis + if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { + return true; + } + + psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping image quality"); + return true; + } + psVector *FWHM_MAJOR = psVectorAllocEmpty(100, PS_TYPE_F32); psVector *FWHM_MINOR = psVectorAllocEmpty(100, PS_TYPE_F32); Index: /branches/eam_branches/psphot.stack.20100120/src/psphotLoadPSF.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotLoadPSF.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotLoadPSF.c (revision 26681) @@ -2,13 +2,17 @@ // load an externally supplied psf model -pmPSF *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe) { +pmPSF *psphotLoadPSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (readout, "missing file?"); // find the currently selected chip - pmChip *chip = pmFPAfileThisChip (config->files, view, "PSPHOT.PSF.LOAD"); - if (!chip) return NULL; + pmChip *chip = pmFPAviewThisChip (view, file->fpa); + if (!chip) return false; // find the currently selected readout - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.PSF.LOAD"); - if (!readout) return NULL; + pmReadout *readout = pmFPAfileThisReadout (view, file->fpa); + if (!readout) return false; // check if a PSF model is supplied by the user @@ -16,13 +20,43 @@ if (psf == NULL) { psLogMsg ("psphot", 3, "no psf supplied for this chip"); - return NULL; + return true; } - if (!psphotPSFstats (readout, recipe, psf)) psAbort("cannot measure PSF shape terms"); + if (!psphotPSFstats (readout, recipe, psf)) { + psAbort("cannot measure PSF shape terms"); + } psLogMsg ("psphot", 3, "using externally supplied PSF model for this readout"); - // we return a psf which can be free (and which may be removed from the metadata) - psMemIncrRefCounter (psf); - return psf; + // save PSF on readout->analysis + if (!psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_PTR, "psphot psf model", psf)) { + psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout"); + return false; + } + + return true; } + +bool psphotLoadPSF (pmConfig *config, const pmFPAview *view) { + + bool status = false; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // XXX PSPHOT.PSF.LOAD vs PSPHOT.INPUT?? + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // 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, "PSPHOT.PSF.LOAD", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to load PSF model for PSPHOT.PSF.LOAD entry %d", i); + return false; + } + } + return true; +} Index: /branches/eam_branches/psphot.stack.20100120/src/psphotMaskReadout.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotMaskReadout.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotMaskReadout.c (revision 26681) @@ -1,8 +1,41 @@ # include "psphotInternal.h" +bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view) { + + bool status = false; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // 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 (!psphotSetMaskAndVarianceReadout (config, readout, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for PSPHOT.INPUT entry %d", i); + return false; + } + + // display the image, weight, mask (ch 1,2,3) + psphotVisualShowImage (readout); + } + return true; +} + // generate mask and variance if not defined, additional mask for restricted subregion -bool psphotSetMaskAndVarianceReadout (pmConfig *config, pmReadout *readout, psMetadata *recipe) { +bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; + + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest + psAssert (file, "missing file?"); + + // find the currently selected readout + pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); + psAssert (readout, "missing readout?"); // ** Interpret the mask values: @@ -73,33 +106,7 @@ } + // display the image, weight, mask (ch 1,2,3) + psphotVisualShowImage (readout); + return true; } - -bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view, psMetadata *recipe) { - - bool status = false; - - int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); - psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest - PS_ASSERT_PTR_NON_NULL (file, false); - - // find the currently selected readout - pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); - PS_ASSERT_PTR_NON_NULL (readout, false); - - // Generate the mask and weight images, including the user-defined analysis region of interest - if (!psphotSetMaskAndVarianceReadout (config, readout, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for PSPHOT.INPUT entry %d", i); - return false; - } - - // display the image, weight, mask (ch 1,2,3) - psphotVisualShowImage (readout); - } - return true; -} Index: /branches/eam_branches/psphot.stack.20100120/src/psphotModelBackground.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotModelBackground.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotModelBackground.c (revision 26681) @@ -356,8 +356,9 @@ // find the currently selected readout pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest - PS_ASSERT_PTR_NON_NULL (file, false); + psAssert (file, "missing file?"); pmFPA *inFPA = file->fpa; pmReadout *readout = pmFPAviewThisReadout(view, inFPA); + psAssert (readout, "missing readout?"); psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters Index: /branches/eam_branches/psphot.stack.20100120/src/psphotOutput.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotOutput.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotOutput.c (revision 26681) @@ -161,5 +161,5 @@ } -bool psphotSetHeaderNstars (psMetadata *recipe, psArray *sources) { +bool psphotSetHeaderNstars (psMetadata *header, psArray *sources) { int nSrc = 0; @@ -190,14 +190,14 @@ } - psMetadataAddS32 (recipe, PS_LIST_TAIL, "NSTARS", PS_META_REPLACE, "Number of sources", nSrc); - psMetadataAddS32 (recipe, PS_LIST_TAIL, "NDET_EXT", PS_META_REPLACE, "Number of extended sources", nEXT); - psMetadataAddS32 (recipe, PS_LIST_TAIL, "NDET_CR", PS_META_REPLACE, "Number of cosmic rays", nCR); - return true; -} - -// these values are saved in an output header stub - they are added to either the -// PHU header (CMP) or the MEF table header (CMF) -// XXX these are currently carried by the recipe -- this will not work in a Stack context -// XXX they should be place in the readout->analysis of the given image + psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", PS_META_REPLACE, "Number of sources", nSrc); + psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_EXT", PS_META_REPLACE, "Number of extended sources", nEXT); + psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_CR", PS_META_REPLACE, "Number of cosmic rays", nCR); + return true; +} + +// these values are saved in an output header stub - they are added to either the PHU header +// (CMP) or the MEF table header (CMF) XXX these are currently carried by the recipe -- this +// will not work in a Stack context XXX they should be place in the readout->analysis of the +// given image or directly on the output header psMetadata *psphotDefineHeader (psMetadata *recipe) { Index: /branches/eam_branches/psphot.stack.20100120/src/psphotReadout.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotReadout.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotReadout.c (revision 26681) @@ -3,5 +3,4 @@ // this should be called by every program that links against libpsphot bool psphotInit (void) { - psphotErrorRegister(); // register our error codes/messages psphotModelClassInit (); // load implementation-specific models @@ -12,8 +11,6 @@ bool psphotReadout(pmConfig *config, const pmFPAview *view) { - // measure the total elapsed time in psphotReadout. XXX the current threading plan - // for psphot envisions threading within psphotReadout, not multiple threads calling - // the same psphotReadout. In the current plan, this dtime is the elapsed time used - // jointly by the multiple threads, not the total time used by all threads. + // measure the total elapsed time in psphotReadout. dtime is the elapsed time used jointly + // by the multiple threads, not the total time used by all threads. psTimerStart ("psphotReadout"); @@ -26,4 +23,7 @@ return false; } + // optional break-point for processing + char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); + psAssert (breakPt, "configuration error: set BREAK_POINT"); // set the photcode for this image @@ -33,122 +33,102 @@ } - // find the currently selected readout - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - PS_ASSERT_PTR_NON_NULL (readout, 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, recipe); + psphotSetMaskAndVariance (config, view); if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL); - } - - // display the image, weight, mask (ch 1,2,3) - // XXX move this into the loop above - psphotVisualShowImage (readout); + return psphotReadoutCleanup(config, view); + } // generate a background model (median, smoothed image) if (!psphotModelBackground (config, view)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } if (!psphotSubtractBackground (config, view)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); - } - - // display the backsub and backgnd images - // move this inthe the subtract background loop - psphotVisualShowBackground (config, view, readout); + return psphotReadoutCleanup (config, view); + } // run a single-model test if desired (exits from here if test is run) - psphotModelTest (config, view, recipe); + // XXX drop this and only keep the stand-alone program? + // psphotModelTest (config, view, recipe); // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe - pmPSF *psf = psphotLoadPSF (config, view, recipe); - - // several functions below behave differently if we have a PSF model already - bool havePSF = (psf != NULL); - + // XXX this needs to save the PSF on the chip->analysis or readout->analysis + if (!psphotLoadPSF (config, view, recipe)) { + // this only happens if we had a programming error in psphotLoadPSF + psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); + return psphotReadoutCleanup (config, view); + } + // find the detections (by peak and/or footprint) in the image. if (!psphotFindDetections (config, view)) { // this only happens if we had an error in psphotFindDetections psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL); - } + return psphotReadoutCleanup (config, view); + } + +XXX: // need to handle or ignore this !! if (!detections->peaks->n) { psLogMsg ("psphot", 3, "unable to find detections in this image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL); + return psphotReadoutCleanup (config, view); } // construct sources and measure basic stats if (!psphotSourceStats (config, view, true)) { - // XXX do I need to raise an error here? - return false; + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); + return psphotReadoutCleanup (config, view); } if (!strcasecmp (breakPt, "PEAKS")) { - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup(config, view); } // find blended neighbors of very saturated stars - // XXX merge this with Basic Deblend? if (!psphotDeblendSatstars (config, view)) { - psLogMsg ("psphot", 3, "failed on satstar deblend analysis"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); + return psphotReadoutCleanup (config, view); } // mark blended peaks PS_SOURCE_BLEND if (!psphotBasicDeblend (config, view)) { - psLogMsg ("psphot", 3, "failed on deblend analysis"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - - // classify sources based on moments, brightness - if (!psphotRoughClass (config, view, havePSF)) { - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis"); + return psphotReadoutCleanup (config, view); + } + + // classify sources based on moments, brightness. if a PSF model has been loaded, the PSF + // clump defined for it is used not measured + if (!psphotRoughClass (config, view)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications"); + return psphotReadoutCleanup (config, view); } if (!strcasecmp (breakPt, "MOMENTS")) { - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); - } - - if (!havePSF && !psphotImageQuality (recipe, sources)) { + return psphotReadoutCleanup(config, view); + } + + // if we were not supplied a PSF model, determine the IQ stats here + if (!psphotImageQuality (config, view)) { psLogMsg("psphot", 3, "failed to measure image quality"); - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); - } - - // if we were not supplied a PSF, choose one here - if (psf == NULL) { - // use bright stellar objects to measure PSF - // XXX if we do not have enough stars to generate the PSF, build one - // from the SEEING guess and model class - if (!psphotChoosePSF (config, view)) { - psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - havePSF = true; + return psphotReadoutCleanup(config, view); + } + + // use bright stellar objects to measure PSF if we were supplied a PSF for any input file, + // this step is skipped + if (!psphotChoosePSF (config, view)) { + psLogMsg ("psphot", 3, "failure to construct a psf model"); + return psphotReadoutCleanup (config, view); } if (!strcasecmp (breakPt, "PSFMODEL")) { - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - // move into psphotChoosePSF - psphotVisualShowPSFModel (readout, psf); + return psphotReadoutCleanup (config, view); + } // include externally-supplied sources + // XXX hmm... psphotLoadExtSources (config, view, sources); // construct an initial model for each object, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, readout, sources, psf); + psphotGuessModels (config, view); // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); - - // We have to place these visualizations here because the models are not realized until - // psphotGuessModels or fitted until psphotFitSourcesLinear. - psphotVisualShowPSFStars (recipe, psf, sources); + psphotFitSourcesLinear (config, view, FALSE); // identify CRs and extended sources @@ -157,15 +137,14 @@ goto finish; } - psphotVisualShowSatStars (recipe, psf, sources); // non-linear PSF and EXT fit to brighter sources // replace model flux, adjust mask as needed, fit, subtract the models (full stamp) - psphotBlendFit (config, readout, sources, psf); + psphotBlendFit (config, view); // replace all sources - psphotReplaceAllSources (sources, recipe); + psphotReplaceAllSources (config, view); // linear fit to include all sources (subtract again) - psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); + psphotFitSourcesLinear (config, view, TRUE); // if we only do one pass, skip to extended source analysis @@ -176,33 +155,38 @@ // add noise for subtracted objects - psphotAddNoise (readout, sources, recipe); + psphotAddNoise (config, view); // find fainter sources (pass 2) - detections = psphotFindDetections (detections, readout, recipe); + // XXX need to distinguish old from new sources + psphotFindDetections (config, view); // remove noise for subtracted objects (ie, return to normal noise level) - psphotSubNoise (readout, sources, recipe); + // XXX this needs to operate only on the OLD sources + psphotSubNoise (config, view); // define new sources based on only the new peaks - psArray *newSources = psphotSourceStats (config, readout, detections, false); + // XXX need to distinguish old from new sources + psphotSourceStats (config, view, false); // set source type - if (!psphotRoughClass (config, view, havePSF)) { + if (!psphotRoughClass (config, view)) { psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, view); } // create full input models, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, readout, newSources, psf); + // XXX need to distinguish old from new sources + psphotGuessModels (config, view); // replace all sources so fit below applies to all at once - psphotReplaceAllSources (sources, recipe); + psphotReplaceAllSources (config, view); // merge the newly selected sources into the existing list - psphotMergeSources (sources, newSources); + // XXX old vs new???? +FIX: psphotMergeSources (sources, newSources); psFree (newSources); // linear fit to all sources - psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); + psphotFitSourcesLinear (config, view, TRUE); pass1finish: @@ -211,7 +195,6 @@ psphotSourceSize (config, view); - psphotExtendedSourceAnalysis (readout, sources, recipe); - - psphotExtendedSourceFits (readout, sources, recipe); +FIX: psphotExtendedSourceAnalysis (readout, sources, recipe); +FIX: psphotExtendedSourceFits (readout, sources, recipe); finish: @@ -221,13 +204,13 @@ // measure aperture photometry corrections - if (!psphotApResid (config, readout, sources, psf)) { +FIX: if (!psphotApResid (config, readout, sources, psf)) { psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, view); } // calculate source magnitudes - psphotMagnitudes(config, readout, view, sources, psf); - - if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) { +FIX: psphotMagnitudes(config, readout, view, sources, psf); + +FIX: if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) { psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear(); @@ -238,11 +221,10 @@ // replace background in residual image - psphotSkyReplace (config, view); +FIX: psphotSkyReplace (config, view); // drop the references to the image pixels held by each source - psphotSourceFreePixels (sources); +FIX: psphotSourceFreePixels (sources); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup(config, view); } - Index: /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutCleanup.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutCleanup.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutCleanup.c (revision 26681) @@ -1,8 +1,8 @@ # include "psphotInternal.h" -// psphotReadoutCleanup is called on exit from psphotReadout. If the last raised error is -// not a DATA error, then there was a serious problem. Only in this case, or if the fail -// on the stats measurement, do we return false -bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmDetections *detections, pmPSF *psf, psArray *sources) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view) +{ + bool status = true; // remove internal pmFPAfiles, if created @@ -12,9 +12,48 @@ } if (psErrorCodeLast() != PS_ERR_NONE) { - psFree (psf); - psFree (sources); - psFree (detections); return false; } + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotReadoutCleanupReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on psphotReadoutCleanup for PSPHOT.INPUT entry %d", i); + return false; + } + } + + // XXX move this to top of loop? + pmKapaClose (); + + return true; +} + +// psphotReadoutCleanup is called on exit from psphotReadout. If the last raised error is +// not a DATA error, then there was a serious problem. Only in this case, or if the fail +// on the stats measurement, do we return false +bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { + + pmReadout *readout, psMetadata *recipe, pmDetections *detections, pmPSF *psf, psArray *sources; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (readout, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // when psphotReadoutCleanup is called, these are not necessarily defined + psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + + // XXX where do we free these, in here (psMetadataRemove?) // use the psf-model to measure FWHM stats @@ -22,7 +61,4 @@ if (!psphotPSFstats (readout, recipe, psf)) { psError(PSPHOT_ERR_PROG, false, "Failed to measure PSF shape parameters"); - psFree (psf); - psFree (sources); - psFree (detections); return false; } @@ -32,7 +68,4 @@ if (!psphotMomentsStats (readout, recipe, sources)) { psError(PSPHOT_ERR_PROG, false, "Failed to measure Moment shape parameters"); - psFree (psf); - psFree (sources); - psFree (detections); return false; } @@ -40,4 +73,5 @@ // Check to see if the image quality was measured + // XXX not sure we want / need this test if (!psf) { bool mdok; // Status of MD lookup @@ -45,22 +79,23 @@ if (!mdok || nIQ <= 0) { psError(PSPHOT_ERR_DATA, false, "Unable to measure image quality"); - psFree (psf); - psFree (sources); - psFree (detections); return false; } } + // create an output header with stats results XXX this needs to be reworked : + // psphotImageQuality and others need to write these values to a header. that needs to be + // stored on the readout->analysis + psMetadata *header = psphotDefineHeader (recipe); // write NSTARS to the image header - psphotSetHeaderNstars (recipe, sources); - - // create an output header with stats results - psMetadata *header = psphotDefineHeader (recipe); + psphotSetHeaderNstars (header, sources); // save the results of the analysis + // this should happen way up stream (when needed?) psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.HEADER", PS_DATA_METADATA, "header stats", header); if (psf) { + // XXX this seems a little silly : we saved the psf on readout->analysis above, but now + // we are moving it to chip->analysis. // save the psf for possible output. if there was already an entry, it was loaded from external sources // the new one may have been updated or modified, so replace the existing entry. We @@ -77,12 +112,5 @@ } - // XXX move this to top of loop? - pmKapaClose (); - - psFree (detections); - psFree (psf); psFree (header); - psFree (sources); - return true; } Index: /branches/eam_branches/psphot.stack.20100120/src/psphotReplaceUnfit.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotReplaceUnfit.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotReplaceUnfit.c (revision 26681) @@ -22,5 +22,27 @@ } -bool psphotReplaceAllSources (psArray *sources, psMetadata *recipe) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotReplaceAllSourcesReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; @@ -29,7 +51,17 @@ psTimerStart ("psphot.replace"); + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (readout, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); + psAssert (sources, "missing sources?"); + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels - assert (maskVal); + psAssert (maskVal, "missing mask value?"); for (int i = 0; i < sources->n; i++) { Index: /branches/eam_branches/psphot.stack.20100120/src/psphotRoughClass.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotRoughClass.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotRoughClass.c (revision 26681) @@ -7,8 +7,27 @@ } } -bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *target, psMetadata *recipe, const bool havePSF); +// for now, let's store the detections on the readout->analysis for each readout +bool psphotRoughClass (pmConfig *config, const pmFPAview *view) +{ + bool status = true; -// 2006.02.02 : no leaks -bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, const bool havePSF) { + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; @@ -23,10 +42,17 @@ psAssert (readout, "missing readout?"); + // if we have a PSF, use the existing PSF clump region below + bool havePSF = false; + if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { + havePSF = true; + } + psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); psAssert (sources, "missing sources?"); - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - psAssert (recipe, "missing recipe?"); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification"); + return true; + } // we make this measurement on a NxM grid of regions across the readout @@ -125,20 +151,2 @@ return true; } - -// for now, let's store the detections on the readout->analysis for each readout -bool psphotRoughClass (pmConfig *config, const pmFPAview *view, const bool havePSF) -{ - bool status = true; - - int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); - psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, havePSF)) { - psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); - return false; - } - } - return true; -} Index: /branches/eam_branches/psphot.stack.20100120/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotSourceSize.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotSourceSize.c (revision 26681) @@ -15,10 +15,12 @@ } psphotSourceSizeOptions; +// local functions: +bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index); bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf); bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options); -bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); -bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); +bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal); +bool psphotMaskCosmicRayFootprintCheck (psArray *sources); // we need to call this function after sources have been fitted to the PSF model and @@ -33,4 +35,8 @@ bool status = true; + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); @@ -38,5 +44,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i)) { + if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); return false; @@ -47,6 +53,5 @@ // XXX use an internal flag to mark sources which have already been measured -// how to track the sources already measured? -bool psphotSourceSize(pmConfig *config, const pmFPAview *view, const char *filename, int index) +bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; @@ -65,10 +70,11 @@ psAssert (sources, "missing sources?"); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + return true; + } + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); psAssert (psf, "missing psf?"); - - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - psAssert (recipe, "missing recipe?"); // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -106,5 +112,5 @@ // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness - // of and object. We need to model this distribution for the PSF stars before we can test + // of an object. We need to model this distribution for the PSF stars before we can test // the significance for a specific object // XXX move this to the code that generates the PSF? @@ -122,9 +128,497 @@ psphotVisualShowSourceSize (readout, sources); psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr); - - return true; -} - -bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { + psphotVisualShowSatStars (recipe, psf, sources); + + return true; +} + +// model the apmifit distribution for the psf stars: +bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) { + + // select stats from the psf stars + psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32); + + psImageMaskType maskVal = options->maskVal | options->markVal; + + // XXX why PHOT_WEIGHT?? + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal); + } + + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); + + float apMag = -2.5*log10(source->moments->Sum); + float dMag = source->psfMag - apMag; + + psVectorAppend (Ap, 100, dMag); + psVectorAppend (ApErr, 100, source->errMag); + } + + // model the distribution as a mean or median value and a systematic error from that value: + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + psVectorStats (stats, Ap, NULL, NULL, 0); + + psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32); + for (int i = 0; i < Ap->n; i++) { + dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian; + } + + options->ApResid = stats->robustMedian; + options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05); + psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr); + + psFree (Ap); + psFree (ApErr); + psFree (stats); + psFree (dAp); + + return true; +} + +// classify sources based on the combination of psf-mag, Mxx, Myy +bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) { + + bool status; + pmPSFClump psfClump; + char regionName[64]; + + psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip"); + + int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS"); + for (int i = 0; i < nRegions; i ++) { + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); + psMetadata *regionMD = psMetadataLookupPtr (&status, readout->analysis, regionName); + psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion"); + + psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION"); + psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion"); + + // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y + psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); psAssert (status, "missing PSF.CLUMP.X"); + psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); psAssert (status, "missing PSF.CLUMP.Y"); + psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); psAssert (status, "missing PSF.CLUMP.DX"); + psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); psAssert (status, "missing PSF.CLUMP.DY"); + + if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) { + psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + continue; + } + + if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) { + psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + continue; + } + } + + return true; +} + +bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) { + + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(recipe, false); + + int Nsat = 0; + int Next = 0; + int Npsf = 0; + int Ncr = 0; + int Nmiss = 0; + int Nskip = 0; + + pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN; + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + + psImageMaskType maskVal = options->maskVal | options->markVal; + + for (psS32 i = 0 ; i < sources->n ; i++) { + + pmSource *source = (pmSource *) sources->data[i]; + + // psfClumps are found for image subregions: + // skip sources not in this region + if (source->peak->x < region->x0) continue; + if (source->peak->x >= region->x1) continue; + if (source->peak->y < region->y0) continue; + if (source->peak->y >= region->y1) continue; + + // 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; + } + + // we are basically classifying by moments; use the default if not found + psAssert (source->moments, "why is this source missing moments?"); + if (source->mode & noMoments) { + Nskip ++; + continue; + } + + psF32 Mxx = source->moments->Mxx; + psF32 Myy = source->moments->Myy; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal); + } + + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); + + float apMag = -2.5*log10(source->moments->Sum); + float dMag = source->psfMag - apMag; + float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr); + + source->extNsigma = nSigma; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + + // Anything within this region is a probably PSF-like object. Saturated stars may land + // in this region, but are detected elsewhere on the basis of their peak value. + bool isPSF = ((fabs(nSigma) < options->nSigmaApResid) && + (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && + (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY)); + if (isPSF) { + psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g PSF\t%g %g\n", + source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma, + options->nSigmaApResid,options->nSigmaMoments); + Npsf ++; + continue; + } + + // Defects may not always match CRs from peak curvature analysis + // Defects may also be marked as SATSTAR -- XXX deactivate this flag? + // XXX this rule is not great + if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) { + + psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g CR\t%g %g\n", + source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma, + options->nSigmaApResid,options->nSigmaMoments); + source->mode |= PM_SOURCE_MODE_DEFECT; + Ncr ++; + continue; + } + + // saturated star (determined in PSF fit). These may also be saturated galaxies, or + // just large saturated regions. + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g SAT\t%g %g\n", + source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma, + options->nSigmaApResid,options->nSigmaMoments); + + Nsat ++; + continue; + } + + // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)? + bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y)); + if (isEXT) { + psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g Ext\t%g %g\n", + source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma, + options->nSigmaApResid,options->nSigmaMoments); + + source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + Next ++; + continue; + } + psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g Unk\t%g %g\n", + source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma, + options->nSigmaApResid,options->nSigmaMoments); + + psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma); + Nmiss ++; + } + + psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip); + + return true; +} + +// given an object suspected to be a defect, generate a pixel mask using the Lapacian transform +// if enough of the object is detected as 'sharp', consider the object a cosmic ray +bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) { + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // 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"); + // psTrace("psphot.czw", 2, "Not calculating crNsigma due to edge\n"); + continue; + } + + // this source is thought to be a cosmic ray. flag the detection and mask the pixels + if (source->mode & PM_SOURCE_MODE_DEFECT) { + // XXX this is running slowly and is too agressive, but it more-or-less works + // psphotMaskCosmicRayCZW(readout, source, options->crMask); + + // XXX these are the old versions which we are not using + // psphotMaskCosmicRay (readout->mask, source, maskVal, crMask); + // 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; + } + + // XXX test : save the mask image + if (0) { + psphotSaveImage (NULL, readout->mask, "mask.fits"); + } + return true; +} + +// Comments by CZW 20091209 : Mechanics of how to identify CR pixels taken from "Cosmic-Ray +// Rejection by Laplacian Edge Detection" by Pieter van Dokkum, arXiv:astro-ph/0108003. This +// does no repair or recovery of the CR pixels, it only masks them out. My test code can be +// found at /data/ipp031.0/watersc1/psphot.20091209/algo_check.c +bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal) { + + // Get the actual images and information about the peak. + psImage *mask = readout->mask; + pmPeak *peak = source->peak; + pmFootprint *footprint = peak->footprint; + + int xm = footprint->bbox.x0; + int xM = footprint->bbox.x1; + int ym = footprint->bbox.y0; + int yM = footprint->bbox.y1; + + if (xm < 0) { xm = 0; } + if (ym < 0) { ym = 0; } + if (xM > mask->numCols) { xM = mask->numCols; } + if (yM > mask->numRows) { yM = mask->numRows; } + int dx = xM - xm; + int dy = yM - ym; + + // Bounding boxes are inclusive of final pixel + // XXX ensure dx,dy do not become too large here + dx++; + dy++; + + psImage *image= readout->image; + psImage *variance = readout->variance; + + int binning = 1; + float sigma_thresh = 2.0; + int iteration = 0; + int max_iter = 5; + float noise_factor = 5.0 / 4.0; // Intrinsic to the Laplacian making noise spikes spikier. + + // Temporary images. + psImage *mypix = psImageAlloc(dx,dy,image->type.type); + psImage *myvar = psImageAlloc(dx,dy,image->type.type); + psImage *binned = psImageAlloc(dx * binning,dy * binning,image->type.type); + psImage *conved = psImageAlloc(dx * binning,dy * binning,image->type.type); + psImage *edges = psImageAlloc(dx,dy,image->type.type); + psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK); + + int x,y; + // Load my copy of things. + for (x = 0; x < dx; x++) { + for (y = 0; y < dy; y++) { + psImageSet(mypix,x,y,psImageGet(image,x+xm,y+ym)); + psImageSet(myvar,x,y,psImageGet(variance,x+xm,y+ym)); + mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00; + } + } + // Mask so I can see on the output image where the footprint is. + for (int i = 0; i < footprint->spans->n; i++) { + pmSpan *sp = footprint->spans->data[i]; + for (int j = sp->x0; j <= sp->x1; j++) { + y = sp->y - ym; + x = j - xm; + mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01; + } + } + + int CRpix_count = 0; + do { + CRpix_count = 0; + // Zero out my temp images. + for (x = 0; x < binning * dx; x++) { + for (y = 0; y < binning * dy; y++) { + psImageSet(binned,x,y,0.0); + psImageSet(conved,x,y,0.0); + psImageSet(edges,x/binning,y/ binning,0.0); + } + } + // Make subsampled image. Maybe this should be called "unbinned" or something + for (x = 0; x < binning * dx; x++) { + for (y = 0; y < binning * dy; y++) { + psImageSet(binned,x,y,psImageGet(mypix,x / binning,y / binning)); + } + } + // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero + for (x = 1; x < dx - 1; x++) { + for (y = 1; y < dy - 1; y++) { + psImageSet(conved,x,y,psImageGet(binned,x,y) - 0.25 * + (psImageGet(binned,x-1,y) + psImageGet(binned,x+1,y) + + psImageGet(binned,x,y-1) + psImageGet(binned,x,y+1))); + if (psImageGet(conved,x,y) < 0.0) { + psImageSet(conved,x,y,0.0); + } + } + } + // Create an edge map by rebinning + for (x = 0; x < binning * dx; x++) { + for (y = 0; y < binning * dy; y++) { + psImageSet(edges,x / binning, y / binning, + psImageGet(edges, x / binning, y / binning) + + psImageGet(conved,x,y)); + } + } + // Modify my mask if we're above the significance threshold + for (x = 0; x < dx; x++) { + for (y = 0; y < dy; y++) { + if ( psImageGet(edges,x,y) / (binning * sqrt(noise_factor * psImageGet(myvar,x,y))) > sigma_thresh ) { + if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01) { + mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40; + CRpix_count++; + } + } + } + } + + // "Repair" Masked pixels for the next round. + for (x = 1; x < dx - 1; x++) { + for (y = 1; y < dy - 1; y++) { + if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) { + psImageSet(mypix,x,y,0.25 * + (psImageGet(mypix,x-1,y) + psImageGet(mypix,x+1,y) + + psImageGet(mypix,x,y-1) + psImageGet(mypix,x,y+1))); + } + } + } + +# if 0 + if ((psTraceGetLevel("psphot.czw") >= 2)&&(abs(xm - 2770) < 5)&&(abs(ym - 2581) < 5)&&(iteration == 0)) { + psTrace("psphot.czw",2,"TRACEMOTRON %d %d %d %d %d\n",xm,ym,dx,dy,iteration); + psphotSaveImage (NULL, mypix, "czw.pix.fits"); + psphotSaveImage (NULL, myvar, "czw.var.fits"); + psphotSaveImage (NULL, binned, "czw.binn.fits"); + psphotSaveImage (NULL, conved, "czw.conv.fits"); + psphotSaveImage (NULL, edges, "czw.edge.fits"); + psphotSaveImage (NULL, mymask, "czw.mask.fits"); + } +# endif + + psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration,CRpix_count); + iteration++; + } while ((iteration < max_iter)&&(CRpix_count > 0)); + + // A solitary masked pixel is likely a lie. Remove those. + for (x = 0; x < dx; x++) { + for (y = 0; y < dy; y++) { + if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) { + if ((x-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) { + continue; + } + if ((y-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) { + continue; + } + if ((x+1 < dx)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) { + continue; + } + if ((y+1 < dy)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) { + continue; + } + mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40; + } + } + } + + // transfer temporary mask to real mask + for (x = 0; x < dx; x++) { + for (y = 0; y < dy; y++) { + if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= maskVal; + } + } + } + + // XXX if we decide this REALLY is a cosmic ray, set the CR_LIMIT bit + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + + psFree(mypix); + psFree(myvar); + psFree(binned); + psFree(conved); + psFree(edges); + psFree(mymask); + + return true; +} + +bool psphotMaskCosmicRayFootprintCheck (psArray *sources) { + + 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"); + } + } + return true; +} + +/**** ------ old versions of cosmic ray masking ----- ****/ + +// 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 @@ -139,4 +633,6 @@ 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); @@ -145,21 +641,23 @@ 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; - } - } + 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; } @@ -228,223 +726,4 @@ } } - return true; -} - -// model the apmifit distribution for the psf stars: -bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) { - - // select stats from the psf stars - psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32); - psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32); - - psImageMaskType maskVal = options->maskVal | options->markVal; - - // XXX why PHOT_WEIGHT?? - pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; - - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; - - // replace object in image - if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { - pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal); - } - - // clear the mask bit and set the circular mask pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); - - // XXX can we test if psfMag is set and calculate only if needed? - pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal - - // clear the mask bit - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); - - // re-subtract the object, leave local sky - pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); - - float apMag = -2.5*log10(source->moments->Sum); - float dMag = source->psfMag - apMag; - - psVectorAppend (Ap, dMag); - psVectorAppend (ApErr, source->errMag); - } - - // model the distribution as a mean or median value and a systematic error from that value: - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - psVectorStats (stats, Ap, NULL, NULL, 0); - - psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32); - for (int i = 0; i < Ap->n; i++) { - dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian; - } - - options->ApResid = stats->robustMedian; - options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05); - // XXX this is quite arbitrary... - if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.01; - psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr); - - psFree (Ap); - psFree (ApErr); - psFree (stats); - psFree (dAp); - - return true; -} - -// classify sources based on the combination of psf-mag, Mxx, Myy -bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) { - - bool status; - pmPSFClump psfClump; - char regionName[64]; - - psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip"); - - int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS"); - for (int i = 0; i < nRegions; i ++) { - snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); - psMetadata *regionMD = psMetadataLookupPtr (&status, readout->analysis, regionName); - psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion"); - - psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION"); - psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion"); - - // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y - psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); psAssert (status, "missing PSF.CLUMP.X"); - psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); psAssert (status, "missing PSF.CLUMP.Y"); - psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); psAssert (status, "missing PSF.CLUMP.DX"); - psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); psAssert (status, "missing PSF.CLUMP.DY"); - - if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) { - psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); - continue; - } - - if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) { - psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); - continue; - } - } - - return true; -} - -bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) { - - PS_ASSERT_PTR_NON_NULL(sources, false); - PS_ASSERT_PTR_NON_NULL(recipe, false); - - int Nsat = 0; - int Next = 0; - int Npsf = 0; - int Ncr = 0; - int Nmiss = 0; - int Nskip = 0; - - pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN; - pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; - - psImageMaskType maskVal = options->maskVal | options->markVal; - - for (psS32 i = 0 ; i < sources->n ; i++) { - - pmSource *source = (pmSource *) sources->data[i]; - - // psfClumps are found for image subregions: - // skip sources not in this region - if (source->peak->x < region->x0) continue; - if (source->peak->x >= region->x1) continue; - if (source->peak->y < region->y0) continue; - if (source->peak->y >= region->y1) continue; - - // 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; - } - - // we are basically classifying by moments; use the default if not found - psAssert (source->moments, "why is this source missing moments?"); - if (source->mode & noMoments) { - Nskip ++; - continue; - } - - psF32 Mxx = source->moments->Mxx; - psF32 Myy = source->moments->Myy; - - // replace object in image - if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { - pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal); - } - - // clear the mask bit and set the circular mask pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); - - // XXX can we test if psfMag is set and calculate only if needed? - pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal - - // clear the mask bit - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); - - // re-subtract the object, leave local sky - pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); - - float apMag = -2.5*log10(source->moments->Sum); - float dMag = source->psfMag - apMag; - float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr); - - source->extNsigma = nSigma; - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; - - // Anything within this region is a probably PSF-like object. Saturated stars may land - // in this region, but are detected elsewhere on the basis of their peak value. - bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY); - if (isPSF) { - Npsf ++; - continue; - } - - // Defects may not always match CRs from peak curvature analysis - // Defects may also be marked as SATSTAR -- XXX deactivate this flag? - // XXX this rule is not great - if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) { - source->mode |= PM_SOURCE_MODE_DEFECT; - Ncr ++; - continue; - } - - // saturated star (determined in PSF fit). These may also be saturated galaxies, or - // just large saturated regions. - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - Nsat ++; - continue; - } - - // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)? - bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y)); - if (isEXT) { - source->mode |= PM_SOURCE_MODE_EXT_LIMIT; - Next ++; - continue; - } - - psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma); - Nmiss ++; - } - - psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip); - return true; } @@ -524,5 +803,6 @@ } -bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) { +// 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?) @@ -656,2 +936,3 @@ return true; } + Index: /branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c (revision 26681) @@ -1,26 +1,42 @@ # include "psphotInternal.h" -bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources); - -bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool setWindow) { - - bool status = false; - psArray *sources = NULL; - - psTimerStart ("psphot.stats"); - - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest - psAssert (readout, "missing file?"); - - pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); - psAssert (readout, "missing readout?"); - - pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psAssert (detections, "missing detections?"); +// convert detections to sources and measure their basic properties (moments, local sky, sky variance) +bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow) +{ + bool status = true; // select the appropriate recipe information psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); psAssert (recipe, "missing recipe?"); + + int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotSourceStatsReadout (config, view, "PSPHOT.INPUT", i, recipe, setWindow)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool setWindow) { + + bool status = false; + psArray *sources = NULL; + + psTimerStart ("psphot.stats"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // 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?"); // determine the number of allowed threads @@ -38,5 +54,5 @@ char *breakPt = psMetadataLookupStr (&status, recipe, "BREAK_POINT"); - if (!status) return false; + psAssert (status, "missing BREAK_POINT?"); psArray *peaks = detections->peaks; @@ -48,4 +64,16 @@ // generate the array of sources, define the associated pixel sources = psArrayAllocEmpty (peaks->n); + + // if there are no peaks, we save the empty source array and return + if (!peaks->n) { + // save sources on the readout->analysis + if (!psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources)) { + psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout"); + psFree(sources); + return false; + } + psFree(sources); + return true; + } for (int i = 0; i < peaks->n; i++) { @@ -485,19 +513,2 @@ // otherwise it only contains the new peaks. -// for now, let's store the detections on the readout->analysis for each readout -bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow) -{ - bool status = true; - - int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); - psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - if (!psphotFindDetectionsReadout (config, view, "PSPHOT.INPUT", i, setWindow)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i); - return false; - } - } - return true; -} Index: /branches/eam_branches/psphot.stack.20100120/src/psphotSubtractBackground.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotSubtractBackground.c (revision 26680) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotSubtractBackground.c (revision 26681) @@ -4,5 +4,5 @@ // generate the median in NxN boxes, clipping heavily // linear interpolation to generate full-scale model -bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) +bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status = true; @@ -25,8 +25,4 @@ pmReadout *model = pmFPAviewThisReadout (view, modelFile->fpa); assert (model); - - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - assert (recipe); // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -114,4 +110,8 @@ // the pmReadout selected in this function are all view on entries in config->files + // display the backsub and backgnd images + // move this inthe the subtract background loop + psphotVisualShowBackground (config, view, readout); + npass ++; return true; @@ -122,4 +122,8 @@ bool status = false; + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); @@ -127,5 +131,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotSubtractBackgroundReadout (config, view, "PSPHOT.INPUT", i)) { + if (!psphotSubtractBackgroundReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to subtract background for PSPHOT.INPUT entry %d", i); return false;