1 //! Quadratic Bézier segments.
2
3 use std::ops::{Mul, Range};
4
5 use arrayvec::ArrayVec;
6
7 use crate::common::solve_cubic;
8 use crate::MAX_EXTREMA;
9 use crate::{
10 Affine, CubicBez, Line, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveCurvature,
11 ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, Point,
12 };
13
14 /// A single quadratic Bézier segment.
15 #[derive(Clone, Copy, Debug, PartialEq)]
16 #[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
17 #[allow(missing_docs)]
18 pub struct QuadBez {
19 pub p0: Point,
20 pub p1: Point,
21 pub p2: Point,
22 }
23
24 impl QuadBez {
25 /// Create a new quadratic Bézier segment.
26 #[inline]
new<V: Into<Point>>(p0: V, p1: V, p2: V) -> QuadBez27 pub fn new<V: Into<Point>>(p0: V, p1: V, p2: V) -> QuadBez {
28 QuadBez {
29 p0: p0.into(),
30 p1: p1.into(),
31 p2: p2.into(),
32 }
33 }
34
35 /// Raise the order by 1.
36 ///
37 /// Returns a cubic Bézier segment that exactly represents this quadratic.
38 #[inline]
raise(&self) -> CubicBez39 pub fn raise(&self) -> CubicBez {
40 CubicBez::new(
41 self.p0,
42 self.p0 + (2.0 / 3.0) * (self.p1 - self.p0),
43 self.p2 + (2.0 / 3.0) * (self.p1 - self.p2),
44 self.p2,
45 )
46 }
47
48 /// Estimate the number of subdivisions for flattening.
estimate_subdiv(&self, sqrt_tol: f64) -> FlattenParams49 pub(crate) fn estimate_subdiv(&self, sqrt_tol: f64) -> FlattenParams {
50 // Determine transformation to $y = x^2$ parabola.
51 let d01 = self.p1 - self.p0;
52 let d12 = self.p2 - self.p1;
53 let dd = d01 - d12;
54 let cross = (self.p2 - self.p0).cross(dd);
55 let x0 = d01.dot(dd) * cross.recip();
56 let x2 = d12.dot(dd) * cross.recip();
57 let scale = (cross / (dd.hypot() * (x2 - x0))).abs();
58
59 // Compute number of subdivisions needed.
60 let a0 = approx_parabola_integral(x0);
61 let a2 = approx_parabola_integral(x2);
62 let val = if scale.is_finite() {
63 let da = (a2 - a0).abs();
64 let sqrt_scale = scale.sqrt();
65 if x0.signum() == x2.signum() {
66 da * sqrt_scale
67 } else {
68 // Handle cusp case (segment contains curvature maximum)
69 let xmin = sqrt_tol / sqrt_scale;
70 sqrt_tol * da / approx_parabola_integral(xmin)
71 }
72 } else {
73 0.0
74 };
75 let u0 = approx_parabola_inv_integral(a0);
76 let u2 = approx_parabola_inv_integral(a2);
77 let uscale = (u2 - u0).recip();
78 FlattenParams {
79 a0,
80 a2,
81 u0,
82 uscale,
83 val,
84 }
85 }
86
87 // Maps a value from 0..1 to 0..1.
determine_subdiv_t(&self, params: &FlattenParams, x: f64) -> f6488 pub(crate) fn determine_subdiv_t(&self, params: &FlattenParams, x: f64) -> f64 {
89 let a = params.a0 + (params.a2 - params.a0) * x;
90 let u = approx_parabola_inv_integral(a);
91 (u - params.u0) * params.uscale
92 }
93 }
94
95 pub(crate) struct FlattenParams {
96 a0: f64,
97 a2: f64,
98 u0: f64,
99 uscale: f64,
100 /// The number of subdivisions * 2 * sqrt_tol.
101 pub(crate) val: f64,
102 }
103
104 /// An approximation to $\int (1 + 4x^2) ^ -0.25 dx$
105 ///
106 /// This is used for flattening curves.
approx_parabola_integral(x: f64) -> f64107 fn approx_parabola_integral(x: f64) -> f64 {
108 const D: f64 = 0.67;
109 x / (1.0 - D + (D.powi(4) + 0.25 * x * x).sqrt().sqrt())
110 }
111
112 /// An approximation to the inverse parabola integral.
approx_parabola_inv_integral(x: f64) -> f64113 fn approx_parabola_inv_integral(x: f64) -> f64 {
114 const B: f64 = 0.39;
115 x * (1.0 - B + (B * B + 0.25 * x * x).sqrt())
116 }
117
118 impl ParamCurve for QuadBez {
119 #[inline]
eval(&self, t: f64) -> Point120 fn eval(&self, t: f64) -> Point {
121 let mt = 1.0 - t;
122 (self.p0.to_vec2() * (mt * mt)
123 + (self.p1.to_vec2() * (mt * 2.0) + self.p2.to_vec2() * t) * t)
124 .to_point()
125 }
126
127 #[inline]
start(&self) -> Point128 fn start(&self) -> Point {
129 self.p0
130 }
131
132 #[inline]
end(&self) -> Point133 fn end(&self) -> Point {
134 self.p2
135 }
136
137 /// Subdivide into halves, using de Casteljau.
138 #[inline]
subdivide(&self) -> (QuadBez, QuadBez)139 fn subdivide(&self) -> (QuadBez, QuadBez) {
140 let pm = self.eval(0.5);
141 (
142 QuadBez::new(self.p0, self.p0.midpoint(self.p1), pm),
143 QuadBez::new(pm, self.p1.midpoint(self.p2), self.p2),
144 )
145 }
146
subsegment(&self, range: Range<f64>) -> QuadBez147 fn subsegment(&self, range: Range<f64>) -> QuadBez {
148 let (t0, t1) = (range.start, range.end);
149 let p0 = self.eval(t0);
150 let p2 = self.eval(t1);
151 let p1 = p0 + (self.p1 - self.p0).lerp(self.p2 - self.p1, t0) * (t1 - t0);
152 QuadBez { p0, p1, p2 }
153 }
154 }
155
156 impl ParamCurveDeriv for QuadBez {
157 type DerivResult = Line;
158
159 #[inline]
deriv(&self) -> Line160 fn deriv(&self) -> Line {
161 Line::new(
162 (2.0 * (self.p1.to_vec2() - self.p0.to_vec2())).to_point(),
163 (2.0 * (self.p2.to_vec2() - self.p1.to_vec2())).to_point(),
164 )
165 }
166 }
167
168 impl ParamCurveArclen for QuadBez {
169 /// Arclength of a quadratic Bézier segment.
170 ///
171 /// This computation is based on an analytical formula. Since that formula suffers
172 /// from numerical instability when the curve is very close to a straight line, we
173 /// detect that case and fall back to Legendre-Gauss quadrature.
174 ///
175 /// Accuracy should be better than 1e-13 over the entire range.
176 ///
177 /// Adapted from <http://www.malczak.linuxpl.com/blog/quadratic-bezier-curve-length/>
178 /// with permission.
arclen(&self, _accuracy: f64) -> f64179 fn arclen(&self, _accuracy: f64) -> f64 {
180 let d2 = self.p0.to_vec2() - 2.0 * self.p1.to_vec2() + self.p2.to_vec2();
181 let a = d2.hypot2();
182 let d1 = self.p1 - self.p0;
183 let c = d1.hypot2();
184 if a < 5e-4 * c {
185 // This case happens for nearly straight Béziers.
186 //
187 // Calculate arclength using Legendre-Gauss quadrature using formula from Behdad
188 // in https://github.com/Pomax/BezierInfo-2/issues/77
189 let v0 = (-0.492943519233745 * self.p0.to_vec2()
190 + 0.430331482911935 * self.p1.to_vec2()
191 + 0.0626120363218102 * self.p2.to_vec2())
192 .hypot();
193 let v1 = ((self.p2 - self.p0) * 0.4444444444444444).hypot();
194 let v2 = (-0.0626120363218102 * self.p0.to_vec2()
195 - 0.430331482911935 * self.p1.to_vec2()
196 + 0.492943519233745 * self.p2.to_vec2())
197 .hypot();
198 return v0 + v1 + v2;
199 }
200 let b = 2.0 * d2.dot(d1);
201
202 let sabc = (a + b + c).sqrt();
203 let a2 = a.powf(-0.5);
204 let a32 = a2.powi(3);
205 let c2 = 2.0 * c.sqrt();
206 let ba_c2 = b * a2 + c2;
207
208 let v0 = 0.25 * a2 * a2 * b * (2.0 * sabc - c2) + sabc;
209 // TODO: justify and fine-tune this exact constant.
210 if ba_c2 < 1e-13 {
211 // This case happens for Béziers with a sharp kink.
212 v0
213 } else {
214 v0 + 0.25
215 * a32
216 * (4.0 * c * a - b * b)
217 * (((2.0 * a + b) * a2 + 2.0 * sabc) / ba_c2).ln()
218 }
219 }
220 }
221
222 impl ParamCurveArea for QuadBez {
223 #[inline]
signed_area(&self) -> f64224 fn signed_area(&self) -> f64 {
225 (self.p0.x * (2.0 * self.p1.y + self.p2.y) + 2.0 * self.p1.x * (self.p2.y - self.p0.y)
226 - self.p2.x * (self.p0.y + 2.0 * self.p1.y))
227 * (1.0 / 6.0)
228 }
229 }
230
231 impl ParamCurveNearest for QuadBez {
232 /// Find nearest point, using analytical algorithm based on cubic root finding.
nearest(&self, p: Point, _accuracy: f64) -> (f64, f64)233 fn nearest(&self, p: Point, _accuracy: f64) -> (f64, f64) {
234 fn eval_t(p: Point, t_best: &mut f64, r_best: &mut Option<f64>, t: f64, p0: Point) {
235 let r = (p0 - p).hypot2();
236 if r_best.map(|r_best| r < r_best).unwrap_or(true) {
237 *r_best = Some(r);
238 *t_best = t;
239 }
240 }
241 fn try_t(
242 q: &QuadBez,
243 p: Point,
244 t_best: &mut f64,
245 r_best: &mut Option<f64>,
246 t: f64,
247 ) -> bool {
248 if t < 0.0 || t > 1.0 {
249 return true;
250 }
251 eval_t(p, t_best, r_best, t, q.eval(t));
252 false
253 }
254 let d0 = self.p1 - self.p0;
255 let d1 = self.p0.to_vec2() + self.p2.to_vec2() - 2.0 * self.p1.to_vec2();
256 let d = self.p0 - p;
257 let c0 = d.dot(d0);
258 let c1 = 2.0 * d0.hypot2() + d.dot(d1);
259 let c2 = 3.0 * d1.dot(d0);
260 let c3 = d1.hypot2();
261 let roots = solve_cubic(c0, c1, c2, c3);
262 let mut r_best = None;
263 let mut t_best = 0.0;
264 let mut need_ends = false;
265 for &t in &roots {
266 need_ends |= try_t(self, p, &mut t_best, &mut r_best, t);
267 }
268 if need_ends {
269 eval_t(p, &mut t_best, &mut r_best, 0.0, self.p0);
270 eval_t(p, &mut t_best, &mut r_best, 1.0, self.p2);
271 }
272 (t_best, r_best.unwrap())
273 }
274 }
275
276 impl ParamCurveCurvature for QuadBez {}
277
278 impl ParamCurveExtrema for QuadBez {
extrema(&self) -> ArrayVec<[f64; MAX_EXTREMA]>279 fn extrema(&self) -> ArrayVec<[f64; MAX_EXTREMA]> {
280 let mut result = ArrayVec::new();
281 let d0 = self.p1 - self.p0;
282 let d1 = self.p2 - self.p1;
283 let dd = d1 - d0;
284 if dd.x != 0.0 {
285 let t = -d0.x / dd.x;
286 if t > 0.0 && t < 1.0 {
287 result.push(t);
288 }
289 }
290 if dd.y != 0.0 {
291 let t = -d0.y / dd.y;
292 if t > 0.0 && t < 1.0 {
293 result.push(t);
294 if result.len() == 2 && result[0] > t {
295 result.swap(0, 1);
296 }
297 }
298 }
299 result
300 }
301 }
302
303 impl Mul<QuadBez> for Affine {
304 type Output = QuadBez;
305
306 #[inline]
mul(self, other: QuadBez) -> QuadBez307 fn mul(self, other: QuadBez) -> QuadBez {
308 QuadBez {
309 p0: self * other.p0,
310 p1: self * other.p1,
311 p2: self * other.p2,
312 }
313 }
314 }
315
316 #[cfg(test)]
317 mod tests {
318 use crate::{
319 Affine, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveDeriv, ParamCurveExtrema,
320 ParamCurveNearest, Point, QuadBez,
321 };
322
assert_near(p0: Point, p1: Point, epsilon: f64)323 fn assert_near(p0: Point, p1: Point, epsilon: f64) {
324 assert!((p1 - p0).hypot() < epsilon, "{:?} != {:?}", p0, p1);
325 }
326
327 #[test]
quadbez_deriv()328 fn quadbez_deriv() {
329 let q = QuadBez::new((0.0, 0.0), (0.0, 0.5), (1.0, 1.0));
330 let deriv = q.deriv();
331
332 let n = 10;
333 for i in 0..=n {
334 let t = (i as f64) * (n as f64).recip();
335 let delta = 1e-6;
336 let p = q.eval(t);
337 let p1 = q.eval(t + delta);
338 let d_approx = (p1 - p) * delta.recip();
339 let d = deriv.eval(t).to_vec2();
340 assert!((d - d_approx).hypot() < delta * 2.0);
341 }
342 }
343
344 #[test]
quadbez_arclen()345 fn quadbez_arclen() {
346 let q = QuadBez::new((0.0, 0.0), (0.0, 0.5), (1.0, 1.0));
347 let true_arclen = 0.5 * 5.0f64.sqrt() + 0.25 * (2.0 + 5.0f64.sqrt()).ln();
348 for i in 0..12 {
349 let accuracy = 0.1f64.powi(i);
350 let est = q.arclen(accuracy);
351 let error = est - true_arclen;
352 assert!(error.abs() < accuracy, "{} != {}", est, true_arclen);
353 }
354 }
355
356 #[test]
quadbez_arclen_pathological()357 fn quadbez_arclen_pathological() {
358 let q = QuadBez::new((-1.0, 0.0), (1.03, 0.0), (1.0, 0.0));
359 let true_arclen = 2.0008737864167325; // A rough empirical calculation
360 let accuracy = 1e-11;
361 let est = q.arclen(accuracy);
362 assert!(
363 (est - true_arclen).abs() < accuracy,
364 "{} != {}",
365 est,
366 true_arclen
367 );
368 }
369
370 #[test]
quadbez_subsegment()371 fn quadbez_subsegment() {
372 let q = QuadBez::new((3.1, 4.1), (5.9, 2.6), (5.3, 5.8));
373 let t0 = 0.1;
374 let t1 = 0.8;
375 let qs = q.subsegment(t0..t1);
376 let epsilon = 1e-12;
377 let n = 10;
378 for i in 0..=n {
379 let t = (i as f64) * (n as f64).recip();
380 let ts = t0 + t * (t1 - t0);
381 assert_near(q.eval(ts), qs.eval(t), epsilon);
382 }
383 }
384
385 #[test]
quadbez_raise()386 fn quadbez_raise() {
387 let q = QuadBez::new((3.1, 4.1), (5.9, 2.6), (5.3, 5.8));
388 let c = q.raise();
389 let qd = q.deriv();
390 let cd = c.deriv();
391 let epsilon = 1e-12;
392 let n = 10;
393 for i in 0..=n {
394 let t = (i as f64) * (n as f64).recip();
395 assert_near(q.eval(t), c.eval(t), epsilon);
396 assert_near(qd.eval(t), cd.eval(t), epsilon);
397 }
398 }
399
400 #[test]
quadbez_signed_area()401 fn quadbez_signed_area() {
402 // y = 1 - x^2
403 let q = QuadBez::new((1.0, 0.0), (0.5, 1.0), (0.0, 1.0));
404 let epsilon = 1e-12;
405 assert!((q.signed_area() - 2.0 / 3.0).abs() < epsilon);
406 assert!(((Affine::rotate(0.5) * q).signed_area() - 2.0 / 3.0).abs() < epsilon);
407 assert!(((Affine::translate((0.0, 1.0)) * q).signed_area() - 3.5 / 3.0).abs() < epsilon);
408 assert!(((Affine::translate((1.0, 0.0)) * q).signed_area() - 3.5 / 3.0).abs() < epsilon);
409 }
410
411 #[test]
quadbez_nearest()412 fn quadbez_nearest() {
413 fn verify(result: (f64, f64), expected: f64) {
414 assert!(
415 (result.0 - expected).abs() < 1e-6,
416 "got {:?} expected {}",
417 result,
418 expected
419 );
420 }
421 // y = x^2
422 let q = QuadBez::new((-1.0, 1.0), (0.0, -1.0), (1.0, 1.0));
423 verify(q.nearest((0.0, 0.0).into(), 1e-3), 0.5);
424 verify(q.nearest((0.0, 0.1).into(), 1e-3), 0.5);
425 verify(q.nearest((0.0, -0.1).into(), 1e-3), 0.5);
426 verify(q.nearest((0.5, 0.25).into(), 1e-3), 0.75);
427 verify(q.nearest((1.0, 1.0).into(), 1e-3), 1.0);
428 verify(q.nearest((1.1, 1.1).into(), 1e-3), 1.0);
429 verify(q.nearest((-1.1, 1.1).into(), 1e-3), 0.0);
430 let a = Affine::rotate(0.5);
431 verify((a * q).nearest(a * Point::new(0.5, 0.25), 1e-3), 0.75);
432 }
433
434 // This test exposes a degenerate case in the solver used internally
435 // by the "nearest" calculation - the cubic term is zero.
436 #[test]
quadbez_nearest_low_order()437 fn quadbez_nearest_low_order() {
438 fn verify(result: (f64, f64), expected: f64) {
439 assert!(
440 (result.0 - expected).abs() < 1e-6,
441 "got {:?} expected {}",
442 result,
443 expected
444 );
445 }
446
447 let q = QuadBez::new((-1.0, 0.0), (0.0, 0.0), (1.0, 0.0));
448
449 verify(q.nearest((0.0, 0.0).into(), 1e-3), 0.5);
450 verify(q.nearest((0.0, 1.0).into(), 1e-3), 0.5);
451 }
452
453 #[test]
quadbez_extrema()454 fn quadbez_extrema() {
455 // y = x^2
456 let q = QuadBez::new((-1.0, 1.0), (0.0, -1.0), (1.0, 1.0));
457 let extrema = q.extrema();
458 assert_eq!(extrema.len(), 1);
459 assert!((extrema[0] - 0.5).abs() < 1e-6);
460
461 let q = QuadBez::new((0.0, 0.5), (1.0, 1.0), (0.5, 0.0));
462 let extrema = q.extrema();
463 assert_eq!(extrema.len(), 2);
464 assert!((extrema[0] - 1.0 / 3.0).abs() < 1e-6);
465 assert!((extrema[1] - 2.0 / 3.0).abs() < 1e-6);
466
467 // Reverse direction
468 let q = QuadBez::new((0.5, 0.0), (1.0, 1.0), (0.0, 0.5));
469 let extrema = q.extrema();
470 assert_eq!(extrema.len(), 2);
471 assert!((extrema[0] - 1.0 / 3.0).abs() < 1e-6);
472 assert!((extrema[1] - 2.0 / 3.0).abs() < 1e-6);
473 }
474 }
475