'''Hough.py
An implementation of the Hough transform.
This program supports "Pixels, Numbers and Programs"
(chapter on Image Analysis)
1. read in an image.
A demo image containing strong lines has been included:
"floor.jpg"
To avoid overly long compute times, use the 100 x 67 version.
Otherwise, use the 200 x 134 version
2. compute edge strength.
We use the Roberts cross operator.
3. take Hough transform.
Note that the parameter space is based on a polar representation
of lines where the pole is at the center of the image. This means
a translation (shift) by half the width and half the height when using
PixelMath (x,y) coordinates in the line formula.
Instead of the standard rho = x cos theta + y sin theta
we use rho = (x - w/2) theta + (y - h/2) sin theta
4. suppress secondary maxima.
Two simple schemes are provided: (a) scan the Hough array and whenever
a value is max in its neighborhood, lower the neighbors by a factor of 0.5.
This has a disadvantage in making lots of local maxes by virtue of one
peak's neighborhood suppression makes it easier for the next non-neighbor
to win, because his neighbors have been suppressed.
(b) scan the Hough array and whenever a value is max in its neighborhood,
flag, but do not lower the values of, the neighbors. This is accomplished
with a separate array called NONPEAK whose values are initially all False,
and some of which get changed to True. After this is done, a local max is
a peak only if it is also not a nonpeak.
5. detect peaks.
Works in conjunction with the previous step.
6. report dominant lines
Indicates both the indexes (theta_idx and rho_idx) of the peak in the Hough
array, and the corresponding values of theta and rho.
Values of theta range from 0 to 2 pi (HN - 1)/HN, and
values of rho range from 0 to half_diag, which is the distance, in pixel
widths, from the center of the original image to a corner of it.
7. create composite of original with detected lines.
This is performed using a different version of the original, for aesthetic
reasons. Its resolution is double, and it is darkened by a factor of 0.5,
so that the white lines plotted on it will show up well.
Characteristics:
flexible input resolution (automatically detects width and height)
Hough array dimensions specified as global variables.
Thresholds for secondary max suppression specified as globals.
Hough array represented as Python array (?)
'''
ORIG = pmOpenTwo("floor100x67-monochrome.jpg")
EDGES = ORIG+1
pmSetDestination(EDGES)
pmSetFormula("sqrt(sqr(s1(x,y)-s1(x-1,y-1))+sqr(s1(x-1,y)-s1(x,y-1)))")
pmCompute()
# Dimensions of Hough array:
HM = 128 # number of different rho values
HN = 128 # number of different theta values
HRow = HM*[0.0] # an empty row of the Hough array.
H = [HRow[:] for i in range(HN)] # The whole, empty Hough array
EDGE_THRESH = 50 # Only pixels with edge value greater get to vote
HIT_THRESH = 1 # Pixels must lie within distance of 1 (pixel width)
# of a line for that line to receive the vote.
RAD = 3 # "radius" of neighborhood used in peak det. and supp.
# The neighborhoods, will be of size (2*RAD + 1)^2
# except at the borders of the image where they're smaller.
SUPRESSION_FACTOR = 0.5 # How much to reduce off-peak values.
HALF_DIAG = None # global that indicates range of rho values.
import math
def houghXform():
global EDGES, HALF_DIAG
w = pmGetImageWidth(EDGES)
h = pmGetImageHeight(EDGES)
diag = math.sqrt(w*w + h*h)
HALF_DIAG = diag / 2.0
for x_idx in range(1,w):
print ".", # show some progress.
x = x_idx - w/2.0
for y_idx in range(1,h):
y = y_idx - h/2.0
(r,g,b) = pmGetPixel(EDGES,x_idx,y_idx)
#monochr = (r+g+b)/3.0 # Use only if color orig.
monochr = r
if monochr > EDGE_THRESH:
for rho_idx in range(HM):
rho = HALF_DIAG * rho_idx/HM
for theta_idx in range(HN):
theta = math.pi * 2 * theta_idx / HN
if abs(x * math.cos(theta) + y * math.sin(theta) - rho)\
< HIT_THRESH:
H[theta_idx][rho_idx] += monochr
def showHough():
global H
temp = pmNewImage(0, "Hough array", HM, HN, 0, 0, 0)
pmPositionWindow(temp, 500, 500, 400, 400)
for i in range(2): pmZoom(temp, 200, 200)
for rho_idx in range(HM):
for theta_idx in range(HN):
r = g = b = int(H[theta_idx][rho_idx]/100)
pmSetPixel(temp, theta_idx, rho_idx, r, g, b)
HRow = HM*[False]
NOTPEAK = [HRow[:] for i in range(HN)]
def suppressSecondary():
for rho_idx in range(HM):
for theta_idx in range(HN):
if isLocalMax(rho_idx, theta_idx):
suppress_neighbors(rho_idx, theta_idx)
def isLocalMax(ri, ti):
global NOTPEAK
if NOTPEAK[ti][ri]: return False
startRI = max(0, ri-RAD)
endRI = min(HM, ri+RAD)
startTI = max(0, ti-RAD)
endTI = min(HN, ti+RAD)
refValue = H[ti][ri]
for rii in range(startRI, endRI):
for tii in range(startTI, endTI):
testValue = H[tii][rii]
if testValue > refValue: return False
return True
def suppress_neighbors(ri, ti):
startRI = max(0, ri-RAD)
endRI = min(HM, ri+RAD)
startTI = max(0, ti-RAD)
endTI = min(HN, ti+RAD)
refValue = H[ti][ri]
for rii in range(startRI, endRI):
for tii in range(startTI, endTI):
testValue = H[tii][rii]
if testValue < refValue:
NOTPEAK[tii][rii] = True
def findPeaks():
global HALF_DIAG, NOTPEAK
for rho_idx in range(HM):
for theta_idx in range(HN):
if NOTPEAK[theta_idx][rho_idx]: continue
if H[theta_idx][rho_idx] > 7000:
if isLocalMax(rho_idx, theta_idx):
rho = (rho_idx * HALF_DIAG) / HM
theta = theta_idx * 2 * math.pi / HN
print "Peak of strength "+str(int(H[theta_idx][rho_idx]/100))+\
" found at ("+str(theta_idx)+","+str(rho_idx)+");"+\
" theta = "+str(theta)+"; rho = "+str(rho)+"."
plotline(rho, theta)
PLOTWIN = None
BACKGDWIN = None
def plotline(rho, theta):
global PLOTWIN, BACKGDWIN, HALF_DIAG
if not PLOTWIN:
BACKGDWIN = pmOpenImage(0, "floor200x134-monochrome.jpg")
PLOTWIN = pmOpenImage(0, "floor200x134-monochrome.jpg")
HALF_DIAG = math.sqrt(pmGetImageWidth(PLOTWIN)**2 +\
pmGetImageHeight(PLOTWIN)**2)/2
pmSetSource1(BACKGDWIN)
pmSetDestination(PLOTWIN)
pmSetFormula("S1(x,y)/2")
pmCompute()
pmSetSource1(BACKGDWIN)
pmSetDestination(PLOTWIN)
pmSetFormula("dest(x,y)+if abs((x-w/2) * cos("+str(theta)+")+ "+\
"(y-h/2) * sin("+str(theta)+") - "+\
str(rho)+") < 0.8 then 255 else 0")
pmCompute()
houghXform()
showHough()
suppressSecondary()
findPeaks()
showHough()