From 2edb33f9500cfe25133b84ef7203b89d30003a3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 21:17:21 +0200 Subject: [PATCH] Make moravec standalone cython function and fix doc string --- skimage/feature/__init__.py | 1 + skimage/feature/_interest.pyx | 39 --------------- skimage/feature/interest.py | 55 +-------------------- skimage/feature/interest_cy.pyx | 85 +++++++++++++++++++++++++++++++++ skimage/feature/setup.py | 4 +- 5 files changed, 89 insertions(+), 95 deletions(-) delete mode 100644 skimage/feature/_interest.pyx create mode 100644 skimage/feature/interest_cy.pyx diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 4beca93e..233cb2f5 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -2,4 +2,5 @@ from ._hog import hog from .texture import greycomatrix, greycoprops, local_binary_pattern from .peak import peak_local_max from .interest import harris +from .interest_cy import moravec from .template import match_template diff --git a/skimage/feature/_interest.pyx b/skimage/feature/_interest.pyx deleted file mode 100644 index 7af5aa7e..00000000 --- a/skimage/feature/_interest.pyx +++ /dev/null @@ -1,39 +0,0 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False -import numpy as np -cimport numpy as cnp -from libc.float cimport DBL_MAX - - -def moravec(cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] image, - int block_size): - cdef int rows = image.shape[0] - cdef int cols = image.shape[1] - - cdef cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] out = \ - np.zeros_like(image) - - cdef double* image_data = image.data - cdef double* out_data = out.data - - cdef double msum, min_msum - cdef int r, c, br, bc, mr, mc, a, b - for r in range(2 * block_size, rows - 2 * block_size): - for c in range(2 * block_size, cols - 2 * block_size): - min_msum = DBL_MAX - for br in range(r - block_size, r + block_size + 1): - for bc in range(c - block_size, c + block_size + 1): - if br != r and bc != c: - msum = 0 - for mr in range(- block_size, block_size + 1): - for mc in range(- block_size, block_size + 1): - a = (r + mr) * cols + c + mc - b = (br + mr) * cols + bc + mc - msum += (image_data[a] - image_data[b]) ** 2 - min_msum = min(msum, min_msum) - - out_data[r * cols + c] = min_msum - - return out \ No newline at end of file diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index b3855c3a..1ce19ce6 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -1,59 +1,6 @@ import numpy as np from scipy import ndimage -from skimage.color import rgb2grey -from skimage.util import img_as_float -from . import peak, _interest - - -def moravec(image, block_size=3, mode='constant', cval=0): - """Compute Moravec response image. - - This interest operator is comparatively fast but not rotation invariant. - - Parameters - ---------- - image : ndarray - Input image. - block_size : int, optional - Block size for mean filtering the squared gradients. - mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional - The mode parameter determines how the array borders are handled, where - cval is the value when mode is equal to 'constant'. - cval : double, optional - Constant value to use for constant mode. - - Returns - ------- - coordinates : (N, 2) array - `(row, column)` coordinates of interest points. - - Examples - ------- - >>> from skimage.feature import moravec, peak_local_max - >>> square = np.zeros([10, 10]) - >>> square[1:9,1:9] = 1 - >>> square - array([[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.], - [ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.], - [ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.], - [ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.], - [ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.], - [ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]) - >>> moravec(square)), square.shape) - (2, 6) - - """ - - if image.ndim == 3: - image = rgb2grey(image) - - image = np.ascontiguousarray(img_as_float(image)) - - return _interest._moravec(image, block_size) +from . import peak def harris(image, eps=1e-6, gaussian_deviation=1): diff --git a/skimage/feature/interest_cy.pyx b/skimage/feature/interest_cy.pyx new file mode 100644 index 00000000..e3fa72db --- /dev/null +++ b/skimage/feature/interest_cy.pyx @@ -0,0 +1,85 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False +import numpy as np +cimport numpy as cnp +from libc.float cimport DBL_MAX + +from skimage.color import rgb2grey +from skimage.util import img_as_float + + +def moravec(image, int block_size=1): + """Compute Moravec response image. + + This interest operator is comparatively fast but not rotation invariant. + + Parameters + ---------- + image : ndarray + Input image. + block_size : int, optional + Block size for mean filtering the squared gradients. + + Returns + ------- + coordinates : (N, 2) array + `(row, column)` coordinates of interest points. + + Examples + ------- + >>> from skimage.feature import moravec, peak_local_max + >>> square = np.zeros([7, 7]) + >>> square[3, 3] = 1 + >>> square + array([[ 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 1., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0.]]) + >>> moravec(square) + array([[ 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 1., 1., 1., 0., 0.], + [ 0., 0., 1., 2., 1., 0., 0.], + [ 0., 0., 1., 1., 1., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0.]]) + """ + + cdef int rows = image.shape[0] + cdef int cols = image.shape[1] + + cdef cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] cimage, out + + if image.ndim == 3: + cimage = rgb2grey(image) + cimage = np.ascontiguousarray(img_as_float(image)) + + out = np.zeros_like(image) + + cdef double* image_data = cimage.data + cdef double* out_data = out.data + + cdef double msum, min_msum + cdef int r, c, br, bc, mr, mc, a, b + for r in range(2 * block_size, rows - 2 * block_size): + for c in range(2 * block_size, cols - 2 * block_size): + min_msum = DBL_MAX + for br in range(r - block_size, r + block_size + 1): + for bc in range(c - block_size, c + block_size + 1): + if br != r and bc != c: + msum = 0 + for mr in range(- block_size, block_size + 1): + for mc in range(- block_size, block_size + 1): + a = (r + mr) * cols + c + mc + b = (br + mr) * cols + bc + mc + msum += (image_data[a] - image_data[b]) ** 2 + min_msum = min(msum, min_msum) + + out_data[r * cols + c] = min_msum + + return out diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 905f5b99..20195500 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -12,11 +12,11 @@ def configuration(parent_package='', top_path=None): config = Configuration('feature', parent_package, top_path) config.add_data_dir('tests') - cython(['_interest.pyx'], working_path=base_path) + cython(['interest_cy.pyx'], working_path=base_path) cython(['_texture.pyx'], working_path=base_path) cython(['_template.pyx'], working_path=base_path) - config.add_extension('_interest', sources=['_interest.c'], + config.add_extension('interest_cy', sources=['interest_cy.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('_texture', sources=['_texture.c'], include_dirs=[get_numpy_include_dirs(), '../_shared'])