Initializations

> with(linalg):
with(plots):
Digits := 10: # Controls the precision of floating point arithmetic

Warning, new definition for norm

Warning, new definition for trace

The Burns Line Extraction Algorithm

This worksheet demonstrates the complete Brian Burns line segment extraction algorithm. The intent is to provide a complete blue print which explains the mathematics and provides a working, all be it, very slow working implementation.

Ross Beveridge March 24, 1999

Read a Raw PGM Image

Given an file name as an argument, this procedure returns an array representing an image stored in PGM raw format. Note this procedure is not robust. It assumes there are no comments in the file and that the raw data is stored one-byte-per-pixel.

> readImage := proc(imageFile)
local fd, A, i, j, type, rs, cs, bits:
fd := fopen(imageFile,READ);
type := fscanf(fd,"%s")[1];
rs := fscanf(fd,"%d")[1];
cs := fscanf(fd,"%d")[1];
bits := fscanf(fd,"%d")[1];
print("In Read Image ", type, rs, cs, bits):
A := matrix(cs,rs):
for i from 1 to cs do
for j from 1 to rs do
A[i,j] := readbytes(fd,1)[1]:
od:
od:
fclose(fd);
eval(A);
end:

>

Create a Test Image

Build a simple test image

> testImage := proc()
local i, j, width, height, A, noise:
width := 12: height := 12:
A := matrix(width,height):
noise := rand(0..10):

> for i from 1 to width do
for j from 1 to height do
A[i,j] := 100.0 + noise():
od:
od:
for i from 4 to 8 do
for j from 4 to 8 do
A[i,j] := A[i,j] + 100.0:
od:
od:
eval(A):
end:

Gradient estimation

The gradient masks

> Sh := matrix(3,3, [-1,-2,-1,0,0,0,1,2,1]):
Sv := transpose(Sh):
[Sh = eval(Sh), Sv = eval(Sv)];

The procedure to find Gradient Mag and Orientation

The procedure below convolves IM with the two Sobel Masks, Sh and Sv. An intermediate result is the approximation to the gradient in the horizontal and vertical directions, dx and dy. These are then used to compute the gradient magnitude and gradient orientation. Note here the true gradient magnitude is used. Often a city block approximation is used to save time.

> GradImages := proc(IM,Sh,Sv)
local rs, cs, Gh, Gv, Gm, Go, dx, dy, Ghf, Gvf, i, j, r, c, Pf;
description "Takes an image and two gradient masks, returns the gradient images.";
rs := rowdim(IM); cs := coldim(IM);
Pf := convert(Pi,float):
Gm := matrix(rs,cs,0):
Go := matrix(rs,cs,0):
Ghf := unapply(sum(sum(IM[r+i,c+j]*Sh[i+2,j+2],j=-1..1),i=-1..1),r,c);
Gvf := unapply(sum(sum(IM[r+i,c+j]*Sv[i+2,j+2],j=-1..1),i=-1..1),r,c):
for r from 2 to (rs-1) do
for c from 2 to (cs-1) do
dx := convert(Ghf(r,c),float):
dy := convert(Gvf(r,c),float):
Gm[r,c] := sqrt((dx*dx)+(dy*dy)):
Go[r,c] := arctan(dy,dx) + Pf:
od:
od:
[eval(Gm),eval(Go)]:
end:

Orientation Bucket Labels

Quantize the orientation into eight buckets. Adjacent pixels with the same bucket labels will subsequently be grouped into edge support regions.

Bucket Labeling

The bucket labels are integer between 1 and where is typically 8. They represent a quantization of the orientations. There is an offset that shifts the starting boundary of the first bucket off radian. The label is special, and indicates a pixel where the gradient magnitude is too small reliably yield a gradient orientation estimate. This is controlled by the gradient magnitude threshold .

> bucket := proc(o,g,k,t,magThresh)
local b;
if (g < magThresh) then
b := 0:
else
b := floor(((o+t)/(2*convert(Pi,float)))*k)+1:
fi:
# Special Case when o is 2PI instead of zero
if (b > k) then b := 1 fi:
b:
end:

