(* spline.sml
 *
 * COPYRIGHT (c) 1991 by AT&T Bell Laboratories  See COPYRIGHT file for details.
 *)

structure Spline : SPLINE =
  struct

    structure G = Geometry
    open G

    fun round x = if x > 0.0 then floor (x+0.5) else ~1*floor(~x+0.5)

    fun addSeg ([], x0,y0,x1,y1) = 
          [PT{x = round x0, y = round y0}, PT{x = round x1, y = round y1}]
      | addSeg (l, x0, y0, _,_) = (PT{x = round x0, y = round y0})::l

  (* isFlat:
   * Returns true if the polygon determined by the four points
   * is flat enough. Flatness is measured by the maximum distance
   * of (x1,y1) and (x2,y2) from the line determined by (x0,y0)
   * and (x3,y3). In addition, check that p1, p2 are close to the 
   * line segment. To do this, make sure they are roughly within 
   * the circle with center (p0+p3)/2 and radius = |p3-p0|/2+flatness
   *)
    fun isFlat {x0, y0, x1, y1, x2, y2, x3, y3} = let
	  fun sqr x = x * x
	  val dx = x3 - x0 and dy = y3 - y0
	  val midx = 0.5*dx and midy = 0.5*dy and dist2 = sqr dy + sqr dx
	  val flatness2 = sqr 1.0 * dist2 and halfd2 = 0.25*dist2
	  fun inFlatRange (x,y) = sqr(dy * x - dx * y) <= flatness2
		andalso let val d = sqr(x-midx) + sqr(y-midy)
			in d<=halfd2 orelse sqr(d-halfd2)<=flatness2 end
	  in
	    inFlatRange (x1-x0, y1-y0) andalso inFlatRange (x2-x0, y2-y0)
	  end

  (* bezier:
   * Recursively compute a Bezier cubic section. If the points
   * determine a polygon flat enough to be represented as a line
   * segment, the segment is added to the list. Otherwise, the
   * the curve is bisected and each part is recursively computed,
   * with the lists concatenated.
   *
   * From "The Beta2-split: A special case of the Beta-spline Curve and
   * Surface Representation." B. A. Barsky and A. D. DeRose. IEEE, 1985,
   * as adapted by Crispin Goswell for xps.
   *)
    fun bezier (arg as {x0,y0,x1,y1,x2,y2,x3,y3}, l) = if (isFlat arg)
	  then addSeg(l, x0, y0, x3, y3)
	  else let
	    val mid_x = (x0 + x3) / 8.0 + 3.0 * (x1 + x2) / 8.0
	    val mid_y = (y0 + y3) / 8.0 + 3.0 * (y1 + y2) / 8.0
	    val l' = bezier ({
		    x0 = mid_x,                     y0 = mid_y,
		    x1 = (x1+x3) / 4.0 + x2 / 2.0,  y1 = (y1+y3) / 4.0 + y2 / 2.0,
		    x2 = (x2+x3) / 2.0,             y2 = (y2+y3) / 2.0,
		    x3 = x3,                        y3 = y3
		  }, l)
	    in
	      bezier ({
		  x0 = x0,                        y0 = y0,
		  x1 = (x0+x1) / 2.0,             y1 = (y0+y1) / 2.0,
		  x2 = (x0+x2) / 4.0 + x1 / 2.0,  y2 = (y0+y2) / 4.0 + y1 / 2.0,
		  x3 = mid_x,                     y3 = mid_y
		}, l')
	    end

  (* curve:
   * Given four points [p0,p1,p2,p3], return a list of points corresponding to 
   * to a Bezier cubic section, starting at p0, ending at p3, with p1, p2 as
   * control points.
   *
   *)
    fun curve (PT{x=x0,y=y0}, PT{x=x1,y=y1}, PT{x=x2,y=y2}, PT{x=x3,y=y3}) =
	  bezier ({
	      x0 = real x0, y0 = real y0, 
	      x1 = real x1, y1 = real y1, 
	      x2 = real x2, y2 = real y2, 
	      x3 = real x3, y3 = real y3
	    }, [])

  (* doSpline:
   * Given four points (p0,p1,p2,p3), return a list of points corresponding to 
   * to the B-spline curve section, accumulating the results on the argument list.
   * We compute the curve by determining the corresponding Bezier control points,
   * and then use the Bezier routines above.
   *
   *)
    fun doSpline (PT{x=x0,y=y0}, PT{x=x1,y=y1}, PT{x=x2,y=y2}, PT{x=x3,y=y3}, l) = let
	  val x0 = real x0 and y0 = real y0 
	  val x1 = real x1 and y1 = real y1 
	  val x2 = real x2 and y2 = real y2 
	  val x3 = real x3 and y3 = real y3
	  in
	    bezier ({
		x0 = (x0 + 4.0*x1 + x2)/6.0, y0 = (y0+ 4.0*y1 + y2)/6.0, 
		x1 = (2.0*x1 + x2)/3.0,      y1 = (2.0*y1 + y2)/3.0, 
		x2 = (x1 + 2.0*x2)/3.0,      y2 = (y1 + 2.0*y2)/3.0, 
		x3 = (x1 + 4.0*x2 + x3)/6.0, y3 = (y1 + 4.0*y2 + y3)/6.0
	      }, l)
	  end

  (* loopSpline:
   * Given a list of spline control points, generate the corresponding
   * spline. Since we accumulate on the head of the list, we assume
   * the calling function has reversed the list of control points.
   * The loop continues as long as there are 4 control points left.
   *)
    fun loopSpline (p3, p2, p1, p0::tl, l) = 
          loopSpline (p2, p1, p0, tl, doSpline(p0, p1, p2, p3, l))
      | loopSpline (_, _, _, [], l) = l 

  (* simpleBSpline:
   * Compute a simple B-spline with the given control points.
   *)
    fun simpleBSpline arg = (case (rev arg)
	   of (p3::p2::p1::tl) => loopSpline (p3,p2,p1,tl,[])
	    | _ => arg
	  (* end case *))
          
  (* bSpline:
   * Compute a B-spline using the given control points.
   * In addition, we constrain the resultant spline to connect the
   * first and last points by adding copies of these points.
   *)
    fun bSpline (arg as (p0::_::_::_)) = let
	  val (pn::tl) = rev (p0::p0::arg)
	  in
	    loopSpline (pn, pn, pn, tl, [])
	  end
      | bSpline l = l
    
  (* closedBSpline:
   * Compute a closed B-spline. This is done by repeating the first
   * three points at the end of the list.
   * Note that the first and last points of the result are the same.
   *)
    fun closedBSpline (arg as (p0::p1::p2::_)) = loopSpline(p2,p1,p0,rev arg,[])
      | closedBSpline l = l

  end (* Spline *)