diff --git a/doc/examples/plot_corner.py b/doc/examples/plot_corner.py new file mode 100644 index 00000000..4d01d177 --- /dev/null +++ b/doc/examples/plot_corner.py @@ -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() diff --git a/doc/examples/plot_harris.py b/doc/examples/plot_harris.py deleted file mode 100644 index 6212ac4a..00000000 --- a/doc/examples/plot_harris.py +++ /dev/null @@ -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() diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 1597e9a8..e0756bfe 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -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 diff --git a/skimage/feature/_harris.py b/skimage/feature/_harris.py deleted file mode 100644 index ae30a29e..00000000 --- a/skimage/feature/_harris.py +++ /dev/null @@ -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 diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py new file mode 100644 index 00000000..867acb5d --- /dev/null +++ b/skimage/feature/corner.py @@ -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 diff --git a/skimage/feature/corner_cy.pyx b/skimage/feature/corner_cy.pyx new file mode 100644 index 00000000..2f001ea4 --- /dev/null +++ b/skimage/feature/corner_cy.pyx @@ -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 = cimage.data + cdef double* out_data = 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 diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 6f820163..e769621d 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -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'], diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py new file mode 100644 index 00000000..46e5466f --- /dev/null +++ b/skimage/feature/tests/test_corner.py @@ -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() diff --git a/skimage/feature/tests/test_harris.py b/skimage/feature/tests/test_harris.py deleted file mode 100644 index 43bf28a3..00000000 --- a/skimage/feature/tests/test_harris.py +++ /dev/null @@ -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()