Changeset 12802 for trunk/psLib

Timestamp:

Apr 11, 2007, 9:54:47 AM (19 years ago)

Author:

Paul Price

Message:

It's not the flux contribution of bad pixels that matters for marking
a pixel as "poor" or "bad", but the kernel contribution relative to
the total kernel.

File:

• trunk/psLib/src/imageops/psImageInterpolate.c (modified) (7 diffs)

trunk/psLib/src/imageops/psImageInterpolate.c

@@ -7 +7 @@
  *  @author Paul Price, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-10 03:59:13 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-11 19:54:47 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -236 +236 @@
             for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
                 for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
-                    value += values[j][i] = kernel[j][i] * image->data.TYPE[yPix][xPix]; \
+                    *imageValue += kernel[j][i] * image->data.TYPE[yPix][xPix]; \
                 } \
             } \
@@ -243 +243 @@
 
     // Calculate the value for the image
-    double value = 0.0;                 // Value to return
-    double values[yNum][xNum];          // Values in image
-    switch (image->type.type) {
-        KERNEL_IMAGE_CASE(U8);
-        KERNEL_IMAGE_CASE(U16);
-        KERNEL_IMAGE_CASE(U32);
-        KERNEL_IMAGE_CASE(U64);
-        KERNEL_IMAGE_CASE(S8);
-        KERNEL_IMAGE_CASE(S16);
-        KERNEL_IMAGE_CASE(S32);
-        KERNEL_IMAGE_CASE(S64);
-        KERNEL_IMAGE_CASE(F32);
-        KERNEL_IMAGE_CASE(F64);
-      default:
-        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
-                image->type.type);
-        return false;
-    }
     if (imageValue) {
-        *imageValue = value;
+        *imageValue = 0.0;
+        switch (image->type.type) {
+            KERNEL_IMAGE_CASE(U8);
+            KERNEL_IMAGE_CASE(U16);
+            KERNEL_IMAGE_CASE(U32);
+            KERNEL_IMAGE_CASE(U64);
+            KERNEL_IMAGE_CASE(S8);
+            KERNEL_IMAGE_CASE(S16);
+            KERNEL_IMAGE_CASE(S32);
+            KERNEL_IMAGE_CASE(S64);
+            KERNEL_IMAGE_CASE(F32);
+            KERNEL_IMAGE_CASE(F64);
+          default:
+            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
+                    image->type.type);
+            return false;
+        }
     }
 
@@ -268 +266 @@
     const psImage *mask = options->mask; // Image mask
     if (maskValue && mask) {
+        *maskValue = 0;
         psMaskType maskVal = options->maskVal; // Mask value
-        double badValue = 0.0;          // Amount of flux in bad pixels
-        int badPix = 0;                 // Number of bad pixels
-        *maskValue = 0;
+        double badContrib = 0.0;        // Amount of kernel on bad pixels
+        double totContrib = 0.0;        // Total amount of kernel
+        bool badPix = false;            // Is there a bad pixel?
         for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) {
             for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) {
                 if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) {
-                    badValue += values[j][i];
-                    badPix++;
+                    badContrib += fabs(kernel[j][i]);
+                    badPix = true;
                 }
                 *maskValue |= mask->data.PS_TYPE_MASK_DATA[yPix][xPix];
+                totContrib += fabs(kernel[j][i]);
             }
         }
 
-        if (badPix > 0) {
-            if (fabs(badValue) >= fabs(value) * options->poorFrac) {
+        if (badPix) {
+            if (badContrib / totContrib >= options->poorFrac) {
                 *maskValue |= options->badMask;
             } else {
@@ -453 +453 @@
             double xInterpValue = 0.0; /* Interpolation in x */ \
             for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
-                xInterpValue += values[j][i] = xKernel[i] * image->data.TYPE[yPix][xPix]; \
+                xInterpValue += xKernel[i] * image->data.TYPE[yPix][xPix]; \
             } \
-            value += yKernel[j] * xInterpValue; /* Interpolating in y */ \
+            *imageValue += yKernel[j] * xInterpValue; /* Interpolating in y */ \
         } \
         break; \
@@ -461 +461 @@
 
     // Calculate the value for the image
-    double value = 0.0;                 // Value to return
-    double values[yNum][xNum];          // Values of interest
-    switch (image->type.type) {
-        SEPARATE_IMAGE_CASE(U8);
-        SEPARATE_IMAGE_CASE(U16);
-        SEPARATE_IMAGE_CASE(U32);
-        SEPARATE_IMAGE_CASE(U64);
-        SEPARATE_IMAGE_CASE(S8);
-        SEPARATE_IMAGE_CASE(S16);
-        SEPARATE_IMAGE_CASE(S32);
-        SEPARATE_IMAGE_CASE(S64);
-        SEPARATE_IMAGE_CASE(F32);
-        SEPARATE_IMAGE_CASE(F64);
-      default:
-        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
-                image->type.type);
-        return false;
-    }
     if (imageValue) {
-        *imageValue = value;
+        *imageValue = 0.0;
+        switch (image->type.type) {
+            SEPARATE_IMAGE_CASE(U8);
+            SEPARATE_IMAGE_CASE(U16);
+            SEPARATE_IMAGE_CASE(U32);
+            SEPARATE_IMAGE_CASE(U64);
+            SEPARATE_IMAGE_CASE(S8);
+            SEPARATE_IMAGE_CASE(S16);
+            SEPARATE_IMAGE_CASE(S32);
+            SEPARATE_IMAGE_CASE(S64);
+            SEPARATE_IMAGE_CASE(F32);
+            SEPARATE_IMAGE_CASE(F64);
+          default:
+            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
+                    image->type.type);
+            return false;
+        }
     }
 
@@ -487 +485 @@
     if (maskValue && mask) {
         psMaskType maskVal = options->maskVal; // Mask value
-        double badValue = 0.0;          // Amount of flux in bad pixels
-        int badPix = 0;                 // Number of bad pixels
+        double badContrib = 0.0;        // Amount of kernel on bad pixels
+        double totContrib = 0.0;        // Total amount of kernel
+        bool badPix = false;            // Number of bad pixels
         *maskValue = 0;
         for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) {
             // Interpolation in x
-            double xInterpValue = 0.0;
+            double xBadContrib = 0.0;   // Amount of x kernel on bad pixels
+            double xTotContrib = 0.0;   // Total amount of x kernel
             for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) {
                 if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) {
-                    xInterpValue += values[j][i];
-                    badPix++;
+                    xBadContrib += fabs(xKernel[i]);
+                    badPix = true;
                 }
                 *maskValue |= mask->data.PS_TYPE_MASK_DATA[yPix][xPix];
+                xTotContrib += fabs(xKernel[i]);
             }
             // Interpolating in y
-            badValue += yKernel[j] * xInterpValue;
-        }
-
-        if (badPix > 0) {
-            if (fabs(badValue) >= fabs(value) * options->poorFrac) {
+            badContrib += fabs(yKernel[j]) * xBadContrib;
+            totContrib += fabs(yKernel[j]) * xTotContrib;
+        }
+
+        if (badPix) {
+            if (badContrib / totContrib >= options->poorFrac) {
                 *maskValue |= options->badMask;
             } else {