Implementing the Hough Transform
Originally posted by Dan on 06:00 Thu 24 February 2005, last modified 19:43 Tue 5 June 2007.
File under: 3rd year project image processing programming python
The Hough Transform (HT) is an image processing operation which enables the extraction of shapes, essentially lines from images. The principle is that there are an infinite number of line that pass though any point, each at a different orientation. The HT aims to determine which of those theoretical lines pass though most features in an image.
We represent using a polar co-ordinate system, therefore requiring two parameters r, and theta. These represent the length and angle of the normal to the theoretical line from the origin of the image (taken to be upper left corner). We determine which theoretical lines are actually most likely in the image, by systematically building up an array (2D) of 'votes'. The maximum value of this array will reveal the angle and length of the normal to the most likely line. This array is called the accumulator. Therefore the dimensions of the accumulator will be the maximum length, given by Pythagoras (and rounded to the nearest int), by 180 (for the range of theta).
The input to the HT is not a raw image, but one where some form of edge detection has operated on the original image. This operation can be done easily by standard modules in the Python Image Library or ImageMagick.
So the HT indexes every pixel in this edge-detected image. Assuming edges are white, the algorithm checks to see if the pixel intensity is 255. Then it proceeds to calculate all the possible lines in any orientation at that point. For each theoretical line the length is calculated, and if it is less than the maximum length and greater than zero, then the value of the accumulator for that r and angle value is incrememted.
After the accumulator array has been filled it's just a matter of finding the indices of the maximum value. These are our values of r and theta for the most likely edge in the image.
From here some geometry is required, so that the lines (normal, and edge) can be added to the edge data. The proceedure is shown in the 3 images below, and the code below them.
width, height = im.size
rmax = int(round(math.sqrt(height**2 + width**2)))
print rmax
acc = Numeric.zeros((rmax, 180))
for x,y in [(x,y) for x in range(1, width) for y in range(1, height)]:
if im.getpixel((x,y)) == 255:
for m in range(1, 180):
arg = (m*math.pi) / 180
r = int(round( (x * math.cos(arg)) + (y * math.sin(arg)) ))
if 0 < r < rmax:
acc[r][m] += 1
max = 0
for x,y in [(x,y) for x in range(rmax) for y in range(180)]:
if acc[x][y] > max:
max = acc[x][y];
r, theta = (x,y)
piTheta180 = (math.pi * theta)/180
cosTheta = math.cos(piTheta180)
sinTheta = math.sin(piTheta180)
if (r/sinTheta) > height:
m = -1 * math.tan(piTheta180)
c = int(round(r/cosTheta))
upper = (c,0)
lower = (int(round(abs((m*512)+c))),512)
else:
m = (-1.0) / math.tan(piTheta180)
c = int(round(r/sinTheta))
c2 = int(round(r/cosTheta))
upper = (c2,0)
lower = (0,c)
return (upper, lower), (r*cosTheta, r*sinTheta)