From 0ca7933a7d1cf573d12fcf15f72b9f03afeb84f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 16 Aug 2012 19:08:13 +0200 Subject: [PATCH] Share bilinear interpolation function for other code Make bilinear_interpolation function callable by other cython code and change fast_homography to use this function. Also, improve performance of some functions by making them inline and fix broken test of fast_homography. --- skimage/transform/_project.pxd | 12 +++ skimage/transform/_project.pyx | 112 +++++++++++--------------- skimage/transform/tests/test_warps.py | 23 +++--- 3 files changed, 72 insertions(+), 75 deletions(-) create mode 100644 skimage/transform/_project.pxd diff --git a/skimage/transform/_project.pxd b/skimage/transform/_project.pxd new file mode 100644 index 00000000..36a32d2d --- /dev/null +++ b/skimage/transform/_project.pxd @@ -0,0 +1,12 @@ +cimport numpy as np +import numpy as np + + +cdef inline double bilinear_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=*) + +cdef inline int coord_map(int dim, int coord, char mode) \ No newline at end of file diff --git a/skimage/transform/_project.pyx b/skimage/transform/_project.pyx index 734d6bff..963670c6 100644 --- a/skimage/transform/_project.pyx +++ b/skimage/transform/_project.pyx @@ -1,38 +1,50 @@ -#cython: cdivison=True boundscheck=False +#cython: cdivison=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False -__all__ = ['homography'] - -cimport cython cimport numpy as np - import numpy as np -import cython - from cython.operator import dereference +from libc.math cimport ceil, floor -np.import_array() -cdef extern from "math.h": - double floor(double) - double fmod(double, double) +cdef inline double bilinear_interpolation(double* image, int rows, int cols, + double r, double c, char mode, + double cval=0): + cdef double dr, dc + cdef int minr, minc, maxr, maxc -cdef double get_pixel(double *image, int rows, int cols, - int r, int c, char mode, double cval=0): + minr = floor(r) + minc = floor(c) + maxr = ceil(r) + maxc = ceil(c) + dr = r - minr + dc = c - minc + top = (1 - dc) * get_pixel(image, rows, cols, minr, minc, mode, cval) \ + + dc * get_pixel(image, rows, cols, minr, maxc, mode, cval) + bottom = (1 - dc) * get_pixel(image, rows, cols, maxr, minc, mode, cval) \ + + dc * get_pixel(image, rows, cols, maxr, maxc, mode, cval) + return (1 - dr) * top + dr * bottom + + +cdef inline double get_pixel(double* image, int rows, int cols, int r, int c, + char mode, double cval=0): """Get a pixel from the image, taking wrapping mode into consideration. Parameters ---------- - image : *double + image : array of dtype double Input image. - rows, cols : int - Dimensions of image. + rows, cols: int + Shape of image. r, c : int Position at which to get the pixel. mode : {'C', 'W', 'M'} - Wrapping mode. Constant, Wrap or Mirror. + Wrapping mode. Constant, Wrap or Mirror. cval : double - Constant value to use for mode constant. - + Constant value to use for constant mode. + """ if mode == 'C': if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1): @@ -40,13 +52,13 @@ cdef double get_pixel(double *image, int rows, int cols, else: return image[r * cols + c] else: - return image[coord_map(rows, r, mode) * cols + - coord_map(cols, c, mode)] + return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)] -cdef int coord_map(int dim, int coord, char mode): + +cdef inline int coord_map(int dim, int coord, char mode): """ - Wrap a coordinate, according to a given dimension and mode. - + Wrap a coordinate, according to a given mode. + Parameters ---------- dim : int @@ -56,7 +68,7 @@ cdef int coord_map(int dim, int coord, char mode): mode : {'W', 'M'} Whether to wrap or mirror the coordinate if it falls outside [0, dim). - + """ dim = dim - 1 if mode == 'M': # mirror @@ -79,7 +91,8 @@ cdef int coord_map(int dim, int coord, char mode): return coord -cdef tf(double x, double y, double* H, double *x_, double *y_): + +cdef inline tf(double x, double y, double* H, double *x_, double *y_): """Apply a homography to a coordinate. Parameters @@ -98,18 +111,15 @@ cdef tf(double x, double y, double* H, double *x_, double *y_): yy = H[3] * x + H[4] * y + H[5] zz = H[6] * x + H[7] * y + H[8] - xx = xx / zz - yy = yy / zz + x_[0] = xx / zz + y_[0] = yy / zz - x_[0] = xx - y_[0] = yy -@cython.boundscheck(False) def homography(np.ndarray image, np.ndarray H, output_shape=None, mode='constant', double cval=0): """ Projective transformation (homography). - + Perform a projective transformation (homography) of a floating point image, using bi-linear interpolation. @@ -140,8 +150,6 @@ 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 : int - Order of splines used in interpolation. mode : {'constant', 'mirror', 'wrap'} How to handle values outside the image borders. cval : string @@ -150,8 +158,7 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, """ - cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] img = \ - np.ascontiguousarray(image, dtype=np.double) + cdef np.ndarray[dtype=np.double_t, ndim=2] img = image.astype(np.double) cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] M = \ np.ascontiguousarray(np.linalg.inv(H)) @@ -165,7 +172,6 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, elif mode == 'mirror': mode_c = ord('M') - cdef int out_r, out_c, columns, rows if output_shape is None: out_r = img.shape[0] out_c = img.shape[1] @@ -173,37 +179,17 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, out_r = output_shape[0] out_c = output_shape[1] - rows = img.shape[0] - columns = img.shape[1] - cdef np.ndarray[dtype=np.double_t, ndim=2] out = \ np.zeros((out_r, out_c), dtype=np.double) - - cdef int tfr, tfc, r_int, c_int - cdef double y0, y1, y2, y3 - cdef double r, c, z, t, u + + cdef int tfr, tfc + cdef double r, c + cdef int rows = img.shape[0] + cdef int cols = img.shape[1] for tfr in range(out_r): for tfc in range(out_c): tf(tfc, tfr, M.data, &c, &r) - r_int = floor(r) - c_int = floor(c) - - t = r - r_int - u = c - c_int - - y0 = get_pixel(img.data, rows, columns, - r_int, c_int, mode_c) - y1 = get_pixel(img.data, rows, columns, - r_int + 1, c_int, mode_c) - y2 = get_pixel(img.data, rows, columns, - r_int + 1, c_int + 1, mode_c) - y3 = get_pixel(img.data, rows, columns, - r_int, c_int + 1, mode_c) - - out[tfr, tfc] = \ - (1 - t) * (1 - u) * y0 + \ - t * (1 - u) * y1 + \ - t * u * y2 + (1 - t) * u * y3; + 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 7c2e52f2..f35b7f57 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -2,7 +2,7 @@ from numpy.testing import assert_array_almost_equal, run_module_suite import numpy as np from skimage.transform import (warp, homography, fast_homography, - SimilarityTransform) + SimilarityTransform, ProjectiveTransform) from skimage import transform as tf, data, img_as_float from skimage.color import rgb2gray @@ -34,7 +34,7 @@ def test_homography(): def test_fast_homography(): - img = rgb2gray(data.lena()).astype(np.uint8) + img = rgb2gray(data.lena()) img = img[:, :100] theta = np.deg2rad(30) @@ -49,20 +49,19 @@ def test_fast_homography(): H[:2, 2] = [tx, ty] for mode in ('constant', 'mirror', 'wrap'): - p0 = homography(img, H, mode=mode, order=1) + p0 = warp(img, ProjectiveTransform(H).inverse, mode=mode, order=1) p1 = fast_homography(img, H, mode=mode) - p1 = np.round(p1) - ## 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() + # 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.2 + assert d < 0.001 def test_swirl():