Index: trunk/ppSub/src/Makefile.am =================================================================== --- trunk/ppSub/src/Makefile.am (revision 24833) +++ trunk/ppSub/src/Makefile.am (revision 24862) @@ -45,5 +45,6 @@ ppSubReadoutStats.c \ ppSubReadoutSubtract.c \ - ppSubSetMasks.c + ppSubSetMasks.c \ + ppSubThreshold.c ppSubKernel_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSUB_CFLAGS) Index: trunk/ppSub/src/ppSub.h =================================================================== --- trunk/ppSub/src/ppSub.h (revision 24833) +++ trunk/ppSub/src/ppSub.h (revision 24862) @@ -76,4 +76,8 @@ /// looking up the parameters in the recipe and supplying them to the function pmSubtractionMatch() bool ppSubMatchPSFs(ppSubData *data ///< Processing data + ); + +/// Threshold low pixels in image +bool ppSubLowThreshold(ppSubData *data ///< Processing data ); Index: trunk/ppSub/src/ppSubLoop.c =================================================================== --- trunk/ppSub/src/ppSubLoop.c (revision 24833) +++ trunk/ppSub/src/ppSubLoop.c (revision 24862) @@ -62,4 +62,9 @@ if (!ppSubFilesIterateUp(config, PPSUB_FILES_INPUT)) { psError(PPSUB_ERR_IO, false, "Unable to close input files."); + return false; + } + + if (!ppSubLowThreshold(data)) { + psError(PS_ERR_UNKNOWN, false, "Unable to threshold images."); return false; } Index: trunk/ppSub/src/ppSubMatchPSFs.c =================================================================== --- trunk/ppSub/src/ppSubMatchPSFs.c (revision 24833) +++ trunk/ppSub/src/ppSubMatchPSFs.c (revision 24862) @@ -18,5 +18,4 @@ #include #include -#include #include "ppSub.h" Index: trunk/ppSub/src/ppSubReadout.c =================================================================== --- trunk/ppSub/src/ppSubReadout.c (revision 24833) +++ (revision ) @@ -1,67 +1,0 @@ -/** @file ppSubReadout.c - * - * @brief - * - * @ingroup ppSub - * - * @author IfA - * @version $Revision: 1.113 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 00:31:20 $ - * Copyright 2009 Institute for Astronomy, University of Hawaii - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include - -#include "ppSub.h" - -bool ppSubReadout(const char *name, bool reverse, ppSubData *data, const pmFPAview *view) -{ - psAssert(data, "Require processing data"); - pmConfig *config = data->config; // Configuration - psAssert(config, "Require configuration"); - - psAssert(name, "Require name"); - psAssert(view, "Require view"); - - if (!ppSubDefineOutput(name, config, data, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define output."); - return false; - } - - if (!data->quality && !ppSubMakePSF(name, config, data, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate PSF."); - return false; - } - - if (!ppSubReadoutSubtract(name, reverse, config, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to subtract images."); - return false; - } - - // Higher order background subtraction using psphot - if (!ppSubBackground(name, config, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); - return false; - } - - if (!data->quality && data->!ppSubReadoutPhotometry(name, config, data, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry."); - return false; - } - - if (!ppSubReadoutUpdate(name, config, data, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to update."); - return false; - } - - - - - return true; -} Index: trunk/ppSub/src/ppSubReadoutSubtract.c =================================================================== --- trunk/ppSub/src/ppSubReadoutSubtract.c (revision 24833) +++ trunk/ppSub/src/ppSubReadoutSubtract.c (revision 24862) @@ -44,20 +44,4 @@ subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); } - -// XXX this should be replaced by SAVE options to the inConv and refConf pmFPAfiles -#ifdef TESTING - { - pmReadoutMaskApply(minuend, maskVal); - psFits *fits = psFitsOpen("minuend.fits", "w"); - psFitsWriteImage(fits, NULL, minuend->image, 0, NULL); - psFitsClose(fits); - } - { - pmReadoutMaskApply(subtrahend, maskVal); - psFits *fits = psFitsOpen("subtrahend.fits", "w"); - psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL); - psFitsClose(fits); - } -#endif // Do the actual subtraction Index: trunk/ppSub/src/ppSubThreshold.c =================================================================== --- trunk/ppSub/src/ppSubThreshold.c (revision 24862) +++ trunk/ppSub/src/ppSubThreshold.c (revision 24862) @@ -0,0 +1,114 @@ +/** @file ppSubThreshold.c + * + * @brief Mask pixels below threshold + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.3 $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppSub.h" + +// Apply low threshold to image +static bool lowThreshold( + pmReadout *ro, // Readout to threshold + float thresh, // Threshold to apply + psImageMaskType maskIgnore, // Ignore pixels with this mask + psImageMaskType maskThresh, // Give pixels this mask if below threshold + const char *description // Description of image + ) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + PM_ASSERT_READOUT_MASK(ro, false); + PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false); + + psImage *image = ro->image; // Image + psImage *mask = ro->mask; // Mask + int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(stats, NULL, image, mask, maskIgnore, rng)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine threshold."); + psFree(rng); + psFree(stats); + return false; + } + psFree(rng); + + float threshold = stats->robustMedian - thresh * stats->robustStdev; // Threshold below which to clip + psFree(stats); + psLogMsg("ppSub", PS_LOG_INFO, "Masking pixels below %f in %s", threshold, description); + + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskIgnore) { + continue; + } + if (image->data.F32[y][x] < threshold) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskThresh; + } + } + } + + return true; +} + + +bool ppSubLowThreshold(ppSubData *data) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psImageMaskType maskThresh = pmConfigMaskGet("LOW", config); // Bits to mask for low pixels + + float thresh = psMetadataLookupF32(NULL, recipe, "LOW.THRESHOLD"); // Threshold for masking + if (!isfinite(thresh)) { + return true; + } + + pmFPAview *view = ppSubViewReadout(); // View to readout + + // Input images + pmReadout *in = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input image + if (!in) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find readout."); + return false; + } + if (!lowThreshold(in, thresh, maskVal, maskThresh, "input convolved image")) { + psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image."); + return false; + } + + pmReadout *ref = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference image + if (!ref) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find readout."); + return false; + } + if (!lowThreshold(ref, thresh, maskVal, maskThresh, "reference convolved image")) { + psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image."); + return false; + } + + psFree(view); + + return true; +}