From c6f16dbff5c42ca35993c52f9380a9706b7b4a9b Mon Sep 17 00:00:00 2001 From: Nelle Varoquaux Date: Wed, 21 Dec 2011 19:09:40 +0100 Subject: [PATCH 1/7] Harris corner detector --- doc/examples/plot_harris.py | 30 +++++++++++ skimage/filter/__init__.py | 1 + skimage/filter/harris.py | 103 ++++++++++++++++++++++++++++++++++++ 3 files changed, 134 insertions(+) create mode 100644 doc/examples/plot_harris.py create mode 100644 skimage/filter/harris.py diff --git a/doc/examples/plot_harris.py b/doc/examples/plot_harris.py new file mode 100644 index 00000000..258925ed --- /dev/null +++ b/doc/examples/plot_harris.py @@ -0,0 +1,30 @@ +""" +=============================================================================== +Harris Corner detector +=============================================================================== + +The Harris corner filter detects interest points using edge detection in many +direction. +""" +from matplotlib import pyplot as plt +from matplotlib import cm + +from skimage import data +from skimage.filter import harris + + +def plot_harris_points(image, filtered_coords): + """ plots corners found in image""" + + plt.subplot(111) + plt.imshow(image, cmap=cm.gray) + plt.plot([p[1] for p in filtered_coords], + [p[0] for p in filtered_coords], + 'b.') + plt.axis('off') + plt.show() + + +im = data.lena().astype(float) +filtered_coords = harris.harris_corner_detector(im, 6) +plot_harris_points(im, filtered_coords) diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index 1acc33f2..6075291f 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -5,3 +5,4 @@ from edges import sobel, hsobel, vsobel, hprewitt, vprewitt, prewitt from tv_denoise import tv_denoise from rank_order import rank_order from thresholding import threshold_otsu +from harris import harris_corner_detector diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py new file mode 100644 index 00000000..53b58b01 --- /dev/null +++ b/skimage/filter/harris.py @@ -0,0 +1,103 @@ +# +# Harris detector +# +# http://www.janeriksolem.net/2009/01/harris-corner-detector-in-python.html + +import numpy as np +from scipy import ndimage + + +def _compute_harris_response(image, eps=1e-6): + """Compute the Harris corner detector response function + for each pixel in the image + + Params + ------- + image: ndarray + + eps: float, optional, default: 1e-6 + normalisation factor + + Returns + -------- + ndarray + """ + if len(image.shape) == 3: + image = image.mean(axis=2) + + # derivatives + image = ndimage.gaussian_filter(image, 1) + 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 + harris = Wdet / (Wtr + eps) + + # Non maximum filter of size 3 + harris_max = ndimage.maximum_filter(harris, 3, mode='constant') + harris *= harris == harris_max + # Remove the image corners + harris[:3] = 0 + harris[-3:] = 0 + harris[:, :3] = 0 + harris[:, -3:] = 0 + + return harris + + +def harris_corner_detector(image, min_distance=10, threshold=0.1, eps=1e-6): + """Return corners from a Harris response image + + params + ------- + harrisim: ndarray + + min_distance: int, optional, default: 10 + minimum number of pixels separating corners and image boundary + + threshold: float, optional, default: 0.1 + + eps: float, optional, default: 1e-6 + + returns: + -------- + array: coordinates + """ + harrisim = _compute_harris_response(image, eps=eps) + corner_threshold = np.max(harrisim.ravel()) * threshold + # find top corner candidates above a threshold + # corner_threshold = max(harrisim.ravel()) * threshold + harrisim_t = (harrisim >= corner_threshold) * 1 + + # get coordinates of candidates + candidates = harrisim_t.nonzero() + coords = [(candidates[0][c], candidates[1][c]) for c + in range(len(candidates[0]))] + + # ...and their values + candidate_values = [harrisim[c[0]][c[1]] for c in coords] + + # sort candidates + index = np.argsort(candidate_values) + + # store allowed point locations in array + allowed_locations = np.zeros(harrisim.shape) + allowed_locations[min_distance:-min_distance, + min_distance:-min_distance] = 1 + + # select the best points taking min_distance into account + filtered_coords = [] + for i in index: + if allowed_locations[coords[i][0]][coords[i][1]] == 1: + filtered_coords.append(coords[i]) + allowed_locations[ + (coords[i][0] - min_distance):(coords[i][0] + min_distance), + (coords[i][1] - min_distance):(coords[i][1] + min_distance)] = 0 + + return np.array(filtered_coords) From 5d8110bab201b80191de8907ddfe78ff7dfb68dd Mon Sep 17 00:00:00 2001 From: Nelle Varoquaux Date: Thu, 22 Dec 2011 12:27:06 +0100 Subject: [PATCH 2/7] Added tests --- doc/examples/plot_harris.py | 5 +++-- skimage/filter/harris.py | 3 ++- skimage/filter/tests/test_harris.py | 23 +++++++++++++++++++++++ 3 files changed, 28 insertions(+), 3 deletions(-) create mode 100644 skimage/filter/tests/test_harris.py diff --git a/doc/examples/plot_harris.py b/doc/examples/plot_harris.py index 258925ed..6cb02ba1 100644 --- a/doc/examples/plot_harris.py +++ b/doc/examples/plot_harris.py @@ -6,11 +6,12 @@ Harris Corner detector The Harris corner filter detects interest points using edge detection in many direction. """ + from matplotlib import pyplot as plt from matplotlib import cm from skimage import data -from skimage.filter import harris +from skimage.filter import harris_corner_detector def plot_harris_points(image, filtered_coords): @@ -26,5 +27,5 @@ def plot_harris_points(image, filtered_coords): im = data.lena().astype(float) -filtered_coords = harris.harris_corner_detector(im, 6) +filtered_coords = harris_corner_detector(im, 6) plot_harris_points(im, filtered_coords) diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py index 53b58b01..b5abe7a3 100644 --- a/skimage/filter/harris.py +++ b/skimage/filter/harris.py @@ -1,6 +1,7 @@ # # Harris detector # +# Inspired from Solem's implementation # http://www.janeriksolem.net/2009/01/harris-corner-detector-in-python.html import numpy as np @@ -56,7 +57,7 @@ def harris_corner_detector(image, min_distance=10, threshold=0.1, eps=1e-6): params ------- - harrisim: ndarray + harrisim: ndarray of floats min_distance: int, optional, default: 10 minimum number of pixels separating corners and image boundary diff --git a/skimage/filter/tests/test_harris.py b/skimage/filter/tests/test_harris.py new file mode 100644 index 00000000..19d8cd3b --- /dev/null +++ b/skimage/filter/tests/test_harris.py @@ -0,0 +1,23 @@ +import numpy as np + +import unittest + +from skimage.filter import harris_corner_detector + + +class TestHarris(unittest.TestCase): + + def test_square_image(self): + im = np.zeros((50, 50)).astype(float) + im[:25, :25] = 1. + results = harris_corner_detector(im) + self.assertTrue(results.any()) + self.assertTrue(len(results) == 1) + + def test_noisy_square_image(self): + im = np.zeros((50, 50)).astype(float) + im[:25, :25] = 1. + im = im + np.random.uniform(size=im.shape) * .5 + results = harris_corner_detector(im) + assert results.any() + assert len(results) == 1 From 7ff98597e2cf1a92e4a1b957926a43fdfbb29e1e Mon Sep 17 00:00:00 2001 From: Nelle Varoquaux Date: Thu, 22 Dec 2011 22:08:27 +0100 Subject: [PATCH 3/7] Small fixes on the Harris corner detector - documentation, method renaming etc --- CONTRIBUTORS.txt | 1 + doc/examples/plot_harris.py | 21 ++++++------- skimage/filter/__init__.py | 2 +- skimage/filter/harris.py | 48 +++++++++++++++++------------ skimage/filter/tests/test_harris.py | 23 ++++++++------ 5 files changed, 54 insertions(+), 41 deletions(-) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 007f267a..2e97d7e1 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -83,6 +83,7 @@ - Nelle Varoquaux Renaming of the package to ``skimage``. + Harris corner detector - W. Randolph Franklin Point in polygon test. diff --git a/doc/examples/plot_harris.py b/doc/examples/plot_harris.py index 6cb02ba1..9c499a3b 100644 --- a/doc/examples/plot_harris.py +++ b/doc/examples/plot_harris.py @@ -3,29 +3,28 @@ Harris Corner detector =============================================================================== -The Harris corner filter detects interest points using edge detection in many -direction. +The Harris corner filter detects interest points using edge detection in +multiple direction. """ from matplotlib import pyplot as plt -from matplotlib import cm -from skimage import data -from skimage.filter import harris_corner_detector +from skimage import data, img_as_float +from skimage.filter import harris def plot_harris_points(image, filtered_coords): """ plots corners found in image""" - plt.subplot(111) - plt.imshow(image, cmap=cm.gray) + plt.plot() + plt.imshow(image) plt.plot([p[1] for p in filtered_coords], - [p[0] for p in filtered_coords], - 'b.') + [p[0] for p in filtered_coords], + 'b.') plt.axis('off') plt.show() -im = data.lena().astype(float) -filtered_coords = harris_corner_detector(im, 6) +im = img_as_float(data.lena()) +filtered_coords = harris(im, 6) plot_harris_points(im, filtered_coords) diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index 6075291f..529331f0 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -5,4 +5,4 @@ from edges import sobel, hsobel, vsobel, hprewitt, vprewitt, prewitt from tv_denoise import tv_denoise from rank_order import rank_order from thresholding import threshold_otsu -from harris import harris_corner_detector +from harris import harris diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py index b5abe7a3..303431dd 100644 --- a/skimage/filter/harris.py +++ b/skimage/filter/harris.py @@ -10,18 +10,20 @@ from scipy import ndimage def _compute_harris_response(image, eps=1e-6): """Compute the Harris corner detector response function - for each pixel in the image + for each pixel in the image - Params - ------- - image: ndarray + Parameters + ---------- + image: ndarray of floats + input image - eps: float, optional, default: 1e-6 - normalisation factor + eps: float, optional + normalisation factor Returns -------- - ndarray + features: (M, 2) ndarray + Harris image response """ if len(image.shape) == 3: image = image.mean(axis=2) @@ -42,8 +44,10 @@ def _compute_harris_response(image, eps=1e-6): # Non maximum filter of size 3 harris_max = ndimage.maximum_filter(harris, 3, mode='constant') - harris *= harris == harris_max - # Remove the image corners + mask = (harris == harris_max) + harris *= mask + + # Remove the image borders harris[:3] = 0 harris[-3:] = 0 harris[:, :3] = 0 @@ -52,23 +56,26 @@ def _compute_harris_response(image, eps=1e-6): return harris -def harris_corner_detector(image, min_distance=10, threshold=0.1, eps=1e-6): +def harris(image, min_distance=10, threshold=0.1, eps=1e-6): """Return corners from a Harris response image - params - ------- - harrisim: ndarray of floats + Parameters + ---------- + image: ndarray of floats + Input image - min_distance: int, optional, default: 10 - minimum number of pixels separating corners and image boundary + min_distance: int, optional + minimum number of pixels separating interest points and image boundary - threshold: float, optional, default: 0.1 + threshold: float, optional + relative threshold impacting the number of interest points. - eps: float, optional, default: 1e-6 + eps: float, optional + Normalisation factor returns: -------- - array: coordinates + array: coordinates of interest points """ harrisim = _compute_harris_response(image, eps=eps) corner_threshold = np.max(harrisim.ravel()) * threshold @@ -78,8 +85,9 @@ def harris_corner_detector(image, min_distance=10, threshold=0.1, eps=1e-6): # get coordinates of candidates candidates = harrisim_t.nonzero() - coords = [(candidates[0][c], candidates[1][c]) for c - in range(len(candidates[0]))] + coords = np.concatenate((candidates[0].reshape((len(candidates[0]), 1)), + candidates[1].reshape((len(candidates[0]), 1))), + axis=1) # ...and their values candidate_values = [harrisim[c[0]][c[1]] for c in coords] diff --git a/skimage/filter/tests/test_harris.py b/skimage/filter/tests/test_harris.py index 19d8cd3b..c78d3e26 100644 --- a/skimage/filter/tests/test_harris.py +++ b/skimage/filter/tests/test_harris.py @@ -1,23 +1,28 @@ import numpy as np -import unittest - -from skimage.filter import harris_corner_detector +from skimage.filter import harris +from skimage import img_as_float -class TestHarris(unittest.TestCase): - +class TestHarris(): def test_square_image(self): im = np.zeros((50, 50)).astype(float) im[:25, :25] = 1. - results = harris_corner_detector(im) - self.assertTrue(results.any()) - self.assertTrue(len(results) == 1) + results = harris(im) + assert results.any() + assert len(results) == 1 def test_noisy_square_image(self): im = np.zeros((50, 50)).astype(float) im[:25, :25] = 1. im = im + np.random.uniform(size=im.shape) * .5 - results = harris_corner_detector(im) + results = harris(im) assert results.any() assert len(results) == 1 + + def test_squared_dot(self): + im = np.zeros((50, 50)) + im[4:8, 4:8] = 1 + im = img_as_float(im) + results = harris(im) + assert results == np.array([6, 6]) From bc0f77622e28e4c9eb27a2c0bb58f8261dd3c8e3 Mon Sep 17 00:00:00 2001 From: Nelle Varoquaux Date: Fri, 23 Dec 2011 17:21:55 +0100 Subject: [PATCH 4/7] Added tests for harris, and small fixes --- skimage/filter/harris.py | 10 ++++------ skimage/filter/tests/test_harris.py | 18 ++++++++++++++++-- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py index 303431dd..2d48f7ec 100644 --- a/skimage/filter/harris.py +++ b/skimage/filter/harris.py @@ -22,7 +22,7 @@ def _compute_harris_response(image, eps=1e-6): Returns -------- - features: (M, 2) ndarray + image: (M, N) ndarray Harris image response """ if len(image.shape) == 3: @@ -85,12 +85,10 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6): # get coordinates of candidates candidates = harrisim_t.nonzero() - coords = np.concatenate((candidates[0].reshape((len(candidates[0]), 1)), - candidates[1].reshape((len(candidates[0]), 1))), - axis=1) + coords = np.transpose(candidates) # ...and their values - candidate_values = [harrisim[c[0]][c[1]] for c in coords] + candidate_values = harrisim[candidates] # sort candidates index = np.argsort(candidate_values) @@ -103,7 +101,7 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6): # select the best points taking min_distance into account filtered_coords = [] for i in index: - if allowed_locations[coords[i][0]][coords[i][1]] == 1: + if allowed_locations[tuple(coords[i])] == 1: filtered_coords.append(coords[i]) allowed_locations[ (coords[i][0] - min_distance):(coords[i][0] + min_distance), diff --git a/skimage/filter/tests/test_harris.py b/skimage/filter/tests/test_harris.py index c78d3e26..853a00a0 100644 --- a/skimage/filter/tests/test_harris.py +++ b/skimage/filter/tests/test_harris.py @@ -1,8 +1,10 @@ import numpy as np -from skimage.filter import harris +from skimage import data from skimage import img_as_float +from skimage.filter import harris + class TestHarris(): def test_square_image(self): @@ -24,5 +26,17 @@ class TestHarris(): im = np.zeros((50, 50)) im[4:8, 4:8] = 1 im = img_as_float(im) + results = harris(im, min_distance=3) + print results + assert (results == np.array([[6, 6]])).all() + + def test_rotated_lena(self): + """ + The harris filter should yield the same results with an image and it's + rotation. + """ + im = img_as_float(data.lena().mean(axis=2)) results = harris(im) - assert results == np.array([6, 6]) + im_rotated = im.T + results_rotated = harris(im_rotated) + assert (results[:, 0] == results_rotated[:, 1]).all() From 374c44671517952d3cf0f97c588fddcac05db70a Mon Sep 17 00:00:00 2001 From: Nelle Varoquaux Date: Fri, 23 Dec 2011 17:28:57 +0100 Subject: [PATCH 5/7] Added parameter to Harris corner detector to set the deviation for the gaussian kernel of the harris response computation --- skimage/filter/harris.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py index 2d48f7ec..8e8a46ec 100644 --- a/skimage/filter/harris.py +++ b/skimage/filter/harris.py @@ -8,7 +8,7 @@ import numpy as np from scipy import ndimage -def _compute_harris_response(image, eps=1e-6): +def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1): """Compute the Harris corner detector response function for each pixel in the image @@ -20,6 +20,9 @@ def _compute_harris_response(image, eps=1e-6): eps: float, optional normalisation factor + gaussian_deviation: integer, optional + standard deviation used for the Gaussian kernel + Returns -------- image: (M, N) ndarray @@ -29,7 +32,7 @@ def _compute_harris_response(image, eps=1e-6): image = image.mean(axis=2) # derivatives - image = ndimage.gaussian_filter(image, 1) + image = ndimage.gaussian_filter(image, gaussian_deviation) imx = ndimage.sobel(image, axis=0, mode='constant') imy = ndimage.sobel(image, axis=1, mode='constant') @@ -56,7 +59,8 @@ def _compute_harris_response(image, eps=1e-6): return harris -def harris(image, min_distance=10, threshold=0.1, eps=1e-6): +def harris(image, min_distance=10, threshold=0.1, eps=1e-6, + gaussian_deviation=1): """Return corners from a Harris response image Parameters @@ -73,11 +77,15 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6): eps: float, optional Normalisation factor + gaussian_deviation: integer, optional + standard deviation used for the Gaussian kernel + returns: -------- array: coordinates of interest points """ - harrisim = _compute_harris_response(image, eps=eps) + harrisim = _compute_harris_response(image, eps=eps, + gaussian_deviation=gaussian_deviation) corner_threshold = np.max(harrisim.ravel()) * threshold # find top corner candidates above a threshold # corner_threshold = max(harrisim.ravel()) * threshold From 3c2988b4bc3a7ad83f730c8c78ae5b0e7dd02541 Mon Sep 17 00:00:00 2001 From: Nelle Varoquaux Date: Sun, 25 Dec 2011 09:54:47 +0100 Subject: [PATCH 6/7] small fixes on the Harris' documentation: added capital letters where needed --- skimage/filter/harris.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py index 8e8a46ec..7eeecd30 100644 --- a/skimage/filter/harris.py +++ b/skimage/filter/harris.py @@ -15,13 +15,13 @@ def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1): Parameters ---------- image: ndarray of floats - input image + Input image eps: float, optional - normalisation factor + Normalisation factor gaussian_deviation: integer, optional - standard deviation used for the Gaussian kernel + Standard deviation used for the Gaussian kernel Returns -------- @@ -69,16 +69,16 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6, Input image min_distance: int, optional - minimum number of pixels separating interest points and image boundary + Minimum number of pixels separating interest points and image boundary threshold: float, optional - relative threshold impacting the number of interest points. + Relative threshold impacting the number of interest points. eps: float, optional Normalisation factor gaussian_deviation: integer, optional - standard deviation used for the Gaussian kernel + Standard deviation used for the Gaussian kernel returns: -------- @@ -86,9 +86,9 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6, """ harrisim = _compute_harris_response(image, eps=eps, gaussian_deviation=gaussian_deviation) - corner_threshold = np.max(harrisim.ravel()) * threshold + # find top corner candidates above a threshold - # corner_threshold = max(harrisim.ravel()) * threshold + corner_threshold = np.max(harrisim.ravel()) * threshold harrisim_t = (harrisim >= corner_threshold) * 1 # get coordinates of candidates From 7285cbe29604a5a0d362348ad823d04578039a2a Mon Sep 17 00:00:00 2001 From: Nelle Varoquaux Date: Thu, 5 Jan 2012 18:20:39 +0100 Subject: [PATCH 7/7] =?UTF-8?q?Added=20Ga=C3=ABl=20Varoquaux=20to=20the=20?= =?UTF-8?q?corner=20detector=20contributors?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CONTRIBUTORS.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 2e97d7e1..a5f83b5c 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -87,3 +87,6 @@ - W. Randolph Franklin Point in polygon test. + +- Gaƫl Varoquaux + Harris corner detector