|
5 | 5 | //! well as some code that was copied from kurbo, which is needed to reimplement the |
6 | 6 | //! full `flatten` method. |
7 | 7 |
|
| 8 | +use crate::flatten::TOL_2; |
| 9 | +#[cfg(not(feature = "std"))] |
| 10 | +use crate::kurbo::common::FloatFuncs as _; |
8 | 11 | use crate::kurbo::{CubicBez, ParamCurve, PathEl, Point, QuadBez}; |
9 | 12 | use alloc::vec; |
10 | 13 | use alloc::vec::Vec; |
11 | 14 | use bytemuck::{Pod, Zeroable}; |
12 | 15 | use fearless_simd::*; |
13 | 16 |
|
14 | | -#[cfg(not(feature = "std"))] |
15 | | -use crate::kurbo::common::FloatFuncs as _; |
16 | | - |
17 | 17 | // Unlike kurbo, which takes a closure with a callback for outputting the lines, we use a trait |
18 | 18 | // instead. The reason is that this way the callback can be inlined, which is not possible with |
19 | 19 | // a closure and turned out to have a noticeable overhead. |
@@ -49,38 +49,84 @@ pub(crate) fn flatten<S: Simd>( |
49 | 49 | } |
50 | 50 | PathEl::QuadTo(p1, p2) => { |
51 | 51 | if let Some(p0) = last_pt { |
52 | | - let q = QuadBez::new(p0, p1, p2); |
53 | | - let params = q.estimate_subdiv(sqrt_tol); |
54 | | - let n = ((0.5 * params.val / sqrt_tol).ceil() as usize).max(1); |
55 | | - let step = 1.0 / (n as f64); |
56 | | - for i in 1..n { |
57 | | - let u = (i as f64) * step; |
58 | | - let t = q.determine_subdiv_t(¶ms, u); |
59 | | - let p = q.eval(t); |
60 | | - callback.callback(PathEl::LineTo(p)); |
| 52 | + // An upper bound on the shortest distance of any point on the quadratic Bezier |
| 53 | + // curve to the line segment [p0, p2] is 1/2 of the maximum of the |
| 54 | + // endpoint-to-control-point distances. |
| 55 | + // |
| 56 | + // The derivation is similar to that for the cubic Bezier (see below). In |
| 57 | + // short: |
| 58 | + // |
| 59 | + // q(t) = B0(t) p0 + B1(t) p1 + B2(t) p2 |
| 60 | + // dist(q(t), [p0, p1]) <= B1(t) dist(p1, [p0, p1]) |
| 61 | + // = 2 (1-t)t dist(p1, [p0, p1]). |
| 62 | + // |
| 63 | + // The maximum occurs at t=1/2, hence |
| 64 | + // max(dist(q(t), [p0, p1] <= 1/2 dist(p1, [p0, p1])). |
| 65 | + // |
| 66 | + // A cheap upper bound for dist(p1, [p0, p1]) is max(dist(p1, p0), dist(p1, p2)). |
| 67 | + // |
| 68 | + // The following takes the square to elide the square root of the Euclidean |
| 69 | + // distance. |
| 70 | + if f64::max((p1 - p0).hypot2(), (p1 - p2).hypot2()) <= 4. * TOL_2 { |
| 71 | + callback.callback(PathEl::LineTo(p2)); |
| 72 | + } else { |
| 73 | + let q = QuadBez::new(p0, p1, p2); |
| 74 | + let params = q.estimate_subdiv(sqrt_tol); |
| 75 | + let n = ((0.5 * params.val / sqrt_tol).ceil() as usize).max(1); |
| 76 | + let step = 1.0 / (n as f64); |
| 77 | + for i in 1..n { |
| 78 | + let u = (i as f64) * step; |
| 79 | + let t = q.determine_subdiv_t(¶ms, u); |
| 80 | + let p = q.eval(t); |
| 81 | + callback.callback(PathEl::LineTo(p)); |
| 82 | + } |
| 83 | + callback.callback(PathEl::LineTo(p2)); |
61 | 84 | } |
62 | | - callback.callback(PathEl::LineTo(p2)); |
63 | 85 | } |
64 | 86 | last_pt = Some(p2); |
65 | 87 | } |
66 | 88 | PathEl::CurveTo(p1, p2, p3) => { |
67 | 89 | if let Some(p0) = last_pt { |
68 | | - let c = CubicBez::new(p0, p1, p2, p3); |
69 | | - let max = simd.vectorize( |
70 | | - #[inline(always)] |
71 | | - || { |
72 | | - flatten_cubic_simd( |
73 | | - simd, |
74 | | - c, |
75 | | - flatten_ctx, |
76 | | - tolerance as f32, |
77 | | - &mut flattened_cubics, |
78 | | - ) |
79 | | - }, |
80 | | - ); |
81 | | - |
82 | | - for p in &flattened_cubics[1..max] { |
83 | | - callback.callback(PathEl::LineTo(Point::new(p.x as f64, p.y as f64))); |
| 90 | + // An upper bound on the shortest distance of any point on the cubic Bezier |
| 91 | + // curve to the line segment [p0, p3] is 3/4 of the maximum of the |
| 92 | + // endpoint-to-control-point distances. |
| 93 | + // |
| 94 | + // With Bernstein weights Bi(t), we have |
| 95 | + // c(t) = B0(t) p0 + B1(t) p1 + B2(t) p2 + B3(t) p3 |
| 96 | + // with t from 0 to 1 (inclusive). |
| 97 | + // |
| 98 | + // Through convexivity of the Euclidean distance function and the line segment, |
| 99 | + // we have |
| 100 | + // dist(c(t), [p0, p3]) <= B1(t) dist(p1, [p0, p3]) + B2(t) dist(p2, [p0, p3]) |
| 101 | + // <= B1(t) ||p1-p0|| + B2(t) ||p2-p3|| |
| 102 | + // <= (B1(t) + B2(t)) max(||p1-p0||, ||p2-p3|||) |
| 103 | + // = 3 ((1-t)t^2 + (1-t)^2t) max(||p1-p0||, ||p2-p3||). |
| 104 | + // |
| 105 | + // The inner polynomial has its maximum of 1/4 at t=1/2, hence |
| 106 | + // max(dist(c(t), [p0, p3])) <= 3/4 max(||p1-p0||, ||p2-p3||). |
| 107 | + // |
| 108 | + // The following takes the square to elide the square root of the Euclidean |
| 109 | + // distance. |
| 110 | + if f64::max((p0 - p1).hypot2(), (p3 - p2).hypot2()) <= 16. / 9. * TOL_2 { |
| 111 | + callback.callback(PathEl::LineTo(p3)); |
| 112 | + } else { |
| 113 | + let c = CubicBez::new(p0, p1, p2, p3); |
| 114 | + let max = simd.vectorize( |
| 115 | + #[inline(always)] |
| 116 | + || { |
| 117 | + flatten_cubic_simd( |
| 118 | + simd, |
| 119 | + c, |
| 120 | + flatten_ctx, |
| 121 | + tolerance as f32, |
| 122 | + &mut flattened_cubics, |
| 123 | + ) |
| 124 | + }, |
| 125 | + ); |
| 126 | + |
| 127 | + for p in &flattened_cubics[1..max] { |
| 128 | + callback.callback(PathEl::LineTo(Point::new(p.x as f64, p.y as f64))); |
| 129 | + } |
84 | 130 | } |
85 | 131 | } |
86 | 132 | last_pt = Some(p3); |
|
0 commit comments