summaryrefslogtreecommitdiff
path: root/curves
diff options
context:
space:
mode:
authorSébastien Dailly <sebastien@chimrod.com>2020-12-16 14:39:42 +0100
committerSébastien Dailly <sebastien@chimrod.com>2020-12-16 14:39:42 +0100
commit4f262d6540281487f79870aff589ca92f5d2f6c6 (patch)
tree940e59d943715366d1aa72bb93f248dcd65ab992 /curves
Initial commit
Diffstat (limited to 'curves')
-rwxr-xr-xcurves/bezier.ml197
-rwxr-xr-xcurves/bezier.mli40
-rwxr-xr-xcurves/bspline.ml149
-rwxr-xr-xcurves/bspline.mli24
-rwxr-xr-xcurves/dd_splines.pdfbin0 -> 184248 bytes
-rwxr-xr-xcurves/dune7
6 files changed, 417 insertions, 0 deletions
diff --git a/curves/bezier.ml b/curves/bezier.ml
new file mode 100755
index 0000000..3dedc70
--- /dev/null
+++ b/curves/bezier.ml
@@ -0,0 +1,197 @@
+(**
+
+ 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 rec get_closest_point
+ : Gg.v2 -> t -> Gg.v2
+ = fun point t ->
+ (* First devide the curve in two *)
+ let seq_0, seq_1 = slice 0.5 t 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
+ p01
+ else if (norm (point - center0)) < (norm (point - center1)) then
+ get_closest_point point seq_0
+ else
+ get_closest_point point seq_1
+
+let reverse
+ : t -> t
+ = fun bezier ->
+ {
+ p0 = bezier.p1
+ ; p1 = bezier.p0
+ ; ctrl0 = bezier.ctrl1
+ ; ctrl1 = bezier.ctrl0 }
diff --git a/curves/bezier.mli b/curves/bezier.mli
new file mode 100755
index 0000000..e90163c
--- /dev/null
+++ b/curves/bezier.mli
@@ -0,0 +1,40 @@
+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 *)
+
+type quadratic
+
+(**
+ 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
+*)
+val three_points_quadratic
+ : Gg.v2 -> Gg.v2 -> Gg.v2 -> quadratic
+
+(**
+ Create a curve from three points.
+
+ This is an implementation for
+ https://pomax.github.io/bezierinfo/#pointcurves
+
+*)
+val three_points_cubic
+ : float -> Gg.v2 -> Gg.v2 -> Gg.v2 -> t
+
+val quadratic_to_cubic
+ : quadratic -> t
+
+(** Split a bezier curve in two at a given position *)
+val slice
+ : float -> t -> t * t
+
+(** Return the closest point to the curve by approximation *)
+val get_closest_point
+ : Gg.v2 -> t -> Gg.v2
+
+val reverse: t -> t
diff --git a/curves/bspline.ml b/curves/bspline.ml
new file mode 100755
index 0000000..bb60227
--- /dev/null
+++ b/curves/bspline.ml
@@ -0,0 +1,149 @@
+open StdLabels
+
+type err = [`InvalidPath ]
+
+module M = Matrix.MakeMatrix (struct
+
+ type t = Float.t
+
+ let compare a b =
+
+ let v = Float.compare a b in
+ if v = 0 then Matrix.Order.Equal
+ else if v > 0 then Matrix.Order.Greater
+ else Matrix.Order.Less
+
+ let zero = Float.zero
+ let one = Float.one
+
+ let divide = (/.)
+ let multiply = ( *. )
+ let add = (+.)
+ let subtract = (-.)
+ exception NonElt
+
+
+ end)
+
+type t = Gg.v2 list
+
+let from_points
+ : Gg.v2 array -> (Gg.v2 array, [> `InvalidPath]) Result.t
+ = fun points ->
+
+ let n = (Array.length points - 2) in
+
+ if n <= 1 then
+ Result.error `InvalidPath
+ else
+
+ (* Create the initial matrix.
+
+ The matrix is augmented with two additionals columns, which will be
+ populated with the points from the path.
+ *)
+ let arr = Array.init n ~f:(fun line ->
+ Array.init (n +2) ~f:(fun row ->
+ match row - line with
+ | (-1) -> 1.
+ | 0 -> 4.
+ | 1 -> 1.
+ | _ -> 0.
+ )
+ ) in
+ let matrix = M.from_array arr in
+
+ (* Add the points from the augmented matrix *)
+ let points_array = points in
+ for line = 0 to (n -1) do
+
+ let point =
+ if line = 0 then
+ let p0 = points_array.(0)
+ and p1 = points_array.(1) in
+ Gg.V2.(6. * p1 - p0)
+ else if (line + 1) = n then
+ let p_n_2 = points_array.(n)
+ and p_n_1 = points_array.(n + 1) in
+ Gg.V2.(6. * p_n_2 - p_n_1)
+ else
+ let n' = line + 1 in
+ Gg.V2.(6. * points_array.(n'))
+ in
+ let x = (Gg.V2.x point)
+ and y = (Gg.V2.y point) in
+
+
+ M.set_elt matrix (line + 1, n + 1) x;
+ M.set_elt matrix (line + 1, n + 2) y;
+ done;
+
+ (* Resolve the matrix *)
+ let res' = M.row_reduce matrix in
+
+ (* Extract the result as points *)
+ let _, col_x = M.get_column res' (n + 1)
+ and _, col_y = M.get_column res' (n + 2) in
+
+ (* Build the result *)
+ let res = Array.make (n + 2) (Array.get points_array (n + 1) ) in
+ for i = 1 to n do
+ let point = Gg.V2.v col_x.(i - 1) col_y.(i - 1) in
+ Array.set res i point;
+ done;
+ Array.set res 0 (Array.get points_array 0);
+ Result.ok res
+
+let (let*) = Result.bind
+
+(** Build a continue curve from path
+
+ see https://www.math.ucla.edu/~baker/149.1.02w/handouts/dd_splines.pdf
+*)
+let to_bezier
+ : ?connexion0:Gg.v2 -> ?connexion1:Gg.v2 -> t -> (Bezier.t array, [> `InvalidPath]) Result.t
+ = fun ?connexion0 ?connexion1 points ->
+
+ let points' = match connexion0 with
+ | None -> points
+ | Some pt -> pt::points in
+
+ let arr_points = match connexion1 with
+ | None -> Array.of_list points'
+ | Some pt ->
+ let arr = Array.make (1 + (List.length points')) pt in
+ List.iteri points'
+ ~f:(fun i value -> Array.set arr i value);
+ arr in
+
+ let* bspline_points = from_points arr_points in
+
+ let start = match connexion0 with
+ | None -> 1
+ | Some _ -> 2
+ and end_ = match connexion1 with
+ | None -> (Array.length bspline_points) - 1
+ | Some _ -> (Array.length bspline_points) - 2 in
+
+ let result = Array.init (end_ - start + 1) ~f:(fun i ->
+
+ let i = i + start in
+
+ let prev_b = Array.get bspline_points (i - 1)
+ and bpoint = Array.get bspline_points i
+ and prev_p = Array.get arr_points (i - 1)
+ and point = Array.get arr_points i in
+ let ctrl0 = Gg.V2.(mix prev_b bpoint (1. /. 3.))
+ and ctrl1 = Gg.V2.(mix prev_b bpoint (2. /. 3.)) in
+
+ let bezier =
+ { Bezier.p0 = prev_p
+ ; Bezier.p1 = point
+ ; Bezier.ctrl0
+ ; Bezier.ctrl1 } in
+
+ bezier
+
+ ) in
+ Result.Ok result
+
diff --git a/curves/bspline.mli b/curves/bspline.mli
new file mode 100755
index 0000000..074658d
--- /dev/null
+++ b/curves/bspline.mli
@@ -0,0 +1,24 @@
+type t
+
+type err =
+ [ `InvalidPath (* Too few points in the path for building the curve *)
+ ]
+
+(** Convert a list of points into a beziers curves.
+
+ At least 4 points are required for building the path.
+
+ [to_bezier ~connexion points] create a list of beziers segments joining all
+ the points together.
+
+ [connexion0] add a virtual point in the begining for helping to get the
+ appropriate tangent when connecting path together
+
+ [connexion1] does the same at the end
+
+*)
+val to_bezier
+ : ?connexion0:Gg.v2
+ -> ?connexion1:Gg.v2
+ -> Gg.v2 list
+ -> (Bezier.t array, [> err]) Result.t
diff --git a/curves/dd_splines.pdf b/curves/dd_splines.pdf
new file mode 100755
index 0000000..2618162
--- /dev/null
+++ b/curves/dd_splines.pdf
Binary files differ
diff --git a/curves/dune b/curves/dune
new file mode 100755
index 0000000..34c87fb
--- /dev/null
+++ b/curves/dune
@@ -0,0 +1,7 @@
+(library
+ (name curves)
+ (libraries
+ tools
+ matrix
+ )
+ )