> bucketLabel := proc(gIms,magThresh)
local Go, Gm, L1, L2, t1, t2, rs, cs, i, j:
# DEBUG;
Gm := gIms[1]: Go := gIms[2]:
rs := rowdim(Go): cs := coldim(Go):
L1 := matrix(rs,cs):
t1 := 0.0:
L2 := matrix(rs,cs):
t2 := (convert(Pi,float)*2.0)/16.0:
for i from 1 to rs do
for j from 1 to cs do
L1[i,j] := bucket(Go[i,j],Gm[i,j],8,t1,magThresh):
L2[i,j] := bucket(Go[i,j],Gm[i,j],8,t2,magThresh):
od:
od:
[eval(L1), eval(L2)]:
end:

Finding the Connected Components in the Bucket Orientation Images

In this step, each 4 connected edge support region is given a unique integer label. The labels run from 1 to . The algorithm used does not use any recursion. It is based upon repeated sweeps through the data propogating integer labels to adjacent pixels in each region. This algorithm can be slow under certain pathological configurations. However, these involve regions that form spirals, etc. These are nearly impossible configurations given the way that the edge support regions depend upon the underlying gradient orientation.

The first step builds two data structures. The first is a region label image where each pixel with a non-zero gradient orientation direction gets a unique integer label. The second is a vector of labels used to propogate the labels. The indices to this vector represent the original unique region labels. The value stored is the current label to be used for this region. So, if regions 1 and 2 are collapsed, position 2 will get the label 1. Initially, this vector just contains sequential integers starting at 1.

The second step is to take four sweeps through the image, going from left to right, top to bottom, right to left, and finally bottom to top, propogating region labels. If any labels change, set a flag indicating that a change has been made. This algorithm is then used repeatedly to update labels until the result is stable: i.e. no more changes happen.

There will be holes in the final set of labels. In other words, the integer labels will not be sequential starting at 1. To create regions with sequential labels, a pass is taken through the region label vector checking to see which ones have not been changed. In other words, which ones have a value equal to the position in the vector (index). It turns out that these are the final labels in the label image. Therefore, in going through the region label vector, a similar vector is made in which those entries which are actually used in the region labeling are assigned sequential integers. Then, this vector is used to update the actual region labels.

Initialize the label plane and region label sequence

> initLabels := proc(L)
local i, j, rs, cs, R, ri:
rs := rowdim(L): cs := coldim(L):
# Give each non zero pixel in label plane a unique region label.
R := matrix(rs,cs):
ri := 1:
for i from 1 to rs do
for j from 1 to cs do
if (L[i,j] > 0) then
R[i,j] := ri:
ri := ri + 1:
else
R[i,j] := 0:
fi:
od:
od:
eval(R)
end:

Test Adjoining Neighbors and Update

This procedure is passed in a pair of coordinates representing a pixel (a,b) and a pixel to the left or down (c,d). The values in the label plane are first checked to see they are both non zero. Next, they are checked to see if they have the same orientation label and different region labels. If all these tests are true, than these pixels are not currently in the same region but they should be in the same region. Thus, the values of the labels in the two pixels is checked and the lower of the two is propogated to the other.

> testNeighbors := proc(L,R,a,b,c,d)
local flag:
#print("Enter Test Neighbors", a, b, c, d):
#print("L and R", L, R):
#print(R[a,b],R[c,d],L[a,b],L[c,d]):
flag := false:
if (L[a,b] <> 0) and (L[c,d] <> 0) then
if (L[a,b] = L[c,d]) and (R[a,b] <> R[c,d]) then
flag := true:
if (R[a,b] < R[c,d]) then
R[c,d] := R[a,b]:
#print("Propogating", R[a,b], a,b,"to", R[c,d], c,d):
else
R[a,b] := R[c,d]:
#print("Propogating", R[c,d], c,d, "to", R[a,b], a,b):
fi:
fi:
fi:
flag:
end:

Propogate the labels along row and columns

