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