summaryrefslogtreecommitdiff
path: root/shapes/bezier.ml
blob: f5f288c4fbc01d86fd79ba197cbcbcc2c2937e71 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
(**

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