Changeset 21456 for trunk/psModules/src/detrend/pmFlatField.c

Timestamp:

Feb 12, 2009, 9:25:52 AM (17 years ago)

Author:

Paul Price

Message:

Adding modification of variance map by the flat-field. Removing multiple type cases for flat-field, since all images are of type F32.

File:

• trunk/psModules/src/detrend/pmFlatField.c (modified) (4 diffs)

trunk/psModules/src/detrend/pmFlatField.c

@@ -17 +17 @@
     PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
 
-    psImage *inImage   = job->args->data[0]; // Input image
-    psImage *inMask    = job->args->data[1]; // Input mask
-    const psImage *flatImage = job->args->data[2]; // Flat-field image
-    const psImage *flatMask  = job->args->data[3]; // Flat-field mask
+    psImage *inImage = job->args->data[0]; // Input image
+    psImage *inMask  = job->args->data[1]; // Input mask
+    psImage *inVar   = job->args->data[2]; // Input variance
+    const psImage *flatImage = job->args->data[3]; // Flat-field image
+    const psImage *flatMask  = job->args->data[4]; // Flat-field mask
 
-    psImageMaskType badFlat = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
-    int xOffset        = PS_SCALAR_VALUE(job->args->data[5],S32);
-    int yOffset        = PS_SCALAR_VALUE(job->args->data[6],S32);
-    int rowStart       = PS_SCALAR_VALUE(job->args->data[7],S32);
-    int rowStop        = PS_SCALAR_VALUE(job->args->data[8],S32);
-    return pmFlatFieldScan(inImage, inMask, flatImage, flatMask, badFlat,
+    psImageMaskType badFlat = PS_SCALAR_VALUE(job->args->data[5], PS_TYPE_IMAGE_MASK_DATA);
+    int xOffset        = PS_SCALAR_VALUE(job->args->data[6], S32);
+    int yOffset        = PS_SCALAR_VALUE(job->args->data[7], S32);
+    int rowStart       = PS_SCALAR_VALUE(job->args->data[8], S32);
+    int rowStop        = PS_SCALAR_VALUE(job->args->data[9], S32);
+
+    return pmFlatFieldScan(inImage, inMask, inVar, flatImage, flatMask, badFlat,
                            xOffset, yOffset, rowStart, rowStop);
 }
 
-// Macro for all PS types
-#define FLAT_DIVISION_CASE(TYPE, SPECIAL)       \
-    case PS_TYPE_##TYPE:                        \
-    for (int j = rowStart; j < rowStop; j++) { \
-        for (int i = 0; i < inImage->numCols; i++) { \
-            ps##TYPE flatValue = flatImage->data.TYPE[j + yOffset][i + xOffset]; \
-            if (!isfinite(flatValue) || flatValue <= 0.0 || \
-                (flatMask && flatMask->data.PS_TYPE_IMAGE_MASK_DATA[j + yOffset][i + xOffset])) { \
-                if (inMask) { \
-                    inMask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= badFlat; \
-                } \
-                inImage->data.TYPE[j][i] = SPECIAL; \
-            } else { \
-                inImage->data.TYPE[j][i] /= flatValue; \
-            } \
-        } \
-    } \
-    break;
 
