From 167b43fc33a6fb1f5cf25cd3ac1f2707c6073247 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sun, 25 Sep 2011 16:45:52 -0700 Subject: [PATCH] ENH: fast_homography: 20x speed improvement. --- scikits/image/transform/_project.pyx | 71 ++++++++++++++++------------ 1 file changed, 42 insertions(+), 29 deletions(-) diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx index 5cbbb608..06bc53cb 100644 --- a/scikits/image/transform/_project.pyx +++ b/scikits/image/transform/_project.pyx @@ -6,26 +6,24 @@ cimport numpy as np import numpy as np import cython +from cython.operator import dereference + np.import_array() cdef extern from "math.h": double floor(double) double fmod(double, double) -cdef double get_pixel(np.ndarray image, int r, int c, char mode, - double cval=0): - cdef np.ndarray[dtype=np.double_t, ndim=2] img = image - cdef int rows = img.shape[0] - cdef int cols = img.shape[1] - +cdef double get_pixel(double *image, int rows, int cols, + int r, int c, char mode, double cval=0): if mode == 'C': if (r < 0) or (r >= cols) or (c < 0) or (c >= cols): return cval else: - return img[r, c] + return image[r * cols + c] else: - return img[coord_map(rows, r, mode), - 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): dim = dim - 1 @@ -42,22 +40,28 @@ cdef int coord_map(int dim, int coord, char mode): return coord -cdef tf(double x, double y, H): - cdef np.ndarray[np.double_t, ndim=2] M = H +cdef tf(double x, double y, double* H, double *x_, double *y_): cdef double xx, yy, zz - xx = M[0, 0] * x + M[0, 1] * y + M[0, 2] - yy = M[1, 0] * x + M[1, 1] * y + M[1, 2] - zz = M[2, 0] * x + M[2, 1] * y + M[2, 2] + ## print + ## print H[0], H[1], H[2] + ## print H[3], H[4], H[5] + ## print H[6], H[7], H[8] + ## print - xx /= zz - yy /= zz + xx = H[0] * x + H[1] * y + H[2] + yy = H[3] * x + H[4] * y + H[5] + zz = H[6] * x + H[7] * y + H[8] - return xx, yy + xx = xx / zz + yy = yy / zz + + x_[0] = xx + y_[0] = yy @cython.boundscheck(False) def homography(np.ndarray image, np.ndarray H, output_shape=None, - mode='C', double cval=0): + mode='constant', double cval=0): """ Projective transformation (homography). @@ -93,9 +97,8 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, Shape of the output image generated. order : int Order of splines used in interpolation. - mode : {'C', 'M', 'W'} + mode : {'constant', 'mirror', 'wrap'} How to handle values outside the image borders. - Constant, Mirror or Wrap. cval : string Used in conjunction with mode 'C' (constant), the value outside the image boundaries. @@ -103,12 +106,18 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, """ cdef np.ndarray[dtype=np.double_t, ndim=2] img = image - cdef np.ndarray[dtype=np.double_t, ndim=2] M = np.linalg.inv(H) + cdef np.ndarray[dtype=np.double_t, ndim=2] M = \ + np.ascontiguousarray(np.linalg.inv(H)) - if mode not in ('C', 'W', 'M'): + if mode not in ('constant', 'wrap', 'mirror'): raise ValueError("Invalid mode specified. Please use " - "C [constant], W [wrap] or M [mirror].") - cdef char mode_c = ord(mode[0]) + "`constant`, `wrap` or `mirror`.") + if mode == 'constant': + mode_c = ord('C') + elif mode == 'wrap': + mode_c = ord('W') + elif mode == 'mirror': + mode_c = ord('M') cdef int out_r, out_c, columns, rows if output_shape is None: @@ -130,17 +139,21 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, for tfr in range(out_r): for tfc in range(out_c): - c, r = tf(tfc, tfr, M) + 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, r_int, c_int, mode_c) - y1 = get_pixel(img, r_int + 1, c_int, mode_c) - y2 = get_pixel(img, r_int + 1, c_int + 1, mode_c) - y3 = get_pixel(img, r_int, c_int + 1, mode_c) + 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 + \