Sweep from left to right, top to bottom, right to left and bottom to top, propogating the smaller of the two integer labels whenever a pair of adjacent pixels is found to be connected. Return false if no labels are propogated, otherwise return true. The intent is that this procedure is called repeatedly until it returns false, i.e. no more labels need to be propogated and the connected components are complete.

> upLabels := proc(L,R)
local i, j, rs, cs, flag:
rs := rowdim(R): cs := coldim(R):
flag := false:
# Go from the left to the right collapsing labels
for i from 1 to rs do
for j from 1 to (cs-1) do
if testNeighbors(L,R,i,j,i,j+1) then flag := true fi:
od:
od:
# Go from the top to the bottom collapsing labels
for j from 1 to cs do
for i from 2 to (rs-1) do
if testNeighbors(L,R,i,j,i+1,j) then flag := true fi:
od:
od:
# Go from the right to the left collapsing labels
for i from 1 to rs do
for j from cs by -1 to 2 do
if testNeighbors(L,R,i,j,i,j-1) then flag := true fi:
od:
od:
# Go from the bottom to the topt collapsing labels
for j from 1 to cs do
for i from (rs-1) by -1 to 2 do
if testNeighbors(L,R,i,j,i-1,j) then flag := true fi:
od:
od:
flag:
end:

Resequence the Region Labels

After labels have been propogated, there will be holes in the numbering. In other words, there may be a connected region labeled "1" and a region labeled "3", but no region labelled "2". The following procedure goes through a renumbers the regions sequentially starting with label "1".

> seqLabels := proc(R)
local C, N, k, i, j, rs, cs, ri, l:
rs := rowdim(R):
cs := coldim(R):
ri := rs * cs:
C := vector(ri,0):
N := vector(ri):
# Count how many times each non-zero label is used.
for i from 1 to rs do
for j from 1 to cs do
if (R[i,j] <> 0) then
l := R[i,j]:
C[l] := C[l] + 1:
fi:
od:
od:
# Old region label i goes to new region label k.
# Observe that the default zeroes in C are used to
# indicate an old label i which is not in use.
k := 1:
for i from 1 to ri do
if (C[i] > 0) then
N[i] := k:
k := k + 1:
fi:
od:
# Map in the new region labels by modifying the region label image R.
for i from 1 to rs do
for j from 1 to cs do
if (R[i,j] <> 0) then
R[i,j] := N[R[i,j]]:
fi:
od:
od:
eval(R):
end:

Connected Components Procdure

This procedure uses the procedure defined above. It first places initial labels in R. Next it propogates labels until all connected pixels have the same label. Finally, it resequences the labels so they begin at 1 and run up to k. There is a maxIter term to prevent endless looping. It is set to be the row dimension of the region label image: this should be enought to guarantee convergence of the region labels.

> connectedComponents := proc(L)
local R, flag, maxIter, i:
maxIter := rowdim(L):
i := 0:
R := initLabels(L):
flag := true:
while (flag and (i < maxIter)) do
flag := upLabels(L, R):
i := i + 1:
od:
seqLabels(R):
end:

Build One List of Edge Support Regions

There is a one-to-one correspondence between the final straight line segments produced by the Burns algorithm and the regions in the two edge support region images. One caveat, regions with fewer than a specified number of pixels are usually not considered further. Thus, in this step a master list of regions is created. What is stored in this list is the index of the region, its source image (either R1 or R2), the extents of the region and the size of the region.

Find the max value in an image

> imMaxVal := proc(A)
local i, j, max:
max := 0:
for i from 1 to rowdim(A) do
for j from 1 to coldim(A) do
if (A[i,j] > max) then max := A[i,j] fi:
od:
od:
max:
end:

Remove Small Regions

Go through the image once and determine which regions are smallerr than a specified threshold. Use a Region Label Vector to map from the old labels to the new labels, where the old regions below the specified size are mapped to the background - zero - region. Also, include an offset which is normally zero. However, it will allow the Burns algorithm to begin the region indexes for the the second label plane at one plus the max label of the first. It will become apparent that this simplifies later processing.

