Add nearest neighbour interpolation

skimage/_shared/interpolation.pxd

@@ -1,4 +1,8 @@
 
+cdef inline double nearest_neighbour(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=*)

skimage/_shared/interpolation.pyx

@@ -2,7 +2,31 @@
 #cython: boundscheck=False
 #cython: nonecheck=False
 #cython: wraparound=False
-from libc.math cimport ceil, floor
+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=0):
+    """Nearest neighbour 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', 'M'}
+        Wrapping mode. Constant, Wrap or Mirror.
+    cval : double
+        Constant value to use for constant mode.
+
+    """
+
+    return get_pixel(image, rows, cols, <int>round(r), <int>round(c),
+                     mode, cval)
 
 
 cdef inline double bilinear_interpolation(double* image, int rows, int cols,

skimage/transform/_project.pyx

@@ -5,7 +5,8 @@
 
 cimport numpy as np
 import numpy as np
-from skimage._shared.interpolation cimport bilinear_interpolation
+from skimage._shared.interpolation cimport (nearest_neighbour,
+                                            bilinear_interpolation)
 
 
 cdef inline _matrix_transform(double x, double y, double* H, double *x_,
@@ -32,7 +33,7 @@ cdef inline _matrix_transform(double x, double y, double* H, double *x_,
     y_[0] = yy / zz
 
 
-def homography(np.ndarray image, np.ndarray H, output_shape=None,
+def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
                mode='constant', double cval=0):
     """
     Projective transformation (homography).
@@ -67,6 +68,10 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
         Transformation matrix H that defines the homography.
     output_shape : tuple (rows, cols)
         Shape of the output image generated.
+    order : {0, 1}
+        Order of interpolation::
+        * 0: Nearest-neighbour interpolation.
+        * 1: Bilinear interpolation (default).
     mode : {'constant', 'mirror', 'wrap'}
         How to handle values outside the image borders.
     cval : string
@@ -104,10 +109,17 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
     cdef int rows = img.shape[0]
     cdef int cols = img.shape[1]
 
-    for tfr in range(out_r):
-        for tfc in range(out_c):
-            _matrix_transform(tfc, tfr, <double*>M.data, &c, &r)
-            out[tfr, tfc] = bilinear_interpolation(<double*>img.data, rows,
-                                                   cols, r, c, mode_c)
+    if order == 0:
+        for tfr in range(out_r):
+            for tfc in range(out_c):
+                _matrix_transform(tfc, tfr, <double*>M.data, &c, &r)
+                out[tfr, tfc] = nearest_neighbour(<double*>img.data, rows,
+                                                  cols, r, c, mode_c)
+    elif order == 1:
+        for tfr in range(out_r):
+            for tfc in range(out_c):
+                _matrix_transform(tfc, tfr, <double*>M.data, &c, &r)
+                out[tfr, tfc] = bilinear_interpolation(<double*>img.data, rows,
+                                                       cols, r, c, mode_c)
 
     return out

skimage/transform/tests/test_warps.py

@@ -48,20 +48,23 @@ def test_fast_homography():
     H[:2, :2] = [[C, -S], [S, C]]
     H[:2, 2] = [tx, ty]
 
-    for mode in ('constant', 'mirror', 'wrap'):
-        p0 = warp(img, ProjectiveTransform(H).inverse, mode=mode, order=1)
-        p1 = fast_homography(img, H, mode=mode)
+    tform = ProjectiveTransform(H)
 
-        # import matplotlib.pyplot as plt
-        # f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
-        # ax0.imshow(img)
-        # ax1.imshow(p0, cmap=plt.cm.gray)
-        # ax2.imshow(p1, cmap=plt.cm.gray)
-        # ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
-        # plt.show()
+    for order in range(2):
+        for mode in ('constant', 'mirror', 'wrap'):
+            p0 = warp(img, tform.inverse, mode=mode, order=order)
+            p1 = fast_homography(img, H, mode=mode, order=order)
 
-        d = np.mean(np.abs(p0 - p1))
-        assert d < 0.001
+            # import matplotlib.pyplot as plt
+            # f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
+            # ax0.imshow(img)
+            # ax1.imshow(p0, cmap=plt.cm.gray)
+            # ax2.imshow(p1, cmap=plt.cm.gray)
+            # ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
+            # plt.show()
+
+            d = np.mean(np.abs(p0 - p1))
+            assert d < 0.001
 
 
 def test_swirl():