From 3b227e226d053685a5cc55b6d6d1f9e48cd826f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 21 Aug 2012 09:20:59 +0200 Subject: [PATCH] Add nearest neighbour interpolation --- skimage/_shared/interpolation.pxd | 4 ++++ skimage/_shared/interpolation.pyx | 26 +++++++++++++++++++++++++- skimage/transform/_project.pyx | 26 +++++++++++++++++++------- skimage/transform/tests/test_warps.py | 27 +++++++++++++++------------ 4 files changed, 63 insertions(+), 20 deletions(-) diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd index 7ad5d223..c883d00d 100644 --- a/skimage/_shared/interpolation.pxd +++ b/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=*) diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index defa2e9b..71852ace 100644 --- a/skimage/_shared/interpolation.pyx +++ b/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, round(r), round(c), + mode, cval) cdef inline double bilinear_interpolation(double* image, int rows, int cols, diff --git a/skimage/transform/_project.pyx b/skimage/transform/_project.pyx index 300ba543..4f9a4704 100644 --- a/skimage/transform/_project.pyx +++ b/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, M.data, &c, &r) - out[tfr, tfc] = bilinear_interpolation(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, M.data, &c, &r) + out[tfr, tfc] = nearest_neighbour(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, M.data, &c, &r) + out[tfr, tfc] = bilinear_interpolation(img.data, rows, + cols, r, c, mode_c) return out diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index f35b7f57..8c0d81d3 100644 --- a/skimage/transform/tests/test_warps.py +++ b/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():