Changeset 19140

Timestamp:

Aug 20, 2008, 3:12:44 PM (18 years ago)

Author:

Paul Price

Message:

Adding option of allowing 'exact' shifts when doing interpolation. This cuts down on a bit of the work in the event that the shift is integral. This can be turned off via a new boolean in psImageInterpolateOptions, so that a real interpolation kernel (not just with 0 and 1, but a real sinc kernel, for example) is returned. Hopefully this allows the use of psImageInterpolate to interpolate over bad pixels.

Location:

trunk/psLib/src/imageops

Files:

• psImageInterpolate.c (modified) (19 diffs)
• psImageInterpolate.h (modified) (2 diffs)

trunk/psLib/src/imageops/psImageInterpolate.c

@@ -7 +7 @@
  *  @author Paul Price, IfA
  *
- *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-08-20 04:06:31 $
+ *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-08-21 01:12:44 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -60 +60 @@
     options->poorMask = poorMask;
     options->poorFrac = poorFrac;
+    options->shifting = true;
 
     return options;
@@ -166 +167 @@
 static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel
                                              double kernel[xNum][yNum], // Kernel, to be set
+                                             bool *xExact, bool *yExact, // Exact shift?
                                              int xCentral, int yCentral, // Central pixel of convolution
                                              float x, float y, // Coordinates of interest
-                                             psImageInterpolateMode mode // Mode for interpolation
+                                             psImageInterpolateMode mode, // Mode for interpolation
+                                             bool allowExact // Allow exact shifts?
                                              )
 {
+    double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
+    double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
+    *xExact = (allowExact && fabs(xFrac) < DBL_EPSILON);
+    *yExact = (allowExact && fabs(yFrac) < DBL_EPSILON);
+
     switch (mode) {
       case PS_INTERPOLATE_BILINEAR: {
-          double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
-          double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
           kernel[0][0] = (1.0 - xFrac) * (1.0 - yFrac);
           kernel[0][1] = xFrac * (1.0 - yFrac);
@@ -220 +226 @@
         psAbort("Invalid interpolation mode.");
     }
-}
+    return;
+}
+
+// Can't interpolate: kernel is off the image
+#define INTERPOLATE_OFF() { \
+        *imageValue = options->badImage; \
+        if (varianceValue) { \
+            *varianceValue = options->badVariance; \
+        } \
+        if (maskValue) { \
+            *maskValue = options->badMask; \
+        } \
+        return PS_INTERPOLATE_STATUS_OFF; \
+    }
 
 // Set up and check interpolation bounds
@@ -231 +250 @@
     if (xStart > xLast || xStop < 0 || yStart > yLast || yStop < 0) { \
         /* Interpolation kernel is entirely off the image */ \
-        *imageValue = options->badImage; \
-        if (varianceValue) { \
-            *varianceValue = options->badVariance; \
-        } \
-        if (maskValue) { \
-            *maskValue = options->badMask; \
-        } \
-        return PS_INTERPOLATE_STATUS_OFF; \
+        INTERPOLATE_OFF(); \
     } \
     int xMin, xMax, yMin, yMax; /* Minimum and maximum valid pixels on image */ \
@@ -276 +288 @@
     }
 
+// Refine limits if the interpolation is exact
+#define INTERPOLATE_EXACT() { \
+    if (xExact) { \
+        if (iMin > 0 || iMax < 1) { \
+            INTERPOLATE_OFF(); /* Returns */ \
+        } \
+        iMin = (xNum - 1) / 2; \
+        iMax = iMin + 1; \
+    } \
+    if (yExact) { \
+        if (jMin > 0 || jMax < 1) { \
+            INTERPOLATE_OFF(); /* Returns */ \
+        } \
+        jMin = (yNum - 1) / 2; \
+        jMax = jMin + 1; \
+    } \
+    if (xExact && yExact) { \
+        return interpolateFlat(imageValue, varianceValue, maskValue, x, y, options); \
+    } \
+}
+
 // Interpolation engine using interpolation kernel
 static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue,
@@ -294 +327 @@
 
     double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
-    interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode);
+    bool xExact, yExact;                // Exact shifts?
+    interpolateKernelGenerate(xNum, yNum, kernel, &xExact, &yExact, xCentral, yCentral,
+                              x, y, options->mode, options->shifting);
 
     float sumKernel = 0.0;              // Sum of the kernel