> remSmallRegions := proc(R,thresh,offset)
local i, j, max, C, RV, k:
# print("Regions Entering Small Regions", R):
max := imMaxVal(R): # The largest label used.
C := vector(max,0): # Count, how many pixels use a label
RV := vector(max,0): # New labels for regions above thresh in size
# Count the number of pixel in each region
for i from 1 to rowdim(R) do
for j from 1 to coldim(R) do
if (R[i,j] <> 0) then C[R[i,j]] := C[R[i,j]] + 1: fi:
od:
od:
# Initialize the new count, k, at the offset plus 1 and begin entering new labels for large regions.
# It is important that all entries for small regions already contain the label "0".
k := offset + 1:
for i from 1 to max do
if C[i] > thresh then
RV[i] := k:
k := k + 1:
fi:
od:
# Go back through giving new labels to all non zero labelled pixels.
for i from 1 to rowdim(R) do
for j from 1 to coldim(R) do
if (R[i,j] <> 0) then R[i,j] := RV[R[i,j]]: fi:
od:
od:
# print("Regions Leaving Small Regions", R):
R:
end:

Make the Regions List

Note that since this procedure modifies the region label planes, it also returns them. Specifically, this function returns a list with three elements. The first is the list of regions, i.e. the region index, source image, extents and size. Next is the first region label plane followed by the second region label plane.

> makeRegList := proc(Regs,minPixels)
local nRegs, R1, R2, R1n, R2n, RV, RR, RS, rs, cs, i, j:
rs := rowdim(Regs[1]): cs := coldim(Regs[1]):
R1 := remSmallRegions(Regs[1],minPixels,0):
R1n := imMaxVal(R1):
R2 := remSmallRegions(Regs[2],minPixels,R1n):
R2n := imMaxVal(R2):
RV := vector(R2n): # The extents of the region
RR := vector(R2n): # The label plane for region, 1 or 2
RS := vector(R2n,0): # The size of each region
for i from 1 to R2n do
RV[i] := [rs, cs, 1, 1]:
od:
# Create the label plane index, used later as a back pointer, for each region.
for i from 1 to R1n do
RR[i] := 1:
od:
for i from R1n + 1 to R2n do
RR[i] := 2:
od:
# Computer to size of the regions
for i from 1 to rs do
for j from 1 to cs do
if (R1[i,j] <> 0) then RS[R1[i,j]] := RS[R1[i,j]] + 1: fi:
if (R2[i,j] <> 0) then RS[R2[i,j]] := RS[R2[i,j]] + 1: fi:
od:
od:
# Go through all pixels in parallel, updating extents.
for i from 1 to rs do
for j from 1 to cs do
if R1[i,j] > 0 then
if RV[R1[i,j]][1] > i then RV[R1[i,j]][1] := i: fi:
if RV[R1[i,j]][2] > j then RV[R1[i,j]][2] := j: fi:
if RV[R1[i,j]][3] < i then RV[R1[i,j]][3] := i: fi:
if RV[R1[i,j]][4] < j then RV[R1[i,j]][4] := j: fi:
fi:
if R2[i,j] > 0 then
if RV[R2[i,j]][1] > i then RV[R2[i,j]][1] := i: fi:
if RV[R2[i,j]][2] > j then RV[R2[i,j]][2] := j: fi:
if RV[R2[i,j]][3] < i then RV[R2[i,j]][3] := i: fi:
if RV[R2[i,j]][4] < j then RV[R2[i,j]][4] := j: fi:
fi:
od:
od:
[[eval(RR), eval(RV), eval(RS)], eval(R1), eval(R2)]:
end:

>

Find Line for Each Region

For each region, fit a plane to the intensity surface of the image weighted by the gradient magnitude. Intersect this plane with a plane of constant intensity placed at the average intensity for the edge support region. Then clip the line to the extents box of the region.

Plane Fitting, the Mathematics

The plane fitting is done by using the explicit equation for the plane.

> eq3 := z = alpha*x + beta*y + gamma:
eq3;

>

