Index: branches/pap/ppSub/src/Makefile.am =================================================================== --- branches/pap/ppSub/src/Makefile.am (revision 23948) +++ branches/pap/ppSub/src/Makefile.am (revision 25027) @@ -30,4 +30,5 @@ ppSubVersion.c \ ppSubBackground.c \ + ppSubVarianceRescale.c \ ppSubCamera.c \ ppSubData.c \ @@ -44,5 +45,6 @@ ppSubReadoutStats.c \ ppSubReadoutSubtract.c \ - ppSubSetMasks.c + ppSubSetMasks.c \ + ppSubThreshold.c ppSubKernel_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSUB_CFLAGS) Index: branches/pap/ppSub/src/ppSub.h =================================================================== --- branches/pap/ppSub/src/ppSub.h (revision 23948) +++ branches/pap/ppSub/src/ppSub.h (revision 25027) @@ -24,4 +24,5 @@ #define PPSUB_RECIPE "PPSUB" /// Name of the recipe to use +#define WCS_TOLERANCE 0.001 // Tolerance for WCS // Output files, for activation/deactivation @@ -77,4 +78,8 @@ ); +/// Threshold low pixels in image +bool ppSubLowThreshold(ppSubData *data ///< Processing data + ); + /// Generate the output readout and pass the kernel info to the header bool ppSubDefineOutput(const char *name,///< Name of output to define @@ -98,4 +103,8 @@ /// Higher-order background subtraction bool ppSubBackground(pmConfig *config ///< Configuration + ); + +/// Perform Variance correction (rescale within a modest range) +bool ppSubVarianceRescale(pmConfig *config ///< Configuration ); Index: branches/pap/ppSub/src/ppSubArguments.c =================================================================== --- branches/pap/ppSub/src/ppSubArguments.c (revision 23948) +++ branches/pap/ppSub/src/ppSubArguments.c (revision 25027) @@ -84,4 +84,5 @@ psMetadataAddS32(arguments, PS_LIST_TAIL, "-threads", 0, "Number of threads", 0); psMetadataAddStr(arguments, PS_LIST_TAIL, "-dumpconfig", 0, "file to dump configuration to", NULL); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-convolve", 0, "Image to convolve [1 or 2]", 0); psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Perform photometry?", NULL); psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", NULL); Index: branches/pap/ppSub/src/ppSubCamera.c =================================================================== --- branches/pap/ppSub/src/ppSubCamera.c (revision 23948) +++ branches/pap/ppSub/src/ppSubCamera.c (revision 25027) @@ -163,4 +163,25 @@ + // Now that the camera has been determined, we can read the recipe + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + if (!recipe) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSUB_RECIPE); + return false; + } + if (psMetadataLookupBool(NULL, config->arguments, "-photometry")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "PHOTOMETRY", PS_META_REPLACE, + "Perform photometry?", true); + } + if (psMetadataLookupBool(NULL, config->arguments, "-inverse")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "INVERSE", PS_META_REPLACE, + "Generate inverse subtractions?", true); + } + + data->inverse = psMetadataLookupBool(NULL, recipe, "INVERSE"); + data->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); + + bool mdok; // Status of MD lookup + bool saveConv = psMetadataLookupBool(&mdok, recipe, "SAVE.CONVOLVED"); // Save convolved images? + // Convolved input image pmFPAfile *inConvImage = defineOutputFile(config, input, true, "PPSUB.INPUT.CONV", PM_FPA_FILE_IMAGE); @@ -171,4 +192,8 @@ return false; } + if (saveConv) { + inConvImage->save = true; + inConvMask->save = true; + } if (inVar) { pmFPAfile *inConvVar = defineOutputFile(config, inConvImage, false, "PPSUB.INPUT.CONV.VARIANCE", @@ -177,4 +202,7 @@ psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); return false; + } + if (saveConv) { + inConvVar->save = true; } } @@ -188,4 +216,8 @@ return false; } + if (saveConv) { + refConvImage->save = true; + refConvMask->save = true; + } if (refVar) { pmFPAfile *refConvVar = defineOutputFile(config, refConvImage, false, "PPSUB.REF.CONV.VARIANCE", @@ -195,25 +227,8 @@ return false; } - } - - - // Now that the camera has been determined, we can read the recipe - psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim - if (!recipe) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSUB_RECIPE); - return false; - } - if (psMetadataLookupBool(NULL, config->arguments, "-photometry")) { - psMetadataAddBool(recipe, PS_LIST_TAIL, "PHOTOMETRY", PS_META_REPLACE, - "Perform photometry?", true); - } - if (psMetadataLookupBool(NULL, config->arguments, "-inverse")) { - psMetadataAddBool(recipe, PS_LIST_TAIL, "INVERSE", PS_META_REPLACE, - "Generate inverse subtractions?", true); - } - - data->inverse = psMetadataLookupBool(NULL, recipe, "INVERSE"); - data->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); - + if (saveConv) { + refConvVar->save = true; + } + } // Output image Index: branches/pap/ppSub/src/ppSubKernel.c =================================================================== --- branches/pap/ppSub/src/ppSubKernel.c (revision 23948) +++ branches/pap/ppSub/src/ppSubKernel.c (revision 25027) @@ -29,5 +29,5 @@ } - psTraceSetLevel("psModules.imcombine", 7); + (void) psTraceSetLevel("psModules.imcombine", 7); const char *inName = argv[1]; // Input file name Index: branches/pap/ppSub/src/ppSubLoop.c =================================================================== --- branches/pap/ppSub/src/ppSubLoop.c (revision 23948) +++ branches/pap/ppSub/src/ppSubLoop.c (revision 25027) @@ -65,4 +65,9 @@ } + if (!ppSubLowThreshold(data)) { + psError(PS_ERR_UNKNOWN, false, "Unable to threshold images."); + return false; + } + // Set up subtraction files if (!ppSubFilesIterateDown(config, PPSUB_FILES_SUB)) { @@ -95,4 +100,10 @@ if (!ppSubBackground(config)) { psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); + return false; + } + + // Perform Variance correction (rescale within a modest range) + if (!ppSubVarianceRescale(config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to rescale variance."); return false; } Index: branches/pap/ppSub/src/ppSubMatchPSFs.c =================================================================== --- branches/pap/ppSub/src/ppSubMatchPSFs.c (revision 23948) +++ branches/pap/ppSub/src/ppSubMatchPSFs.c (revision 25027) @@ -18,7 +18,23 @@ #include #include -#include #include "ppSub.h" + +// Normalise a region on an image +static void normaliseRegion(psImage *image, // Image to normalise + const psRegion *region, // Region of image to normalise + float norm // Normalisation + ) +{ + if (!image) { + return; + } + psAssert(region, "Expect region"); + psImage *subImage = psImageSubset(image, *region); // Sub-image + psBinaryOp(subImage, subImage, "*", psScalarAlloc(norm, PS_TYPE_F32)); + psFree(subImage); + return; +} + bool ppSubMatchPSFs(ppSubData *data) @@ -60,6 +76,8 @@ pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); - psArray *inSources = psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.SOURCES"); // Source list - psArray *refSources = psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.SOURCES"); // Source list + psArray *inSources = inSourceRO ? psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.SOURCES") : + NULL; // Source list from input image + psArray *refSources = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.SOURCES") : + NULL ; // Source list from reference image psArray *sources = NULL; // Merged list of sources @@ -92,5 +110,4 @@ float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps - const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps const char *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type @@ -135,5 +152,24 @@ bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution? - pmSubtractionMode subMode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtracn mode + pmSubtractionMode subMode; // Subtraction mode + if (dual) { + subMode = PM_SUBTRACTION_MODE_DUAL; + } else { + int convolve = psMetadataLookupS32(NULL, config->arguments, "-convolve"); // Image number to convolve + switch (convolve) { + case 0: + subMode = PM_SUBTRACTION_MODE_UNSURE; + break; + case 1: + subMode = PM_SUBTRACTION_MODE_1; + break; + case 2: + subMode = PM_SUBTRACTION_MODE_2; + break; + default: + psErrorStackPrint(stderr, "Invalid value for -convolve"); + return false; + } + } int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads @@ -149,5 +185,5 @@ } else { success = pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, - spacing, threshold, sources, stampsName, type, size, order, + spacing, threshold, sources, data->stamps, type, size, order, widths, orders, inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, sys, maskVal, @@ -165,4 +201,9 @@ ppSubDataQuality(data, error, PPSUB_FILES_ALL); return true; + } else if (error == PM_ERR_SMALL_AREA) { + psErrorStackPrint(stderr, "Insufficient area for PSF matching"); + psWarning("Insufficient area for PSF matching --- suspect bad data quality."); + ppSubDataQuality(data, error, PPSUB_FILES_ALL); + return true; } else { psError(PS_ERR_UNKNOWN, false, "Unable to match images."); @@ -171,4 +212,49 @@ } + // Need to be careful with the normalisation + // We will normalise everything to the normalisation of the *input* image + { + // Since the entries are MULTI, we have to retrieve them differently + psMetadataIterator *regIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, + "^" PM_SUBTRACTION_ANALYSIS_REGION "$"); + psMetadataIterator *modeIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, + "^" PM_SUBTRACTION_ANALYSIS_MODE "$"); + psMetadataIterator *normIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, + "^" PM_SUBTRACTION_ANALYSIS_NORM "$"); + psMetadataItem *regItem; // Item with region + while ((regItem = psMetadataGetAndIncrement(regIter))) { + psAssert(regItem->type == PS_DATA_REGION && regItem->data.V, "Expect region type"); + psRegion *region = regItem->data.V; // Region of interest + psMetadataItem *modeItem = psMetadataGetAndIncrement(modeIter); // Item with mode + psAssert(modeItem && modeItem->type == PS_TYPE_S32, "Expect subtraction mode"); + pmSubtractionMode mode = modeItem->data.S32; // Subtraction mode + psMetadataItem *normItem = psMetadataGetAndIncrement(normIter); // Item with normalisation + psAssert(normItem && normItem->type == PS_TYPE_F32 && isfinite(normItem->data.F32), + "Expect normalisation"); + float norm = normItem->data.F32; // Normalisation + + switch (mode) { + case PM_SUBTRACTION_MODE_1: // Convolved the input to match template + case PM_SUBTRACTION_MODE_DUAL: // Convolved both; template should have flux conserved + psLogMsg("ppSub", PS_LOG_INFO, "Correcting image for normalisation of %f\n", norm); + normaliseRegion(inConv->image, region, 1.0 / norm); + normaliseRegion(refConv->image, region, 1.0 / norm); + normaliseRegion(inConv->variance, region, 1.0 / PS_SQR(norm)); + normaliseRegion(refConv->variance, region, 1.0 / PS_SQR(norm)); + break; + case PM_SUBTRACTION_MODE_2: // Convolved the template to match input + // We're already happy! + psLogMsg("ppSub", PS_LOG_INFO, "Image normalisation is correct\n"); + break; + default: + psAbort("Invalid subtraction mode: %x", mode); + } + } + psFree(regIter); + psFree(modeIter); + psFree(normIter); + } + + pmConceptsCopyFPA(inConv->parent->parent->parent, inRO->parent->parent->parent, true, true); pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true); Index: branches/pap/ppSub/src/ppSubReadout.c =================================================================== --- branches/pap/ppSub/src/ppSubReadout.c (revision 23948) +++ (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: branches/pap/ppSub/src/ppSubReadoutInverse.c =================================================================== --- branches/pap/ppSub/src/ppSubReadoutInverse.c (revision 23948) +++ branches/pap/ppSub/src/ppSubReadoutInverse.c (revision 25027) @@ -7,5 +7,4 @@ #include "ppSub.h" - bool ppSubReadoutInverse(pmConfig *config) @@ -24,5 +23,23 @@ invRO->data_exists = invRO->parent->data_exists = invRO->parent->parent->data_exists = true; - pmConceptsCopyFPA(invRO->parent->parent->parent, outRO->parent->parent->parent, true, true); + // Get concepts from reference + pmFPAfile *refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF.CONV"); // File with concepts + pmFPA *invFPA = invRO->parent->parent->parent; // Inverse FPA + pmConceptsCopyFPA(invFPA, refFile->fpa, true, true); + + // Get astrometry from (forward) subtraction + pmChip *outChip = outRO->parent->parent; // Output chip + pmFPA *outFPA = outChip->parent; // Output FPA + pmChip *invChip = invRO->parent->parent; // Inverse chip + pmHDU *invHDU = invFPA->hdu; // Inverse HDU + if (!pmAstromWriteWCS(invHDU->header, outFPA, outChip, WCS_TOLERANCE)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to PPSUB.INVERSE."); + return false; + } + // Read from newly written astrometry so that it exists in the "inverse" FPA (for sources) + if (!pmAstromReadWCS(invFPA, invChip, invHDU->header, 1.0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry."); + return false; + } return true; Index: branches/pap/ppSub/src/ppSubReadoutPhotometry.c =================================================================== --- branches/pap/ppSub/src/ppSubReadoutPhotometry.c (revision 23948) +++ branches/pap/ppSub/src/ppSubReadoutPhotometry.c (revision 25027) @@ -65,4 +65,7 @@ // equivalent to the minuend image. pmReadout *inRO = pmFPAfileThisReadout(config->files, view, name); // Readout with image and sources + if (psMetadataLookup(inRO->analysis, "PSPHOT.SOURCES")) { + psMetadataRemoveKey(inRO->analysis, "PSPHOT.SOURCES"); + } pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file @@ -92,6 +95,11 @@ } + // If no sources were found, there's no error, but we want to trigger 'bad quality' + psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources + if (!sources) { + ppSubDataQuality(data, PSPHOT_ERR_DATA, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); + } + if (data->stats) { - psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources bool mdok; int numSources = psMetadataLookupS32(&mdok, data->stats, "NUM_SOURCES"); // Number of sources Index: branches/pap/ppSub/src/ppSubReadoutSubtract.c =================================================================== --- branches/pap/ppSub/src/ppSubReadoutSubtract.c (revision 23948) +++ branches/pap/ppSub/src/ppSubReadoutSubtract.c (revision 25027) @@ -20,6 +20,4 @@ #include "ppSub.h" - -#define WCS_TOLERANCE 0.001 // Tolerance for WCS bool ppSubReadoutSubtract(pmConfig *config) @@ -46,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: branches/pap/ppSub/src/ppSubThreshold.c =================================================================== --- branches/pap/ppSub/src/ppSubThreshold.c (revision 25027) +++ branches/pap/ppSub/src/ppSubThreshold.c (revision 25027) @@ -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; +} Index: branches/pap/ppSub/src/ppSubVarianceRescale.c =================================================================== --- branches/pap/ppSub/src/ppSubVarianceRescale.c (revision 25027) +++ branches/pap/ppSub/src/ppSubVarianceRescale.c (revision 25027) @@ -0,0 +1,205 @@ +/** @file ppSubVarianceRescale.c + * + * @brief measure the background signal/noise distribution and rescale the variance if needed + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.7 $ $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 + +#include "ppSub.h" + +bool ppSubVarianceRescale(pmConfig *config) +{ + psAssert(config, "Require configuration"); + + bool mdok; // Status of metadata lookups + + psMetadata *ppSubRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub + psAssert(ppSubRecipe, "Need PPSUB recipe"); + + if (!psMetadataLookupBool(&mdok, ppSubRecipe, "VARIANCE.RESCALE")) return true; + + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask + + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image + + psImage *image = outRO->image; // Image of interest + psImage *variance = outRO->variance; // Variance image + psImage *mask = outRO->mask; // Mask of interest + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + + psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Generating the signal-to-noise image"); + + // generate an image representing the signal/noise: + psImage *SN = psImageCopy (NULL, outRO->image, PS_TYPE_F32); + + int nx = image->numCols; // Size of image + int ny = image->numRows; // Size of image + + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskBad) { + SN->data.F32[y][x] = 0.0; + continue; + } + if (!isfinite(image->data.F32[y][x])) { + SN->data.F32[y][x] = 0.0; + continue; + } + if (!isfinite(variance->data.F32[y][x]) || (variance->data.F32[y][x] < 0.0)) { + SN->data.F32[y][x] = 0.0; + continue; + } + SN->data.F32[y][x] = image->data.F32[y][x] / sqrt(variance->data.F32[y][x]); + } + } + + // these should not be chosen by the user: only a ROBUST version makes sense + psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + + int Nsubset = psMetadataLookupS32(&mdok, ppSubRecipe, "VARIANCE.RESCALE.NSAMPLE"); + psAssert (mdok, "missing VARIANCE.RESCALE.NSAMPLE in ppSub recipe"); + + const int Npixels = nx*ny; // Total number of pixels + Nsubset = PS_MIN(Nsubset, Npixels); // Number of pixels in subset + psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Searching for %d pixels in S/N image", Nsubset); + + psVector *values = psVectorAllocEmpty(Nsubset, PS_TYPE_F32); // Vector containing subsample + + // Minimum and maximum values + float min = +PS_MAX_F32; + float max = -PS_MAX_F32; + + // select a subset of the image pixels to measure the stats + long n = 0; // Number of actual pixels in subset + if (Nsubset >= Npixels) { + // if we have an image smaller than Nsubset, just loop over the image pixels + for (int iy = 0; iy < ny; iy++) { + for (int ix = 0; ix < nx; ix++) { + if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskBad)) { + continue; + } + + float value = SN->data.F32[iy][ix]; + min = PS_MIN(value, min); + max = PS_MAX(value, max); + values->data.F32[n] = value; + n++; + } + } + } else { + for (long i = 0; i < Nsubset; i++) { + double frnd = psRandomUniform(rng); + int pixel = Npixels * frnd; + int ix = pixel % nx; + int iy = pixel / nx; + + if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskBad)) { + continue; + } + + float value = SN->data.F32[iy][ix]; + min = PS_MIN(value, min); + max = PS_MAX(value, max); + values->data.F32[n] = value; + n++; + } + } + if (n < 0.01*Nsubset) { + psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld)", n); + psFree(values); + return false; + } + values->n = n; + psFree (SN); + + psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Using for %ld pixels in S/N image", values->n); + + if (!psVectorStats (stats, values, NULL, NULL, 0)) { + if (psTraceGetLevel("psLib.imageops") >= 5) { + FILE *f = fopen ("vector.dat", "w"); + int fd = fileno(f); + p_psVectorPrint (fd, values, "values"); + fclose (f); + } + psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background " + "(%dx%d, (row0,col0) = (%d,%d)", + image->numRows, image->numCols, image->row0, image->col0); + psFree(values); + return false; + } + if (psTraceGetLevel("psLib.imageops") >= 6) { + FILE *f = fopen ("vector.dat", "w"); + int fd = fileno(f); + p_psVectorPrint (fd, values, "values"); + fclose (f); + } + psFree (values); + psLogMsg ("ppSub", PS_LOG_INFO, "background stdev %f\n", stats->robustStdev); + + float minValid = psMetadataLookupF32(&mdok, ppSubRecipe, "VARIANCE.RESCALE.MIN.VALID"); + psAssert (mdok, "missing VARIANCE.RESCALE.MIN.VALID in ppSub recipe"); + + float maxValid = psMetadataLookupF32(&mdok, ppSubRecipe, "VARIANCE.RESCALE.MAX.VALID"); + psAssert (mdok, "missing VARIANCE.RESCALE.MAX.VALID in ppSub recipe"); + + psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Stdev of S/N image: %f (mean: %f)", stats->robustStdev, stats->robustMedian); + + // valid range of correction; 0.66 to 1.5 (skip if in range 0.98 to 1.02) + if ((stats->robustStdev < minValid) || (stats->robustStdev > maxValid)) { + psWarning ("background stdev (%f) is out of allowed range; not correcting\n", stats->robustStdev); + // XXX set quality to poor + psFree (stats); + return true; + } + + // valid range of correction; 0.66 to 1.5 (skip if in range 0.98 to 1.02) + // if ((stats->robustStdev > 0.98) && (stats->robustStdev > 1.02)) { + // psLogMsg ("ppSub", PS_LOG_INFO, "background stdev (%f) is close to 1.0; do not modify variance\n", stats->robustStdev); + // // XXX set quality to poor + // psFree (stats); + // return true; + // } + + float varianceFactor = PS_SQR(stats->robustStdev); + psLogMsg ("ppSub", PS_LOG_INFO, "background stdev is not 1.0: multiplying variance by %f\n", varianceFactor); + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskBad) { + continue; + } + if (!isfinite(image->data.F32[y][x])) { + continue; + } + if (!isfinite(variance->data.F32[y][x]) || (variance->data.F32[y][x] < 0.0)) { + continue; + } + variance->data.F32[y][x] *= varianceFactor; + } + } + + pmFPA *outFPA = outRO->parent->parent->parent; + pmHDU *outHDU = outFPA->hdu; // Output HDU + + char history[80]; + snprintf (history, 80, "Rescaled variance by %6.4f (stdev by %6.4f)", varianceFactor, stats->robustStdev); + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history); + + psFree (stats); + return true; +}