hough transform for circles

skimage/transform/_hough_transform.pyx

@@ -7,7 +7,7 @@ import numpy as np
 cimport numpy as np
 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 +17,90 @@ 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(np.ndarray img, \
+                  np.ndarray[ndim=1, dtype=np.npy_intp] radius, \
+                  normalize=True):
+
+    if img.ndim != 2:
+        raise ValueError('The input image must be 2D.')
+
+    # compute the nonzero indexes
+    cdef np.ndarray[ndim=1, dtype=np.npy_intp] x, y
+    y, x = np.PyArray_Nonzero(img)
+
+    # Offset the image
+    max_radius = radius.max()
+    x = x + max_radius
+    y = y + max_radius
+
+    cdef list H = list()
+
+    for rad in radius:
+        # Accumulator
+        out = np.zeros((img.shape[0] + 2 * max_radius, img.shape[1] + 2 * max_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
+            out[tx, ty] += 1
+
+        if normalize:
+            out = out / len(circle_x)
+
+        H.append(out)
+    return np.array(H)
+
+
+@cython.boundscheck(False)
+def _hough_circle(np.ndarray img, \
+                  np.ndarray[ndim=1, dtype=np.npy_intp] radius, \
+                  normalize=True):
+
+    if img.ndim != 2:
+        raise ValueError('The input image must be 2D.')
+
+    # compute the nonzero indexes
+    cdef np.ndarray[ndim=1, dtype=np.npy_intp] x, y
+    y, x = np.PyArray_Nonzero(img)
+
+    # Offset the image
+    max_radius = radius.max()
+    x = x + max_radius
+    y = y + max_radius
+
+    cdef list H = list()
+
+    for rad in radius:
+        # Accumulator
+        out = np.zeros((img.shape[0] + 2 * max_radius, img.shape[1] + 2 * max_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
+            out[tx, ty] += 1
+
+        if normalize:
+            out = out / len(circle_x)
+
+        H.append(out)
+    return np.array(H)
+
+
 def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.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,29 @@ def hough(img, theta=None):
     """
     return _hough(img, theta)
 
+def hough_circle(img, radius, normalize=True):
+    """Perform a circle 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
+        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
+
+    Examples
+    --------
+
+    """
+    return _hough_circle(img, radius, normalize)
 
 def hough_peaks(hspace, angles, dists, min_distance=10, min_angle=10,
                 threshold=None, num_peaks=np.inf):