Changeset 20325 for trunk/psLib/src/imageops/psImageInterpolate.c

Timestamp:

Oct 22, 2008, 8:40:36 AM (18 years ago)

Author:

Paul Price

Message:

Fix interpolation with separable kernels:

• sumKernel2 is now for the entire kernel, not just the good pixels.
• was using the wrong values of the kernel when near the edge
• was using the wrong kernelValue in one case
• can optimise bad contributions since kernel^{2 sums are known}

File:

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

trunk/psLib/src/imageops/psImageInterpolate.c

@@ -7 +7 @@
  *  @author Paul Price, IfA
  *
- *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-10-22 03:53:10 $
+ *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-10-22 18:40:36 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -429 +429 @@
     *imageValue = sumImage / sumKernel; \
     if (wantVariance) { \
-        *varianceValue = sumVariance / sumKernel2; \
+        *varianceValue = sumVariance / (sumKernel2 - sumBad); \
     } \
     if (sumBad == 0) { \
         /* Completely good pixel */ \
         status = PS_INTERPOLATE_STATUS_GOOD; \
-    } else if (sumBad < PS_SQR(interp->poorFrac) * (sumBad + sumKernel2) ) { \
+    } else if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) { \
         /* Some pixels masked: poor pixel */ \
         if (haveMask && maskValue) { \
@@ -535 +535 @@
 
     if (offImage) {
-    // Add contributions in an area outside the image
-#define INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL() \
-    float xSumBad = 0.0;
-#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL() { \
-        xSumBad += xKernel2[i]; \
-    }
-#define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW() { \
-        sumBad += yKernel2[j] * xSumBad; \
-    }
-
         // Bottom rows
-        if (!yExact) {
-            for (int j = 0; j < jMin; j++) {
-                INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
-                for (int i = 0; i < size; i++) {
-                    INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
-                }
-                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
+        for (int j = 0; j < jMin; j++) {
+            sumBad += yKernel2[j] * xSumKernel2;
+        }
+        // Two sides
+        for (int j = jMin; j < jMax; j++) {
+            float xSumBad = 0.0;
+            for (int i = 0; i < iMin; i++) {
+                xSumBad += xKernel2[i];
             }
-        }
-        // Two sides
-        if (!xExact) {
-            for (int j = jMin; j < jMax; j++) {
-                INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
-                for (int i = 0; i < iMin; i++) {
-                    INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
-                }
-                for (int i = iMax; i < size; i++) {
-                    INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
-                }
-                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
+            for (int i = iMax; i < size; i++) {
+                xSumBad += xKernel2[i];
             }
+            sumBad += ySumKernel2 * xSumBad;
         }
         // Top rows
-        if (!yExact) {
-            for (int j = jMax; j < size; j++) {
-                if (!xExact) {
-                    INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
-                    for (int i = 0; i < size; i++) {
-                        INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
-                    }
-                    INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
-                }
-            }
-        }
-    }
+        for (int j = jMax; j < size; j++) {
+            sumBad += yKernel2[j] * xSumKernel2;
+        }
+    }
+
+    yKernel += jMin;
+    yKernel2 += jMin;
+    xKernel += iMin;
+    xKernel2 += iMin;
 
 #define INTERPOLATE_SEPARATE_CASE(TYPE) \
@@ -600 +579 @@
                   const psF32 *xKernelData = xKernel; \
                   const psF32 *xKernel2Data = xKernel2; \
-                  for (int i = iMin, xPix = xMin; i < iMax; \
-                           i++, xPix++, imageData++, varianceData++, maskData++, \
-                           xKernelData++, xKernel2Data++) { \
+                  for (int i = iMin; i < iMax; i++, imageData++, varianceData++, maskData++, \
+                                               xKernelData++, xKernel2Data++) { \
                       float kernelValue = *xKernelData; /* Value of kernel in x */ \
                       float kernelValue2 = *xKernel2Data; /* Square of kernel in x */ \
@@ -613 +591 @@
                       } \
                   } \
-                  float kernelValue = *yKernel; /* Value of kernel in y */ \
+                  float kernelValue = *yKernelData; /* Value of kernel in y */ \
                   float kernelValue2 = *yKernel2Data; /* Value of kernel-squared in y */ \
                   sumImage += kernelValue * xSumImage; \
@@ -633 +611 @@
                   const psF32 *xKernelData = xKernel; \
                   const psF32 *xKernel2Data = xKernel2; \
-                  for (int i = iMin, xPix = xMin; i < iMax; \
-                           i++, xPix++, imageData++, varianceData++, xKernelData++, xKernel2Data++) { \
+                  for (int i = iMin; i < iMax; i++, imageData++, varianceData++, \
+                                               xKernelData++, xKernel2Data++) { \
                       float kernelValue = *xKernelData; /* Value of kernel */ \
                       float kernelValue2 = *xKernel2Data; /* Square of kernel */ \
@@ -661 +639 @@
               const psF32 *xKernelData = xKernel; \
               const psF32 *xKernel2Data = xKernel2; \
-              for (int i = iMin, xPix = xMin; i < iMax; \
-                       i++, xPix++, imageData++, maskData++, xKernelData++, xKernel2Data++) { \
+              for (int i = iMin; i < iMax; i++, imageData++, maskData++, xKernelData++, xKernel2Data++) { \
                   float kernelValue = *xKernelData; /* Value of kernel */ \
                   float kernelValue2 = *xKernel2Data; /* Value of kernel-squared */ \
@@ -687 +664 @@
               const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \
               const psF32 *xKernelData = xKernel; \
-              for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++, imageData++, xKernelData++) { \
+              for (int i = iMin; i < iMax; i++, imageData++, xKernelData++) { \
                   float kernelValue = *xKernelData; /* Value of kernel */ \
                   xSumImage += kernelValue * *imageData; \