diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd index 7d212969..088be1ef 100644 --- a/skimage/_shared/interpolation.pxd +++ b/skimage/_shared/interpolation.pxd @@ -9,7 +9,8 @@ cdef inline Py_ssize_t round(double r): return ((r + 0.5) if (r > 0.0) else (r - 0.5)) -cdef inline double nearest_neighbour_interpolation(double* image, Py_ssize_t rows, +cdef inline double nearest_neighbour_interpolation(double* image, + Py_ssize_t rows, Py_ssize_t cols, double r, double c, char mode, double cval): @@ -63,30 +64,33 @@ cdef inline double bilinear_interpolation(double* image, Py_ssize_t rows, """ cdef double dr, dc - cdef Py_ssize_t minr, minc, maxr, maxc + cdef long minr, minc, maxr, maxc - minr = floor(r) - minc = floor(c) - maxr = ceil(r) - maxc = ceil(c) + minr = floor(r) + minc = floor(c) + maxr = ceil(r) + maxc = ceil(c) dr = r - minr dc = c - minc top = (1 - dc) * get_pixel2d(image, rows, cols, minr, minc, mode, cval) \ + dc * get_pixel2d(image, rows, cols, minr, maxc, mode, cval) - bottom = (1 - dc) * get_pixel2d(image, rows, cols, maxr, minc, mode, cval) \ + bottom = (1 - dc) * get_pixel2d(image, rows, cols, maxr, minc, mode, + cval) \ + dc * get_pixel2d(image, rows, cols, maxr, maxc, mode, cval) return (1 - dr) * top + dr * bottom cdef inline double quadratic_interpolation(double x, double[3] f): - """Quadratic interpolation. + """WARNING: Do not use, not implemented correctly. + + Quadratic interpolation. Parameters ---------- x : double - Position in the interval [-1, 1]. - f : double[4] - Function values at positions [-1, 0, 1]. + Position in the interval [0, 2]. + f : double[3] + Function values at positions [0, 2]. Returns ------- @@ -94,13 +98,17 @@ cdef inline double quadratic_interpolation(double x, double[3] f): Interpolated value. """ - return f[1] - 0.25 * (f[0] - f[2]) * x + return (x * f[2] * (x - 1)) / 2 - \ + x * f[1] * (x - 2) + \ + (f[0] * (x - 1) * (x - 2)) / 2 cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows, - Py_ssize_t cols, double r, double c, - char mode, double cval): - """Biquadratic interpolation at a given position in the image. + Py_ssize_t cols, double r, + double c, char mode, double cval): + """WARNING: Do not use, not implemented correctly. + + Biquadratic interpolation at a given position in the image. Parameters ---------- @@ -122,30 +130,23 @@ cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows, """ - cdef Py_ssize_t r0 = round(r) - cdef Py_ssize_t c0 = round(c) - if r < 0: - r0 -= 1 - if c < 0: - c0 -= 1 - # scale position to range [-1, 1] - cdef double xr = (r - r0) - 1 - cdef double xc = (c - c0) - 1 - if r == r0: - xr += 1 - if c == c0: - xc += 1 + cdef long r0 = round(r) - 1 + cdef long c0 = round(c) - 1 + + cdef double xr = r - r0 + cdef double xc = c - c0 cdef double fc[3] cdef double fr[3] - cdef Py_ssize_t pr, pc + cdef long pr, pc # row-wise cubic interpolation - for pr in range(r0, r0 + 3): - for pc in range(c0, c0 + 3): - fc[pc - c0] = get_pixel2d(image, rows, cols, pr, pc, mode, cval) - fr[pr - r0] = quadratic_interpolation(xc, fc) + for pr in range(3): + for pc in range(3): + fc[pc] = get_pixel2d(image, rows, cols, + r0 + pr, c0 + pc, mode, cval) + fr[pr] = quadratic_interpolation(xc, fc) # cubic interpolation for interpolated values of each row return quadratic_interpolation(xr, fr) @@ -159,7 +160,7 @@ cdef inline double cubic_interpolation(double x, double[4] f): x : double Position in the interval [0, 1]. f : double[4] - Function values at positions [0, 1/3, 2/3, 1]. + Function values at positions [-1, 0, 1, 2]. Returns ------- @@ -179,6 +180,9 @@ cdef inline double bicubic_interpolation(double* image, Py_ssize_t rows, char mode, double cval): """Bicubic interpolation at a given position in the image. + Interpolation using Catmull-Rom splines, based on the bicubic convolution + algorithm described in [1]_. + Parameters ---------- image : double array @@ -197,35 +201,42 @@ cdef inline double bicubic_interpolation(double* image, Py_ssize_t rows, value : double Interpolated value. + References + ---------- + .. [1] R. Keys, (1981). "Cubic convolution interpolation for digital image + processing". IEEE Transactions on Signal Processing, Acoustics, + Speech, and Signal Processing 29 (6): 1153–1160. + """ - cdef Py_ssize_t r0 = r - 1 - cdef Py_ssize_t c0 = c - 1 - if r < 0: - r0 -= 1 - if c < 0: - c0 -= 1 + cdef long r0 = floor(r) + cdef long c0 = floor(c) + # scale position to range [0, 1] - cdef double xr = (r - r0) / 3 - cdef double xc = (c - c0) / 3 + cdef double xr = r - r0 + cdef double xc = c - c0 + + r0 -= 1 + c0 -= 1 cdef double fc[4] cdef double fr[4] - cdef Py_ssize_t pr, pc + cdef long pr, pc # row-wise cubic interpolation - for pr in range(r0, r0 + 4): - for pc in range(c0, c0 + 4): - fc[pc - c0] = get_pixel2d(image, rows, cols, pr, pc, mode, cval) - fr[pr - r0] = cubic_interpolation(xc, fc) + for pr in range(4): + for pc in range(4): + fc[pc] = get_pixel2d(image, rows, cols, pr + r0, pc + c0, mode, cval) + fr[pr] = cubic_interpolation(xc, fc) # cubic interpolation for interpolated values of each row return cubic_interpolation(xr, fr) cdef inline double get_pixel2d(double* image, Py_ssize_t rows, Py_ssize_t cols, - Py_ssize_t r, Py_ssize_t c, char mode, double cval): + long r, long c, char mode, + double cval): """Get a pixel from the image, taking wrapping mode into consideration. Parameters @@ -248,7 +259,7 @@ cdef inline double get_pixel2d(double* image, Py_ssize_t rows, Py_ssize_t cols, """ if mode == 'C': - if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1): + if (r < 0) or (r >= rows) or (c < 0) or (c >= cols): return cval else: return image[r * cols + c] @@ -257,8 +268,8 @@ cdef inline double get_pixel2d(double* image, Py_ssize_t rows, Py_ssize_t cols, cdef inline double get_pixel3d(double* image, Py_ssize_t rows, Py_ssize_t cols, - Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c, Py_ssize_t d, - char mode, double cval): + Py_ssize_t dims, long r, long c, + long d, char mode, double cval): """Get a pixel from the image, taking wrapping mode into consideration. Parameters @@ -281,19 +292,18 @@ cdef inline double get_pixel3d(double* image, Py_ssize_t rows, Py_ssize_t cols, """ if mode == 'C': - if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1): + if (r < 0) or (r >= rows) or (c < 0) or (c >= cols): return cval else: return image[r * cols * dims + c * dims + d] else: return image[coord_map(rows, r, mode) * cols * dims + coord_map(cols, c, mode) * dims - + d] + + coord_map(dims, d, mode)] -cdef inline Py_ssize_t coord_map(Py_ssize_t dim, Py_ssize_t coord, char mode): - """ - Wrap a coordinate, according to a given mode. +cdef inline Py_ssize_t coord_map(Py_ssize_t dim, long coord, char mode): + """Wrap a coordinate, according to a given mode. Parameters ---------- diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 78387cf1..456570f2 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -1134,7 +1134,18 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, out = None - if order in range(4) and not map_args: + if order == 2: + # When fixing this issue, make sure to fix the branches further + # below in this function + warnings.warn("Bi-quadratic interpolation behavior has changed due " + "to a bug in the implementation of scikit-image. " + "The new version now serves as a wrapper " + "around SciPy's interpolation functions, which itself " + "is not verified to be a correct implementation. Until " + "skimage's implementation is fixed, we recommend " + "to use bi-linear or bi-cubic interpolation instead.") + + if order in (0, 1, 3) and not map_args: # use fast Cython version for specific interpolation orders and input matrix = None diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index b5206414..a2d6265f 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -1,4 +1,4 @@ -from numpy.testing import (assert_array_almost_equal, run_module_suite, +from numpy.testing import (assert_almost_equal, run_module_suite, assert_array_equal, assert_raises) import numpy as np from scipy.ndimage import map_coordinates @@ -22,10 +22,10 @@ def test_warp_tform(): tform = SimilarityTransform(scale=1, rotation=theta, translation=(0, 4)) x90 = warp(x, tform, order=1) - assert_array_almost_equal(x90, np.rot90(x)) + assert_almost_equal(x90, np.rot90(x)) x90 = warp(x, tform.inverse, order=1) - assert_array_almost_equal(x90, np.rot90(x)) + assert_almost_equal(x90, np.rot90(x)) def test_warp_callable(): @@ -37,7 +37,7 @@ def test_warp_callable(): shift = lambda xy: xy + 1 outx = warp(x, shift, order=1) - assert_array_almost_equal(outx, refx) + assert_almost_equal(outx, refx) def test_warp_matrix(): @@ -50,7 +50,7 @@ def test_warp_matrix(): # _warp_fast outx = warp(x, matrix, order=1) - assert_array_almost_equal(outx, refx) + assert_almost_equal(outx, refx) # check for ndimage.map_coordinates outx = warp(x, matrix, order=5) @@ -71,7 +71,7 @@ def test_warp_nd(): outx = warp(x, coords, order=0, cval=0) - assert_array_almost_equal(outx, refx) + assert_almost_equal(outx, refx) def test_warp_clip(): @@ -79,10 +79,10 @@ def test_warp_clip(): matrix = np.eye(3) outx = warp(x, matrix, order=0, clip=False) - assert_array_almost_equal(x, outx) + assert_almost_equal(x, outx) outx = warp(x, matrix, order=0, clip=True) - assert_array_almost_equal(x / 2, outx) + assert_almost_equal(x / 2, outx) def test_homography(): @@ -96,49 +96,14 @@ def test_homography(): x90 = warp(x, inverse_map=ProjectiveTransform(M).inverse, order=1) - assert_array_almost_equal(x90, np.rot90(x)) - - -def test_fast_homography(): - img = rgb2gray(data.astronaut()).astype(np.uint8) - img = img[:, :100] - - theta = np.deg2rad(30) - scale = 0.5 - tx, ty = 50, 50 - - H = np.eye(3) - S = scale * np.sin(theta) - C = scale * np.cos(theta) - - H[:2, :2] = [[C, -S], [S, C]] - H[:2, 2] = [tx, ty] - - tform = ProjectiveTransform(H) - coords = warp_coords(tform.inverse, (img.shape[0], img.shape[1])) - - for order in range(4): - for mode in ('constant', 'reflect', 'wrap', 'nearest'): - p0 = map_coordinates(img, coords, mode=mode, order=order) - p1 = warp(img, tform, mode=mode, order=order) - - # 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 + assert_almost_equal(x90, np.rot90(x)) def test_rotate(): x = np.zeros((5, 5), dtype=np.double) x[1, 1] = 1 x90 = rotate(x, 90) - assert_array_almost_equal(x90, np.rot90(x)) + assert_almost_equal(x90, np.rot90(x)) def test_rotate_resize(): @@ -158,9 +123,9 @@ def test_rotate_center(): refx = np.zeros((10, 10), dtype=np.double) refx[2, 5] = 1 x20 = rotate(x, 20, order=0, center=(0, 0)) - assert_array_almost_equal(x20, refx) + assert_almost_equal(x20, refx) x0 = rotate(x20, -20, order=0, center=(0, 0)) - assert_array_almost_equal(x0, x) + assert_almost_equal(x0, x) def test_rescale(): @@ -170,7 +135,7 @@ def test_rescale(): scaled = rescale(x, 2, order=0) ref = np.zeros((10, 10)) ref[2:4, 2:4] = 1 - assert_array_almost_equal(scaled, ref) + assert_almost_equal(scaled, ref) # different scale factors x = np.zeros((5, 5), dtype=np.double) @@ -178,7 +143,7 @@ def test_rescale(): scaled = rescale(x, (2, 1), order=0) ref = np.zeros((10, 5)) ref[2:4, 1] = 1 - assert_array_almost_equal(scaled, ref) + assert_almost_equal(scaled, ref) def test_resize2d(): @@ -187,7 +152,7 @@ def test_resize2d(): resized = resize(x, (10, 10), order=0) ref = np.zeros((10, 10)) ref[2:4, 2:4] = 1 - assert_array_almost_equal(resized, ref) + assert_almost_equal(resized, ref) def test_resize3d_keep(): @@ -197,9 +162,9 @@ def test_resize3d_keep(): resized = resize(x, (10, 10), order=0) ref = np.zeros((10, 10, 3)) ref[2:4, 2:4, :] = 1 - assert_array_almost_equal(resized, ref) + assert_almost_equal(resized, ref) resized = resize(x, (10, 10, 3), order=0) - assert_array_almost_equal(resized, ref) + assert_almost_equal(resized, ref) def test_resize3d_resize(): @@ -209,7 +174,7 @@ def test_resize3d_resize(): resized = resize(x, (10, 10, 1), order=0) ref = np.zeros((10, 10, 1)) ref[2:4, 2:4] = 1 - assert_array_almost_equal(resized, ref) + assert_almost_equal(resized, ref) def test_resize3d_bilinear(): @@ -223,7 +188,7 @@ def test_resize3d_bilinear(): ref[1:5, 2:4, :] = 0.09375 ref[2:4, 1:5, :] = 0.09375 ref[2:4, 2:4, :] = 0.28125 - assert_array_almost_equal(resized, ref) + assert_almost_equal(resized, ref) def test_swirl():