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

Merge pull request #12 from holtzhau/hough

Browse Source 
Probabilistic Hough transform.
This commit is contained in:
Stefan van der Walt
2011-10-09 16:03:03 -07:00
parent 463e1d2f2d ed011f95d1
commit 12ab1a0aa4
4 changed files with 371 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/source/tutorials/hough_transform.txt
+133
View File
@@ -0,0 +1,133 @@
***************
Hough transform
***************
The Hough transform in its simplest form is a method to detect straight lines.
http://en.wikipedia.org/wiki/Hough_transform
As a first example we construct a line intersection.
.. ipython::
In [1]: import numpy as np
In [2]: from scikits.image.transform import hough, probabilistic_hough
In [3]: import matplotlib.pyplot as plt
In [4]: from matplotlib.lines import Line2D
In [5]: image = np.zeros((100, 100))
In [6]: for i in range(25, 75):
...: image[100 - i, i] = 255
...: image[i, i] = 255
...:
In [7]: plt.imshow(image)
@savefig hough_original.png width=4in
In [8]: plt.show()
The Hough transform converts the image into a parameter space that represents
lines. A line can be represented by the distance r of its closest point to the
origin and by the angle theta of this vector.
Every non-zero pixel of the image votes for potential line candidates, and the
local maxima represents the parameters of probable lines.
.. ipython::
In [9]: h, theta, d = hough(image)
In [10]: plt.figure()
In [10]: plt.title("hough transform")
In [10]: plt.xlabel("degrees")
In [10]: plt.ylabel("distance")
In [11]: plt.imshow(h)
@savefig hough_transform.png width=4in
In [12]: plt.show()
As can be seen, the maxima occur at 45 and 135 degrees, corresponding to the
normal vector angles of each line.
Another method is to use the function probabilistic_hough, an implementation
based on the Progressive Probabilistic Hough Transform [1]. It states that a
random subset of voting points give good enough results, and that lines can
be extracted during the voting process by walking along connected components.
This returns the beginning and end of line segments, which are useful.
The function has three parameters: a general threshold that is applied to
the Hough accumulator, a minimum line length and the line gap that influences
line merging.
.. ipython::
In [13]: lines = probabilistic_hough(image, threshold=10, line_length=10, line_gap=1)
In [14]: plt.figure()
In [15]: for line in lines:
....: p0, p1 = line
....: plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
....:
@savefig hough_probabilistic1.png width=4in
In [16]: plt.show()
The Hough transform are often used on edge detected images.
.. ipython::
In [17]: from scikits.image.io import imread
In [18]: from scikits.image import data_dir
In [19]: from scikits.image.filter import canny
In [20]: image = imread(data_dir + "/camera.png")
In [21]: edges = canny(image, 2, 1, 25)
In [22]: plt.imshow(edges)
@savefig hough_edge_detected.png width=4in
In [23]: plt.show()
Apply the Probabilistic Hough Transform and find lines longer than 10 with a
gap less than 3 pixels.
.. ipython::
In [24]: plt.figure()
In [25]: plt.imshow(np.zeros(edges.shape))
In [26]: lines = probabilistic_hough(edges, threshold=10, line_length=5, line_gap=3)
In [27]: for line in lines:
....: p0, p1 = line
....: plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
....:
@savefig hough_lines.png width=4in
In [28]: plt.show()
References
----------
.. [1] C. Galamhos, J. Matas and J. Kittler,"Progressive probabilistic Hough
transform for line detection", in IEEE Computer Society Conference on
Computer Vision and Pattern Recognition, 1999.
[2] Duda, R. O. and P. E. Hart, "Use of the Hough Transformation to Detect
Lines and Curves in Pictures," Comm. ACM, Vol. 15, pp. 11–15 (January,
1972)
scikits/image/transform/_hough_transform.pyx
+176 -12
View File
@@ -1,12 +1,15 @@
cimport cython
import numpy as np
cimport numpy as np
from random import randint
np.import_array()
cdef extern from "stdlib.h":
int rand()
cdef extern from "math.h":
int abs(int)
double fabs(double)
double sqrt(double)
double ceil(double)
double floor(double)
@@ -17,7 +20,6 @@ cdef double round(double val):
cdef double PI_2 = 1.5707963267948966
cdef double NEG_PI_2 = -PI_2
@cython.boundscheck(False)
def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
@@ -34,14 +36,14 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
ctheta = np.cos(theta)
stheta = np.sin(theta)
# compute the bins and allocate the output array
cdef np.ndarray[ndim=2, dtype=np.uint64_t] out
# compute the bins and allocate the accumulator array
cdef np.ndarray[ndim=2, dtype=np.uint64_t] accum
cdef np.ndarray[ndim=1, dtype=np.double_t] bins
cdef int max_distance, offset
max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] +
img.shape[1] * img.shape[1])))
out = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64)
accum = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64)
bins = np.linspace(-max_distance / 2.0, max_distance / 2.0, max_distance)
offset = max_distance / 2
@@ -49,17 +51,179 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
cdef np.ndarray[ndim=1, dtype=np.int_t] x_idxs, y_idxs
y_idxs, x_idxs = np.PyArray_Nonzero(img)
# finally, run the transform
cdef int nidxs, nthetas, i, j, x, y, out_idx
cdef int nidxs, nthetas, i, j, x, y, accum_idx
nidxs = y_idxs.shape[0] # x and y are the same shape
nthetas = theta.shape[0]
for i in range(nidxs):
x = x_idxs[i]
y = y_idxs[i]
y = y_idxs[i]
for j in range(nthetas):
out_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
out[out_idx, j] += 1
accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
accum[accum_idx, j] += 1
return accum, theta, bins
return out, theta, bins
import math
@cython.cdivision(True)
@cython.boundscheck(False)
def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \
int line_gap, np.ndarray[ndim=1, dtype=np.double_t] theta=None):
if img.ndim != 2:
raise ValueError('The input image must be 2D.')
# compute the array of angles and their sine and cosine
cdef np.ndarray[ndim=1, dtype=np.double_t] ctheta
cdef np.ndarray[ndim=1, dtype=np.double_t] stheta
# calculate thetas if none specified
if theta is None:
theta = np.linspace(math.pi/2, -math.pi/2, 180)
theta = math.pi/2-np.arange(180)/180.0* math.pi
ctheta = np.cos(theta)
stheta = np.sin(theta)
cdef int height = img.shape[0]
cdef int width = img.shape[1]
# compute the bins and allocate the accumulator array
cdef np.ndarray[ndim=2, dtype=np.int64_t] accum
cdef np.ndarray[ndim=2, dtype=np.uint8_t] mask = np.zeros((height, width), dtype=np.uint8)
cdef np.ndarray[ndim=2, dtype=np.int32_t] line_end = np.zeros((2, 2), dtype=np.int32)
cdef int max_distance, offset, num_indexes, index
cdef double a, b
cdef int nidxs, nthetas, i, j, x, y, px, py, accum_idx, value, max_value, max_theta
cdef int shift = 16
# maximum line number cutoff
cdef int lines_max = 2 ** 15
cdef int xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, good_line, count
max_distance = 2 * <int>ceil((sqrt(img.shape[0] * img.shape[0] +
img.shape[1] * img.shape[1])))
accum = np.zeros((max_distance, theta.shape[0]), dtype=np.int64)
offset = max_distance / 2
# find the nonzero indexes
cdef np.ndarray[ndim=1, dtype=np.int_t] x_idxs, y_idxs
y_idxs, x_idxs = np.nonzero(img)
num_indexes = y_idxs.shape[0] # x and y are the same shape
nthetas = theta.shape[0]
points = []
for i in range(num_indexes):
points.append((x_idxs[i], y_idxs[i]))
lines = []
# create mask of all non-zero indexes
for i in range(num_indexes):
mask[y_idxs[i], x_idxs[i]] = 1
while 1:
# select random non-zero point
count = len(points)
if count == 0:
break
index = rand() % (count)
x = points[index][0]
y = points[index][1]
del points[index]
# if previously eliminated, skip
if not mask[y, x]:
continue
value = 0
max_value = value_threshold-1
max_theta = -1
# apply hough transform on point
for j in range(nthetas):
accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
accum[accum_idx, j] += 1
value = accum[accum_idx, j]
if value > max_value:
max_value = value
max_theta = j
if max_value < value_threshold:
continue
# from the random point walk in opposite directions and find line beginning and end
a = -stheta[max_theta]
b = ctheta[max_theta]
x0 = x
y0 = y
# calculate gradient of walks using fixed point math
xflag = fabs(a) > fabs(b)
if xflag:
if a > 0:
dx0 = 1
else:
dx0 = -1
dy0 = <int>round(b * (1 << shift) / fabs(a))
y0 = (y0 << shift) + (1 << (shift - 1))
else:
if b > 0:
dy0 = 1
else:
dy0 = -1
dx0 = <int>round(a * (1 << shift) / fabs(b))
x0 = (x0 << shift) + (1 << (shift - 1))
# pass 1: walk the line, merging lines less than specified gap length
for k in range(2):
gap = 0
px = x0
py = y0
dx = dx0
dy = dy0
if k > 0:
dx = -dx
dy = -dy
while 1:
if xflag:
x1 = px
y1 = py >> shift
else:
x1 = px >> shift
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;
line_end[k, 1] = y1
line_end[k, 0] = x1
# if gap to this point was too large, end the line
elif gap > line_gap:
break
px += dx
py += dy
# confirm line length is sufficient
good_line = abs(line_end[1, 1] - line_end[0, 1]) >= line_length or \
abs(line_end[1, 0] - line_end[0, 0]) >= line_length
# pass 2: walk the line again and reset accumulator and mask
for k in range(2):
px = x0
py = y0
dx = dx0
dy = dy0
if k > 0:
dx = -dx
dy = -dy
while 1:
if xflag:
x1 = px
y1 = py >> shift
else:
x1 = px >> shift
y1 = py
# if non-zero point found, continue the line
if mask[y1, x1]:
if good_line:
accum_idx = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset
accum[accum_idx, max_theta] -= 1
mask[y1, x1] = 0
# exit when the point is the line end
if x1 == line_end[k, 0] and y1 == line_end[k, 1]:
break
px += dx
py += dy
# add line to the result
if good_line:
lines.append(((line_end[0, 0], line_end[0, 1]), (line_end[1, 0], line_end[1, 1])))
if len(lines) > lines_max:
return lines
return lines
scikits/image/transform/hough_transform.py
+38 -2
View File
@@ -1,7 +1,8 @@
__all__ = ['hough']
__all__ = ['hough', 'probabilistic_hough']
from itertools import izip
import numpy as np
from _hough_transform import _probabilistic_hough
def _hough(img, theta=None):
if img.ndim != 2:
@@ -58,6 +59,39 @@ except ImportError:
pass
def probabilistic_hough(img, threshold=10, line_length=50, line_gap=10, theta=None):
"""Performs a progressive probabilistic line Hough transform and returns the detected lines.
Parameters
----------
img : (M, N) ndarray
Input image with nonzero values representing edges.
value_threshold: int
Threshold
line_length: int, optional (default 50)
Minimum accepted length of detected lines.
Increase the parameter to extract longer lines.
line_gap: int, optional, (default 10)
Maximum gap between pixels to still form a line.
Increase the parameter to merge broken lines more aggresively.
theta :1D ndarray, dtype=double, optional, default (-pi/2 .. pi/2)
Angles at which to compute the transform, in radians.
Returns
-------
lines : list
List of lines identified, lines in format ((x0, y0), (x1, y0)), indicating
line start and end.
References
----------
.. [1] C. Galamhos, J. Matas and J. Kittler,"Progressive probabilistic Hough
transform for line detection", in IEEE Computer Society Conference on
Computer Vision and Pattern Recognition, 1999.
"""
return _probabilistic_hough(img, threshold, line_length, line_gap, theta)
def hough(img, theta=None):
"""Perform a straight line Hough transform.
@@ -67,7 +101,7 @@ def hough(img, theta=None):
Input image with nonzero values representing edges.
theta : 1D ndarray of double
Angles at which to compute the transform, in radians.
Defaults to -pi/2 - pi/2
Defaults to -pi/2 .. pi/2
Returns
-------
@@ -106,3 +140,5 @@ def hough(img, theta=None):
"""
return _hough(img, theta)
scikits/image/transform/tests/test_hough_transform.py
+24
View File
@@ -3,6 +3,7 @@ from numpy.testing import *
import scikits.image.transform as tf
import scikits.image.transform.hough_transform as ht
from scikits.image.transform import probabilistic_hough
def append_desc(func, description):
"""Append the test function ``func`` and append
@@ -12,6 +13,8 @@ def append_desc(func, description):
return func
from scikits.image.transform import *
def test_hough():
# Generate a test image
img = np.zeros((100, 100), dtype=int)
@@ -27,6 +30,7 @@ def test_hough():
assert_equal(dist > 70, dist < 72)
assert_equal(theta > 0.78, theta < 0.79)
def test_hough_angles():
img = np.zeros((10, 10))
img[0, 0] = 1
@@ -43,6 +47,26 @@ def test_py_hough():
tf._hough = fast_hough
def test_probabilistic_hough():
# Generate a test image
img = np.zeros((100, 100), dtype=int)
for i in range(25, 75):
img[100 - i, i] = 100
img[i, i] = 100
# decrease default theta sampling because similar orientations may confuse
# as mentioned in article of Galambos et al
theta=np.linspace(0, np.pi, 45)
lines = probabilistic_hough(img, theta=theta, threshold=10, line_length=10, line_gap=1)
# sort the lines according to the x-axis
sorted_lines = []
for line in lines:
line = list(line)
line.sort(lambda x,y: cmp(x[0], y[0]))
sorted_lines.append(line)
assert([(25, 75), (74, 26)] in sorted_lines)
assert([(25, 25), (74, 74)] in sorted_lines)
if __name__ == "__main__":
run_module_suite()