(****************************************************************************** ** BulldozerHull.sml ** sml ** ** Hal Burch ** Bulldozer Convex Hull ******************************************************************************) structure Queue = struct datatype 'a WaitQueue = Empty | Leaf of int * 'a | Node of int * 'a WaitQueue * 'a WaitQueue exception EmptyQueue val qcnt = ref 0 fun queueSize Empty = 0 | queueSize (Leaf(c,_)) = c | queueSize (Node(c,_,_)) = c fun enqueue Empty simp cnt = Leaf(cnt,simp) | enqueue (Leaf(c,s)) simp cnt = Node(c+cnt,Leaf(c,s),Leaf(cnt,simp)) | enqueue (Node(c,l,r)) simp cnt = let in if (queueSize(l) > queueSize(r)) then Node(c+cnt,l,enqueue r simp cnt) else Node(c+cnt,enqueue l simp cnt,r) end fun findPos Empty num = raise EmptyQueue | findPos (Leaf(c,s)) num = (s,Empty,num) | findPos (Node(c,l,r)) num = let in if (num < queueSize(l)) then let val (s,nt,ind) = findPos l num in case nt of Empty => (s, r, ind) | Leaf(c,_) => (s, Node(c + queueSize r,nt,r), ind) | Node(c,_,_) => (s, Node(c + queueSize r,nt,r), ind) end (* if num < queueSize(l) let *) else let val (s,nt,ind) = findPos r (num - queueSize(l)) in case nt of Empty => (s, l, ind) | Leaf(c,_) => (s, Node(c + queueSize l,l,nt), ind) | Node(c,_,_) => (s, Node(c + queueSize l,l,nt), ind) end (* if num < queueSize(l) else let *) end (* findPos Node num *) fun dequeue Empty seed = raise EmptyQueue | dequeue (Leaf(c,s)) seed = let val (num, seed) = RC4UnitReal.nextSeed seed val num = Real.trunc(num * real(c)) in (s,Empty,seed,num) end (* dequeue Leaf seed let *) | dequeue (Node(c,l,r)) seed = let val (num, seed) = RC4UnitReal.nextSeed seed val num = Real.trunc(num * real(c)) val (s,nt,ind) = findPos (Node(c,l,r)) num in (s,nt,seed,ind) end (* dequeue Node seed let *) end (* Queue *) functor BulldozerHull(structure Geometry : GEOMETRY structure Complex : SIMPLICIAL_COMPLEX where type Vertex.point = Geometry.Point.t) : HULL where Geometry = Geometry and Complex = Complex = struct structure Geometry = Geometry structure Complex = Complex structure Number = Geometry.Number structure Simplex = Complex.Simplex structure Vertex = Complex.Vertex val makePlaneCount = ref 0 val lineSideCount = ref 0 val searchCount = ref 0 val addCount = ref 0 val remCount = ref 0 fun clearCount () = (lineSideCount := 0; makePlaneCount := 0; searchCount := 0; addCount := 0; remCount := 0) datatype 'a group = EMPTY | ELEMENT of 'a * 'a group | SPLIT of 'a group * 'a group type faceData = Vertex.vertex group fun groupFoldr (f:('a * 'b -> 'b)) (base:'b) (EMPTY:'a group) = base | groupFoldr f base (ELEMENT(a,b)) = groupFoldr f (f (a,base)) b | groupFoldr f base (SPLIT(a,b)) = groupFoldr f (groupFoldr f base a) b fun groupJoin (EMPTY,b) = b | groupJoin (a,EMPTY) = a | groupJoin (a,b) = SPLIT(a,b) fun groupAdd (a,b) = ELEMENT(a,b) fun makeGroup list = List.foldr groupAdd EMPTY list fun groupSize0 EMPTY = true | groupSize0 (ELEMENT(_,b)) = false | groupSize0 (SPLIT(a,b)) = case (groupSize0 a) of true => groupSize0 b | false => false fun groupSize1 EMPTY = false | groupSize1 (ELEMENT(_,b)) = groupSize0 b | groupSize1 (SPLIT(a,b)) = case (groupSize1 a) of true => groupSize0 b | false => case (groupSize0 a) of false => false | true => groupSize1 b fun groupSize EMPTY = 0 | groupSize (ELEMENT(_,b)) = 1 + (groupSize b) | groupSize (SPLIT(a,b)) = (groupSize a) + (groupSize b) exception TooLarge (* debugging stuff *) fun pointToString simp = Geometry.Point.toString(Vertex.location(simp)) fun countPoints hull = let val simps = Complex.simplices hull 2 fun addpts (simp : Simplex.simplex, cnt : int) = case Complex.getData(hull,simp) of NONE => cnt | SOME(a) => cnt + (length(a)) in Sequence.foldl addpts 0 simps end (* countPoints *) fun status hull = Int.toString(countPoints hull) ^ " (" ^ Int.toString(countPoints hull) ^ ")" (******************************************************************* ******************************************************************** proc: lightFace Description: Curried version for finding light side of a face Input: a triangle Output: a function that takes a point and returns true if the point is above the triangle Algorithm: this partially evaluates the plane-side test based on the thre points of a triangle without knowing the fourth point. ********************************************************************) exception lightFaceBadInput val LFcnt = ref 0 val SUcnt = ref 0 val Fcnt = ref 0 (* val flag = ref true val CurList = ref (nil:Vertex.vertex vector list) *) fun distance points = if not(Sequence.length(points) = 3) then raise lightFaceBadInput else let val locations = Sequence.map Vertex.location points val halfspace = Geometry.HalfSpace.fromPoints(locations) in (Geometry.HalfSpace.pointDist halfspace) o Vertex.location end fun lightFace points = if not(Sequence.length(points) = 3) then raise lightFaceBadInput else let val locations = Sequence.map Vertex.location points val halfspace = Geometry.HalfSpace.fromPoints(locations) open Number in (fn pt => (Geometry.HalfSpace.pointDist halfspace (Vertex.location pt)) > zero) end exception NoPoints fun furthest (simplex,vertices:Vertex.vertex group) = let exception NoMatch fun CalcFirst EMPTY = raise NoMatch | CalcFirst (ELEMENT(a,b)) = a | CalcFirst (SPLIT(a,b)) = CalcFirst a handle NoMatch => CalcFirst b in if (groupSize1 vertices) then (CalcFirst vertices, EMPTY) else let val CalcDist = Number.toReal o (distance (Simplex.vertices simplex)) fun findBest (EMPTY,p,d) = (p,d,EMPTY) | findBest (ELEMENT(p1,g),NONE,d) = let val d1 = CalcDist p1 in findBest(g,SOME(p1),d1) end | findBest (ELEMENT(p1,g),SOME(p),d) = let val d1 = CalcDist p1 val (p,d,p1) = if (d1 > d) then (p1,d1,p) else (p,d,p1) val (p',d,g) = findBest(g,SOME(p),d) in (p',d,ELEMENT(p1,g)) end | findBest (SPLIT(EMPTY,g2),p,d) = findBest (g2,p,d) | findBest (SPLIT(g1,EMPTY),p,d) = findBest (g1,p,d) | findBest (SPLIT(g1,g2),p,d) = let val (p,d,g1) = findBest (g1,p,d) val (p,d,g2) = findBest (g2,p,d) in (p,d,SPLIT(g1,g2)) end in case findBest (vertices,NONE,0.0-1.0) of (NONE,_,_) => raise NoMatch | (SOME(p),_,g) => (p,g) end end (******************************************************************* ******************************************************************** Proc: lightOrDark Input: a) a simplex b) a list of vertices Output: a) A list of the vertices visible from the face of the simplex b) A list of the vertices invisible from the face of the simplex ********************************************************************) val LDcnt = ref 0 val LDcntpts = ref 0 fun lightOrDark (simplex, EMPTY) = (EMPTY,EMPTY) | lightOrDark (simplex, vertices:Vertex.vertex group) = let val sideTest = lightFace (Simplex.vertices simplex) (* val _ = (LDcnt := (!LDcnt) + 1) val _ = (LDcntpts := (!LDcntpts) + (groupSize vertices)) *) fun select (point, (light, dark)) = if sideTest(point) then (groupAdd(point,light),dark) else (light,groupAdd(point,dark)) in groupFoldr select (EMPTY,EMPTY) vertices end fun seperateCount sideTest (simplex, vertices:Vertex.vertex group) = let fun select (point, (light, dark,cnt)) = if sideTest(point) then (groupAdd(point,light),dark, cnt+1) else (light,groupAdd(point,dark), cnt) in groupFoldr select (EMPTY,EMPTY,0) vertices end fun seperate (test, vertices:Vertex.vertex group) = let fun select (point, (light, dark)) = if test(point) then (groupAdd(point,light), dark) else (light, groupAdd(point,dark)) in groupFoldr select (EMPTY,EMPTY) vertices end (******************************************************************* ******************************************************************** Proc: start Description: Intialized the Hull with two simplices Input: a) The list of all the vertices to put into the hull Output: a) An initial convex hull consisting of the first d+1 points from the input inserted as a simplex in both orientations. All the other vertices are appropriately assigned to one of the two simplices based on which one they are visible from. b) A list containing the two simplices. Algorithm: Extract first d+1 points, add them and their inverse to the hull ********************************************************************) exception TooFiewVertices fun start(vertices, dimension) = if length(vertices) < dimension then raise TooFiewVertices else let val head = List.take(vertices, dimension) val rest = List.drop(vertices, dimension) val t1 = Simplex.fromSeq(Sequence.fromList(head)) val t2 = Simplex.flip t1 (* Split into the light and dark vertices *) val (d1, d2) = lightOrDark(t1, makeGroup rest) val hull = Complex.new(dimension - 1) val _ = (addCount := !addCount + 2) val hull = Complex.add(hull,t1) val hull = Complex.addData(hull,t1,d1) val hull = Complex.add(hull,t2); val hull = Complex.addData(hull,t2,d2) open Geometry fun sum (vert,vec) = Point.toVec(Point.translate(Vertex.location vert, vec)) val vec = List.foldl sum Vec.zero head val origin = Point.fromVec(Vec.scale(vec, Number.fromReal(1.0/3.0))) val origin = Vertex.new(origin) val queue = Queue.Empty val queue = Queue.enqueue queue t1 (groupSize d1) val queue = Queue.enqueue queue t2 (groupSize d2) in (origin, hull, queue) end exception NoMatch fun addPoint origin hull light visface queue list = let fun PushPoints hull light (p1,p2) face pts queue = case face of NONE => raise NoMatch | SOME(simplex) => let val panSeq = Sequence.fromList [light,p2,p1]; val faceTest = lightFace panSeq val #[v1,v2,v3] = Simplex.vertices simplex val ov = case (Vertex.compare(v1,p1)) of EQUAL => (case (Vertex.compare(v2,p2)) of EQUAL => v3 | _ => v2) | _ => (case (Vertex.compare(v1,p2)) of EQUAL => (case (Vertex.compare(v2,p1)) of EQUAL => v3 | _ => v2) | _ => v1) in if (faceTest ov) then case pts of EMPTY => (hull, queue, true) | _ => let val (newpts,_) = (pts,pts) (* this doesn't seem to be helpful val (newpts,_) = seperate(faceTest, pts) *) val hull = Complex.update (hull, simplex, fn SOME(a) => SOME(groupJoin(newpts,a))) in (hull, queue, true) end (* if let *) else let (* Add Panel *) (* val Fcnt = (Fcnt := (!Fcnt) + 1) *) val newfc = Simplex.fromSeq(panSeq) val _ = (addCount := !addCount + 1) val (pts,_,cnt) = seperateCount faceTest (newfc, pts) val queue = case cnt of 0 => queue | _ => Queue.enqueue queue newfc cnt val hull = Complex.add(hull, newfc) val hull = Complex.addData(hull, newfc, pts) in (hull, queue, false) end (* else let *) end (* case let *) (* recurse these points across the arc (p1,p2) to the given face *) fun RecurseArc hull light (p1,p2) face queue = case face of NONE => raise NoMatch | SOME(simplex) => traverse (hull,simplex,light,(p1,p2),queue) and traverse (hull,curface,light,(p1,p2),queue) = let val verts = Simplex.vertices curface val pts = Complex.getData (hull, curface) val pts = case pts of SOME(d) => d | NONE => raise NoMatch val #[v1,v2,v3] = verts (* find which of these three vertices is not p1 and not p2 *) (* the `other vertex' *) val ov = case (Vertex.compare(v1,p1)) of EQUAL => (case (Vertex.compare(v2,p2)) of EQUAL => v3 | _ => v2) | _ => (case (Vertex.compare(v1,p2)) of EQUAL => (case (Vertex.compare(v2,p1)) of EQUAL => v3 | _ => v2) | _ => v1) val docheck = lightFace (Sequence.fromList([light,origin,ov])) (* val _ = (LDcnt := (!LDcnt) + 1) val _ = (LDcntpts := (!LDcntpts) + 2) *) val check1 = docheck p1 val check2 = case check1 of true => true | _ => docheck p2 in if (check1 = check2) then (* in-degree 2 *) let val otherarc = if (check1) then Simplex.fromSeq(Sequence.fromList([p1,ov])) else Simplex.fromSeq(Sequence.fromList([ov,p2])) val _ = (searchCount := !searchCount + 1) val otherfacedeleted = case (Complex.find(hull, otherarc)) of NONE => true | SOME(_) => false (* val _ = print (Bool.toString(otherfacedeleted) ^ "\n") val _ = print (Bool.toString(check1) ^ "\n") *) in if (otherfacedeleted) then let val (nextarc, arc) = if (check1) then (Simplex.fromSeq(Sequence.fromList([ov,p2])), (ov,p2)) else (Simplex.fromSeq(Sequence.fromList([p1,ov])), (p1,ov)) val _ = (searchCount := !searchCount + 1) val nextface = Complex.find (hull, nextarc) val hull = Complex.rem (hull, curface) val _ = (remCount := !remCount + 1) val (hull, queue, FaceVisible) = PushPoints hull light arc nextface pts queue val (hull, queue) = if (FaceVisible) then RecurseArc hull light arc nextface queue else (hull, queue) in (hull,queue) end else let (* val _ = (LDcnt := (!LDcnt) + 1) val _ = (LDcntpts := (!LDcntpts) + 1) *) in (hull, queue) end end else (* in-degree 1 *) let (* val _ = (LDcntpts := (!LDcntpts) + (groupSize pts)) *) val (rp, lp) = case pts of EMPTY => (EMPTY,EMPTY) | _ => seperate(docheck, pts) val _ = (searchCount := !searchCount + 2) val f1 = Complex.find(hull, Simplex.fromSeq(Sequence.fromList([ov,p2]))) and f2 = Complex.find(hull, Simplex.fromSeq(Sequence.fromList([p1,ov]))) val _ = (remCount := !remCount + 1) val hull = Complex.rem (hull, curface) val (hull,queue,FV1) = PushPoints hull light (ov,p2) f1 rp queue val (hull,queue,FV2) = PushPoints hull light (p1,ov) f2 lp queue val (hull,queue) = if (FV1) then RecurseArc hull light (ov,p2) f1 queue else (hull, queue) val _ = if (FV1 andalso FV2) then searchCount := !searchCount + 1 else () val (hull,queue) = if (FV2) then let (* this could've been removed by the previous RecurseArc call *) val f2 = if (FV1) then Complex.find(hull, Simplex.fromSeq(Sequence.fromList([p1,ov]))) else f2 in case f2 of NONE => (hull, queue) | SOME(_) => RecurseArc hull light (p1,ov) f2 queue end else (hull,queue) in (hull,queue) end end fun startTraverse p1 p2 p3 queue pts = let val (d1,d2) = lightOrDark(Simplex.fromSeq(Sequence.fromList [light,origin,p1]), pts) val (d3,d4) = lightOrDark(Simplex.fromSeq(Sequence.fromList [light,origin,p3]), d1) and (d5,d6) = lightOrDark(Simplex.fromSeq(Sequence.fromList [light,origin,p2]), d2) val _ = (remCount := !remCount + 1) val hull = Complex.rem(hull, visface) val _ = (searchCount := !searchCount + 3) val f1 = Complex.find(hull, Simplex.fromSeq(Sequence.fromList([p2,p1]))) and f2 = Complex.find(hull, Simplex.fromSeq(Sequence.fromList([p3,p2]))) and f3 = Complex.find(hull, Simplex.fromSeq(Sequence.fromList([p1,p3]))) val (hull,queue,FV1) = PushPoints hull light (p2,p1) f1 d5 queue val (hull,queue,FV2) = PushPoints hull light (p3,p2) f2 (groupJoin(d3,d6)) queue val (hull,queue,FV3) = PushPoints hull light (p1,p3) f3 d4 queue val (hull,queue) = if (FV1) then RecurseArc hull light (p2,p1) f1 queue else (hull, queue) (* this could've been removed *) val (hull,queue) = if (FV2) then let (* this could've been removed by the previous RecurseArc call *) val _ = if (FV1) then (searchCount := !searchCount + 1) else () val f2 = if (FV1) then Complex.find(hull, Simplex.fromSeq(Sequence.fromList([p3,p2]))) else f2 in case f2 of NONE => (hull, queue) | SOME(_) => RecurseArc hull light (p3,p2) f2 queue end else (hull,queue) (* this could've been removed *) val _ = (searchCount := !searchCount + 1) val f3 = Complex.find(hull, Simplex.fromSeq(Sequence.fromList([p1,p3]))) val (hull,queue) = if (FV3) then case f3 of SOME(_) => RecurseArc hull light (p1,p3) f3 queue | NONE => (hull,queue) else (hull,queue) in (hull, queue) end val vertices = Simplex.vertices visface val p1 = Sequence.nth (Simplex.vertices visface) 0 and p2 = Sequence.nth (Simplex.vertices visface) 1 and p3 = Sequence.nth (Simplex.vertices visface) 2 in startTraverse p1 p2 p3 queue list end (******************************************************************* ******************************************************************** Proc: hull Input: a) The list of all the vertices to put into the hull Output: a) The convex hull as a simplicial complex Algorithm: Variant of incremental convex hull. Keeps the invariant at each step that a) the current hull is a hull of a subset of the points b) the remaining points that are not contained within the hull are assigned to one of the faces visible from the point Each step a) removes one of the points from a face (currently it picks an arbitrary vertex) b) finds all the faces visible from that point, removes them, and collects their vertices c) builds a new tent into the hull with the point as the apex and the removed boundary as the rim, and assigns each collected vertex to one of the faces from which it is visible. It drops any vertices that are not visible from any face. ********************************************************************) fun hull (points:Vertex.point list) = let val seed = RC4UnitReal.State.start "3DHull" val points = List.map Vertex.new points val dimension = 3 val (origin, hull, queue) = start(points, dimension) fun ch hull Queue.Empty seed = hull | ch (hull:Vertex.vertex group Complex.complex) queue seed = let val (simplex,queue,seed,ind) = Queue.dequeue queue seed in case Complex.getData(hull,simplex) of NONE => ch hull queue seed | SOME(EMPTY) => ch hull queue seed | SOME(list) => let val (light,list) = furthest (simplex, list) (* val light = List.nth (list,ind) *) val (newHull,newQueue) = addPoint origin hull light simplex queue list in ch newHull newQueue seed end end in ch hull queue seed end end (* BulldozerHull *)