diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 38e2ea32..fd0656db 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -132,6 +132,7 @@ - François Boulogne Andres Method for circle perimeter, ellipse perimeter drawing. + Circular Hough Transform - Thouis Jones Vectorized operators for arrays of 16-bit ints. diff --git a/doc/examples/plot_circular_hough_transform.py b/doc/examples/plot_circular_hough_transform.py new file mode 100755 index 00000000..2574a1ac --- /dev/null +++ b/doc/examples/plot_circular_hough_transform.py @@ -0,0 +1,71 @@ +""" +======================== +Circular Hough Transform +======================== + +The Hough transform in its simplest form is a `method to detect +straight lines `__ +but it can also be used to detect circles. + +In the following example, the Hough transform is used to detect +coin positions and match their edges. We provide a range of +plausible radii. For each radius, two circles are extracted and +we finally keep the five most prominent candidates. +The result shows that coin positions are well-detected. + + +Algorithm overview +------------------ + +Given a black circle on a white background, we first guess its +radius (or a range of radii) to construct a new circle. +This circle is applied on each black pixel of the original picture +and the coordinates of this circle are voting in an accumulator. +From this geometrical construction, the original circle center +position receives the highest score. + +Note that the accumulator size is built to be larger than the +original picture in order to detect centers outside the frame. +Its size is extended by two times the larger radius. + +""" + + +import numpy as np +import matplotlib.pyplot as plt + +from skimage import data, filter +from skimage.transform import hough_circle +from skimage.feature import peak_local_max +from skimage.draw import circle_perimeter + +# Load picture and detect edges +image = data.coins()[0:95, 70:370] +edges = filter.canny(image, sigma=3, low_threshold=10, high_threshold=50) + +fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6)) + +# Detect two radii +hough_radii = np.arange(15, 30, 2) +hough_res = hough_circle(edges, hough_radii) + +centers = [] +accums = [] +radii = [] + +for radius, h in zip(hough_radii, hough_res): + # For each radius, extract two circles + peaks = peak_local_max(h, num_peaks=2) + centers.extend(peaks - hough_radii.max()) + accums.extend(h[peaks[:, 0], peaks[:, 1]]) + radii.extend([radius, radius]) + +# Draw the most prominent 5 circles +for idx in np.argsort(accums)[::-1][:5]: + center_x, center_y = centers[idx] + radius = radii[idx] + cx, cy = circle_perimeter(center_y, center_x, radius) + image[cy, cx] = 0 + +ax.imshow(image, cmap=plt.cm.gray) +plt.show() diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx index 5d1f27af..29a38eb7 100644 --- a/skimage/transform/_hough_transform.pyx +++ b/skimage/transform/_hough_transform.pyx @@ -5,9 +5,12 @@ import numpy as np cimport numpy as cnp +cimport cython + from libc.math cimport abs, fabs, sqrt, ceil from libc.stdlib cimport rand +from skimage.draw import circle_perimeter cdef double PI_2 = 1.5707963267948966 cdef double NEG_PI_2 = -PI_2 @@ -17,6 +20,65 @@ cdef inline Py_ssize_t round(double r): return ((r + 0.5) if (r > 0.0) else (r - 0.5)) +@cython.boundscheck(False) +def _hough_circle(cnp.ndarray img, \ + cnp.ndarray[ndim=1, dtype=cnp.intp_t] radius, \ + normalize=True): + """Perform a circular Hough transform. + + Parameters + ---------- + img : (M, N) ndarray + Input image with nonzero values representing edges. + radius : ndarray + Radii at which to compute the Hough transform. + normalize : boolean, optional + Normalize the accumulator with the number + of pixels used to draw the radius + + Returns + ------- + H : 3D ndarray (radius index, (M, N) ndarray) + Hough transform accumulator for each radius + + """ + if img.ndim != 2: + raise ValueError('The input image must be 2D.') + + # compute the nonzero indexes + cdef cnp.ndarray[ndim=1, dtype=cnp.intp_t] x, y + x, y = np.nonzero(img) + + # Offset the image + cdef int max_radius = radius.max() + x = x + max_radius + y = y + max_radius + + cdef int px, py + cdef cnp.ndarray[ndim=1, dtype=cnp.intp_t] tx, ty, circle_x, circle_y + cdef cnp.ndarray acc = np.zeros((radius.size, + img.shape[0] + 2 * max_radius, + img.shape[1] + 2 * max_radius)) + + for i, rad in enumerate(radius): + # Store in memory the circle of given radius + # centered at (0,0) + circle_x, circle_y = circle_perimeter(0, 0, rad) + + # For each non zero pixel + for (px, py) in zip(x, y): + # Plug the circle at (px, py), + # its coordinates are (tx, ty) + tx = circle_x + px + ty = circle_y + py + acc[i, tx, ty] += 1 + + if normalize: + acc[i] = acc[i] / len(circle_x) + + return acc + + def _hough(cnp.ndarray img, cnp.ndarray[ndim=1, dtype=cnp.double_t] theta=None): if img.ndim != 2: diff --git a/skimage/transform/hough_transform.py b/skimage/transform/hough_transform.py index 05d5cfbe..17b76826 100644 --- a/skimage/transform/hough_transform.py +++ b/skimage/transform/hough_transform.py @@ -1,4 +1,4 @@ -__all__ = ['hough', 'hough_peaks', 'probabilistic_hough'] +__all__ = ['hough', 'hough_line', 'hough_circle', 'hough_peaks', 'probabilistic_hough'] from itertools import izip as zip @@ -96,8 +96,15 @@ def probabilistic_hough(img, threshold=10, line_length=50, line_gap=10, """ return _probabilistic_hough(img, threshold, line_length, line_gap, theta) +from skimage._shared.utils import deprecated +@deprecated('hough_line') def hough(img, theta=None): + return hough_line(img, theta) + +from ._hough_transform import _hough_circle + +def hough_line(img, theta=None): """Perform a straight line Hough transform. Parameters @@ -138,6 +145,26 @@ def hough(img, theta=None): """ return _hough(img, theta) +def hough_circle(img, radius, normalize=True): + """Perform a circular Hough transform. + + Parameters + ---------- + img : (M, N) ndarray + Input image with nonzero values representing edges. + radius : ndarray + Radii at which to compute the Hough transform. + normalize : boolean, optional + Normalize the accumulator with the number + of pixels used to draw the radius + + Returns + ------- + H : 3D ndarray (radius index, (M, N) ndarray) + Hough transform accumulator for each radius + + """ + return _hough_circle(img, radius, normalize) def hough_peaks(hspace, angles, dists, min_distance=10, min_angle=10, threshold=None, num_peaks=np.inf): diff --git a/skimage/transform/tests/test_hough_transform.py b/skimage/transform/tests/test_hough_transform.py index a38e6e76..00427332 100644 --- a/skimage/transform/tests/test_hough_transform.py +++ b/skimage/transform/tests/test_hough_transform.py @@ -4,6 +4,7 @@ from numpy.testing import * import skimage.transform as tf import skimage.transform.hough_transform as ht from skimage.transform import probabilistic_hough +from skimage.draw import circle_perimeter def append_desc(func, description): @@ -14,8 +15,6 @@ def append_desc(func, description): return func -from skimage.transform import * - def test_hough(): # Generate a test image @@ -23,7 +22,7 @@ def test_hough(): for i in range(25, 75): img[100 - i, i] = 1 - out, angles, d = tf.hough(img) + out, angles, d = tf.hough_line(img) y, x = np.where(out == out.max()) dist = d[y[0]] @@ -37,7 +36,7 @@ def test_hough_angles(): img = np.zeros((10, 10)) img[0, 0] = 1 - out, angles, d = tf.hough(img, np.linspace(0, 360, 10)) + out, angles, d = tf.hough_line(img, np.linspace(0, 360, 10)) assert_equal(len(angles), 10) @@ -76,7 +75,7 @@ def test_hough_peaks_dist(): img = np.zeros((100, 100), dtype=np.bool_) img[:, 30] = True img[:, 40] = True - hspace, angles, dists = tf.hough(img) + hspace, angles, dists = tf.hough_line(img) assert len(tf.hough_peaks(hspace, angles, dists, min_distance=5)[0]) == 2 assert len(tf.hough_peaks(hspace, angles, dists, min_distance=15)[0]) == 1 @@ -86,17 +85,17 @@ def test_hough_peaks_angle(): img[:, 0] = True img[0, :] = True - hspace, angles, dists = tf.hough(img) + hspace, angles, dists = tf.hough_line(img) assert len(tf.hough_peaks(hspace, angles, dists, min_angle=45)[0]) == 2 assert len(tf.hough_peaks(hspace, angles, dists, min_angle=90)[0]) == 1 theta = np.linspace(0, np.pi, 100) - hspace, angles, dists = tf.hough(img, theta) + hspace, angles, dists = tf.hough_line(img, theta) assert len(tf.hough_peaks(hspace, angles, dists, min_angle=45)[0]) == 2 assert len(tf.hough_peaks(hspace, angles, dists, min_angle=90)[0]) == 1 theta = np.linspace(np.pi / 3, 4. / 3 * np.pi, 100) - hspace, angles, dists = tf.hough(img, theta) + hspace, angles, dists = tf.hough_line(img, theta) assert len(tf.hough_peaks(hspace, angles, dists, min_angle=45)[0]) == 2 assert len(tf.hough_peaks(hspace, angles, dists, min_angle=90)[0]) == 1 @@ -105,10 +104,25 @@ def test_hough_peaks_num(): img = np.zeros((100, 100), dtype=np.bool_) img[:, 30] = True img[:, 40] = True - hspace, angles, dists = tf.hough(img) + hspace, angles, dists = tf.hough_line(img) assert len(tf.hough_peaks(hspace, angles, dists, min_distance=0, min_angle=0, num_peaks=1)[0]) == 1 +def test_houghcircle(): + # Prepare picture + img = np.zeros((120, 100), dtype=int) + radius = 20 + x_0, y_0 = (99, 50) + x, y = circle_perimeter(y_0, x_0, radius) + img[y, x] = 1 + + out = tf.hough_circle(img, np.array([radius])) + + x, y = np.where(out[0] == out[0].max()) + # Offset for x_0, y_0 + assert_equal(x[0], x_0 + radius) + assert_equal(y[0], y_0 + radius) + if __name__ == "__main__": run_module_suite()