1// Copyright 2021 the Kurbo Authors
2// SPDX-License-Identifier: Apache-2.0 OR MIT
3
4//! Minimum distance between two Bézier curves
5//!
6//! This implements the algorithm in "Computing the minimum distance between
7//! two Bézier curves", Chen et al., *Journal of Computational and Applied
8//! Mathematics* 229(2009), 294-301
9
10use crate::Vec2;
11use core::cmp::Ordering;
12
13#[cfg(not(feature = "std"))]
14use crate::common::FloatFuncs;
15
16pub(crate) fn min_dist_param(
17 bez1: &[Vec2],
18 bez2: &[Vec2],
19 u: (f64, f64),
20 v: (f64, f64),
21 epsilon: f64,
22 best_alpha: Option<f64>,
23) -> (f64, f64, f64) {
24 assert!(!bez1.is_empty() && !bez2.is_empty());
25 let n = bez1.len() - 1;
26 let m = bez2.len() - 1;
27 let (umin, umax) = u;
28 let (vmin, vmax) = v;
29 let umid = (umin + umax) / 2.0;
30 let vmid = (vmin + vmax) / 2.0;
31 let svalues: [(f64, f64, f64); 4] = [
32 (S(umin, vmin, bez1, bez2), umin, vmin),
33 (S(umin, vmax, bez1, bez2), umin, vmax),
34 (S(umax, vmin, bez1, bez2), umax, vmin),
35 (S(umax, vmax, bez1, bez2), umax, vmax),
36 ];
37 let alpha: f64 = svalues.iter().map(|(a, _, _)| *a).reduce(f64::min).unwrap();
38 if let Some(best) = best_alpha {
39 if alpha > best {
40 return (alpha, umid, vmid);
41 }
42 }
43
44 if (umax - umin).abs() < epsilon || (vmax - vmin).abs() < epsilon {
45 return (alpha, umid, vmid);
46 }
47
48 // Property one: D(r>k) > alpha
49 let mut is_outside = true;
50 let mut min_drk = None;
51 let mut min_ij = None;
52 for r in 0..(2 * n) {
53 for k in 0..(2 * m) {
54 let d_rk = D_rk(r, k, bez1, bez2);
55 if d_rk < alpha {
56 is_outside = false;
57 }
58 if min_drk.is_none() || Some(d_rk) < min_drk {
59 min_drk = Some(d_rk);
60 min_ij = Some((r, k));
61 }
62 }
63 }
64 if is_outside {
65 return (alpha, umid, vmid);
66 }
67
68 // Property two: boundary check
69 let mut at_boundary0on_bez1 = true;
70 let mut at_boundary1on_bez1 = true;
71 let mut at_boundary0on_bez2 = true;
72 let mut at_boundary1on_bez2 = true;
73 for i in 0..2 * n {
74 for j in 0..2 * m {
75 let dij = D_rk(i, j, bez1, bez2);
76 let dkj = D_rk(0, j, bez1, bez2);
77 if dij < dkj {
78 at_boundary0on_bez1 = false
79 }
80 let dkj = D_rk(2 * n, j, bez1, bez2);
81 if dij < dkj {
82 at_boundary1on_bez1 = false
83 }
84 let dkj = D_rk(i, 0, bez1, bez2);
85 if dij < dkj {
86 at_boundary0on_bez2 = false
87 }
88 let dkj = D_rk(i, 2 * n, bez1, bez2);
89 if dij < dkj {
90 at_boundary1on_bez2 = false
91 }
92 }
93 }
94 if at_boundary0on_bez1 && at_boundary0on_bez2 {
95 return svalues[0];
96 }
97 if at_boundary0on_bez1 && at_boundary1on_bez2 {
98 return svalues[1];
99 }
100 if at_boundary1on_bez1 && at_boundary0on_bez2 {
101 return svalues[2];
102 }
103 if at_boundary1on_bez1 && at_boundary1on_bez2 {
104 return svalues[3];
105 }
106
107 let (min_i, min_j) = min_ij.unwrap();
108 let new_umid = umin + (umax - umin) * (min_i as f64 / (2 * n) as f64);
109 let new_vmid = vmin + (vmax - vmin) * (min_j as f64 / (2 * m) as f64);
110
111 // Subdivide
112 let results = [
113 min_dist_param(
114 bez1,
115 bez2,
116 (umin, new_umid),
117 (vmin, new_vmid),
118 epsilon,
119 Some(alpha),
120 ),
121 min_dist_param(
122 bez1,
123 bez2,
124 (umin, new_umid),
125 (new_vmid, vmax),
126 epsilon,
127 Some(alpha),
128 ),
129 min_dist_param(
130 bez1,
131 bez2,
132 (new_umid, umax),
133 (vmin, new_vmid),
134 epsilon,
135 Some(alpha),
136 ),
137 min_dist_param(
138 bez1,
139 bez2,
140 (new_umid, umax),
141 (new_vmid, vmax),
142 epsilon,
143 Some(alpha),
144 ),
145 ];
146
147 *results
148 .iter()
149 .min_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(Ordering::Less))
150 .unwrap()
151}
152
153// $ S(u,v)=\sum_{r=0}^{2n} \sum _{k=0}^{2m} D_{r,k} B_{2n}^r(u) B_{2m}^k(v) $
154
155#[allow(non_snake_case)]
156fn S(u: f64, v: f64, bez1: &[Vec2], bez2: &[Vec2]) -> f64 {
157 let n: usize = bez1.len() - 1;
158 let m: usize = bez2.len() - 1;
159 let mut summand: f64 = 0.0;
160 for r: usize in 0..=2 * n {
161 for k: usize in 0..=2 * m {
162 summand +=
163 D_rk(r, k, bez1, bez2) * basis_function(n:2 * n, i:r, u) * basis_function(n:2 * m, i:k, u:v)
164 }
165 }
166 summand
167}
168
169// $ C_{r,k} = ( \sum_{i=\theta}^\upsilon P_i C_n^i C_n^{r-i} / C_{2n}^r ) \cdot (\sum_{j=\sigma}^\varsigma Q_j C_m^j C_m^{k-j} / C_{2m}^k ) $
170#[allow(non_snake_case)]
171fn C_rk(r: usize, k: usize, bez1: &[Vec2], bez2: &[Vec2]) -> f64 {
172 let n: usize = bez1.len() - 1;
173 let upsilon: usize = r.min(n);
174 let theta: usize = r - n.min(r);
175 let mut left: Vec2 = Vec2::ZERO;
176 for (i: usize, &item: Vec2) in bez1.iter().enumerate().take(upsilon + 1).skip(theta) {
177 left += item * (choose(n, k:i) * choose(n, k:r - i)) as f64 / (choose(n:2 * n, k:r) as f64);
178 }
179
180 let m: usize = bez2.len() - 1;
181 let varsigma: usize = k.min(m);
182 let sigma: usize = k - m.min(k);
183 let mut right: Vec2 = Vec2::ZERO;
184 for (j: usize, &item: Vec2) in bez2.iter().enumerate().take(varsigma + 1).skip(sigma) {
185 right += item * (choose(n:m, k:j) * choose(n:m, k:k - j)) as f64 / (choose(n:2 * m, k) as f64);
186 }
187 left.dot(right)
188}
189
190// $ A_r=\sum _{i=\theta} ^\upsilon (P_i \cdot P_{r-i}) C_n^i C_n^{r-i} / C_{2n}^r $
191// ($ B_k $ is just the same as $ A_r $ but for the other curve.)
192
193#[allow(non_snake_case)]
194fn A_r(r: usize, p: &[Vec2]) -> f64 {
195 let n: usize = p.len() - 1;
196 let upsilon: usize = r.min(n);
197 let theta: usize = r - n.min(r);
198 (theta..=upsilon)
199 .map(|i: usize| {
200 let dot: f64 = p[i].dot(p[r - i]); // These are bounds checked by the sum limits
201 let factor: f64 = (choose(n, k:i) * choose(n, k:r - i)) as f64 / (choose(n:2 * n, k:r) as f64);
202 dot * factor
203 })
204 .sum()
205}
206
207#[allow(non_snake_case)]
208fn D_rk(r: usize, k: usize, bez1: &[Vec2], bez2: &[Vec2]) -> f64 {
209 // In the paper, B_k is used for the second factor, but it's the same thing
210 A_r(r, p:bez1) + A_r(r:k, p:bez2) - 2.0 * C_rk(r, k, bez1, bez2)
211}
212
213// Bezier basis function
214fn basis_function(n: usize, i: usize, u: f64) -> f64 {
215 choose(n, k:i) as f64 * (1.0 - u).powi((n - i) as i32) * u.powi(i as i32)
216}
217
218// Binomial co-efficient, but returning zeros for values outside of domain
219fn choose(n: usize, k: usize) -> u32 {
220 let mut n: usize = n;
221 if k > n {
222 return 0;
223 }
224 let mut p: usize = 1;
225 for i: usize in 1..=(n - k) {
226 p *= n;
227 p /= i;
228 n -= 1;
229 }
230 p as u32
231}
232
233#[cfg(test)]
234mod tests {
235 use crate::mindist::A_r;
236 use crate::mindist::{choose, D_rk};
237 use crate::{CubicBez, Line, PathSeg, Vec2};
238
239 #[test]
240 fn test_choose() {
241 assert_eq!(choose(6, 0), 1);
242 assert_eq!(choose(6, 1), 6);
243 assert_eq!(choose(6, 2), 15);
244 }
245
246 #[test]
247 fn test_d_rk() {
248 let bez1 = vec![
249 Vec2::new(129.0, 139.0),
250 Vec2::new(190.0, 139.0),
251 Vec2::new(201.0, 364.0),
252 Vec2::new(90.0, 364.0),
253 ];
254 let bez2 = vec![
255 Vec2::new(309.0, 159.0),
256 Vec2::new(178.0, 159.0),
257 Vec2::new(215.0, 408.0),
258 Vec2::new(309.0, 408.0),
259 ];
260 let b = A_r(1, &bez2);
261 assert!((b - 80283.0).abs() < 0.005, "B_1(Q)={b:?}");
262 let d = D_rk(0, 1, &bez1, &bez2);
263 assert!((d - 9220.0).abs() < 0.005, "D={d:?}");
264 }
265
266 #[test]
267 fn test_mindist() {
268 let bez1 = PathSeg::Cubic(CubicBez::new(
269 (129.0, 139.0),
270 (190.0, 139.0),
271 (201.0, 364.0),
272 (90.0, 364.0),
273 ));
274 let bez2 = PathSeg::Cubic(CubicBez::new(
275 (309.0, 159.0),
276 (178.0, 159.0),
277 (215.0, 408.0),
278 (309.0, 408.0),
279 ));
280 let mindist = bez1.min_dist(bez2, 0.001);
281 assert!((mindist.distance - 50.9966).abs() < 0.5);
282 }
283
284 #[test]
285 fn test_overflow() {
286 let bez1 = PathSeg::Cubic(CubicBez::new(
287 (232.0, 126.0),
288 (134.0, 126.0),
289 (139.0, 232.0),
290 (141.0, 301.0),
291 ));
292 let bez2 = PathSeg::Line(Line::new((359.0, 416.0), (367.0, 755.0)));
293 let mindist = bez1.min_dist(bez2, 0.001);
294 assert!((mindist.distance - 246.4731222669117).abs() < 0.5);
295 }
296
297 #[test]
298 fn test_out_of_order() {
299 let bez1 = PathSeg::Cubic(CubicBez::new(
300 (287.0, 182.0),
301 (346.0, 277.0),
302 (356.0, 299.0),
303 (359.0, 416.0),
304 ));
305 let bez2 = PathSeg::Line(Line::new((141.0, 301.0), (152.0, 709.0)));
306 let mindist1 = bez1.min_dist(bez2, 0.5);
307 let mindist2 = bez2.min_dist(bez1, 0.5);
308 assert!((mindist1.distance - mindist2.distance).abs() < 0.5);
309 }
310}
311