diff --git a/doc/examples/plot_harris.py b/doc/examples/plot_harris.py index 9c499a3b..0e0222cb 100644 --- a/doc/examples/plot_harris.py +++ b/doc/examples/plot_harris.py @@ -3,14 +3,17 @@ 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 from skimage import data, img_as_float -from skimage.filter import harris +from skimage.feature import harris def plot_harris_points(image, filtered_coords): @@ -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/feature/__init__.py b/skimage/feature/__init__.py index 6b3b7014..49721c73 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,2 +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/feature/harris.py b/skimage/feature/harris.py new file mode 100644 index 00000000..1a7e3b6d --- /dev/null +++ b/skimage/feature/harris.py @@ -0,0 +1,85 @@ +""" +Harris corner detector + +Inspired from Solem's implementation +http://www.janeriksolem.net/2009/01/harris-corner-detector-in-python.html +""" +from scipy import ndimage + +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 + + 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. + """ + harrisim = _compute_harris_response(image, eps=eps, + gaussian_deviation=gaussian_deviation) + coordinates = peak.peak_local_max(harrisim, min_distance=min_distance, + threshold=threshold) + return coordinates + diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py new file mode 100644 index 00000000..225f42cd --- /dev/null +++ b/skimage/feature/peak.py @@ -0,0 +1,48 @@ +import numpy as np +from scipy import ndimage + + +def peak_local_max(image, min_distance=10, threshold=0.1): + """Return coordinates of peaks in an image. + + Peaks are the local maxima in a region of `2 * min_distance + 1` + (i.e. peaks are separated by at least `min_distance`). + + 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 + 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[: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 peaks + coordinates = np.transpose(image_t.nonzero()) + + return coordinates + diff --git a/skimage/feature/tests/test_harris.py b/skimage/feature/tests/test_harris.py new file mode 100644 index 00000000..4715ee25 --- /dev/null +++ b/skimage/feature/tests/test_harris.py @@ -0,0 +1,48 @@ +import numpy as np + +from skimage import data +from skimage import img_as_float + +from skimage.feature import harris + + +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(): + 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(): + 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 (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 + testing.run_module_suite() + diff --git a/skimage/feature/tests/test_peak.py b/skimage/feature/tests/test_peak.py new file mode 100644 index 00000000..19f9e95f --- /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_local_max(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/__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 diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py deleted file mode 100644 index 7eeecd30..00000000 --- a/skimage/filter/harris.py +++ /dev/null @@ -1,118 +0,0 @@ -# -# Harris 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 - - -def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1): - """Compute the Harris corner detector response function - for each pixel in the 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 - 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 - - -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: - -------- - array: coordinates of interest points - """ - 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) - - # 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) diff --git a/skimage/filter/tests/test_harris.py b/skimage/filter/tests/test_harris.py deleted file mode 100644 index 853a00a0..00000000 --- a/skimage/filter/tests/test_harris.py +++ /dev/null @@ -1,42 +0,0 @@ -import numpy as np - -from skimage import data -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_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_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_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()