(******************************************************************************
 ** BoxFromGeometryPrims.sml
 ** sml
 **
 ** Franklin Chen and Guy Blelloch 
 ** Generates a BOX given a collection of Geometry Primitives
 ******************************************************************************)

functor BoxFromGeometryPrims(structure GP : GEOMETRY_PRIMS
			     structure Point : POINT where type point = GP.Vec.t)
    : BOX =
    struct

      (* This should not be needed any more:
       * exception NYI*)

	structure Point = Point
	structure Number = Point.Number
	structure Vec = Point.Vec

	(* Abbreviations. *)
	structure Seq = Sequence

	(*
	 * Corner and offsets.
	 * This will also be OK in n dimensions.
	 *)
	type box = Point.t * Vec.t
	type t = box

	fun ==((p1,v1),(p2,v2)) = Point.==(p1,p2) andalso Vec.==(v1,v2)
	fun !=((p1,v1),(p2,v2)) = Point.!=(p1,p2) orelse Vec.!=(v1,v2)

	val empty : t = (Point.origin, Vec.zero)

	fun toString((p,v) : t) =
	    ("Box(" ^ Point.toString(p) ^ "," ^ Vec.toString(v) ^ ")")

	fun translate ((p,v) : t ,tr : Vec.t) =
	    (Point.translate(p,tr),v)

		fun fromCorners(p1 : Point.t ,p2 : Point.t) =
		    let
			val lower = GP.minCoords(p1,p2)
			val upper = GP.maxCoords(p1,p2)
		    in
			(lower, Point.-(upper,lower))
		    end

		fun center((p,v) : t) =
		    Point.translate(p,Vec.scale(v,Number.fromReal(0.5)))

		fun radius((p,v) : t) =
		    Number.halve(Vec.norm(v))

		fun toVec((p,v) : t) = v

		fun intersection((p1,v1) : t,(p2,v2) : t) = 
		  let
		    val lower = GP.maxCoords(p1,p2);
		    val upper = GP.minCoords(Point.translate(p1,v1),
					     Point.translate(p2,v2));
		  in
		    if Point.==(lower, GP.minCoords(lower,upper))
		      then SOME(lower, Point.-(upper,lower))
		    else NONE
		  end

                (*
                 * Fold with min and max to get bounds in each dimension.
                 *
                 * Note:  zipping with self because reducel, unlike foldl,
                 * requires use of only one type 'a.
                 *
                 * FMC:  degeneracy (case of 0 or 1 element) may be a problem.
                 *)
                fun fromPoints (ps : Point.t Seq.seq) : t =
                    let
                        open Number

			val posInf = Point.fromVec(Vec.infinity)
			val negInf = Point.fromVec(Vec.~ Vec.infinity)

                        val (p, p') =
                            Seq.reduce (fn ((min1,max1),(min2,max2)) =>
                                         (GP.minCoords(min1,min2),
					  GP.maxCoords(max1,max2)))
                            (posInf,negInf)
                            (Seq.zip (ps, ps))

                    in
                        (p, Point.- (p', p))
                    end

                (*
                 * Check within bounds of each dimension.
                 * SHOULD VERIFY THE DIRECTION
                 *)
                fun pointIn (p,v) p' = 
		  let
		    val lower = p
		    val upper = Point.translate(p,v)
		  in
		    if Point.==(p',GP.maxCoords(p',lower)) andalso
		      Point.==(p',GP.minCoords(p',upper)) then
		      let
			val lo = Point.toSeq lower
			val hi = Point.toSeq upper
			val mid = Point.toSeq p'
		      in
			(* Converts each point to a sequence representing
			 * its position in each dimension; if the comparison
			 * point shares any dimension with either the upper
			 * or lower point of the box, return On *)
			if
			  Seq.foldr (fn (a,b)=>a orelse b) false
			  (Seq.map Number.== (Seq.zip(mid,lo)))
			  orelse
			  Seq.foldr (fn (a,b)=>a orelse b) false
			  (Seq.map Number.== (Seq.zip(mid,hi)))
			  then Side.On else Side.In
		      end
		    else Side.Out
		  end

                (*
                 * Size of the box = maximum width in a dimension
                 *)
                fun size ((_, dp) : t) : Number.t = 
		  Seq.reduce Number.max (Number.fromReal 0.0) (Vec.toSeq dp)

		(* Note from cce@andrew: This function was (as is) in the
		 * code given to me, but it isn't in the signature and
		 * doesn't seem to be needed.
		 *
                fun minDist _ _ =
                    raise NYI
		      *)

		(* Here I (cce@andrew) will use the definition of scatter
		 * annotated by RH,GB as it is unambiguous as to the order
		 * of the 'a seqs*)
		fun split n (p,v) =
		  let
		    val vs = Vec.toSeq v
		    val dim = Seq.sequence(1,(n,Number.halve(Seq.nth vs n)))
		    val (d1,d2) = Seq.split (fn (a,b)=>(a,b)) dim
		    (* New vector for halved boxes *)
		    val vs' = Seq.scatter (Vec.toSeq v, d1,d2)
		    val dp = Vec.fromSeq(Seq.scatter(Vec.toSeq Vec.zero,d1,d2))
		    val v' = Vec.fromSeq vs'
		  in
		    ((p,v'),(Point.translate(p,dp),v'))
		  end
		  
                fun bisect (p,v) = 
		  let
		    val vs = Vec.toSeq v
		    val vs' =
		      Seq.zip(vs,Seq.tabulate (fn i=>i) (Seq.length vs))
		    val (_,idx) =
		      Seq.foldr (fn ((d,i),(dmax,imax)) =>
				 if Number.>(d,dmax) then (d,i)
				 else (dmax,imax))
		      (Number.fromReal 0.0,~1) vs'
		  in
		    split idx (p,v)
		  end

                fun divide b = 
		  let
		    fun append' (a,b) = Seq.append a b
		    (* val div' : t seq -> int -> t seq *)
		    fun div' n s = 
		      (div' (n+1)
		       ((append' o Seq.unzip) (Seq.map (split n) s)))
		      handle Seq.InvalidIndex => s
		  in
		    div' 0 (Seq.sequence(1,b))
		  end
    end