Changeset 13870

Timestamp:

Jun 18, 2007, 5:40:48 PM (19 years ago)

Author:

Paul Price

Message:

Splitting out shutter correction so that it can run as we read in bits and pieces.

Location:

trunk/psModules/src/detrend

Files:

• pmShutterCorrection.c (modified) (5 diffs)
• pmShutterCorrection.h (modified) (2 diffs)

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -13 +13 @@
 #include "psVectorBracket.h"
 #include "pmConceptsAverage.h"
+#include "pmReadoutStack.h"
 
 #include "pmShutterCorrection.h"
@@ -58 +59 @@
 
 
+#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 +352 @@
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-#define MEASURE_SAMPLES 4
 
 bool pmShutterCorrectionMeasure(pmReadout *output, const psArray *readouts, int size, psStatsOptions meanStat,
@@ -680 +684 @@
 
     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 +716 @@
     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;
+}

trunk/psModules/src/detrend/pmShutterCorrection.h

@@ -5 +5 @@
  * @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 +128 @@
                              );
 
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// 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