(******************************************************************************
 ** 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 *)