Changeset 26836 for branches/eam_branches/20091201/psModules/src/detrend/pmPattern.c

Timestamp:

Feb 9, 2010, 8:18:37 PM (16 years ago)

Author:

Paul Price

Message:

Adding I/O for pattern correction (both row and cell, in the same parts).

File:

• branches/eam_branches/20091201/psModules/src/detrend/pmPattern.c (modified) (11 diffs)

branches/eam_branches/20091201/psModules/src/detrend/pmPattern.c

@@ -1 +1 @@
 #include <stdio.h>
+#include <string.h>
 #include <pslib.h>
 
@@ -9 +10 @@
 // Mask a row as bad
 static void patternMaskRow(pmReadout *ro, // Readout to mask
-                           int y,         // Row to mask
+                           int y,       // Row to mask
                            psImageMaskType bad // Mask value to give
                            )
@@ -17 +18 @@
     psAssert(y < image->numRows, "Row not in image");
 
-    int numCols = image->numCols;       // Size of image
+    int numCols = image->numCols;       // Number of columns
     for (int x = 0; x < numCols; x++) {
         image->data.F32[y][x] = NAN;
@@ -29 +30 @@
     return;
 }
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Measurement and application
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
 bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh,
@@ -80 +85 @@
     psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, order); // Polynomial to fit
     psVector *data = psVectorAlloc(numCols, PS_TYPE_F32); // Data to fit
+    psArray *corrs = psArrayAlloc(numRows); // Correction parameters
+
+    psImage *corr = psImageAlloc(order + 1, numRows, PS_TYPE_F64); // Corrections applied
 
     for (int y = 0; y < numRows; y++) {
@@ -105 +113 @@
             continue;
         }
+
+        memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        corr->data.F64[y][0] -= background;
+
         psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector
         if (!solution) {
@@ -118 +130 @@
         psFree(solution);
     }
+
+    psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE,
+                       "Pattern row correction", corr);
+    psFree(corrs);
 
     psFree(indices);
@@ -128 +144 @@
 }
 
-
+bool pmPatternRowApply(pmReadout *ro, psImageMaskType maskBad)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PM_ASSERT_READOUT_IMAGE(ro, false);
+
+    bool mdok;                          // Status of MD lookup
+    psImage *corr = psMetadataLookupPtr(&mdok, ro->analysis, PM_PATTERN_ROW_CORRECTION); // Correction
+    if (!mdok) {
+        // No correction to apply
+        return true;
+    }
+
+    psImage *image = ro->image; // Image of interest
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, numCols - 1); // Polynomial to apply
+    psVector *indices = psVectorAlloc(numCols, PS_TYPE_F32); // Indices for polynomial
+    float norm = 2.0 / (float)numCols;  // Normalisation for indices
+    for (int x = 0; x < numCols; x++) {
+        indices->data.F32[x] = x * norm - 1.0;
+    }
+
+    for (int y = 0; y < numRows; y++) {
+        memcpy(poly->coeff, corr->data.F64[y], numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector
+        if (!solution) {
+            psWarning("Unable to evaluate polynomial for row %d", y);
+            psErrorClear();
+            patternMaskRow(ro, y, maskBad);
+            continue;
+        }
+
+        for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] -= solution->data.F32[x];
+        }
+        psFree(solution);
+    }
+
+    psFree(poly);
+    psFree(indices);
+
+    return true;
+}
 
 
@@ -255 +313 @@
             psImageInit(ro->image, NAN);
             psBinaryOp(ro->mask, ro->mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK));
+            psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE,
+                             "Pattern cell correction solution", NAN);
             continue;
         }
@@ -262 +322 @@
                  cellName, correction);
         psBinaryOp(ro->image, ro->image, "+", psScalarAlloc(correction, PS_TYPE_F32));
+        psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE,
+                         "Pattern cell correction solution", correction);
     }
 
@@ -269 +331 @@
     return true;
 }
+
+bool pmPatternCellApply(pmReadout *ro, psImageMaskType maskBad)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PM_ASSERT_READOUT_IMAGE(ro, false);
+
+    bool mdok;                          // Status of MD lookup
+    float corr = psMetadataLookupF32(&mdok, ro->analysis, PM_PATTERN_CELL_CORRECTION); // Correction to apply
+    if (!mdok) {
+        // No correction to apply
+        return true;
+    }
+
+    psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    if (!isfinite(corr)) {
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                image->data.F32[y][x] = NAN;
+            }
+        }
+        if (mask) {
+            for (int y = 0; y < numRows; y++) {
+                for (int x = 0; x < numCols; x++) {
+                    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad;
+                }
+            }
+        }
+    } else {
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                image->data.F32[y][x] += corr;
+            }
+        }
+    }
+
+    return true;
+}
+
+