From 83dc4a1d40e4145d74d3ba82c22f3ed6627cdc10 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 9 Jan 2012 22:28:56 -0500 Subject: [PATCH] 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