MODULE Bezier;

(* Bezier curve predicates and procedures. *)

IMPORT R2, Geometry;

(* A Bezier curve is defined by a univariate cubic polynomial 
   "F(t)": 

|  F(t) = p0 + t*p1 + t^2*p2 + t^3*p3,

   where "t" is a scalar, and the coefficients "p_i" are vectors 
   (pairs of real numbers). Here, "*" denotes the multiplication 
   of a vector by a scalar, and "+" denotes vector addition. 

   Alternatively, a Bezier curve "a", "b", "c", "d" may be 
   specified by its two endpoints "a" and "d" and its two 
   control points "b" and "c". Since "t = 0" corresponds to the 
   first endpoint "a", and "F(0) = p0", we have that "a = p0". 
   The "MapRep" predicate below defines the mapping between the 
   two representations. *)

(* MAPPING BETWEEN REPRESENTATIONS *)

PRED MapRep(a, b, c, d, p1, p2, p3) IS 
  (E ba, bc, b2, ab2, bc3, ab2c, bc3d :: 
    ba = R2.Minus(b, a) AND 
    bc = R2.Minus(b, c) AND 
    b2 = R2.Times(2, b) AND 
    ab2 = R2.Minus(a, b2) AND 
    bc3 = R2.Times(3, bc) AND 
    ab2c = R2.Plus(ab2, c) AND 
    bc3d = R2.Plus(bc3, d) AND 
    p1 = R2.Times(3, ba) AND 
    p2 = R2.Times(3, ab2c) AND 
    p3 = R2.Minus(bc3d, a)) 
END;

(* The curve defined by the points "a", "b", "c", and "d" 
   corresponds to the univariate cubic polynomial "F" with 
   coefficients "a", "p1", "p2", and "p3". *)

(* COEFFICIENT REPRESENTATION *)

PRIVATE FUNC r = Sum3(a, b, c) IS 
  r = R2.Plus(a, R2.Plus(b, c)) 
END;

/* "r" is the vector sum of "a", "b", and "c", which are 2D 
   points. */

PRIVATE FUNC r = Sum4(a, b, c, d) IS 
  r = R2.Plus(a, R2.Plus(b, R2.Plus(c, d))) 
END;

/* "r" is the vector sum of "a", "b", "c", and "d", which are 2D 
   points. */

FUNC p = AtTFromCoeffs(p0, p1, p2, p3, t) IS 
  (E t2, t3, s1, s2, s3 :: 
    t2 = t * t AND 
    t3 = t * t2 AND 
    s1 = R2.Times(t, p1) AND 
    s2 = R2.Times(t2, p2) AND 
    s3 = R2.Times(t3, p3) AND 
    p = Sum4(p0, s1, s2, s3)) 
END;

(* If "F" is the Bezier curve defined by the four coefficients 
   "p0", "p1", "p2", and "p3", then "p = F(t)". *)

FUNC p = PrimeAtTFromCoeffs(p0, p1, p2, p3, t) IS 
  p = Sum3(p1, R2.Times(2 * t, p2), R2.Times(3 * t * t, p3)) 
END;

(* If "F" is the Bezier curve defined by the four coefficients 
   "p0", "p1", "p2", and "p3", then "p = F'(t)", where "F'" 
   denotes the first derivative of the polynomial "F" with 
   respect to "t". *)

(* POINT REPRESENTATION *)

FUNC p = AtT(a, b, c, d, t) IS 
  (E p1, p2, p3 :: 
    MapRep(a, b, c, d, p1, p2, p3) AND 
    p = AtTFromCoeffs(a, p1, p2, p3, t)) 
END;

(* If "F" is the Bezier curve defined by the four points "a", 
   "b", "c", and "d", then "p = F(t)". *)

FUNC p = PrimeAtT(a, b, c, d, t) IS 
  (E p1, p2, p3 :: 
    MapRep(a, b, c, d, p1, p2, p3) AND 
    p = PrimeAtTFromCoeffs(a, p1, p2, p3, t)) 
END;

(* If "F" is the Bezier curve defined by the four points "a", 
   "b", "c", and "d", then "p = F'(t)", where "F'" denotes the 
   first derivative of the polynomial "F" with respect to "t". *)

PRIVATE PROC p := EvalAtT(p0, p1, p2, p3, t) IS 
  VAR t2, t3, s1, s2, s3 IN 
    t2 := t * t; 
    t3 := t * t2; 
    s1 := R2.Times(t, p1); 
    s2 := R2.Times(t2, p2); 
    s3 := R2.Times(t3, p3); 
    p := Sum4(p0, s1, s2, s3) 
  END 
