(******************************************************************************
** IncrementalHull.sml
** sml
**
** Guy Blelloch
** Simple Incremental Convex Hull
******************************************************************************)
functor IncrementalHull (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
type faceData = Geometry.HalfSpace.t * Vertex.vertex list
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)
(*******************************************************************
********************************************************************
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
fun makePlane points =
let
val _ = makePlaneCount := !makePlaneCount + 1
in
Geometry.HalfSpace.fromPoints(Sequence.map Vertex.location points)
end
fun planeDist(plane,point) =
let
val _ = lineSideCount := !lineSideCount + 1
val loc = Vertex.location(point)
in
Geometry.HalfSpace.pointDist plane loc
end
fun planeSide(plane,point) =
let
val _ = lineSideCount := !lineSideCount + 1
val loc = Vertex.location(point)
open Number
in
(Geometry.HalfSpace.pointDist plane loc) > zero
end
(*******************************************************************
********************************************************************
Proc: searchFirst()
Description: Search the hull for simplices that are visible from a given
vertex and remove them.
Input:
a) The current convex hull
b) A simplex on the hull
c) A light source (vertex) visible from the simplex
Output:
a) A new convex hull with the visible region deleted
b) The boundary (horizon) of the visible region
c) A list of all the datavertices taken from the visible region
Algorithm:
Depth first search across the visible region.
Whenever reaching a visible simplex
1) the data vertices are taken off of it and added to running list of data
2) the simplex is removed
3) the other outgoing edges are searched
When reaching a boundary (nonvisible simplex) the boundary is added to
the running list of boundary edges
********************************************************************)
fun searchFirst (hull, simplex, light) =
let
(* Recursive depth first search to find visible simplices *)
fun search (arc, (hull, boundary, data)) =
let
val _ = searchCount := !searchCount + 1
in
case Complex.find(hull, Simplex.flip(arc)) of
NONE => (hull, boundary, data)
| SOME(simplex) =>
let
val SOME(plane,newdata) = Complex.getData(hull, simplex)
in
if planeSide(plane,light) then
let
val data = newdata @ data
val _ = remCount := !remCount + 1
val hull = Complex.rem(hull,simplex)
val neighborArcs = Simplex.faces(simplex)
(* the following optimization depends on the
order of Simplex.faces and is commented out *)
(* val neighborArcs = Sequence.drop 1 neighborArcs *)
in
Sequence.foldr search (hull, boundary,data) neighborArcs
end
else
(hull, arc::boundary, data)
end
end
val neighborArcs = Simplex.faces(simplex)
val _ = remCount := !remCount + 1
val hull = Complex.rem(hull,simplex)
in
Sequence.foldr search (hull, nil, nil) neighborArcs
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
********************************************************************)
fun lightOrDark (plane, vertices) =
let
fun select (point, (maxLight, light, maxDark, dark)) =
let
val dist = planeDist(plane,point)
open Number
in
if dist > zero then
case maxLight of
NONE => (SOME(dist,point), light, maxDark, dark)
| SOME(d,pt) => (if dist > d then
(SOME(dist,point), pt::light, maxDark, dark)
else
(maxLight, point::light, maxDark, dark))
else
case maxDark of
NONE => (maxLight, light, SOME(dist,point), dark)
| SOME(d,pt) => (if dist < d then
(maxLight, light, SOME(dist,point), pt::dark)
else
(maxLight, light, maxDark, point::dark))
end
val (maxLight,light,maxDark,dark) = List.foldr select (NONE,nil,NONE,nil) vertices
fun addExtreme(NONE,pts) = pts
| addExtreme(SOME(_,pt),pts) = pt::pts
in
(addExtreme(maxLight,light),addExtreme(maxDark,dark))
end
(*******************************************************************
********************************************************************
Proc: tent
Description: Build a new tent on the hull around a given apex
Input:
a) The convex hull
b) The boundary of a region that is missing (a hole) -- a list of its edges
this boundary will form the rim of the tent
c) a vertex used as the apex of the tent
d) The data -- a list of vertices to be added to each simplex of the tent if
they are visible from that simplex
Output:
a) A new convex hull with the visible region deleted
b) A list of the vertices (data) that were not assigned to any simplex of the tent
c) A list of all the simplices of the tent that have new data on them
Algorithm:
Circles around the boundary adding one simplex at a time
********************************************************************)
fun tent (hull, boundary, apex, data) =
let
(* adds one simplex to the tent *)
fun addsimplex(arc, (hull, data, q)) =
let
(* adds the apex to one of the boundary edges *)
val simplex = Simplex.join(apex, arc)
(* Split into the light and dark vertices *)
val plane = makePlane(Simplex.vertices(simplex))
val (light, dark) = lightOrDark(plane, data)
val _ = addCount := !addCount + 1
(* add the simplex to the hull using filtered data *)
val hull = Complex.add(hull, simplex)
val hull = Complex.addData(hull, simplex, (plane,light))
in
case light of
nil => (hull, dark, q)
| _ => (hull, dark, simplex::q)
end
in
List.foldr addsimplex (hull, data, nil) boundary
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 orders. 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
val p1 = makePlane(Simplex.vertices(t1))
val p2 = makePlane(Simplex.vertices(t2))
(* Split into the light and dark vertices *)
val (d1, d2) = lightOrDark(p1, rest)
val hull = Complex.new(dimension - 1)
val hull = Complex.add(hull,t1)
val hull = Complex.addData(hull,t1,(p1,d1))
val hull = Complex.add(hull,t2)
val hull = Complex.addData(hull,t2,(p2,d2))
in
(hull, [t1,t2])
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 safeGetData(hull,simplex) =
Complex.getData(hull,simplex) handle
SimplexNotFound => NONE
fun hull points =
let
fun ch(hull, nil) = hull
| ch(hull, simplex::queue) =
case safeGetData(hull,simplex) of
NONE => ch(hull, queue)
| SOME(plane,nil) => ch(hull, queue)
| SOME(plane,light::rest) =>
let
val (hull, boundary, data) = searchFirst(hull, simplex, light)
val data = rest @ data
val (hull, dark, q) = tent(hull, boundary, light, data)
in
ch(hull, queue @ q)
end
val _ = clearCount()
val points = List.map Vertex.new points
val dimension = 3
val (hull, queue) = start(points, dimension)
in
ch(hull, queue)
end
end (* IncrementalHull *)