Add (bi-)cubic interpolation

skimage/_shared/interpolation.pxd

@@ -1,12 +1,19 @@
 
-cdef inline double nearest_neighbour(double* image, int rows, int cols,
-                                     double r, double c, char mode,
-                                     double cval)
+cdef inline double nearest_neighbour_interpolation(double* image, int rows,
+                                                   int cols, double r,
+                                                   double c, char mode,
+                                                   double cval)
 
 cdef inline double bilinear_interpolation(double* image, int rows, int cols,
                                           double r, double c, char mode,
                                           double cval)
 
+cdef inline double cubic_interpolation(double x, double[4] f)
+
+cdef inline double bicubic_interpolation(double* image, int rows, int cols,
+                                         double r, double c, char mode,
+                                         double cval)
+
 cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
                              char mode, double cval)
 

skimage/_shared/interpolation.pyx

@@ -5,16 +5,17 @@
 from libc.math cimport ceil, floor, round
 
 
-cdef inline double nearest_neighbour(double* image, int rows, int cols,
-                                     double r, double c, char mode,
-                                     double cval):
+cdef inline double nearest_neighbour_interpolation(double* image, int rows,
+                                                   int cols, double r,
+                                                   double c, char mode,
+                                                   double cval):
     """Nearest neighbour interpolation at a given position in the image.
 
     Parameters
     ----------
     image : double array
         Input image.
-    rows, cols: int
+    rows, cols : int
         Shape of image.
     r, c : int
         Position at which to interpolate.
@@ -23,6 +24,11 @@ cdef inline double nearest_neighbour(double* image, int rows, int cols,
     cval : double
         Constant value to use for constant mode.
 
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
     """
 
     return get_pixel(image, rows, cols, <int>round(r), <int>round(c),
@@ -38,7 +44,7 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols,
     ----------
     image : double array
         Input image.
-    rows, cols: int
+    rows, cols : int
         Shape of image.
     r, c : int
         Position at which to interpolate.
@@ -47,6 +53,11 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols,
     cval : double
         Constant value to use for constant mode.
 
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
     """
     cdef double dr, dc
     cdef int minr, minc, maxr, maxc
@@ -64,6 +75,79 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols,
     return (1 - dr) * top + dr * bottom
 
 
+cdef inline double cubic_interpolation(double x, double[4] f):
+    """Ccubic interpolation.
+
+    Parameters
+    ----------
+    x : double
+        Position in the interval [0, 1].
+    f : double[4]
+        Function values at positions [0, 1/3, 2/3, 1].
+
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
+    """
+    return \
+        f[1] + 0.5 * x * \
+            (f[2] - f[0] + x * \
+                (2.0 * f[0] - 5.0 * f[1] + 4.0 * f[2] - f[3] + x * \
+                    (3.0 * (f[1] - f[2]) + f[3] - f[0])))
+
+
+cdef inline double bicubic_interpolation(double* image, int rows, int cols,
+                                         double r, double c, char mode,
+                                         double cval):
+    """Bicubic interpolation at a given position in the image.
+
+    Parameters
+    ----------
+    image : double array
+        Input image.
+    rows, cols : int
+        Shape of image.
+    r, c : int
+        Position at which to interpolate.
+    mode : {'C', 'W', 'R'}
+        Wrapping mode. Constant, Wrap or Reflect.
+    cval : double
+        Constant value to use for constant mode.
+
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
+    """
+
+    cdef int r0 = <int>r
+    cdef int c0 = <int>c
+    if r < 0:
+        r0 -= 1
+    if c < 0:
+        c0 -= 1
+    # scale position to range [0, 1]
+    cdef double xr = (r - r0 + 1) / 3
+    cdef double xc = (c - c0 + 1) / 3
+
+    cdef double fc[4], fr[4]
+
+    cdef int pr, pc
+
+    for pr in range(r0 - 1, r0 + 3):
+
+        # do row-wise cubic interpolation
+        for pc in range(c0 - 1, c0 + 3):
+            fc[pc + 1 - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval)
+        fr[pr + 1 - r0] = cubic_interpolation(xc, fc)
+
+    # do cubic interpolation for interpolated values of each row
+    return cubic_interpolation(xr, fr)
+
+
 cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
                              char mode, double cval):
     """Get a pixel from the image, taking wrapping mode into consideration.
@@ -72,7 +156,7 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
     ----------
     image : double array
         Input image.
-    rows, cols: int
+    rows, cols : int
         Shape of image.
     r, c : int
         Position at which to get the pixel.
@@ -81,6 +165,11 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
     cval : double
         Constant value to use for constant mode.
 
+    Returns
+    -------
+    value : double
+        Pixel value at given position.
+
     """
     if mode == 'C':
         if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):

skimage/transform/_geometric.py

@@ -834,7 +834,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
     # use fast Cython version for specific parameters
     fast_modes = ('constant', 'reflect', 'wrap')
-    if order in (0, 1) and mode in fast_modes and not map_args:
+    if order in (0, 1, 3) and mode in fast_modes and not map_args:
         matrix = None
         if isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS):
             matrix = inverse_map._matrix

skimage/transform/_warps_cy.pyx

@@ -5,8 +5,9 @@
 
 cimport numpy as np
 import numpy as np
-from skimage._shared.interpolation cimport (nearest_neighbour,
-                                            bilinear_interpolation)
+from skimage._shared.interpolation cimport (nearest_neighbour_interpolation,
+                                            bilinear_interpolation,
+                                            bicubic_interpolation)
 
 
 cdef inline void _matrix_transform(double x, double y, double* H, double *x_,
@@ -108,9 +109,11 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
     cdef double (*interp_func)(double*, int, int, double, double,
                                char, double)
     if order == 0:
-        interp_func = nearest_neighbour
+        interp_func = nearest_neighbour_interpolation
     elif order == 1:
         interp_func = bilinear_interpolation
+    elif order == 3:
+        interp_func = bicubic_interpolation
 
     for tfr in range(out_r):
         for tfc in range(out_c):

skimage/transform/tests/test_warps.py

@@ -57,7 +57,7 @@ def test_fast_homography():
     tform = ProjectiveTransform(H)
     coords = warp_coords(tform.inverse, (img.shape[0], img.shape[1]))
 
-    for order in range(2):
+    for order in range(4):
         for mode in ('constant', 'reflect', 'wrap'):
             p0 = map_coordinates(img, coords, mode=mode, order=order)
             p1 = warp(img, tform, mode=mode, order=order)