The values of , and must minimize the sum of squared distances, measured along the z axis, between the height of the pixels and the z as a function of x and y. Here let A be the source image. Let G be the gradient magnitude image used to weight the contribution of the different pixels. Note that we will assume G is zero outside of each region. Thus, this equation is valid for the entire bounding box around the region.
The plane must minimize the following error function:

> unassign('i','j','x','y','err','A','G'):
err := G[i,j]*(alpha*i+beta*j+gamma-A[i,j])^2:
eq4 := E = Sum(Sum(err,j=y[min]..y[max]),i=x[min]..x[max]):
eq4;

Take the partial derivatives of this error with respect to the three plane parameters and set them equal to zero.

> DM := matrix(3,1, [diff(rhs(eq4),alpha),diff(rhs(eq4),beta),diff(rhs(eq4),gamma)]):
eq5 := eval(DM) = matrix(3,1, [0,0,0]):
eval(eq5);

Expand the resulting equation, taking the summations inside the matrix and moving the parameters outside.

> DM2 := matrix(3,1):
DM2[1,1] := alpha*Sum(Sum(G[i,j]*i*i,j=y[min]..y[max]),i=x[min]..x[max])
+ beta*Sum(Sum(G[i,j]*j*i,j=y[min]..y[max]),i=x[min]..x[max])
+ gamma*Sum(Sum(G[i,j]*i,j=y[min]..y[max]),i=x[min]..x[max])
- Sum(Sum(G[i,j]*A[i,j]*i,j=y[min]..y[max]),i=x[min]..x[max]):
DM2[2,1] := alpha*Sum(Sum(G[i,j]*i*j,j=y[min]..y[max]),i=x[min]..x[max])
+ beta*Sum(Sum(G[i,j]*j*j,j=y[min]..y[max]),i=x[min]..x[max])
+ gamma*Sum(Sum(G[i,j]*j,j=y[min]..y[max]),i=x[min]..x[max])
- Sum(Sum(G[i,j]*A[i,j]*j,j=y[min]..y[max]),i=x[min]..x[max]):
DM2[3,1] := alpha*Sum(Sum(G[i,j]*i,j=y[min]..y[max]),i=x[min]..x[max])
+ beta*Sum(Sum(G[i,j]*j,j=y[min]..y[max]),i=x[min]..x[max])
+ gamma*Sum(Sum(G[i,j],j=y[min]..y[max]),i=x[min]..x[max])
- Sum(Sum(G[i,j]*A[i,j],j=y[min]..y[max]),i=x[min]..x[max]):
eval(DM2);
MMij := matrix(3,3, [G[i,j]*i*i, G[i,j]*i*j, G[i,j]*i,
G[i,j]*i*j, G[i,j]*j*j, G[i,j]*j,
G[i,j]*i, G[i,j]*j, G[i,j]]):
VVij := matrix(3,1, [G[i,j]*A[i,j]*i, G[i,j]*A[i,j]*j, G[i,j]*A[i,j]]):
MM := Sum(Sum(eval(MMij),j=y[min]..y[max]),i=x[min]..x[max]):
VV := Sum(Sum(eval(VVij),j=y[min]..y[max]),i=x[min]..x[max]):
eval(MM) * matrix(3,1, [alpha,beta,gamma]) = eval(VV);
MMs := matrix(3,3, [a, b, c, b, d, e, c, e, h]):
VVs := matrix(3,1, [m, n, o]):
XXs := matrix(3,1, [alpha,beta,gamma]):
eval(MMs) * eval(XXs) = eval(VVs);

To solve for the plane parameters, invert the matrix. This can be done symbolically.

> MMsi := inverse(MMs):
Minv = eval(MMsi);
MMsin := scalarmul(MMsi,(-a*d*h+a*e^2+b^2*h-2*b*c*e+c^2*d)):
MMsid := -a*d*h+a*e^2+b^2*h-2*b*c*e+c^2*d:
Minv = eval(MMsin) / MMsid;

The x and y values given z and the other.

The following will come in handy below.

>

Solving for the plane

Plane parameters from M and B

Use the formula above to find the equation for the plane, i.e. solve for alpha, beta and gamma.

