mirror of
https://github.com/wassname/scikit-image.git
synced 2026-08-08 11:26:12 +08:00
Merge pull request #110 from tonysyu/peak-detection
ENH: Add peak detection.
This commit is contained in:
Stefan van der Walt
@@ -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)
|
||||
|
||||
|
||||
@@ -1,2 +1,4 @@
|
||||
from hog import hog
|
||||
from greycomatrix import greycomatrix, greycoprops
|
||||
from peak import peak_local_max
|
||||
from harris import harris
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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()
|
||||
|
||||
@@ -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()
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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)
|
||||
@@ -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()
|
||||
Reference in New Issue
Block a user