kurbo/
cubicbez.rs

Help
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(always)]
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    #[inline(always)]
495    fn area(&self) -> f64 {
496        0.0
497    }
498
499    #[inline]
500    fn perimeter(&self, accuracy: f64) -> f64 {
501        self.arclen(accuracy)
502    }
503
504    #[inline(always)]
505    fn winding(&self, _pt: Point) -> i32 {
506        0
507    }
508
509    #[inline]
510    fn bounding_box(&self) -> Rect {
511        ParamCurveExtrema::bounding_box(self)
512    }
513}
514
515impl Iterator for CubicBezIter {
516    type Item = PathEl;
517
518    fn next(&mut self) -> Option<PathEl> {
519        self.ix += 1;
520        match self.ix {
521            1 => Some(PathEl::MoveTo(self.cubic.p0)),
522            2 => Some(PathEl::CurveTo(self.cubic.p1, self.cubic.p2, self.cubic.p3)),
523            _ => None,
524        }
525    }
526}
527
528impl ParamCurve for CubicBez {
529    #[inline]
530    fn eval(&self, t: f64) -> Point {
531        let mt = 1.0 - t;
532        let v = self.p0.to_vec2() * (mt * mt * mt)
533            + (self.p1.to_vec2() * (mt * mt * 3.0)
534                + (self.p2.to_vec2() * (mt * 3.0) + self.p3.to_vec2() * t) * t)
535                * t;
536        v.to_point()
537    }
538
539    fn subsegment(&self, range: Range<f64>) -> CubicBez {
540        let (t0, t1) = (range.start, range.end);
541        let p0 = self.eval(t0);
542        let p3 = self.eval(t1);
543        let d = self.deriv();
544        let scale = (t1 - t0) * (1.0 / 3.0);
545        let p1 = p0 + scale * d.eval(t0).to_vec2();
546        let p2 = p3 - scale * d.eval(t1).to_vec2();
547        CubicBez { p0, p1, p2, p3 }
548    }
549
550    /// Subdivide into halves, using de Casteljau.
551    #[inline]
552    fn subdivide(&self) -> (CubicBez, CubicBez) {
553        let pm = self.eval(0.5);
554        (
555            CubicBez::new(
556                self.p0,
557                self.p0.midpoint(self.p1),
558                ((self.p0.to_vec2() + self.p1.to_vec2() * 2.0 + self.p2.to_vec2()) * 0.25)
559                    .to_point(),
560                pm,
561            ),
562            CubicBez::new(
563                pm,
564                ((self.p1.to_vec2() + self.p2.to_vec2() * 2.0 + self.p3.to_vec2()) * 0.25)
565                    .to_point(),
566                self.p2.midpoint(self.p3),
567                self.p3,
568            ),
569        )
570    }
571
572    #[inline(always)]
573    fn start(&self) -> Point {
574        self.p0
575    }
576
577    #[inline(always)]
578    fn end(&self) -> Point {
579        self.p3
580    }
581}
582
583impl ParamCurveDeriv for CubicBez {
584    type DerivResult = QuadBez;
585
586    #[inline]
587    fn deriv(&self) -> QuadBez {
588        QuadBez::new(
589            (3.0 * (self.p1 - self.p0)).to_point(),
590            (3.0 * (self.p2 - self.p1)).to_point(),
591            (3.0 * (self.p3 - self.p2)).to_point(),
592        )
593    }
594}
595
596fn arclen_quadrature_core(coeffs: &[(f64, f64)], dm: Vec2, dm1: Vec2, dm2: Vec2) -> f64 {
597    coeffs
598        .iter()
599        .map(|&(wi, xi)| {
600            let d = dm + dm2 * (xi * xi);
601            let dpx = (d + dm1 * xi).hypot();
602            let dmx = (d - dm1 * xi).hypot();
603            (2.25f64.sqrt() * wi) * (dpx + dmx)
604        })
605        .sum::<f64>()
606}
607
608fn arclen_rec(c: &CubicBez, accuracy: f64, depth: usize) -> f64 {
609    let d03 = c.p3 - c.p0;
610    let d01 = c.p1 - c.p0;
611    let d12 = c.p2 - c.p1;
612    let d23 = c.p3 - c.p2;
613    let lp_lc = d01.hypot() + d12.hypot() + d23.hypot() - d03.hypot();
614    let dd1 = d12 - d01;
615    let dd2 = d23 - d12;
616    // It might be faster to do direct multiplies, the data dependencies would be shorter.
617    // The following values don't have the factor of 3 for first deriv
618    let dm = 0.25 * (d01 + d23) + 0.5 * d12; // first derivative at midpoint
619    let dm1 = 0.5 * (dd2 + dd1); // second derivative at midpoint
620    let dm2 = 0.25 * (dd2 - dd1); // 0.5 * (third derivative at midpoint)
621
622    let est = GAUSS_LEGENDRE_COEFFS_8
623        .iter()
624        .map(|&(wi, xi)| {
625            wi * {
626                let d_norm2 = (dm + dm1 * xi + dm2 * (xi * xi)).hypot2();
627                let dd_norm2 = (dm1 + dm2 * (2.0 * xi)).hypot2();
628                dd_norm2 / d_norm2
629            }
630        })
631        .sum::<f64>();
632    let est_gauss8_error = (est.powi(3) * 2.5e-6).min(3e-2) * lp_lc;
633    if est_gauss8_error < accuracy {
634        return arclen_quadrature_core(GAUSS_LEGENDRE_COEFFS_8_HALF, dm, dm1, dm2);
635    }
636    let est_gauss16_error = (est.powi(6) * 1.5e-11).min(9e-3) * lp_lc;
637    if est_gauss16_error < accuracy {
638        return arclen_quadrature_core(GAUSS_LEGENDRE_COEFFS_16_HALF, dm, dm1, dm2);
639    }
640    let est_gauss24_error = (est.powi(9) * 3.5e-16).min(3.5e-3) * lp_lc;
641    if est_gauss24_error < accuracy || depth >= 20 {
642        return arclen_quadrature_core(GAUSS_LEGENDRE_COEFFS_24_HALF, dm, dm1, dm2);
643    }
644    let (c0, c1) = c.subdivide();
645    arclen_rec(&c0, accuracy * 0.5, depth + 1) + arclen_rec(&c1, accuracy * 0.5, depth + 1)
646}
647
648impl ParamCurveArclen for CubicBez {
649    /// Arclength of a cubic Bézier segment.
650    ///
651    /// This is an adaptive subdivision approach using Legendre-Gauss quadrature
652    /// in the base case, and an error estimate to decide when to subdivide.
653    fn arclen(&self, accuracy: f64) -> f64 {
654        arclen_rec(self, accuracy, 0)
655    }
656}
657
658impl ParamCurveArea for CubicBez {
659    #[inline]
660    fn signed_area(&self) -> f64 {
661        (self.p0.x * (6.0 * self.p1.y + 3.0 * self.p2.y + self.p3.y)
662            + 3.0
663                * (self.p1.x * (-2.0 * self.p0.y + self.p2.y + self.p3.y)
664                    - self.p2.x * (self.p0.y + self.p1.y - 2.0 * self.p3.y))
665            - self.p3.x * (self.p0.y + 3.0 * self.p1.y + 6.0 * self.p2.y))
666            * (1.0 / 20.0)
667    }
668}
669
670impl ParamCurveNearest for CubicBez {
671    /// Find the nearest point, using subdivision.
672    fn nearest(&self, p: Point, accuracy: f64) -> Nearest {
673        let mut best_r = None;
674        let mut best_t = 0.0;
675        for (t0, t1, q) in self.to_quads(accuracy) {
676            let nearest = q.nearest(p, accuracy);
677            if best_r
678                .map(|best_r| nearest.distance_sq < best_r)
679                .unwrap_or(true)
680            {
681                best_t = t0 + nearest.t * (t1 - t0);
682                best_r = Some(nearest.distance_sq);
683            }
684        }
685        Nearest {
686            t: best_t,
687            distance_sq: best_r.unwrap(),
688        }
689    }
690}
691
692impl ParamCurveCurvature for CubicBez {}
693
694impl ParamCurveExtrema for CubicBez {
695    fn extrema(&self) -> ArrayVec<f64, MAX_EXTREMA> {
696        fn one_coord(result: &mut ArrayVec<f64, MAX_EXTREMA>, d0: f64, d1: f64, d2: f64) {
697            let a = d0 - 2.0 * d1 + d2;
698            let b = 2.0 * (d1 - d0);
699            let c = d0;
700            let roots = solve_quadratic(c, b, a);
701            for &t in &roots {
702                if t > 0.0 && t < 1.0 {
703                    result.push(t);
704                }
705            }
706        }
707        let mut result = ArrayVec::new();
708        let d0 = self.p1 - self.p0;
709        let d1 = self.p2 - self.p1;
710        let d2 = self.p3 - self.p2;
711        one_coord(&mut result, d0.x, d1.x, d2.x);
712        one_coord(&mut result, d0.y, d1.y, d2.y);
713        result.sort_by(|a, b| a.partial_cmp(b).unwrap());
714        result
715    }
716}
717
718impl Mul<CubicBez> for Affine {
719    type Output = CubicBez;
720
721    #[inline]
722    fn mul(self, c: CubicBez) -> CubicBez {
723        CubicBez {
724            p0: self * c.p0,
725            p1: self * c.p1,
726            p2: self * c.p2,
727            p3: self * c.p3,
728        }
729    }
730}
731
732impl Iterator for ToQuads {
733    type Item = (f64, f64, QuadBez);
734
735    fn next(&mut self) -> Option<(f64, f64, QuadBez)> {
736        if self.i == self.n {
737            return None;
738        }
739        let t0 = self.i as f64 / self.n as f64;
740        let t1 = (self.i + 1) as f64 / self.n as f64;
741        let seg = self.c.subsegment(t0..t1);
742        let p1x2 = 3.0 * seg.p1.to_vec2() - seg.p0.to_vec2();
743        let p2x2 = 3.0 * seg.p2.to_vec2() - seg.p3.to_vec2();
744        let result = QuadBez::new(seg.p0, ((p1x2 + p2x2) / 4.0).to_point(), seg.p3);
745        self.i += 1;
746        Some((t0, t1, result))
747    }
748
749    fn size_hint(&self) -> (usize, Option<usize>) {
750        let remaining = self.n - self.i;
751        (remaining, Some(remaining))
752    }
753}
754
755/// Convert multiple cubic Bézier curves to quadratic splines.
756///
757/// Ensures that the resulting splines have the same number of control points.
758///
759/// Rust port of cu2qu [cubic_approx_quadratic](https://github.com/fonttools/fonttools/blob/3b9a73ff8379ab49d3ce35aaaaf04b3a7d9d1655/Lib/fontTools/cu2qu/cu2qu.py#L322)
760pub fn cubics_to_quadratic_splines(curves: &[CubicBez], accuracy: f64) -> Option<Vec<QuadSpline>> {
761    let mut result = Vec::new();
762    let mut split_order = 0;
763
764    while split_order <= MAX_SPLINE_SPLIT {
765        split_order += 1;
766        result.clear();
767
768        for curve in curves {
769            match curve.approx_spline_n(split_order, accuracy) {
770                Some(spline) => result.push(spline),
771                None => break,
772            }
773        }
774
775        if result.len() == curves.len() {
776            return Some(result);
777        }
778    }
779    None
780}
781#[cfg(test)]
782mod tests {
783    use crate::{
784        cubics_to_quadratic_splines, Affine, CubicBez, Nearest, ParamCurve, ParamCurveArclen,
785        ParamCurveArea, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, Point, QuadBez,
786        QuadSpline,
787    };
788
789    #[test]
790    fn cubicbez_deriv() {
791        // y = x^2
792        let c = CubicBez::new(
793            (0.0, 0.0),
794            (1.0 / 3.0, 0.0),
795            (2.0 / 3.0, 1.0 / 3.0),
796            (1.0, 1.0),
797        );
798        let deriv = c.deriv();
799
800        let n = 10;
801        for i in 0..=n {
802            let t = (i as f64) * (n as f64).recip();
803            let delta = 1e-6;
804            let p = c.eval(t);
805            let p1 = c.eval(t + delta);
806            let d_approx = (p1 - p) * delta.recip();
807            let d = deriv.eval(t).to_vec2();
808            assert!((d - d_approx).hypot() < delta * 2.0);
809        }
810    }
811
812    #[test]
813    fn cubicbez_arclen() {
814        // y = x^2
815        let c = CubicBez::new(
816            (0.0, 0.0),
817            (1.0 / 3.0, 0.0),
818            (2.0 / 3.0, 1.0 / 3.0),
819            (1.0, 1.0),
820        );
821        let true_arclen = 0.5 * 5.0f64.sqrt() + 0.25 * (2.0 + 5.0f64.sqrt()).ln();
822        for i in 0..12 {
823            let accuracy = 0.1f64.powi(i);
824            let error = c.arclen(accuracy) - true_arclen;
825            assert!(error.abs() < accuracy);
826        }
827    }
828
829    #[test]
830    fn cubicbez_inv_arclen() {
831        // y = x^2 / 100
832        let c = CubicBez::new(
833            (0.0, 0.0),
834            (100.0 / 3.0, 0.0),
835            (200.0 / 3.0, 100.0 / 3.0),
836            (100.0, 100.0),
837        );
838        let true_arclen = 100.0 * (0.5 * 5.0f64.sqrt() + 0.25 * (2.0 + 5.0f64.sqrt()).ln());
839        for i in 0..12 {
840            let accuracy = 0.1f64.powi(i);
841            let n = 10;
842            for j in 0..=n {
843                let arc = (j as f64) * ((n as f64).recip() * true_arclen);
844                let t = c.inv_arclen(arc, accuracy * 0.5);
845                let actual_arc = c.subsegment(0.0..t).arclen(accuracy * 0.5);
846                assert!(
847                    (arc - actual_arc).abs() < accuracy,
848                    "at accuracy {accuracy:e}, wanted {actual_arc} got {arc}"
849                );
850            }
851        }
852        // corner case: user passes accuracy larger than total arc length
853        let accuracy = true_arclen * 1.1;
854        let arc = true_arclen * 0.5;
855        let t = c.inv_arclen(arc, accuracy);
856        let actual_arc = c.subsegment(0.0..t).arclen(accuracy);
857        assert!(
858            (arc - actual_arc).abs() < 2.0 * accuracy,
859            "at accuracy {accuracy:e}, want {actual_arc} got {arc}"
860        );
861    }
862
863    #[test]
864    fn cubicbez_inv_arclen_accuracy() {
865        let c = CubicBez::new((0.2, 0.73), (0.35, 1.08), (0.85, 1.08), (1.0, 0.73));
866        let true_t = c.inv_arclen(0.5, 1e-12);
867        for i in 1..12 {
868            let accuracy = (0.1f64).powi(i);
869            let approx_t = c.inv_arclen(0.5, accuracy);
870            assert!((approx_t - true_t).abs() <= accuracy);
871        }
872    }
873
874    #[test]
875    #[allow(clippy::float_cmp)]
876    fn cubicbez_signed_area_linear() {
877        // y = 1 - x
878        let c = CubicBez::new(
879            (1.0, 0.0),
880            (2.0 / 3.0, 1.0 / 3.0),
881            (1.0 / 3.0, 2.0 / 3.0),
882            (0.0, 1.0),
883        );
884        let epsilon = 1e-12;
885        assert_eq!((Affine::rotate(0.5) * c).signed_area(), 0.5);
886        assert!(((Affine::rotate(0.5) * c).signed_area() - 0.5).abs() < epsilon);
887        assert!(((Affine::translate((0.0, 1.0)) * c).signed_area() - 1.0).abs() < epsilon);
888        assert!(((Affine::translate((1.0, 0.0)) * c).signed_area() - 1.0).abs() < epsilon);
889    }
890
891    #[test]
892    fn cubicbez_signed_area() {
893        // y = 1 - x^3
894        let c = CubicBez::new((1.0, 0.0), (2.0 / 3.0, 1.0), (1.0 / 3.0, 1.0), (0.0, 1.0));
895        let epsilon = 1e-12;
896        assert!((c.signed_area() - 0.75).abs() < epsilon);
897        assert!(((Affine::rotate(0.5) * c).signed_area() - 0.75).abs() < epsilon);
898        assert!(((Affine::translate((0.0, 1.0)) * c).signed_area() - 1.25).abs() < epsilon);
899        assert!(((Affine::translate((1.0, 0.0)) * c).signed_area() - 1.25).abs() < epsilon);
900    }
901
902    #[test]
903    fn cubicbez_nearest() {
904        fn verify(result: Nearest, expected: f64) {
905            assert!(
906                (result.t - expected).abs() < 1e-6,
907                "got {result:?} expected {expected}"
908            );
909        }
910        // y = x^3
911        let c = CubicBez::new((0.0, 0.0), (1.0 / 3.0, 0.0), (2.0 / 3.0, 0.0), (1.0, 1.0));
912        verify(c.nearest((0.1, 0.001).into(), 1e-6), 0.1);
913        verify(c.nearest((0.2, 0.008).into(), 1e-6), 0.2);
914        verify(c.nearest((0.3, 0.027).into(), 1e-6), 0.3);
915        verify(c.nearest((0.4, 0.064).into(), 1e-6), 0.4);
916        verify(c.nearest((0.5, 0.125).into(), 1e-6), 0.5);
917        verify(c.nearest((0.6, 0.216).into(), 1e-6), 0.6);
918        verify(c.nearest((0.7, 0.343).into(), 1e-6), 0.7);
919        verify(c.nearest((0.8, 0.512).into(), 1e-6), 0.8);
920        verify(c.nearest((0.9, 0.729).into(), 1e-6), 0.9);
921        verify(c.nearest((1.0, 1.0).into(), 1e-6), 1.0);
922        verify(c.nearest((1.1, 1.1).into(), 1e-6), 1.0);
923        verify(c.nearest((-0.1, 0.0).into(), 1e-6), 0.0);
924        let a = Affine::rotate(0.5);
925        verify((a * c).nearest(a * Point::new(0.1, 0.001), 1e-6), 0.1);
926    }
927
928    // ensure to_quads returns something given collinear points
929    #[test]
930    fn degenerate_to_quads() {
931        let c = CubicBez::new((0., 9.), (6., 6.), (12., 3.0), (18., 0.0));
932        let quads = c.to_quads(1e-6).collect::<Vec<_>>();
933        assert_eq!(quads.len(), 1, "{:?}", &quads);
934    }
935
936    #[test]
937    fn cubicbez_extrema() {
938        // y = x^2
939        let q = CubicBez::new((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0));
940        let extrema = q.extrema();
941        assert_eq!(extrema.len(), 1);
942        assert!((extrema[0] - 0.5).abs() < 1e-6);
943
944        let q = CubicBez::new((0.4, 0.5), (0.0, 1.0), (1.0, 0.0), (0.5, 0.4));
945        let extrema = q.extrema();
946        assert_eq!(extrema.len(), 4);
947    }
948
949    #[test]
950    fn cubicbez_toquads() {
951        // y = x^3
952        let c = CubicBez::new((0.0, 0.0), (1.0 / 3.0, 0.0), (2.0 / 3.0, 0.0), (1.0, 1.0));
953        for i in 0..10 {
954            let accuracy = 0.1f64.powi(i);
955            let mut worst: f64 = 0.0;
956            for (t0, t1, q) in c.to_quads(accuracy) {
957                let epsilon = 1e-12;
958                assert!((q.start() - c.eval(t0)).hypot() < epsilon);
959                assert!((q.end() - c.eval(t1)).hypot() < epsilon);
960                let n = 4;
961                for j in 0..=n {
962                    let t = (j as f64) * (n as f64).recip();
963                    let p = q.eval(t);
964                    let err = (p.y - p.x.powi(3)).abs();
965                    worst = worst.max(err);
966                    assert!(err < accuracy, "got {err} wanted {accuracy}");
967                }
968            }
969        }
970    }
971
972    #[test]
973    fn cubicbez_approx_spline() {
974        let c1 = CubicBez::new(
975            (550.0, 258.0),
976            (1044.0, 482.0),
977            (2029.0, 1841.0),
978            (1934.0, 1554.0),
979        );
980
981        let quad = c1.try_approx_quadratic(344.0);
982        let expected = QuadBez::new(
983            Point::new(550.0, 258.0),
984            Point::new(1673.665720592873, 767.5164401068898),
985            Point::new(1934.0, 1554.0),
986        );
987        assert!(quad.is_some());
988        assert_eq!(quad.unwrap(), expected);
989
990        let quad = c1.try_approx_quadratic(343.0);
991        assert!(quad.is_none());
992
993        let spline = c1.approx_spline_n(2, 343.0);
994        assert!(spline.is_some());
995        let spline = spline.unwrap();
996        let expected = [
997            Point::new(550.0, 258.0),
998            Point::new(920.5, 426.0),
999            Point::new(2005.25, 1769.25),
1000            Point::new(1934.0, 1554.0),
1001        ];
1002        assert_eq!(spline.points().len(), expected.len());
1003        for (got, &wanted) in spline.points().iter().zip(expected.iter()) {
1004            assert!(got.distance(wanted) < 5.0);
1005        }
1006
1007        let spline = c1.approx_spline(5.0);
1008        let expected = [
1009            Point::new(550.0, 258.0),
1010            Point::new(673.5, 314.0),
1011            Point::new(984.8777777777776, 584.2666666666667),
1012            Point::new(1312.6305555555557, 927.825),
1013            Point::new(1613.1194444444443, 1267.425),
1014            Point::new(1842.7055555555555, 1525.8166666666666),
1015            Point::new(1957.75, 1625.75),
1016            Point::new(1934.0, 1554.0),
1017        ];
1018        assert!(spline.is_some());
1019        let spline = spline.unwrap();
1020        assert_eq!(spline.points().len(), expected.len());
1021        for (got, &wanted) in spline.points().iter().zip(expected.iter()) {
1022            assert!(got.distance(wanted) < 5.0);
1023        }
1024    }
1025
1026    #[test]
1027    fn cubicbez_cubics_to_quadratic_splines() {
1028        let curves = vec![
1029            CubicBez::new(
1030                (550.0, 258.0),
1031                (1044.0, 482.0),
1032                (2029.0, 1841.0),
1033                (1934.0, 1554.0),
1034            ),
1035            CubicBez::new(
1036                (859.0, 384.0),
1037                (1998.0, 116.0),
1038                (1596.0, 1772.0),
1039                (8.0, 1824.0),
1040            ),
1041            CubicBez::new(
1042                (1090.0, 937.0),
1043                (418.0, 1300.0),
1044                (125.0, 91.0),
1045                (104.0, 37.0),
1046            ),
1047        ];
1048        let converted = cubics_to_quadratic_splines(&curves, 5.0);
1049        assert!(converted.is_some());
1050        let converted = converted.unwrap();
1051        assert_eq!(converted[0].points().len(), 8);
1052        assert_eq!(converted[1].points().len(), 8);
1053        assert_eq!(converted[2].points().len(), 8);
1054        assert!(converted[0].points()[1].distance(Point::new(673.5, 314.0)) < 0.0001);
1055        assert!(
1056            converted[0].points()[2].distance(Point::new(88639.0 / 90.0, 52584.0 / 90.0)) < 0.0001
1057        );
1058    }
1059
1060    #[test]
1061    fn cubicbez_approx_spline_div_exact() {
1062        // Ensure rounding behavior for division matches fonttools
1063        // cu2qu.
1064        // See <https://github.com/linebender/kurbo/issues/272>
1065        let cubic = CubicBez::new(
1066            Point::new(408.0, 321.0),
1067            Point::new(408.0, 452.0),
1068            Point::new(342.0, 560.0),
1069            Point::new(260.0, 560.0),
1070        );
1071        let spline = cubic.approx_spline(1.0).unwrap();
1072        assert_eq!(
1073            spline.points(),
1074            &[
1075                Point::new(408.0, 321.0),
1076                // Previous behavior produced 386.49999999999994 for the
1077                // y coordinate leading to inconsistent rounding.
1078                Point::new(408.0, 386.5),
1079                Point::new(368.16666666666663, 495.0833333333333),
1080                Point::new(301.0, 560.0),
1081                Point::new(260.0, 560.0)
1082            ]
1083        );
1084    }
1085
1086    #[test]
1087    fn cubicbez_inflections() {
1088        let c = CubicBez::new((0., 0.), (0.8, 1.), (0.2, 1.), (1., 0.));
1089        let inflections = c.inflections();
1090        assert_eq!(inflections.len(), 2);
1091        assert!((inflections[0] - 0.311018).abs() < 1e-6);
1092        assert!((inflections[1] - 0.688982).abs() < 1e-6);
1093        let c = CubicBez::new((0., 0.), (1., 1.), (2., -1.), (3., 0.));
1094        let inflections = c.inflections();
1095        assert_eq!(inflections.len(), 1);
1096        assert!((inflections[0] - 0.5).abs() < 1e-6);
1097        let c = CubicBez::new((0., 0.), (1., 1.), (2., 1.), (3., 0.));
1098        let inflections = c.inflections();
1099        assert_eq!(inflections.len(), 0);
1100    }
1101
1102    #[test]
1103    fn cubic_to_quadratic_matches_python() {
1104        // from https://github.com/googlefonts/fontmake-rs/issues/217
1105        let cubic = CubicBez {
1106            p0: (796.0, 319.0).into(),
1107            p1: (727.0, 314.0).into(),
1108            p2: (242.0, 303.0).into(),
1109            p3: (106.0, 303.0).into(),
1110        };
1111
1112        // FontTools can approximate this curve successfully in 7 splits, we can too
1113        assert!(cubic.approx_spline_n(7, 1.0).is_some());
1114
1115        // FontTools can solve this with accuracy 0.001, we can too
1116        assert!(cubics_to_quadratic_splines(&[cubic], 0.001).is_some());
1117    }
1118
1119    #[test]
1120    fn cubics_to_quadratic_splines_matches_python() {
1121        // https://github.com/linebender/kurbo/pull/273
1122        let light = CubicBez::new((378., 608.), (378., 524.), (355., 455.), (266., 455.));
1123        let regular = CubicBez::new((367., 607.), (367., 511.), (338., 472.), (243., 472.));
1124        let bold = CubicBez::new(
1125            (372.425, 593.05),
1126            (372.425, 524.95),
1127            (355.05, 485.95),
1128            (274., 485.95),
1129        );
1130        let qsplines = cubics_to_quadratic_splines(&[light, regular, bold], 1.0).unwrap();
1131        assert_eq!(
1132            qsplines,
1133            [
1134                QuadSpline::new(vec![
1135                    (378.0, 608.0).into(),
1136                    (378.0, 566.0).into(),
1137                    (359.0833333333333, 496.5).into(),
1138                    (310.5, 455.0).into(),
1139                    (266.0, 455.0).into(),
1140                ]),
1141                QuadSpline::new(vec![
1142                    (367.0, 607.0).into(),
1143                    (367.0, 559.0).into(),
1144                    // Previous behavior produced 496.5 for the y coordinate
1145                    (344.5833333333333, 499.49999999999994).into(),
1146                    (290.5, 472.0).into(),
1147                    (243.0, 472.0).into(),
1148                ]),
1149                QuadSpline::new(vec![
1150                    (372.425, 593.05).into(),
1151                    (372.425, 559.0).into(),
1152                    (356.98333333333335, 511.125).into(),
1153                    (314.525, 485.95).into(),
1154                    (274.0, 485.95).into(),
1155                ]),
1156            ]
1157        );
1158    }
1159}
kurbo/cubicbez.rs

kurbo/
cubicbez.rs