-bool pmFlatFieldScan(psImage *inImage, psImage *inMask, const psImage *flatImage, const psImage *flatMask,
-                     psImageMaskType badFlat, int xOffset, int yOffset, int rowStart, int rowStop)
+bool pmFlatFieldScan(psImage *inImage, psImage *inMask, psImage *inVar, const psImage *flatImage,
+                     const psImage *flatMask, psImageMaskType badFlat,
+                     int xOffset, int yOffset, int rowStart, int rowStop)
 {
-    switch (inImage->type.type) {
-        FLAT_DIVISION_CASE(U8,  0);
-        FLAT_DIVISION_CASE(U16, 0);
-        FLAT_DIVISION_CASE(U32, 0);
-        FLAT_DIVISION_CASE(U64, 0);
-        FLAT_DIVISION_CASE(S8,  0);
-        FLAT_DIVISION_CASE(S16, 0);
-        FLAT_DIVISION_CASE(S32, 0);
-        FLAT_DIVISION_CASE(S64, 0);
-        FLAT_DIVISION_CASE(F32, NAN);
-        FLAT_DIVISION_CASE(F64, NAN);
-    default:
-        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unsupported type for input image: %x\n",
-                inImage->type.type);
-        return false;
+    // Neglecting asserts, because inputs should have been checked already
+
+    int numCols = inImage->numCols;     // Number of columns
+
+    // Using i,j for image; x,y for flat.
+    for (int j = rowStart, y = yOffset; j < rowStop; j++, y++) {
+        for (int i = 0, x = xOffset; i < numCols; i++, x++) {
+            float flatValue = 1.0 / flatImage->data.F32[y][x];
+            if (!isfinite(flatValue) || flatValue <= 0.0 ||
+                (flatMask && flatMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x])) {
+                if (inMask) {
+                    inMask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= badFlat;
+                }
+                inImage->data.F32[j][i] = NAN;
+                if (inVar) {
+                    inVar->data.F32[j][i] = NAN;
+                }
+            } else {
+                inImage->data.F32[j][i] *= flatValue;
+                if (inVar) {
+                    inVar->data.F32[j][i] *= PS_SQR(flatValue);
+                }
+            }
+        }
     }
+
     return true;
 }
@@ -74 +69 @@
 bool pmFlatField(pmReadout *in, const pmReadout *flat, psImageMaskType badFlat)
 {
-    PS_ASSERT_PTR_NON_NULL(in, false);
-    PS_ASSERT_PTR_NON_NULL(in->image, false);
-    PS_ASSERT_IMAGE_NON_EMPTY(in->image, false);
-    PS_ASSERT_PTR_NON_NULL(flat, false);
-    PS_ASSERT_PTR_NON_NULL(flat->image, false);
-    PS_ASSERT_IMAGE_NON_EMPTY(flat->image, false);
+    PM_ASSERT_READOUT_NON_NULL(in, false);
+    PM_ASSERT_READOUT_IMAGE(in, false);
+    PM_ASSERT_READOUT_NON_NULL(flat, false);
+    PM_ASSERT_READOUT_IMAGE(flat, false);
     if (in->mask) {
-        PS_ASSERT_IMAGE_TYPE(in->mask, PS_TYPE_IMAGE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(in->mask, in->image, false);
+        PM_ASSERT_READOUT_MASK(in, false);
     }
-    PS_ASSERT_IMAGE_TYPE(flat->image, in->image->type.type, false);
+    if (in->variance) {
+        PM_ASSERT_READOUT_VARIANCE(in, false);
+    }
     if (flat->mask) {
-        PS_ASSERT_IMAGE_TYPE(flat->mask, PS_TYPE_IMAGE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(flat->mask, flat->image, false);
+        PM_ASSERT_READOUT_MASK(flat, false);
     }
 
     psImage *inImage   = in->image;     // Input image
     psImage *inMask    = in->mask;      // Mask for input image
+    psImage *inVar     = in->variance;  // Variance for input image
     psImage *flatImage = flat->image;   // Flat-field image
     psImage *flatMask  = flat->mask;    // Mask for flat-field image
@@ -146 +140 @@
           psArrayAdd(job->args, 1, inImage);
           psArrayAdd(job->args, 1, inMask);
+          psArrayAdd(job->args, 1, inVar);
           psArrayAdd(job->args, 1, flatImage);
           psArrayAdd(job->args, 1, flatMask);
@@ -159 +154 @@
           }
           psFree(job);
-      } else if (!pmFlatFieldScan(inImage, inMask, flatImage, flatMask, badFlat,
+      } else if (!pmFlatFieldScan(inImage, inMask, inVar, flatImage, flatMask, badFlat,
                                   xOffset, yOffset, rowStart, rowStop)) {
           psError(PS_ERR_UNKNOWN, false, "Unable to flat-field image.");