wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-29 11:26:57 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #675 from sciunto/ellipseht_heapq_fix

Browse Source 
Hough transform ellipse parameter fix.
This commit is contained in:
Stefan van der Walt
2013-10-12 09:08:35 -07:00
parent 5dcb2d9fcc 6c67738808
commit c6a08d6a9c
4 changed files with 274 additions and 90 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/examples/plot_circular_elliptical_hough_transform.py
+21 -20
View File
@@ -74,7 +74,6 @@ for idx in np.argsort(accums)[::-1][:5]:
image[cy, cx] = (220, 20, 20)
ax.imshow(image, cmap=plt.cm.gray)
plt.show()
"""
@@ -96,13 +95,13 @@ an ellipse passes to them. A good match corresponds to high accumulator values.
A full description of the algorithm can be found in reference [1]_.
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
"""
import matplotlib.pyplot as plt
from skimage import data, filter, color
@@ -110,7 +109,7 @@ from skimage.transform import hough_ellipse
from skimage.draw import ellipse_perimeter
# Load picture, convert to grayscale and detect edges
image_rgb = data.load('coffee.png')[0:220, 100:450]
image_rgb = data.coffee()[0:220, 160:420]
image_gray = color.rgb2gray(image_rgb)
edges = filter.canny(image_gray, sigma=2.0,
low_threshold=0.55, high_threshold=0.8)
@@ -119,29 +118,31 @@ edges = filter.canny(image_gray, sigma=2.0,
# The accuracy corresponds to the bin size of a major axis.
# The value is chosen in order to get a single high accumulator.
# The threshold eliminates low accumulators
accum = hough_ellipse(edges, accuracy=10, threshold=170, min_size=50)
accum.sort(key=lambda x:x[5])
result = hough_ellipse(edges, accuracy=20, threshold=250,
min_size=100, max_size=120)
result.sort(order='accumulator')
# Estimated parameters for the ellipse
center_y = int(accum[-1][0])
center_x = int(accum[-1][1])
xradius = int(accum[-1][2])
yradius = int(accum[-1][3])
angle = np.pi - accum[-1][4]
best = result[-1]
yc = int(best[1])
xc = int(best[2])
a = int(best[3])
b = int(best[4])
orientation = best[5]
# Draw the ellipse on the original image
cx, cy = ellipse_perimeter(center_y, center_x,
yradius, xradius, orientation=angle)
image_rgb[cy, cx] = (0, 0, 1)
cy, cx = ellipse_perimeter(yc, xc, a, b, orientation)
image_rgb[cy, cx] = (0, 0, 255)
# Draw the edge (white) and the resulting ellipse (red)
edges = color.gray2rgb(edges)
edges[cy, cx] = (250, 0, 0)
fig = plt.subplots(figsize=(10, 6))
plt.subplot(1, 2, 1)
plt.title('Original picture')
plt.imshow(image_rgb)
plt.subplot(1, 2, 2)
plt.title('Edge (white) and result (red)')
plt.imshow(edges)
fig2, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(10, 6))
ax1.set_title('Original picture')
ax1.imshow(image_rgb)
ax2.set_title('Edge (white) and result (red)')
ax2.imshow(edges)
plt.show()
skimage/draw/_draw.pyx
+2 -2
View File
@@ -449,9 +449,9 @@ def ellipse_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t yradius,
----------
cy, cx : int
Centre coordinate of ellipse.
yradius, xradius: int
yradius, xradius : int
Minor and major semi-axes. ``(x/xradius)**2 + (y/yradius)**2 = 1``.
orientation: double, optional (default 0)
orientation : double, optional (default 0)
Major axis orientation in clockwise direction as radians.
Returns
skimage/transform/_hough_transform.pyx
+46 -37
View File
@@ -7,7 +7,7 @@ import numpy as np
cimport numpy as cnp
cimport cython
from libc.math cimport abs, fabs, sqrt, ceil
from libc.math cimport abs, fabs, sqrt, ceil, atan2, M_PI
from libc.stdlib cimport rand
from skimage.draw import circle_perimeter
@@ -122,17 +122,18 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1,
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.
result : ndarray with fields [(accumulator, y0, x0, a, b, orientation)]
Where ``(yc, xc)`` is the center, ``(a, b)`` the major and minor
axes, respectively. The `orientation` value follows
`skimage.draw.ellipse_perimeter` convention.
Examples
--------
>>> img = np.zeros((25, 25), dtype=int)
>>> rr, cc = draw.ellipse_perimeter(10, 10, 6, 8)
>>> img = np.zeros((25, 25), dtype=np.uint8)
>>> rr, cc = ellipse_perimeter(10, 10, 6, 8)
>>> img[cc, rr] = 1
>>> result = hough_ellipse(img, threshold=8)
[(10.0, 10.0, 8.0, 6.0, 0.0, 10)]
[(10, 10.0, 8.0, 6.0, 0.0, 10.0)]
Notes
-----
@@ -149,47 +150,47 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1,
if img.ndim != 2:
raise ValueError('The input image must be 2D.')
cdef Py_ssize_t[:, :] pixels = np.transpose(np.nonzero(img))
cdef Py_ssize_t num_pixels = pixels.shape[0]
cdef Py_ssize_t[:, ::1] pixels = np.row_stack(np.nonzero(img))
cdef Py_ssize_t num_pixels = pixels.shape[1]
cdef list acc = list()
cdef list results = list()
cdef bin_size = accuracy**2
cdef double 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
max_b_squared = np.round(0.5 * img.shape[0]) ** 2
else:
max_b_squared = np.round(0.5 * img.shape[1])**2
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
cdef double xc, yc, a, b, d, k
cdef double cos_tau_squared, b_squared, f_squared, orientation
for p1 in range(num_pixels):
p1x = pixels[p1, 1]
p1y = pixels[p1, 0]
p1x = pixels[1, p1]
p1y = pixels[0, p1]
for p2 in range(p1):
p2x = pixels[p2, 1]
p2y = pixels[p2, 0]
p2x = pixels[1, p2]
p2y = pixels[0, p2]
# Candidate: center (x0, y0) and main axis a
# Candidate: center (xc, yc) 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)
xc = 0.5 * (p1x + p2x)
yc = 0.5 * (p1y + p2y)
for p3 in range(num_pixels):
p3x = pixels[p3, 1]
p3y = pixels[p3, 0]
p3x = pixels[1, p3]
p3y = pixels[0, p3]
d = sqrt((p3x - x0)**2 + (p3y - y0)**2)
d = sqrt((p3x - xc)**2 + (p3y - yc)**2)
if d > min_size:
f_squared = (p3x - p1x)**2 + (p3y - p1y)**2
cos_tau_squared = ((a**2 + d**2 - f_squared) \
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
@@ -205,21 +206,29 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1,
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
orientation = atan2(p1x - p2x, p1y - p2y)
b = sqrt(bin_edges[hist.argmax()])
results.append((x0,
y0,
a,
b,
angle,
hist_max, # Accumulator
))
# to keep ellipse_perimeter() convention
if orientation != 0:
orientation = M_PI - orientation
# When orientation is not in [-pi:pi]
# it would mean in ellipse_perimeter()
# that a < b. But we keep a > b.
if orientation > M_PI:
orientation = orientation - M_PI / 2.
a, b = b, a
results.append((hist_max, # Accumulator
yc, xc,
a, b,
orientation))
acc = []
return results
return np.array(results, dtype=[('accumulator', np.intp),
('yc', np.double),
('xc', np.double),
('a', np.double),
('b', np.double),
('orientation', np.double)])
def hough_line(cnp.ndarray img,
skimage/transform/tests/test_hough_transform.py
+205 -31
View File
@@ -1,7 +1,5 @@
import numpy as np
from numpy.testing import (assert_almost_equal,
assert_equal,
)
from numpy.testing import assert_almost_equal, assert_equal
import skimage.transform as tf
from skimage.draw import line, circle_perimeter, ellipse_perimeter
@@ -81,8 +79,10 @@ def test_hough_line_peaks_dist():
img[:, 30] = True
img[:, 40] = True
hspace, angles, dists = tf.hough_line(img)
assert len(tf.hough_line_peaks(hspace, angles, dists, min_distance=5)[0]) == 2
assert len(tf.hough_line_peaks(hspace, angles, dists, min_distance=15)[0]) == 1
assert len(tf.hough_line_peaks(hspace, angles, dists,
min_distance=5)[0]) == 2
assert len(tf.hough_line_peaks(hspace, angles, dists,
min_distance=15)[0]) == 1
def test_hough_line_peaks_angle():
@@ -91,18 +91,24 @@ def test_hough_line_peaks_angle():
img[0, :] = True
hspace, angles, dists = tf.hough_line(img)
assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=45)[0]) == 2
assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=90)[0]) == 1
assert len(tf.hough_line_peaks(hspace, angles, dists,
min_angle=45)[0]) == 2
assert len(tf.hough_line_peaks(hspace, angles, dists,
min_angle=90)[0]) == 1
theta = np.linspace(0, np.pi, 100)
hspace, angles, dists = tf.hough_line(img, theta)
assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=45)[0]) == 2
assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=90)[0]) == 1
assert len(tf.hough_line_peaks(hspace, angles, dists,
min_angle=45)[0]) == 2
assert len(tf.hough_line_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_line(img, theta)
assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=45)[0]) == 2
assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=90)[0]) == 1
assert len(tf.hough_line_peaks(hspace, angles, dists,
min_angle=45)[0]) == 2
assert len(tf.hough_line_peaks(hspace, angles, dists,
min_angle=90)[0]) == 1
def test_hough_line_peaks_num():
@@ -149,36 +155,204 @@ def test_hough_circle_extended():
def test_hough_ellipse_zero_angle():
img = np.zeros((25, 25), dtype=int)
a = 6
b = 8
rx = 6
ry = 8
x0 = 12
y0 = 12
y0 = 15
angle = 0
rr, cc = ellipse_perimeter(x0, x0, b, a)
rr, cc = ellipse_perimeter(y0, x0, ry, rx)
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)
best = result[-1]
assert_equal(best[1], y0)
assert_equal(best[2], x0)
assert_almost_equal(best[3], ry, decimal=1)
assert_almost_equal(best[4], rx, decimal=1)
assert_equal(best[5], angle)
# Check if I re-draw the ellipse, points are the same!
# ie check API compatibility between hough_ellipse and ellipse_perimeter
rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]),
orientation=best[5])
assert_equal(rr, rr2)
assert_equal(cc, cc2)
def test_hough_ellipse_non_zero_angle():
img = np.zeros((20, 20), dtype=int)
a = 6
b = 9
def test_hough_ellipse_non_zero_posangle1():
# ry > rx, angle in [0:pi/2]
img = np.zeros((30, 24), dtype=int)
rx = 6
ry = 12
x0 = 10
y0 = 10
y0 = 15
angle = np.pi / 1.35
rr, cc = ellipse_perimeter(x0, x0, b, a, orientation=angle)
rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
img[rr, cc] = 1
result = tf.hough_ellipse(img, threshold=15, accuracy=3)
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)
result.sort(order='accumulator')
best = result[-1]
assert_almost_equal(best[1] / 100., y0 / 100., decimal=1)
assert_almost_equal(best[2] / 100., x0 / 100., decimal=1)
assert_almost_equal(best[3] / 10., ry / 10., decimal=1)
assert_almost_equal(best[4] / 100., rx / 100., decimal=1)
assert_almost_equal(best[5], angle, decimal=1)
# Check if I re-draw the ellipse, points are the same!
# ie check API compatibility between hough_ellipse and ellipse_perimeter
rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]),
orientation=best[5])
assert_equal(rr, rr2)
assert_equal(cc, cc2)
def test_hough_ellipse_non_zero_posangle2():
# ry < rx, angle in [0:pi/2]
img = np.zeros((30, 24), dtype=int)
rx = 12
ry = 6
x0 = 10
y0 = 15
angle = np.pi / 1.35
rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
img[rr, cc] = 1
result = tf.hough_ellipse(img, threshold=15, accuracy=3)
result.sort(order='accumulator')
best = result[-1]
assert_almost_equal(best[1] / 100., y0 / 100., decimal=1)
assert_almost_equal(best[2] / 100., x0 / 100., decimal=1)
assert_almost_equal(best[3] / 10., ry / 10., decimal=1)
assert_almost_equal(best[4] / 100., rx / 100., decimal=1)
assert_almost_equal(best[5], angle, decimal=1)
# Check if I re-draw the ellipse, points are the same!
# ie check API compatibility between hough_ellipse and ellipse_perimeter
rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]),
orientation=best[5])
assert_equal(rr, rr2)
assert_equal(cc, cc2)
def test_hough_ellipse_non_zero_posangle3():
# ry < rx, angle in [pi/2:pi]
img = np.zeros((30, 24), dtype=int)
rx = 12
ry = 6
x0 = 10
y0 = 15
angle = np.pi / 1.35 + np.pi / 2.
rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
img[rr, cc] = 1
result = tf.hough_ellipse(img, threshold=15, accuracy=3)
result.sort(order='accumulator')
best = result[-1]
# Check if I re-draw the ellipse, points are the same!
# ie check API compatibility between hough_ellipse and ellipse_perimeter
rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]),
orientation=best[5])
assert_equal(rr, rr2)
assert_equal(cc, cc2)
def test_hough_ellipse_non_zero_posangle4():
# ry < rx, angle in [pi:3pi/4]
img = np.zeros((30, 24), dtype=int)
rx = 12
ry = 6
x0 = 10
y0 = 15
angle = np.pi / 1.35 + np.pi
rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
img[rr, cc] = 1
result = tf.hough_ellipse(img, threshold=15, accuracy=3)
result.sort(order='accumulator')
best = result[-1]
# Check if I re-draw the ellipse, points are the same!
# ie check API compatibility between hough_ellipse and ellipse_perimeter
rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]),
orientation=best[5])
assert_equal(rr, rr2)
assert_equal(cc, cc2)
def test_hough_ellipse_non_zero_negangle1():
# ry > rx, angle in [0:-pi/2]
img = np.zeros((30, 24), dtype=int)
rx = 6
ry = 12
x0 = 10
y0 = 15
angle = - np.pi / 1.35
rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
img[rr, cc] = 1
result = tf.hough_ellipse(img, threshold=15, accuracy=3)
result.sort(order='accumulator')
best = result[-1]
# Check if I re-draw the ellipse, points are the same!
# ie check API compatibility between hough_ellipse and ellipse_perimeter
rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]),
orientation=best[5])
assert_equal(rr, rr2)
assert_equal(cc, cc2)
def test_hough_ellipse_non_zero_negangle2():
# ry < rx, angle in [0:-pi/2]
img = np.zeros((30, 24), dtype=int)
rx = 12
ry = 6
x0 = 10
y0 = 15
angle = - np.pi / 1.35
rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
img[rr, cc] = 1
result = tf.hough_ellipse(img, threshold=15, accuracy=3)
result.sort(order='accumulator')
best = result[-1]
# Check if I re-draw the ellipse, points are the same!
# ie check API compatibility between hough_ellipse and ellipse_perimeter
rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]),
orientation=best[5])
assert_equal(rr, rr2)
assert_equal(cc, cc2)
def test_hough_ellipse_non_zero_negangle3():
# ry < rx, angle in [-pi/2:-pi]
img = np.zeros((30, 24), dtype=int)
rx = 12
ry = 6
x0 = 10
y0 = 15
angle = - np.pi / 1.35 - np.pi / 2.
rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
img[rr, cc] = 1
result = tf.hough_ellipse(img, threshold=15, accuracy=3)
result.sort(order='accumulator')
best = result[-1]
# Check if I re-draw the ellipse, points are the same!
# ie check API compatibility between hough_ellipse and ellipse_perimeter
rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]),
orientation=best[5])
assert_equal(rr, rr2)
assert_equal(cc, cc2)
def test_hough_ellipse_non_zero_negangle4():
# ry < rx, angle in [-pi:-3pi/4]
img = np.zeros((30, 24), dtype=int)
rx = 12
ry = 6
x0 = 10
y0 = 15
angle = - np.pi / 1.35 - np.pi
rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle)
img[rr, cc] = 1
result = tf.hough_ellipse(img, threshold=15, accuracy=3)
result.sort(order='accumulator')
best = result[-1]
# Check if I re-draw the ellipse, points are the same!
# ie check API compatibility between hough_ellipse and ellipse_perimeter
rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]),
orientation=best[5])
assert_equal(rr, rr2)
assert_equal(cc, cc2)
if __name__ == "__main__":