Refactor peak detection algorithm from Harris detector.

skimage/feature/__init__.py

@@ -1,2 +1,3 @@
 from hog import hog
 from greycomatrix import greycomatrix, greycoprops
+from peak import peak_min_dist

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)
+

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()
+

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