1// Copyright 2018 the Kurbo Authors
2// SPDX-License-Identifier: Apache-2.0 OR MIT
3
4//! Cubic Bézier segments.
5
6use alloc::vec;
7use alloc::vec::Vec;
8use core::ops::{Mul, Range};
9
10use crate::MAX_EXTREMA;
11use crate::{Line, QuadSpline, Vec2};
12use arrayvec::ArrayVec;
13
14use crate::common::{
15 solve_quadratic, solve_quartic, GAUSS_LEGENDRE_COEFFS_16_HALF, GAUSS_LEGENDRE_COEFFS_24_HALF,
16 GAUSS_LEGENDRE_COEFFS_8, GAUSS_LEGENDRE_COEFFS_8_HALF,
17};
18use crate::{
19 Affine, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveCurvature,
20 ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point, QuadBez, Rect, Shape,
21};
22
23#[cfg(not(feature = "std"))]
24use crate::common::FloatFuncs;
25
26const MAX_SPLINE_SPLIT: usize = 100;
27
28/// A single cubic Bézier segment.
29#[derive(Clone, Copy, Debug, PartialEq)]
30#[cfg_attr(feature = "schemars", derive(schemars::JsonSchema))]
31#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
32#[allow(missing_docs)]
33pub struct CubicBez {
34 pub p0: Point,
35 pub p1: Point,
36 pub p2: Point,
37 pub p3: Point,
38}
39
40/// An iterator which produces quadratic Bézier segments.
41struct ToQuads {
42 c: CubicBez,
43 i: usize,
44 n: usize,
45}
46
47/// Classification result for cusp detection.
48#[derive(Debug)]
49pub enum CuspType {
50 /// Cusp is a loop.
51 Loop,
52 /// Cusp has two closely spaced inflection points.
53 DoubleInflection,
54}
55
56impl CubicBez {
57 /// Create a new cubic Bézier segment.
58 #[inline]
59 pub fn new<P: Into<Point>>(p0: P, p1: P, p2: P, p3: P) -> CubicBez {
60 CubicBez {
61 p0: p0.into(),
62 p1: p1.into(),
63 p2: p2.into(),
64 p3: p3.into(),
65 }
66 }
67
68 /// Convert to quadratic Béziers.
69 ///
70 /// The iterator returns the start and end parameter in the cubic of each quadratic
71 /// segment, along with the quadratic.
72 ///
73 /// Note that the resulting quadratic Béziers are not in general G1 continuous;
74 /// they are optimized for minimizing distance error.
75 ///
76 /// This iterator will always produce at least one `QuadBez`.
77 #[inline]
78 pub fn to_quads(&self, accuracy: f64) -> impl Iterator<Item = (f64, f64, QuadBez)> {
79 // The maximum error, as a vector from the cubic to the best approximating
80 // quadratic, is proportional to the third derivative, which is constant
81 // across the segment. Thus, the error scales down as the third power of
82 // the number of subdivisions. Our strategy then is to subdivide `t` evenly.
83 //
84 // This is an overestimate of the error because only the component
85 // perpendicular to the first derivative is important. But the simplicity is
86 // appealing.
87
88 // This magic number is the square of 36 / sqrt(3).
89 // See: http://caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
90 let max_hypot2 = 432.0 * accuracy * accuracy;
91 let p1x2 = 3.0 * self.p1.to_vec2() - self.p0.to_vec2();
92 let p2x2 = 3.0 * self.p2.to_vec2() - self.p3.to_vec2();
93 let err = (p2x2 - p1x2).hypot2();
94 let n = ((err / max_hypot2).powf(1. / 6.0).ceil() as usize).max(1);
95
96 ToQuads { c: *self, n, i: 0 }
97 }
98
99 /// Return a [`QuadSpline`] approximating this cubic Bézier.
100 ///
101 /// Returns `None` if no suitable approximation is found within the given
102 /// tolerance.
103 pub fn approx_spline(&self, accuracy: f64) -> Option<QuadSpline> {
104 (1..=MAX_SPLINE_SPLIT).find_map(|n| self.approx_spline_n(n, accuracy))
105 }
106
107 // Approximate a cubic curve with a quadratic spline of `n` curves
108 fn approx_spline_n(&self, n: usize, accuracy: f64) -> Option<QuadSpline> {
109 if n == 1 {
110 return self
111 .try_approx_quadratic(accuracy)
112 .map(|quad| QuadSpline::new(vec![quad.p0, quad.p1, quad.p2]));
113 }
114 let mut cubics = self.split_into_n(n);
115
116 // The above function guarantees that the iterator returns n items,
117 // which is why we're unwrapping things with wild abandon.
118 let mut next_cubic = cubics.next().unwrap();
119 let mut next_q1: Point = next_cubic.approx_quad_control(0.0);
120 let mut q2 = self.p0;
121 let mut d1 = Vec2::ZERO;
122 let mut spline = vec![self.p0, next_q1];
123 for i in 1..=n {
124 let current_cubic: CubicBez = next_cubic;
125 let q0 = q2;
126 let q1 = next_q1;
127 q2 = if i < n {
128 next_cubic = cubics.next().unwrap();
129 next_q1 = next_cubic.approx_quad_control(i as f64 / (n - 1) as f64);
130
131 spline.push(next_q1);
132 q1.midpoint(next_q1)
133 } else {
134 current_cubic.p3
135 };
136 let d0 = d1;
137 d1 = q2.to_vec2() - current_cubic.p3.to_vec2();
138
139 if d1.hypot() > accuracy
140 || !CubicBez::new(
141 d0.to_point(),
142 q0.lerp(q1, 2.0 / 3.0) - current_cubic.p1.to_vec2(),
143 q2.lerp(q1, 2.0 / 3.0) - current_cubic.p2.to_vec2(),
144 d1.to_point(),
145 )
146 .fit_inside(accuracy)
147 {
148 return None;
149 }
150 }
151 spline.push(self.p3);
152 Some(QuadSpline::new(spline))
153 }
154
155 fn approx_quad_control(&self, t: f64) -> Point {
156 let p1 = self.p0 + (self.p1 - self.p0) * 1.5;
157 let p2 = self.p3 + (self.p2 - self.p3) * 1.5;
158 p1.lerp(p2, t)
159 }
160
161 /// Approximate a cubic with a single quadratic
162 ///
163 /// Returns a quadratic approximating the given cubic that maintains
164 /// endpoint tangents if that is within tolerance, or `None` otherwise.
165 fn try_approx_quadratic(&self, accuracy: f64) -> Option<QuadBez> {
166 if let Some(q1) = Line::new(self.p0, self.p1).crossing_point(Line::new(self.p2, self.p3)) {
167 let c1 = self.p0.lerp(q1, 2.0 / 3.0);
168 let c2 = self.p3.lerp(q1, 2.0 / 3.0);
169 if !CubicBez::new(
170 Point::ZERO,
171 c1 - self.p1.to_vec2(),
172 c2 - self.p2.to_vec2(),
173 Point::ZERO,
174 )
175 .fit_inside(accuracy)
176 {
177 return None;
178 }
179 return Some(QuadBez::new(self.p0, q1, self.p3));
180 }
181 None
182 }
183
184 fn split_into_n(&self, n: usize) -> impl Iterator<Item = CubicBez> {
185 // for certain small values of `n` we precompute all our values.
186 // if we have six or fewer items we precompute them.
187 let mut storage = ArrayVec::<_, 6>::new();
188
189 match n {
190 1 => storage.push(*self),
191 2 => {
192 let (l, r) = self.subdivide();
193 storage.try_extend_from_slice(&[r, l]).unwrap();
194 }
195 3 => {
196 let (left, mid, right) = self.subdivide_3();
197 storage.try_extend_from_slice(&[right, mid, left]).unwrap();
198 }
199 4 => {
200 let (l, r) = self.subdivide();
201 let (ll, lr) = l.subdivide();
202 let (rl, rr) = r.subdivide();
203 storage.try_extend_from_slice(&[rr, rl, lr, ll]).unwrap();
204 }
205 6 => {
206 let (l, r) = self.subdivide();
207 let (l1, l2, l3) = l.subdivide_3();
208 let (r1, r2, r3) = r.subdivide_3();
209 storage
210 .try_extend_from_slice(&[r3, r2, r1, l3, l2, l1])
211 .unwrap();
212 }
213 _ => (),
214 }
215
216 // a limitation of returning 'impl Trait' is that the implementation
217 // can only return a single concrete type; that is you cannot return
218 // Vec::into_iter() from one branch, and then HashSet::into_iter from
219 // another branch.
220 //
221 // This means we have to get a bit fancy, and have a single concrete
222 // type that represents both of our possible cases.
223
224 let mut storage = if storage.is_empty() {
225 None
226 } else {
227 Some(storage)
228 };
229
230 // used in the fallback case
231 let mut i = 0;
232 let (a, b, c, d) = self.parameters();
233 let dt = 1.0 / n as f64;
234 let delta_2 = dt * dt;
235 let delta_3 = dt * delta_2;
236
237 core::iter::from_fn(move || {
238 // if storage exists, we use it exclusively
239 if let Some(storage) = storage.as_mut() {
240 return storage.pop();
241 }
242
243 // if storage does not exist, we are exclusively working down here.
244 if i >= n {
245 return None;
246 }
247
248 let t1 = i as f64 * dt;
249 let t1_2 = t1 * t1;
250 let a1 = a * delta_3;
251 let b1 = (3.0 * a * t1 + b) * delta_2;
252 let c1 = (2.0 * b * t1 + c + 3.0 * a * t1_2) * dt;
253 let d1 = a * t1 * t1_2 + b * t1_2 + c * t1 + d;
254 let result = CubicBez::from_parameters(a1, b1, c1, d1);
255 i += 1;
256 Some(result)
257 })
258 }
259
260 fn parameters(&self) -> (Vec2, Vec2, Vec2, Vec2) {
261 let c = (self.p1 - self.p0) * 3.0;
262 let b = (self.p2 - self.p1) * 3.0 - c;
263 let d = self.p0.to_vec2();
264 let a = self.p3.to_vec2() - d - c - b;
265 (a, b, c, d)
266 }
267
268 /// Rust port of cu2qu [calc_cubic_points](https://github.com/fonttools/fonttools/blob/3b9a73ff8379ab49d3ce35aaaaf04b3a7d9d1655/Lib/fontTools/cu2qu/cu2qu.py#L63-L68)
269 fn from_parameters(a: Vec2, b: Vec2, c: Vec2, d: Vec2) -> Self {
270 let p0 = d.to_point();
271 let p1 = c.div_exact(3.0).to_point() + d;
272 let p2 = (b + c).div_exact(3.0).to_point() + p1.to_vec2();
273 let p3 = (a + d + c + b).to_point();
274 CubicBez::new(p0, p1, p2, p3)
275 }
276
277 fn subdivide_3(&self) -> (CubicBez, CubicBez, CubicBez) {
278 let (p0, p1, p2, p3) = (
279 self.p0.to_vec2(),
280 self.p1.to_vec2(),
281 self.p2.to_vec2(),
282 self.p3.to_vec2(),
283 );
284 // The original Python cu2qu code here does not use division operator to divide by 27 but
285 // instead uses multiplication by the reciprocal 1 / 27. We want to match it exactly
286 // to avoid any floating point differences, hence in this particular case we do not use div_exact.
287 // I could directly use the Vec2 Div trait (also implemented as multiplication by reciprocal)
288 // but I prefer to be explicit here.
289 // Source: https://github.com/fonttools/fonttools/blob/85c80be/Lib/fontTools/cu2qu/cu2qu.py#L215-L218
290 // See also: https://github.com/linebender/kurbo/issues/272
291 let one_27th = 27.0_f64.recip();
292 let mid1 = ((8.0 * p0 + 12.0 * p1 + 6.0 * p2 + p3) * one_27th).to_point();
293 let deriv1 = (p3 + 3.0 * p2 - 4.0 * p0) * one_27th;
294 let mid2 = ((p0 + 6.0 * p1 + 12.0 * p2 + 8.0 * p3) * one_27th).to_point();
295 let deriv2 = (4.0 * p3 - 3.0 * p1 - p0) * one_27th;
296
297 let left = CubicBez::new(
298 self.p0,
299 (2.0 * p0 + p1).div_exact(3.0).to_point(),
300 mid1 - deriv1,
301 mid1,
302 );
303 let mid = CubicBez::new(mid1, mid1 + deriv1, mid2 - deriv2, mid2);
304 let right = CubicBez::new(
305 mid2,
306 mid2 + deriv2,
307 (p2 + 2.0 * p3).div_exact(3.0).to_point(),
308 self.p3,
309 );
310 (left, mid, right)
311 }
312
313 /// Does this curve fit inside the given distance from the origin?
314 ///
315 /// Rust port of cu2qu [cubic_farthest_fit_inside](https://github.com/fonttools/fonttools/blob/3b9a73ff8379ab49d3ce35aaaaf04b3a7d9d1655/Lib/fontTools/cu2qu/cu2qu.py#L281)
316 fn fit_inside(&self, distance: f64) -> bool {
317 if self.p2.to_vec2().hypot() <= distance && self.p1.to_vec2().hypot() <= distance {
318 return true;
319 }
320 let mid =
321 (self.p0.to_vec2() + 3.0 * (self.p1.to_vec2() + self.p2.to_vec2()) + self.p3.to_vec2())
322 * 0.125;
323 if mid.hypot() > distance {
324 return false;
325 }
326 // Split in two. Note that cu2qu here uses a 3/8 subdivision. I don't know why.
327 let (left, right) = self.subdivide();
328 left.fit_inside(distance) && right.fit_inside(distance)
329 }
330
331 /// Is this cubic Bezier curve finite?
332 #[inline]
333 pub fn is_finite(&self) -> bool {
334 self.p0.is_finite() && self.p1.is_finite() && self.p2.is_finite() && self.p3.is_finite()
335 }
336
337 /// Is this cubic Bezier curve NaN?
338 #[inline]
339 pub fn is_nan(&self) -> bool {
340 self.p0.is_nan() || self.p1.is_nan() || self.p2.is_nan() || self.p3.is_nan()
341 }
342
343 /// Determine the inflection points.
344 ///
345 /// Return value is t parameter for the inflection points of the curve segment.
346 /// There are a maximum of two for a cubic Bézier.
347 ///
348 /// See <https://www.caffeineowl.com/graphics/2d/vectorial/cubic-inflexion.html>
349 /// for the theory.
350 pub fn inflections(&self) -> ArrayVec<f64, 2> {
351 let a = self.p1 - self.p0;
352 let b = (self.p2 - self.p1) - a;
353 let c = (self.p3 - self.p0) - 3. * (self.p2 - self.p1);
354 solve_quadratic(a.cross(b), a.cross(c), b.cross(c))
355 .iter()
356 .copied()
357 .filter(|t| *t >= 0.0 && *t <= 1.0)
358 .collect()
359 }
360
361 /// Find points on the curve where the tangent line passes through the
362 /// given point.
363 ///
364 /// Result is array of t values such that the tangent line from the curve
365 /// evaluated at that point goes through the argument point.
366 pub fn tangents_to_point(&self, p: Point) -> ArrayVec<f64, 4> {
367 let (a, b, c, d_orig) = self.parameters();
368 let d = d_orig - p.to_vec2();
369 // coefficients of x(t) \cross x'(t)
370 let c4 = b.cross(a);
371 let c3 = 2.0 * c.cross(a);
372 let c2 = c.cross(b) + 3.0 * d.cross(a);
373 let c1 = 2.0 * d.cross(b);
374 let c0 = d.cross(c);
375 solve_quartic(c0, c1, c2, c3, c4)
376 .iter()
377 .copied()
378 .filter(|t| *t >= 0.0 && *t <= 1.0)
379 .collect()
380 }
381
382 /// Preprocess a cubic Bézier to ease numerical robustness.
383 ///
384 /// If the cubic Bézier segment has zero or near-zero derivatives, perturb
385 /// the control points to make it easier to process (especially offset and
386 /// stroke), avoiding numerical robustness problems.
387 pub(crate) fn regularize(&self, dimension: f64) -> CubicBez {
388 let mut c = *self;
389 // First step: if control point is too near the endpoint, nudge it away
390 // along the tangent.
391 let dim2 = dimension * dimension;
392 if c.p0.distance_squared(c.p1) < dim2 {
393 let d02 = c.p0.distance_squared(c.p2);
394 if d02 >= dim2 {
395 // TODO: moderate if this would move closer to p3
396 c.p1 = c.p0.lerp(c.p2, (dim2 / d02).sqrt());
397 } else {
398 c.p1 = c.p0.lerp(c.p3, 1.0 / 3.0);
399 c.p2 = c.p3.lerp(c.p0, 1.0 / 3.0);
400 return c;
401 }
402 }
403 if c.p3.distance_squared(c.p2) < dim2 {
404 let d13 = c.p1.distance_squared(c.p2);
405 if d13 >= dim2 {
406 // TODO: moderate if this would move closer to p0
407 c.p2 = c.p3.lerp(c.p1, (dim2 / d13).sqrt());
408 } else {
409 c.p1 = c.p0.lerp(c.p3, 1.0 / 3.0);
410 c.p2 = c.p3.lerp(c.p0, 1.0 / 3.0);
411 return c;
412 }
413 }
414 if let Some(cusp_type) = self.detect_cusp(dimension) {
415 let d01 = c.p1 - c.p0;
416 let d01h = d01.hypot();
417 let d23 = c.p3 - c.p2;
418 let d23h = d23.hypot();
419 match cusp_type {
420 CuspType::Loop => {
421 c.p1 += (dimension / d01h) * d01;
422 c.p2 -= (dimension / d23h) * d23;
423 }
424 CuspType::DoubleInflection => {
425 // Avoid making control distance smaller than dimension
426 if d01h > 2.0 * dimension {
427 c.p1 -= (dimension / d01h) * d01;
428 }
429 if d23h > 2.0 * dimension {
430 c.p2 += (dimension / d23h) * d23;
431 }
432 }
433 }
434 }
435 c
436 }
437
438 /// Detect whether there is a cusp.
439 ///
440 /// Return a cusp classification if there is a cusp with curvature greater than
441 /// the reciprocal of the given dimension.
442 fn detect_cusp(&self, dimension: f64) -> Option<CuspType> {
443 let d01 = self.p1 - self.p0;
444 let d02 = self.p2 - self.p0;
445 let d03 = self.p3 - self.p0;
446 let d12 = self.p2 - self.p1;
447 let d23 = self.p3 - self.p2;
448 let det_012 = d01.cross(d02);
449 let det_123 = d12.cross(d23);
450 let det_013 = d01.cross(d03);
451 let det_023 = d02.cross(d03);
452 if det_012 * det_123 > 0.0 && det_012 * det_013 < 0.0 && det_012 * det_023 < 0.0 {
453 let q = self.deriv();
454 // accuracy isn't used for quadratic nearest
455 let nearest = q.nearest(Point::ORIGIN, 1e-9);
456 // detect whether curvature at minimum derivative exceeds 1/dimension,
457 // without division.
458 let d = q.eval(nearest.t);
459 let d2 = q.deriv().eval(nearest.t);
460 let cross = d.to_vec2().cross(d2.to_vec2());
461 if nearest.distance_sq.powi(3) <= (cross * dimension).powi(2) {
462 let a = 3. * det_012 + det_023 - 2. * det_013;
463 let b = -3. * det_012 + det_013;
464 let c = det_012;
465 let d = b * b - 4. * a * c;
466 if d > 0.0 {
467 return Some(CuspType::DoubleInflection);
468 } else {
469 return Some(CuspType::Loop);
470 }
471 }
472 }
473 None
474 }
475}
476
477/// An iterator for cubic beziers.
478pub struct CubicBezIter {
479 cubic: CubicBez,
480 ix: usize,
481}
482
483impl Shape for CubicBez {
484 type PathElementsIter<'iter> = CubicBezIter;
485
486 #[inline]
487 fn path_elements(&self, _tolerance: f64) -> CubicBezIter {
488 CubicBezIter {
489 cubic: *self,
490 ix: 0,
491 }
492 }
493
494 fn area(&self) -> f64 {
495 0.0
496 }
497
498 #[inline]
499 fn perimeter(&self, accuracy: f64) -> f64 {
500 self.arclen(accuracy)
501 }
502
503 fn winding(&self, _pt: Point) -> i32 {
504 0
505 }
506
507 #[inline]
508 fn bounding_box(&self) -> Rect {
509 ParamCurveExtrema::bounding_box(self)
510 }
511}
512
513impl Iterator for CubicBezIter {
514 type Item = PathEl;
515
516 fn next(&mut self) -> Option<PathEl> {
517 self.ix += 1;
518 match self.ix {
519 1 => Some(PathEl::MoveTo(self.cubic.p0)),
520 2 => Some(PathEl::CurveTo(self.cubic.p1, self.cubic.p2, self.cubic.p3)),
521 _ => None,
522 }
523 }
524}
525
526impl ParamCurve for CubicBez {
527 #[inline]
528 fn eval(&self, t: f64) -> Point {
529 let mt = 1.0 - t;
530 let v = self.p0.to_vec2() * (mt * mt * mt)
531 + (self.p1.to_vec2() * (mt * mt * 3.0)
532 + (self.p2.to_vec2() * (mt * 3.0) + self.p3.to_vec2() * t) * t)
533 * t;
534 v.to_point()
535 }
536
537 fn subsegment(&self, range: Range<f64>) -> CubicBez {
538 let (t0, t1) = (range.start, range.end);
539 let p0 = self.eval(t0);
540 let p3 = self.eval(t1);
541 let d = self.deriv();
542 let scale = (t1 - t0) * (1.0 / 3.0);
543 let p1 = p0 + scale * d.eval(t0).to_vec2();
544 let p2 = p3 - scale * d.eval(t1).to_vec2();
545 CubicBez { p0, p1, p2, p3 }
546 }
547
548 /// Subdivide into halves, using de Casteljau.
549 #[inline]
550 fn subdivide(&self) -> (CubicBez, CubicBez) {
551 let pm = self.eval(0.5);
552 (
553 CubicBez::new(
554 self.p0,
555 self.p0.midpoint(self.p1),
556 ((self.p0.to_vec2() + self.p1.to_vec2() * 2.0 + self.p2.to_vec2()) * 0.25)
557 .to_point(),
558 pm,
559 ),
560 CubicBez::new(
561 pm,
562 ((self.p1.to_vec2() + self.p2.to_vec2() * 2.0 + self.p3.to_vec2()) * 0.25)
563 .to_point(),
564 self.p2.midpoint(self.p3),
565 self.p3,
566 ),
567 )
568 }
569
570 #[inline]
571 fn start(&self) -> Point {
572 self.p0
573 }
574
575 #[inline]
576 fn end(&self) -> Point {
577 self.p3
578 }
579}
580
581impl ParamCurveDeriv for CubicBez {
582 type DerivResult = QuadBez;
583
584 #[inline]
585 fn deriv(&self) -> QuadBez {
586 QuadBez::new(
587 (3.0 * (self.p1 - self.p0)).to_point(),
588 (3.0 * (self.p2 - self.p1)).to_point(),
589 (3.0 * (self.p3 - self.p2)).to_point(),
590 )
591 }
592}
593
594fn arclen_quadrature_core(coeffs: &[(f64, f64)], dm: Vec2, dm1: Vec2, dm2: Vec2) -> f64 {
595 coeffsimpl Iterator
596 .iter()
597 .map(|&(wi: f64, xi: f64)| {
598 let d: Vec2 = dm + dm2 * (xi * xi);
599 let dpx: f64 = (d + dm1 * xi).hypot();
600 let dmx: f64 = (d - dm1 * xi).hypot();
601 (2.25f64.sqrt() * wi) * (dpx + dmx)
602 })
603 .sum::<f64>()
604}
605
606fn arclen_rec(c: &CubicBez, accuracy: f64, depth: usize) -> f64 {
607 let d03 = c.p3 - c.p0;
608 let d01 = c.p1 - c.p0;
609 let d12 = c.p2 - c.p1;
610 let d23 = c.p3 - c.p2;
611 let lp_lc = d01.hypot() + d12.hypot() + d23.hypot() - d03.hypot();
612 let dd1 = d12 - d01;
613 let dd2 = d23 - d12;
614 // It might be faster to do direct multiplies, the data dependencies would be shorter.
615 // The following values don't have the factor of 3 for first deriv
616 let dm = 0.25 * (d01 + d23) + 0.5 * d12; // first derivative at midpoint
617 let dm1 = 0.5 * (dd2 + dd1); // second derivative at midpoint
618 let dm2 = 0.25 * (dd2 - dd1); // 0.5 * (third derivative at midpoint)
619
620 let est = GAUSS_LEGENDRE_COEFFS_8
621 .iter()
622 .map(|&(wi, xi)| {
623 wi * {
624 let d_norm2 = (dm + dm1 * xi + dm2 * (xi * xi)).hypot2();
625 let dd_norm2 = (dm1 + dm2 * (2.0 * xi)).hypot2();
626 dd_norm2 / d_norm2
627 }
628 })
629 .sum::<f64>();
630 let est_gauss8_error = (est.powi(3) * 2.5e-6).min(3e-2) * lp_lc;
631 if est_gauss8_error < accuracy {
632 return arclen_quadrature_core(GAUSS_LEGENDRE_COEFFS_8_HALF, dm, dm1, dm2);
633 }
634 let est_gauss16_error = (est.powi(6) * 1.5e-11).min(9e-3) * lp_lc;
635 if est_gauss16_error < accuracy {
636 return arclen_quadrature_core(GAUSS_LEGENDRE_COEFFS_16_HALF, dm, dm1, dm2);
637 }
638 let est_gauss24_error = (est.powi(9) * 3.5e-16).min(3.5e-3) * lp_lc;
639 if est_gauss24_error < accuracy || depth >= 20 {
640 return arclen_quadrature_core(GAUSS_LEGENDRE_COEFFS_24_HALF, dm, dm1, dm2);
641 }
642 let (c0, c1) = c.subdivide();
643 arclen_rec(&c0, accuracy * 0.5, depth + 1) + arclen_rec(&c1, accuracy * 0.5, depth + 1)
644}
645
646impl ParamCurveArclen for CubicBez {
647 /// Arclength of a cubic Bézier segment.
648 ///
649 /// This is an adaptive subdivision approach using Legendre-Gauss quadrature
650 /// in the base case, and an error estimate to decide when to subdivide.
651 fn arclen(&self, accuracy: f64) -> f64 {
652 arclen_rec(self, accuracy, depth:0)
653 }
654}
655
656impl ParamCurveArea for CubicBez {
657 #[inline]
658 fn signed_area(&self) -> f64 {
659 (self.p0.x * (6.0 * self.p1.y + 3.0 * self.p2.y + self.p3.y)
660 + 3.0
661 * (self.p1.x * (-2.0 * self.p0.y + self.p2.y + self.p3.y)
662 - self.p2.x * (self.p0.y + self.p1.y - 2.0 * self.p3.y))
663 - self.p3.x * (self.p0.y + 3.0 * self.p1.y + 6.0 * self.p2.y))
664 * (1.0 / 20.0)
665 }
666}
667
668impl ParamCurveNearest for CubicBez {
669 /// Find the nearest point, using subdivision.
670 fn nearest(&self, p: Point, accuracy: f64) -> Nearest {
671 let mut best_r: Option = None;
672 let mut best_t: f64 = 0.0;
673 for (t0: f64, t1: f64, q: QuadBez) in self.to_quads(accuracy) {
674 let nearest: Nearest = q.nearest(p, accuracy);
675 if best_r
676 .map(|best_r| nearest.distance_sq < best_r)
677 .unwrap_or(default:true)
678 {
679 best_t = t0 + nearest.t * (t1 - t0);
680 best_r = Some(nearest.distance_sq);
681 }
682 }
683 Nearest {
684 t: best_t,
685 distance_sq: best_r.unwrap(),
686 }
687 }
688}
689
690impl ParamCurveCurvature for CubicBez {}
691
692impl ParamCurveExtrema for CubicBez {
693 fn extrema(&self) -> ArrayVec<f64, MAX_EXTREMA> {
694 fn one_coord(result: &mut ArrayVec<f64, MAX_EXTREMA>, d0: f64, d1: f64, d2: f64) {
695 let a: f64 = d0 - 2.0 * d1 + d2;
696 let b: f64 = 2.0 * (d1 - d0);
697 let c: f64 = d0;
698 let roots: ArrayVec = solve_quadratic(c0:c, c1:b, c2:a);
699 for &t: f64 in &roots {
700 if t > 0.0 && t < 1.0 {
701 result.push(element:t);
702 }
703 }
704 }
705 let mut result: ArrayVec = ArrayVec::new();
706 let d0: Vec2 = self.p1 - self.p0;
707 let d1: Vec2 = self.p2 - self.p1;
708 let d2: Vec2 = self.p3 - self.p2;
709 one_coord(&mut result, d0.x, d1.x, d2.x);
710 one_coord(&mut result, d0.y, d1.y, d2.y);
711 result.sort_by(|a: &f64, b: &f64| a.partial_cmp(b).unwrap());
712 result
713 }
714}
715
716impl Mul<CubicBez> for Affine {
717 type Output = CubicBez;
718
719 #[inline]
720 fn mul(self, c: CubicBez) -> CubicBez {
721 CubicBez {
722 p0: self * c.p0,
723 p1: self * c.p1,
724 p2: self * c.p2,
725 p3: self * c.p3,
726 }
727 }
728}
729
730impl Iterator for ToQuads {
731 type Item = (f64, f64, QuadBez);
732
733 fn next(&mut self) -> Option<(f64, f64, QuadBez)> {
734 if self.i == self.n {
735 return None;
736 }
737 let t0: f64 = self.i as f64 / self.n as f64;
738 let t1: f64 = (self.i + 1) as f64 / self.n as f64;
739 let seg: CubicBez = self.c.subsegment(range:t0..t1);
740 let p1x2: Vec2 = 3.0 * seg.p1.to_vec2() - seg.p0.to_vec2();
741 let p2x2: Vec2 = 3.0 * seg.p2.to_vec2() - seg.p3.to_vec2();
742 let result: QuadBez = QuadBez::new(seg.p0, ((p1x2 + p2x2) / 4.0).to_point(), p2:seg.p3);
743 self.i += 1;
744 Some((t0, t1, result))
745 }
746
747 fn size_hint(&self) -> (usize, Option<usize>) {
748 let remaining: usize = self.n - self.i;
749 (remaining, Some(remaining))
750 }
751}
752
753/// Convert multiple cubic Bézier curves to quadratic splines.
754///
755/// Ensures that the resulting splines have the same number of control points.
756///
757/// Rust port of cu2qu [cubic_approx_quadratic](https://github.com/fonttools/fonttools/blob/3b9a73ff8379ab49d3ce35aaaaf04b3a7d9d1655/Lib/fontTools/cu2qu/cu2qu.py#L322)
758pub fn cubics_to_quadratic_splines(curves: &[CubicBez], accuracy: f64) -> Option<Vec<QuadSpline>> {
759 let mut result: Vec = Vec::new();
760 let mut split_order: usize = 0;
761
762 while split_order <= MAX_SPLINE_SPLIT {
763 split_order += 1;
764 result.clear();
765
766 for curve: &CubicBez in curves {
767 match curve.approx_spline_n(n:split_order, accuracy) {
768 Some(spline: QuadSpline) => result.push(spline),
769 None => break,
770 }
771 }
772
773 if result.len() == curves.len() {
774 return Some(result);
775 }
776 }
777 None
778}
779#[cfg(test)]
780mod tests {
781 use crate::{
782 cubics_to_quadratic_splines, Affine, CubicBez, Nearest, ParamCurve, ParamCurveArclen,
783 ParamCurveArea, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, Point, QuadBez,
784 QuadSpline,
785 };
786
787 #[test]
788 fn cubicbez_deriv() {
789 // y = x^2
790 let c = CubicBez::new(
791 (0.0, 0.0),
792 (1.0 / 3.0, 0.0),
793 (2.0 / 3.0, 1.0 / 3.0),
794 (1.0, 1.0),
795 );
796 let deriv = c.deriv();
797
798 let n = 10;
799 for i in 0..=n {
800 let t = (i as f64) * (n as f64).recip();
801 let delta = 1e-6;
802 let p = c.eval(t);
803 let p1 = c.eval(t + delta);
804 let d_approx = (p1 - p) * delta.recip();
805 let d = deriv.eval(t).to_vec2();
806 assert!((d - d_approx).hypot() < delta * 2.0);
807 }
808 }
809
810 #[test]
811 fn cubicbez_arclen() {
812 // y = x^2
813 let c = CubicBez::new(
814 (0.0, 0.0),
815 (1.0 / 3.0, 0.0),
816 (2.0 / 3.0, 1.0 / 3.0),
817 (1.0, 1.0),
818 );
819 let true_arclen = 0.5 * 5.0f64.sqrt() + 0.25 * (2.0 + 5.0f64.sqrt()).ln();
820 for i in 0..12 {
821 let accuracy = 0.1f64.powi(i);
822 let error = c.arclen(accuracy) - true_arclen;
823 assert!(error.abs() < accuracy);
824 }
825 }
826
827 #[test]
828 fn cubicbez_inv_arclen() {
829 // y = x^2 / 100
830 let c = CubicBez::new(
831 (0.0, 0.0),
832 (100.0 / 3.0, 0.0),
833 (200.0 / 3.0, 100.0 / 3.0),
834 (100.0, 100.0),
835 );
836 let true_arclen = 100.0 * (0.5 * 5.0f64.sqrt() + 0.25 * (2.0 + 5.0f64.sqrt()).ln());
837 for i in 0..12 {
838 let accuracy = 0.1f64.powi(i);
839 let n = 10;
840 for j in 0..=n {
841 let arc = (j as f64) * ((n as f64).recip() * true_arclen);
842 let t = c.inv_arclen(arc, accuracy * 0.5);
843 let actual_arc = c.subsegment(0.0..t).arclen(accuracy * 0.5);
844 assert!(
845 (arc - actual_arc).abs() < accuracy,
846 "at accuracy {accuracy:e}, wanted {actual_arc} got {arc}"
847 );
848 }
849 }
850 // corner case: user passes accuracy larger than total arc length
851 let accuracy = true_arclen * 1.1;
852 let arc = true_arclen * 0.5;
853 let t = c.inv_arclen(arc, accuracy);
854 let actual_arc = c.subsegment(0.0..t).arclen(accuracy);
855 assert!(
856 (arc - actual_arc).abs() < 2.0 * accuracy,
857 "at accuracy {accuracy:e}, want {actual_arc} got {arc}"
858 );
859 }
860
861 #[test]
862 fn cubicbez_inv_arclen_accuracy() {
863 let c = CubicBez::new((0.2, 0.73), (0.35, 1.08), (0.85, 1.08), (1.0, 0.73));
864 let true_t = c.inv_arclen(0.5, 1e-12);
865 for i in 1..12 {
866 let accuracy = (0.1f64).powi(i);
867 let approx_t = c.inv_arclen(0.5, accuracy);
868 assert!((approx_t - true_t).abs() <= accuracy);
869 }
870 }
871
872 #[test]
873 #[allow(clippy::float_cmp)]
874 fn cubicbez_signed_area_linear() {
875 // y = 1 - x
876 let c = CubicBez::new(
877 (1.0, 0.0),
878 (2.0 / 3.0, 1.0 / 3.0),
879 (1.0 / 3.0, 2.0 / 3.0),
880 (0.0, 1.0),
881 );
882 let epsilon = 1e-12;
883 assert_eq!((Affine::rotate(0.5) * c).signed_area(), 0.5);
884 assert!(((Affine::rotate(0.5) * c).signed_area() - 0.5).abs() < epsilon);
885 assert!(((Affine::translate((0.0, 1.0)) * c).signed_area() - 1.0).abs() < epsilon);
886 assert!(((Affine::translate((1.0, 0.0)) * c).signed_area() - 1.0).abs() < epsilon);
887 }
888
889 #[test]
890 fn cubicbez_signed_area() {
891 // y = 1 - x^3
892 let c = CubicBez::new((1.0, 0.0), (2.0 / 3.0, 1.0), (1.0 / 3.0, 1.0), (0.0, 1.0));
893 let epsilon = 1e-12;
894 assert!((c.signed_area() - 0.75).abs() < epsilon);
895 assert!(((Affine::rotate(0.5) * c).signed_area() - 0.75).abs() < epsilon);
896 assert!(((Affine::translate((0.0, 1.0)) * c).signed_area() - 1.25).abs() < epsilon);
897 assert!(((Affine::translate((1.0, 0.0)) * c).signed_area() - 1.25).abs() < epsilon);
898 }
899
900 #[test]
901 fn cubicbez_nearest() {
902 fn verify(result: Nearest, expected: f64) {
903 assert!(
904 (result.t - expected).abs() < 1e-6,
905 "got {result:?} expected {expected}"
906 );
907 }
908 // y = x^3
909 let c = CubicBez::new((0.0, 0.0), (1.0 / 3.0, 0.0), (2.0 / 3.0, 0.0), (1.0, 1.0));
910 verify(c.nearest((0.1, 0.001).into(), 1e-6), 0.1);
911 verify(c.nearest((0.2, 0.008).into(), 1e-6), 0.2);
912 verify(c.nearest((0.3, 0.027).into(), 1e-6), 0.3);
913 verify(c.nearest((0.4, 0.064).into(), 1e-6), 0.4);
914 verify(c.nearest((0.5, 0.125).into(), 1e-6), 0.5);
915 verify(c.nearest((0.6, 0.216).into(), 1e-6), 0.6);
916 verify(c.nearest((0.7, 0.343).into(), 1e-6), 0.7);
917 verify(c.nearest((0.8, 0.512).into(), 1e-6), 0.8);
918 verify(c.nearest((0.9, 0.729).into(), 1e-6), 0.9);
919 verify(c.nearest((1.0, 1.0).into(), 1e-6), 1.0);
920 verify(c.nearest((1.1, 1.1).into(), 1e-6), 1.0);
921 verify(c.nearest((-0.1, 0.0).into(), 1e-6), 0.0);
922 let a = Affine::rotate(0.5);
923 verify((a * c).nearest(a * Point::new(0.1, 0.001), 1e-6), 0.1);
924 }
925
926 // ensure to_quads returns something given colinear points
927 #[test]
928 fn degenerate_to_quads() {
929 let c = CubicBez::new((0., 9.), (6., 6.), (12., 3.0), (18., 0.0));
930 let quads = c.to_quads(1e-6).collect::<Vec<_>>();
931 assert_eq!(quads.len(), 1, "{:?}", &quads);
932 }
933
934 #[test]
935 fn cubicbez_extrema() {
936 // y = x^2
937 let q = CubicBez::new((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0));
938 let extrema = q.extrema();
939 assert_eq!(extrema.len(), 1);
940 assert!((extrema[0] - 0.5).abs() < 1e-6);
941
942 let q = CubicBez::new((0.4, 0.5), (0.0, 1.0), (1.0, 0.0), (0.5, 0.4));
943 let extrema = q.extrema();
944 assert_eq!(extrema.len(), 4);
945 }
946
947 #[test]
948 fn cubicbez_toquads() {
949 // y = x^3
950 let c = CubicBez::new((0.0, 0.0), (1.0 / 3.0, 0.0), (2.0 / 3.0, 0.0), (1.0, 1.0));
951 for i in 0..10 {
952 let accuracy = 0.1f64.powi(i);
953 let mut worst: f64 = 0.0;
954 for (t0, t1, q) in c.to_quads(accuracy) {
955 let epsilon = 1e-12;
956 assert!((q.start() - c.eval(t0)).hypot() < epsilon);
957 assert!((q.end() - c.eval(t1)).hypot() < epsilon);
958 let n = 4;
959 for j in 0..=n {
960 let t = (j as f64) * (n as f64).recip();
961 let p = q.eval(t);
962 let err = (p.y - p.x.powi(3)).abs();
963 worst = worst.max(err);
964 assert!(err < accuracy, "got {err} wanted {accuracy}");
965 }
966 }
967 }
968 }
969
970 #[test]
971 fn cubicbez_approx_spline() {
972 let c1 = CubicBez::new(
973 (550.0, 258.0),
974 (1044.0, 482.0),
975 (2029.0, 1841.0),
976 (1934.0, 1554.0),
977 );
978
979 let quad = c1.try_approx_quadratic(344.0);
980 let expected = QuadBez::new(
981 Point::new(550.0, 258.0),
982 Point::new(1673.665720592873, 767.5164401068898),
983 Point::new(1934.0, 1554.0),
984 );
985 assert!(quad.is_some());
986 assert_eq!(quad.unwrap(), expected);
987
988 let quad = c1.try_approx_quadratic(343.0);
989 assert!(quad.is_none());
990
991 let spline = c1.approx_spline_n(2, 343.0);
992 assert!(spline.is_some());
993 let spline = spline.unwrap();
994 let expected = [
995 Point::new(550.0, 258.0),
996 Point::new(920.5, 426.0),
997 Point::new(2005.25, 1769.25),
998 Point::new(1934.0, 1554.0),
999 ];
1000 assert_eq!(spline.points().len(), expected.len());
1001 for (got, &wanted) in spline.points().iter().zip(expected.iter()) {
1002 assert!(got.distance(wanted) < 5.0);
1003 }
1004
1005 let spline = c1.approx_spline(5.0);
1006 let expected = [
1007 Point::new(550.0, 258.0),
1008 Point::new(673.5, 314.0),
1009 Point::new(984.8777777777776, 584.2666666666667),
1010 Point::new(1312.6305555555557, 927.825),
1011 Point::new(1613.1194444444443, 1267.425),
1012 Point::new(1842.7055555555555, 1525.8166666666666),
1013 Point::new(1957.75, 1625.75),
1014 Point::new(1934.0, 1554.0),
1015 ];
1016 assert!(spline.is_some());
1017 let spline = spline.unwrap();
1018 assert_eq!(spline.points().len(), expected.len());
1019 for (got, &wanted) in spline.points().iter().zip(expected.iter()) {
1020 assert!(got.distance(wanted) < 5.0);
1021 }
1022 }
1023
1024 #[test]
1025 fn cubicbez_cubics_to_quadratic_splines() {
1026 let curves = vec![
1027 CubicBez::new(
1028 (550.0, 258.0),
1029 (1044.0, 482.0),
1030 (2029.0, 1841.0),
1031 (1934.0, 1554.0),
1032 ),
1033 CubicBez::new(
1034 (859.0, 384.0),
1035 (1998.0, 116.0),
1036 (1596.0, 1772.0),
1037 (8.0, 1824.0),
1038 ),
1039 CubicBez::new(
1040 (1090.0, 937.0),
1041 (418.0, 1300.0),
1042 (125.0, 91.0),
1043 (104.0, 37.0),
1044 ),
1045 ];
1046 let converted = cubics_to_quadratic_splines(&curves, 5.0);
1047 assert!(converted.is_some());
1048 let converted = converted.unwrap();
1049 assert_eq!(converted[0].points().len(), 8);
1050 assert_eq!(converted[1].points().len(), 8);
1051 assert_eq!(converted[2].points().len(), 8);
1052 assert!(converted[0].points()[1].distance(Point::new(673.5, 314.0)) < 0.0001);
1053 assert!(
1054 converted[0].points()[2].distance(Point::new(88639.0 / 90.0, 52584.0 / 90.0)) < 0.0001
1055 );
1056 }
1057
1058 #[test]
1059 fn cubicbez_approx_spline_div_exact() {
1060 // Ensure rounding behavior for division matches fonttools
1061 // cu2qu.
1062 // See <https://github.com/linebender/kurbo/issues/272>
1063 let cubic = CubicBez::new(
1064 Point::new(408.0, 321.0),
1065 Point::new(408.0, 452.0),
1066 Point::new(342.0, 560.0),
1067 Point::new(260.0, 560.0),
1068 );
1069 let spline = cubic.approx_spline(1.0).unwrap();
1070 assert_eq!(
1071 spline.points(),
1072 &[
1073 Point::new(408.0, 321.0),
1074 // Previous behavior produced 386.49999999999994 for the
1075 // y coordinate leading to inconsistent rounding.
1076 Point::new(408.0, 386.5),
1077 Point::new(368.16666666666663, 495.0833333333333),
1078 Point::new(301.0, 560.0),
1079 Point::new(260.0, 560.0)
1080 ]
1081 );
1082 }
1083
1084 #[test]
1085 fn cubicbez_inflections() {
1086 let c = CubicBez::new((0., 0.), (0.8, 1.), (0.2, 1.), (1., 0.));
1087 let inflections = c.inflections();
1088 assert_eq!(inflections.len(), 2);
1089 assert!((inflections[0] - 0.311018).abs() < 1e-6);
1090 assert!((inflections[1] - 0.688982).abs() < 1e-6);
1091 let c = CubicBez::new((0., 0.), (1., 1.), (2., -1.), (3., 0.));
1092 let inflections = c.inflections();
1093 assert_eq!(inflections.len(), 1);
1094 assert!((inflections[0] - 0.5).abs() < 1e-6);
1095 let c = CubicBez::new((0., 0.), (1., 1.), (2., 1.), (3., 0.));
1096 let inflections = c.inflections();
1097 assert_eq!(inflections.len(), 0);
1098 }
1099
1100 #[test]
1101 fn cubic_to_quadratic_matches_python() {
1102 // from https://github.com/googlefonts/fontmake-rs/issues/217
1103 let cubic = CubicBez {
1104 p0: (796.0, 319.0).into(),
1105 p1: (727.0, 314.0).into(),
1106 p2: (242.0, 303.0).into(),
1107 p3: (106.0, 303.0).into(),
1108 };
1109
1110 // FontTools can approximate this curve successfully in 7 splits, we can too
1111 assert!(cubic.approx_spline_n(7, 1.0).is_some());
1112
1113 // FontTools can solve this with accuracy 0.001, we can too
1114 assert!(cubics_to_quadratic_splines(&[cubic], 0.001).is_some());
1115 }
1116
1117 #[test]
1118 fn cubics_to_quadratic_splines_matches_python() {
1119 // https://github.com/linebender/kurbo/pull/273
1120 let light = CubicBez::new((378., 608.), (378., 524.), (355., 455.), (266., 455.));
1121 let regular = CubicBez::new((367., 607.), (367., 511.), (338., 472.), (243., 472.));
1122 let bold = CubicBez::new(
1123 (372.425, 593.05),
1124 (372.425, 524.95),
1125 (355.05, 485.95),
1126 (274., 485.95),
1127 );
1128 let qsplines = cubics_to_quadratic_splines(&[light, regular, bold], 1.0).unwrap();
1129 assert_eq!(
1130 qsplines,
1131 [
1132 QuadSpline::new(vec![
1133 (378.0, 608.0).into(),
1134 (378.0, 566.0).into(),
1135 (359.0833333333333, 496.5).into(),
1136 (310.5, 455.0).into(),
1137 (266.0, 455.0).into(),
1138 ]),
1139 QuadSpline::new(vec![
1140 (367.0, 607.0).into(),
1141 (367.0, 559.0).into(),
1142 // Previous behavior produced 496.5 for the y coordinate
1143 (344.5833333333333, 499.49999999999994).into(),
1144 (290.5, 472.0).into(),
1145 (243.0, 472.0).into(),
1146 ]),
1147 QuadSpline::new(vec![
1148 (372.425, 593.05).into(),
1149 (372.425, 559.0).into(),
1150 (356.98333333333335, 511.125).into(),
1151 (314.525, 485.95).into(),
1152 (274.0, 485.95).into(),
1153 ]),
1154 ]
1155 );
1156 }
1157}
1158