diff --git a/skimage/feature/_interest.pyx b/skimage/feature/_interest.pyx new file mode 100644 index 00000000..7af5aa7e --- /dev/null +++ b/skimage/feature/_interest.pyx @@ -0,0 +1,39 @@ +#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 87163247..cddd6039 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -2,7 +2,58 @@ import numpy as np from scipy import ndimage from skimage.color import rgb2grey from skimage.util import img_as_float -from . import peak +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 _corner._moravec(image, block_size) def harris(image, eps=1e-6, gaussian_deviation=1): diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 6f820163..905f5b99 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -12,9 +12,12 @@ 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(['_texture.pyx'], working_path=base_path) cython(['_template.pyx'], working_path=base_path) + config.add_extension('_interest', sources=['_interest.c'], + include_dirs=[get_numpy_include_dirs()]) config.add_extension('_texture', sources=['_texture.c'], include_dirs=[get_numpy_include_dirs(), '../_shared']) config.add_extension('_template', sources=['_template.c'],