Index: /trunk/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 13869) +++ /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 13870) @@ -13,4 +13,5 @@ #include "psVectorBracket.h" #include "pmConceptsAverage.h" +#include "pmReadoutStack.h" #include "pmShutterCorrection.h" @@ -58,4 +59,9 @@ +#define MEASURE_SAMPLES 4 // Number of samples to make over the image. This should only be + // changed with great caution, since assumptions on its value are in + // the code (see pmShutterCorrectionDataAlloc). + + static void pmShutterCorrectionFree(pmShutterCorrection *pars) { @@ -346,6 +352,4 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -#define MEASURE_SAMPLES 4 bool pmShutterCorrectionMeasure(pmReadout *output, const psArray *readouts, int size, psStatsOptions meanStat, @@ -680,22 +684,22 @@ if (exptime <= 0.0) { - // In the extreme case that we have exptime <= 0.0, we correct the image to - // counts-per-second, rather than counts in the nominal exposure time - for (int y = 0; y < image->numRows; y++) { - for (int x = 0; x < image->numCols; x++) { - image->data.F32[y][x] *= 1.0 / (exptime + shutterImage->data.F32[y][x]); - } - } - psMetadataAddF32 (cell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "exposure time re-normalized to 1.0", 1.0); // Exposure time - psString line = NULL; - psStringAppend (&line, "extreme exposure time %f, re-normalized to 1.0", exptime); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, line, ""); - psFree (line); + // In the extreme case that we have exptime <= 0.0, we correct the image to + // counts-per-second, rather than counts in the nominal exposure time + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + image->data.F32[y][x] *= 1.0 / (exptime + shutterImage->data.F32[y][x]); + } + } + psMetadataAddF32 (cell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "exposure time re-normalized to 1.0", 1.0); // Exposure time + psString line = NULL; + psStringAppend (&line, "extreme exposure time %f, re-normalized to 1.0", exptime); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, line, ""); + psFree (line); } else { - for (int y = 0; y < image->numRows; y++) { - for (int x = 0; x < image->numCols; x++) { - image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]); - } - } + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]); + } + } } psFree(shutterImage); @@ -712,2 +716,272 @@ return true; } + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + +#define IMAGES_BUFFER 10 // Allocate space for this many images at a time + +static void shutterCorrectionDataFree(pmShutterCorrectionData *data) +{ + psFree(data->regions); + psFree(data->mean); + psFree(data->stdev); + + psFree(data->exptimes); + psFree(data->refs); + + return; +} + +pmShutterCorrectionData *pmShutterCorrectionDataAlloc(int numCols, int numRows, int size) +{ + pmShutterCorrectionData *data = psAlloc(sizeof(pmShutterCorrectionData)); + psMemSetDeallocator(data, (psFreeFunc)shutterCorrectionDataFree); + + data->num = 0; + data->numCols = 0; + data->numRows = 0; + + data->regions = psArrayAlloc(MEASURE_SAMPLES); + for (int j = 0; j < MEASURE_SAMPLES; j++) { + int x = (j % 2) ? size : numCols - size - 1; + int y = (j > 1) ? size : numRows - size - 1; + psRegion region = psRegionForSquare(x, y, size); + data->regions->data[j] = psRegionAlloc(region.x0, region.x1, region.y0, region.y1); + } + + data->mean = psArrayAlloc(MEASURE_SAMPLES); + data->stdev = psArrayAlloc(MEASURE_SAMPLES); + for (int i = 0; i < MEASURE_SAMPLES; i++) { + data->mean->data[i] = psVectorAllocEmpty(IMAGES_BUFFER, PS_TYPE_F32); + data->stdev->data[i] = psVectorAllocEmpty(IMAGES_BUFFER, PS_TYPE_F32); + } + + data->exptimes = psVectorAllocEmpty(IMAGES_BUFFER, PS_TYPE_F32); + data->refs = psVectorAllocEmpty(IMAGES_BUFFER, PS_TYPE_F32); + + return data; +} + +bool pmShutterCorrectionAddReadout(pmShutterCorrectionData *data, + const pmReadout *readout, ///< Readout to add + psStatsOptions meanStat, ///< Statistic to use for mean + psStatsOptions stdevStat, ///< Statistic to use for stdev + psMaskType maskVal, ///< Mask value + psRandom *rng ///< Random number generator + ) +{ + PS_ASSERT_PTR_NON_NULL(data, NULL); + PS_ASSERT_PTR_NON_NULL(readout, NULL); + PS_ASSERT_IMAGE_NON_NULL(readout->image, NULL); + PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, NULL); + if (data->num == 0) { + data->numCols = readout->image->numCols; + data->numRows = readout->image->numRows; + } else { + PS_ASSERT_IMAGE_SIZE(readout->image, data->numCols, data->numRows, NULL); + } + if (readout->mask) { + PS_ASSERT_IMAGE_NON_NULL(readout->mask, NULL); + PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, NULL); + PS_ASSERT_IMAGE_SIZE(readout->mask, data->numCols, data->numRows, NULL); + } + if (readout->weight) { + PS_ASSERT_IMAGE_NON_NULL(readout->weight, NULL); + PS_ASSERT_IMAGE_TYPE(readout->weight, PS_TYPE_F32, NULL); + PS_ASSERT_IMAGE_SIZE(readout->weight, data->numCols, data->numRows, NULL); + } + + // Add the exposure time + float exptime = psMetadataLookupF32(NULL, readout->parent->concepts, "CELL.EXPOSURE"); // Exp. time + if (!isfinite(exptime)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure time is not set."); + return false; + } + data->exptimes->data.F32[data->exptimes->n] = exptime; + data->exptimes = psVectorExtend(data->exptimes, IMAGES_BUFFER, 1); + + // Add the statistics + + // Add the reference value + psStats *stats = psStatsAlloc(meanStat | stdevStat); // Statistics to apply + if (!rng) { + rng = psRandomAlloc(PS_RANDOM_TAUS, 0); + } else { + psMemIncrRefCounter(rng); + } + if (!psImageBackground(stats, readout->image, readout->mask, maskVal, rng)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure reference statistics.\n"); + psFree(stats); + psFree(rng); + return false; + } + psFree(rng); + float refValue = psStatsGetValue(stats, meanStat); // Reference value + psTrace("psModules.detrend", 3, "Reference value for shutter image = %f\n", refValue); + if (refValue <= 0.0) { + psError(PS_ERR_UNKNOWN, true, "Measured non-positive reference value.\n"); + psFree(stats); + return false; + } + refValue = 1.0 / refValue; + data->refs->data.F32[data->refs->n] = refValue; + data->refs = psVectorExtend(data->refs, IMAGES_BUFFER, 1); + + // Add the region statistics + for (int j = 0; j < MEASURE_SAMPLES; j++) { + psRegion *region = data->regions->data[j]; // Region of interest + psImage *subImage = psImageSubset(readout->image, *region); // Sub-image + psImage *subMask = NULL; // Sub-image of mask + if (readout->mask) { + subMask = psImageSubset(readout->mask, *region); + } + if (!psImageStats(stats, subImage, subMask, maskVal)) { + psString regionString = psRegionToString(*region); + psWarning("Unable to measure sample statistics at %s in image.\n", + regionString); + psFree(regionString); + } + psFree(subImage); + psFree(subMask); + + psVector *mean = data->mean->data[j]; // Vector of means for this region + psVector *stdev = data->stdev->data[j]; // Vector of standard deviations for this region + + mean->data.F32[mean->n] = psStatsGetValue(stats, meanStat) * refValue; + stdev->data.F32[stdev->n] = psStatsGetValue(stats, stdevStat) * refValue; + data->mean->data[j] = psVectorExtend(mean, IMAGES_BUFFER, 1); + data->stdev->data[j] = psVectorExtend(stdev, IMAGES_BUFFER, 1); + } + + data->num++; + + return true; +} + +float pmShutterCorrectionReference(const pmShutterCorrectionData *data) +{ + PS_ASSERT_PTR_NON_NULL(data, NAN); + PS_ASSERT_INT_POSITIVE(data->num, NAN); + + float meanRef = 0.0; // Mean reference offset + int numGood = 0; // Number of good measurements + for (int i = 0; i < MEASURE_SAMPLES; i++) { + psVector *counts = data->mean->data[i]; // Mean for each image + psVector *errors = data->stdev->data[i]; // Stdev for each image + pmShutterCorrection *guess = pmShutterCorrectionGuess(data->exptimes, counts); // Guess at correction + pmShutterCorrection *corr = pmShutterCorrectionFullFit(data->exptimes, counts, + errors, guess); // The actual correction + psFree(guess); + if (corr && isfinite(corr->offref)) { + psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref); + meanRef += corr->offref; + numGood++; + } + psFree(corr); + } + + if (numGood == 0) { + psError(PS_ERR_UNKNOWN, true, "Unable to measure mean reference offset.\n"); + return false; + } + meanRef /= (float)numGood; + psTrace("psModules.detrend", 3, "Mean reference value: %f\n", meanRef); + return meanRef; +} + +bool pmShutterCorrectionGenerate(pmReadout *shutter, pmReadout *pattern, const psArray *inputs, + float reference, const pmShutterCorrectionData *data, + int nIter, float rej, psMaskType maskVal) +{ + PS_ASSERT_PTR_NON_NULL(shutter, false); + PS_ASSERT_ARRAY_NON_NULL(inputs, false); + PS_ASSERT_INT_EQUAL(data->num, inputs->n, false); + + int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine + int xSize, ySize; // Size of the output image + if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, + inputs)) { + psError(PS_ERR_UNKNOWN, false, "No valid input readouts."); + return false; + } + + pmReadoutUpdateSize(shutter, minInputCols, minInputRows, xSize, ySize, false); + if (pattern) { + pmReadoutUpdateSize(pattern, minInputCols, minInputRows, xSize, ySize, false); + } + + psImage *shutterImage = shutter->image; // Shutter correction image + psImage *patternImage; // Illumination pattern + if (pattern) { + patternImage = psMemIncrRefCounter(pattern->image); + } else { + patternImage = psImageAlloc(xSize, ySize, PS_TYPE_F32); + } + + int num = data->num; // Number of images + psVector *counts = psVectorAlloc(num, PS_TYPE_F32); // Counts in each image + psVector *errors = psVectorAlloc(num, PS_TYPE_F32); // Counts in each image + psVector *mask = psVectorAlloc(num, PS_TYPE_MASK); // Mask for each image + psTrace("psModules.detrend", 2, "Performing linear fit on individual pixels...\n"); + for (int i = minInputRows; i < maxInputRows; i++) { + int yOut = i - shutter->row0; // y position on output readout + for (int j = minInputCols; j < maxInputCols; j++) { + int xOut = j - shutter->col0; // x position on output readout + + psVectorInit(mask, 0); + for (int r = 0; r < num; r++) { + pmReadout *readout = inputs->data[r]; // Readout of interest + int yIn = i - readout->row0; // y position on input readout + int xIn = j - readout->col0; // x position on input readout + psImage *image = readout->image; // Image of interest + float ref = data->refs->data.F32[r]; // (Inverse) reference value + counts->data.F32[r] = image->data.F32[yIn][xIn] * ref; + if (readout->mask) { + mask->data.PS_TYPE_MASK_DATA[r] = readout->mask->data.PS_TYPE_MASK_DATA[yIn][xIn]; + } + if (readout->weight) { + errors->data.F32[r] = sqrtf(readout->weight->data.F32[yIn][xIn]) * ref; + } else { + errors->data.F32[r] = sqrtf(image->data.F32[yIn][xIn]) * ref; + } + } + + pmShutterCorrection *corr = pmShutterCorrectionLinFit(data->exptimes, counts, errors, mask, + reference, nIter, rej, maskVal); + shutterImage->data.F32[yOut][xOut] = corr->offset; + patternImage->data.F32[yOut][xOut] = corr->scale; + psFree(corr); + } + } + psFree(mask); + psFree(errors); + psFree(counts); + + // Update the "concepts" + psList *inputCells = psListAlloc(NULL); // List of cells + for (long i = 0; i < inputs->n; i++) { + pmReadout *readout = inputs->data[i]; // Readout of interest + psListAdd(inputCells, PS_LIST_TAIL, readout->parent); + } + bool success = pmConceptsAverageCells(shutter->parent, inputCells, NULL, NULL, true); + psFree(inputCells); + + // Correct the exposure times --- they don't make sense any more. + psMetadataItem *item = psMetadataLookup(shutter->parent->concepts, "CELL.EXPOSURE"); + item->data.F32 = NAN; + item = psMetadataLookup(shutter->parent->concepts, "CELL.DARKTIME"); + item->data.F32 = NAN; + + shutter->data_exists = true; + shutter->parent->data_exists = true; + shutter->parent->parent->data_exists = true; + + if (pattern) { + pattern->data_exists = true; + pattern->parent->data_exists = true; + pattern->parent->parent->data_exists = true; + } + + return success; +} Index: /trunk/psModules/src/detrend/pmShutterCorrection.h =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.h (revision 13869) +++ /trunk/psModules/src/detrend/pmShutterCorrection.h (revision 13870) @@ -5,6 +5,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-04-24 21:17:19 $ + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-06-19 03:40:48 $ * Copyright 2006 Institute for Astronomy, University of Hawaii */ @@ -128,4 +128,59 @@ ); +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Functions for doing the shutter correction piece-meal (don't have to read entire image stack into memory at +// once). A single read run through the stack is required, calling pmShutterCorrectionAddReadout on each. +// Then pmShutterCorrectionReference provides the required reference shutter time, so that +// pmShutterCorrectionGenerate can generate a shutter correction piece by piece as overlapping pixels from +// each input are read in. + + +/// Data for measuring the shutter correction +typedef struct { + int num; ///< Number of images + int numCols, numRows; ///< Size of images + psArray *regions; ///< Regions at which to measure statistics + psArray *mean; ///< Vector of means at each region + psArray *stdev; ///< Vector of standard deviations at each region + psVector *exptimes; ///< Exposure times for each image + psVector *refs; ///< Reference fluxes +} pmShutterCorrectionData; + +/// Allocator for pmShutterCorrectionData +pmShutterCorrectionData *pmShutterCorrectionDataAlloc(int numCols, int numRows, ///< Size of images + int size ///< Size of regions + ); + +/// Add a readout to the correction data +/// +/// Performs statistics on the readout, recording the data +bool pmShutterCorrectionAddReadout(pmShutterCorrectionData *data, ///< Correction data + const pmReadout *readout, ///< Readout to add + psStatsOptions meanStat, ///< Statistic to use for mean + psStatsOptions stdevStat, ///< Statistic to use for stdev + psMaskType maskVal, ///< Mask value + psRandom *rng ///< Random number generator + ); + +/// Calculate the reference shutter time from the correction data +float pmShutterCorrectionReference(const pmShutterCorrectionData *data ///< Correction data + ); + +/// Generate a shutter correction +/// +/// Performs the linear fit to each pixel in the stack. +bool pmShutterCorrectionGenerate(pmReadout *shutter, ///< Shutter correction + pmReadout *pattern, ///< Background pattern (or NULL) + const psArray *inputs, ///< Stack of input pmReadouts + float reference, ///< Reference shutter time (from pmShutterCorrectionRef) + const pmShutterCorrectionData *data, ///< Correction data + int nIter, ///< Number of iterations + float rej, ///< Rejection threshold (sigma) + psMaskType maskVal ///< Mask value + ); + + + /// @} #endif