diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx index 662b3b88..6164f304 100644 --- a/skimage/transform/_hough_transform.pyx +++ b/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 ((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: diff --git a/skimage/transform/hough_transform.py b/skimage/transform/hough_transform.py index 05d5cfbe..88ef4b83 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,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):