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

Timestamp:

Mar 28, 2008, 5:10:17 PM (18 years ago)

Author:

Paul Price

Message:

Merging pap_branch_080328 so we can use the modernised version of ppMerge.

File:

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

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -78 +78 @@
     corr->offset = 0.0;
     corr->offref = 0.0;
+    corr->num = 0;
+    corr->stdev = NAN;
 
     return corr;
@@ -190 +192 @@
 
 // linear fit to the counts and exptime, given a value for offref
-pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime,
-        const psVector *counts,
-        const psVector *cntError,
-        const psVector *mask,
-        float offref,
-        int nIter,
-        float rej,
-        psMaskType maskVal
-                                              )
+pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, const psVector *counts,
+                                               const psVector *cntError, const psVector *mask, float offref,
+                                               int nIter, float rej, psMaskType maskVal)
 {
     PS_ASSERT_VECTOR_NON_NULL(exptime, NULL);
@@ -249 +245 @@
         return NULL;
     }
-    psFree(stats);
 
     pmShutterCorrection *corr = pmShutterCorrectionAlloc();
@@ -255 +250 @@
     corr->scale  = line->coeff[1][0];
     corr->offset = line->coeff[0][1] / line->coeff[1][0];
-
+    corr->num = stats->clippedNvalues;
+    corr->stdev = stats->clippedStdev;
+
+    psFree(stats);
     psFree(x);
     psFree(y);
@@ -263 +261 @@
 }
 
-static psF32 pmShutterCorrectionModel(psVector *deriv,
-                                      const psVector *params,
-                                      const psVector *x)
+static psF32 pmShutterCorrectionModel(psVector *deriv, const psVector *params, const psVector *x)
 {
     // This is in a tight loop, so we won't assert here.
@@ -283 +279 @@
 
 // non-linear fit to the counts and exptime, given a guess for the three parameters
-pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime,
-        const psVector *counts,
-        const psVector *cntError,
-        const pmShutterCorrection *guess
-                                               )
+pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, const psVector *counts,
+                                                const psVector *cntError, const pmShutterCorrection *guess)
 {
     PS_ASSERT_VECTOR_NON_NULL(exptime, NULL);
@@ -638 +631 @@
 
 
-bool pmShutterCorrectionApply(pmReadout *readout, const pmReadout *shutter)
+bool pmShutterCorrectionApply(pmReadout *readout, const pmReadout *shutter, psMaskType blank)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
@@ -682 +675 @@
     }
     psImage *image = readout->image;    // Image to correct
+    psImage *mask = readout->mask;      // Corresponding mask
 
     if (exptime <= 0.0) {
@@ -688 +682 @@
         for (int y = 0; y < image->numRows; y++) {
             for (int x = 0; x < image->numCols; x++) {
+                if (mask && !isfinite(shutterImage->data.F32[y][x])) {
+                    mask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
+                    image->data.F32[y][x] = NAN;
+                    continue;
+                }
                 image->data.F32[y][x] *= 1.0 / (exptime + shutterImage->data.F32[y][x]);
             }
@@ -699 +698 @@
         for (int y = 0; y < image->numRows; y++) {
             for (int x = 0; x < image->numCols; x++) {
+                if (mask && !isfinite(shutterImage->data.F32[y][x])) {
+                    mask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
+                    image->data.F32[y][x] = NAN;
+                    continue;
+                }
                 image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]);
             }
@@ -903 +907 @@
     PS_ASSERT_ARRAY_NON_NULL(inputs, false);
     PS_ASSERT_INT_EQUAL(data->num, inputs->n, false);
+    PS_ASSERT_INT_NONNEGATIVE(nIter, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
 
     int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
@@ -915 +921 @@
     if (pattern) {
         pmReadoutUpdateSize(pattern, minInputCols, minInputRows, xSize, ySize, false, false, maskVal);
+    }
+
+    psImage *nums = pmReadoutAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_COUNT, xSize, ySize,
+                                           PS_TYPE_U16, 0); // Image with number fitted per pixel
+    if (!nums) {
+        return false;
+    }
+    psImage *sigma = pmReadoutAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_SIGMA, xSize, ySize,
+                                            PS_TYPE_F32, NAN); // Image with stdev per pixel
+    if (!sigma) {
+        psFree(nums);
+        return false;
     }
 
@@ -955 +973 @@
             pmShutterCorrection *corr = pmShutterCorrectionLinFit(data->exptimes, counts, errors, mask,
                                                                   reference, nIter, rej, maskVal);
+            if (!corr) {
+                // Nothing we can do about it
+                psErrorClear();
+                shutterImage->data.F32[yOut][xOut] = NAN;
+                patternImage->data.F32[yOut][xOut] = NAN;
+                nums->data.U16[yOut][xOut] = 0;
+                sigma->data.F32[yOut][xOut] = NAN;
+                continue;
+            }
             shutterImage->data.F32[yOut][xOut] = corr->offset;
             patternImage->data.F32[yOut][xOut] = corr->scale;
+            nums->data.U16[yOut][xOut] = corr->num;
+            sigma->data.F32[yOut][xOut] = corr->stdev;
             psFree(corr);
         }
@@ -963 +992 @@
     psFree(errors);
     psFree(counts);
+    psFree(nums);
+    psFree(sigma);
 
     // Update the "concepts"