Merge pull request #431 from sciunto/houghcircle

ENH: Add circular Hough transform.

Conflicts:
	skimage/transform/_hough_transform.pyx

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.

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 <http://en.wikipedia.org/wiki/Hough_transform>`__
+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()

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 <Py_ssize_t>((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:

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

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