> planeFromMandB := proc(M,B)
local a, b, c, d, e, h, m, n, o, K, Mi, X, holdDigits:
holdDigits := Digits:
a := M[1,1]:
b := M[1,2]:
c := M[1,3]:
d := M[2,2]:
e := M[2,3]:
h := M[3,3]:
m := B[1,1]:
n := B[2,1]:
o := B[3,1]:
Mi := matrix(3,3):
K := -a*d*h+a*e^2+b^2*h-2*b*c*e+c^2*d:
# print("K in planeFromMandB ", K):
Mi[1,1] := (-d*h+e^2)/K:
Mi[1,2] := (b*h-c*e)/K:
Mi[1,3] := (-b*e+c*d)/K:
Mi[2,1] := Mi[1,2]:
Mi[2,2] := (-a*h+c^2)/K:
Mi[2,3] := (a*e-b*c)/K:
Mi[3,1] := Mi[1,3]:
Mi[3,2] := Mi[2,3]:
Mi[3,3] := (-a*d+b^2)/K:
Digits := holdDigits:
X := evalm(Mi &* B):
# print("Test Plane Solution: ", linsolve(M,B), "X ", X):
#X:
X:
end:

M and B matrix from regions

Compute the M matrix and the B column vector for the region indicated by the index ri

> regMandB := proc(RL, ri, Rs, Gs, A)
local M, B, lab, ext, L, G, i, j:
M := matrix(3,3, 0):
B := matrix(3,1, 0):
lab := RL[1][ri]: # the index of the label and gradient plane
ext := RL[2][ri]: # ext is the extends of the region
# print("Region ", ri, "Label ", lab, "Extents ", ext, "Size", RL[3][ri]):
for i from ext[1] to ext[3] do
for j from ext[2] to ext[4] do
L := Rs[lab]: # the label plane
G := Gs[1]: # the gradient plane
# This had been a test for greater than zero - an error.
if (L[i,j] = ri) then
M[1,1] := M[1,1] + (G[i,j] * i * i):
M[1,2] := M[1,2] + (G[i,j] * i * j):
M[1,3] := M[1,3] + (G[i,j] * i):
M[2,1] := M[1,2]:
M[2,2] := M[2,2] + (G[i,j] * j * j):
M[2,3] := M[2,3] + (G[i,j] * j):
M[3,1] := M[1,3]:
M[3,2] := M[2,3]:
M[3,3] := M[3,3] + G[i,j]:
B[1,1] := B[1,1] + (G[i,j] * A[i,j] * i):
B[2,1] := B[2,1] + (G[i,j] * A[i,j] * j):
B[3,1] := B[3,1] + (G[i,j] * A[i,j]):
fi:
od: # End of colum within extents
od: # End of row within extents
[eval(M), eval(B)]:
end:

Line from Plane & Extents

Represent the line in with the implicit form and then compute the for points where the infinite line crosses boundaries of the extents box.

Equations for x in terms of y and z, y in terms of x and z.

> x = solve(eq3,x);
y = solve(eq3,y);

Procedures for xofy, etc.

> xofy := proc(y,X,avg)
-(y*X[2,1]+X[3,1]-avg)/X[1,1]
end:
yofx := proc(x,X,avg)
-(x*X[1,1]+X[3,1]-avg)/X[2,1]
end:
zofxy := proc(x,y,X)
X[1,1] * x + X[2,1] * y + X[3,1]:
end:

3D Coords for Corners of Fit Plane

> planeCorners := proc(X,ext)
local xmin, ymin, xmax, ymax:
xmin := ext[1]: ymin := ext[2]: xmax := ext[3]: ymax := ext[4]:
[[xmin,ymin,zofxy(xmin,ymin,X)],
[xmin,ymax,zofxy(xmin,ymax,X)],
[xmax,ymax,zofxy(xmax,ymax,X)],
[xmax,ymin,zofxy(xmax,ymin,X)]]:
end:

Line Segment from Plane

