(**

   Bezier curve
*)

module Utils = Tools.Utils

type quadratic =
  { p0:Gg.v2     (* The starting point *)
  ; p1:Gg.v2     (* The end point *)
  ; ctrl:Gg.v2 } (* The control point *)


type t =
  { p0:Gg.v2      (* The starting point *)
  ; p1:Gg.v2      (* The end point *)
  ; ctrl0:Gg.v2   (* The control point *)
  ; ctrl1:Gg.v2 } (* The control point *)


(**
   Build a control point for a quadratic curve for passing throuht 3 points.
   taken from https://xuhehuan.com/2608.html


   also look to https://pomax.github.io/bezierinfo/#pointcurves
*)
let three_points_quadratic
  : Gg.v2 -> Gg.v2 -> Gg.v2 -> quadratic
  = fun p0 c1 p1 ->

    let open Gg.V2 in

    let vect_1 = p0 - c1
    and vect_2 = p1 - c1 in
    let norm1 = norm vect_1
    and norm2 = norm vect_2 in
    let v = (Float.sqrt (norm1 *. norm2)) /. 2. in

    let ctrl = c1 - v * (( vect_1 / norm1) + (vect_2 / norm2)) in
    {p0; p1; ctrl}

(**

   Convert a cubic bezier curve into a quadratic one

*)
let quadratic_to_cubic
  : quadratic -> t
  = fun {p0; p1; ctrl} ->

    let coef = 2. /. 3. in

    let open Gg.V2 in
    { p0
    ; p1
    ; ctrl0 = mix p0 ctrl coef
    ; ctrl1 = mix p1 ctrl coef }



let abc_ratio
  : int -> float -> float
  = fun n t ->
    let n' = Float.of_int n in
    let bottom = (Float.pow t n') +. (Float.pow  (1. -. t) n') in
    let top = bottom -. 1. in
    Float.abs (top /. bottom)

let half_cubic_ratio = abc_ratio 3 0.5

exception Not_found

(**

   https://pomax.github.io/bezierinfo/#pointcurves

*)
let three_points_cubic
  : float -> Gg.v2 -> Gg.v2 -> Gg.v2 -> t
  = fun f p0 p1 p2 ->

    let open Gg.V2 in

    let c = half ( p0 + p2) in
    let a = p1 + ((p1 - c) / half_cubic_ratio) in

    let vect1_0 = p1 - p0 in
    let vect2_0 = p2 - p0 in

    let d1 = norm vect1_0
    and d2 = norm (p2 - p1) in
    let t = d1 /. (d1 +. d2) in

    let angle_1_0 = angle vect1_0
    and angle_2_0 = angle vect2_0 in

    (* get our e1-e2 distances *)
    let angle = mod_float
        (Gg.Float.two_pi
         +. angle_2_0
         -. angle_1_0)
        Gg.Float.two_pi in

    let distance = (norm vect2_0) *. f in

    let bc =
      if angle < 0. || angle > Gg.Float.pi then
        Float.(neg distance)
      else
        distance in
    let de1  = t *. bc
    and de2 = (1. -. t) *. bc in

    (* get the circle-aligned slope as normalized dx/dy *)
    let center = Utils.center p0 p1 p2 in
    match center with
    | None -> raise Not_found
    | Some center ->
      let t' = p1 - center in
      let tangent0 = v
          ((x p1) -. (y t'))
          ((y p1) +. (x t'))
      and tangent1 = v
          ((x p1) +. (y t'))
          ((y p1) -. (x t')) in

      let d = unit (tangent1 - tangent0) in

      (* then set up an e1 and e2 parallel to the baseline *)
      let e1 = p1 + de1 * d
      and e2 = p1 - de2 * d in

      (* then use those e1/e2 to derive the new hull coordinates *)
      let v1 = a + (e1 - a) / (1. -. t)
      and v2 = a + (e2 - a) / t in

      let ctrl0 = p0 + (v1 - p0) / t
      and ctrl1 = p2 + (v2 -p2) / (1. -. t) in

      {p0; p1 = p2; ctrl0; ctrl1}

(** Split a bezier curve in two at a given position *)
let slice
  : float -> t -> t * t
  = fun t {p0; p1; ctrl0; ctrl1} ->

    let mix p1 p2 = Gg.V2.mix p1 p2 t in

    let p12 = mix p0 ctrl0
    and p23 = mix ctrl0 ctrl1
    and p34 = mix ctrl1 p1 in

    let p123 = mix p12 p23
    and p234 = mix p23 p34 in

    let p1234 = mix p123 p234 in

    ( { p0
      ; ctrl0 = p12
      ; ctrl1 = p123
      ; p1 = p1234 }
    , { p0 = p1234
      ; ctrl0 = p234
      ; ctrl1 = p34
      ; p1 } )


let get_closest_point
  : Gg.v2 -> t -> float * Gg.v2
  = fun point t ->

    let rec f min max t =

      (* First devide the curve in two *)
      let seq_0, seq_1 = slice 0.5 t in
      let avg = (min +. max) /. 2. in

      let p0 = t.p0
      and p1 = t.p1
      and p01 = seq_0.p1 in (* seq_0.p1 = seq_1.p0 *)

      let open Gg.V2 in
      let center0 = mix p0 p01 0.5
      and center1 = mix p01 p1 0.5 in

      if Tools.Utils.equal_point 0.001 p0 p1 then
        avg, p01
      else if (norm (point - center0)) < (norm (point - center1)) then
        f min avg seq_0
      else
        f avg max seq_1

    in f 0. 1. t

let reverse
  : t -> t
  = fun bezier ->
    {
      p0 = bezier.p1
    ; p1 = bezier.p0
    ; ctrl0 = bezier.ctrl1
    ; ctrl1 = bezier.ctrl0 }

(**

   see https://github.com/Pomax/BezierInfo-2/blob/master/docs/js/graphics-element/lib/bezierjs/bezier.js#L504

   let root
   : t -> 'a
   = fun bezier ->

    let accept
      : float -> bool
      = fun t ->
        0. <= t && t <= 1. in

    let cuberoot v =
      if v < 0.  then
        Float.pow (Float.neg v) ( 1. /. 3.)
        |> Float.neg
      else Float.pow v (1. /. 3.) in




    failwith "Non implemented"
*)