@@ -306 +341 @@
     bool haveMask = mask && maskVal; // Does the user want the variance value?
 
+    INTERPOLATE_EXACT();
+
     // Add contributions in an area outside the image
 #define INTERPOLATE_KERNEL_ADD_OFFIMAGE() { \
@@ -316 +353 @@
     // Measure kernel contribution from outside image
     if (offImage) {
-        // Bottom rows
-        for (int j = 0; j < jMin; j++) {
-            for (int i = 0; i < xNum; i++) {
-                INTERPOLATE_KERNEL_ADD_OFFIMAGE();
+        if (!yExact) {
+            // Bottom rows
+            for (int j = 0; j < jMin; j++) {
+                for (int i = 0; i < xNum; i++) {
+                    INTERPOLATE_KERNEL_ADD_OFFIMAGE();
+                }
             }
         }
-        // Two sides
-        for (int j = jMin; j < jMax; j++) {
-            for (int i = 0; i < iMin; i++) {
-                INTERPOLATE_KERNEL_ADD_OFFIMAGE();
+        if (!xExact) {
+            // Two sides
+            for (int j = jMin; j < jMax; j++) {
+                for (int i = 0; i < iMin; i++) {
+                    INTERPOLATE_KERNEL_ADD_OFFIMAGE();
+                }
+                for (int i = iMax + 1; i < xNum; i++) {
+                    INTERPOLATE_KERNEL_ADD_OFFIMAGE();
+                }
             }
-            for (int i = iMax + 1; i < xNum; i++) {
-                INTERPOLATE_KERNEL_ADD_OFFIMAGE();
-            }
-        }
-        // Top rows
-        for (int j = jMax + 1; j < yNum; j++) {
-            for (int i = 0; i < xNum; i++) {
-                INTERPOLATE_KERNEL_ADD_OFFIMAGE();
+        }
+        if (!yExact) {
+            // Top rows
+            for (int j = jMax + 1; j < yNum; j++) {
+                for (int i = 0; i < xNum; i++) {
+                    INTERPOLATE_KERNEL_ADD_OFFIMAGE();
+                }
             }
         }
@@ -445 +488 @@
 // Generate Lanczos interpolation kernels
 static void lanczos(double values[],    // Interpolation kernel to generate
+                    bool *exact,        // Exact shift?
                     int num,            // Number of values in the kernel
-                    float frac          // Sub-pixel position
+                    float frac,         // Sub-pixel position
+                    bool allowExact     // Allow exact shift?
     )
 {
-    // XXX: Instead of generating a convolution kernel that does no shifting, need to set a boolean that says
-    // we can do an exact shift.
-    if (fabs(frac) < DBL_EPSILON) {
+    if (allowExact && fabs(frac) < DBL_EPSILON) {
+        *exact = true;
         // No real shift
         for (int i = 0; i < (num - 1) / 2; i++) {
@@ -461 +505 @@
         }
     } else {
-        double norm4 = 0.0;             // Normalisation
+        *exact = false;
         double norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos
         double norm2 = M_PI * 4.0 / (float)num; // Normalisation for sinc function 1
@@ -467 +511 @@
         double pos = - (num - 1)/2 - frac;  // Position of interest
         for (int i = 0; i < num; i++, pos += 1.0) {
-            norm4 += values[i] = norm1 * sin(pos * norm2) * sin(pos * norm3) / PS_SQR(pos);
-        }
-        norm4 = 1.0 / norm4;
-        for (int i = 0; i < num; i++) {
-            values[i] *= norm4;
+            if (fabs(pos) < DBL_EPSILON) {
+                values[i] = 1.0;
+            } else {
+                values[i] = norm1 * sin(pos * norm2) * sin(pos * norm3) / PS_SQR(pos);
+            }
         }
     }
@@ -510 +554 @@
 static inline void interpolateSeparateGenerate(int xNum, int yNum, // Size of interpolation kernel
                                                double xKernel[xNum], double yKernel[yNum], // Kernels
+                                               bool *xExact, bool *yExact, // Exact shifts?
                                                int xCentral, int yCentral, // Central pixel of convolution
                                                float x, float y, // Coordinates of interest
-                                               psImageInterpolateMode mode // Mode for interpolation
+                                               psImageInterpolateMode mode, // Mode for interpolation
+                                               bool allowExact // Allow an exact shift?
                                                )
 {
@@ -521 +567 @@
       case PS_INTERPOLATE_LANCZOS4: {
           double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
-          lanczos(xKernel, xNum, xFrac);
+          lanczos(xKernel, xExact, xNum, xFrac, allowExact);
           double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
-          lanczos(yKernel, yNum, yFrac);
+          lanczos(yKernel, yExact, yNum, yFrac, allowExact);
           break;
       }
@@ -553 +599 @@
 
     double xKernel[xNum], yKernel[yNum]; // Interpolation kernels
-    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode);
+    bool xExact, yExact;                // Exact shift?
+    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, &xExact, &yExact,
+                                xCentral, yCentral, x, y, options->mode, options->shifting);
 
     float sumKernel = 0.0;              // Sum of the kernel
@@ -564 +612 @@
     bool wantVariance = variance && varianceValue; // Does the user want the variance value?
     bool haveMask = mask && maskVal; // Does the user want the variance value?
+
+    INTERPOLATE_EXACT();
 
     // Add contributions in an area outside the image
@@ -585 +635 @@
     if (offImage) {
         // Bottom rows
-        for (int j = 0; j < jMin; j++) {
-            INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
-            for (int i = 0; i < xNum; i++) {
-                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
+        if (!yExact) {
+            for (int j = 0; j < jMin; j++) {
+                INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
+                for (int i = 0; i < xNum; i++) {
+                    INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
+                }
+                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
             }
-            INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
         }
         // Two sides
-        for (int j = jMin; j < jMax; j++) {
-            INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
-            for (int i = 0; i < iMin; i++) {
-                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
+        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 + 1; i < xNum; i++) {
+                    INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
+                }
+                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
             }
-            for (int i = iMax + 1; i < xNum; i++) {
-                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
+        }
+        // Top rows
+        if (!yExact) {
+            for (int j = jMax + 1; j < yNum; j++) {
+                if (!xExact) {
+                    INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
+                    for (int i = 0; i < xNum; i++) {
+                        INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
+                    }
+                    INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
+                }
             }
-            INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
-        }
-        // Top rows
-        for (int j = jMax + 1; j < yNum; j++) {
-            INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
-            for (int i = 0; i < xNum; i++) {
-                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
-            }
-            INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
         }
     }
@@ -826 +884 @@
 
     double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
-    interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode);
+    bool xExact, yExact;                // Exact shift?
+    // Pretend we're not shifting pixels (i.e., we don't care about any exact shifting) because on the off
+    // chance that a shift is integral, we don't want to get a trivial kernel back, but something that sort of
+    // reflects what's going on.
+    interpolateKernelGenerate(xNum, yNum, kernel, &xExact, &yExact, xCentral, yCentral,
+                              x, y, options->mode, false);
 
     double sumKernel2 = 0.0;            // Sum of kernel squares
@@ -855 +918 @@
 
     double xKernel[xNum], yKernel[yNum]; // Interpolation kernels for separable interpolation
-    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode);
+    bool xExact, yExact;                // Exact shift?
+    // Pretend we're not shifting pixels (i.e., we don't care about any exact shifting) because on the off
+    // chance that a shift is integral, we don't want to get a trivial kernel back, but something that sort of
+    // reflects what's going on.
+    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, &xExact, &yExact,
+                                xCentral, yCentral, x, y, options->mode, false);
 
     double ySumKernel2 = 0.0;           // Sum of kernel squared in y

trunk/psLib/src/imageops/psImageInterpolate.h

@@ -7 +7 @@
  * @author Paul Price, Institute for Astronomy
  *
- * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-06-17 21:24:58 $
+ * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2008-08-21 01:12:44 $
  * Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  */
@@ -53 +53 @@
     psMaskType poorMask;                ///< Mask value to give poor pixels
     float poorFrac;                     ///< Fraction of flux in bad pixels before output is marked bad
+    bool shifting;                      ///< Shifting images? Don't interpolate if the shift is exact.
 } psImageInterpolateOptions;