Check first for vertical or horizontal special case. Otherwise, intersect infinite line with extents box and return the endpoints. Note, this code does not check for possible ties and thus intersection with more than two sides.

> lineFromPlane := proc(X,ext,avg)
local Line, xt, yt, xmin, ymin, xmax, ymax:
xmin := ext[1]: ymin := ext[2]: xmax := ext[3]: ymax := ext[4]:
if (X[1,1] = 0) then # Vertical line
Line := [[xmin,yofx(xmin,X,avg)], [xmax,yofx(xmax,X,avg)]]:
elif (X[2,1] = 0) then # Horizontal line
Line := [[xofy(ymin,X,avg),ymin],[xofy(ymax,X,avg)]]:
else
Line := []:
# Test if line crosses left boundary at ymin
xt := xofy(ymin,X,avg):
if ((xt > xmin) and (xt < xmax)) then
Line := [op(Line), [xt,ymin]]:
fi:
# Test if line crosses right boundary at ymax
xt := xofy(ymax,X,avg):
if ((xt > xmin) and (xt < xmax)) then
Line := [op(Line), [xt,ymax]]:
fi:
# Test if line crosses top boundary at xmin
yt := yofx(xmin,X,avg):
if ((yt > ymin) and (yt < ymax)) then
Line := [op(Line), [xmin,yt]]:
fi:
# Test if line crosses botom boundary at xmax
yt := yofx(xmax,X,avg):
if ((yt > ymin) and (yt < ymax)) then
Line := [op(Line), [xmax,yt]]:
fi:
fi:
if (nops(Line) <> 2) then
printf("Line to Extents Intersection Failed for Line\n"):
Line := [[0.0, 0.0], [0.0, 0.0]]:
fi:
Line:
end:

Voting for support

Each line segment, or equivalently, each edge support region, can be said to be supported by some pixels within it. To see what this means, realize that every pixel in the original image belongs to two edge support regions, one for each of the two label planes. A pixel "supports" the larger of the two regions to which it belongs. Thus, the support is the ratio of the number of pixels voting for a region over the total number of pixels in that region. To compute support, first the size of each region must be computed.

> lineSupport := proc(Size,len,Rs)
local i, j, rs, cs, r1, r2, s1, s2, Votes:
Votes := vector(len,0.0):
rs := rowdim(Rs[1]):
cs := coldim(Rs[1]):
for i from 1 to rs do
for j from 1 to cs do
r1 := Rs[1][i,j]: r2 := Rs[2][i,j]:
if (r1 > 0) or (r2 > 0) then
if (r1 > 0) then s1 := Size[r1]: else s1 := 0: fi:
if (r2 > 0) then s2 := Size[r2]: else s2 := 0: fi:
if (s1 > s2) then
Votes[r1] := Votes[r1] + 1:
else
Votes[r2] := Votes[r2] + 1:
fi:
fi:
od:
od:
for i from 1 to len do
Votes[i] := Votes[i] / Size[i]:
od:
Votes:
end:

>

>

Map Regions to Lines

This function takes in the list of edge support regions, the two region label planes, the gradient images and the source image. I generates a list of line segments. Each element in the list contains the 2D endpoints as clipped to the extents box, the support as computed by voting, the 3D line segment and 3D plane over the extents box for the purposes of visualization.

> line3D := proc(Line,avg)
# print("Entering Line 3D", Line, avg):
[[Line[1][1], Line[1][2], avg], [Line[2][1], Line[2][2], avg]]:
end:

> regToLines := proc(RL, Rs, Gs, A)
local len, ri, i, j, avg, Supp, M, B, MB, PL, PC, LN, RES:
# Setup to solve M*X=B, where X is the a, b and c for the plane
len := rhs(op(2,eval(RL[2]))):
# Accumulate results in vector RES
RES := vector(len):
# Compute support for regions
Supp := lineSupport(RL[3],len,Rs):
# For each of the regions in the Label Planes List
for ri from 1 to len do
MB := regMandB(RL, ri, Rs, Gs, A):
M := MB[1]: B := MB[2]: avg := B[3,1] / M[3,3]:
PL := planeFromMandB(M,B):
PC := planeCorners(PL,RL[2][ri]):
LN := lineFromPlane(PL,RL[2][ri],avg):
RES[ri] := [ri, LN, line3D(LN,avg), RL[2][ri], PC, RL[3][ri], Supp[ri]]:
od: #End of iteration through regions
RES:
end:

