diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd index 3e8e74e9..40df6524 100644 --- a/skimage/_shared/interpolation.pxd +++ b/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) diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index 63f5b70c..bcebd78d 100644 --- a/skimage/_shared/interpolation.pyx +++ b/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, round(r), 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 = r + cdef int c0 = 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): diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index f11f1b44..6f895fa7 100644 --- a/skimage/transform/_geometric.py +++ b/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 diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index 18b1a5b0..37a577f7 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/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): diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index d5ffff59..41d1a46a 100644 --- a/skimage/transform/tests/test_warps.py +++ b/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)