END;

PROC cl := Closure(a, b, c, d) IS 
  IF 
    VAR p1, p2, p3 IN 
      MapRep(a, b, c, d, p1, p2, p3) -> 
        cl := CLOSE(EvalAtT, a, p1, p2, p3) 
    END 
  FI 
END;

(* Return a closure taking a single argument "t" and returning 
   the point "F(t)", where "F" is the Bezier curve defined by 
   the points "a", "b", "c", and "d". *)

PRIVATE FUNC t = EstimateT(p, a, b, c, d) IS 
  (E ab = (0.5, 0) REL (a, b), bc = (0.5, 0) REL (b, c), 
     cd = (0.5, 0) REL (c, d), abc = (0.5, 0) REL (ab, bc), 
     bcd = (0.5, 0) REL (bc, cd), 
     abcd = (0.5, 0) REL (abc, bcd), x, y :: 
    p = (x, y) REL (abcd, bcd) AND t = 0.5 + x / 6) 
END;

PRED On(p, a, b, c, d) IS 
  (E t ~ EstimateT(p, a, b, c, d), ab ~ (t, 0) REL (a, b), 
     bc ~ (t, 0) REL (b, c), cd ~ (t, 0) REL (c, d), 
     abc ~ (t, 0) REL (ab, bc), bcd ~ (t, 0) REL (bc, cd) :: 
    ab = (t, 0) REL (a, b) AND 
    bc = (t, 0) REL (b, c) AND 
    cd = (t, 0) REL (c, d) AND 
    abc = (t, 0) REL (ab, bc) AND 
    bcd = (t, 0) REL (bc, cd) AND 
    p = (t, 0) REL (abc, bcd)) 
END;

UI PointTool(On);

(* The point "p" lies on the Bezier curve "a", "b", "c", "d". *)

PROC p1, p2, p3, p4, p5 := Split(a, b, c, d) IS 
  IF 
    VAR 
      ab = Geometry.Mid(a, b), bc = Geometry.Mid(b, c), 
      cd = Geometry.Mid(c, d), abc = Geometry.Mid(ab, bc), 
      bcd = Geometry.Mid(bc, cd), abcd = Geometry.Mid(abc, bcd) 
    IN 
      p1, p2, p3, p4, p5 := ab, abc, abcd, bcd, cd 
    END 
  FI 
END;

(* Set the output parameters so that 

|    PS.MoveTo(a);
|    PS.CurveTo(p1, p2, p3);
|    PS.CurveTo(p4, p5, d)

   produces a path equivalent to 

|    PS.MoveTo(a);
|    PS.CurveTo(b, c, d).

   The point "p3" is the midpoint of the curve "a", b", "c", 
   "d". That is, "p3 = AtT(a, b, c, d, 0.5)". *)

PRIVATE CONST MinDist = 8, MaxRatio = 0.004;

PRIVATE PROC Ck, Pk := LengthAux(a, b, c, d, depth, max) IS 
  VAR split = 1 IN 
    IF 
      depth < max -> SKIP 
    | Ck, Pk := 
        Geometry.Dist(a, d), 
        Geometry.Dist(a, b) + Geometry.Dist(b, c) + 
        Geometry.Dist(c, d); 
      IF 
        Pk < MinDist OR (Pk - Ck) / Pk < MaxRatio -> split := NIL 
      | SKIP 
      FI 
    FI; 
    IF 
      split # NIL -> 
        VAR p1, p2, p3, p4, p5, Ck2, Pk2 IN 
          p1, p2, p3, p4, p5 := Split(a, b, c, d); 
          Ck, Pk := LengthAux(a, p1, p2, p3, depth + 1, max); 
          Ck2, Pk2 := LengthAux(p3, p4, p5, d, depth + 1, max); 
          Ck, Pk := Ck + Ck2, Pk + Pk2 
        END | SKIP 
    FI 
  END 
END;

/* Sets "Ck" and "Pk" to a lower- and upper-bound, respectively, 
   on the length of the Bezier curve defined by the points "a", 
   "b", "c", and "d". Recursion proceeds so long as the relative 
   difference between the upper and lower bound is greater than 
   the constant "Epsilon". */

PROC len := Length(a, b, c, d) IS 
  VAR Ck, Pk IN 
    Ck, Pk := LengthAux(a, b, c, d, 0, 3); len := (Ck + Pk) / 2 
  END 
END;

(* Set "len" to an approximation of the length (in points) of 
   the Bezier curve defined by the four points "a", "b", "c", 
   "d". *)