mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-07 07:43:56 +08:00
Merge pull request #321 from ahojnnes/interest-points
ENH: Improved corner detection.
This commit is contained in:
Stefan van der Walt
@@ -0,0 +1,38 @@
|
||||
"""
|
||||
================
|
||||
Corner detection
|
||||
================
|
||||
|
||||
Detect corner points using the Harris corner detector and determine subpixel
|
||||
position of corners.
|
||||
|
||||
.. [1] http://en.wikipedia.org/wiki/Corner_detection
|
||||
.. [2] http://en.wikipedia.org/wiki/Interest_point_detection
|
||||
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
from skimage import data
|
||||
from skimage.feature import corner_harris, corner_subpix, corner_peaks
|
||||
from skimage.transform import warp, AffineTransform
|
||||
from skimage.draw import ellipse
|
||||
|
||||
tform = AffineTransform(scale=(1.3, 1.1), rotation=1, shear=0.7,
|
||||
translation=(210, 50))
|
||||
image = warp(data.checkerboard(), tform.inverse, output_shape=(350, 350))
|
||||
rr, cc = ellipse(310, 175, 10, 100)
|
||||
image[rr, cc] = 1
|
||||
image[180:230, 10:60] = 1
|
||||
image[230:280, 60:110] = 1
|
||||
|
||||
coords = corner_peaks(corner_harris(image), min_distance=5)
|
||||
coords_subpix = corner_subpix(image, coords, window_size=13)
|
||||
|
||||
plt.gray()
|
||||
plt.imshow(image, interpolation='nearest')
|
||||
plt.plot(coords[:, 1], coords[:, 0], '.b', markersize=3)
|
||||
plt.plot(coords_subpix[:, 1], coords_subpix[:, 0], '+r', markersize=15)
|
||||
plt.axis((0, 350, 350, 0))
|
||||
plt.show()
|
||||
@@ -1,42 +0,0 @@
|
||||
"""
|
||||
===============================================================================
|
||||
Harris Corner detector
|
||||
===============================================================================
|
||||
|
||||
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
|
||||
"""
|
||||
import numpy as np
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
from skimage import data, img_as_float
|
||||
from skimage.feature import harris
|
||||
|
||||
|
||||
def plot_harris_points(image, filtered_coords):
|
||||
""" plots corners found in image"""
|
||||
|
||||
plt.imshow(image)
|
||||
y, x = np.transpose(filtered_coords)
|
||||
plt.plot(x, y, 'b.')
|
||||
plt.axis('off')
|
||||
|
||||
# display results
|
||||
plt.figure(figsize=(8, 3))
|
||||
im_lena = img_as_float(data.lena())
|
||||
im_text = img_as_float(data.text())
|
||||
|
||||
filtered_coords = harris(im_lena, min_distance=4)
|
||||
|
||||
plt.axes([0, 0, 0.3, 0.95])
|
||||
plot_harris_points(im_lena, filtered_coords)
|
||||
|
||||
filtered_coords = harris(im_text, min_distance=4)
|
||||
|
||||
plt.axes([0.2, 0, 0.77, 1])
|
||||
plot_harris_points(im_text, filtered_coords)
|
||||
|
||||
plt.show()
|
||||
@@ -1,5 +1,7 @@
|
||||
from ._hog import hog
|
||||
from .texture import greycomatrix, greycoprops, local_binary_pattern
|
||||
from .peak import peak_local_max
|
||||
from ._harris import harris
|
||||
from .corner import (corner_kitchen_rosenfeld, corner_harris, corner_shi_tomasi,
|
||||
corner_foerstner, corner_subpix, corner_peaks)
|
||||
from .corner_cy import corner_moravec
|
||||
from .template import match_template
|
||||
|
||||
@@ -1,109 +0,0 @@
|
||||
"""
|
||||
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.
|
||||
|
||||
Examples
|
||||
-------
|
||||
>>> square = np.zeros([10,10])
|
||||
>>> square[2:8,2:8] = 1
|
||||
>>> square
|
||||
array([[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.],
|
||||
[ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.],
|
||||
[ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.],
|
||||
[ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.],
|
||||
[ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.],
|
||||
[ 0., 0., 1., 1., 1., 1., 1., 1., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
|
||||
>>> harris(square, min_distance=1)
|
||||
|
||||
Corners of the square
|
||||
|
||||
array([[3, 3],
|
||||
[3, 6],
|
||||
[6, 3],
|
||||
[6, 6]])
|
||||
"""
|
||||
|
||||
harrisim = _compute_harris_response(image, eps=eps,
|
||||
gaussian_deviation=gaussian_deviation)
|
||||
coordinates = peak.peak_local_max(harrisim, min_distance=min_distance,
|
||||
threshold_rel=threshold)
|
||||
return coordinates
|
||||
@@ -0,0 +1,505 @@
|
||||
import numpy as np
|
||||
from scipy import ndimage
|
||||
from scipy import stats
|
||||
from skimage.color import rgb2grey
|
||||
from skimage.util import img_as_float
|
||||
from skimage.feature import peak_local_max
|
||||
|
||||
|
||||
def _compute_derivatives(image):
|
||||
"""Compute derivatives in x and y direction using the Sobel operator.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input image.
|
||||
|
||||
Returns
|
||||
-------
|
||||
imx : ndarray
|
||||
Derivative in x-direction.
|
||||
imy : ndarray
|
||||
Derivative in y-direction.
|
||||
|
||||
"""
|
||||
|
||||
imy = ndimage.sobel(image, axis=0, mode='constant', cval=0)
|
||||
imx = ndimage.sobel(image, axis=1, mode='constant', cval=0)
|
||||
|
||||
return imx, imy
|
||||
|
||||
|
||||
def _compute_auto_correlation(image, sigma):
|
||||
"""Compute auto-correlation matrix using sum of squared differences.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input image.
|
||||
sigma : float
|
||||
Standard deviation used for the Gaussian kernel, which is used as
|
||||
weighting function for the auto-correlation matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
Axx : ndarray
|
||||
Element of the auto-correlation matrix for each pixel in input image.
|
||||
Axy : ndarray
|
||||
Element of the auto-correlation matrix for each pixel in input image.
|
||||
Ayy : ndarray
|
||||
Element of the auto-correlation matrix for each pixel in input image.
|
||||
|
||||
"""
|
||||
|
||||
if image.ndim == 3:
|
||||
image = img_as_float(rgb2grey(image))
|
||||
|
||||
imx, imy = _compute_derivatives(image)
|
||||
|
||||
# structure tensore
|
||||
Axx = ndimage.gaussian_filter(imx * imx, sigma, mode='constant', cval=0)
|
||||
Axy = ndimage.gaussian_filter(imx * imy, sigma, mode='constant', cval=0)
|
||||
Ayy = ndimage.gaussian_filter(imy * imy, sigma, mode='constant', cval=0)
|
||||
|
||||
return Axx, Axy, Ayy
|
||||
|
||||
|
||||
def corner_kitchen_rosenfeld(image):
|
||||
"""Compute Kitchen and Rosenfeld corner measure response image.
|
||||
|
||||
The corner measure is calculated as follows::
|
||||
|
||||
(imxx * imy**2 + imyy * imx**2 - 2 * imxy * imx * imy)
|
||||
------------------------------------------------------
|
||||
(imx**2 + imy**2)
|
||||
|
||||
Where imx and imy are the first and imxx, imxy, imyy the second derivatives.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input image.
|
||||
|
||||
Returns
|
||||
-------
|
||||
response : ndarray
|
||||
Kitchen and Rosenfeld response image.
|
||||
|
||||
"""
|
||||
|
||||
imx, imy = _compute_derivatives(image)
|
||||
imxx, imxy = _compute_derivatives(imx)
|
||||
imyx, imyy = _compute_derivatives(imy)
|
||||
|
||||
response = (imxx * imy**2 + imyy * imx**2 - 2 * imxy * imx * imy) \
|
||||
/ (imx**2 + imy**2)
|
||||
|
||||
return response
|
||||
|
||||
|
||||
def corner_harris(image, method='k', k=0.05, eps=1e-6, sigma=1):
|
||||
"""Compute Harris corner measure response image.
|
||||
|
||||
This corner detector uses information from the auto-correlation matrix A::
|
||||
|
||||
A = [(imx**2) (imx*imy)] = [Axx Axy]
|
||||
[(imx*imy) (imy**2)] [Axy Ayy]
|
||||
|
||||
Where imx and imy are the first derivatives averaged with a gaussian filter.
|
||||
The corner measure is then defined as::
|
||||
|
||||
det(A) - k * trace(A)**2
|
||||
|
||||
or::
|
||||
|
||||
2 * det(A) / (trace(A) + eps)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input image.
|
||||
method : {'k', 'eps'}, optional
|
||||
Method to compute the response image from the auto-correlation matrix.
|
||||
k : float, optional
|
||||
Sensitivity factor to separate corners from edges, typically in range
|
||||
`[0, 0.2]`. Small values of k result in detection of sharp corners.
|
||||
eps : float, optional
|
||||
Normalisation factor (Noble's corner measure).
|
||||
sigma : float, optional
|
||||
Standard deviation used for the Gaussian kernel, which is used as
|
||||
weighting function for the auto-correlation matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
response : ndarray
|
||||
Harris response image.
|
||||
|
||||
References
|
||||
----------
|
||||
..[1] http://kiwi.cs.dal.ca/~dparks/CornerDetection/harris.htm
|
||||
..[2] http://en.wikipedia.org/wiki/Corner_detection
|
||||
|
||||
Examples
|
||||
-------
|
||||
>>> from skimage.feature import corner_harris, corner_peaks
|
||||
>>> square = np.zeros([10, 10])
|
||||
>>> square[2:8, 2:8] = 1
|
||||
>>> square
|
||||
array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
|
||||
>>> corner_peaks(corner_harris(square), min_distance=1)
|
||||
array([[2, 2],
|
||||
[2, 7],
|
||||
[7, 2],
|
||||
[7, 7]])
|
||||
|
||||
"""
|
||||
|
||||
Axx, Axy, Ayy = _compute_auto_correlation(image, sigma)
|
||||
|
||||
# determinant
|
||||
detA = Axx * Ayy - Axy**2
|
||||
# trace
|
||||
traceA = Axx + Ayy
|
||||
|
||||
if method == 'k':
|
||||
response = detA - k * traceA**2
|
||||
else:
|
||||
response = 2 * detA / (traceA + eps)
|
||||
|
||||
return response
|
||||
|
||||
|
||||
def corner_shi_tomasi(image, sigma=1):
|
||||
"""Compute Shi-Tomasi (Kanade-Tomasi) corner measure response image.
|
||||
|
||||
This corner detector uses information from the auto-correlation matrix A::
|
||||
|
||||
A = [(imx**2) (imx*imy)] = [Axx Axy]
|
||||
[(imx*imy) (imy**2)] [Axy Ayy]
|
||||
|
||||
Where imx and imy are the first derivatives averaged with a gaussian filter.
|
||||
The corner measure is then defined as the smaller eigenvalue of A::
|
||||
|
||||
((Axx + Ayy) - sqrt((Axx - Ayy)**2 + 4 * Axy**2)) / 2
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input image.
|
||||
sigma : float, optional
|
||||
Standard deviation used for the Gaussian kernel, which is used as
|
||||
weighting function for the auto-correlation matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
response : ndarray
|
||||
Shi-Tomasi response image.
|
||||
|
||||
References
|
||||
----------
|
||||
..[1] http://kiwi.cs.dal.ca/~dparks/CornerDetection/harris.htm
|
||||
..[2] http://en.wikipedia.org/wiki/Corner_detection
|
||||
|
||||
Examples
|
||||
-------
|
||||
>>> from skimage.feature import corner_shi_tomasi, corner_peaks
|
||||
>>> square = np.zeros([10, 10])
|
||||
>>> square[2:8, 2:8] = 1
|
||||
>>> square
|
||||
array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
|
||||
>>> corner_peaks(corner_shi_tomasi(square), min_distance=1)
|
||||
array([[2, 2],
|
||||
[2, 7],
|
||||
[7, 2],
|
||||
[7, 7]])
|
||||
|
||||
"""
|
||||
|
||||
Axx, Axy, Ayy = _compute_auto_correlation(image, sigma)
|
||||
|
||||
# minimum eigenvalue of A
|
||||
response = ((Axx + Ayy) - np.sqrt((Axx - Ayy)**2 + 4 * Axy**2)) / 2
|
||||
|
||||
return response
|
||||
|
||||
|
||||
def corner_foerstner(image, sigma=1):
|
||||
"""Compute Foerstner corner measure response image.
|
||||
|
||||
This corner detector uses information from the auto-correlation matrix A::
|
||||
|
||||
A = [(imx**2) (imx*imy)] = [Axx Axy]
|
||||
[(imx*imy) (imy**2)] [Axy Ayy]
|
||||
|
||||
Where imx and imy are the first derivatives averaged with a gaussian filter.
|
||||
The corner measure is then defined as::
|
||||
|
||||
w = det(A) / trace(A) (size of error ellipse)
|
||||
q = 4 * det(A) / trace(A)**2 (roundness of error ellipse)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input image.
|
||||
sigma : float, optional
|
||||
Standard deviation used for the Gaussian kernel, which is used as
|
||||
weighting function for the auto-correlation matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
w : ndarray
|
||||
Error ellipse sizes.
|
||||
q : ndarray
|
||||
Roundness of error ellipse.
|
||||
|
||||
References
|
||||
----------
|
||||
..[1] http://www.ipb.uni-bonn.de/uploads/tx_ikgpublication/\
|
||||
foerstner87.fast.pdf
|
||||
..[2] http://en.wikipedia.org/wiki/Corner_detection
|
||||
|
||||
Examples
|
||||
-------
|
||||
>>> from skimage.feature import corner_foerstner, corner_peaks
|
||||
>>> square = np.zeros([10, 10])
|
||||
>>> square[2:8, 2:8] = 1
|
||||
>>> square
|
||||
array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
|
||||
>>> w, q = corner_foerstner(square)
|
||||
>>> accuracy_thresh = 0.5
|
||||
>>> roundness_thresh = 0.3
|
||||
>>> foerstner = (q > roundness_thresh) * (w > accuracy_thresh) * w
|
||||
>>> corner_peaks(foerstner, min_distance=1)
|
||||
array([[2, 2],
|
||||
[2, 7],
|
||||
[7, 2],
|
||||
[7, 7]])
|
||||
|
||||
"""
|
||||
|
||||
Axx, Axy, Ayy = _compute_auto_correlation(image, sigma)
|
||||
|
||||
# determinant
|
||||
detA = Axx * Ayy - Axy**2
|
||||
# trace
|
||||
traceA = Axx + Ayy
|
||||
|
||||
w = detA / traceA
|
||||
q = 4 * detA / traceA**2
|
||||
|
||||
return w, q
|
||||
|
||||
|
||||
def corner_subpix(image, corners, window_size=11, alpha=0.99):
|
||||
"""Determine subpixel position of corners.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input image.
|
||||
corners : (N, 2) ndarray
|
||||
Corner coordinates `(row, col)`.
|
||||
window_size : int, optional
|
||||
Search window size for subpixel estimation.
|
||||
alpha : float, optional
|
||||
Significance level for point classification.
|
||||
|
||||
Returns
|
||||
-------
|
||||
positions : (N, 2) ndarray
|
||||
Subpixel corner positions. NaN for "not classified" corners.
|
||||
|
||||
References
|
||||
----------
|
||||
..[1] http://www.ipb.uni-bonn.de/uploads/tx_ikgpublication/\
|
||||
foerstner87.fast.pdf
|
||||
..[2] http://en.wikipedia.org/wiki/Corner_detection
|
||||
|
||||
"""
|
||||
|
||||
# window extent in one direction
|
||||
wext = (window_size - 1) / 2
|
||||
|
||||
# normal equation arrays
|
||||
N_dot = np.zeros((2, 2), dtype=np.double)
|
||||
N_edge = np.zeros((2, 2), dtype=np.double)
|
||||
b_dot = np.zeros((2, ), dtype=np.double)
|
||||
b_edge = np.zeros((2, ), dtype=np.double)
|
||||
|
||||
# critical statistical test values
|
||||
redundancy = window_size**2 - 2
|
||||
t_crit_dot = stats.f.isf(1 - alpha, redundancy, redundancy)
|
||||
t_crit_edge = stats.f.isf(alpha, redundancy, redundancy)
|
||||
|
||||
# coordinates of pixels within window
|
||||
y, x = np.mgrid[- wext:wext + 1, - wext:wext + 1]
|
||||
|
||||
corners_subpix = np.zeros_like(corners, dtype=np.double)
|
||||
|
||||
for i, (y0, x0) in enumerate(corners):
|
||||
|
||||
# crop window around corner + border for sobel operator
|
||||
miny = y0 - wext - 1
|
||||
maxy = y0 + wext + 2
|
||||
minx = x0 - wext - 1
|
||||
maxx = x0 + wext + 2
|
||||
window = image[miny:maxy, minx:maxx]
|
||||
|
||||
winx, winy = _compute_derivatives(window)
|
||||
|
||||
# compute gradient suares and remove border
|
||||
winx_winx = (winx * winx)[1:-1, 1:-1]
|
||||
winx_winy = (winx * winy)[1:-1, 1:-1]
|
||||
winy_winy = (winy * winy)[1:-1, 1:-1]
|
||||
|
||||
# sum of squared differences (mean instead of gaussian filter)
|
||||
Axx = np.sum(winx_winx)
|
||||
Axy = np.sum(winx_winy)
|
||||
Ayy = np.sum(winy_winy)
|
||||
|
||||
# sum of squared differences weighted with coordinates
|
||||
# (mean instead of gaussian filter)
|
||||
bxx_x = np.sum(winx_winx * x)
|
||||
bxx_y = np.sum(winx_winx * y)
|
||||
bxy_x = np.sum(winx_winy * x)
|
||||
bxy_y = np.sum(winx_winy * y)
|
||||
byy_x = np.sum(winy_winy * x)
|
||||
byy_y = np.sum(winy_winy * y)
|
||||
|
||||
# normal equations for subpixel position
|
||||
N_dot[0, 0] = Axx
|
||||
N_dot[0, 1] = N_dot[1, 0] = - Axy
|
||||
N_dot[1, 1] = Ayy
|
||||
|
||||
N_edge[0, 0] = Ayy
|
||||
N_edge[0, 1] = N_edge[1, 0] = Axy
|
||||
N_edge[1, 1] = Axx
|
||||
|
||||
b_dot[:] = bxx_y - bxy_x, byy_x - bxy_y
|
||||
b_edge[:] = byy_y + bxy_x, bxx_x + bxy_y
|
||||
|
||||
# estimated positions
|
||||
est_dot = np.linalg.solve(N_dot, b_dot)
|
||||
est_edge = np.linalg.solve(N_edge, b_edge)
|
||||
|
||||
# residuals
|
||||
ry_dot = y - est_dot[0]
|
||||
rx_dot = x - est_dot[1]
|
||||
ry_edge = y - est_edge[0]
|
||||
rx_edge = x - est_edge[1]
|
||||
# squared residuals
|
||||
rxx_dot = rx_dot * rx_dot
|
||||
rxy_dot = rx_dot * ry_dot
|
||||
ryy_dot = ry_dot * ry_dot
|
||||
rxx_edge = rx_edge * rx_edge
|
||||
rxy_edge = rx_edge * ry_edge
|
||||
ryy_edge = ry_edge * ry_edge
|
||||
|
||||
# determine corner class (dot or edge)
|
||||
# variance for different models
|
||||
var_dot = np.sum(winx_winx * ryy_dot - 2 * winx_winy * rxy_dot \
|
||||
+ winy_winy * rxx_dot)
|
||||
var_edge = np.sum(winy_winy * ryy_edge + 2 * winx_winy * rxy_edge \
|
||||
+ winx_winx * rxx_edge)
|
||||
# test value (F-distributed)
|
||||
t = var_edge / var_dot
|
||||
# 1 for edge, -1 for dot, 0 for "not classified"
|
||||
corner_class = (t < t_crit_edge) - (t > t_crit_dot)
|
||||
|
||||
if corner_class == - 1:
|
||||
corners_subpix[i, :] = y0 + est_dot[0], x0 + est_dot[1]
|
||||
elif corner_class == 0:
|
||||
corners_subpix[i, :] = np.nan, np.nan
|
||||
elif corner_class == 1:
|
||||
corners_subpix[i, :] = y0 + est_edge[0], x0 + est_edge[1]
|
||||
|
||||
return corners_subpix
|
||||
|
||||
|
||||
def corner_peaks(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
|
||||
exclude_border=True, indices=True, num_peaks=np.inf,
|
||||
footprint=None, labels=None):
|
||||
"""Find corners in corner measure response image.
|
||||
|
||||
This differs from `skimage.feature.peak_local_max` in that it suppresses
|
||||
multiple connected peaks with the same accumulator value.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
See `skimage.feature.peak_local_max`.
|
||||
|
||||
Returns
|
||||
-------
|
||||
See `skimage.feature.peak_local_max`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from skimage.feature import peak_local_max, corner_peaks
|
||||
>>> response = np.zeros((5, 5))
|
||||
>>> response[2:4, 2:4] = 1
|
||||
>>> response
|
||||
array([[ 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 1., 1., 0.],
|
||||
[ 0., 0., 1., 1., 0.],
|
||||
[ 0., 0., 0., 0., 0.]])
|
||||
>>> peak_local_max(response, exclude_border=False)
|
||||
array([[2, 2],
|
||||
[2, 3],
|
||||
[3, 2],
|
||||
[3, 3]])
|
||||
>>> corner_peaks(response, exclude_border=False)
|
||||
array([[2, 2]])
|
||||
>>> corner_peaks(response, exclude_border=False, min_distance=0)
|
||||
array([[2, 2],
|
||||
[2, 3],
|
||||
[3, 2],
|
||||
[3, 3]])
|
||||
|
||||
"""
|
||||
|
||||
peaks = peak_local_max(image, min_distance=min_distance,
|
||||
threshold_abs=threshold_abs,
|
||||
threshold_rel=threshold_rel,
|
||||
exclude_border=exclude_border,
|
||||
indices=False, num_peaks=np.inf,
|
||||
footprint=footprint, labels=labels)
|
||||
if min_distance > 0:
|
||||
coords = np.transpose(peaks.nonzero())
|
||||
for r, c in coords:
|
||||
if peaks[r, c]:
|
||||
peaks[r - min_distance:r + min_distance + 1,
|
||||
c - min_distance:c + min_distance + 1] = False
|
||||
peaks[r, c] = True
|
||||
|
||||
if indices is True:
|
||||
return np.transpose(peaks.nonzero())
|
||||
else:
|
||||
return peaks
|
||||
@@ -0,0 +1,91 @@
|
||||
#cython: cdivision=True
|
||||
#cython: boundscheck=False
|
||||
#cython: nonecheck=False
|
||||
#cython: wraparound=False
|
||||
import numpy as np
|
||||
cimport numpy as cnp
|
||||
from libc.float cimport DBL_MAX
|
||||
|
||||
from skimage.color import rgb2grey
|
||||
from skimage.util import img_as_float
|
||||
|
||||
|
||||
def corner_moravec(image, int window_size=1):
|
||||
"""Compute Moravec corner measure response image.
|
||||
|
||||
This is one of the simplest corner detectors and is comparatively fast but
|
||||
has several limitations (e.g. not rotation invariant).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input image.
|
||||
window_size : int, optional
|
||||
Window size.
|
||||
|
||||
Returns
|
||||
-------
|
||||
response : ndarray
|
||||
Moravec response image.
|
||||
|
||||
References
|
||||
----------
|
||||
..[1] http://kiwi.cs.dal.ca/~dparks/CornerDetection/moravec.htm
|
||||
..[2] http://en.wikipedia.org/wiki/Corner_detection
|
||||
|
||||
Examples
|
||||
-------
|
||||
>>> from skimage.feature import moravec, peak_local_max
|
||||
>>> square = np.zeros([7, 7])
|
||||
>>> square[3, 3] = 1
|
||||
>>> square
|
||||
array([[ 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 0., 1., 0., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0.]])
|
||||
>>> moravec(square)
|
||||
array([[ 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 1., 1., 1., 0., 0.],
|
||||
[ 0., 0., 1., 2., 1., 0., 0.],
|
||||
[ 0., 0., 1., 1., 1., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0.],
|
||||
[ 0., 0., 0., 0., 0., 0., 0.]])
|
||||
"""
|
||||
|
||||
cdef int rows = image.shape[0]
|
||||
cdef int cols = image.shape[1]
|
||||
|
||||
cdef cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] cimage, out
|
||||
|
||||
if image.ndim == 3:
|
||||
cimage = rgb2grey(image)
|
||||
cimage = np.ascontiguousarray(img_as_float(image))
|
||||
|
||||
out = np.zeros(image.shape, dtype=np.double)
|
||||
|
||||
cdef double* image_data = <double*>cimage.data
|
||||
cdef double* out_data = <double*>out.data
|
||||
|
||||
cdef double msum, min_msum
|
||||
cdef int r, c, br, bc, mr, mc, a, b
|
||||
for r in range(2 * window_size, rows - 2 * window_size):
|
||||
for c in range(2 * window_size, cols - 2 * window_size):
|
||||
min_msum = DBL_MAX
|
||||
for br in range(r - window_size, r + window_size + 1):
|
||||
for bc in range(c - window_size, c + window_size + 1):
|
||||
if br != r and bc != c:
|
||||
msum = 0
|
||||
for mr in range(- window_size, window_size + 1):
|
||||
for mc in range(- window_size, window_size + 1):
|
||||
a = (r + mr) * cols + c + mc
|
||||
b = (br + mr) * cols + bc + mc
|
||||
msum += (image_data[a] - image_data[b]) ** 2
|
||||
min_msum = min(msum, min_msum)
|
||||
|
||||
out_data[r * cols + c] = min_msum
|
||||
|
||||
return out
|
||||
@@ -12,9 +12,12 @@ def configuration(parent_package='', top_path=None):
|
||||
config = Configuration('feature', parent_package, top_path)
|
||||
config.add_data_dir('tests')
|
||||
|
||||
cython(['corner_cy.pyx'], working_path=base_path)
|
||||
cython(['_texture.pyx'], working_path=base_path)
|
||||
cython(['_template.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('corner_cy', sources=['corner_cy.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_texture', sources=['_texture.c'],
|
||||
include_dirs=[get_numpy_include_dirs(), '../_shared'])
|
||||
config.add_extension('_template', sources=['_template.c'],
|
||||
|
||||
@@ -0,0 +1,115 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_equal
|
||||
|
||||
from skimage import data
|
||||
from skimage import img_as_float
|
||||
|
||||
from skimage.feature import (corner_moravec, corner_harris, corner_shi_tomasi,
|
||||
corner_subpix, peak_local_max, corner_peaks)
|
||||
|
||||
|
||||
def test_square_image():
|
||||
im = np.zeros((50, 50)).astype(float)
|
||||
im[:25, :25] = 1.
|
||||
|
||||
# Moravec
|
||||
results = peak_local_max(corner_moravec(im))
|
||||
# interest points along edge
|
||||
assert len(results) == 57
|
||||
|
||||
# Harris
|
||||
results = peak_local_max(corner_harris(im))
|
||||
# interest at corner
|
||||
assert len(results) == 1
|
||||
|
||||
# Shi-Tomasi
|
||||
results = peak_local_max(corner_shi_tomasi(im))
|
||||
# interest at corner
|
||||
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) * .2
|
||||
|
||||
# Moravec
|
||||
results = peak_local_max(corner_moravec(im))
|
||||
# undefined number of interest points
|
||||
assert results.any()
|
||||
|
||||
# Harris
|
||||
results = peak_local_max(corner_harris(im, sigma=1.5))
|
||||
assert len(results) == 1
|
||||
|
||||
# Shi-Tomasi
|
||||
results = peak_local_max(corner_shi_tomasi(im, sigma=1.5))
|
||||
assert len(results) == 1
|
||||
|
||||
|
||||
def test_squared_dot():
|
||||
im = np.zeros((50, 50))
|
||||
im[4:8, 4:8] = 1
|
||||
im = img_as_float(im)
|
||||
|
||||
# Moravec fails
|
||||
|
||||
# Harris
|
||||
results = peak_local_max(corner_harris(im))
|
||||
assert (results == np.array([[6, 6]])).all()
|
||||
|
||||
# Shi-Tomasi
|
||||
results = peak_local_max(corner_shi_tomasi(im))
|
||||
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))
|
||||
im_rotated = im.T
|
||||
|
||||
# Moravec
|
||||
results = peak_local_max(corner_moravec(im))
|
||||
results_rotated = peak_local_max(corner_moravec(im_rotated))
|
||||
assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all()
|
||||
assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all()
|
||||
|
||||
# Harris
|
||||
results = peak_local_max(corner_harris(im))
|
||||
results_rotated = peak_local_max(corner_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()
|
||||
|
||||
# Shi-Tomasi
|
||||
results = peak_local_max(corner_shi_tomasi(im))
|
||||
results_rotated = peak_local_max(corner_shi_tomasi(im_rotated))
|
||||
assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all()
|
||||
assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all()
|
||||
|
||||
|
||||
def test_subpix():
|
||||
img = np.zeros((50, 50))
|
||||
img[:25,:25] = 255
|
||||
img[25:,25:] = 255
|
||||
corner = peak_local_max(corner_harris(img), num_peaks=1)
|
||||
subpix = corner_subpix(img, corner)
|
||||
assert_array_equal(subpix[0], (24.5, 24.5))
|
||||
|
||||
|
||||
def test_corner_peaks():
|
||||
response = np.zeros((5, 5))
|
||||
response[2:4, 2:4] = 1
|
||||
|
||||
corners = corner_peaks(response, exclude_border=False)
|
||||
assert len(corners) == 1
|
||||
|
||||
corners = corner_peaks(response, exclude_border=False, min_distance=0)
|
||||
assert len(corners) == 4
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy import testing
|
||||
testing.run_module_suite()
|
||||
@@ -1,49 +0,0 @@
|
||||
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)
|
||||
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()
|
||||
Reference in New Issue
Block a user