From d87ed28d8eb1197f0f443d208a9853ff1be9bd94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 21 Aug 2012 08:20:57 +0200 Subject: [PATCH 01/13] Add new package for shared code --- skimage/_shared/__init__.py | 0 skimage/_shared/setup.py | 28 ++++++++++++++++++++++++++++ 2 files changed, 28 insertions(+) create mode 100644 skimage/_shared/__init__.py create mode 100644 skimage/_shared/setup.py diff --git a/skimage/_shared/__init__.py b/skimage/_shared/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/skimage/_shared/setup.py b/skimage/_shared/setup.py new file mode 100644 index 00000000..913a020c --- /dev/null +++ b/skimage/_shared/setup.py @@ -0,0 +1,28 @@ +#!/usr/bin/env python + +import os + +from skimage._build import cython + +base_path = os.path.abspath(os.path.dirname(__file__)) + + +def configuration(parent_package='', top_path=None): + from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs + + config = Configuration('_shared', parent_package, top_path) + config.add_data_dir('tests') + + return config + + +if __name__ == '__main__': + from numpy.distutils.core import setup + setup(maintainer='Scikits-image Developers', + author='Scikits-image Developers', + maintainer_email='scikits-image@googlegroups.com', + description='Transforms', + url='https://github.com/scikits-image/scikits-image', + license='SciPy License (BSD Style)', + **(configuration(top_path='').todict()) + ) From 824997af0a4c0070c5ddaade9c21efaaaac61e9c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 21 Aug 2012 08:32:20 +0200 Subject: [PATCH 02/13] Move bilinear interpolation code to shared package --- .../interpolation.pxd} | 8 +- skimage/_shared/interpolation.pyx | 104 ++++++++++++++++++ skimage/_shared/setup.py | 5 + skimage/feature/_texture.pyx | 2 +- skimage/setup.py | 1 + skimage/transform/_project.pyx | 102 +---------------- skimage/transform/setup.py | 2 +- 7 files changed, 118 insertions(+), 106 deletions(-) rename skimage/{transform/_project.pxd => _shared/interpolation.pxd} (52%) create mode 100644 skimage/_shared/interpolation.pyx diff --git a/skimage/transform/_project.pxd b/skimage/_shared/interpolation.pxd similarity index 52% rename from skimage/transform/_project.pxd rename to skimage/_shared/interpolation.pxd index dee51026..7ad5d223 100644 --- a/skimage/transform/_project.pxd +++ b/skimage/_shared/interpolation.pxd @@ -1,7 +1,9 @@ -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) diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx new file mode 100644 index 00000000..defa2e9b --- /dev/null +++ b/skimage/_shared/interpolation.pyx @@ -0,0 +1,104 @@ +#cython: cdivison=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False +from libc.math cimport ceil, floor + + +cdef inline double bilinear_interpolation(double* image, int rows, int cols, + double r, double c, char mode, + double cval=0): + """Bilinear 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. + + """ + cdef double dr, dc + cdef int minr, minc, maxr, maxc + + 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 array + Input 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. + cval : double + Constant value to use for constant mode. + + """ + if mode == 'C': + if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1): + return cval + else: + return image[r * cols + c] + else: + return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)] + + +cdef inline int coord_map(int dim, int coord, char mode): + """ + Wrap a coordinate, according to a given mode. + + Parameters + ---------- + dim : int + Maximum coordinate. + coord : int + Coord provided by user. May be < 0 or > dim. + mode : {'W', 'M'} + Whether to wrap or mirror the coordinate if it + falls outside [0, dim). + + """ + dim = dim - 1 + if mode == 'M': # mirror + if (coord < 0): + # How many times times does the coordinate wrap? + if ((-coord / dim) % 2 != 0): + return dim - (-coord % dim) + else: + return (-coord % dim) + elif (coord > dim): + if ((coord / dim) % 2 != 0): + return (dim - (coord % dim)) + else: + return (coord % dim) + elif mode == 'W': # wrap + if (coord < 0): + return (dim - (-coord % dim)) + elif (coord > dim): + return (coord % dim) + + return coord diff --git a/skimage/_shared/setup.py b/skimage/_shared/setup.py index 913a020c..6e4d1b6b 100644 --- a/skimage/_shared/setup.py +++ b/skimage/_shared/setup.py @@ -13,6 +13,11 @@ def configuration(parent_package='', top_path=None): config = Configuration('_shared', parent_package, top_path) config.add_data_dir('tests') + cython(['interpolation.pyx'], working_path=base_path) + + config.add_extension('interpolation', sources=['interpolation.c'], + include_dirs=[get_numpy_include_dirs()]) + return config diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index ca51e177..20b61513 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -5,7 +5,7 @@ import numpy as np cimport numpy as np from libc.math cimport sin, cos, abs -from skimage.transform._project cimport bilinear_interpolation +from skimage._shared.interpolation cimport bilinear_interpolation def _glcm_loop(np.ndarray[dtype=np.uint8_t, ndim=2, diff --git a/skimage/setup.py b/skimage/setup.py index 02c0b52d..1082ba07 100644 --- a/skimage/setup.py +++ b/skimage/setup.py @@ -6,6 +6,7 @@ def configuration(parent_package='', top_path=None): config = Configuration('skimage', parent_package, top_path) + config.add_subpackage('_shared') config.add_subpackage('color') config.add_subpackage('data') config.add_subpackage('draw') diff --git a/skimage/transform/_project.pyx b/skimage/transform/_project.pyx index 3f31f229..300ba543 100644 --- a/skimage/transform/_project.pyx +++ b/skimage/transform/_project.pyx @@ -5,107 +5,7 @@ cimport numpy as np import numpy as np -from cython.operator import dereference -from libc.math cimport ceil, floor - - -cdef inline double bilinear_interpolation(double* image, int rows, int cols, - double r, double c, char mode, - double cval=0): - """Bilinear 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. - - """ - cdef double dr, dc - cdef int minr, minc, maxr, maxc - - 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 array - Input 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. - cval : double - Constant value to use for constant mode. - - """ - if mode == 'C': - if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1): - return cval - else: - return image[r * cols + c] - else: - return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)] - - -cdef inline int coord_map(int dim, int coord, char mode): - """ - Wrap a coordinate, according to a given mode. - - Parameters - ---------- - dim : int - Maximum coordinate. - coord : int - Coord provided by user. May be < 0 or > dim. - mode : {'W', 'M'} - Whether to wrap or mirror the coordinate if it - falls outside [0, dim). - - """ - dim = dim - 1 - if mode == 'M': # mirror - if (coord < 0): - # How many times times does the coordinate wrap? - if ((-coord / dim) % 2 != 0): - return dim - (-coord % dim) - else: - return (-coord % dim) - elif (coord > dim): - if ((coord / dim) % 2 != 0): - return (dim - (coord % dim)) - else: - return (coord % dim) - elif mode == 'W': # wrap - if (coord < 0): - return (dim - (-coord % dim)) - elif (coord > dim): - return (coord % dim) - - return coord +from skimage._shared.interpolation cimport bilinear_interpolation cdef inline _matrix_transform(double x, double y, double* H, double *x_, diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py index 4ecb6ba4..75210b54 100644 --- a/skimage/transform/setup.py +++ b/skimage/transform/setup.py @@ -20,7 +20,7 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('_project', sources=['_project.c'], - include_dirs=[get_numpy_include_dirs()]) + include_dirs=[get_numpy_include_dirs(), '../_shared']) return config From 72473848824f3d837efed4c63ba64dac1513e07e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 21 Aug 2012 08:59:29 +0200 Subject: [PATCH 03/13] Add integrate function to shared package --- skimage/_shared/setup.py | 3 +++ skimage/_shared/transform.pxd | 6 +++++ skimage/_shared/transform.pyx | 45 +++++++++++++++++++++++++++++++++ skimage/feature/_template.pyx | 47 ++--------------------------------- skimage/feature/setup.py | 5 ++-- 5 files changed, 58 insertions(+), 48 deletions(-) create mode 100644 skimage/_shared/transform.pxd create mode 100644 skimage/_shared/transform.pyx diff --git a/skimage/_shared/setup.py b/skimage/_shared/setup.py index 6e4d1b6b..ed48916e 100644 --- a/skimage/_shared/setup.py +++ b/skimage/_shared/setup.py @@ -14,9 +14,12 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('tests') cython(['interpolation.pyx'], working_path=base_path) + cython(['transform.pyx'], working_path=base_path) config.add_extension('interpolation', sources=['interpolation.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('transform', sources=['transform.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/_shared/transform.pxd b/skimage/_shared/transform.pxd new file mode 100644 index 00000000..2953ac20 --- /dev/null +++ b/skimage/_shared/transform.pxd @@ -0,0 +1,6 @@ +cimport numpy as cnp +import numpy as np + + +cdef float integrate(cnp.ndarray[float, ndim=2, mode="c"] sat, + int r0, int c0, int r1, int c1) diff --git a/skimage/_shared/transform.pyx b/skimage/_shared/transform.pyx new file mode 100644 index 00000000..b6649852 --- /dev/null +++ b/skimage/_shared/transform.pyx @@ -0,0 +1,45 @@ +#cython: cdivison=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False +cimport numpy as cnp +import numpy as np + + +cdef float integrate(cnp.ndarray[float, ndim=2, mode="c"] sat, + int r0, int c0, int r1, int c1): + """ + Using a summed area table / integral image, calculate the sum + over a given window. + + This function is the same as the `integrate` function in + `skimage.transform.integrate`, but this Cython version significantly + speeds up the code. + + Parameters + ---------- + sat : ndarray of float + Summed area table / integral image. + r0, c0 : int + Top-left corner of block to be summed. + r1, c1 : int + Bottom-right corner of block to be summed. + + Returns + ------- + S : int + Sum over the given window. + """ + cdef float S = 0 + + S += sat[r1, c1] + + if (r0 - 1 >= 0) and (c0 - 1 >= 0): + S += sat[r0 - 1, c0 - 1] + + if (r0 - 1 >= 0): + S -= sat[r0 - 1, c1] + + if (c0 - 1 >= 0): + S -= sat[r1, c0 - 1] + return S diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index b83761a8..58d48524 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -35,51 +35,8 @@ cimport numpy as np import numpy as np from scipy.signal import fftconvolve from skimage.transform import integral - - -cdef extern from "math.h": - float sqrt(float x) - float fabs(float x) - - -@cython.boundscheck(False) -cdef float integrate(np.ndarray[float, ndim=2, mode="c"] sat, - int r0, int c0, int r1, int c1): - """ - Using a summed area table / integral image, calculate the sum - over a given window. - - This function is the same as the `integrate` function in - `skimage.transform.integrate`, but this Cython version significantly - speeds up the code. - - Parameters - ---------- - sat : ndarray of float - Summed area table / integral image. - r0, c0 : int - Top-left corner of block to be summed. - r1, c1 : int - Bottom-right corner of block to be summed. - - Returns - ------- - S : int - Sum over the given window. - """ - cdef float S = 0 - - S += sat[r1, c1] - - if (r0 - 1 >= 0) and (c0 - 1 >= 0): - S += sat[r0 - 1, c0 - 1] - - if (r0 - 1 >= 0): - S -= sat[r0 - 1, c1] - - if (c0 - 1 >= 0): - S -= sat[r1, c0 - 1] - return S +from libc.math cimport sqrt, fabs +from skimage._shared.transform cimport integrate @cython.boundscheck(False) diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 0b0b80bd..9c9074be 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -16,10 +16,9 @@ def configuration(parent_package='', top_path=None): cython(['_template.pyx'], working_path=base_path) config.add_extension('_texture', sources=['_texture.c'], - include_dirs=[get_numpy_include_dirs(), - '../transform']) + include_dirs=[get_numpy_include_dirs(), '../_shared']) config.add_extension('_template', sources=['_template.c'], - include_dirs=[get_numpy_include_dirs()]) + include_dirs=[get_numpy_include_dirs(), '../_shared']) return config From ba171937a134d8c5d4fedd807ca548f114feb942 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 21 Aug 2012 09:05:18 +0200 Subject: [PATCH 04/13] Remove unused imports --- skimage/_shared/transform.pxd | 1 - skimage/_shared/transform.pyx | 1 - 2 files changed, 2 deletions(-) diff --git a/skimage/_shared/transform.pxd b/skimage/_shared/transform.pxd index 2953ac20..0edc22a4 100644 --- a/skimage/_shared/transform.pxd +++ b/skimage/_shared/transform.pxd @@ -1,5 +1,4 @@ cimport numpy as cnp -import numpy as np cdef float integrate(cnp.ndarray[float, ndim=2, mode="c"] sat, diff --git a/skimage/_shared/transform.pyx b/skimage/_shared/transform.pyx index b6649852..e4bbfa74 100644 --- a/skimage/_shared/transform.pyx +++ b/skimage/_shared/transform.pyx @@ -3,7 +3,6 @@ #cython: nonecheck=False #cython: wraparound=False cimport numpy as cnp -import numpy as np cdef float integrate(cnp.ndarray[float, ndim=2, mode="c"] sat, 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 05/13] 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(): From a08779e06a72210d061498f0ff2f8201ec1e705b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 21 Aug 2012 15:15:27 +0200 Subject: [PATCH 06/13] Use predefined header files from Cython --- skimage/graph/_spath.pyx | 3 +- skimage/io/_plugins/_colormixer.pyx | 19 +++-------- skimage/segmentation/_quickshift.pyx | 6 +--- skimage/transform/_hough_transform.pyx | 45 +++++++++++--------------- 4 files changed, 25 insertions(+), 48 deletions(-) diff --git a/skimage/graph/_spath.pyx b/skimage/graph/_spath.pyx index f342021d..1624bc7d 100644 --- a/skimage/graph/_spath.pyx +++ b/skimage/graph/_spath.pyx @@ -2,9 +2,8 @@ import _mcp cimport _mcp +from libc.math cimport fabs -cdef extern from "math.h": - double fabs(double f) cdef class MCP_Diff(_mcp.MCP): """MCP_Diff(costs, offsets=None, fully_connected=True) diff --git a/skimage/io/_plugins/_colormixer.pyx b/skimage/io/_plugins/_colormixer.pyx index 5d0dd1b9..ace45c94 100644 --- a/skimage/io/_plugins/_colormixer.pyx +++ b/skimage/io/_plugins/_colormixer.pyx @@ -8,15 +8,10 @@ integers, so currently the only way to clip results efficiently one. """ - +import cython import numpy as np cimport numpy as np - -import cython - -cdef extern from "math.h": - float exp(float) nogil - float pow(float, float) nogil +from libc.math cimport exp, pow @cython.boundscheck(False) @@ -189,7 +184,6 @@ def sigmoid_gamma(np.ndarray[np.uint8_t, ndim=3] img, img[i,j,2] = lut[stateimg[i,j,2]] - @cython.boundscheck(False) def gamma(np.ndarray[np.uint8_t, ndim=3] img, np.ndarray[np.uint8_t, ndim=3] stateimg, @@ -219,7 +213,6 @@ def gamma(np.ndarray[np.uint8_t, ndim=3] img, img[i,j,2] = lut[stateimg[i,j,2]] - @cython.cdivision(True) cdef void rgb_2_hsv(float* RGB, float* HSV) nogil: cdef float R, G, B, H, S, V, MAX, MIN @@ -283,6 +276,7 @@ cdef void rgb_2_hsv(float* RGB, float* HSV) nogil: HSV[1] = S HSV[2] = V + @cython.cdivision(True) cdef void hsv_2_rgb(float* HSV, float* RGB) nogil: cdef float H, S, V @@ -388,6 +382,7 @@ def py_hsv_2_rgb(H, S, V): return (R, G, B) + def py_rgb_2_hsv(R, G, B): '''Convert an HSV value to RGB. @@ -561,9 +556,3 @@ def hsv_multiply(np.ndarray[np.uint8_t, ndim=3] img, img[i, j, 0] = RGB[0] img[i, j, 1] = RGB[1] img[i, j, 2] = RGB[2] - - - - - - diff --git a/skimage/segmentation/_quickshift.pyx b/skimage/segmentation/_quickshift.pyx index 57009ad1..57be1040 100644 --- a/skimage/segmentation/_quickshift.pyx +++ b/skimage/segmentation/_quickshift.pyx @@ -1,6 +1,7 @@ import numpy as np cimport numpy as np cimport cython +from libc.math cimport exp, sqrt from itertools import product from scipy import ndimage @@ -9,11 +10,6 @@ from ..util import img_as_float from ..color import rgb2lab -cdef extern from "math.h": - double exp(double) - double sqrt(double) - - @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx index b34b28e2..906e4464 100644 --- a/skimage/transform/_hough_transform.pyx +++ b/skimage/transform/_hough_transform.pyx @@ -2,27 +2,20 @@ cimport cython import numpy as np cimport numpy as np from random import randint +from libc.math cimport abs, fabs, sqrt, ceil, floor, round +from libc.stdlib cimport rand + + np.import_array() -cdef extern from "stdlib.h": - int rand() - -cdef extern from "math.h": - int abs(int) - double fabs(double) - double sqrt(double) - double ceil(double) - double floor(double) - -cdef double round(double val): - return floor(val + 0.5); cdef double PI_2 = 1.5707963267948966 cdef double NEG_PI_2 = -PI_2 + @cython.boundscheck(False) def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None): - + if img.ndim != 2: raise ValueError('The input image must be 2D.') @@ -31,7 +24,7 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None): cdef np.ndarray[ndim=1, dtype=np.double_t] stheta if theta is None: - theta = np.linspace(PI_2, NEG_PI_2, 180) + theta = np.linspace(PI_2, NEG_PI_2, 180) ctheta = np.cos(theta) stheta = np.sin(theta) @@ -39,14 +32,14 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None): # compute the bins and allocate the accumulator array cdef np.ndarray[ndim=2, dtype=np.uint64_t] accum cdef np.ndarray[ndim=1, dtype=np.double_t] bins - cdef int max_distance, offset + cdef int max_distance, offset - max_distance = 2 * ceil((sqrt(img.shape[0] * img.shape[0] + + max_distance = 2 * ceil((sqrt(img.shape[0] * img.shape[0] + img.shape[1] * img.shape[1]))) accum = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64) bins = np.linspace(-max_distance / 2.0, max_distance / 2.0, max_distance) offset = max_distance / 2 - + # compute the nonzero indexes cdef np.ndarray[ndim=1, dtype=np.npy_intp] x_idxs, y_idxs y_idxs, x_idxs = np.PyArray_Nonzero(img) @@ -58,7 +51,7 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None): nthetas = theta.shape[0] for i in range(nidxs): x = x_idxs[i] - y = y_idxs[i] + y = y_idxs[i] for j in range(nthetas): accum_idx = round((ctheta[j] * x + stheta[j] * y)) + offset accum[accum_idx, j] += 1 @@ -94,7 +87,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \ # maximum line number cutoff cdef int lines_max = 2 ** 15 cdef int xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, good_line, count - max_distance = 2 * ceil((sqrt(img.shape[0] * img.shape[0] + + max_distance = 2 * ceil((sqrt(img.shape[0] * img.shape[0] + img.shape[1] * img.shape[1]))) accum = np.zeros((max_distance, theta.shape[0]), dtype=np.int64) offset = max_distance / 2 @@ -114,11 +107,11 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \ # select random non-zero point count = len(points) if count == 0: - break + break index = rand() % (count) x = points[index][0] y = points[index][1] - del points[index] + del points[index] # if previously eliminated, skip if not mask[y, x]: continue @@ -147,7 +140,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \ dx0 = 1 else: dx0 = -1 - dy0 = round(b * (1 << shift) / fabs(a)) + dy0 = round(b * (1 << shift) / fabs(a)) y0 = (y0 << shift) + (1 << (shift - 1)) else: if b > 0: @@ -156,7 +149,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \ dy0 = -1 dx0 = round(a * (1 << shift) / fabs(b)) x0 = (x0 << shift) + (1 << (shift - 1)) - + # pass 1: walk the line, merging lines less than specified gap length for k in range(2): gap = 0 @@ -208,9 +201,9 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \ x1 = px >> shift y1 = py # if non-zero point found, continue the line - if mask[y1, x1]: - if good_line: - accum_idx = round((ctheta[j] * x1 + stheta[j] * y1)) + offset + if mask[y1, x1]: + if good_line: + accum_idx = round((ctheta[j] * x1 + stheta[j] * y1)) + offset accum[accum_idx, max_theta] -= 1 mask[y1, x1] = 0 # exit when the point is the line end From c471d475eb5f62b930bf10ec8e5da08682a33751 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 21 Aug 2012 15:43:26 +0200 Subject: [PATCH 07/13] Refactor pnpoly function in Cython and add to shared package --- skimage/_shared/geometry.pxd | 7 ++++ skimage/_shared/geometry.pyx | 27 +++++++++++++ skimage/_shared/setup.py | 2 + skimage/draw/_draw.pyx | 8 +--- skimage/draw/setup.py | 2 +- skimage/morphology/_pnpoly.h | 72 ---------------------------------- skimage/morphology/_pnpoly.pyx | 17 +++----- skimage/morphology/setup.py | 2 +- 8 files changed, 45 insertions(+), 92 deletions(-) create mode 100644 skimage/_shared/geometry.pxd create mode 100644 skimage/_shared/geometry.pyx delete mode 100644 skimage/morphology/_pnpoly.h diff --git a/skimage/_shared/geometry.pxd b/skimage/_shared/geometry.pxd new file mode 100644 index 00000000..b228b23a --- /dev/null +++ b/skimage/_shared/geometry.pxd @@ -0,0 +1,7 @@ + +cdef inline unsigned char point_in_polygon(int nr_verts, double *xp, double *yp, + double x, double y) + +cdef void points_in_polygon(int nr_verts, double *xp, double *yp, + int nr_points, double *x, double *y, + unsigned char *result) diff --git a/skimage/_shared/geometry.pyx b/skimage/_shared/geometry.pyx new file mode 100644 index 00000000..de1ed4f4 --- /dev/null +++ b/skimage/_shared/geometry.pyx @@ -0,0 +1,27 @@ +#cython: cdivison=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + + +cdef inline unsigned char point_in_polygon(int nr_verts, double *xp, double *yp, + double x, double y): + cdef int i + cdef unsigned char c = 0 + cdef int j = nr_verts - 1 + for i in range(nr_verts): + if ( + (((yp[i] <= y) and (y < yp[j])) or + ((yp[j] <= y) and (y < yp[i]))) + and (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]) + ): + c = not c + j = i + return c + +cdef void points_in_polygon(int nr_verts, double *xp, double *yp, + int nr_points, double *x, double *y, + unsigned char *result): + cdef int n + for n in range(nr_points): + result[n] = point_in_polygon(nr_verts, xp, yp, x[n], y[n]) diff --git a/skimage/_shared/setup.py b/skimage/_shared/setup.py index ed48916e..81113656 100644 --- a/skimage/_shared/setup.py +++ b/skimage/_shared/setup.py @@ -13,9 +13,11 @@ def configuration(parent_package='', top_path=None): config = Configuration('_shared', parent_package, top_path) config.add_data_dir('tests') + cython(['geometry.pyx'], working_path=base_path) cython(['interpolation.pyx'], working_path=base_path) cython(['transform.pyx'], working_path=base_path) + config.add_extension('geometry', sources=['geometry.c']) config.add_extension('interpolation', sources=['interpolation.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('transform', sources=['transform.c'], diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx index 9dc1d307..fbd525fa 100644 --- a/skimage/draw/_draw.pyx +++ b/skimage/draw/_draw.pyx @@ -3,11 +3,7 @@ import math from libc.math cimport sqrt cimport numpy as np cimport cython - - -cdef extern from "../morphology/_pnpoly.h": - int pnpoly(int nr_verts, double *xp, double *yp, - double x, double y) +from skimage._shared.geometry cimport point_in_polygon @cython.boundscheck(False) @@ -119,7 +115,7 @@ def polygon(y, x, shape=None): for r in range(minr, maxr+1): for c in range(minc, maxc+1): - if pnpoly(nr_verts, cptr, rptr, c, r): + if point_in_polygon(nr_verts, cptr, rptr, c, r): rr.append(r) cc.append(c) diff --git a/skimage/draw/setup.py b/skimage/draw/setup.py index 4503d664..59067f9a 100644 --- a/skimage/draw/setup.py +++ b/skimage/draw/setup.py @@ -15,7 +15,7 @@ def configuration(parent_package='', top_path=None): cython(['_draw.pyx'], working_path=base_path) config.add_extension('_draw', sources=['_draw.c'], - include_dirs=[get_numpy_include_dirs()]) + include_dirs=[get_numpy_include_dirs(), '../shared']) return config diff --git a/skimage/morphology/_pnpoly.h b/skimage/morphology/_pnpoly.h deleted file mode 100644 index 95c89bcb..00000000 --- a/skimage/morphology/_pnpoly.h +++ /dev/null @@ -1,72 +0,0 @@ -/* `pnpoly` is from - - http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html - - Copyright (c) 1970-2003, Wm. Randolph Franklin - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - 1. Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimers. - 2. Redistributions in binary form must reproduce the above - copyright notice in the documentation and/or other materials - provided with the distribution. - 3. The name of W. Randolph Franklin may not be used to endorse or - promote products derived from this Software without specific - prior written permission. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS - BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN - ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN - CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. */ - -#ifdef __cplusplus -extern "C" { -#endif - -unsigned char pnpoly(int nr_verts, double *xp, double *yp, double x, double y) -{ - int i, j; - unsigned char c = 0; - for (i = 0, j = nr_verts-1; i < nr_verts; j = i++) { - if ((((yp[i]<=y) && (yvx.data, vy.data, m, n) + out[m, n] = point_in_polygon(V, vx.data, vy.data, m, n) return out.view(bool) - + def points_inside_poly(points, verts): """Test whether points lie inside a polygon. @@ -84,8 +77,8 @@ def points_inside_poly(points, verts): cdef np.ndarray[np.uint8_t, ndim=1] out = \ np.zeros(x.shape[0], dtype=np.uint8) - - npnpoly(vx.shape[0], vx.data, vy.data, + + points_in_polygon(vx.shape[0], vx.data, vy.data, x.shape[0], x.data, y.data, out.data) diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py index 2dd4a378..bbd176e4 100644 --- a/skimage/morphology/setup.py +++ b/skimage/morphology/setup.py @@ -28,7 +28,7 @@ def configuration(parent_package='', top_path=None): config.add_extension('_skeletonize_cy', sources=['_skeletonize_cy.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('_pnpoly', sources=['_pnpoly.c'], - include_dirs=[get_numpy_include_dirs()]) + include_dirs=[get_numpy_include_dirs(), '../shared']) config.add_extension('_convex_hull', sources=['_convex_hull.c'], include_dirs=[get_numpy_include_dirs()]) From 878554ac35ab5cb8013c7b0b02971a9381a03c79 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 21 Aug 2012 15:45:18 +0200 Subject: [PATCH 08/13] Fix indentation --- skimage/morphology/_pnpoly.pyx | 52 +++++++++++++++++----------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/skimage/morphology/_pnpoly.pyx b/skimage/morphology/_pnpoly.pyx index 9c897d48..7deb6a50 100644 --- a/skimage/morphology/_pnpoly.pyx +++ b/skimage/morphology/_pnpoly.pyx @@ -48,39 +48,39 @@ def grid_points_inside_poly(shape, verts): def points_inside_poly(points, verts): - """Test whether points lie inside a polygon. + """Test whether points lie inside a polygon. - Parameters - ---------- - points : (N, 2) array - Input points, ``(x, y)``. - verts : (M, 2) array - Vertices of the polygon, sorted either clockwise or anti-clockwise. - The first point may (but does not need to be) duplicated. + Parameters + ---------- + points : (N, 2) array + Input points, ``(x, y)``. + verts : (M, 2) array + Vertices of the polygon, sorted either clockwise or anti-clockwise. + The first point may (but does not need to be) duplicated. - Returns - ------- - mask : (N,) array of bool - True if corresponding point is inside the polygon. + Returns + ------- + mask : (N,) array of bool + True if corresponding point is inside the polygon. - """ - cdef np.ndarray[np.double_t, ndim=1, mode="c"] x, y, vx, vy + """ + cdef np.ndarray[np.double_t, ndim=1, mode="c"] x, y, vx, vy - points = np.asarray(points) - verts = np.asarray(verts) + points = np.asarray(points) + verts = np.asarray(verts) - x = points[:, 0].astype(np.double) - y = points[:, 1].astype(np.double) + x = points[:, 0].astype(np.double) + y = points[:, 1].astype(np.double) - vx = verts[:, 0].astype(np.double) - vy = verts[:, 1].astype(np.double) + vx = verts[:, 0].astype(np.double) + vy = verts[:, 1].astype(np.double) - cdef np.ndarray[np.uint8_t, ndim=1] out = \ - np.zeros(x.shape[0], dtype=np.uint8) + cdef np.ndarray[np.uint8_t, ndim=1] out = \ + np.zeros(x.shape[0], dtype=np.uint8) - points_in_polygon(vx.shape[0], vx.data, vy.data, - x.shape[0], x.data, y.data, - out.data) + points_in_polygon(vx.shape[0], vx.data, vy.data, + x.shape[0], x.data, y.data, + out.data) - return out.astype(bool) + return out.astype(bool) From 20afa7edae4dc0726424140625529e94e881e658 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 21 Aug 2012 15:49:18 +0200 Subject: [PATCH 09/13] Add doc string to point in polygon functions --- skimage/_shared/geometry.pyx | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/skimage/_shared/geometry.pyx b/skimage/_shared/geometry.pyx index de1ed4f4..e739f988 100644 --- a/skimage/_shared/geometry.pyx +++ b/skimage/_shared/geometry.pyx @@ -6,6 +6,17 @@ cdef inline unsigned char point_in_polygon(int nr_verts, double *xp, double *yp, double x, double y): + """Test whether point lies inside a polygon. + + Parameters + ---------- + nr_verts : int + Number of vertices of polygon. + xp, yp : double array + Coordinates of polygon with length nr_verts. + x, y : double + Coordinates of point. + """ cdef int i cdef unsigned char c = 0 cdef int j = nr_verts - 1 @@ -22,6 +33,21 @@ cdef inline unsigned char point_in_polygon(int nr_verts, double *xp, double *yp, cdef void points_in_polygon(int nr_verts, double *xp, double *yp, int nr_points, double *x, double *y, unsigned char *result): + """Test whether points lie inside a polygon. + + Parameters + ---------- + nr_verts : int + Number of vertices of polygon. + xp, yp : double array + Coordinates of polygon with length nr_verts. + nr_points : int + Number of points to test. + x, y : double array + Coordinates of points. + result : unsigned char array + Test results for each point. + """ cdef int n for n in range(nr_points): result[n] = point_in_polygon(nr_verts, xp, yp, x[n], y[n]) From dc91a0ee2f911c88c190a737726ead6608afec14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 21 Aug 2012 15:50:03 +0200 Subject: [PATCH 10/13] Add empty line between functions --- skimage/_shared/geometry.pyx | 1 + 1 file changed, 1 insertion(+) diff --git a/skimage/_shared/geometry.pyx b/skimage/_shared/geometry.pyx index e739f988..6f7de4cd 100644 --- a/skimage/_shared/geometry.pyx +++ b/skimage/_shared/geometry.pyx @@ -30,6 +30,7 @@ cdef inline unsigned char point_in_polygon(int nr_verts, double *xp, double *yp, j = i return c + cdef void points_in_polygon(int nr_verts, double *xp, double *yp, int nr_points, double *x, double *y, unsigned char *result): From 9731d4c21bfdd9c3578e6435543b22f528a982c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 23 Aug 2012 17:29:23 +0200 Subject: [PATCH 11/13] Simplify code by removing duplicate loops --- skimage/transform/_project.pyx | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/skimage/transform/_project.pyx b/skimage/transform/_project.pyx index 4f9a4704..dbbdfc7d 100644 --- a/skimage/transform/_project.pyx +++ b/skimage/transform/_project.pyx @@ -109,16 +109,13 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1, cdef int rows = img.shape[0] cdef int cols = img.shape[1] - if order == 0: - for tfr in range(out_r): - for tfc in range(out_c): - _matrix_transform(tfc, tfr, M.data, &c, &r) + for tfr in range(out_r): + for tfc in range(out_c): + _matrix_transform(tfc, tfr, M.data, &c, &r) + if order == 0: 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) + elif order == 1: out[tfr, tfc] = bilinear_interpolation(img.data, rows, cols, r, c, mode_c) From af708d2e00b8ec41681a593dc33905285100adb6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 25 Aug 2012 08:43:42 +0200 Subject: [PATCH 12/13] Fix typo in Cython compiler directives --- skimage/_shared/geometry.pyx | 2 +- skimage/_shared/interpolation.pyx | 2 +- skimage/_shared/transform.pyx | 2 +- skimage/feature/_texture.pyx | 2 +- skimage/transform/_project.pyx | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/skimage/_shared/geometry.pyx b/skimage/_shared/geometry.pyx index 6f7de4cd..3f4850b0 100644 --- a/skimage/_shared/geometry.pyx +++ b/skimage/_shared/geometry.pyx @@ -1,4 +1,4 @@ -#cython: cdivison=True +#cython: cdivision=True #cython: boundscheck=False #cython: nonecheck=False #cython: wraparound=False diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index 71852ace..137004e5 100644 --- a/skimage/_shared/interpolation.pyx +++ b/skimage/_shared/interpolation.pyx @@ -1,4 +1,4 @@ -#cython: cdivison=True +#cython: cdivision=True #cython: boundscheck=False #cython: nonecheck=False #cython: wraparound=False diff --git a/skimage/_shared/transform.pyx b/skimage/_shared/transform.pyx index e4bbfa74..ba0efc71 100644 --- a/skimage/_shared/transform.pyx +++ b/skimage/_shared/transform.pyx @@ -1,4 +1,4 @@ -#cython: cdivison=True +#cython: cdivision=True #cython: boundscheck=False #cython: nonecheck=False #cython: wraparound=False diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index 20b61513..0e82c306 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -1,4 +1,4 @@ -#cython: cdivison=True +#cython: cdivision=True #cython: boundscheck=False #cython: nonecheck=False #cython: wraparound=False diff --git a/skimage/transform/_project.pyx b/skimage/transform/_project.pyx index dbbdfc7d..3e9a9295 100644 --- a/skimage/transform/_project.pyx +++ b/skimage/transform/_project.pyx @@ -1,4 +1,4 @@ -#cython: cdivison=True +#cython: cdivision=True #cython: boundscheck=False #cython: nonecheck=False #cython: wraparound=False From 50d184374ad4260ce2c1cdde26aed5c4a030a03d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 25 Aug 2012 13:57:32 +0200 Subject: [PATCH 13/13] Optimize fast homography --- skimage/transform/_project.pyx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/skimage/transform/_project.pyx b/skimage/transform/_project.pyx index 3e9a9295..4df6e001 100644 --- a/skimage/transform/_project.pyx +++ b/skimage/transform/_project.pyx @@ -80,13 +80,15 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1, """ - 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"] img = \ + np.ascontiguousarray(image, dtype=np.double) cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] M = \ np.ascontiguousarray(np.linalg.inv(H)) if mode not in ('constant', 'wrap', 'mirror'): raise ValueError("Invalid mode specified. Please use " "`constant`, `wrap` or `mirror`.") + cdef char mode_c if mode == 'constant': mode_c = ord('C') elif mode == 'wrap': @@ -94,6 +96,7 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1, elif mode == 'mirror': mode_c = ord('M') + cdef int out_r, out_c if output_shape is None: out_r = img.shape[0] out_c = img.shape[1]