From fd0e88b9867b1f398912c8a6cd309c26d0d566db Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 9 Jan 2012 20:52:05 -0500 Subject: [PATCH 1/7] Minor modifications to Harris corner detector * Add some references. * Make keyword argument explicit in example. * Remove test class in favor of functions since no setup was required. * Clean up docstrings. --- doc/examples/plot_harris.py | 10 +++-- skimage/filter/harris.py | 53 ++++++++++++----------- skimage/filter/tests/test_harris.py | 67 ++++++++++++++++------------- 3 files changed, 72 insertions(+), 58 deletions(-) diff --git a/doc/examples/plot_harris.py b/doc/examples/plot_harris.py index 9c499a3b..2d818e41 100644 --- a/doc/examples/plot_harris.py +++ b/doc/examples/plot_harris.py @@ -3,8 +3,11 @@ Harris Corner detector =============================================================================== -The Harris corner filter detects interest points using edge detection in -multiple direction. +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 """ from matplotlib import pyplot as plt @@ -26,5 +29,6 @@ def plot_harris_points(image, filtered_coords): im = img_as_float(data.lena()) -filtered_coords = harris(im, 6) +filtered_coords = harris(im, min_distance=6) plot_harris_points(im, filtered_coords) + diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py index 7eeecd30..b0add2cd 100644 --- a/skimage/filter/harris.py +++ b/skimage/filter/harris.py @@ -1,8 +1,9 @@ -# -# Harris detector -# -# Inspired from Solem's implementation -# http://www.janeriksolem.net/2009/01/harris-corner-detector-in-python.html +""" +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 @@ -14,18 +15,18 @@ def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1): Parameters ---------- - image: ndarray of floats - Input image + image : ndarray of floats + Input image. - eps: float, optional - Normalisation factor + eps : float, optional + Normalisation factor. - gaussian_deviation: integer, optional - Standard deviation used for the Gaussian kernel + gaussian_deviation : integer, optional + Standard deviation used for the Gaussian kernel. Returns -------- - image: (M, N) ndarray + image : (M, N) ndarray Harris image response """ if len(image.shape) == 3: @@ -43,6 +44,8 @@ def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1): # 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) # Non maximum filter of size 3 @@ -65,24 +68,25 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6, Parameters ---------- - image: ndarray of floats - Input image + image : ndarray of floats + Input image. - min_distance: int, optional - Minimum number of pixels separating interest points and image boundary + min_distance : int, optional + Minimum number of pixels separating interest points and image boundary. - threshold: float, optional + threshold : float, optional Relative threshold impacting the number of interest points. - eps: float, optional - Normalisation factor + eps : float, optional + Normalisation factor. - gaussian_deviation: integer, optional - Standard deviation used for the Gaussian kernel + gaussian_deviation : integer, optional + Standard deviation used for the Gaussian kernel. - returns: - -------- - array: coordinates of interest points + Returns + ------- + coordinates : (N, 2) array + (row, column) coordinates of interest points. """ harrisim = _compute_harris_response(image, eps=eps, gaussian_deviation=gaussian_deviation) @@ -116,3 +120,4 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6, (coords[i][1] - min_distance):(coords[i][1] + min_distance)] = 0 return np.array(filtered_coords) + diff --git a/skimage/filter/tests/test_harris.py b/skimage/filter/tests/test_harris.py index 853a00a0..295ab358 100644 --- a/skimage/filter/tests/test_harris.py +++ b/skimage/filter/tests/test_harris.py @@ -6,37 +6,42 @@ from skimage import img_as_float from skimage.filter import harris -class TestHarris(): - def test_square_image(self): - im = np.zeros((50, 50)).astype(float) - im[:25, :25] = 1. - results = harris(im) - assert results.any() - assert len(results) == 1 +def test_square_image(): + im = np.zeros((50, 50)).astype(float) + im[:25, :25] = 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(im) - assert results.any() - 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 = 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, min_distance=3) - print results - assert (results == np.array([[6, 6]])).all() +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) + print results + assert (results == np.array([[6, 6]])).all() + +def test_rotated_lena(): + """ + 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) + im_rotated = im.T + results_rotated = harris(im_rotated) + assert (results[:, 0] == results_rotated[:, 1]).all() + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() - 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) - im_rotated = im.T - results_rotated = harris(im_rotated) - assert (results[:, 0] == results_rotated[:, 1]).all() From 63f17344e59f0fa27920a9221ddb28fe64dec0b7 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 9 Jan 2012 22:24:56 -0500 Subject: [PATCH 2/7] Fix peak detection algorithm. Reverse sorted values so that max peak is checked first. --- skimage/filter/harris.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py index b0add2cd..2f5b2198 100644 --- a/skimage/filter/harris.py +++ b/skimage/filter/harris.py @@ -103,7 +103,7 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6, candidate_values = harrisim[candidates] # sort candidates - index = np.argsort(candidate_values) + index = np.argsort(candidate_values)[::-1] # store allowed point locations in array allowed_locations = np.zeros(harrisim.shape) From 83dc4a1d40e4145d74d3ba82c22f3ed6627cdc10 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 9 Jan 2012 22:28:56 -0500 Subject: [PATCH 3/7] Refactor peak detection algorithm from Harris detector. --- skimage/feature/__init__.py | 1 + skimage/feature/peak.py | 69 ++++++++++++++++++++++++++++++ skimage/feature/tests/test_peak.py | 24 +++++++++++ skimage/filter/harris.py | 48 +++------------------ 4 files changed, 99 insertions(+), 43 deletions(-) create mode 100644 skimage/feature/peak.py create mode 100644 skimage/feature/tests/test_peak.py diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 6b3b7014..b7046358 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,2 +1,3 @@ from hog import hog from greycomatrix import greycomatrix, greycoprops +from peak import peak_min_dist diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py new file mode 100644 index 00000000..317622d6 --- /dev/null +++ b/skimage/feature/peak.py @@ -0,0 +1,69 @@ +import numpy as np +from scipy import ndimage + + +def peak_min_dist(image, min_distance=10, threshold=0.1): + """Return coordinates of peaks in an image. + + Candidate peaks are determined by a relative `threshold`, and peaks that + are too close (as determined by `min_distance`) to larger peaks are + rejected. + + Parameters + ---------- + image: ndarray of floats + Input image. + + min_distance: int, optional + Minimum number of pixels separating peaks and image boundary. + + threshold: float, optional + Candidate peaks are calculated as `max(image) * threshold`. + + Returns + ------- + coordinates : (N, 2) array + (row, column) coordinates of peaks. + """ + image = image.copy() + # Non maximum filter of size 3 + image_max = ndimage.maximum_filter(image, 3, mode='constant') + mask = (image == image_max) + image *= mask + + # Remove the image borders + image[:3] = 0 + image[-3:] = 0 + image[:, :3] = 0 + image[:, -3:] = 0 + + # find top corner candidates above a threshold + corner_threshold = np.max(image.ravel()) * threshold + image_t = (image >= corner_threshold) * 1 + + # get coordinates of candidates + candidates = image_t.nonzero() + coords = np.transpose(candidates) + + # ...and their values + candidate_values = image[candidates] + + # sort candidates + index = np.argsort(candidate_values)[::-1] + + # store allowed point locations in array + allowed_locations = np.zeros(image.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[tuple(coords[i])] == 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) + diff --git a/skimage/feature/tests/test_peak.py b/skimage/feature/tests/test_peak.py new file mode 100644 index 00000000..99b9bc27 --- /dev/null +++ b/skimage/feature/tests/test_peak.py @@ -0,0 +1,24 @@ +import numpy as np + +from skimage import feature + + +def test_noisy_peaks(): + peak_locations = [(7, 7), (7, 13), (13, 7), (13, 13)] + + # image with noise of amplitude 0.8 and peaks of amplitude 1 + image = 0.8 * np.random.random((20, 20)) + for r, c in peak_locations: + image[r, c] = 1 + + peaks_detected = feature.peak_min_dist(image, min_distance=5) + + assert len(peaks_detected) == len(peak_locations) + for loc in peaks_detected: + assert tuple(loc) in peak_locations + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() + diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py index 2f5b2198..0298ddfe 100644 --- a/skimage/filter/harris.py +++ b/skimage/filter/harris.py @@ -4,10 +4,10 @@ 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 import feature + def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1): """Compute the Harris corner detector response function @@ -48,17 +48,6 @@ def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1): # Alison Noble, "Descriptions of Image Surfaces", PhD thesis (1989) harris = Wdet / (Wtr + eps) - # Non maximum filter of size 3 - harris_max = ndimage.maximum_filter(harris, 3, mode='constant') - mask = (harris == harris_max) - harris *= mask - - # Remove the image borders - harris[:3] = 0 - harris[-3:] = 0 - harris[:, :3] = 0 - harris[:, -3:] = 0 - return harris @@ -90,34 +79,7 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6, """ harrisim = _compute_harris_response(image, eps=eps, gaussian_deviation=gaussian_deviation) - - # find top corner candidates above a threshold - corner_threshold = np.max(harrisim.ravel()) * threshold - harrisim_t = (harrisim >= corner_threshold) * 1 - - # get coordinates of candidates - candidates = harrisim_t.nonzero() - coords = np.transpose(candidates) - - # ...and their values - candidate_values = harrisim[candidates] - - # sort candidates - index = np.argsort(candidate_values)[::-1] - - # 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[tuple(coords[i])] == 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) + coordinates = feature.peak_min_dist(harrisim, min_distance=min_distance, + threshold=threshold) + return coordinates From 4831185b4cde4c3450b4c543c08028c841f8e6f8 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 9 Jan 2012 22:43:52 -0500 Subject: [PATCH 4/7] Make test more robust. Test shouldn't depend on the order of the detected corners. --- skimage/filter/tests/test_harris.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/skimage/filter/tests/test_harris.py b/skimage/filter/tests/test_harris.py index 295ab358..16204e0f 100644 --- a/skimage/filter/tests/test_harris.py +++ b/skimage/filter/tests/test_harris.py @@ -38,7 +38,8 @@ def test_rotated_lena(): results = harris(im) im_rotated = im.T results_rotated = harris(im_rotated) - assert (results[:, 0] == results_rotated[:, 1]).all() + 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 8bb8d93552a5773383558ba033b2e7689043555e Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 9 Jan 2012 22:45:21 -0500 Subject: [PATCH 5/7] Use ndimage.maximum_filter to simplify peak detection. This change should give identitical results to the previous implementation. --- skimage/feature/peak.py | 40 ++++++++++------------------------------ 1 file changed, 10 insertions(+), 30 deletions(-) diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py index 317622d6..55f6010c 100644 --- a/skimage/feature/peak.py +++ b/skimage/feature/peak.py @@ -26,44 +26,24 @@ def peak_min_dist(image, min_distance=10, threshold=0.1): (row, column) coordinates of peaks. """ image = image.copy() - # Non maximum filter of size 3 - image_max = ndimage.maximum_filter(image, 3, mode='constant') + # Non maximum filter + size = 2 * min_distance + 1 + image_max = ndimage.maximum_filter(image, size=size, mode='constant') mask = (image == image_max) image *= mask # Remove the image borders - image[:3] = 0 - image[-3:] = 0 - image[:, :3] = 0 - image[:, -3:] = 0 + image[:min_distance] = 0 + image[-min_distance:] = 0 + image[:, :min_distance] = 0 + image[:, -min_distance:] = 0 # find top corner candidates above a threshold corner_threshold = np.max(image.ravel()) * threshold image_t = (image >= corner_threshold) * 1 - # get coordinates of candidates - candidates = image_t.nonzero() - coords = np.transpose(candidates) + # get coordinates of peaks + coordinates = np.transpose(image_t.nonzero()) - # ...and their values - candidate_values = image[candidates] - - # sort candidates - index = np.argsort(candidate_values)[::-1] - - # store allowed point locations in array - allowed_locations = np.zeros(image.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[tuple(coords[i])] == 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) + return coordinates From 1c9b340fe6de94ddb1ffa32afe729ef416921107 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 9 Jan 2012 22:54:12 -0500 Subject: [PATCH 6/7] Rename peak_min_dist to peak_local_max. --- skimage/feature/__init__.py | 2 +- skimage/feature/peak.py | 7 +++---- skimage/feature/tests/test_peak.py | 2 +- skimage/filter/harris.py | 2 +- 4 files changed, 6 insertions(+), 7 deletions(-) diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index b7046358..c63a773d 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,3 +1,3 @@ from hog import hog from greycomatrix import greycomatrix, greycoprops -from peak import peak_min_dist +from peak import peak_local_max diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py index 55f6010c..225f42cd 100644 --- a/skimage/feature/peak.py +++ b/skimage/feature/peak.py @@ -2,12 +2,11 @@ import numpy as np from scipy import ndimage -def peak_min_dist(image, min_distance=10, threshold=0.1): +def peak_local_max(image, min_distance=10, threshold=0.1): """Return coordinates of peaks in an image. - Candidate peaks are determined by a relative `threshold`, and peaks that - are too close (as determined by `min_distance`) to larger peaks are - rejected. + Peaks are the local maxima in a region of `2 * min_distance + 1` + (i.e. peaks are separated by at least `min_distance`). Parameters ---------- diff --git a/skimage/feature/tests/test_peak.py b/skimage/feature/tests/test_peak.py index 99b9bc27..19f9e95f 100644 --- a/skimage/feature/tests/test_peak.py +++ b/skimage/feature/tests/test_peak.py @@ -11,7 +11,7 @@ def test_noisy_peaks(): for r, c in peak_locations: image[r, c] = 1 - peaks_detected = feature.peak_min_dist(image, min_distance=5) + peaks_detected = feature.peak_local_max(image, min_distance=5) assert len(peaks_detected) == len(peak_locations) for loc in peaks_detected: diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py index 0298ddfe..feb848ac 100644 --- a/skimage/filter/harris.py +++ b/skimage/filter/harris.py @@ -79,7 +79,7 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6, """ harrisim = _compute_harris_response(image, eps=eps, gaussian_deviation=gaussian_deviation) - coordinates = feature.peak_min_dist(harrisim, min_distance=min_distance, + coordinates = feature.peak_local_max(harrisim, min_distance=min_distance, threshold=threshold) return coordinates From 96ba7eb4bddd5e96e0a85b48a899ca1e624bea96 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 2 Feb 2012 22:45:23 -0500 Subject: [PATCH 7/7] Move harris corner detection from filter to feature subpackage. --- doc/examples/plot_harris.py | 2 +- skimage/feature/__init__.py | 1 + skimage/{filter => feature}/harris.py | 4 ++-- skimage/{filter => feature}/tests/test_harris.py | 2 +- skimage/filter/__init__.py | 1 - 5 files changed, 5 insertions(+), 5 deletions(-) rename skimage/{filter => feature}/harris.py (95%) rename skimage/{filter => feature}/tests/test_harris.py (97%) diff --git a/doc/examples/plot_harris.py b/doc/examples/plot_harris.py index 2d818e41..0e0222cb 100644 --- a/doc/examples/plot_harris.py +++ b/doc/examples/plot_harris.py @@ -13,7 +13,7 @@ detection in multiple directions. from matplotlib import pyplot as plt from skimage import data, img_as_float -from skimage.filter import harris +from skimage.feature import harris def plot_harris_points(image, filtered_coords): diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index c63a773d..49721c73 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,3 +1,4 @@ from hog import hog from greycomatrix import greycomatrix, greycoprops from peak import peak_local_max +from harris import harris diff --git a/skimage/filter/harris.py b/skimage/feature/harris.py similarity index 95% rename from skimage/filter/harris.py rename to skimage/feature/harris.py index feb848ac..1a7e3b6d 100644 --- a/skimage/filter/harris.py +++ b/skimage/feature/harris.py @@ -6,7 +6,7 @@ http://www.janeriksolem.net/2009/01/harris-corner-detector-in-python.html """ from scipy import ndimage -from skimage import feature +from . import peak def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1): @@ -79,7 +79,7 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6, """ harrisim = _compute_harris_response(image, eps=eps, gaussian_deviation=gaussian_deviation) - coordinates = feature.peak_local_max(harrisim, min_distance=min_distance, + coordinates = peak.peak_local_max(harrisim, min_distance=min_distance, threshold=threshold) return coordinates diff --git a/skimage/filter/tests/test_harris.py b/skimage/feature/tests/test_harris.py similarity index 97% rename from skimage/filter/tests/test_harris.py rename to skimage/feature/tests/test_harris.py index 16204e0f..4715ee25 100644 --- a/skimage/filter/tests/test_harris.py +++ b/skimage/feature/tests/test_harris.py @@ -3,7 +3,7 @@ import numpy as np from skimage import data from skimage import img_as_float -from skimage.filter import harris +from skimage.feature import harris def test_square_image(): diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index 529331f0..1acc33f2 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -5,4 +5,3 @@ 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