The Complete Algorithm

Here are the pieces put together. The input is a single image array. The output is a sequence of 2D line segments.

> burns := proc(A,magThresh,minRegSize)
local Sh, Sv, Gs, Ls, Rs,Rs1, Rs2, RGL:
Sh := matrix(3,3, [-1,-2,-1,0,0,0,1,2,1]):
Sv := transpose(Sh):
print("Computing Gradient Images"):
Gs := GradImages(A,Sh,Sv):
#print(A):
#print(Gs[1]):
print("Computing Bucket Labels"):
Ls := bucketLabel(Gs,magThresh):
print("Computing Connected Components on Labels 1"):
Rs1 := connectedComponents(Ls[1]):
print("Computing Connected Components on Labels 2"):
Rs2 := connectedComponents(Ls[2]):
Rs := [eval(Rs1), eval(Rs2)]:
print("Computing Edge Support Regions List"):
RGL := makeRegList(Rs,minRegSize):
print("Computing Line Segments for Regions");
Ls := regToLines(RGL[1], [RGL[2], RGL[3]], Gs,A):
eval(Ls):
end:

Test the Algorithm

> A := readImage(`CS 510, Spring 99:CS510 Current:Assignments:BurnsTest.pgm`):
# A := testImage():
# A := readImage(`CS 510, Spring 99:CS510 Current:Assignments:CagedSphereGrayscale.pgm`):

> LS := burns(A,5.0,5):

Write Line Segments to a File

After running the Burns Algorithm, the following procedure can read the sequence of line segments and write the endpoints to a text file

> writeLines := proc(LS,linesFile)
local i, len, fd:
len := rhs(op(2,eval(LS))):
fd := fopen(linesFile,WRITE);
for i from 1 to len do
if (LS[i][7] > 1/2) then
fprintf(fd, "%6.2f %6.2f %6.2f %6.2f\n", LS[i][2][1][1], LS[i][2][1][2], LS[i][2][2][1], LS[i][2][2][2]):
fi:
od:
fclose(fd):
end:

> # writeLines(LS,`CS 510, Spring 99:CS510 Current:Assignments:BurnsTest.dat`);

Plot Line Segments

Run the Burns Algorithm and draw the resulting segments

> drawLines2D := proc(LS,A)
local objs, len, i, rs, cs:
rs := rowdim(A): cs := coldim(A):
objs := {}:
len := rhs(op(2,eval(LS))):
for i from 1 to len do
if (LS[i][7] > 1/2) then
objs := {op(objs),CURVES(LS[i][2],THICKNESS(3),COLOR(HUE,0.4))}:
fi:
od:
display(objs,axes=FRAME,view=[0..rs,0..cs],scaling=CONSTRAINED,font=[HELVETICA,BOLD,14],labels=["Rows","Columns"]);
end:

> drawLines3D := proc(LS,A)
local objs, zVals, len, i, j, rs, cs:
rs := rowdim(A): cs := coldim(A):
objs := {}:
zVals := []:
len := rhs(op(2,eval(LS))):
for i from 1 to rs do
zVals := [op(zVals),[]]:
for j from 1 to cs do
zVals[i] := [op(zVals[i]),A[i,j]]:
od:
od:
objs := {op(objs),GRID(1..rs,1..cs,zVals,COLOR(ZGREYSCALE))}:
for i from 1 to len do
if (LS[i][7] > 1/2) then
objs := {op(objs),CURVES(LS[i][3],THICKNESS(10),COLOR(HUE,0.4))}:
fi:
od:
#print(objs):
display(objs,axes=FRAME,view=[0..rs,0..cs,0..255],font=[HELVETICA,BOLD,14]);
end:

> drawLines2D(LS,A);

> drawLines3D(LS,A);