Index: trunk/ppSub/src/Makefile.am =================================================================== --- trunk/ppSub/src/Makefile.am (revision 24334) +++ trunk/ppSub/src/Makefile.am (revision 24672) @@ -30,4 +30,5 @@ ppSubVersion.c \ ppSubBackground.c \ + ppSubVarianceRescale.c \ ppSubCamera.c \ ppSubData.c \ Index: trunk/ppSub/src/ppSub.h =================================================================== --- trunk/ppSub/src/ppSub.h (revision 24334) +++ trunk/ppSub/src/ppSub.h (revision 24672) @@ -101,4 +101,8 @@ ); +/// Perform Variance correction (rescale within a modest range) +bool ppSubVarianceRescale(pmConfig *config ///< Configuration + ); + /// Put the program version information into a header bool ppSubVersionHeader(psMetadata *header ///< Header to populate Index: trunk/ppSub/src/ppSubLoop.c =================================================================== --- trunk/ppSub/src/ppSubLoop.c (revision 24334) +++ trunk/ppSub/src/ppSubLoop.c (revision 24672) @@ -98,4 +98,10 @@ } + // Perform Variance correction (rescale within a modest range) + if (!ppSubVarianceRescale(config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to rescale variance."); + return false; + } + if (!ppSubFilesIterateDown(config, PPSUB_FILES_PHOT_SUB)) { psError(PPSUB_ERR_IO, false, "Unable to set up photometry files."); Index: trunk/ppSub/src/ppSubVarianceRescale.c =================================================================== --- trunk/ppSub/src/ppSubVarianceRescale.c (revision 24672) +++ trunk/ppSub/src/ppSubVarianceRescale.c (revision 24672) @@ -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; +}