diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index b8501258..f490eec7 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -132,7 +132,7 @@ - François Boulogne Drawing: Andres Method for circle perimeter, ellipse perimeter drawing, Bezier curve. - Circular Hough Transform + Circular and elliptical Hough Transforms Various fixes - Thouis Jones diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index 859c7515..5aab2700 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -1,4 +1,4 @@ -from ._hough_transform import (hough_circle, hough_line, +from ._hough_transform import (hough_circle, hough_ellipse, hough_line, probabilistic_hough_line) from .hough_transform import (hough, probabilistic_hough, hough_peaks, hough_line_peaks) @@ -15,6 +15,7 @@ from .pyramids import (pyramid_reduce, pyramid_expand, __all__ = ['hough_circle', + 'hough_ellipse', 'hough_line', 'probabilistic_hough_line', 'hough', diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx index d8a0b1bf..d158a84d 100644 --- a/skimage/transform/_hough_transform.pyx +++ b/skimage/transform/_hough_transform.pyx @@ -101,7 +101,129 @@ def hough_circle(cnp.ndarray img, return acc -def hough_line(cnp.ndarray img, cnp.ndarray[ndim=1, dtype=cnp.double_t] theta=None): +def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1, + int min_size=4, max_size=None): + """Perform an elliptical Hough transform. + + Parameters + ---------- + img : (M, N) ndarray + Input image with nonzero values representing edges. + threshold: int, optional (default 4) + Accumulator threshold value. + accuracy : double, optional (default 1) + Bin size on the minor axis used in the accumulator. + min_size : int, optional (default 4) + Minimal major axis length. + max_size : int, optional + Maximal minor axis length. (default None) + If None, the value is set to the half of the smaller + image dimension. + + Returns + ------- + res : list of tuples [(x0, y0, a, b, angle, accumulator)] + Where (x0, y0) is the center, (a, b) major and minor axis. + The angle value follows `draw.ellipse_perimeter()` convention. + + Examples + -------- + >>> img = np.zeros((25, 25), dtype=int) + >>> rr, cc = draw.ellipse_perimeter(10, 10, 6, 8) + >>> img[rr, cc] = 1 + >>> result = hough_ellipse(img, threshold=6) + [(10.0, 10.0, 8.0, 6.0474292058692187, 0.0, 8)] + + Notes + ----- + The accuracy must be chosen to produce a peak in the accumulator + distribution. In other words, a flat accumulator distribution with low + values may be caused by a too low bin size. + + References + ---------- + .. [1] Xie, Yonghong, and Qiang Ji. "A new efficient ellipse detection + method." Pattern Recognition, 2002. Proceedings. 16th International + Conference on. Vol. 2. IEEE, 2002 + """ + if img.ndim != 2: + raise ValueError('The input image must be 2D.') + + cdef long[:, :] pixels = np.transpose(np.nonzero(img)) + cdef Py_ssize_t num_pixels = pixels.shape[0] + cdef list acc = list() + cdef list results = list() + cdef bin_size = accuracy**2 + + cdef int max_b_squared + if max_size is None: + if img.shape[0] < img.shape[1]: + max_b_squared = np.round(0.5 * img.shape[0])**2 + else: + max_b_squared = np.round(0.5 * img.shape[1])**2 + else: + max_b_squared = max_size**2 + + cdef Py_ssize_t p1, p2, p3, p1x, p1y, p2x, p2y, p3x, p3y + cdef double x0, y0, a, b, d, k + cdef double cos_tau_squared, b_squared, f_squared, angle + + for p1 in range(num_pixels): + p1x = pixels[p1, 1] + p1y = pixels[p1, 0] + + for p2 in range(p1): + p2x = pixels[p2, 1] + p2y = pixels[p2, 0] + + # Candidate: center (x0, y0) and main axis a + a = 0.5 * sqrt((p1x - p2x)**2 + (p1y - p2y)**2) + if a > 0.5 * min_size: + x0 = 0.5 * (p1x + p2x) + y0 = 0.5 * (p1y + p2y) + + for p3 in range(num_pixels): + p3x = pixels[p3, 1] + p3y = pixels[p3, 0] + + d = sqrt((p3x - x0)**2 + (p3y - y0)**2) + if d > min_size: + f_squared = (p3x - p1x)**2 + (p3y - p1y)**2 + cos_tau_squared = ((a**2 + d**2 - f_squared) \ + / (2 * a * d))**2 + # Consider b2 > 0 and avoid division by zero + k = a**2 - d**2 * cos_tau_squared + if k > 0 and cos_tau_squared < 1: + b_squared = a**2 * d**2 * (1 - cos_tau_squared) / k + # b2 range is limited to avoid histogram memory + # overflow + if b_squared <= max_b_squared: + acc.append(b_squared) + + if len(acc) > 0: + bins = np.arange(0, np.max(acc) + bin_size, bin_size) + hist, bin_edges = np.histogram(acc, bins=bins) + hist_max = np.max(hist) + if hist_max > threshold: + angle = np.arctan2(p1x - p2x, p1y - p2y) + # pi - angle to keep ellipse_perimeter() convention + if angle != 0: + angle = np.pi - angle + b = sqrt(bin_edges[hist.argmax()]) + results.append((x0, + y0, + a, + b, + angle, + hist_max, # Accumulator + )) + acc = [] + + return results + + +def hough_line(cnp.ndarray img, + cnp.ndarray[ndim=1, dtype=cnp.double_t] theta=None): """Perform a straight line Hough transform. Parameters @@ -211,8 +333,8 @@ def probabilistic_hough_line(cnp.ndarray img, int threshold=10, Returns ------- lines : list - List of lines identified, lines in format ((x0, y0), (x1, y0)), indicating - line start and end. + List of lines identified, lines in format ((x0, y0), (x1, y0)), + indicating line start and end. References ---------- @@ -334,14 +456,14 @@ def probabilistic_hough_line(cnp.ndarray img, int threshold=10, y1 = py >> shift else: x1 = px >> shift - y1 = py; + y1 = py # check when line exits image boundary if x1 < 0 or x1 >= width or y1 < 0 or y1 >= height: break gap += 1 # if non-zero point found, continue the line if mask[y1, x1]: - gap = 0; + gap = 0 line_end[k, 1] = y1 line_end[k, 0] = x1 # if gap to this point was too large, end the line diff --git a/skimage/transform/tests/test_hough_transform.py b/skimage/transform/tests/test_hough_transform.py index 5dd47b15..344dbea0 100644 --- a/skimage/transform/tests/test_hough_transform.py +++ b/skimage/transform/tests/test_hough_transform.py @@ -2,7 +2,7 @@ import numpy as np from numpy.testing import * import skimage.transform as tf -from skimage.draw import circle_perimeter, line +from skimage.draw import line, circle_perimeter, ellipse_perimeter def append_desc(func, description): @@ -126,6 +126,7 @@ def test_hough_circle(): assert_equal(x[0], x_0) assert_equal(y[0], y_0) + def test_hough_circle_extended(): # Prepare picture # The circle center is outside the image @@ -133,7 +134,7 @@ def test_hough_circle_extended(): radius = 20 x_0, y_0 = (-5, 50) y, x = circle_perimeter(y_0, x_0, radius) - img[x[np.where(x>0)], y[np.where(x>0)]] = 1 + img[x[np.where(x > 0)], y[np.where(x > 0)]] = 1 out = tf.hough_circle(img, np.array([radius]), full_output=True) @@ -142,5 +143,40 @@ def test_hough_circle_extended(): assert_equal(x[0], x_0 + radius) assert_equal(y[0], y_0 + radius) + +def test_hough_ellipse_zero_angle(): + img = np.zeros((25, 25), dtype=int) + a = 6 + b = 8 + x0 = 12 + y0 = 12 + angle = 0 + rr, cc = ellipse_perimeter(x0, x0, b, a) + img[rr, cc] = 1 + result = tf.hough_ellipse(img, threshold=9) + assert_equal(result[0][0], x0) + assert_equal(result[0][1], y0) + assert_almost_equal(result[0][2], b, decimal=1) + assert_almost_equal(result[0][3], a, decimal=1) + assert_equal(result[0][4], angle) + + +def test_hough_ellipse_non_zero_angle(): + img = np.zeros((20, 20), dtype=int) + a = 6 + b = 9 + x0 = 10 + y0 = 10 + angle = np.pi/1.35 + rr, cc = ellipse_perimeter(x0, x0, b, a, orientation=angle) + img[rr, cc] = 1 + result = tf.hough_ellipse(img, threshold=15, accuracy=3) + print(result) + assert_almost_equal(result[0][0]/100., x0/100., decimal=1) + assert_almost_equal(result[0][1]/100., y0/100., decimal=1) + assert_almost_equal(result[0][2]/100., b/100., decimal=1) + assert_almost_equal(result[0][3]/100., a/100., decimal=1) + assert_almost_equal(result[0][4], angle, decimal=1) + if __name__ == "__main__": run_module_suite()