From 8c4ba00b12133494dc7bca45d5275ee9ca770342 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 20:59:19 +0200 Subject: [PATCH 01/37] Move harris file to a combined interest point file --- skimage/feature/__init__.py | 2 +- skimage/feature/{_harris.py => interest.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename skimage/feature/{_harris.py => interest.py} (100%) diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 1597e9a8..4beca93e 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,5 +1,5 @@ from ._hog import hog from .texture import greycomatrix, greycoprops, local_binary_pattern from .peak import peak_local_max -from ._harris import harris +from .interest import harris from .template import match_template diff --git a/skimage/feature/_harris.py b/skimage/feature/interest.py similarity index 100% rename from skimage/feature/_harris.py rename to skimage/feature/interest.py From 32a8442d613665f3c0cffdac0edc00fe2ee3c5db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 21:01:03 +0200 Subject: [PATCH 02/37] Refactor harris corner detection, so it returns harris response image --- skimage/feature/interest.py | 103 +++++++++++------------------------- 1 file changed, 32 insertions(+), 71 deletions(-) diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index ae30a29e..87163247 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -1,85 +1,31 @@ -""" -Harris corner detector - -Inspired from Solem's implementation -http://www.janeriksolem.net/2009/01/harris-corner-detector-in-python.html -""" +import numpy as np from scipy import ndimage - +from skimage.color import rgb2grey +from skimage.util import img_as_float from . import peak -def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1): - """Compute the Harris corner detector response function - for each pixel in the image +def harris(image, eps=1e-6, gaussian_deviation=1): + """Compute Harris response image. Parameters ---------- image : ndarray of floats Input image. - eps : float, optional Normalisation factor. - - gaussian_deviation : integer, optional - Standard deviation used for the Gaussian kernel. - - Returns - -------- - image : (M, N) ndarray - Harris image response - """ - if len(image.shape) == 3: - image = image.mean(axis=2) - - # derivatives - image = ndimage.gaussian_filter(image, gaussian_deviation) - imx = ndimage.sobel(image, axis=0, mode='constant') - imy = ndimage.sobel(image, axis=1, mode='constant') - - Wxx = ndimage.gaussian_filter(imx * imx, 1.5, mode='constant') - Wxy = ndimage.gaussian_filter(imx * imy, 1.5, mode='constant') - Wyy = ndimage.gaussian_filter(imy * imy, 1.5, mode='constant') - - # determinant and trace - Wdet = Wxx * Wyy - Wxy**2 - Wtr = Wxx + Wyy - # Alternate formula for Harris response. - # Alison Noble, "Descriptions of Image Surfaces", PhD thesis (1989) - harris = Wdet / (Wtr + eps) - - return harris - - -def harris(image, min_distance=10, threshold=0.1, eps=1e-6, - gaussian_deviation=1): - """Return corners from a Harris response image - - Parameters - ---------- - image : ndarray of floats - Input image. - - min_distance : int, optional - Minimum number of pixels separating interest points and image boundary. - - threshold : float, optional - Relative threshold impacting the number of interest points. - - eps : float, optional - Normalisation factor. - gaussian_deviation : integer, optional Standard deviation used for the Gaussian kernel. Returns ------- coordinates : (N, 2) array - (row, column) coordinates of interest points. + `(row, column)` coordinates of interest points. Examples ------- - >>> square = np.zeros([10,10]) + >>> from skimage.feature import harris, peak_local_max + >>> square = np.zeros([10, 10]) >>> square[2:8,2:8] = 1 >>> square array([[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], @@ -92,18 +38,33 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6, [ 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.]]) - >>> harris(square, min_distance=1) - - Corners of the square - + >>> peak_local_max(harris(square), min_distance=1) array([[3, 3], [3, 6], [6, 3], [6, 6]]) + """ - harrisim = _compute_harris_response(image, eps=eps, - gaussian_deviation=gaussian_deviation) - coordinates = peak.peak_local_max(harrisim, min_distance=min_distance, - threshold_rel=threshold) - return coordinates + if image.ndim == 3: + image = rgb2grey(image) + + # derivatives + image = ndimage.gaussian_filter(image, gaussian_deviation, mode='constant', + cval=0) + imx = ndimage.sobel(image, axis=0, mode='constant', cval=0) + imy = ndimage.sobel(image, axis=1, mode='constant', cval=0) + + Wxx = ndimage.gaussian_filter(imx * imx, 1.5, mode='constant', cval=0) + Wxy = ndimage.gaussian_filter(imx * imy, 1.5, mode='constant', cval=0) + Wyy = ndimage.gaussian_filter(imy * imy, 1.5, mode='constant', cval=0) + + # determinant and trace + Wdet = Wxx * Wyy - Wxy**2 + Wtr = Wxx + Wyy + + # Alternate formula for Harris response. + # Alison Noble, "Descriptions of Image Surfaces", PhD thesis (1989) + harris = Wdet / (Wtr + eps) + + return harris From 4991acf34ead450778303777bb93b35affb85b4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 21:02:45 +0200 Subject: [PATCH 03/37] Add Moravec interest point detection --- skimage/feature/_interest.pyx | 39 ++++++++++++++++++++++++++ skimage/feature/interest.py | 53 ++++++++++++++++++++++++++++++++++- skimage/feature/setup.py | 3 ++ 3 files changed, 94 insertions(+), 1 deletion(-) create mode 100644 skimage/feature/_interest.pyx 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'], From 1d96b29d775c3c11bde519aad44d70736da51218 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 21:03:46 +0200 Subject: [PATCH 04/37] Rename harris test case file --- skimage/feature/interest.py | 2 +- skimage/feature/tests/{test_harris.py => test_interest.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename skimage/feature/tests/{test_harris.py => test_interest.py} (100%) diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index cddd6039..b3855c3a 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -53,7 +53,7 @@ def moravec(image, block_size=3, mode='constant', cval=0): image = np.ascontiguousarray(img_as_float(image)) - return _corner._moravec(image, block_size) + return _interest._moravec(image, block_size) def harris(image, eps=1e-6, gaussian_deviation=1): diff --git a/skimage/feature/tests/test_harris.py b/skimage/feature/tests/test_interest.py similarity index 100% rename from skimage/feature/tests/test_harris.py rename to skimage/feature/tests/test_interest.py 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 05/37] 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']) From ecd2687226716074308c8ab155c97f0e2dde364c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 21:18:30 +0200 Subject: [PATCH 06/37] Fix return types in doc string of interest point detectors --- skimage/feature/interest.py | 6 +++--- skimage/feature/interest_cy.pyx | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index 1ce19ce6..d1a9a76c 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -8,7 +8,7 @@ def harris(image, eps=1e-6, gaussian_deviation=1): Parameters ---------- - image : ndarray of floats + image : ndarray Input image. eps : float, optional Normalisation factor. @@ -17,8 +17,8 @@ def harris(image, eps=1e-6, gaussian_deviation=1): Returns ------- - coordinates : (N, 2) array - `(row, column)` coordinates of interest points. + response : ndarray + Moravec response image. Examples ------- diff --git a/skimage/feature/interest_cy.pyx b/skimage/feature/interest_cy.pyx index e3fa72db..8d6867ed 100644 --- a/skimage/feature/interest_cy.pyx +++ b/skimage/feature/interest_cy.pyx @@ -24,8 +24,8 @@ def moravec(image, int block_size=1): Returns ------- - coordinates : (N, 2) array - `(row, column)` coordinates of interest points. + response : ndarray + Moravec response image. Examples ------- From f692055a25e3b746e8667cf8edddfd4cc658f541 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 21:19:46 +0200 Subject: [PATCH 07/37] Add reference to moravec operator --- skimage/feature/interest_cy.pyx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/skimage/feature/interest_cy.pyx b/skimage/feature/interest_cy.pyx index 8d6867ed..1974f8d7 100644 --- a/skimage/feature/interest_cy.pyx +++ b/skimage/feature/interest_cy.pyx @@ -27,6 +27,10 @@ def moravec(image, int block_size=1): response : ndarray Moravec response image. + References + ---------- + ..[1] http://kiwi.cs.dal.ca/~dparks/CornerDetection/moravec.htm + Examples ------- >>> from skimage.feature import moravec, peak_local_max From 3529afe0d2a7256cad89a25d09733ad9d59fad4a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 21:21:06 +0200 Subject: [PATCH 08/37] Rename window size parameter of moravec operator --- skimage/feature/interest_cy.pyx | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/skimage/feature/interest_cy.pyx b/skimage/feature/interest_cy.pyx index 1974f8d7..04a0e180 100644 --- a/skimage/feature/interest_cy.pyx +++ b/skimage/feature/interest_cy.pyx @@ -10,7 +10,7 @@ from skimage.color import rgb2grey from skimage.util import img_as_float -def moravec(image, int block_size=1): +def moravec(image, int window_size=1): """Compute Moravec response image. This interest operator is comparatively fast but not rotation invariant. @@ -19,8 +19,8 @@ def moravec(image, int block_size=1): ---------- image : ndarray Input image. - block_size : int, optional - Block size for mean filtering the squared gradients. + window_size : int, optional + Window size. Returns ------- @@ -70,15 +70,15 @@ def moravec(image, int block_size=1): 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): + for r in range(2 * window_size, rows - 2 * window_size): + for c in range(2 * window_size, cols - 2 * window_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): + for br in range(r - window_size, r + window_size + 1): + for bc in range(c - window_size, c + window_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): + for mr in range(- window_size, window_size + 1): + for mc in range(- window_size, window_size + 1): a = (r + mr) * cols + c + mc b = (br + mr) * cols + bc + mc msum += (image_data[a] - image_data[b]) ** 2 From 2ae23b427edfc0b1f3ddc1200b5256cca93ae0be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 21:47:32 +0200 Subject: [PATCH 09/37] Fix missing sigma for gaussian filter in harris operator --- skimage/feature/interest.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index d1a9a76c..15b249bc 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -3,7 +3,7 @@ from scipy import ndimage from . import peak -def harris(image, eps=1e-6, gaussian_deviation=1): +def harris(image, eps=1e-6, sigma=1): """Compute Harris response image. Parameters @@ -12,7 +12,7 @@ def harris(image, eps=1e-6, gaussian_deviation=1): Input image. eps : float, optional Normalisation factor. - gaussian_deviation : integer, optional + sigma : float, optional Standard deviation used for the Gaussian kernel. Returns @@ -48,14 +48,16 @@ def harris(image, eps=1e-6, gaussian_deviation=1): image = rgb2grey(image) # derivatives - image = ndimage.gaussian_filter(image, gaussian_deviation, mode='constant', - cval=0) + image = ndimage.gaussian_filter(image, sigma, mode='constant', cval=0) imx = ndimage.sobel(image, axis=0, mode='constant', cval=0) imy = ndimage.sobel(image, axis=1, mode='constant', cval=0) - Wxx = ndimage.gaussian_filter(imx * imx, 1.5, mode='constant', cval=0) - Wxy = ndimage.gaussian_filter(imx * imy, 1.5, mode='constant', cval=0) - Wyy = ndimage.gaussian_filter(imy * imy, 1.5, mode='constant', cval=0) + Wxx = ndimage.gaussian_filter(imx * imx, sigma, + mode='constant', cval=0) + Wxy = ndimage.gaussian_filter(imx * imy, sigma, + mode='constant', cval=0) + Wyy = ndimage.gaussian_filter(imy * imy, sigma, + mode='constant', cval=0) # determinant and trace Wdet = Wxx * Wyy - Wxy**2 From c4ca726566f391b010bbb91ed79c3bd4d56d99a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 22:20:35 +0200 Subject: [PATCH 10/37] Fix harris corner bug and add alternative corner measure param --- doc/examples/plot_harris.py | 6 ++--- skimage/feature/interest.py | 36 ++++++++++++++++---------- skimage/feature/tests/test_interest.py | 12 ++++----- 3 files changed, 31 insertions(+), 23 deletions(-) diff --git a/doc/examples/plot_harris.py b/doc/examples/plot_harris.py index 6212ac4a..39c1f29e 100644 --- a/doc/examples/plot_harris.py +++ b/doc/examples/plot_harris.py @@ -13,7 +13,7 @@ import numpy as np from matplotlib import pyplot as plt from skimage import data, img_as_float -from skimage.feature import harris +from skimage.feature import harris, peak_local_max def plot_harris_points(image, filtered_coords): @@ -29,12 +29,12 @@ plt.figure(figsize=(8, 3)) im_lena = img_as_float(data.lena()) im_text = img_as_float(data.text()) -filtered_coords = harris(im_lena, min_distance=4) +filtered_coords = peak_local_max(harris(im_lena), min_distance=4) plt.axes([0, 0, 0.3, 0.95]) plot_harris_points(im_lena, filtered_coords) -filtered_coords = harris(im_text, min_distance=4) +filtered_coords = peak_local_max(harris(im_text), min_distance=4) plt.axes([0.2, 0, 0.77, 1]) plot_harris_points(im_text, filtered_coords) diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index 15b249bc..88280a27 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -1,24 +1,31 @@ import numpy as np from scipy import ndimage +from skimage.color import rgb2grey from . import peak -def harris(image, eps=1e-6, sigma=1): +def harris(image, method='k', k=0.05, eps=1e-6, sigma=1): """Compute Harris response image. Parameters ---------- image : ndarray Input image. + method : {'k', 'eps'}, optional + Method to + k : float, optional + Sensitivity factor to separate corners from edges, typically in range + `[0, 0.2]`. Small values of k result in detection of sharp corners. eps : float, optional - Normalisation factor. + Normalisation factor (Noble's corner measure). sigma : float, optional - Standard deviation used for the Gaussian kernel. + Standard deviation used for the Gaussian kernel, which is used as + weighting function for the auto-correlation matrix. Returns ------- response : ndarray - Moravec response image. + Harris response image. Examples ------- @@ -48,23 +55,24 @@ def harris(image, eps=1e-6, sigma=1): image = rgb2grey(image) # derivatives - image = ndimage.gaussian_filter(image, sigma, mode='constant', cval=0) imx = ndimage.sobel(image, axis=0, mode='constant', cval=0) imy = ndimage.sobel(image, axis=1, mode='constant', cval=0) - Wxx = ndimage.gaussian_filter(imx * imx, sigma, + Axx = ndimage.gaussian_filter(imx * imx, sigma, mode='constant', cval=0) - Wxy = ndimage.gaussian_filter(imx * imy, sigma, + Axy = ndimage.gaussian_filter(imx * imy, sigma, mode='constant', cval=0) - Wyy = ndimage.gaussian_filter(imy * imy, sigma, + Ayy = ndimage.gaussian_filter(imy * imy, sigma, mode='constant', cval=0) - # determinant and trace - Wdet = Wxx * Wyy - Wxy**2 - Wtr = Wxx + Wyy + # determinant + detA = Axx * Ayy - Axy**2 + # trace + traceA = Axx + Ayy - # Alternate formula for Harris response. - # Alison Noble, "Descriptions of Image Surfaces", PhD thesis (1989) - harris = Wdet / (Wtr + eps) + if method == 'k': + harris = detA - k * traceA**2 + else: + harris = 2 * detA / (traceA + eps) return harris diff --git a/skimage/feature/tests/test_interest.py b/skimage/feature/tests/test_interest.py index 43bf28a3..22de33c3 100644 --- a/skimage/feature/tests/test_interest.py +++ b/skimage/feature/tests/test_interest.py @@ -3,13 +3,13 @@ import numpy as np from skimage import data from skimage import img_as_float -from skimage.feature import harris +from skimage.feature import moravec, harris, peak_local_max def test_square_image(): im = np.zeros((50, 50)).astype(float) im[:25, :25] = 1. - results = harris(im) + results = peak_local_max(harris(im)) assert results.any() assert len(results) == 1 @@ -18,7 +18,7 @@ def test_noisy_square_image(): im = np.zeros((50, 50)).astype(float) im[:25, :25] = 1. im = im + np.random.uniform(size=im.shape) * .5 - results = harris(im) + results = peak_local_max(harris(im)) assert results.any() assert len(results) == 1 @@ -27,7 +27,7 @@ def test_squared_dot(): im = np.zeros((50, 50)) im[4:8, 4:8] = 1 im = img_as_float(im) - results = harris(im, min_distance=3) + results = peak_local_max(harris(im)) assert (results == np.array([[6, 6]])).all() @@ -37,9 +37,9 @@ def test_rotated_lena(): rotation. """ im = img_as_float(data.lena().mean(axis=2)) - results = harris(im) + results = peak_local_max(harris(im)) im_rotated = im.T - results_rotated = harris(im_rotated) + results_rotated = peak_local_max(harris(im_rotated)) assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all() assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() From 43393860b2fad0c92420b6e68fb2b70183da6311 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 10 Sep 2012 22:23:54 +0200 Subject: [PATCH 11/37] Add comment to explain meaning of structure matrix of harris --- skimage/feature/interest.py | 1 + 1 file changed, 1 insertion(+) diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index 88280a27..cf9b3ed0 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -58,6 +58,7 @@ def harris(image, method='k', k=0.05, eps=1e-6, sigma=1): imx = ndimage.sobel(image, axis=0, mode='constant', cval=0) imy = ndimage.sobel(image, axis=1, mode='constant', cval=0) + # sum of squared differences / structure tensore Axx = ndimage.gaussian_filter(imx * imx, sigma, mode='constant', cval=0) Axy = ndimage.gaussian_filter(imx * imy, sigma, From 4a75c649e32a59dc15857329e3d6c464fcea07f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 08:19:04 +0200 Subject: [PATCH 12/37] Add Shi-Tomasi corner detector and fix bug in auto-correlation matrix --- skimage/feature/__init__.py | 2 +- skimage/feature/interest.py | 94 ++++++++++++++++++++++++++++++------- 2 files changed, 77 insertions(+), 19 deletions(-) diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 233cb2f5..ba1d8c76 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,6 +1,6 @@ from ._hog import hog from .texture import greycomatrix, greycoprops, local_binary_pattern from .peak import peak_local_max -from .interest import harris +from .interest import harris, shi_tomasi from .interest_cy import moravec from .template import match_template diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index cf9b3ed0..307f8ded 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -4,15 +4,57 @@ from skimage.color import rgb2grey from . import peak +def _compute_auto_correlation(image, sigma): + """Compute auto-correlation matrix using sum of squared differences. + + Parameters + ---------- + image : ndarray + Input image. + sigma : float + Standard deviation used for the Gaussian kernel, which is used as + weighting function for the auto-correlation matrix. + + Returns + ------- + Axx, Axy, Ayy : arrays + Elements of the auto-correlation matrix for each pixel in input image. + + """ + + if image.ndim == 3: + image = rgb2grey(image) + + # derivatives + gradient_weights = np.array([-1, 0, 1]) + imx = ndimage.convolve1d(image, gradient_weights, axis=0, + mode='constant', cval=0) + imy = ndimage.convolve1d(image, gradient_weights, axis=1, + mode='constant', cval=0) + + # structure tensore + Axx = ndimage.gaussian_filter(imx * imx, sigma, + mode='constant', cval=0) + Axy = ndimage.gaussian_filter(imx * imy, sigma, + mode='constant', cval=0) + Ayy = ndimage.gaussian_filter(imy * imy, sigma, + mode='constant', cval=0) + + return Axx, Axy, Ayy + + def harris(image, method='k', k=0.05, eps=1e-6, sigma=1): """Compute Harris response image. + This corner detector uses information in the auto-correlation matrix + (sum of squared differences) to make assumptions about the type of point. + Parameters ---------- image : ndarray Input image. method : {'k', 'eps'}, optional - Method to + Method to compute the response image from the auto-correlation matrix. k : float, optional Sensitivity factor to separate corners from edges, typically in range `[0, 0.2]`. Small values of k result in detection of sharp corners. @@ -51,20 +93,7 @@ def harris(image, method='k', k=0.05, eps=1e-6, sigma=1): """ - if image.ndim == 3: - image = rgb2grey(image) - - # derivatives - imx = ndimage.sobel(image, axis=0, mode='constant', cval=0) - imy = ndimage.sobel(image, axis=1, mode='constant', cval=0) - - # sum of squared differences / structure tensore - Axx = ndimage.gaussian_filter(imx * imx, sigma, - mode='constant', cval=0) - Axy = ndimage.gaussian_filter(imx * imy, sigma, - mode='constant', cval=0) - Ayy = ndimage.gaussian_filter(imy * imy, sigma, - mode='constant', cval=0) + Axx, Axy, Ayy = _compute_auto_correlation(image, sigma) # determinant detA = Axx * Ayy - Axy**2 @@ -72,8 +101,37 @@ def harris(image, method='k', k=0.05, eps=1e-6, sigma=1): traceA = Axx + Ayy if method == 'k': - harris = detA - k * traceA**2 + response = detA - k * traceA**2 else: - harris = 2 * detA / (traceA + eps) + response = 2 * detA / (traceA + eps) - return harris + return response + + +def shi_tomasi(image, sigma=1): + """Compute Shi-Tomasi (Kanade-Tomasi) response image. + + This corner detector uses information in the auto-correlation matrix + (sum of squared differences) to make assumptions about the type of point. + It is computationally more expensive than the harris corner detector as + it directly computes the minimum eigenvalue of the auto-correlation matrix. + + Parameters + ---------- + image : ndarray + Input image. + sigma : float, optional + Standard deviation used for the Gaussian kernel, which is used as + weighting function for the auto-correlation matrix. + + response : ndarray + Shi-Tomasi response image. + + """ + + Axx, Axy, Ayy = _compute_auto_correlation(image, sigma) + + # minimum eigenvalue of A + response = ((Axx + Ayy) - np.sqrt((Axx - Ayy)**2 + 4 * Axy**2)) / 2 + + return response From 82089f7fc8607f71dabb1cf17997181f2f5e280a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 08:37:25 +0200 Subject: [PATCH 13/37] Update test cases for shi-tomasi and moravec operator --- skimage/feature/tests/test_interest.py | 54 +++++++++++++++++++++++--- 1 file changed, 49 insertions(+), 5 deletions(-) diff --git a/skimage/feature/tests/test_interest.py b/skimage/feature/tests/test_interest.py index 22de33c3..2877e81d 100644 --- a/skimage/feature/tests/test_interest.py +++ b/skimage/feature/tests/test_interest.py @@ -3,23 +3,45 @@ import numpy as np from skimage import data from skimage import img_as_float -from skimage.feature import moravec, harris, peak_local_max +from skimage.feature import moravec, harris, shi_tomasi, peak_local_max def test_square_image(): im = np.zeros((50, 50)).astype(float) im[:25, :25] = 1. + + # Moravec + results = peak_local_max(moravec(im)) + # interest points along edge + assert len(results) == 57 + + # Harris results = peak_local_max(harris(im)) - assert results.any() + # interest at corner + assert len(results) == 1 + + # Shi-Tomasi + results = peak_local_max(shi_tomasi(im)) + # interest at corner assert len(results) == 1 def test_noisy_square_image(): im = np.zeros((50, 50)).astype(float) im[:25, :25] = 1. - im = im + np.random.uniform(size=im.shape) * .5 - results = peak_local_max(harris(im)) + im = im + np.random.uniform(size=im.shape) * .2 + + # Moravec + results = peak_local_max(moravec(im)) + # undefined number of interest points assert results.any() + + # Harris + results = peak_local_max(harris(im, sigma=1.5)) + assert len(results) == 1 + + # Shi-Tomasi + results = peak_local_max(shi_tomasi(im, sigma=1.5)) assert len(results) == 1 @@ -27,9 +49,17 @@ def test_squared_dot(): im = np.zeros((50, 50)) im[4:8, 4:8] = 1 im = img_as_float(im) + + # Moravec fails + + # Harris results = peak_local_max(harris(im)) assert (results == np.array([[6, 6]])).all() + # Shi-Tomasi + results = peak_local_max(shi_tomasi(im)) + assert (results == np.array([[6, 6]])).all() + def test_rotated_lena(): """ @@ -37,12 +67,26 @@ def test_rotated_lena(): rotation. """ im = img_as_float(data.lena().mean(axis=2)) - results = peak_local_max(harris(im)) im_rotated = im.T + + # Moravec + results = peak_local_max(moravec(im)) + results_rotated = peak_local_max(moravec(im_rotated)) + assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all() + assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() + + # Harris + results = peak_local_max(harris(im)) results_rotated = peak_local_max(harris(im_rotated)) assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all() assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() + # Shi-Tomasi + results = peak_local_max(shi_tomasi(im)) + results_rotated = peak_local_max(shi_tomasi(im_rotated)) + assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all() + assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() + if __name__ == '__main__': from numpy import testing From 6593e27b4d9ba305a8b0d4ef246fed64b54edaf0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 08:37:45 +0200 Subject: [PATCH 14/37] Fix contiguous bug in moravec operator --- skimage/feature/interest_cy.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/feature/interest_cy.pyx b/skimage/feature/interest_cy.pyx index 04a0e180..b14e3f6a 100644 --- a/skimage/feature/interest_cy.pyx +++ b/skimage/feature/interest_cy.pyx @@ -63,7 +63,7 @@ def moravec(image, int window_size=1): cimage = rgb2grey(image) cimage = np.ascontiguousarray(img_as_float(image)) - out = np.zeros_like(image) + out = np.zeros(image.shape, dtype=np.double) cdef double* image_data = cimage.data cdef double* out_data = out.data From 3508aa90321f2098433bb6f605310f8f35727606 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 08:43:01 +0200 Subject: [PATCH 15/37] Add example to doc string of shi-tomasi --- skimage/feature/interest.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index 307f8ded..1e2a86af 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -127,6 +127,28 @@ def shi_tomasi(image, sigma=1): response : ndarray Shi-Tomasi response image. + Examples + ------- + >>> from skimage.feature import shi_tomasi, peak_local_max + >>> square = np.zeros([10, 10]) + >>> square[2:8,2:8] = 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.]]) + >>> peak_local_max(shi_tomasi(square), min_distance=1) + array([[3, 3], + [3, 6], + [6, 3], + [6, 6]]) + """ Axx, Axy, Ayy = _compute_auto_correlation(image, sigma) From 410280f4a9c6cde82f515d9d54c7bab157a0c674 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 08:52:10 +0200 Subject: [PATCH 16/37] Add some more references to interest point functions --- skimage/feature/interest.py | 10 ++++++++++ skimage/feature/interest_cy.pyx | 1 + 2 files changed, 11 insertions(+) diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index 1e2a86af..ef47e7cf 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -69,6 +69,11 @@ def harris(image, method='k', k=0.05, eps=1e-6, sigma=1): response : ndarray Harris response image. + References + ---------- + ..[1] http://kiwi.cs.dal.ca/~dparks/CornerDetection/harris.htm + ..[2] http://en.wikipedia.org/wiki/Corner_detection + Examples ------- >>> from skimage.feature import harris, peak_local_max @@ -127,6 +132,11 @@ def shi_tomasi(image, sigma=1): response : ndarray Shi-Tomasi response image. + References + ---------- + ..[1] http://kiwi.cs.dal.ca/~dparks/CornerDetection/harris.htm + ..[2] http://en.wikipedia.org/wiki/Corner_detection + Examples ------- >>> from skimage.feature import shi_tomasi, peak_local_max diff --git a/skimage/feature/interest_cy.pyx b/skimage/feature/interest_cy.pyx index b14e3f6a..ae49faf1 100644 --- a/skimage/feature/interest_cy.pyx +++ b/skimage/feature/interest_cy.pyx @@ -30,6 +30,7 @@ def moravec(image, int window_size=1): References ---------- ..[1] http://kiwi.cs.dal.ca/~dparks/CornerDetection/moravec.htm + ..[2] http://en.wikipedia.org/wiki/Corner_detection Examples ------- From d55ef7ede8dbe0163de9bf7c8cf066b441b1b708 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 08:52:51 +0200 Subject: [PATCH 17/37] Add missing doc string return section --- skimage/feature/interest.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index ef47e7cf..f4075d57 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -129,6 +129,8 @@ def shi_tomasi(image, sigma=1): Standard deviation used for the Gaussian kernel, which is used as weighting function for the auto-correlation matrix. + Returns + ------- response : ndarray Shi-Tomasi response image. From 3b9ffe1cd89316f226dd437fdd2b3fa0246769dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 08:53:33 +0200 Subject: [PATCH 18/37] Improve code layout --- skimage/feature/interest.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py index f4075d57..cc63acae 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/interest.py @@ -33,12 +33,9 @@ def _compute_auto_correlation(image, sigma): mode='constant', cval=0) # structure tensore - Axx = ndimage.gaussian_filter(imx * imx, sigma, - mode='constant', cval=0) - Axy = ndimage.gaussian_filter(imx * imy, sigma, - mode='constant', cval=0) - Ayy = ndimage.gaussian_filter(imy * imy, sigma, - mode='constant', cval=0) + Axx = ndimage.gaussian_filter(imx * imx, sigma, mode='constant', cval=0) + Axy = ndimage.gaussian_filter(imx * imy, sigma, mode='constant', cval=0) + Ayy = ndimage.gaussian_filter(imy * imy, sigma, mode='constant', cval=0) return Axx, Axy, Ayy From 16f09d358dc7debc7b122bd8e8bc0ceeb68c24c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 18:58:45 +0200 Subject: [PATCH 19/37] Rename corner operators and file --- skimage/feature/__init__.py | 4 +-- skimage/feature/{interest.py => corner.py} | 4 +-- .../{interest_cy.pyx => corner_cy.pyx} | 2 +- skimage/feature/setup.py | 4 +-- .../{test_interest.py => test_corner.py} | 31 ++++++++++--------- 5 files changed, 23 insertions(+), 22 deletions(-) rename skimage/feature/{interest.py => corner.py} (98%) rename skimage/feature/{interest_cy.pyx => corner_cy.pyx} (98%) rename skimage/feature/tests/{test_interest.py => test_corner.py} (66%) diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index ba1d8c76..5b421c9c 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,6 +1,6 @@ from ._hog import hog from .texture import greycomatrix, greycoprops, local_binary_pattern from .peak import peak_local_max -from .interest import harris, shi_tomasi -from .interest_cy import moravec +from .corner import corner_harris, corner_shi_tomasi +from .corner_cy import corner_moravec from .template import match_template diff --git a/skimage/feature/interest.py b/skimage/feature/corner.py similarity index 98% rename from skimage/feature/interest.py rename to skimage/feature/corner.py index cc63acae..e8f6aa2a 100644 --- a/skimage/feature/interest.py +++ b/skimage/feature/corner.py @@ -40,7 +40,7 @@ def _compute_auto_correlation(image, sigma): return Axx, Axy, Ayy -def harris(image, method='k', k=0.05, eps=1e-6, sigma=1): +def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): """Compute Harris response image. This corner detector uses information in the auto-correlation matrix @@ -110,7 +110,7 @@ def harris(image, method='k', k=0.05, eps=1e-6, sigma=1): return response -def shi_tomasi(image, sigma=1): +def corner_shi_tomasi(image, sigma=1): """Compute Shi-Tomasi (Kanade-Tomasi) response image. This corner detector uses information in the auto-correlation matrix diff --git a/skimage/feature/interest_cy.pyx b/skimage/feature/corner_cy.pyx similarity index 98% rename from skimage/feature/interest_cy.pyx rename to skimage/feature/corner_cy.pyx index ae49faf1..6094d937 100644 --- a/skimage/feature/interest_cy.pyx +++ b/skimage/feature/corner_cy.pyx @@ -10,7 +10,7 @@ from skimage.color import rgb2grey from skimage.util import img_as_float -def moravec(image, int window_size=1): +def corner_moravec(image, int window_size=1): """Compute Moravec response image. This interest operator is comparatively fast but not rotation invariant. diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 20195500..e769621d 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_cy.pyx'], working_path=base_path) + cython(['corner_cy.pyx'], working_path=base_path) cython(['_texture.pyx'], working_path=base_path) cython(['_template.pyx'], working_path=base_path) - config.add_extension('interest_cy', sources=['interest_cy.c'], + config.add_extension('corner_cy', sources=['corner_cy.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('_texture', sources=['_texture.c'], include_dirs=[get_numpy_include_dirs(), '../_shared']) diff --git a/skimage/feature/tests/test_interest.py b/skimage/feature/tests/test_corner.py similarity index 66% rename from skimage/feature/tests/test_interest.py rename to skimage/feature/tests/test_corner.py index 2877e81d..2c75442c 100644 --- a/skimage/feature/tests/test_interest.py +++ b/skimage/feature/tests/test_corner.py @@ -3,7 +3,8 @@ import numpy as np from skimage import data from skimage import img_as_float -from skimage.feature import moravec, harris, shi_tomasi, peak_local_max +from skimage.feature import (corner_moravec, corner_harris, corner_shi_tomasi, + peak_local_max) def test_square_image(): @@ -11,17 +12,17 @@ def test_square_image(): im[:25, :25] = 1. # Moravec - results = peak_local_max(moravec(im)) + results = peak_local_max(corner_moravec(im)) # interest points along edge assert len(results) == 57 # Harris - results = peak_local_max(harris(im)) + results = peak_local_max(corner_harris(im)) # interest at corner assert len(results) == 1 # Shi-Tomasi - results = peak_local_max(shi_tomasi(im)) + results = peak_local_max(corner_shi_tomasi(im)) # interest at corner assert len(results) == 1 @@ -32,16 +33,16 @@ def test_noisy_square_image(): im = im + np.random.uniform(size=im.shape) * .2 # Moravec - results = peak_local_max(moravec(im)) + results = peak_local_max(corner_moravec(im)) # undefined number of interest points assert results.any() # Harris - results = peak_local_max(harris(im, sigma=1.5)) + results = peak_local_max(corner_harris(im, sigma=1.5)) assert len(results) == 1 # Shi-Tomasi - results = peak_local_max(shi_tomasi(im, sigma=1.5)) + results = peak_local_max(corner_shi_tomasi(im, sigma=1.5)) assert len(results) == 1 @@ -53,11 +54,11 @@ def test_squared_dot(): # Moravec fails # Harris - results = peak_local_max(harris(im)) + results = peak_local_max(corner_harris(im)) assert (results == np.array([[6, 6]])).all() # Shi-Tomasi - results = peak_local_max(shi_tomasi(im)) + results = peak_local_max(corner_shi_tomasi(im)) assert (results == np.array([[6, 6]])).all() @@ -70,20 +71,20 @@ def test_rotated_lena(): im_rotated = im.T # Moravec - results = peak_local_max(moravec(im)) - results_rotated = peak_local_max(moravec(im_rotated)) + results = peak_local_max(corner_moravec(im)) + results_rotated = peak_local_max(corner_moravec(im_rotated)) assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all() assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() # Harris - results = peak_local_max(harris(im)) - results_rotated = peak_local_max(harris(im_rotated)) + results = peak_local_max(corner_harris(im)) + results_rotated = peak_local_max(corner_harris(im_rotated)) assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all() assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() # Shi-Tomasi - results = peak_local_max(shi_tomasi(im)) - results_rotated = peak_local_max(shi_tomasi(im_rotated)) + results = peak_local_max(corner_shi_tomasi(im)) + results_rotated = peak_local_max(corner_shi_tomasi(im_rotated)) assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all() assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() From e26dc02ead617887b6a6c29011e7d0d4a95374b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 19:20:28 +0200 Subject: [PATCH 20/37] Add Kitchen and Rosenfeld corner detector --- skimage/feature/__init__.py | 2 +- skimage/feature/corner.py | 59 +++++++++++++++++++++++++++++++++---- 2 files changed, 54 insertions(+), 7 deletions(-) diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 5b421c9c..f74f31a8 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,6 +1,6 @@ from ._hog import hog from .texture import greycomatrix, greycoprops, local_binary_pattern from .peak import peak_local_max -from .corner import corner_harris, corner_shi_tomasi +from .corner import corner_kitchen_rosenfeld, corner_harris, corner_shi_tomasi from .corner_cy import corner_moravec from .template import match_template diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index e8f6aa2a..9ebee9f2 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -4,6 +4,30 @@ from skimage.color import rgb2grey from . import peak +def _compute_derivatives(image): + """Compute derivatives in x and y direction. + + Parameters + ---------- + image : ndarray + Input image. + + Returns + ------- + imx, imy : arrays + Derivatives in x and y direction. + + """ + + gradient_weights = np.array([-1, 0, 1]) + imx = ndimage.convolve1d(image, gradient_weights, axis=0, + mode='constant', cval=0) + imy = ndimage.convolve1d(image, gradient_weights, axis=1, + mode='constant', cval=0) + + return imx, imy + + def _compute_auto_correlation(image, sigma): """Compute auto-correlation matrix using sum of squared differences. @@ -25,12 +49,7 @@ def _compute_auto_correlation(image, sigma): if image.ndim == 3: image = rgb2grey(image) - # derivatives - gradient_weights = np.array([-1, 0, 1]) - imx = ndimage.convolve1d(image, gradient_weights, axis=0, - mode='constant', cval=0) - imy = ndimage.convolve1d(image, gradient_weights, axis=1, - mode='constant', cval=0) + imx, imy = _compute_derivatives(image) # structure tensore Axx = ndimage.gaussian_filter(imx * imx, sigma, mode='constant', cval=0) @@ -40,6 +59,34 @@ def _compute_auto_correlation(image, sigma): return Axx, Axy, Ayy +def corner_kitchen_rosenfeld(image): + """Compute Kitchen and Rosenfeld response image. + + This corner detector uses information in the auto-correlation matrix + (sum of squared differences) to make assumptions about the type of point. + + Parameters + ---------- + image : ndarray + Input image. + + Returns + ------- + response : ndarray + Kitchen and Rosenfeld response image. + + """ + + imx, imy = _compute_derivatives(image) + imxx, imxy = _compute_derivatives(imx) + imyx, imyy = _compute_derivatives(imy) + + response = (imxx * imy**2 + imyy * imx**2 - 2 * imxy * imx * imy) \ + / (imx**2 + imy**2) + + return response + + def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): """Compute Harris response image. From 4366462033fd797210213cf438dbba4205aaa148 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 19:23:15 +0200 Subject: [PATCH 21/37] Use sobel for image derivatives --- skimage/feature/corner.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index 9ebee9f2..c7e64233 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -5,7 +5,7 @@ from . import peak def _compute_derivatives(image): - """Compute derivatives in x and y direction. + """Compute derivatives in x and y direction using the Sobel operator. Parameters ---------- @@ -19,11 +19,8 @@ def _compute_derivatives(image): """ - gradient_weights = np.array([-1, 0, 1]) - imx = ndimage.convolve1d(image, gradient_weights, axis=0, - mode='constant', cval=0) - imy = ndimage.convolve1d(image, gradient_weights, axis=1, - mode='constant', cval=0) + imx = ndimage.sobel(image, axis=0, mode='constant', cval=0) + imy = ndimage.sobel(image, axis=1, mode='constant', cval=0) return imx, imy From 2d3cc8e0a047f23a747cbef991138a04db4b077c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 19:39:27 +0200 Subject: [PATCH 22/37] Add more detailed description of corner detectors in doc strings --- skimage/feature/corner.py | 37 +++++++++++++++++++++++++++-------- skimage/feature/corner_cy.pyx | 3 ++- 2 files changed, 31 insertions(+), 9 deletions(-) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index c7e64233..4302e8d0 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -59,8 +59,13 @@ def _compute_auto_correlation(image, sigma): def corner_kitchen_rosenfeld(image): """Compute Kitchen and Rosenfeld response image. - This corner detector uses information in the auto-correlation matrix - (sum of squared differences) to make assumptions about the type of point. + The corner measure is calculated as follows:: + + (imxx * imy**2 + imyy * imx**2 - 2 * imxy * imx * imy) + ------------------------------------------------------ + (imx**2 + imy**2) + + Where imx and imy are the first and imxx, imxy, imyy the second derivatives. Parameters ---------- @@ -87,8 +92,19 @@ def corner_kitchen_rosenfeld(image): def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): """Compute Harris response image. - This corner detector uses information in the auto-correlation matrix - (sum of squared differences) to make assumptions about the type of point. + This corner detector uses information from the auto-correlation matrix A:: + + A = [(imx**2) (imx*imy)] = [Axx Axy] + [(imx*imy) (imy**2)] [Axy Ayy] + + Where imx and imy are the first derivatives averaged with a gaussian filter. + The corner measure is then defined as:: + + det(A) - k * trace(A)**2 + + or:: + + 2 * det(A) / (trace(A) + eps) Parameters ---------- @@ -157,10 +173,15 @@ def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): def corner_shi_tomasi(image, sigma=1): """Compute Shi-Tomasi (Kanade-Tomasi) response image. - This corner detector uses information in the auto-correlation matrix - (sum of squared differences) to make assumptions about the type of point. - It is computationally more expensive than the harris corner detector as - it directly computes the minimum eigenvalue of the auto-correlation matrix. + This corner detector uses information from the auto-correlation matrix A:: + + A = [(imx**2) (imx*imy)] = [Axx Axy] + [(imx*imy) (imy**2)] [Axy Ayy] + + Where imx and imy are the first derivatives averaged with a gaussian filter. + The corner measure is then defined as the smaller eigenvalue of A:: + + ((Axx + Ayy) - sqrt((Axx - Ayy)**2 + 4 * Axy**2)) / 2 Parameters ---------- diff --git a/skimage/feature/corner_cy.pyx b/skimage/feature/corner_cy.pyx index 6094d937..94d69772 100644 --- a/skimage/feature/corner_cy.pyx +++ b/skimage/feature/corner_cy.pyx @@ -13,7 +13,8 @@ from skimage.util import img_as_float def corner_moravec(image, int window_size=1): """Compute Moravec response image. - This interest operator is comparatively fast but not rotation invariant. + This is one of the simplest corner detectors and is comparatively fast but + has several limitations (e.g. not rotation invariant). Parameters ---------- From 0f639e4989d5856d7aac18159e774e7c8571ce61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 20:06:50 +0200 Subject: [PATCH 23/37] Add foerstner corner detector --- skimage/feature/__init__.py | 3 +- skimage/feature/corner.py | 57 +++++++++++++++++++++++++++++++++++-- 2 files changed, 56 insertions(+), 4 deletions(-) diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index f74f31a8..c1e7f845 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,6 +1,7 @@ from ._hog import hog from .texture import greycomatrix, greycoprops, local_binary_pattern from .peak import peak_local_max -from .corner import corner_kitchen_rosenfeld, corner_harris, corner_shi_tomasi +from .corner import (corner_kitchen_rosenfeld, corner_harris, corner_shi_tomasi, + corner_foerstner) from .corner_cy import corner_moravec from .template import match_template diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index 4302e8d0..0d35c93a 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -57,7 +57,7 @@ def _compute_auto_correlation(image, sigma): def corner_kitchen_rosenfeld(image): - """Compute Kitchen and Rosenfeld response image. + """Compute Kitchen and Rosenfeld corner measure response image. The corner measure is calculated as follows:: @@ -90,7 +90,7 @@ def corner_kitchen_rosenfeld(image): def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): - """Compute Harris response image. + """Compute Harris corner measure response image. This corner detector uses information from the auto-correlation matrix A:: @@ -171,7 +171,7 @@ def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): def corner_shi_tomasi(image, sigma=1): - """Compute Shi-Tomasi (Kanade-Tomasi) response image. + """Compute Shi-Tomasi (Kanade-Tomasi) corner measure response image. This corner detector uses information from the auto-correlation matrix A:: @@ -231,3 +231,54 @@ def corner_shi_tomasi(image, sigma=1): response = ((Axx + Ayy) - np.sqrt((Axx - Ayy)**2 + 4 * Axy**2)) / 2 return response + + +def corner_foerstner(image, sigma=1): + """Compute Foerstner corner measure response image. + + This corner detector uses information from the auto-correlation matrix A:: + + A = [(imx**2) (imx*imy)] = [Axx Axy] + [(imx*imy) (imy**2)] [Axy Ayy] + + Where imx and imy are the first derivatives averaged with a gaussian filter. + The corner measure is then defined as:: + + w = det(A) / trace(A) (size of error ellipse) + q = 4 * det(A) / trace(A)**2 (roundness of error ellipse) + w * q (corner measure) + + Parameters + ---------- + image : ndarray + Input image. + sigma : float, optional + Standard deviation used for the Gaussian kernel, which is used as + weighting function for the auto-correlation matrix. + + Returns + ------- + response : ndarray + Foerstner response image. + + References + ---------- + ..[1] http://www.ipb.uni-bonn.de/uploads/tx_ikgpublication/\ + foerstner87.fast.pdf + ..[2] http://en.wikipedia.org/wiki/Corner_detection + + """ + + Axx, Axy, Ayy = _compute_auto_correlation(image, sigma) + + # determinant + detA = Axx * Ayy - Axy**2 + # trace + traceA = Axx + Ayy + + w = detA / traceA + q = 4 * detA / traceA**2 + + response = w * q + + return response From 6e193175f92857e1585cf068f7089471dc12e1f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 11 Sep 2012 20:19:24 +0200 Subject: [PATCH 24/37] Fix short description of moravec corner detector --- skimage/feature/corner_cy.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/feature/corner_cy.pyx b/skimage/feature/corner_cy.pyx index 94d69772..2f001ea4 100644 --- a/skimage/feature/corner_cy.pyx +++ b/skimage/feature/corner_cy.pyx @@ -11,7 +11,7 @@ from skimage.util import img_as_float def corner_moravec(image, int window_size=1): - """Compute Moravec response image. + """Compute Moravec corner measure response image. This is one of the simplest corner detectors and is comparatively fast but has several limitations (e.g. not rotation invariant). From bf3d34ccf27e251c2d36a5bf005218454add5bd5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 12 Sep 2012 08:24:12 +0200 Subject: [PATCH 25/37] Let foerstner return both corner measures --- skimage/feature/corner.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index 0d35c93a..4c1dd89c 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -246,7 +246,6 @@ def corner_foerstner(image, sigma=1): w = det(A) / trace(A) (size of error ellipse) q = 4 * det(A) / trace(A)**2 (roundness of error ellipse) - w * q (corner measure) Parameters ---------- @@ -258,8 +257,8 @@ def corner_foerstner(image, sigma=1): Returns ------- - response : ndarray - Foerstner response image. + w, q : ndarray + Foerstner response images. References ---------- @@ -279,6 +278,4 @@ def corner_foerstner(image, sigma=1): w = detA / traceA q = 4 * detA / traceA**2 - response = w * q - - return response + return w, q From cea889c99761ee612f9e9d1f3a79f6fb68fd74cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 16 Sep 2012 10:48:18 +0200 Subject: [PATCH 26/37] Add subpixel corner detection --- skimage/feature/__init__.py | 2 +- skimage/feature/corner.py | 128 +++++++++++++++++++++++++++++++++++- 2 files changed, 126 insertions(+), 4 deletions(-) diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index c1e7f845..ec81849f 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -2,6 +2,6 @@ from ._hog import hog from .texture import greycomatrix, greycoprops, local_binary_pattern from .peak import peak_local_max from .corner import (corner_kitchen_rosenfeld, corner_harris, corner_shi_tomasi, - corner_foerstner) + corner_foerstner, corner_subpix) from .corner_cy import corner_moravec from .template import match_template diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index 4c1dd89c..9c7a2fdf 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -1,6 +1,8 @@ import numpy as np from scipy import ndimage +from scipy import stats from skimage.color import rgb2grey +from skimage.util import img_as_float from . import peak @@ -19,8 +21,8 @@ def _compute_derivatives(image): """ - imx = ndimage.sobel(image, axis=0, mode='constant', cval=0) - imy = ndimage.sobel(image, axis=1, mode='constant', cval=0) + imy = ndimage.sobel(image, axis=0, mode='constant', cval=0) + imx = ndimage.sobel(image, axis=1, mode='constant', cval=0) return imx, imy @@ -44,7 +46,7 @@ def _compute_auto_correlation(image, sigma): """ if image.ndim == 3: - image = rgb2grey(image) + image = img_as_float(rgb2grey(image)) imx, imy = _compute_derivatives(image) @@ -279,3 +281,123 @@ def corner_foerstner(image, sigma=1): q = 4 * detA / traceA**2 return w, q + + +def corner_subpix(image, corners, window_size=11, alpha=0.99): + """Determine subpixel position of corners. + + Parameters + ---------- + image : ndarray + Input image. + corners : (N, 2) ndarray + Corner coordinates `(row, cols)`. + window_size : int, optional + Search window size for subpixel estimation. + alpha : float, optional + Significance level for point classification. + + Returns + ------- + positions : (N, 2) ndarray + Subpixel corner positions. NaN for "not classified" corners. + + """ + + # window extent in one direction + wext = (window_size - 1) / 2 + + # normal equation arrays + N_dot = np.zeros((2, 2), dtype=np.double) + N_edge = np.zeros((2, 2), dtype=np.double) + b_dot = np.zeros((2, ), dtype=np.double) + b_edge = np.zeros((2, ), dtype=np.double) + + # critical statistical test values + redundancy = window_size**2 - 2 + t_crit_dot = stats.f.isf(1 - alpha, redundancy, redundancy) + t_crit_edge = stats.f.isf(alpha, redundancy, redundancy) + + # coordinates of pixels within window + y, x = np.mgrid[- wext:wext + 1, - wext:wext + 1] + + # output arrays + corners_subpix = np.zeros_like(corners, dtype=np.double) + corners_class = np.zeros(corners.shape[0], dtype=np.int8) + + for i, (y0, x0) in enumerate(corners): + + # crop window around corner + border for sobel operator + miny = y0 - wext - 1 + maxy = y0 + wext + 2 + minx = x0 - wext - 1 + maxx = x0 + wext + 2 + window = image[miny:maxy, minx:maxx] + + winx, winy = _compute_derivatives(window) + + # compute gradient suares and remove border + winx_winx = (winx * winx)[1:-1, 1:-1] + winx_winy = (winx * winy)[1:-1, 1:-1] + winy_winy = (winy * winy)[1:-1, 1:-1] + + # sum of squared differences (mean instead of gaussian filter) + Axx = np.sum(winx_winx) + Axy = np.sum(winx_winy) + Ayy = np.sum(winy_winy) + + # sum of squared differences weighted with coordinates + # (mean instead of gaussian filter) + bxx_x = np.sum(winx_winx * x) + bxx_y = np.sum(winx_winx * y) + bxy_x = np.sum(winx_winy * x) + bxy_y = np.sum(winx_winy * y) + byy_x = np.sum(winy_winy * x) + byy_y = np.sum(winy_winy * y) + + # normal equations for subpixel position + N_dot[0, 0] = Axx + N_dot[0, 1] = N_dot[1, 0] = - Axy + N_dot[1, 1] = Ayy + + N_edge[0, 0] = Ayy + N_edge[0, 1] = N_edge[1, 0] = Axy + N_edge[1, 1] = Axx + + b_dot[:] = bxx_y - bxy_x, byy_x - bxy_y + b_edge[:] = byy_y + bxy_x, bxx_x + bxy_y + + # estimated positions + est_dot = np.linalg.solve(N_dot, b_dot) + est_edge = np.linalg.solve(N_edge, b_edge) + + # residuals + ry_dot = y - est_dot[0] + rx_dot = x - est_dot[1] + ry_edge = y - est_edge[0] + rx_edge = x - est_edge[1] + # squared residuals + rxx_dot = rx_dot * rx_dot + rxy_dot = rx_dot * ry_dot + ryy_dot = ry_dot * ry_dot + rxx_edge = rx_edge * rx_edge + rxy_edge = rx_edge * ry_edge + ryy_edge = ry_edge * ry_edge + + # determine corner class (dot or edge) + w_dot = np.sum(winx_winx * ryy_dot - 2 * winx_winy * rxy_dot \ + + winy_winy * rxx_dot) + w_edge = np.sum(winy_winy * ryy_edge + 2 * winx_winy * rxy_edge \ + + winx_winx * rxx_edge) + t = w_edge / w_dot + # 1 for edge, -1 for dot, 0 for "not classified" + corner_class = (t < t_crit_edge) - (t > t_crit_dot) + + if corner_class == - 1: + corners_subpix[i, :] = y0 + est_dot[0], x0 + est_dot[1] + elif corner_class == 0: + corners_subpix[i, :] = np.nan, np.nan + elif corner_class == 1: + corners_subpix[i, :] = y0 + est_edge[0], x0 + est_edge[1] + + return corners_subpix From 2273b73d99f39f010216726afe3b5ca9ebf62efe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 16 Sep 2012 11:04:19 +0200 Subject: [PATCH 27/37] Add references for subpixel corner detection --- skimage/feature/corner.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index 9c7a2fdf..aa5211bf 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -302,6 +302,12 @@ def corner_subpix(image, corners, window_size=11, alpha=0.99): positions : (N, 2) ndarray Subpixel corner positions. NaN for "not classified" corners. + References + ---------- + ..[1] http://www.ipb.uni-bonn.de/uploads/tx_ikgpublication/\ + foerstner87.fast.pdf + ..[2] http://en.wikipedia.org/wiki/Corner_detection + """ # window extent in one direction From 8af43f2d4940044cc98f206ca735a91b9343d9bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 16 Sep 2012 11:13:31 +0200 Subject: [PATCH 28/37] Add test case for subpixel corner detection --- skimage/feature/tests/test_corner.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py index 2c75442c..600d743e 100644 --- a/skimage/feature/tests/test_corner.py +++ b/skimage/feature/tests/test_corner.py @@ -1,10 +1,11 @@ import numpy as np +from numpy.testing import assert_array_equal from skimage import data from skimage import img_as_float from skimage.feature import (corner_moravec, corner_harris, corner_shi_tomasi, - peak_local_max) + corner_subpix, peak_local_max) def test_square_image(): @@ -89,6 +90,15 @@ def test_rotated_lena(): assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all() +def test_subpix(): + img = np.zeros((50, 50)) + img[:25,:25] = 255 + img[25:,25:] = 255 + corner = peak_local_max(corner_harris(img), num_peaks=1) + subpix = corner_subpix(img, corner) + assert_array_equal(subpix[0], (24.5, 24.5)) + + if __name__ == '__main__': from numpy import testing testing.run_module_suite() From c1f9336c2ac25bcd91f4d4f0ffad6b220f015192 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 16 Sep 2012 11:22:27 +0200 Subject: [PATCH 29/37] Update examples of corner detection in doc string --- skimage/feature/corner.py | 78 ++++++++++++++++++++++++++------------- 1 file changed, 52 insertions(+), 26 deletions(-) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index aa5211bf..69511dc7 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -135,21 +135,21 @@ def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): Examples ------- - >>> from skimage.feature import harris, peak_local_max + >>> from skimage.feature import corner_harris, peak_local_max >>> square = np.zeros([10, 10]) - >>> square[2:8,2:8] = 1 + >>> square[2:8, 2:8] = 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.]]) - >>> peak_local_max(harris(square), min_distance=1) + 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]]) + >>> peak_local_max(corner_harris(square), min_distance=1) array([[3, 3], [3, 6], [6, 3], @@ -205,21 +205,21 @@ def corner_shi_tomasi(image, sigma=1): Examples ------- - >>> from skimage.feature import shi_tomasi, peak_local_max + >>> from skimage.feature import corner_shi_tomasi, peak_local_max >>> square = np.zeros([10, 10]) - >>> square[2:8,2:8] = 1 + >>> square[2:8, 2:8] = 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.]]) - >>> peak_local_max(shi_tomasi(square), min_distance=1) + 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]]) + >>> peak_local_max(corner_shi_tomasi(square), min_distance=1) array([[3, 3], [3, 6], [6, 3], @@ -268,6 +268,32 @@ def corner_foerstner(image, sigma=1): foerstner87.fast.pdf ..[2] http://en.wikipedia.org/wiki/Corner_detection + Examples + ------- + >>> from skimage.feature import corner_foerstner, peak_local_max + >>> square = np.zeros([10, 10]) + >>> square[2:8, 2:8] = 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]]) + >>> w, q = corner_foerstner(square) + >>> accuracy_thresh = 0.5 + >>> roundness_thresh = 0.3 + >>> foerstner = (q > roundness_thresh) * (w > accuracy_thresh) * w + >>> peak_local_max(foerstner, min_distance=1) + array([[2, 2], + [2, 7], + [7, 2], + [7, 7]]) + """ Axx, Axy, Ayy = _compute_auto_correlation(image, sigma) From f95e93669658d586eccaf9dc00f144a8d48c036f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 16 Sep 2012 15:39:25 +0200 Subject: [PATCH 30/37] Refactor corner detection example --- doc/examples/plot_corner.py | 38 +++++++++++++++++++++++++++++++++ doc/examples/plot_harris.py | 42 ------------------------------------- 2 files changed, 38 insertions(+), 42 deletions(-) create mode 100644 doc/examples/plot_corner.py delete mode 100644 doc/examples/plot_harris.py diff --git a/doc/examples/plot_corner.py b/doc/examples/plot_corner.py new file mode 100644 index 00000000..7c9caa34 --- /dev/null +++ b/doc/examples/plot_corner.py @@ -0,0 +1,38 @@ +""" +================ +Corner detection +================ + +Detect corner points using the Harris corner detector and determine subpixel +position of corners. + +.. [1] http://en.wikipedia.org/wiki/Corner_detection +.. [2] http://en.wikipedia.org/wiki/Interest_point_detection + +""" + +import numpy as np +from matplotlib import pyplot as plt + +from skimage import data +from skimage.feature import corner_harris, corner_subpix, peak_local_max +from skimage.transform import warp, AffineTransform +from skimage.draw import ellipse + +tform = AffineTransform(scale=(1.3, 1.1), rotation=1, shear=0.7, + translation=(210, 50)) +image = warp(data.checkerboard(), tform.inverse, output_shape=(350, 350)) +rr, cc = ellipse(310, 175, 10, 100) +image[rr, cc] = 1 +image[180:230, 10:60] = 1 +image[230:280, 60:110] = 1 + +coords = peak_local_max(corner_harris(image), min_distance=5) +coords_subpix = corner_subpix(image, coords, window_size=13) + +plt.gray() +plt.imshow(image, interpolation='nearest') +plt.plot(coords[:, 1], coords[:, 0], '.b', markersize=3) +plt.plot(coords_subpix[:, 1], coords_subpix[:, 0], '+r', markersize=15) +plt.axis((0, 350, 350, 0)) +plt.show() diff --git a/doc/examples/plot_harris.py b/doc/examples/plot_harris.py deleted file mode 100644 index 39c1f29e..00000000 --- a/doc/examples/plot_harris.py +++ /dev/null @@ -1,42 +0,0 @@ -""" -=============================================================================== -Harris Corner detector -=============================================================================== - -The Harris corner filter [1]_ detects "interest points" [2]_ using edge -detection in multiple directions. - -.. [1] http://en.wikipedia.org/wiki/Corner_detection -.. [2] http://en.wikipedia.org/wiki/Interest_point_detection -""" -import numpy as np -from matplotlib import pyplot as plt - -from skimage import data, img_as_float -from skimage.feature import harris, peak_local_max - - -def plot_harris_points(image, filtered_coords): - """ plots corners found in image""" - - plt.imshow(image) - y, x = np.transpose(filtered_coords) - plt.plot(x, y, 'b.') - plt.axis('off') - -# display results -plt.figure(figsize=(8, 3)) -im_lena = img_as_float(data.lena()) -im_text = img_as_float(data.text()) - -filtered_coords = peak_local_max(harris(im_lena), min_distance=4) - -plt.axes([0, 0, 0.3, 0.95]) -plot_harris_points(im_lena, filtered_coords) - -filtered_coords = peak_local_max(harris(im_text), min_distance=4) - -plt.axes([0.2, 0, 0.77, 1]) -plot_harris_points(im_text, filtered_coords) - -plt.show() From 5b72c5d52c1ce3e4cd8f810a1d4adbb863fc6846 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 16 Sep 2012 20:10:00 +0200 Subject: [PATCH 31/37] Fix doc string for multiple return values --- skimage/feature/corner.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index 69511dc7..0363ee55 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -16,8 +16,10 @@ def _compute_derivatives(image): Returns ------- - imx, imy : arrays - Derivatives in x and y direction. + imx : ndarray + Derivative in x-direction. + imy : ndarray + Derivative in y-direction. """ @@ -40,8 +42,12 @@ def _compute_auto_correlation(image, sigma): Returns ------- - Axx, Axy, Ayy : arrays - Elements of the auto-correlation matrix for each pixel in input image. + Axx : ndarray + Element of the auto-correlation matrix for each pixel in input image. + Axy : ndarray + Element of the auto-correlation matrix for each pixel in input image. + Ayy : ndarray + Element of the auto-correlation matrix for each pixel in input image. """ @@ -259,8 +265,10 @@ def corner_foerstner(image, sigma=1): Returns ------- - w, q : ndarray - Foerstner response images. + w : ndarray + Error ellipse sizes. + q : ndarray + Roundness of error ellipse. References ---------- From 32cbf8e22609d97d801ab9967e307231107d2a68 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 16 Sep 2012 20:10:46 +0200 Subject: [PATCH 32/37] Fix typo in doc string --- skimage/feature/corner.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index 0363ee55..a62dcc73 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -325,7 +325,7 @@ def corner_subpix(image, corners, window_size=11, alpha=0.99): image : ndarray Input image. corners : (N, 2) ndarray - Corner coordinates `(row, cols)`. + Corner coordinates `(row, col)`. window_size : int, optional Search window size for subpixel estimation. alpha : float, optional From b22d570ff988a34c9b0008a2d9c9d26414b4b6fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 16 Sep 2012 23:18:16 +0200 Subject: [PATCH 33/37] Add comments and rename variables for better readability --- skimage/feature/corner.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index a62dcc73..64ccfb63 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -425,11 +425,13 @@ def corner_subpix(image, corners, window_size=11, alpha=0.99): ryy_edge = ry_edge * ry_edge # determine corner class (dot or edge) - w_dot = np.sum(winx_winx * ryy_dot - 2 * winx_winy * rxy_dot \ - + winy_winy * rxx_dot) - w_edge = np.sum(winy_winy * ryy_edge + 2 * winx_winy * rxy_edge \ - + winx_winx * rxx_edge) - t = w_edge / w_dot + # variance for different models + var_dot = np.sum(winx_winx * ryy_dot - 2 * winx_winy * rxy_dot \ + + winy_winy * rxx_dot) + var_edge = np.sum(winy_winy * ryy_edge + 2 * winx_winy * rxy_edge \ + + winx_winx * rxx_edge) + # test value (F-distributed) + t = var_edge / var_dot # 1 for edge, -1 for dot, 0 for "not classified" corner_class = (t < t_crit_edge) - (t > t_crit_dot) From 04fbee189e8b2eefaba1ea346561bb4601537911 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 16 Sep 2012 23:19:32 +0200 Subject: [PATCH 34/37] Remove unused array --- skimage/feature/corner.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index 64ccfb63..e48a0b41 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -361,9 +361,7 @@ def corner_subpix(image, corners, window_size=11, alpha=0.99): # coordinates of pixels within window y, x = np.mgrid[- wext:wext + 1, - wext:wext + 1] - # output arrays corners_subpix = np.zeros_like(corners, dtype=np.double) - corners_class = np.zeros(corners.shape[0], dtype=np.int8) for i, (y0, x0) in enumerate(corners): From 60c8e95905d11e1b2b1c852fb61d662826958f71 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 9 Dec 2012 18:14:36 +0100 Subject: [PATCH 35/37] Add function to detect corner peaks as a wrapper to peak_local_max --- doc/examples/plot_corner.py | 4 +- skimage/feature/__init__.py | 2 +- skimage/feature/corner.py | 87 ++++++++++++++++++++++++++++++------- 3 files changed, 75 insertions(+), 18 deletions(-) diff --git a/doc/examples/plot_corner.py b/doc/examples/plot_corner.py index 7c9caa34..4d01d177 100644 --- a/doc/examples/plot_corner.py +++ b/doc/examples/plot_corner.py @@ -15,7 +15,7 @@ import numpy as np from matplotlib import pyplot as plt from skimage import data -from skimage.feature import corner_harris, corner_subpix, peak_local_max +from skimage.feature import corner_harris, corner_subpix, corner_peaks from skimage.transform import warp, AffineTransform from skimage.draw import ellipse @@ -27,7 +27,7 @@ image[rr, cc] = 1 image[180:230, 10:60] = 1 image[230:280, 60:110] = 1 -coords = peak_local_max(corner_harris(image), min_distance=5) +coords = corner_peaks(corner_harris(image), min_distance=5) coords_subpix = corner_subpix(image, coords, window_size=13) plt.gray() diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index ec81849f..e0756bfe 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -2,6 +2,6 @@ from ._hog import hog from .texture import greycomatrix, greycoprops, local_binary_pattern from .peak import peak_local_max from .corner import (corner_kitchen_rosenfeld, corner_harris, corner_shi_tomasi, - corner_foerstner, corner_subpix) + corner_foerstner, corner_subpix, corner_peaks) from .corner_cy import corner_moravec from .template import match_template diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index e48a0b41..cbdbf881 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -3,7 +3,7 @@ from scipy import ndimage from scipy import stats from skimage.color import rgb2grey from skimage.util import img_as_float -from . import peak +from skimage.feature import peak_local_max def _compute_derivatives(image): @@ -141,7 +141,7 @@ def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): Examples ------- - >>> from skimage.feature import corner_harris, peak_local_max + >>> from skimage.feature import corner_harris, corner_peaks >>> square = np.zeros([10, 10]) >>> square[2:8, 2:8] = 1 >>> square @@ -155,11 +155,11 @@ def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1): [ 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]]) - >>> peak_local_max(corner_harris(square), min_distance=1) - array([[3, 3], - [3, 6], - [6, 3], - [6, 6]]) + >>> corner_peaks(corner_harris(square), min_distance=1) + array([[2, 2], + [2, 7], + [7, 2], + [7, 7]]) """ @@ -211,7 +211,7 @@ def corner_shi_tomasi(image, sigma=1): Examples ------- - >>> from skimage.feature import corner_shi_tomasi, peak_local_max + >>> from skimage.feature import corner_shi_tomasi, corner_peaks >>> square = np.zeros([10, 10]) >>> square[2:8, 2:8] = 1 >>> square @@ -225,11 +225,11 @@ def corner_shi_tomasi(image, sigma=1): [ 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]]) - >>> peak_local_max(corner_shi_tomasi(square), min_distance=1) - array([[3, 3], - [3, 6], - [6, 3], - [6, 6]]) + >>> corner_peaks(corner_shi_tomasi(square), min_distance=1) + array([[2, 2], + [2, 7], + [7, 2], + [7, 7]]) """ @@ -278,7 +278,7 @@ def corner_foerstner(image, sigma=1): Examples ------- - >>> from skimage.feature import corner_foerstner, peak_local_max + >>> from skimage.feature import corner_foerstner, corner_peaks >>> square = np.zeros([10, 10]) >>> square[2:8, 2:8] = 1 >>> square @@ -296,7 +296,7 @@ def corner_foerstner(image, sigma=1): >>> accuracy_thresh = 0.5 >>> roundness_thresh = 0.3 >>> foerstner = (q > roundness_thresh) * (w > accuracy_thresh) * w - >>> peak_local_max(foerstner, min_distance=1) + >>> corner_peaks(foerstner, min_distance=1) array([[2, 2], [2, 7], [7, 2], @@ -441,3 +441,60 @@ def corner_subpix(image, corners, window_size=11, alpha=0.99): corners_subpix[i, :] = y0 + est_edge[0], x0 + est_edge[1] return corners_subpix + + +def corner_peaks(image, min_distance=10, threshold_abs=0, threshold_rel=0.1, + exclude_border=True, indices=True, num_peaks=np.inf, + footprint=None, labels=None): + """Find corners in corner measure response image. + + This differs from `skimage.feature.peak_local_max` in that it suppresses + multiple connected peaks with the same accumulator value. + + Parameters + ---------- + See `skimage.feature.peak_local_max`. + + Returns + ------- + See `skimage.feature.peak_local_max`. + + Examples + -------- + >>> from skimage.feature import peak_local_max, corner_peaks + >>> response = np.zeros((5, 5)) + >>> response[2:4, 2:4] = 1 + >>> response + array([[ 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0.], + [ 0., 0., 1., 1., 0.], + [ 0., 0., 1., 1., 0.], + [ 0., 0., 0., 0., 0.]]) + >>> peak_local_max(response, exclude_border=False) + array([[2, 2], + [2, 3], + [3, 2], + [3, 3]]) + >>> corner_peaks(response, exclude_border=False) + array([[2, 2]]) + + """ + + peaks = peak_local_max(image, min_distance=min_distance, + threshold_abs=threshold_abs, + threshold_rel=threshold_rel, + exclude_border=exclude_border, + indices=False, num_peaks=np.inf, + footprint=footprint, labels=labels) + if min_distance > 0: + coords = np.transpose(peaks.nonzero()) + for r, c in coords: + if peaks[r, c]: + peaks[r - min_distance:r + min_distance + 1, + c - min_distance:c + min_distance + 1] = False + peaks[r, c] = True + + if indices is True: + return np.transpose(peaks.nonzero()) + else: + return peaks From a7411bc67853aa6573a06765210891f6025cede0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 9 Dec 2012 18:21:31 +0100 Subject: [PATCH 36/37] Extend example of corner_peaks to highlight difference to peak_local_max --- skimage/feature/corner.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index cbdbf881..867acb5d 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -477,6 +477,11 @@ def corner_peaks(image, min_distance=10, threshold_abs=0, threshold_rel=0.1, [3, 3]]) >>> corner_peaks(response, exclude_border=False) array([[2, 2]]) + >>> corner_peaks(response, exclude_border=False, min_distance=0) + array([[2, 2], + [2, 3], + [3, 2], + [3, 3]]) """ From 377ae3b968ef860ede7158919600d7cd614d1ff4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 9 Dec 2012 18:26:09 +0100 Subject: [PATCH 37/37] Add test case for corner_peaks function --- skimage/feature/tests/test_corner.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py index 600d743e..46e5466f 100644 --- a/skimage/feature/tests/test_corner.py +++ b/skimage/feature/tests/test_corner.py @@ -5,7 +5,7 @@ from skimage import data from skimage import img_as_float from skimage.feature import (corner_moravec, corner_harris, corner_shi_tomasi, - corner_subpix, peak_local_max) + corner_subpix, peak_local_max, corner_peaks) def test_square_image(): @@ -99,6 +99,17 @@ def test_subpix(): assert_array_equal(subpix[0], (24.5, 24.5)) +def test_corner_peaks(): + response = np.zeros((5, 5)) + response[2:4, 2:4] = 1 + + corners = corner_peaks(response, exclude_border=False) + assert len(corners) == 1 + + corners = corner_peaks(response, exclude_border=False, min_distance=0) + assert len(corners) == 4 + + if __name__ == '__main__': from numpy import testing testing.run_module_suite()