{VERSION 3 0 "APPLE_PPC_MAC" "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 "" 11 12 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }1 0 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 1 {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" }}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 19 " Create a 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(widt h,height):\n noise := rand(0..10):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 232 " for i from 1 to 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 from 4 to 8 do\n A[i,j] := A[i,j] + 100. 0:\n od:\n od:\n eval(A):\nend:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "A := te stImage();" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"AG-%'matrixG6#7.7.$ \"%!4\"!\"\"$\"%S5F,$\"%55F,F-F-F/$\"%I5F,$\"%!3\"F,F*F*$\"%?5F,F57.F* F-$\"%]5F,$\"%+6F,$\"%g5F,F:F1F3$\"%+5F,F>F$\"% q5F,F8F<7.FF87.F3F>FF3F>F AF/F87.F3F>F*F:F-FAFAF1F/F " 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.F 1F7" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 50 "The procedure to find Gr adient Mag and Orientation" }}{PARA 0 "" 0 "" {TEXT -1 368 " The proc edure below convolves IM with the two Sobel Masks, Sh and Sv. An inter mediate 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 gradi ent magnitude is used. Often a city block approximation is used to sav e time. " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 688 "GradImages := p roc(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, r eturns the gradient images.\";\n rs := rowdim(IM); cs := coldim(IM); \n Pf := convert(Pi,float):\n Gm := matrix(rs,cs,0):\n Go := mat rix(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 fro m 2 to (cs-1) do\n dx := Ghf(r,c):\n dy := Gvf(r,c):\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 Bucket Labels" }} {PARA 0 "" 0 "" {TEXT -1 145 "Quantize the orientation into eight buck ets. 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 "The bucket labe ls 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 typical ly 8. They represent a quantization of the orientations. There is an \+ offset " }{XPPEDIT 18 0 "t;" "6#%\"tG" }{TEXT -1 59 " that 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 m agnitude is too small reliably yield a gradient orientation 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 228 "bucket := proc(o,g,k,t,m)\n local b;\n if (g < m ) then \n b := 0:\n else\n b := floor(((o+t)/(2*convert(Pi ,float)))*k)+1:\n fi:\n # Special Case when o is 2PI instead of ze ro\n if (b > k) then b := 1 fi:\n b:\nend:" }}}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 470 "bucketLabel := proc(gIms)\n local Go, Gm, L 1, L2, t1, t2, rs, cs, i, j:\n # DEBUG;\n Gm := gIms[1]: Go := \+ gIms[2]:\n rs := rowdim(Go): cs := coldim(Go):\n L1 := matrix(rs,c s):\n t1 := 0.0:\n L2 := matrix(rs,cs):\n t2 := (convert(Pi,floa t)*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,1.0):\n L2[i ,j] := bucket(Go[i,j],Gm[i,j],8,t2,1.0):\n od:\n od:\n [eval( L1), eval(L2)]:\nend:" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 65 "Findin g the Connected Components in the Bucket Orientation Images" }}{PARA 0 "" 0 "" {TEXT -1 109 "In this step, each 4 connected edge support re gion is given a unique integer 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 da ta propogating integer labels to adjacent pixels in each region. This \+ algorithm can be slow under certain pathological configurations. Howev er, these involve regions that form spirals, etc. These are nearly imp ossible configurations given the way that the edge support regions dep end upon the underlying gradient orientation. " }}{PARA 0 "" 0 "" {TEXT -1 514 "The first step builds two data structures. The first is \+ a region label image where each pixel with a non-zero gradient orienta tion direction gets a unique integer label. The second is a vector of labels used to propogate the labels. The indices to this vector repre sent the original unique region labels. The value stored is the curren t label to be used for this region. So, if regions 1 and 2 are collaps ed, position 2 will get the label 1. Initially, this vector just conta ins sequential integers starting at 1." }}{PARA 0 "" 0 "" {TEXT -1 354 "The second step is to take four sweeps through the image, going f rom left to right, top to bottom, right to left, and finally bottom to top, propogating region labels. If any labels change, set a flag indi cating that a change has been made. This algorithm is then used repeat edly to update labels until the result 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 se quential starting at 1. To create regions with sequential labels, a pa ss 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 regi on label vector, a similar vector is made in which those entries which are actually used in the region labeling are assigned sequential inte gers. 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 "in itLabels := proc(L)\n local i, j, rs, cs, R, ri:\n rs := rowdim(L) : cs := coldim(L):\n # Give each non zero pixel in label plane a uni que region label.\n R := matrix(rs,cs):\n ri := 1:\n for i from \+ 1 to rs do\n for j from 1 to 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 Neigh bors and Update" }}{PARA 0 "" 0 "" {TEXT -1 530 "This procedure is pas sed 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 check ed to see they are both non zero. Next, they are checked to see if the y 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 propog ated to the other. " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 598 "test Neighbors := 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) a nd (L[c,d] <> 0) then\n if (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 a nd columns" }}{PARA 0 "" 0 "" {TEXT -1 414 "Sweep from left to right, \+ top to bottom, right to left and bottom to top, propogating the smalle r of the two integer labels whenever a pair of adjacent pixels is foun d to be connected. Return false if no labels are propogated, otherwis e 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 fr om the left to the right collapsing 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 fr om 2 to (rs-1) do\n if testNeighbors(L,R,i,j,i+1,j) then flag \+ := true fi:\n od:\n od:\n # Go from the right to the left col lapsing labels\n for i from 1 to rs do\n for j from cs by -1 to 2 do\n if testNeighbors(L,R,i,j,i,j-1) then flag := true fi: \n od:\n od:\n # Go from the bottom to the topt collapsing la bels\n for j from 1 to cs do\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 Labels" }}{PARA 0 "" 0 "" {TEXT -1 287 " Af ter labels have been propogated, there will be holes in the numbering. In other words, there may be a connected region labeled \"1\" and a r egion labeled \"3\", but no region labelled \"2\". The following proce dure goes through a renumbers the regions sequentially starting with l abel \"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 how many times each non-zero label is us ed.\n for i from 1 to rs do\n for j from 1 to cs do\n i f (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 # Observe 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 regio n labels by modifying the region label image R.\n for i from 1 to rs do\n for j from 1 to cs do\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 Compone nts Procdure" }}{PARA 0 "" 0 "" {TEXT -1 416 "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. Fin ally, it resequences the labels so they begin at 1 and run up to k. T here 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 gu arantee convergence of the region labels. " }}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 250 "connectedComponents := proc(L)\n local R, flag, \+ maxIter, i:\n maxIter := rowdim(L):\n i := 0:\n R := initLabels( L):\n flag := true:\n while (flag and (i < maxIter)) do\n fla g := 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 Su pport Regions" }}{PARA 0 "" 0 "" {TEXT -1 469 "There is a one-to-one c orrespondence 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 r egions is created. What is stored in this list is the index of the re gion, its source image (either R1 or R2), the extents of the region an d the size of the region. " }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 30 "Fin d the max value in an image" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 199 "imMaxVal := proc(A)\n local i, j, max:\n max := 0:\n for i \+ from 1 to rowdim(A) do\n for j from 1 to coldim(A) do\n i f (A[i,j] > max) then max := A[i,j] fi:\n od:\n od:\n max:\ne nd:" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 20 "Remove Small Regions" }} {PARA 0 "" 0 "" {TEXT -1 512 "Go through the image once and determine \+ which regions are smallerr than a specified threshold. Use a Region La bel 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 be come apparent that this simplifies later processing." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1109 "remSmallRegions := proc(R,thresh,offset )\n local i, j, max, C, RV, k:\n # print(\"Regions Entering Small \+ Regions\", R):\n max := imMaxVal(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 from 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, at the offset plus 1 and begin enterin g new labels for large regions.\n # It is important that all entries for small regions already contain 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 thr ough giving new labels to all non zero labelled pixels.\n for i from 1 to rowdim(R) do\n for j from 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 # pr int(\"Regions Leaving Small Regions\", R):\n R:\nend:\n" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 21 "Make the Regions List" }}{PARA 0 "" 0 " " {TEXT -1 319 "Note that since this procedure modifies the region lab el 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." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1746 "makeRegList := proc(Regs,minPixel s)\n local nRegs, R1, R2, R1n, R2n, RV, RR, RS, rs, cs, i, j:\n rs := rowdim(Regs[1]): cs := coldim(Regs[1]):\n R1 := remSmallRegion s(Regs[1],minPixels,0):\n R1n := imMaxVal(R1):\n R2 := remSmallRe gions(Regs[2],minPixels,R1n):\n R2n := imMaxVal(R2):\n RV := vect or(R2n): # The extents of the region\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 to R2n do\n RV[i] := [rs, cs , 1, 1]:\n od:\n # Create the label plane index, used later as a b ack pointer, for each region.\n for i from 1 to R1n do\n RR[i] \+ := 1:\n od:\n for i from R1n + 1 to R2n do\n RR[i] := 2:\n \+ od:\n # Computer to size of the regions\n for i from 1 to rs do\n \+ for j from 1 to cs do\n if (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 pixel s in parallel, updating extents.\n for 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 th en RV[R1[i,j]][4] := j: fi:\n fi:\n if R2[i,j] > 0 the n\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(RV), 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 wi th a plane of constant intensity placed at the average intensity for t he edge support region. Then clip the line to the extents box of the r egion. " }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 30 "Plane Fitting, the Mat hematics" }}{PARA 0 "" 0 "" {TEXT -1 71 "The plane fitting is done by \+ using the explicit equation for the 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(*&%%betaGF(%\"yGF(F(%&gammaGF(" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 14 "The values of " }{XPPEDIT 18 0 "al pha;" "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 a nd 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 th at we will assume G is zero outside of each region. Thus, this equatio n is valid for the entire bounding box around the region. \nThe plane \+ must minimize the following error function:" }}}{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[mi n]..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 "Ta ke the partial derivatives of this error with respect to the three pla ne parameters and set them equal to zero." }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 134 "DM := matrix(3,1, [diff(rhs(eq4),alpha),diff(rhs(e q4),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*&%%betaG F5F4F5F5%&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#\"\"!FjnFjn" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 105 "Expand the resulting equation, taking the summations inside th e matrix and moving the parameters outside." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1465 "DM2 := matrix(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(Su m(G[i,j]*j*i,j=y[min]..y[max]),i=x[min]..x[max])\n + gamma*Su m(Sum(G[i,j]*i,j=y[min]..y[max]),i=x[min]..x[max])\n - Sum(Su m(G[i,j]*A[i,j]*i,j=y[min]..y[max]),i=x[min]..x[max]):\nDM2[2,1] := al pha*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[min]..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]):\nD M2[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(M Mij),j=y[min]..y[max]),i=x[min]..x[max]):\nVV := Sum(Sum(eval(VVij), j=y[min]..y[max]),i=x[min]..x[max]):\neval(MM) * matrix(3,1, [alpha,be ta,gamma]) = eval(VV);\nMMs := matrix(3,3, [a, b, c, b, d, e, c, e, h] ):\nVVs := matrix(3,1, [m, n, o]):\nXXs := matrix(3,1, [alpha,beta,gam ma]):\neval(MMs) * eval(XXs) = eval(VVs);\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG6#7%7#,**&%&alphaG\"\"\"-%$SumG6$-F-6$*&&%\"G G6$%\"iG%\"jGF+)F5\"\"#\"\"\"/F6;&%\"yG6#%$minG&F=6#%$maxG/F5;&%\"xGF> &FFFAF+F+*&%%betaGF+-F-6$-F-6$*(F2F9F6F+F5F+F:FCF+F+*&%&gammaGF+-F-6$- F-6$*&F2F9F5F9F:FCF+F+-F-6$-F-6$*(F2F9&%\"AGF4F+F5F9F:FC!\"\"7#,**&F*F 9FJF9F+*&FIF9-F-6$-F-6$*&F2F9)F6F8F9F:FCF+F+*&FPF9-F-6$-F-6$*&F2F9F6F9 F:FCF+F+-F-6$-F-6$*(F2F9FenF9F6F9F:FCFgn7#,**&F*F9FQF9F+*&FIF9FcoF9F+* &FPF9-F-6$-F-6$F2F:FCF+F+-F-6$-F-6$*&F2F9FenF9F:FCFgn" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/*&-%$SumG6$-F&6$-%'matrixG6#7%7%*&&%\"GG6$%\"iG% \"jG\"\"\")F3\"\"#\"\"\"*(F0F8F4F5F3F5*&F0F8F3F87%F9*&F0F8)F4F7F8*&F0F 8F4F87%F:F>F0/F4;&%\"yG6#%$minG&FC6#%$maxG/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 solve for the plane parameters, invert the matrix. This can be done symbolically." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 172 "MMsi := inverse(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\"\"#\"\"\"!\"\"F5,,*(%\"aGF/F.F5F0F5F6*&F9F5F2 F5F/*&)%\"bGF4F5F0F5F/*(F=F/%\"cGF/F3F/!\"#*&)F?F4F5F.F5F/!\"\"F6*&,&* &F=F5F0F5F/*&F?F5F3F5F6F5F7FC,$*&,&*&F=F5F3F5F/*&F?F5F.F5F6F5F7FCF67%F D,$*&,&*&F9F5F0F5F/*$FBF5F6F5F7FCF6,$*&,&*&F9F5F3F5F6*&F=F5F?F5F/F5F7F CF67%FHFS*&,&*&F9F5F.F5F6*$F " 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 1 {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 1132 "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 if (L[i,j] > 0) 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 with in 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 & \+ Extents" }}{PARA 0 "" 0 "" {TEXT -1 45 "Represent the line in with the implicit form " }{XPPEDIT 18 0 "alpha*x+beta*y = rho;" "6#/,&*&%&alph aG\"\"\"%\"xGF'F'*&%%betaGF'%\"yGF'F'%$rhoG" }{TEXT -1 95 " and then c ompute the for points where the infinite line crosses boundaries of th e 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#/ %\"yG,$*&,(%\"zG!\"\"*&%&alphaG\"\"\"%\"xGF,F,%&gammaGF,\"\"\"%%betaG! \"\"F)" }}}}{SECT 0 {PARA 20 "" 0 "" {TEXT -1 25 "Procedures for xofy, etc." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 183 "xofy := proc(y,X,a vg)\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 0 {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, y min, xmax, ymax:\n xmin := ext[1]: ymin := ext[2]: xmax := ext[3]: \+ ymax := ext[4]:\n [[xmin,ymin,zofxy(xmin,ymin,X)],\n [xmin,ymax,z ofxy(xmin,ymax,X)],\n [xmax,ymax,zofxy(xmax,ymax,X)],\n [xmax,ym in,zofxy(xmax,ymin,X)]]:\nend:\n" }}}}{SECT 0 {PARA 20 "" 0 "" {TEXT -1 23 "Line Segment from Plane" }}{PARA 0 "" 0 "" {TEXT -1 227 "Check \+ first for vertical or horizontal special case. Otherwise, intersect in finite 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." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1106 "lineFromPla ne := proc(X,ext,avg)\n local Line, xt, yt, xmin, ymin, xmax, ymax: \n xmin := 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,avg)], [xmax,yofx(xmax,X,avg)]]:\n elif (X[2,1] = 0) then # H orizontal line\n Line := [[xofy(ymin,X,avg),ymin],[xofy(ymax,X,avg )]]:\n else\n Line := []:\n # Test if line crosses left bo undary at ymin\n xt := xofy(ymin,X,avg):\n if ((xt > xmin) a nd (xt < xmax)) then\n Line := [op(Line), [xt,ymin]]:\n f i:\n # Test if line crosses right boundary at ymax\n xt := x ofy(ymax,X,avg):\n if ((xt > xmin) and (xt < xmax)) then\n \+ Line := [op(Line), [xt,ymax]]:\n fi:\n # Test if line cros ses top boundary 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)) the n\n Line := [op(Line), [xmax,yt]]:\n fi:\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, \+ each edge support region, can be said to be supported by some pixels w ithin it. To see what this means, realize that every pixel in the orig inal image belongs to two edge support regions, one for each of the tw o 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 pi xels voting for a region over the total number of pixels in that regio n. To compute support, first the size of each region must be computed. \n" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 697 "lineSupport := proc( Size,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) then\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 len 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 0 {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 817 "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 #print(\"Line & Plane\", LN, PC, \" Supp \", Supp[ri]):\n RES[ri] := [ri, LN, line3D(LN,avg), RL[2][r i], PC, RL[3][ri], Supp[ri]]:\n od: #End of iteration through regi ons\n RES:\nend:" }}}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 22 "The Comp lete Algorithm" }}{PARA 0 "" 0 "" {TEXT -1 115 "Here are the pieces pu t together. The input is a single image array. The output is a sequenc e of 2D line segments. " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "A := testImage():" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 422 "burns \+ := proc(A)\n local Sh, Sv, Gs, Ls, Rs, RGL:\n Sh := matrix(3,3, [- 1,-2,-1,0,0,0,1,2,1]):\n Sv := transpose(Sh):\n Gs := GradImages(A ,Sh,Sv):\n #print(A):\n #print(Gs[1]):\n Ls := bucketLabel(Gs): \n Rs := [connectedComponents(Ls[1]), connectedComponents(Ls[2])]:\n #print(Rs):\n RGL := makeRegList(Rs,3): \n #print(RGL[2],RGL[3] ):\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 Algori thm" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "A := testImage():" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "LS := burns(A);" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%#LSG-%'vectorG6#7-7)\"\"\"7$7$\"\"$$\"+*3> Q^$!\"*7$\"\"($\"+sxoyMF07$7%F-F.$\"+\\(\\4a\"!\"(7%F2F3F77&F-F-F2\"\" %7&7%F-F-$\"++>2D5F97%F-F<$\"+!oy!H?F97%F2F<$\"+g2Nk?F97%F2F-$\"+!)RMg 5F9F2#\"\"#F27)FK7$7$$\"+\"f#[0IF0\"\"&7$$\"+/\\:4WF0\"\"*7$7%FOFQ$\"+ p#oAX\"F97%FSFUFX7&F-FQFQFU7&7%F-FQ$\"+g/j\\9F97%F-FU$\"*O#*Gu(F97%FQF U$\"+c$Qlt\"F97%FQFQ$\"+!ez=T#F9FU#FKFU7)F-7$7$F<$\"+sd.R\")F07$F2$\"+ tF`t7%FbpFU$\"&s0)F`t7%FbpFbs$\"(y4E)F`tF2#F-F27)Fbs7$7$F-$!+)=kX:$!# 57$FbsF^u7$7%F-F^u$\"+>4e[5F97%FbsF^uFdu7&F-FKFbsFK7&7%F-FK$\"$;#\"\"! Fiu7%FbsFKFjuF]vF ()F0F^r7$F^r$\"+P8+U))F07$7%F^]lF^r$\"+!zkXR\"F97%F^rFa]lFe]l7&F^rF^rF UFbp7&7%F^rF^r$\"*]&Rw=F`t7%F^rFbp$\")K1>tF`t7%FUFbp$\"(;8='F`t7%FUF^r $\"*M-j?\"F`tF " 0 "" {MPLTEXT 1 0 387 "drawLines2D := 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) the n\n objs := \{op(objs),CURVES(LS[i][2],THICKNESS(3),COLOR(HUE, 0.5))\}:\n fi:\n od:\n display(objs,axes=FRAME,view=[0..rs,0. .cs],scaling=CONSTRAINED,font=[HELVETICA,BOLD,14]);\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 575 "drawLines3D := proc(LS,A)\n loca l objs, len, i, j, rs, cs:\n rs := rowdim(A): cs := coldim(A):\n \+ objs := \{\}:\n len := rhs(op(2,eval(LS))):\n for i from 1 to rs \+ do\n for j from 1 to cs do\n objs := \{op(objs),POINTS([c onvert(i,float),convert(j,float),A[i,j]],COLOR(HUE,0.1),SYMBOL(CIRCLE) )\}:\n od:\n od:\n for i from 1 to len do\n if (LS[i][7] > 1/2) then\n objs := \{op(objs),CURVES(LS[i][3],THICKNESS(3) ,COLOR(HUE,0.7))\}:\n fi:\n od:\n # print(objs):\n display( objs,axes=FRAME,view=[0..rs,0..cs,0..255],font=[HELVETICA,BOLD,14]);\n end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 113 "drawTest := proc() \n local A, LS:\n A := testImage():\n LS := burns(A):\n drawLi nes3D(LS,A);\nend:\ndrawTest();" }}{PARA 13 "" 1 "" {GLPLOT3D 400 300 300 {PLOTDATA 3 "6dt-%'POINTSG6%7%$\"#6\"\"!$\"#7F)$\"%q5!\"\"-%&COLOR G6$%$HUEG$\"\"\"F.-%'SYMBOLG6#%'CIRCLEG-F$6%7%F*$\"\"#F)$\"%I5F.F/F5-% 'CURVESG6%7$7%$\"+/*H\\:)!\"*\"\"%$\"++XU5;!\"(7%$\"+\"4^hH)FG\"\")FI- %*THICKNESSG6#\"\"$-F06$F2$\"\"(F.-F$6%7%$\"\"&F)F*$\"%g5F.F/F5-F$6%7% $\"#5F)$\"\"*F)$\"%!3\"F.F/F5-F$6%7%Fen$FHF)$\"%I?F.F/F5-F$6%7%$\"\"'F )F^oF`oF/F5-F$6%7%FeoF\\o$\"%!4\"F.F/F5-F$6%7%F[pF\\o$\"%?5F.F/F5-F$6% 7%F\\oF<$\"%+6F.F/F5-F$6%7%$FOF)F[p$\"%??F.F/F5-F$6%7%Feo$FWF)$\"%5?F. F/F5-F$6%7%FF/F5-F$6%7%F*F[pFjpF/F5-F$6%7%F*FenF[rF/ F5-F$6%7%F*FeoFfsF/F5-F$6%7%F'F\\oF,F/F5-F$6%7%F'F'FfsF/F5-FA6%7$7%$\" +%)y)*=jFGF_o$\"+kIVw5FK7%FO$\"+_/SG5!\")F^wFPFT-FA6%7$7%$\"+[QHUPFGFH $\"+5W#>e\"FK7%$\"+(*QwxPFGFOFjwFPFT-F$6%7%F*F*$\"%+5F.F/F5-F$6%7%F*F \\oFjpF/F5-F$6%7%F*FeqF[rF/F5-F$6%7%F'F[pFbxF/F5-F$6%7%F'FeqFjpF/F5-F$ 6%7%F'FenF`oF/F5-F$6%7%F'FeoF[rF/F5-F$6%7%F'FF/F5-F$6%7%F^oF_qFbxF/F5-F$6%7%F^oFbsF`rF/F5-F$6%7%F ^oFeoF>F/F5-F$6%7%F^oFenFbxF/F5-F$6%7%F^oF'FbxF/F5-F$6%7%F^oF\\oFfsF/F 5-F$6%7%FeqFF/F5-F$6%7%F_qF_qFfqF/F5-F$6%7%F_qFeqFa]lF/F5-F$6%7%F^o FhrFfsF/F5-F$6%7%F^oFINg;F K7%FO$\"+Ai\\VTFGFjdlFPFT-F$6%7%F[pFeoFfqF/F5-F$6%7%FenF\\oFjpF/F5-F$6 %7%FenF'F>F/F5-F$6%7%F[pFhrF`oF/F5-F$6%7%F[pF[p$\"%+@F.F/F5-F$6%7%F[pF eq$\"%g?F.F/F5-F$6%7%F[pF_qF`qF/F5-F$6%7%F[pFenF`qF/F5-F$6%7%FenF^oF[r F/F5-F$6%7%F[pFbsF`oF/F5-F$6%7%FenFF/F5-F$6%7%FbsF\\oF>F/F5 -F$6%7%FeoF'FgnF/F5-F$6%7%FeoF*F[rF/F5-F$6%7%FeoF^oF[rF/F5-F$6%7%FeoF_ qF`qF/F5-F$6%7%FeoFenFcflF/F5-F$6%7%FeoF[pFdilF/F5-F$6%7%FeoFbsFepF/F5 -F$6%7%F " 0 " " {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}}{MARK "10 4 1 0" 0 }{VIEWOPTS 1 1 0 3 2 1804 }