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
type point = float * float
let cubic ((x0, y0) : point) ((x1, y1) : point) ((x2, y2) : point) ((x3, y3) : point) (t : float) : point =
let u = 1. -. t in
let f a b c d = (u *. u *. u *. a) +. (3. *. u *. u *. t *. b) +. (3. *. u *. t *. t *. c) +. (t *. t *. t *. d) in
(f x0 x1 x2 x3, f y0 y1 y2 y3)
let quadratic ((x0, y0) : point) ((x1, y1) : point) ((x2, y2) : point) (t : float) : point =
let u = 1. -. t in
let f a b c = (u *. u *. a) +. (2. *. u *. t *. b) +. (t *. t *. c) in
(f x0 x1 x2, f y0 y1 y2)
let catmull_rom ((x0, y0) : point) ((x1, y1) : point) ((x2, y2) : point) ((x3, y3) : point) (t : float) : point =
let f a b c d = 0.5 *. ((2. *. b) +. ((c -. a) *. t) +. (((2. *. a) -. (5. *. b) +. (4. *. c) -. d) *. t *. t) +. (((3. *. b) -. a -. (3. *. c) +. d) *. t *. t *. t)) in
(f x0 x1 x2 x3, f y0 y1 y2 y3)
let cubic_of_catmull_rom (p0 : point) (p1 : point) (p2 : point) (p3 : point) : point * point * point * point =
(p1, Vec2.add p1 (Vec2.scale (1. /. 6.) (Vec2.sub p2 p0)), Vec2.sub p2 (Vec2.scale (1. /. 6.) (Vec2.sub p3 p1)), p2)
let distance_to_chord (a : point) (b : point) (p : point) : float =
let chord = Vec2.sub b a in
let n = Vec2.length chord in
if n = 0. then Vec2.length (Vec2.sub p a) else Float.abs (Vec2.cross chord (Vec2.sub p a)) /. n
let flatten ?(tolerance = 0.1) (p0 : point) (p1 : point) (p2 : point) (p3 : point) : point list =
let flat_enough a b c d = distance_to_chord a d b <= tolerance && distance_to_chord a d c <= tolerance in
let rec go depth a b c d acc =
if depth = 0 || flat_enough a b c d then d :: acc
else
let mid = Vec2.scale 0.5 (Vec2.add b c) in
let lb = Vec2.scale 0.5 (Vec2.add a b) in
let rc = Vec2.scale 0.5 (Vec2.add c d) in
let lc = Vec2.scale 0.5 (Vec2.add lb mid) in
let rb = Vec2.scale 0.5 (Vec2.add mid rc) in
let m = Vec2.scale 0.5 (Vec2.add lc rb) in
go (depth - 1) a lb lc m (go (depth - 1) m rb rc d acc)
in
p0 :: go 16 p0 p1 p2 p3 []
let through ?(steps = 16) (points : point list) : point list =
let a = Array.of_list points in
let n = Array.length a in
let get i = a.(max 0 (min (n - 1) i)) in
List.concat (List.init (n - 1) (fun i -> List.init steps (fun k -> catmull_rom (get (i - 1)) (get i) (get (i + 1)) (get (i + 2)) (float_of_int k /. float_of_int steps))))
@ [ a.(n - 1) ]
type t = { pts : point array; lengths : float array }
let measure (polyline : point list) : t =
let pts = Array.of_list polyline in
let lengths = Array.make (Array.length pts) 0. in
for i = 1 to Array.length pts - 1 do
let (x0, y0), (x1, y1) = (pts.(i - 1), pts.(i)) in
lengths.(i) <- lengths.(i - 1) +. Float.hypot (x1 -. x0) (y1 -. y0)
done;
{ pts; lengths }
let length (c : t) : float = c.lengths.(Array.length c.lengths - 1)
let at (c : t) (s : float) : point * float =
let n = Array.length c.pts in
if n < 2 then (c.pts.(0), 0.)
else
let rec find i = if i < n - 1 && c.lengths.(i) < s then find (i + 1) else i in
let i = max 1 (find 1) in
let (x0, y0), (x1, y1) = (c.pts.(i - 1), c.pts.(i)) in
let seg = c.lengths.(i) -. c.lengths.(i - 1) in
let f = if seg > 0. then Float.max 0. (Float.min 1. ((s -. c.lengths.(i - 1)) /. seg)) else 1. in
((x0 +. (f *. (x1 -. x0)), y0 +. (f *. (y1 -. y0))), Float.atan2 (y1 -. y0) (x1 -. x0))