Changeset 9336 for trunk/psModules/src/detrend/pmShutterCorrection.c

Timestamp:

Oct 5, 2006, 5:06:30 PM (20 years ago)

Author:

Paul Price

Message:

Adding iteration and masking to shutter measurement

File:

• trunk/psModules/src/detrend/pmShutterCorrection.c (modified) (12 diffs)

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -60 +60 @@
 #include <pslib.h>
 
+#include "pmFPA.h"
 #include "pmShutterCorrection.h"
 #include "psVectorBracket.h"
@@ -193 +194 @@
         const psVector *counts,
         const psVector *cntError,
-        float offref
+        const psVector *mask,
+        float offref,
+        int nIter,
+        float rej,
+        psMaskType maskVal
                                               )
 {
@@ -232 +237 @@
     line->coeff[1][1] = 0;
 
-    if (!psVectorFitPolynomial2D(line, NULL, 0, counts, cntError, x, y)) {
+    psStats *stats = psStatsAlloc(0);
+    stats->clipSigma = rej;
+    stats->clipIter = nIter;
+
+    if (!psVectorClipFitPolynomial2D(line, stats, mask, maskVal, counts, cntError, x, y)) {
         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit shutter correction.\n");
+        psFree(stats);
         psFree(x);
         psFree(y);
@@ -239 +249 @@
         return NULL;
     }
+    psFree(stats);
 
     pmShutterCorrection *corr = pmShutterCorrectionAlloc();
@@ -355 +366 @@
                                     psStatsOptions meanStat, // Statistic to use for mean
                                     psStatsOptions stdevStat, // Statistic to use for stdev
+                                    int nIter, // Number of iterations per pixel
+                                    float rej, // Rejection limit
                                     psMaskType maskVal // Mask value
                                    )
@@ -382 +395 @@
         if (!image) {
             continue;
+        }
+        if (image->type.type != PS_TYPE_F32) {
+            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Bad type for image: %x\n", image->type.type);
+            goto MEASURE_ERROR;
         }
         if (numRows == 0 || numCols == 0) {
@@ -401 +418 @@
             goto MEASURE_ERROR;
         }
+        if (masks) {
+            psImage *mask = masks->data[i];
+            if (mask) {
+                if (mask->type.type != PS_TYPE_U8) {
+                    psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Bad type for mask: %x\n", mask->type.type);
+                    goto MEASURE_ERROR;
+                }
+                if (mask->numRows != numRows || mask->numCols != numCols) {
+                    psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                            "Mask sizes don't match: %dx%d vs %dx%d\n", mask->numCols, mask->numRows,
+                            numCols, numRows);
+                    goto MEASURE_ERROR;
+                }
+            }
+        }
 
         // Measure statistics
@@ -407 +439 @@
         }
         refs->data.F32[i] = psStatsGetValue(stats, meanStat);
+        psTrace("psModules.detrend", 3, "Reference value for image %ld = %f\n", i, refs->data.F32[i]);
         if (refs->data.F32[i] <= 0.0) {
             psError(PS_ERR_UNKNOWN, true, "Measured non-positive reference value.\n");
@@ -424 +457 @@
             samplesMean->data.F32[j][i] = psStatsGetValue(stats, meanStat) * refs->data.F32[i];
             samplesStdev->data.F32[j][i] = psStatsGetValue(stats, stdevStat) * refs->data.F32[i];
+            psTrace("psModules.detrend", 5, "Image %ld, sample %d: %f +/- %f\n", i, j,
+                    samplesMean->data.F32[j][i], samplesStdev->data.F32[j][i]);
         }
     }
@@ -437 +472 @@
         errors = psImageRow(errors, samplesStdev, i);
         pmShutterCorrection *guess = pmShutterCorrectionGuess(exptimes, counts); // Guess at correction
-        pmShutterCorrection *params = pmShutterCorrectionFullFit(exptimes, counts, errors, guess); // Correc'n
-        if (isfinite(params->offref)) {
-            meanRef += params->offref;
+        pmShutterCorrection *corr = pmShutterCorrectionFullFit(exptimes, counts, errors, guess); // Correct'n
+        psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref);
+        if (isfinite(corr->offref)) {
+            meanRef += corr->offref;
             numGood++;
         }
-        psFree(params);
+        psFree(corr);
         psFree(guess);
     }
@@ -455 +491 @@
     }
     meanRef /= (float)numGood;
+    psTrace("psModules.detrend", 3, "Mean reference value: %f\n", meanRef);
 
     psImage *shutter = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Shutter correction image
     //psImage *pattern = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Illumination pattern
+    psVector *mask = psVectorAlloc(num, PS_TYPE_U8); // Mask for each image
+    mask->n = num;
+    psVectorInit(mask, 0);
+    psTrace("psModules.detrend", 2, "Performing linear fit on individual pixels...\n");
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
             for (int i = 0; i < num; i++) {
                 psImage *image = images->data[i]; // Image of interest
+                psImage *maskImage;     // Mask image
+                if (masks && (maskImage = masks->data[i])) {
+                    mask->data.U8[i] = maskImage->data.U8[y][x];
+                }
                 counts->data.F32[i] = image->data.F32[y][x] * refs->data.F32[i];
                 errors->data.F32[i] = sqrtf(image->data.F32[y][x]) * refs->data.F32[i];
             }
 
-            pmShutterCorrection *corr = pmShutterCorrectionLinFit(exptimes, counts, errors, meanRef);
+            pmShutterCorrection *corr = pmShutterCorrectionLinFit(exptimes, counts, errors, mask, meanRef,
+                                        nIter, rej, maskVal);
             shutter->data.F32[y][x] = corr->offset;
             //pattern->data.F32[y][x] = corr->scale;
@@ -488 +534 @@
     return NULL;
 }
+
+
+bool pmShutterCorrectionApply(pmReadout *readout, // Readout to which to apply shutter correction
+                              const pmReadout *shutter // Shutter correction readout
+                             )
+{
+    PS_ASSERT_PTR_NON_NULL(readout, false);
+    PS_ASSERT_PTR_NON_NULL(shutter, false);
+    PS_ASSERT_IMAGE_NON_NULL(readout->image, false);
+    PS_ASSERT_IMAGE_NON_NULL(shutter->image, false);
+    PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false);
+    PS_ASSERT_IMAGE_TYPE(shutter->image, PS_TYPE_F32, false);
+
+    psRegion region = psRegionSet(readout->col0, readout->col0 + readout->image->numCols,
+                                  readout->row0, readout->row0 + readout->image->numRows); // Detector region
+
+    pmCell *cell = readout->parent;     // Parent cell
+    if (!cell) {
+        psError(PS_ERR_BAD_PARAMETER_NULL, true,
+                "Parent cell is NULL --- unable to determine exposure time.\n");
+        return false;
+    }
+    float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time
+    if (!isfinite(exptime) || exptime < 0.0) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Bad exposure time: %f.\n", exptime);
+        return false;
+    }
+
+    psImage *shutterImage = psImageSubset(shutter->image, region); // Subimage with shutter
+    if (!shutterImage) {
+        psString regionString = psRegionToString(region);
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Size mismatch: %s vs %dx%d\n",
+                regionString, shutter->image->numCols, shutter->image->numRows);
+        psFree(regionString);
+        psFree(shutterImage);
+        return false;
+    }
+    psImage *image = readout->image;    // Image to correct
+    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);
+
+    return true;
+}