{VERSION 3 0 "SUN SPARC SOLARIS" "3.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }{PSTYLE " Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Text Output" -1 2 1 {CSTYLE "" -1 -1 "Courier" 1 10 0 0 255 1 0 0 0 0 0 1 3 0 3 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 8 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading \+ 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 8 2 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 3" 4 5 1 {CSTYLE "" -1 -1 "" 1 12 0 0 0 0 1 0 0 0 0 0 0 0 0 }0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 }{PSTYLE " Warning" 2 7 1 {CSTYLE "" -1 -1 "" 0 1 0 0 255 1 0 0 0 0 0 0 1 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" 0 11 1 {CSTYLE " " -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Plot" 0 13 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 4" 5 20 1 {CSTYLE "" -1 -1 "" 1 10 0 0 0 0 1 0 0 0 0 0 0 0 0 }0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 }} {SECT 0 {SECT 0 {PARA 3 "" 0 "" {TEXT -1 15 "Initializations" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "with(linalg):\nwith(plots): \nDigits := 10: # Controls the precision of floating point arithmeti c" }}{PARA 7 "" 1 "" {TEXT -1 32 "Warning, new definition for norm" }} {PARA 7 "" 1 "" {TEXT -1 33 "Warning, new definition for trace" }}}} {SECT 0 {PARA 3 "" 0 "" {TEXT -1 35 "The Burns Line Extraction Algorit hm" }}{PARA 0 "" 0 "" {TEXT -1 235 "This worksheet demonstrates the co mplete Brian Burns line segment extraction algorithm. The intent is to provide a complete blue print which explains the mathematics and prov ides a working, all be it, very slow working implementation. " }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 37 "Ross Beve ridge March 24, 1999" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 7 "Updates" }}{PARA 0 "" 0 "" {TEXT -1 265 "4/16/99: I have fixed a bug in regMandB which computes the M matrix \+ and B vectors for each region. It was NOT checking that the region lab el matched. This was a serious bug and caused pixels within the extent s bug and outside the region to influence the plane fit." }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 20 "Read a Raw PGM Image" }}{PARA 0 "" 0 "" {TEXT -1 248 "Given an file name as an argument, this procedure return s an array representing an image stored in PGM raw format. Note this p rocedure is not robust. It assumes there are no comments in the file a nd that the raw data is stored one-byte-per-pixel. " }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 442 "readImage := proc(imageFile)\n local fd, A, i, j, type, rs, cs, bits:\n fd \+ := fopen(imageFile,READ); \n type := fscanf(fd,\"%s\")[1];\n rs \+ := fscanf(fd,\"%d\")[1];\n cs := fscanf(fd,\"%d\")[1];\n bits \+ := fscanf(fd,\"%d\")[1];\n print(\"In Read Image \", type, rs, cs, b its):\n A := matrix(cs,rs):\n for i from 1 to cs do\n for j f rom 1 to rs do\n A[i,j] := readbytes(fd,1)[1]:\n od:\n \+ od:\n fclose(fd);\n eval(A);\nend:\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 19 "Create a Te st Image" }}{PARA 0 "" 0 "" {TEXT -1 25 "Build a simple test image" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 144 "testImage := proc()\n local i, j, width, height, A, noise:\n \+ width := 12: height := 12:\n A := matrix(width,height):\n noise := rand(0..10):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 232 " for i from 1 t o width do\n for j from 1 to height do\n A[i,j] := 100.0 \+ + noise():\n od:\n od:\n for i from 4 to 8 do\n for j fr om 4 to 8 do\n A[i,j] := A[i,j] + 100.0:\n od:\n od:\n \+ eval(A):\nend:" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 19 "Gradient es timation" }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 18 "The gradient masks" } }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 95 "Sh := matrix(3,3, [-1,-2,-1 ,0,0,0,1,2,1]):\nSv := transpose(Sh):\n[Sh = eval(Sh), Sv = eval(Sv)]; " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$/%#ShG-%'matrixG6#7%7%!\"\"!\"#F +7%\"\"!F.F.7%\"\"\"\"\"#F0/%#SvG-F'6#7%7%F+F.F07%F,F.F1F7" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 50 "The procedure to find Gradient Mag and \+ Orientation" }}{PARA 0 "" 0 "" {TEXT -1 368 " The procedure below con volves IM with the two Sobel Masks, Sh and Sv. An intermediate result \+ is the approximation to the gradient in the horizontal and vertical di rections, dx and dy. These are then used to compute the gradient magni tude and gradient orientation. Note here the true gradient magnitude i s used. Often a city block approximation is used to save time. " }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 718 "GradImages := proc(IM,Sh,Sv )\n local rs, cs, Gh, Gv, Gm, Go, dx, dy, Ghf, Gvf, i, j, r, c, Pf; \n description \"Takes an image and two gradient masks, returns the \+ gradient images.\";\n rs := rowdim(IM); cs := coldim(IM);\n Pf := \+ convert(Pi,float):\n Gm := matrix(rs,cs,0):\n Go := matrix(rs,cs,0 ):\n Ghf := unapply(sum(sum(IM[r+i,c+j]*Sh[i+2,j+2],j=-1..1),i=-1..1 ),r,c);\n Gvf := unapply(sum(sum(IM[r+i,c+j]*Sv[i+2,j+2],j=-1..1),i= -1..1),r,c):\n for r from 2 to (rs-1) do\n for c from 2 to (cs- 1) do\n dx := convert(Ghf(r,c),float):\n dy := convert (Gvf(r,c),float):\n Gm[r,c] := sqrt((dx*dx)+(dy*dy)):\n \+ Go[r,c] := arctan(dy,dx) + Pf:\n od:\n od:\n [eval(Gm),eval (Go)]:\nend:" }}}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 25 "Orientation Bu cket Labels" }}{PARA 0 "" 0 "" {TEXT -1 145 "Quantize the orientation \+ into eight buckets. Adjacent pixels with the same bucket labels will \+ subsequently be grouped into edge support regions." }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 15 "Bucket Labeling" }}{PARA 0 "" 0 "" {TEXT -1 44 "T he bucket labels are integer between 1 and " }{XPPEDIT 18 0 "k;" "6#% \"kG" }{TEXT -1 7 " where " }{XPPEDIT 18 0 "k;" "6#%\"kG" }{TEXT -1 88 " is typically 8. They represent a quantization of the orientations . There is an offset " }{XPPEDIT 18 0 "t;" "6#%\"tG" }{TEXT -1 59 " t hat shifts the starting boundary of the first bucket off " }{XPPEDIT 18 0 "0;" "6#\"\"!" }{TEXT -1 19 " radian. The label " }{XPPEDIT 18 0 "0;" "6#\"\"!" }{TEXT -1 181 " is special, and indicates a pixel where the gradient magnitude is too small reliably yield a gradient orienta tion estimate. This is controlled by the gradient magnitude threshold " }{XPPEDIT 18 0 "m;" "6#%\"mG" }{TEXT -1 1 "." }}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 244 "bucket := proc(o,g,k,t,magThresh)\n local b ;\n if (g < magThresh) then \n b := 0:\n else\n b := flo or(((o+t)/(2*convert(Pi,float)))*k)+1:\n fi:\n # Special Case when o is 2PI instead of zero\n if (b > k) then b := 1 fi:\n b:\nend: " }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 492 "bucketLabel := proc(g Ims,magThresh)\n local Go, Gm, L1, L2, t1, t2, rs, cs, i, j:\n # D EBUG;\n Gm := gIms[1]: Go := gIms[2]:\n rs := rowdim(Go): cs := coldim(Go):\n L1 := matrix(rs,cs):\n t1 := 0.0:\n L2 := matrix( rs,cs):\n t2 := (convert(Pi,float)*2.0)/16.0:\n for i from 1 to rs do\n for j from 1 to cs do\n L1[i,j] := bucket(Go[i,j], Gm[i,j],8,t1,magThresh):\n L2[i,j] := bucket(Go[i,j],Gm[i,j], 8,t2,magThresh):\n od:\n od:\n [eval(L1), eval(L2)]:\nend:" } }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 65 "Finding the Connected Componen ts in the Bucket Orientation Images" }}{PARA 0 "" 0 "" {TEXT -1 109 "I n this step, each 4 connected edge support region is given a unique in teger label. The labels run from 1 to " }{XPPEDIT 18 0 "k;" "6#%\"kG" }{TEXT -1 422 ". 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 c ertain pathological configurations. However, these involve regions tha t form spirals, etc. These are nearly impossible configurations given \+ the way that the edge support regions depend upon the underlying gradi ent orientation. " }}{PARA 0 "" 0 "" {TEXT -1 514 "The first step buil ds two data structures. The first is a region label image where each p ixel with a non-zero gradient orientation direction gets a unique inte ger label. The second is a vector of labels used to propogate the lab els. The indices to this vector represent the original unique region l abels. The value stored is the current label to be used for this regio n. 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." }}{PARA 0 "" 0 "" {TEXT -1 354 "The second step is to take four sweeps through the image, going from left to right, top to bottom, ri ght 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 resu lt is stable: i.e. no more changes happen. " }}{PARA 0 "" 0 "" {TEXT -1 655 "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 r egions with sequential labels, a pass is taken through the region labe l vector checking to see which ones have not been changed. In other wo rds, which ones have a value equal to the position in the vector (inde x). It turns out that these are the final labels in the label image. T herefore, in going through the region label vector, a similar vector i s made in which those entries which are actually used in the region la beling are assigned sequential integers. Then, this vector is used to \+ update the actual region labels. " }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 52 "Initialize the label plane and region label sequence" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 411 "initLabels := proc(L)\n local i, j, rs, cs, R, ri:\n rs := rowdim(L): cs := coldim(L):\n # Give ea ch non zero pixel in label plane a unique region label.\n R := matri x(rs,cs):\n ri := 1:\n for i from 1 to rs do\n for j from 1 t o cs do\n if (L[i,j] > 0) then\n R[i,j] := ri:\n \+ ri := ri + 1:\n else\n R[i,j] := 0:\n \+ fi:\n od:\n od:\n eval(R)\nend:" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 35 "Test Adjoining Neighbors and Update" }}{PARA 0 "" 0 "" {TEXT -1 530 "This procedure is passed in a pair of coordinates rep resenting a pixel (a,b) and a pixel to the left or down (c,d). The val ues 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 labe l 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 chec ked and the lower of the two is propogated to the other. " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 598 "testNeighbors := proc(L,R,a,b,c,d) \n local flag:\n #print(\"Enter Test Neighbors\", a, b, c, d):\n \+ #print(\"L and R\", L, R):\n #print(R[a,b],R[c,d],L[a,b],L[c,d]):\n flag := false:\n if (L[a,b] <> 0) and (L[c,d] <> 0) then\n i f (L[a,b] = L[c,d]) and (R[a,b] <> R[c,d]) then\n flag := true :\n if (R[a,b] < R[c,d]) then\n R[c,d] := R[a,b]: \+ \n #print(\"Propogating\", R[a,b], a,b,\"to\", R[c,d], c,d) :\n else\n R[a,b] := R[c,d]: \n #print (\"Propogating\", R[c,d], c,d, \"to\", R[a,b], a,b):\n fi:\n \+ fi:\n fi:\n flag:\nend: " }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 42 "Propogate the labels along row and columns" }}{PARA 0 "" 0 "" {TEXT -1 414 "Sweep from left to right, top to bottom, right to l eft and bottom to top, propogating the smaller of the two integer labe ls whenever a pair of adjacent pixels is found to be connected. Retur n 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." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 889 "upLabels : = proc(L,R)\n local i, j, rs, cs, flag:\n rs := rowdim(R): cs := \+ coldim(R):\n flag := false:\n # Go from the left to the right coll apsing labels\n for i from 1 to rs do\n for j from 1 to (cs-1) \+ do\n if testNeighbors(L,R,i,j,i,j+1) then flag := true fi:\n \+ od:\n od:\n # Go from the top to the bottom collapsing labels \n for j from 1 to cs do\n for i from 2 to (rs-1) do\n \+ if testNeighbors(L,R,i,j,i+1,j) then flag := true fi:\n od:\n o d:\n # Go from the right to the left collapsing labels\n for i fro m 1 to rs do\n for j from cs by -1 to 2 do\n if testNeigh bors(L,R,i,j,i,j-1) then flag := true fi:\n od:\n od:\n # Go \+ from the bottom to the topt collapsing labels\n for j from 1 to cs d o\n for i from (rs-1) by -1 to 2 do\n if testNeighbors(L, R,i,j,i-1,j) then flag := true fi:\n od:\n od:\n flag:\nend: " }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 28 "Resequence the Region Label s" }}{PARA 0 "" 0 "" {TEXT -1 287 " 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 reg ion labelled \"2\". The following procedure goes through a renumbers t he regions sequentially starting with label \"1\"." }}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 897 "seqLabels := proc(R)\n local C, N, k, i, \+ j, rs, cs, ri, l:\n rs := rowdim(R): \n cs := coldim(R):\n ri : = rs * cs:\n C := vector(ri,0):\n N := vector(ri):\n # Count h ow many times each non-zero label is used.\n for i from 1 to rs do\n for j from 1 to cs do\n if (R[i,j] <> 0) then\n \+ l := R[i,j]:\n C[l] := C[l] + 1:\n fi:\n od :\n od:\n # Old region label i goes to new region label k.\n # O bserve that the default zeroes in C are used to \n # indicate an old label i which is not in use.\n k := 1:\n for i from 1 to ri do\n if (C[i] > 0) then\n N[i] := k:\n k := k + 1:\n \+ fi:\n od:\n # Map in the new region labels by modifying the regio n label image R.\n for i from 1 to rs do\n for j from 1 to cs d o\n if (R[i,j] <> 0) then\n R[i,j] := N[R[i,j]]:\n \+ fi:\n od:\n od:\n eval(R):\nend:" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 29 "Connected Components Procdure" }}{PARA 0 "" 0 " " {TEXT -1 416 "This procedure uses the procedure defined above. It f irst places initial labels in R. Next it propogates labels until all c onnected pixels have the same label. Finally, it resequences the label s so they begin at 1 and run up to k. There is a maxIter term to prev ent endless looping. It is set to be the row dimension of the region l abel image: this should be enought to guarantee convergence of the reg ion labels. " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 250 "connectedCo mponents := proc(L)\n local R, flag, maxIter, i:\n maxIter := rowd im(L):\n i := 0:\n R := initLabels(L):\n flag := true:\n while (flag and (i < maxIter)) do\n flag := upLabels(L, R):\n i : = i + 1:\n od:\n seqLabels(R):\nend:" }}}}}{SECT 1 {PARA 3 "" 0 " " {TEXT -1 38 "Build One List of Edge Support Regions" }}{PARA 0 "" 0 "" {TEXT -1 469 "There is a one-to-one correspondence between the fina l straight line segments produced by the Burns algorithm and the regio ns in the two edge support region images. One caveat, regions with few er than a specified number of pixels are usually not considered furthe r. Thus, in this step a master list of regions is created. What is st ored 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. " }} {SECT 1 {PARA 4 "" 0 "" {TEXT -1 30 "Find the max value in an image" } }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 199 "imMaxVal := proc(A)\n lo cal i, j, max:\n max := 0:\n for i from 1 to rowdim(A) do\n f or j from 1 to coldim(A) do\n if (A[i,j] > max) then max := A[ i,j] fi:\n od:\n od:\n max:\nend:" }}}}{SECT 1 {PARA 4 "" 0 " " {TEXT -1 20 "Remove Small Regions" }}{PARA 0 "" 0 "" {TEXT -1 512 "G o 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 si ze are mapped to the background - zero - region. Also, include an offs et 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 simplif ies later processing." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1109 "r emSmallRegions := proc(R,thresh,offset)\n local i, j, max, C, RV, k: \n # print(\"Regions Entering Small Regions\", R):\n max := imMaxV al(R): # The largest label used.\n C := vector(max,0): \+ # Count, how many pixels use a label\n RV := vector(max,0): # New labels for regions above thresh in size\n # Count the number of pixel in each region\n for i from 1 to rowdim(R) do\n for j fr om 1 to coldim(R) do\n if (R[i,j] <> 0) then C[R[i,j]] := C[R[ i,j]] + 1: fi:\n od:\n od:\n # Initialize the new count, k, a t the offset plus 1 and begin entering new labels for large regions.\n # It is important that all entries for small regions already contai n the label \"0\".\n k := offset + 1:\n for i from 1 to max do\n \+ if C[i] > thresh then \n RV[i] := k:\n k := k + 1: \n fi:\n od:\n # Go back through giving new labels to all non zero labelled pixels.\n for i from 1 to rowdim(R) do\n for j f rom 1 to coldim(R) do\n if (R[i,j] <> 0) then R[i,j] := RV[R[i ,j]]: fi:\n od:\n od:\n # print(\"Regions Leaving Small Regio ns\", R):\n R:\nend:\n" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 21 "Mak e the Regions List" }}{PARA 0 "" 0 "" {TEXT -1 319 "Note that since th is procedure modifies the region label planes, it also returns them. S pecifically, this function returns a list with three elements. The fir st is the list of regions, i.e. the region index, source image, extent s and size. Next is the first region label plane followed by the secon d region label plane." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1746 "m akeRegList := proc(Regs,minPixels)\n local nRegs, R1, R2, R1n, R2n, \+ RV, RR, RS, rs, cs, i, j:\n rs := rowdim(Regs[1]): cs := coldim(Reg s[1]):\n R1 := remSmallRegions(Regs[1],minPixels,0):\n R1n := imM axVal(R1):\n R2 := remSmallRegions(Regs[2],minPixels,R1n):\n R2n \+ := imMaxVal(R2):\n RV := vector(R2n): # The extents of the regi on\n RR := vector(R2n): # The label plane for region, 1 or 2\n \+ RS := vector(R2n,0): # The size of each region\n for i from 1 t o R2n do\n RV[i] := [rs, cs, 1, 1]:\n od:\n # Create the labe l plane index, used later as a back pointer, for each region.\n for \+ i from 1 to R1n do\n RR[i] := 1:\n od:\n for i from R1n + 1 t o R2n do\n RR[i] := 2:\n od:\n # Computer to size of the regi ons\n for i from 1 to rs do\n for j from 1 to cs do\n i f (R1[i,j] <> 0) then RS[R1[i,j]] := RS[R1[i,j]] + 1: fi:\n if (R2[i,j] <> 0) then RS[R2[i,j]] := RS[R2[i,j]] + 1: fi:\n od:\n \+ od:\n # Go through all pixels in parallel, updating extents.\n f or i from 1 to rs do\n for j from 1 to cs do\n if R1[i,j] > 0 then\n if RV[R1[i,j]][1] > i then RV[R1[i,j]][1] := i: fi:\n if RV[R1[i,j]][2] > j then RV[R1[i,j]][2] := j: fi: \n if RV[R1[i,j]][3] < i then RV[R1[i,j]][3] := i: fi:\n \+ if RV[R1[i,j]][4] < j then RV[R1[i,j]][4] := j: fi:\n \+ fi:\n if R2[i,j] > 0 then\n if RV[R2[i,j]][1] > i \+ then RV[R2[i,j]][1] := i: fi:\n if RV[R2[i,j]][2] > j then \+ RV[R2[i,j]][2] := j: fi:\n if RV[R2[i,j]][3] < i then RV[R2 [i,j]][3] := i: fi:\n if RV[R2[i,j]][4] < j then RV[R2[i,j] ][4] := j: fi:\n fi:\n od:\n od:\n [[eval(RR), eval(R V), eval(RS)], eval(R1), eval(R2)]:\nend:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 25 "Find Line \+ for Each Region" }}{PARA 0 "" 0 "" {TEXT -1 273 "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 pl aced at the average intensity for the edge support region. Then clip t he line to the extents box of the region. " }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 30 "Plane Fitting, the Mathematics" }}{PARA 0 "" 0 "" {TEXT -1 71 "The plane fitting is done by using the explicit equation for th e plane." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "eq3 := z = alpha *x + beta*y + gamma:\neq3;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#/%\"zG,(*&%&alphaG\"\"\"%\"xGF(F(*&%%b etaGF(%\"yGF(F(%&gammaGF(" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 14 "The \+ values of " }{XPPEDIT 18 0 "alpha;" "6#%&alphaG" }{TEXT -1 1 "," } {XPPEDIT 18 0 "beta;" "6#%%betaG" }{TEXT -1 5 " and " }{XPPEDIT 18 0 " gamma;" "6#%&gammaG" }{TEXT -1 457 " 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. L et G be the gradient magnitude image used to weight the contribution o f 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. \nThe plane must minimize the following error funct ion:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 152 "unassign('i','j',' x','y','err','A','G'):\nerr := G[i,j]*(alpha*i+beta*j+gamma-A[i,j])^2: \neq4 := E = Sum(Sum(err,j=y[min]..y[max]),i=x[min]..x[max]):\neq4;" } }{PARA 11 "" 1 "" {XPPMATH 20 "6#/%\"EG-%$SumG6$-F&6$*&&%\"GG6$%\"iG% \"jG\"\"\"),**&%&alphaGF0F.F0F0*&%%betaGF0F/F0F0%&gammaGF0&%\"AGF-!\" \"\"\"#\"\"\"/F/;&%\"yG6#%$minG&F@6#%$maxG/F.;&%\"xGFA&FIFD" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 113 "Take the partial derivatives of t his error with respect to the three plane parameters and set them equa l to zero." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 134 "DM := matrix (3,1, [diff(rhs(eq4),alpha),diff(rhs(eq4),beta),diff(rhs(eq4),gamma)]) :\neq5 := eval(DM) = matrix(3,1, [0,0,0]):\neval(eq5);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/-%'matrixG6#7%7#-%$SumG6$-F*6$,$*(&%\"GG6$%\"iG% \"jG\"\"\",**&%&alphaGF5F3F5F5*&%%betaGF5F4F5F5%&gammaGF5&%\"AGF2!\"\" F5F3\"\"\"\"\"#/F4;&%\"yG6#%$minG&FD6#%$maxG/F3;&%\"xGFE&FMFH7#-F*6$-F *6$,$*(F0F?F6F?F4F?F@FAFJ7#-F*6$-F*6$,$*&F0F?F6F?F@FAFJ-F%6#7%7#\"\"!F jnFjn" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 105 "Expand the resulting eq uation, taking the summations inside the matrix and moving the paramet ers outside." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1465 "DM2 := ma trix(3,1):\nDM2[1,1] := alpha*Sum(Sum(G[i,j]*i*i,j=y[min]..y[max]),i=x [min]..x[max])\n + beta*Sum(Sum(G[i,j]*j*i,j=y[min]..y[max]), i=x[min]..x[max])\n + gamma*Sum(Sum(G[i,j]*i,j=y[min]..y[max] ),i=x[min]..x[max])\n - Sum(Sum(G[i,j]*A[i,j]*i,j=y[min]..y[m ax]),i=x[min]..x[max]):\nDM2[2,1] := alpha*Sum(Sum(G[i,j]*i*j,j=y[min] ..y[max]),i=x[min]..x[max])\n + beta*Sum(Sum(G[i,j]*j*j,j=y[m in]..y[max]),i=x[min]..x[max])\n + gamma*Sum(Sum(G[i,j]*j,j=y [min]..y[max]),i=x[min]..x[max])\n - Sum(Sum(G[i,j]*A[i,j]*j, j=y[min]..y[max]),i=x[min]..x[max]):\nDM2[3,1] := alpha*Sum(Sum(G[i,j] *i,j=y[min]..y[max]),i=x[min]..x[max])\n + beta*Sum(Sum(G[i,j ]*j,j=y[min]..y[max]),i=x[min]..x[max])\n + gamma*Sum(Sum(G[i ,j],j=y[min]..y[max]),i=x[min]..x[max])\n - Sum(Sum(G[i,j]*A[ i,j],j=y[min]..y[max]),i=x[min]..x[max]):\neval(DM2);\nMMij := matrix( 3,3, [G[i,j]*i*i, G[i,j]*i*j, G[i,j]*i,\n G[i,j]*i *j, G[i,j]*j*j, G[i,j]*j,\n G[i,j]*i, G[i,j]*j, \+ G[i,j]]):\nVVij := matrix(3,1, [G[i,j]*A[i,j]*i, G[i,j]*A[i,j]*j, G[ i,j]*A[i,j]]):\nMM := Sum(Sum(eval(MMij),j=y[min]..y[max]),i=x[min]. .x[max]):\nVV := Sum(Sum(eval(VVij),j=y[min]..y[max]),i=x[min]..x[ma x]):\neval(MM) * matrix(3,1, [alpha,beta,gamma]) = eval(VV);\nMMs := m atrix(3,3, [a, b, c, b, d, e, c, e, h]):\nVVs := matrix(3,1, [m, n, o] ):\nXXs := matrix(3,1, [alpha,beta,gamma]):\neval(MMs) * eval(XXs) = e val(VVs);\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG6#7%7#,**&%&a lphaG\"\"\"-%$SumG6$-F-6$*&&%\"GG6$%\"iG%\"jGF+)F5\"\"#\"\"\"/F6;&%\"y G6#%$minG&F=6#%$maxG/F5;&%\"xGF>&FFFAF+F+*&%%betaGF+-F-6$-F-6$*(F2F9F6 F+F5F+F:FCF+F+*&%&gammaGF+-F-6$-F-6$*&F2F9F5F9F:FCF+F+-F-6$-F-6$*(F2F9 &%\"AGF4F+F5F9F:FC!\"\"7#,**&F*F9FJF9F+*&FIF9-F-6$-F-6$*&F2F9)F6F8F9F: FCF+F+*&FPF9-F-6$-F-6$*&F2F9F6F9F:FCF+F+-F-6$-F-6$*(F2F9FenF9F6F9F:FCF gn7#,**&F*F9FQF9F+*&FIF9FcoF9F+*&FPF9-F-6$-F-6$F2F:FCF+F+-F-6$-F-6$*&F 2F9FenF9F:FCFgn" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/*&-%$SumG6$-F&6$-% 'matrixG6#7%7%*&&%\"GG6$%\"iG%\"jG\"\"\")F3\"\"#\"\"\"*(F0F8F4F5F3F5*& F0F8F3F87%F9*&F0F8)F4F7F8*&F0F8F4F87%F:F>F0/F4;&%\"yG6#%$minG&FC6#%$ma xG/F3;&%\"xGFD&FLFGF5-F+6#7%7#%&alphaG7#%%betaG7#%&gammaGF5-F&6$-F&6$- F+6#7%7#*(F0F8&%\"AGF2F5F3F87#*(F0F8FjnF8F4F87#*&F0F8FjnF8F@FI" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#/*&-%'matrixG6#7%7%%\"aG%\"bG%\"cG7%F+ %\"dG%\"eG7%F,F/%\"hG\"\"\"-F&6#7%7#%&alphaG7#%%betaG7#%&gammaGF2-F&6# 7%7#%\"mG7#%\"nG7#%\"oG" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 85 "To sol ve for the plane parameters, invert the matrix. This can be done symb olically." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 172 "MMsi := inver se(MMs):\nMinv = eval(MMsi);\nMMsin := scalarmul(MMsi,(-a*d*h+a*e^2+b^ 2*h-2*b*c*e+c^2*d)):\nMMsid := -a*d*h+a*e^2+b^2*h-2*b*c*e+c^2*d:\nMinv = eval(MMsin) / MMsid;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%%MinvG-%' matrixG6#7%7%*&,&*&%\"dG\"\"\"%\"hGF.F.*$)%\"eG\"\"#\"\"\"!\"\"F4,,*(% \"aGF.F-F4F/F4F.*&F8F4F1F4F5*&)%\"bGF3F4F/F4F5*(F F3F4F-F4F5!\"\",$*&,&*&FF4F2F4F5F4F6FAF5*&,&*&FF4F-F4F5F4F6FA7%FB*&,&*&F8F4F/F4F.*$F@F4F5F4F6FA,$*&,&*&F8F4F2F4F.* &FF4F5F4F6FAF57%FGFP*&,&*&F8F4F-F4F.*$F;F4F5F4F6FA" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#/%%MinvG*&-%'matrixG6#7%7%*&*&,,*(%\"aG\"\"\"%\" dGF0%\"hGF0!\"\"*&F/\"\"\")%\"eG\"\"#F5F0*&)%\"bGF8F5F2F5F0*(F;F0%\"cG F0F7F0!\"#*&)F=F8F5F1F5F0F0,&*&F1F5F2F5F0*$F6F5F3F0F5,,F.F0F4F3F9F3F " 0 " " {MPLTEXT 1 0 0 "" }}}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 22 "Solving \+ for the plane " }}{SECT 1 {PARA 5 "" 0 "" {TEXT -1 29 "Plane parameter s from M and B" }}{PARA 0 "" 0 "" {TEXT -1 95 "Use the formula above t o find the equation for the plane, i.e. solve for alpha, beta and gamm a." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 724 "planeFromMandB := pro c(M,B)\n local a, b, c, d, e, h, m, n, o, K, Mi, X, holdDigits:\n \+ holdDigits := Digits:\n a := M[1,1]: \n b := M[1,2]:\n c := M[ 1,3]:\n d := M[2,2]:\n e := M[2,3]:\n h := M[3,3]:\n m := B[1, 1]:\n n := B[2,1]:\n o := B[3,1]:\n Mi := matrix(3,3):\n K := -a*d*h+a*e^2+b^2*h-2*b*c*e+c^2*d:\n # print(\"K in planeFromMandB \+ \", K):\n Mi[1,1] := (-d*h+e^2)/K:\n Mi[1,2] := (b*h-c*e)/K:\n M i[1,3] := (-b*e+c*d)/K:\n Mi[2,1] := Mi[1,2]:\n Mi[2,2] := (-a*h+c ^2)/K:\n Mi[2,3] := (a*e-b*c)/K:\n Mi[3,1] := Mi[1,3]:\n Mi[3,2] := Mi[2,3]:\n Mi[3,3] := (-a*d+b^2)/K:\n Digits := holdDigits:\n \+ X := evalm(Mi &* B):\n # print(\"Test Plane Solution: \", linsolve (M,B), \"X \", X):\n #X:\n X:\nend:" }{TEXT -1 0 "" }{MPLTEXT 1 0 1 "\n" }}}}{SECT 0 {PARA 5 "" 0 "" {TEXT -1 27 "M and B matrix from re gions" }}{PARA 0 "" 0 "" {TEXT -1 85 "Compute the M matrix and the B c olumn vector for the region indicated by the index ri" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1198 "regMandB := proc(RL, ri, Rs, Gs, A)\n \+ local M, B, lab, ext, L, G, i, j:\n M := matrix(3,3, 0):\n B := matr ix(3,1, 0):\n lab := RL[1][ri]: # the index of the label and gradien t plane\n ext := RL[2][ri]: # ext is the extends of the region\n # \+ print(\"Region \", ri, \"Label \", lab, \"Extents \", ext, \"Size\", R L[3][ri]): \n for i from ext[1] to ext[3] do\n for j from ext[2 ] to ext[4] do\n L := Rs[lab]: # the label plane\n G := Gs[1]: # the gradient plane\n # This had been a test for gr eater than zero - an error.\n if (L[i,j] = ri) then\n \+ M[1,1] := M[1,1] + (G[i,j] * i * i):\n M[1,2] := M[1,2] + ( G[i,j] * i * j):\n M[1,3] := M[1,3] + (G[i,j] * i):\n \+ M[2,1] := M[1,2]:\n M[2,2] := M[2,2] + (G[i,j] * j * j): \n M[2,3] := M[2,3] + (G[i,j] * j):\n M[3,1] := M[ 1,3]:\n M[3,2] := M[2,3]:\n M[3,3] := M[3,3] + G[i ,j]:\n B[1,1] := B[1,1] + (G[i,j] * A[i,j] * i):\n \+ B[2,1] := B[2,1] + (G[i,j] * A[i,j] * j):\n B[3,1] := B[3,1 ] + (G[i,j] * A[i,j]):\n fi:\n od: # End of colum within \+ extents\n od: # End of row within extents\n [eval(M), eval(B)]: \nend:" }}}}{SECT 0 {PARA 5 "" 0 "" {TEXT -1 25 "Line from Plane & Ext ents" }}{PARA 0 "" 0 "" {TEXT -1 45 "Represent the line in with the im plicit form " }{XPPEDIT 18 0 "alpha*x+beta*y = rho;" "6#/,&*&%&alphaG \"\"\"%\"xGF'F'*&%%betaGF'%\"yGF'F'%$rhoG" }{TEXT -1 95 " and then com pute the for points where the infinite line crosses boundaries of the \+ extents box." }}{SECT 0 {PARA 20 "" 0 "" {TEXT -1 59 "Equations for x \+ in terms of y and z, y in terms of x and z." }}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 35 "x = solve(eq3,x);\ny = solve(eq3,y);" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#/%\"xG,$*&,(%\"zG!\"\"*&%%betaG\"\"\"%\"yGF,F,%& gammaGF,\"\"\"%&alphaG!\"\"F)" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%\"y G,$*&,(%\"zG!\"\"*&%&alphaG\"\"\"%\"xGF,F,%&gammaGF,\"\"\"%%betaG!\"\" F)" }}}}{SECT 1 {PARA 20 "" 0 "" {TEXT -1 25 "Procedures for xofy, etc ." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 183 "xofy := proc(y,X,avg) \n -(y*X[2,1]+X[3,1]-avg)/X[1,1]\nend:\nyofx := proc(x,X,avg)\n -( x*X[1,1]+X[3,1]-avg)/X[2,1]\nend:\nzofxy := proc(x,y,X)\n X[1,1] * x + X[2,1] * y + X[3,1]:\nend:\n" }}}}{SECT 1 {PARA 20 "" 0 "" {TEXT -1 34 "3D Coords for Corners of Fit Plane" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 279 "planeCorners := proc(X,ext)\n local xmin, ymin, xm ax, ymax:\n xmin := ext[1]: ymin := ext[2]: xmax := ext[3]: ymax := ext[4]:\n [[xmin,ymin,zofxy(xmin,ymin,X)],\n [xmin,ymax,zofxy(xm in,ymax,X)],\n [xmax,ymax,zofxy(xmax,ymax,X)],\n [xmax,ymin,zofx y(xmax,ymin,X)]]:\nend:\n" }}}}{SECT 1 {PARA 20 "" 0 "" {TEXT -1 23 "L ine Segment from Plane" }}{PARA 0 "" 0 "" {TEXT -1 227 "Check first fo r vertical or horizontal special case. Otherwise, intersect infinite l ine with extents box and return the endpoints. Note, this code does no t check for possible ties and thus intersection with more than two sid es." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1246 "lineFromPlane := pr oc(X,ext,avg)\n local Line, xt, yt, xmin, ymin, xmax, ymax:\n xmi n := ext[1]: ymin := ext[2]: xmax := ext[3]: ymax := ext[4]:\n if ( X[1,1] = 0) then # Vertical line\n Line := [[xmin,yofx(xmin,X,av g)], [xmax,yofx(xmax,X,avg)]]:\n elif (X[2,1] = 0) then # Horizontal line\n Line := [[xofy(ymin,X,avg),ymin],[xofy(ymax,X,avg)]]:\n \+ else\n Line := []:\n # Test if line crosses left boundary at ymin\n xt := xofy(ymin,X,avg):\n if ((xt > xmin) and (xt < \+ xmax)) then\n Line := [op(Line), [xt,ymin]]:\n fi:\n \+ # Test if line crosses right boundary at ymax\n xt := xofy(ymax, X,avg):\n if ((xt > xmin) and (xt < xmax)) then\n Line := [op(Line), [xt,ymax]]:\n fi:\n # Test if line crosses top b oundary at xmin\n yt := yofx(xmin,X,avg):\n if ((yt > ymin) \+ and (yt < ymax)) then\n Line := [op(Line), [xmin,yt]]:\n \+ fi:\n # Test if line crosses botom boundary at xmax\n yt := \+ yofx(xmax,X,avg):\n if ((yt > ymin) and (yt < ymax)) then\n \+ Line := [op(Line), [xmax,yt]]:\n fi:\n fi:\n if (nops(Line ) <> 2) then\n printf(\"Line to Extents Intersection Failed for L ine\\n\"):\n Line := [[0.0, 0.0], [0.0, 0.0]]:\n fi:\n Line: \nend:" }}}}}{SECT 1 {PARA 5 "" 0 "" {TEXT -1 18 "Voting for support" }}{PARA 0 "" 0 "" {TEXT -1 526 "Each line segment, or equivalently, ea ch edge support region, can be said to be supported by some pixels wit hin it. To see what this means, realize that every pixel in the origin al 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 w hich it belongs. Thus, the support is the ratio of the number of pixe ls 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. \+ \n" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 697 "lineSupport := proc(S ize,len,Rs)\n local i, j, rs, cs, r1, r2, s1, s2, Votes:\n Votes : = vector(len,0.0):\n rs := rowdim(Rs[1]):\n cs := coldim(Rs[1]):\n for i from 1 to rs do\n for j from 1 to cs do\n r1 := Rs[1][i,j]: r2 := Rs[2][i,j]:\n if (r1 > 0) or (r2 > 0) t hen\n if (r1 > 0) then s1 := Size[r1]: else s1 := 0: fi:\n if (r2 > 0) then s2 := Size[r2]: else s2 := 0: fi:\n \+ if (s1 > s2) then \n Votes[r1] := Votes[r1] + 1 : \n else\n Votes[r2] := Votes[r2] + 1:\n \+ fi:\n fi:\n od:\n od:\n for i from 1 to l en do\n Votes[i] := Votes[i] / Size[i]:\n od:\n Votes:\nend: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 20 "Map Regions to Line s" }}{PARA 0 "" 0 "" {TEXT -1 370 "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 eleme nt 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." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 139 "line3D := proc(Line,avg)\n # pri nt(\"Entering Line 3D\", Line, avg):\n [[Line[1][1], Line[1][2], avg ], [Line[2][1], Line[2][2], avg]]:\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 760 "regToLines := proc(RL, Rs, Gs, A)\n local len, ri, i, j, avg, Supp, M, B, MB, PL, PC, LN, RES:\n # Setup to solve M*X= B, where X is the a, b and c for the plane\n len := rhs(op(2,eval(RL [2]))):\n # Accumulate results in vector RES\n RES := vector(len): \n # Compute support for regions\n Supp := lineSupport(RL[3],len,R s):\n # For each of the regions in the Label Planes List\n for ri \+ from 1 to len do\n MB := regMandB(RL, ri, Rs, Gs, A):\n M : = MB[1]: B := MB[2]: avg := B[3,1] / M[3,3]:\n PL := planeFromM andB(M,B):\n PC := planeCorners(PL,RL[2][ri]):\n LN := lineF romPlane(PL,RL[2][ri],avg):\n RES[ri] := [ri, LN, line3D(LN,avg), RL[2][ri], PC, RL[3][ri], Supp[ri]]:\n od: #End of iteration thro ugh regions\n RES:\nend:" }}}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 22 " The Complete Algorithm" }}{PARA 0 "" 0 "" {TEXT -1 115 "Here are the p ieces put together. The input is a single image array. The output is a sequence of 2D line segments. " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 758 "burns := proc(A,magThresh,minRegSize)\n local Sh, Sv, Gs, L s, Rs,Rs1, Rs2, RGL:\n Sh := matrix(3,3, [-1,-2,-1,0,0,0,1,2,1]):\n \+ Sv := transpose(Sh):\n print(\"Computing Gradient Images\"):\n G s := GradImages(A,Sh,Sv):\n #print(A):\n #print(Gs[1]):\n print( \"Computing Bucket Labels\"):\n Ls := bucketLabel(Gs,magThresh):\n \+ print(\"Computing Connected Components on Labels 1\"):\n Rs1 := con nectedComponents(Ls[1]):\n print(\"Computing Connected Components on Labels 2\"):\n Rs2 := connectedComponents(Ls[2]):\n Rs := [eval(R s1), eval(Rs2)]:\n print(\"Computing Edge Support Regions List\"):\n RGL := makeRegList(Rs,minRegSize): \n print(\"Computing Line Segm ents for Regions\");\n Ls := regToLines(RGL[1], [RGL[2], RGL[3]], Gs ,A):\n eval(Ls):\nend:\n " }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 18 "Test the Algorithm" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 188 "A := readImage(`/s/bach/b/class/cs510/public_html/Labs/Lab2/BurnsTest.p gm`):\n# A := testImage():\n# A := readImage(`CS 510, Spring 99:CS510 \+ Current:Assignments:CagedSphereGrayscale.pgm`):\n\n" }}{PARA 11 "" 1 " " {XPPMATH 20 "6'Q/In~Read~Image~6\"Q#P56\"\"#I\"#D\"$b#" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "LS := burns(A,5.0,5):" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#Q:Computing~Gradient~Images6\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#Q8Computing~Bucket~Labels6\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#QKComputing~Connected~Components~on~Labels~16\"" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#QKComputing~Connected~Components~on~La bels~26\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#QDComputing~Edge~Support~ Regions~List6\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#QDComputing~Line~Se gments~for~Regions6\"" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 29 "Write \+ Line Segments to a File" }}{PARA 0 "" 0 "" {TEXT -1 136 "After running the Burns Algorithm, the following procedure can read the sequence of line segments and write the endpoints to a text file" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 328 "writeLines := proc(LS,linesFile)\n loc al i, len, fd:\n len := rhs(op(2,eval(LS))):\n fd := fopen(linesF ile,WRITE); \n for i from 1 to len do\n if (LS[i][7] > 1/2) the n\n 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]):\n fi:\n od: \n fclose(fd):\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "# writeLines(LS,`CS 510, Spring 99:CS510 Current:Assignments:BurnsTest. dat`);" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 18 "Plot Line Segments" } }{PARA 0 "" 0 "" {TEXT -1 55 "Run the Burns Algorithm and draw the res ulting segments" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 413 "drawLine s2D := proc(LS,A)\n local objs, len, i, rs, cs:\n rs := rowdim(A): cs := coldim(A):\n objs := \{\}:\n len := rhs(op(2,eval(LS))): \n for i from 1 to len do\n if (LS[i][7] > 1/2) then\n \+ objs := \{op(objs),CURVES(LS[i][2],THICKNESS(3),COLOR(HUE,0.4))\}:\n \+ fi:\n od:\n display(objs,axes=FRAME,view=[0..rs,0..cs],scaling =CONSTRAINED,font=[HELVETICA,BOLD,14],labels=[\"Rows\",\"Columns\"]); \nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 632 "drawLines3D := p roc(LS,A)\n local objs, zVals, len, i, j, rs, cs:\n rs := rowdim(A ): cs := coldim(A):\n objs := \{\}:\n zVals := []:\n len := rhs(op(2,eval(LS))):\n for i from 1 to rs do\n zVals := [op(zV als),[]]:\n for j from 1 to cs do\n zVals[i] := [op(zVal s[i]),A[i,j]]:\n od:\n od:\n objs := \{op(objs),GRID(1..rs,1. .cs,zVals,COLOR(ZGREYSCALE))\}:\n for i from 1 to len do\n if ( LS[i][7] > 1/2) then\n objs := \{op(objs),CURVES(LS[i][3],THIC KNESS(10),COLOR(HUE,0.4))\}:\n fi:\n od:\n #print(objs):\n \+ display(objs,axes=FRAME,view=[0..rs,0..cs,0..255],font=[HELVETICA,BOLD ,14]);\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "drawLines2D (LS,A);\n\n" }}{PARA 13 "" 1 "" {INLPLOT "6*-%'CURVESG6%7$7$$\"+\"G.>v (!\"*\"\"#7$$\"+>\"\\DF\"!\")\"#F-%*THICKNESSG6#\"\"$-%&COLORG6$%$HUEG $\"\"%!\"\"-F$6%7$7$\"#8$\"+%=,\"*f#F/7$\"#C$\"+8?2*f#F/F1F5-F$6%7$7$ \"#9$\"+PY4s@F/7$FD$\"+c\\F " 0 "" {MPLTEXT 1 0 18 "drawLines3D(LS,A);" }}{PARA 13 "" 1 "" {INLPLOT "6)-%'CURVESG6 %7$7%\"#8$\"+%=,\"*f#!\")$\"+R#yd8%F+7%\"#C$\"+8?2*f#F+F,-%*THICKNESSG 6#\"#5-%&COLORG6$%$HUEG$\"\"%!\"\"-F$6%7$7%\"#9$\"+PY4s@F+$\"+`X8<5!\" (7%F/$\"+c\\F\"FgnFgnFgnFhnFhnFhnFhnFhn\"$;\"\"#')FjnFjnFjn\"#cF[oFUFUFUFUFUFU FUFUFUFUFUFUFUFU7@F]oF]oFgnFgnFgnFhnFhnFhnFhnFhnF^oF^o\"$:\"FboFboF_oF `oF`oF`oF`oF[oFUFUFUFUFUFUFUFUFU7@\"$?\"F]oF]oFgnFgnFgnFhnFhnFhnFhnFhn F^oF^oFboFboFbo\"$9\"Feo\"$8\"Ffo\"#&)\"#bFhoFhoFho\"#FFUFUFUFU7@FdoFd oF]oF]oFgnFgnFgnFhnFhnFhnFhnF^oF^oF^oFboFboFboFeoFeoFfoFfoFfo\"$7\"\"# (*Fgo\"#UFUFUFUFU7@FdoFdoFdoF]oF]oFgnFgnFhnFhnFhnFhnFhnF^oF^oF^oFboFbo FeoFeoFeoFfoF\\pFgo\"#\")\"#$)\"#SFUFUFUFU7@\"$@\"FdoFdoF]oF]oF]oFgnFg nFhnFhnFhnFhnFhnF^oF^oF^oFboFboFeo\"#)*F_oF_p\"#%)F`pF`pFapFUFUFUFU7@F cpFcpFdoFdoF]oF]oFgnFgnFgnFhnFhnFhnFhnFhnF^oF^oFboFboFdp\"##)FepFepFep F`pF`pFapFUFUFUFU7@FcpFcpFdoFdoFdoF]oF]oFgnFgnFgnFhnFhnFhnFhnF^oF^o\"# **FinFgpFepFepFepFepFepF`pFapFUFUFUFU7@\"$A\"FcpFcpFdoFdoF]oF]oF]oFgnF gnFhnFhnFhnFhn\"$+\"FinFgpFgoFgoFgoFepFepFepFepFepFapFUFUFUFU7@F[qFcpF cpFcpFdoFdoF]oF]oFgnFgnFgnFhnF\\qFYF`pFgoFgoFgoFgoFgoFgoFepFepFepFepFa pFUFUFUFU7@F[qF[qFcpFcpFdoFdoF]oF]oF]oFgn\"$,\"FYF`pF_oF_oF_oFgoFgoFgo FgoFgoFepFepFepFepFapFUFUFUFU7@F[qF[qFcpFcpFcpFdoFdoF]oF]o\"$-\"FepF_o F_oF_oF_oF_oF_oFgoFgoFgoFgoFgoFepFepFepFapFUFUFUFU7@F[qF[qF[qFcpFcpFdo FdoFaq\"#*)FepFinFinF_oF_oF_oF_oF_oF_oFgoFgoFgoFgoFgoFepFepFapFUFUFUFU 7@F[qF[qF[qFcpFcpFaq\"#!*FepFinFinFinFinFinF_oF_oF_oF_oF_oF_oFgoFgoFgo FgoFepFepFapFUFUFUFU7@F[qF[qF[q\"$.\"FeqFgoFYFYFinFinFinFinFinFinF_oF_ oF_oF_oF_oFgoFgoFgoFgoFgoFepFapFUFUFUFU7@F[q\"$/\"\"#\"*FgoFYFYFYFYFYF inFinFinFinFinFinF_oF_oF_oF_oF_oFgoFgoFgoFgo\"#j\"#>FUFUFUFU-F76#%+ZGR EYSCALEG-F$6%7$7%$\"+\"G.>v(FJ\"\"#$\"+)\\f'RfF+7%$\"+>\"\\DF\"F+FioFg rF2F6-%*AXESSTYLEG6#%&FRAMEG-%%FONTG6%%*HELVETICAG%%BOLDGFA-%%VIEWG6%; FU$FPFU;FU$FRFU;FU$\"$b#FU" 3 668 550 550 2 0 1 0 2 1 0 3 2 1.000000 45.000000 45.000000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 18835 -16399 0 0 0 0 0 0 }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}}{MARK "12 5 1 " 0 }{VIEWOPTS 1 1 0 1 1 1803 }