rasterize/
curve.rs

1//! All the things you need to handle bezier curves
2
3use crate::{
4    BBox, EPSILON, EllipArc, LineCap, LineJoin, Point, Scalar, StrokeStyle, SvgParserError,
5    SvgPathCmd, SvgPathParser, Transform,
6    utils::{
7        ArrayIter, M3x3, M4x4, QUADRATURE_16, QUADRATURE_32, clamp, cubic_solve,
8        integrate_quadrature, quadratic_solve,
9    },
10};
11use std::{fmt, io::Cursor, str::FromStr};
12
13/// Iterator containing curve roots
14pub type CurveRoots = ArrayIter<Scalar, 3>;
15/// Iterator containing curve extremities
16pub type CurveExtremities = ArrayIter<Scalar, 6>;
17
18/// Set of operations common to all bezier curves.
19pub trait Curve: Into<Segment> {
20    /// Convert curve to an iterator over line segments with desired flatness
21    fn flatten(&self, tr: Transform, flatness: Scalar) -> CurveFlattenIter {
22        CurveFlattenIter::new(self.transform(tr), flatness)
23    }
24
25    /// Correspond to maximum deviation of the curve from the straight line
26    /// `f = max |curve(t) - line(curve_start, curve_end)(t)|`. This function
27    /// actually returns `16.0 * f^2` to avoid unneeded division and square root.
28    fn flatness(&self) -> Scalar;
29
30    /// Apply affine transformation to the curve
31    fn transform(&self, tr: Transform) -> Self;
32
33    /// Point at which curve starts
34    fn start(&self) -> Point;
35
36    /// Point at which curve ends
37    fn end(&self) -> Point;
38
39    /// Evaluate curve at parameter value `t` in (0.0..=1.0)
40    fn at(&self, t: Scalar) -> Point;
41
42    /// Optimized version of `Curve::split_at(0.5)`
43    fn split(&self) -> (Self, Self) {
44        self.split_at(0.5)
45    }
46
47    /// Split the curve at parameter value `t`
48    fn split_at(&self, t: Scalar) -> (Self, Self);
49
50    /// Create sub-curve specified starting at parameter value `a` and ending at value `b`
51    fn cut(&self, a: Scalar, b: Scalar) -> Self;
52
53    /// Extend provided `init` bounding box with the bounding box of the curve
54    fn bbox(&self, init: Option<BBox>) -> BBox;
55
56    /// Offset the curve by distance `dist`, result is inserted into `out` container
57    fn offset(&self, dist: Scalar, out: &mut impl Extend<Segment>);
58
59    /// Derivative with respect to t, `deriv(t) = [curve'(t)_x, curve'(t)_y]`
60    fn deriv(&self) -> Segment;
61
62    /// Identical curve but directed from end to start, instead of start to end.
63    fn reverse(&self) -> Self;
64
65    /// Find roots of the equation `curve(t)_y = 0`. Values of the parameter at which curve
66    /// crosses y axis.
67    fn roots(&self) -> CurveRoots;
68
69    /// Find all extremities of the curve `curve'(t)_x = 0 || curve'(t)_y = 0`
70    fn extremities(&self) -> CurveExtremities;
71
72    /// Calculate length of the curve from `t0` to `t1`
73    fn length(&self, t0: Scalar, t1: Scalar) -> Scalar;
74
75    /// Find value of parameter `t` given desired `l` length of the segment
76    ///
77    /// This method is not particularly fast, parameter value is found by solving
78    /// `f(t) = self.length(0.0, t) - l == 0` using Newton's method with fallback
79    /// to bisection if next iteration will produce out of bound value.
80    ///
81    /// Reference: <https://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf>
82    fn param_at_length(&self, l: Scalar, error: Option<Scalar>) -> Scalar {
83        let deriv = self.deriv();
84        let length = self.length(0.0, 1.0);
85        let error = error.unwrap_or(length / 1e4);
86        let l = clamp(l, 0.0, length);
87
88        let f = |t| self.length(0.0, t) - l;
89        let f_deriv = |t| deriv.at(t).length();
90
91        let mut t = l / length; // initial guess
92        let mut low = 0.0;
93        let mut high = 1.0;
94        for _ in 0..16 {
95            let ft = f(t);
96            if ft.abs() < error {
97                break;
98            }
99            let t_next = t - ft / f_deriv(t);
100            if ft > 0.0 {
101                high = t;
102                t = if t_next <= low {
103                    (low + high) / 2.0
104                } else {
105                    t_next
106                }
107            } else {
108                low = t;
109                t = if t_next >= high {
110                    (low + high) / 2.0
111                } else {
112                    t_next
113                }
114            }
115        }
116        t
117    }
118}
119
120/// Iterator over line segments approximating curve segment
121pub struct CurveFlattenIter {
122    flatness: Scalar,
123    stack: Vec<Segment>,
124}
125
126impl CurveFlattenIter {
127    pub fn new(segment: impl Into<Segment>, flatness: Scalar) -> Self {
128        Self {
129            flatness: 16.0 * flatness * flatness,
130            stack: vec![segment.into()],
131        }
132    }
133}
134
135impl Iterator for CurveFlattenIter {
136    type Item = Line;
137
138    fn next(&mut self) -> Option<Self::Item> {
139        loop {
140            match self.stack.pop() {
141                None => {
142                    return None;
143                }
144                Some(segment) => {
145                    if segment.flatness() < self.flatness {
146                        return Some(Line([segment.start(), segment.end()]));
147                    }
148                    let (s0, s1) = segment.split();
149                    self.stack.push(s1);
150                    self.stack.push(s0);
151                }
152            }
153        }
154    }
155}
156
157// -----------------------------------------------------------------------------
158// Line
159// -----------------------------------------------------------------------------
160
161/// Line segment curve
162#[derive(Clone, Copy, PartialEq)]
163pub struct Line(pub [Point; 2]);
164
165impl fmt::Debug for Line {
166    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
167        let Line([p0, p1]) = self;
168        write!(f, "Line {:?} {:?}", p0, p1)
169    }
170}
171
172impl fmt::Display for Line {
173    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
174        let Line([p0, p1]) = self;
175        write!(f, "M{:?} L{:?}", p0, p1)
176    }
177}
178
179impl Line {
180    pub fn new(p0: impl Into<Point>, p1: impl Into<Point>) -> Self {
181        Self([p0.into(), p1.into()])
182    }
183
184    /// Length of the line
185    pub fn length(&self) -> Scalar {
186        let Self([p0, p1]) = self;
187        p0.dist(*p1)
188    }
189
190    /// Start and end points of the line
191    pub fn points(&self) -> [Point; 2] {
192        self.0
193    }
194
195    pub fn ends(&self) -> (Line, Line) {
196        (*self, *self)
197    }
198
199    /// Find intersection of two lines
200    ///
201    /// Returns pair of `t` parameters for this line and the other line.
202    /// Found by solving `self.at(t0) == other.at(t1)`. Actual intersection of
203    /// line segments can be found by making sure that `0.0 <= t0 <= 1.0 && 0.0 <= t1 <= 1.0`
204    pub fn intersect(&self, other: Line) -> Option<(Scalar, Scalar)> {
205        let Line([Point([x1, y1]), Point([x2, y2])]) = *self;
206        let Line([Point([x3, y3]), Point([x4, y4])]) = other;
207        let det = (x4 - x3) * (y1 - y2) - (x1 - x2) * (y4 - y3);
208        if det.abs() < EPSILON {
209            return None;
210        }
211        let t0 = ((y3 - y4) * (x1 - x3) + (x4 - x3) * (y1 - y3)) / det;
212        let t1 = ((y1 - y2) * (x1 - x3) + (x2 - x1) * (y1 - y3)) / det;
213        Some((t0, t1))
214    }
215
216    /// Find intersection point between two line segments
217    pub fn intersect_point(&self, other: Line) -> Option<Point> {
218        let (t0, t1) = self.intersect(other)?;
219        if (0.0..=1.0).contains(&t0) && (0.0..=1.0).contains(&t1) {
220            Some(self.at(t0))
221        } else {
222            None
223        }
224    }
225
226    /// Direction vector associated with the line segment
227    pub fn direction(&self) -> Point {
228        self.end() - self.start()
229    }
230}
231
232impl Curve for Line {
233    fn flatness(&self) -> Scalar {
234        0.0
235    }
236
237    fn transform(&self, tr: Transform) -> Self {
238        let Line([p0, p1]) = self;
239        Self([tr.apply(*p0), tr.apply(*p1)])
240    }
241
242    fn start(&self) -> Point {
243        self.0[0]
244    }
245
246    fn end(&self) -> Point {
247        self.0[1]
248    }
249
250    fn at(&self, t: Scalar) -> Point {
251        let Self([p0, p1]) = self;
252        (1.0 - t) * p0 + t * p1
253    }
254
255    fn deriv(&self) -> Segment {
256        let deriv = self.end() - self.start();
257        Line::new(deriv, deriv).into()
258    }
259
260    fn split_at(&self, t: Scalar) -> (Self, Self) {
261        let Self([p0, p1]) = self;
262        let mid = self.at(t);
263        (Self([*p0, mid]), Self([mid, *p1]))
264    }
265
266    fn cut(&self, a: Scalar, b: Scalar) -> Self {
267        Self([self.at(a), self.at(b)])
268    }
269
270    fn bbox(&self, init: Option<BBox>) -> BBox {
271        let Self([p0, p1]) = *self;
272        BBox::new(p0, p1).union_opt(init)
273    }
274
275    fn offset(&self, dist: Scalar, out: &mut impl Extend<Segment>) {
276        out.extend(line_offset(*self, dist).map(Segment::from));
277    }
278
279    fn reverse(&self) -> Self {
280        let Self([p0, p1]) = *self;
281        Self([p1, p0])
282    }
283
284    fn roots(&self) -> CurveRoots {
285        let mut result = CurveRoots::new();
286        let Self([Point([_, y0]), Point([_, y1])]) = self;
287        if (y0 - y1).abs() > EPSILON {
288            let t = y0 / (y0 - y1);
289            if (0.0..=1.0).contains(&t) {
290                result.push(t);
291            }
292        }
293        result
294    }
295
296    fn extremities(&self) -> CurveExtremities {
297        CurveExtremities::new()
298    }
299
300    fn length(&self, t0: Scalar, t1: Scalar) -> Scalar {
301        self.at(t0).dist(self.at(t1))
302    }
303}
304
305impl FromStr for Line {
306    type Err = SvgParserError;
307
308    fn from_str(text: &str) -> Result<Self, Self::Err> {
309        let segment = Segment::from_str(text)?;
310        segment
311            .to_line()
312            .ok_or(SvgParserError::UnexpectedSegmentType)
313    }
314}
315
316// -----------------------------------------------------------------------------
317// Quadratic bezier curve
318// -----------------------------------------------------------------------------
319
320// Matrix form for quadratic bezier curve
321#[rustfmt::skip]
322const Q: M3x3 = M3x3([
323    1.0,  0.0, 0.0,
324   -2.0,  2.0, 0.0,
325    1.0, -2.0, 1.0,
326]);
327
328// Inverted matrix form for quadratic bezier curve
329#[rustfmt::skip]
330const QI: M3x3 = M3x3([
331    1.0, 0.0, 0.0,
332    1.0, 0.5, 0.0,
333    1.0, 1.0, 1.0,
334]);
335
336/// Quadratic bezier curve
337///
338/// Polynomial form:
339/// `(1 - t) ^ 2 * p0 + 2 * (1 - t) * t * p1 + t ^ 2 * p2`
340/// Matrix from:
341/// ```text
342///             ┌          ┐ ┌    ┐
343/// ┌         ┐ │  1  0  0 │ │ p0 │
344/// │ 1 t t^2 │ │ -2  2  0 │ │ p1 │
345/// └         ┘ │  1 -2  1 │ │ p2 │
346///             └          ┘ └    ┘
347/// ```
348#[derive(Clone, Copy, PartialEq)]
349pub struct Quad(pub [Point; 3]);
350
351impl fmt::Debug for Quad {
352    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
353        let Quad([p0, p1, p2]) = self;
354        write!(f, "Quad {:?} {:?} {:?}", p0, p1, p2)
355    }
356}
357
358impl fmt::Display for Quad {
359    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
360        let Quad([p0, p1, p2]) = self;
361        write!(f, "M{:?} Q{:?} {:?}", p0, p1, p2)
362    }
363}
364
365impl Quad {
366    /// Create new quadratic bezier curve
367    pub fn new(p0: impl Into<Point>, p1: impl Into<Point>, p2: impl Into<Point>) -> Self {
368        Self([p0.into(), p1.into(), p2.into()])
369    }
370
371    /// Points defining quadratic bezier curve
372    pub fn points(&self) -> [Point; 3] {
373        self.0
374    }
375
376    /// Tangent lines at the ends of the quadratic bezier curve
377    pub fn ends(&self) -> (Line, Line) {
378        let Self([p0, p1, p2]) = *self;
379        let start = Line::new(p0, p1);
380        let end = Line::new(p1, p2);
381        if p0.is_close_to(p1) {
382            (end, end)
383        } else if p1.is_close_to(p2) {
384            (start, start)
385        } else {
386            (start, end)
387        }
388    }
389
390    /// Find smooth point used by SVG parser
391    pub fn smooth(&self) -> Point {
392        let Quad([_p0, p1, p2]) = self;
393        2.0 * p2 - *p1
394    }
395}
396
397impl Curve for Quad {
398    /// Flatness criteria for the cubic curve
399    ///
400    /// It is equal to `f = max d(t) where d(t) = |q(t) - l(t)|, l(t) = (1 - t) * p0 + t * p2`
401    /// for q(t) bezier2 curve with p{0..2} control points, in other words maximum distance
402    /// from parametric line to bezier2 curve for the same parameter t.
403    ///
404    /// Line can be represented as bezier2 curve, if `p1 = (p0 + p2) / 2.0`.
405    /// Grouping polynomial coefficients:
406    /// ```text
407    ///     q(t) = t^2 p2 + 2 (1 - t) t p1 + (1 - t)^2 p0
408    ///     l(t) = t^2 p2 + (1 - t) t (p0 + p2) + (1 - t)^2 p0
409    ///     d(t) = |q(t) - l(t)| = (1 - t) t |2 * p1 - p0 - p2|
410    ///     f    = 1 / 4 * | 2 p1 - p0 - p2 |
411    ///     f^2  = 1/16 |2 * p1 - p0 - p2|^2
412    /// ```
413    fn flatness(&self) -> Scalar {
414        let Self([p0, p1, p2]) = *self;
415        let Point([x, y]) = 2.0 * p1 - p0 - p2;
416        x * x + y * y
417    }
418
419    fn transform(&self, tr: Transform) -> Self {
420        let Quad([p0, p1, p2]) = self;
421        Self([tr.apply(*p0), tr.apply(*p1), tr.apply(*p2)])
422    }
423
424    fn start(&self) -> Point {
425        self.0[0]
426    }
427
428    fn end(&self) -> Point {
429        self.0[2]
430    }
431
432    fn at(&self, t: Scalar) -> Point {
433        // at(t) =
434        //   (1 - t) ^ 2 * p0 +
435        //   2 * (1 - t) * t * p1 +
436        //   t ^ 2 * p2
437        let Self([p0, p1, p2]) = self;
438        let (t1, t_1) = (t, 1.0 - t);
439        let (t2, t_2) = (t1 * t1, t_1 * t_1);
440        t_2 * p0 + 2.0 * t1 * t_1 * p1 + t2 * p2
441    }
442
443    fn deriv(&self) -> Segment {
444        let Self([p0, p1, p2]) = *self;
445        Line::new(2.0 * (p1 - p0), 2.0 * (p2 - p1)).into()
446    }
447
448    /// Optimized version of `split_at(0.5)`
449    fn split(&self) -> (Self, Self) {
450        let Self([p0, p1, p2]) = *self;
451        let mid = 0.25 * (p0 + 2.0 * p1 + p2);
452        (
453            Self([p0, 0.5 * (p0 + p1), mid]),
454            Self([mid, 0.5 * (p1 + p2), p2]),
455        )
456    }
457
458    fn split_at(&self, t: Scalar) -> (Self, Self) {
459        // https://pomax.github.io/bezierinfo/#matrixsplit
460        let Self([p0, p1, p2]) = *self;
461        let (t1, t_1) = (t, 1.0 - t);
462        let (t2, t_2) = (t1 * t1, t_1 * t_1);
463        let mid = t_2 * p0 + 2.0 * t1 * t_1 * p1 + t2 * p2;
464        (
465            Self([p0, t_1 * p0 + t * p1, mid]),
466            Self([mid, t_1 * p1 + t * p2, p2]),
467        )
468    }
469
470    fn cut(&self, a: Scalar, b: Scalar) -> Self {
471        // Given quadratic curve as `Q(t) = [1 t t^2] Q P`
472        // where Q is the matrix of the curve, and P is the column vector of the control points
473        // we can change parameter `t in [a, b]` to `s in [0, 1]`, with `t = a + (b - a) * s`
474        //```text
475        //             ┌                         ┐
476        // ┌         ┐ │  1  a       a^2         │
477        // │ 1 s s^2 │ │  0  (b - a) 2*a*(b - a) │ = [1 s s^2 ] T = [1 t t^2]
478        // └         ┘ │  0  0       (b - a)^2   │
479        //             └                         ┘
480        // Q(t in [a, b])
481        //   = [1 s s^2] T Q P
482        //   = [1 s s^2] Q (QI T Q) P
483        //```
484        // => new control points for Q(s) are `P1 = (QI T Q) P`
485        let Self([p0, p1, p2]) = self;
486        let ba = b - a;
487        #[rustfmt::skip]
488        let t = M3x3([
489            1.0, a  , a * a       ,
490            0.0, ba , 2.0 * a * ba,
491            0.0, 0.0, ba * ba     ,
492        ]);
493        #[rustfmt::skip]
494        let M3x3([
495            m00, m01, m02,
496            m10, m11, m12,
497            m20, m21, m22,
498        ]) = QI * t * Q;
499        let q0 = m00 * p0 + m01 * p1 + m02 * p2;
500        let q1 = m10 * p0 + m11 * p1 + m12 * p2;
501        let q2 = m20 * p0 + m21 * p1 + m22 * p2;
502        Self([q0, q1, q2])
503    }
504
505    fn bbox(&self, init: Option<BBox>) -> BBox {
506        let Self([p0, p1, p2]) = self;
507        let bbox = BBox::new(*p0, *p2).union_opt(init);
508        if bbox.contains(*p1) {
509            return bbox;
510        }
511        self.extremities()
512            .fold(bbox, |bbox, t| bbox.extend(self.at(t)))
513    }
514
515    fn offset(&self, dist: Scalar, out: &mut impl Extend<Segment>) {
516        quad_offset_rec(*self, dist, out, 0)
517    }
518
519    fn reverse(&self) -> Self {
520        let Self([p0, p1, p2]) = *self;
521        Self([p2, p1, p0])
522    }
523
524    fn roots(&self) -> CurveRoots {
525        let mut result = CurveRoots::new();
526        // curve(t)_y = 0
527        let Self([Point([_, y0]), Point([_, y1]), Point([_, y2])]) = *self;
528        let a = y0 - 2.0 * y1 + y2;
529        let b = -2.0 * y0 + 2.0 * y1;
530        let c = y0;
531        result.extend(quadratic_solve(a, b, c).filter(|t| (0.0..=1.0).contains(t)));
532        result
533    }
534
535    fn extremities(&self) -> CurveExtremities {
536        let mut result = CurveExtremities::new();
537        let Self([p0, p1, p2]) = self;
538        let Point([a0, a1]) = *p2 - 2.0 * p1 + *p0;
539        let Point([b0, b1]) = *p1 - *p0;
540        // curve'(t)_x = 0
541        if a0.abs() > EPSILON {
542            let t0 = -b0 / a0;
543            if (0.0..=1.0).contains(&t0) {
544                result.push(t0)
545            }
546        }
547        // curve'(t)_y = 0
548        if a1.abs() > EPSILON {
549            let t1 = -b1 / a1;
550            if (0.0..=1.0).contains(&t1) {
551                result.push(t1)
552            }
553        }
554        result
555    }
556
557    fn length(&self, t0: Scalar, t1: Scalar) -> Scalar {
558        // length(t0, t1) = integral(t0, t1, (self'_x ^ 2 + self'_y ^ 2).sqrt())
559        let deriv = self
560            .deriv()
561            .to_line()
562            .expect("derivative of a quad curve must be a line");
563        integrate_quadrature(t0, t1, |t| deriv.at(t).length(), &QUADRATURE_16)
564    }
565}
566
567impl FromStr for Quad {
568    type Err = SvgParserError;
569
570    fn from_str(text: &str) -> Result<Self, Self::Err> {
571        let segment = Segment::from_str(text)?;
572        segment
573            .to_quad()
574            .ok_or(SvgParserError::UnexpectedSegmentType)
575    }
576}
577
578// -----------------------------------------------------------------------------
579// Cubic bezier curve
580// -----------------------------------------------------------------------------
581
582/// Matrix form for cubic bezier curve
583#[rustfmt::skip]
584const C: M4x4 = M4x4([
585    1.0,  0.0,  0.0, 0.0,
586   -3.0,  3.0,  0.0, 0.0,
587    3.0, -6.0,  3.0, 0.0,
588   -1.0,  3.0, -3.0, 1.0,
589]);
590
591/// Inverted matrix form for cubic bezier curve
592#[rustfmt::skip]
593const CI: M4x4 = M4x4([
594    1.0, 0.0      , 0.0      , 0.0,
595    1.0, 1.0 / 3.0, 0.0      , 0.0,
596    1.0, 2.0 / 3.0, 1.0 / 3.0, 0.0,
597    1.0, 1.0      , 1.0      , 1.0,
598]);
599
600/// Cubic bezier curve
601///
602/// Polynomial form:
603/// `(1 - t) ^ 3 * p0 + 3 * (1 - t) ^ 2 * t * p1 + 3 * (1 - t) * t ^ 2 * p2 + t ^ 3 * p3`
604/// Matrix from:
605/// ```text
606///                 ┌             ┐ ┌    ┐
607/// ┌             ┐ │  1  0  0  0 │ │ p0 │
608/// │ 1 t t^2 t^3 │ │ -3  3  0  0 │ │ p1 │
609/// └             ┘ │  3 -6  3  0 │ │ p2 │
610///                 │ -1  3 -3  1 │ │ p3 │
611///                 └             ┘ └    ┘
612/// ```
613#[derive(Clone, Copy, PartialEq)]
614pub struct Cubic(pub [Point; 4]);
615
616impl fmt::Debug for Cubic {
617    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
618        let Cubic([p0, p1, p2, p3]) = self;
619        write!(f, "Cubic {:?} {:?} {:?} {:?}", p0, p1, p2, p3)
620    }
621}
622
623impl fmt::Display for Cubic {
624    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
625        let Cubic([p0, p1, p2, p3]) = self;
626        write!(f, "M{:?} C{:?} {:?} {:?}", p0, p1, p2, p3)
627    }
628}
629
630impl Cubic {
631    /// Create new cubic bezier curve
632    pub fn new(
633        p0: impl Into<Point>,
634        p1: impl Into<Point>,
635        p2: impl Into<Point>,
636        p3: impl Into<Point>,
637    ) -> Self {
638        Self([p0.into(), p1.into(), p2.into(), p3.into()])
639    }
640
641    /// Points defining cubic bezier curve
642    pub fn points(&self) -> [Point; 4] {
643        self.0
644    }
645
646    /// Tangent lines at the ends of cubic bezier curve
647    pub fn ends(&self) -> (Line, Line) {
648        let ps = self.points();
649        let mut start = 0;
650        for i in 0..3 {
651            if !ps[i].is_close_to(ps[i + 1]) {
652                start = i;
653                break;
654            }
655        }
656        let mut end = 0;
657        for i in (1..4).rev() {
658            if !ps[i].is_close_to(ps[i - 1]) {
659                end = i;
660                break;
661            }
662        }
663        (
664            Line::new(ps[start], ps[start + 1]),
665            Line::new(ps[end - 1], ps[end]),
666        )
667    }
668
669    /// Find smooth point used by SVG parser
670    pub fn smooth(&self) -> Point {
671        let Cubic([_p0, _p1, p2, p3]) = self;
672        2.0 * p3 - *p2
673    }
674}
675
676impl Curve for Cubic {
677    /// Flatness criteria for the cubic curve
678    /// This function actually returns `16 * flatness^2`
679    ///
680    /// It is equal to `f = max d(t) where d(t) = |c(t) - l(t)|, l(t) = (1 - t) * c0 + t * c3`
681    /// for c(t) bezier3 curve with c{0..3} control points, in other words maximum distance
682    /// from parametric line to bezier3 curve for the same parameter t. It is shown in the article
683    /// that:
684    ///     f^2 <= 1/16 (max{u_x^2, v_x^2} + max{u_y^2, v_y^2})
685    /// where:
686    ///     u = 3 * b1 - 2 * b0 - b3
687    ///     v = 3 * b2 - b0 - 2 * b3
688    /// `f == 0` means completely flat so estimating upper bound is sufficient as splitting more
689    /// than needed is not a problem for rendering.
690    ///
691    /// [Linear Approximation of Bezier Curve](https://hcklbrrfnn.files.wordpress.com/2012/08/bez.pdf)
692    fn flatness(&self) -> Scalar {
693        let Self([p0, p1, p2, p3]) = *self;
694        let u = 3.0 * p1 - 2.0 * p0 - p3;
695        let v = 3.0 * p2 - p0 - 2.0 * p3;
696        (u.x() * u.x()).max(v.x() * v.x()) + (u.y() * u.y()).max(v.y() * v.y())
697    }
698
699    fn transform(&self, tr: Transform) -> Self {
700        let Cubic([p0, p1, p2, p3]) = self;
701        Self([tr.apply(*p0), tr.apply(*p1), tr.apply(*p2), tr.apply(*p3)])
702    }
703
704    fn start(&self) -> Point {
705        self.0[0]
706    }
707
708    fn end(&self) -> Point {
709        self.0[3]
710    }
711
712    fn at(&self, t: Scalar) -> Point {
713        // at(t) =
714        //   (1 - t) ^ 3 * p0 +
715        //   3 * (1 - t) ^ 2 * t * p1 +
716        //   3 * (1 - t) * t ^ 2 * p2 +
717        //   t ^ 3 * p3
718        let Self([p0, p1, p2, p3]) = self;
719        let (t1, t_1) = (t, 1.0 - t);
720        let (t2, t_2) = (t1 * t1, t_1 * t_1);
721        let (t3, t_3) = (t2 * t1, t_2 * t_1);
722        t_3 * p0 + 3.0 * t1 * t_2 * p1 + 3.0 * t2 * t_1 * p2 + t3 * p3
723    }
724
725    fn deriv(&self) -> Segment {
726        let Self([p0, p1, p2, p3]) = *self;
727        Quad::new(3.0 * (p1 - p0), 3.0 * (p2 - p1), 3.0 * (p3 - p2)).into()
728    }
729
730    /// Optimized version of `split_at(0.5)`
731    fn split(&self) -> (Self, Self) {
732        let Self([p0, p1, p2, p3]) = *self;
733        let mid = 0.125 * p0 + 0.375 * p1 + 0.375 * p2 + 0.125 * p3;
734        let c0 = Self([
735            p0,
736            0.5 * p0 + 0.5 * p1,
737            0.25 * p0 + 0.5 * p1 + 0.25 * p2,
738            mid,
739        ]);
740        let c1 = Self([
741            mid,
742            0.25 * p1 + 0.5 * p2 + 0.25 * p3,
743            0.5 * p2 + 0.5 * p3,
744            p3,
745        ]);
746        (c0, c1)
747    }
748
749    fn split_at(&self, t: Scalar) -> (Self, Self) {
750        // https://pomax.github.io/bezierinfo/#matrixsplit
751        let Self([p0, p1, p2, p3]) = self;
752        let (t1, t_1) = (t, 1.0 - t);
753        let (t2, t_2) = (t1 * t1, t_1 * t_1);
754        let (t3, t_3) = (t2 * t1, t_2 * t_1);
755        let mid = t_3 * p0 + 3.0 * t1 * t_2 * p1 + 3.0 * t2 * t_1 * p2 + t3 * p3;
756        let c0 = Self([
757            *p0,
758            t_1 * p0 + t * p1,
759            t_2 * p0 + 2.0 * t * t_1 * p1 + t2 * p2,
760            mid,
761        ]);
762        let c1 = Self([
763            mid,
764            t_2 * p1 + 2.0 * t * t_1 * p2 + t2 * p3,
765            t_1 * p2 + t * p3,
766            *p3,
767        ]);
768        (c0, c1)
769    }
770
771    fn cut(&self, a: Scalar, b: Scalar) -> Self {
772        // Given cubic curve as `C(t) = [1 t t^2 t^3] C P`
773        // where C is the matrix of the curve, and P is the column vector of the control points
774        // we can change parameter `t in [a, b]` to `s in [0, 1]`, with `t = a + (b - a) * s`
775        // ```text
776        //                                   ┌                                       ┐
777        // ┌             ┐   ┌             ┐ │  1  a       a^2         a^3           │
778        // │ 1 t t^2 t^3 │ = │ 1 s s^2 s^3 │ │  0  (b - a) 2*a*(b - a) 3*a^2*(b - a) │ = [1 s s^2 s^3] T
779        // └             ┘   └             ┘ │  0  0       (b - a)^2   3*a*(b - a)^2 │
780        //                                   │  0  0       0           (b - a)^3     │
781        //                                   └                                       ┘
782        // C(t in [a, b])
783        //   = [1 s s^2 s^3] T C P
784        //   = [1 s s^2 s^3] C (CI T C) P
785        // ```
786        // => new control points for C(s) are P1 = (CI T C) P
787        let Self([p0, p1, p2, p3]) = self;
788        let ba = b - a;
789        #[rustfmt::skip]
790        let t = M4x4([
791            1.0, a  , a * a       , a * a * a        ,
792            0.0, ba , 2.0 * a * ba, 3.0 * a * a * ba ,
793            0.0, 0.0, ba * ba     , 3.0 * a * ba * ba,
794            0.0, 0.0, 0.0         , ba * ba * ba     ,
795        ]);
796        #[rustfmt::skip]
797        let M4x4([
798            m00, m01, m02, m03,
799            m10, m11, m12, m13,
800            m20, m21, m22, m23,
801            m30, m31, m32, m33,
802        ]) = CI * t * C;
803        let c0 = m00 * p0 + m01 * p1 + m02 * p2 + m03 * p3;
804        let c1 = m10 * p0 + m11 * p1 + m12 * p2 + m13 * p3;
805        let c2 = m20 * p0 + m21 * p1 + m22 * p2 + m23 * p3;
806        let c3 = m30 * p0 + m31 * p1 + m32 * p2 + m33 * p3;
807        Self([c0, c1, c2, c3])
808    }
809
810    fn bbox(&self, init: Option<BBox>) -> BBox {
811        let Self([p0, p1, p2, p3]) = self;
812        let bbox = BBox::new(*p0, *p3).union_opt(init);
813        if bbox.contains(*p1) && bbox.contains(*p2) {
814            return bbox;
815        }
816        self.extremities()
817            .fold(bbox, |bbox, t| bbox.extend(self.at(t)))
818    }
819
820    /// Offset cubic bezier curve with a list of cubic curves
821    ///
822    /// Offset bezier curve using Tiller-Hanson method. In short, it will just offset
823    /// line segment corresponding to control points, then find intersection of this
824    /// lines and treat them as new control points.
825    fn offset(&self, dist: Scalar, out: &mut impl Extend<Segment>) {
826        cubic_offset_rec(*self, None, dist, out, 0);
827    }
828
829    fn reverse(&self) -> Self {
830        let Self([p0, p1, p2, p3]) = *self;
831        Self([p3, p2, p1, p0])
832    }
833
834    fn roots(&self) -> CurveRoots {
835        let mut result = CurveRoots::new();
836        // curve(t)_y = 0
837        let Self(
838            [
839                Point([_, y0]),
840                Point([_, y1]),
841                Point([_, y2]),
842                Point([_, y3]),
843            ],
844        ) = *self;
845        let a = -y0 + 3.0 * y1 - 3.0 * y2 + y3;
846        let b = 3.0 * y0 - 6.0 * y1 + 3.0 * y2;
847        let c = -3.0 * y0 + 3.0 * y1;
848        let d = y0;
849        result.extend(cubic_solve(a, b, c, d).filter(|t| (0.0..=1.0).contains(t)));
850        result
851    }
852
853    fn extremities(&self) -> CurveExtremities {
854        let Self([p0, p1, p2, p3]) = *self;
855        let Point([a0, a1]) = -1.0 * p0 + 3.0 * p1 - 3.0 * p2 + 1.0 * p3;
856        let Point([b0, b1]) = 2.0 * p0 - 4.0 * p1 + 2.0 * p2;
857        let Point([c0, c1]) = -1.0 * p0 + p1;
858
859        // Solve for `curve'(t)_x = 0 || curve'(t)_y = 0`
860        quadratic_solve(a0, b0, c0)
861            .chain(quadratic_solve(a1, b1, c1))
862            .filter(|t| *t >= 0.0 && *t <= 1.0)
863            .collect::<CurveExtremities>()
864    }
865
866    fn length(&self, t0: Scalar, t1: Scalar) -> Scalar {
867        // length(t0, t1) = integral(t0, t1, (self'_x ^ 2 + self'_y ^ 2).sqrt())
868        let deriv = self
869            .deriv()
870            .to_quad()
871            .expect("derivative of a cubic curve must be a quad curve");
872        integrate_quadrature(t0, t1, |t| deriv.at(t).length(), &QUADRATURE_32)
873    }
874}
875
876impl From<Quad> for Cubic {
877    fn from(quad: Quad) -> Self {
878        let Quad([p0, p1, p2]) = quad;
879        Self([
880            p0,
881            (1.0 / 3.0) * p0 + (2.0 / 3.0) * p1,
882            (2.0 / 3.0) * p1 + (1.0 / 3.0) * p2,
883            p2,
884        ])
885    }
886}
887
888impl FromStr for Cubic {
889    type Err = SvgParserError;
890
891    fn from_str(text: &str) -> Result<Self, Self::Err> {
892        let segment = Segment::from_str(text)?;
893        segment
894            .to_cubic()
895            .ok_or(SvgParserError::UnexpectedSegmentType)
896    }
897}
898
899// -----------------------------------------------------------------------------
900// Segment
901// -----------------------------------------------------------------------------
902
903/// `Segment` is an enum of either [`Line`], [`Quad`] or [`Cubic`]
904#[derive(Clone, Copy, PartialEq)]
905pub enum Segment {
906    Line(Line),
907    Quad(Quad),
908    Cubic(Cubic),
909}
910
911impl Segment {
912    pub fn ends(&self) -> (Line, Line) {
913        match self {
914            Segment::Line(line) => line.ends(),
915            Segment::Quad(quad) => quad.ends(),
916            Segment::Cubic(cubic) => cubic.ends(),
917        }
918    }
919
920    /// Find intersection between two segments
921    ///
922    /// This might not be the fastest method possible but works for any two curves.
923    /// Divide curves as long as there is intersection between bounding boxes, if
924    /// the intersection is smaller than tolerance we can treat it as an intersection point.
925    ///
926    /// NOTE: This produces duplicate results which are close to each other. Why?
927    pub fn intersect(self, other: impl Into<Segment>, tolerance: Scalar) -> Vec<Point> {
928        let mut queue = vec![(self, other.into())];
929        let mut result: Vec<Point> = Vec::new();
930        while let Some((s0, s1)) = queue.pop() {
931            let b0 = s0.bbox(None);
932            let b1 = s1.bbox(None);
933            if let Some(b) = b0.intersect(b1) {
934                if b.diag().length() < tolerance {
935                    let p_new = b.diag().at(0.5);
936                    if result.iter().all(|p| p.dist(p_new) > tolerance) {
937                        result.push(p_new);
938                    }
939                } else {
940                    // TODO: can be optimized by splitting only curves with large bounding box
941                    let (s00, s01) = s0.split();
942                    let (s10, s11) = s1.split();
943                    queue.push((s00, s10));
944                    queue.push((s00, s11));
945                    queue.push((s01, s10));
946                    queue.push((s01, s11));
947                }
948            }
949        }
950        result
951    }
952
953    /// Convert to line if it is a line variant of the segment
954    pub fn to_line(&self) -> Option<Line> {
955        match self {
956            Segment::Line(line) => Some(*line),
957            _ => None,
958        }
959    }
960
961    /// Convert to quad if it is a quad variant of the segment
962    pub fn to_quad(&self) -> Option<Quad> {
963        match self {
964            Segment::Quad(quad) => Some(*quad),
965            _ => None,
966        }
967    }
968
969    /// Convert to cubic if it is a cubic variant of the segment
970    pub fn to_cubic(&self) -> Option<Cubic> {
971        match self {
972            Segment::Cubic(cubic) => Some(*cubic),
973            _ => None,
974        }
975    }
976
977    /// Produce iterator over segments that join to segments with the specified method.
978    pub fn line_join(
979        self,
980        other: Segment,
981        stroke_style: StrokeStyle,
982    ) -> impl Iterator<Item = Self> {
983        let mut result = ArrayIter::<Segment, 4>::new();
984        if self.end().is_close_to(other.start()) {
985            return result;
986        }
987        let bevel = Line::new(self.end(), other.start());
988        // https://www.w3.org/TR/SVG2/painting.html#LineJoin
989        match stroke_style.line_join {
990            LineJoin::Bevel => {
991                result.push(bevel.into());
992            }
993            LineJoin::Miter(miter_limit) => {
994                let (_, start) = self.ends();
995                let (end, _) = other.ends();
996                match start.intersect(end) {
997                    Some((t0, t1)) if (0.0..=1.0).contains(&t0) && (0.0..=1.0).contains(&t1) => {
998                        // ends intersect
999                        result.push(bevel.into());
1000                    }
1001                    None => result.push(bevel.into()),
1002                    Some((t, _)) => {
1003                        let p0 = start.end() - start.start();
1004                        let p1 = end.start() - end.end();
1005                        // miter_length = stroke_width / sin(a / 2)
1006                        // sin(a / 2) = +/- ((1 - cos(a)) / 2).sqrt()
1007                        let miter_length = p0
1008                            .cos_between(p1)
1009                            .map(|c| stroke_style.width / ((1.0 - c) / 2.0).sqrt());
1010                        match miter_length {
1011                            Some(miter_length) if miter_length < miter_limit => {
1012                                let p = start.at(t);
1013                                result.push(Line::new(start.end(), p).into());
1014                                result.push(Line::new(p, end.start()).into());
1015                            }
1016                            _ => result.push(bevel.into()),
1017                        }
1018                    }
1019                }
1020            }
1021            LineJoin::Round => {
1022                let (_, start) = self.ends();
1023                let (end, _) = other.ends();
1024                match start.intersect_point(end) {
1025                    Some(_) => result.push(bevel.into()),
1026                    None => {
1027                        let sweep_flag = start.direction().cross(bevel.direction()) >= 0.0;
1028                        let radius = stroke_style.width / 2.0;
1029                        let arc = EllipArc::new_param(
1030                            start.end(),
1031                            end.start(),
1032                            radius,
1033                            radius,
1034                            0.0,
1035                            false,
1036                            sweep_flag,
1037                        );
1038                        match arc {
1039                            Some(arc) => result.extend(arc.to_cubics().map(Segment::from)),
1040                            None => result.push(bevel.into()),
1041                        }
1042                    }
1043                }
1044            }
1045        }
1046        result
1047    }
1048
1049    /// Produce and iterator over segments that adds caps between two segments
1050    pub fn line_cap(self, other: Segment, stroke_style: StrokeStyle) -> impl Iterator<Item = Self> {
1051        let mut result = ArrayIter::<Segment, 4>::new();
1052        if self.end().is_close_to(other.start()) {
1053            return result;
1054        }
1055        let butt = Line::new(self.end(), other.start());
1056        match stroke_style.line_cap {
1057            LineCap::Butt => result.push(butt.into()),
1058            LineCap::Square => {
1059                let (_, from) = self.ends();
1060                if let Some(tang) = from.direction().normalize() {
1061                    let l0 = Line::new(self.end(), self.end() + stroke_style.width / 2.0 * tang);
1062                    result.push(l0.into());
1063                    let l1 = Line::new(l0.end(), l0.end() + butt.direction());
1064                    result.push(l1.into());
1065                    let l2 = Line::new(l1.end(), other.start());
1066                    result.push(l2.into());
1067                }
1068            }
1069            LineCap::Round => {
1070                let stroke_style = StrokeStyle {
1071                    line_join: LineJoin::Round,
1072                    ..stroke_style
1073                };
1074                result.extend(self.line_join(other, stroke_style));
1075            }
1076        }
1077        result
1078    }
1079
1080    /// Whether curve has any NaN values
1081    pub fn has_nans(&self) -> bool {
1082        match self {
1083            Segment::Line(line) => line.0.as_slice(),
1084            Segment::Quad(quad) => quad.0.as_slice(),
1085            Segment::Cubic(cubic) => cubic.0.as_slice(),
1086        }
1087        .iter()
1088        .any(|point| point.0.iter().any(|f| f.is_nan()))
1089    }
1090}
1091
1092impl fmt::Debug for Segment {
1093    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1094        match self {
1095            Segment::Line(line) => line.fmt(f),
1096            Segment::Quad(quad) => quad.fmt(f),
1097            Segment::Cubic(cubic) => cubic.fmt(f),
1098        }
1099    }
1100}
1101
1102impl fmt::Display for Segment {
1103    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1104        match self {
1105            Segment::Line(line) => line.fmt(f),
1106            Segment::Quad(quad) => quad.fmt(f),
1107            Segment::Cubic(cubic) => cubic.fmt(f),
1108        }
1109    }
1110}
1111
1112impl FromStr for Segment {
1113    type Err = SvgParserError;
1114
1115    fn from_str(text: &str) -> Result<Self, Self::Err> {
1116        use SvgPathCmd::*;
1117        let mut iter = SvgPathParser::new(Cursor::new(text));
1118        match [
1119            iter.next().transpose()?,
1120            iter.next().transpose()?,
1121            iter.next().transpose()?,
1122        ] {
1123            [Some(MoveTo(p0)), Some(curve), None] => {
1124                let segment: Segment = match curve {
1125                    LineTo(p1) => Line::new(p0, p1).into(),
1126                    QuadTo(p1, p2) => Quad::new(p0, p1, p2).into(),
1127                    CubicTo(p1, p2, p3) => Cubic::new(p0, p1, p2, p3).into(),
1128                    _ => return Err(SvgParserError::UnexpectedSegmentType),
1129                };
1130                Ok(segment)
1131            }
1132            _ => Err(SvgParserError::UnexpectedSegmentType),
1133        }
1134    }
1135}
1136
1137impl Curve for Segment {
1138    fn flatness(&self) -> Scalar {
1139        match self {
1140            Segment::Line(line) => line.flatness(),
1141            Segment::Quad(quad) => quad.flatness(),
1142            Segment::Cubic(cubic) => cubic.flatness(),
1143        }
1144    }
1145
1146    fn transform(&self, tr: Transform) -> Self {
1147        match self {
1148            Segment::Line(line) => line.transform(tr).into(),
1149            Segment::Quad(quad) => quad.transform(tr).into(),
1150            Segment::Cubic(cubic) => cubic.transform(tr).into(),
1151        }
1152    }
1153
1154    fn start(&self) -> Point {
1155        match self {
1156            Segment::Line(line) => line.start(),
1157            Segment::Quad(quad) => quad.start(),
1158            Segment::Cubic(cubic) => cubic.start(),
1159        }
1160    }
1161
1162    fn end(&self) -> Point {
1163        match self {
1164            Segment::Line(line) => line.end(),
1165            Segment::Quad(quad) => quad.end(),
1166            Segment::Cubic(cubic) => cubic.end(),
1167        }
1168    }
1169
1170    fn at(&self, t: Scalar) -> Point {
1171        match self {
1172            Segment::Line(line) => line.at(t),
1173            Segment::Quad(quad) => quad.at(t),
1174            Segment::Cubic(cubic) => cubic.at(t),
1175        }
1176    }
1177
1178    fn deriv(&self) -> Segment {
1179        match self {
1180            Segment::Line(line) => line.deriv(),
1181            Segment::Quad(quad) => quad.deriv(),
1182            Segment::Cubic(cubic) => cubic.deriv(),
1183        }
1184    }
1185
1186    fn split_at(&self, t: Scalar) -> (Self, Self) {
1187        match self {
1188            Segment::Line(line) => {
1189                let (l0, l1) = line.split_at(t);
1190                (l0.into(), l1.into())
1191            }
1192            Segment::Quad(quad) => {
1193                let (q0, q1) = quad.split_at(t);
1194                (q0.into(), q1.into())
1195            }
1196            Segment::Cubic(cubic) => {
1197                let (c0, c1) = cubic.split_at(t);
1198                (c0.into(), c1.into())
1199            }
1200        }
1201    }
1202
1203    fn cut(&self, a: Scalar, b: Scalar) -> Self {
1204        match self {
1205            Segment::Line(line) => line.cut(a, b).into(),
1206            Segment::Quad(quad) => quad.cut(a, b).into(),
1207            Segment::Cubic(cubic) => cubic.cut(a, b).into(),
1208        }
1209    }
1210
1211    fn bbox(&self, init: Option<BBox>) -> BBox {
1212        match self {
1213            Segment::Line(line) => line.bbox(init),
1214            Segment::Quad(quad) => quad.bbox(init),
1215            Segment::Cubic(cubic) => cubic.bbox(init),
1216        }
1217    }
1218
1219    fn offset(&self, dist: Scalar, out: &mut impl Extend<Segment>) {
1220        match self {
1221            Segment::Line(line) => line.offset(dist, out),
1222            Segment::Quad(quad) => quad.offset(dist, out),
1223            Segment::Cubic(cubic) => cubic.offset(dist, out),
1224        }
1225    }
1226
1227    fn reverse(&self) -> Self {
1228        match self {
1229            Segment::Line(line) => line.reverse().into(),
1230            Segment::Quad(quad) => quad.reverse().into(),
1231            Segment::Cubic(cubic) => cubic.reverse().into(),
1232        }
1233    }
1234
1235    fn roots(&self) -> CurveRoots {
1236        match self {
1237            Segment::Line(line) => line.roots(),
1238            Segment::Quad(quad) => quad.roots(),
1239            Segment::Cubic(cubic) => cubic.roots(),
1240        }
1241    }
1242
1243    fn extremities(&self) -> CurveExtremities {
1244        match self {
1245            Segment::Line(line) => line.extremities(),
1246            Segment::Quad(quad) => quad.extremities(),
1247            Segment::Cubic(cubic) => cubic.extremities(),
1248        }
1249    }
1250
1251    fn length(&self, t0: Scalar, t1: Scalar) -> Scalar {
1252        match self {
1253            Segment::Line(line) => Curve::length(line, t0, t1),
1254            Segment::Quad(quad) => quad.length(t0, t1),
1255            Segment::Cubic(cubic) => cubic.length(t0, t1),
1256        }
1257    }
1258}
1259
1260impl From<Line> for Segment {
1261    fn from(line: Line) -> Self {
1262        Self::Line(line)
1263    }
1264}
1265
1266impl From<Quad> for Segment {
1267    fn from(quad: Quad) -> Self {
1268        Self::Quad(quad)
1269    }
1270}
1271
1272impl From<Cubic> for Segment {
1273    fn from(cubic: Cubic) -> Self {
1274        Self::Cubic(cubic)
1275    }
1276}
1277
1278// -----------------------------------------------------------------------------
1279// Bezier curve offsetting helpers
1280// -----------------------------------------------------------------------------
1281
1282/// Offset line to the distance.
1283pub(crate) fn line_offset(line: Line, dist: Scalar) -> Option<Line> {
1284    let Line([p0, p1]) = line;
1285    let offset = dist * (p1 - p0).normal().normalize()?;
1286    Some(Line::new(p0 + offset, p1 + offset))
1287}
1288
1289/// Offset poly-line specified by points `ps`.
1290///
1291/// Implementation correctly handles repeated points.
1292/// False result indicates that all points are close to each other
1293/// and it is not possible to offset them.
1294fn polyline_offset(ps: &mut [Point], dist: Scalar) -> bool {
1295    if ps.is_empty() {
1296        return true;
1297    }
1298    let mut prev: Option<Line> = None;
1299    let mut index = 0;
1300    loop {
1301        // find identical points repeats
1302        let mut repeats = 1;
1303        for i in index..ps.len() - 1 {
1304            if !ps[i].is_close_to(ps[i + 1]) {
1305                break;
1306            }
1307            repeats += 1;
1308        }
1309        if index + repeats >= ps.len() {
1310            break;
1311        }
1312        index += repeats;
1313        // offset line
1314        let next = line_offset(Line::new(ps[index - 1], ps[index]), dist)
1315            .expect("polyline implementation error");
1316        // find where to move repeated points
1317        let point = match prev {
1318            None => next.start(),
1319            Some(prev) => match prev.intersect(next) {
1320                Some((t, _)) => prev.at(t),
1321                None => {
1322                    // TODO: not sure what to do especially for up/down move
1323                    next.start()
1324                }
1325            },
1326        };
1327        // move repeats
1328        for p in ps.iter_mut().take(index).skip(index - repeats) {
1329            *p = point;
1330        }
1331        prev = Some(next);
1332    }
1333    // handle tail points
1334    match prev {
1335        None => {
1336            // all points are close to each other, can not offset.
1337            false
1338        }
1339        Some(prev) => {
1340            for p in ps.iter_mut().skip(index) {
1341                *p = prev.end();
1342            }
1343            true
1344        }
1345    }
1346}
1347
1348/// Determine if quad curve needs splitting before offsetting.
1349fn quad_offset_should_split(quad: Quad) -> bool {
1350    let Quad([p0, p1, p2]) = quad;
1351    // split if angle is sharp
1352    if (p0 - p1).dot(p2 - p1) > 0.0 {
1353        return true;
1354    }
1355    // distance between center mass and midpoint of a curve,
1356    // should be bigger then 0.1 of the bounding box diagonal
1357    let c_mass = (p0 + p1 + p2) / 3.0;
1358    let c_mid = quad.at(0.5);
1359    let dist = (c_mass - c_mid).length();
1360    let bbox_diag = quad.bbox(None).diag().length();
1361    bbox_diag * 0.1 < dist
1362}
1363
1364/// Recursive quad offset.
1365fn quad_offset_rec(quad: Quad, dist: Scalar, out: &mut impl Extend<Segment>, depth: usize) {
1366    if quad_offset_should_split(quad) && depth < 3 {
1367        let (c0, c1) = quad.split_at(0.5);
1368        quad_offset_rec(c0, dist, out, depth + 1);
1369        quad_offset_rec(c1, dist, out, depth + 1);
1370    } else {
1371        let mut points = quad.points();
1372        if polyline_offset(&mut points, dist) {
1373            out.extend(Some(Quad(points).into()));
1374        }
1375    }
1376}
1377
1378/// Recursive cubic offset.
1379fn cubic_offset_rec(
1380    cubic: Cubic,
1381    last: Option<Segment>,
1382    dist: Scalar,
1383    out: &mut impl Extend<Segment>,
1384    depth: usize,
1385) -> Option<Segment> {
1386    if cubic_offset_should_split(cubic) && depth < 3 {
1387        let (c0, c1) = cubic.split();
1388        let last = cubic_offset_rec(c0, last, dist, out, depth + 1);
1389        return cubic_offset_rec(c1, last, dist, out, depth + 1);
1390    }
1391    let mut points = cubic.points();
1392    if polyline_offset(&mut points, dist) {
1393        let result: Segment = Cubic(points).into();
1394        // there could a disconnect between offset curves.
1395        // For example M0,0 C10,5 0,5 10,0.
1396        if let Some(last) = last
1397            && !last.end().is_close_to(result.start())
1398        {
1399            let stroke_style = StrokeStyle {
1400                width: dist * 2.0,
1401                line_join: LineJoin::Round,
1402                line_cap: LineCap::Round,
1403            };
1404            out.extend(last.line_join(result, stroke_style));
1405        }
1406        out.extend(Some(result));
1407        Some(result)
1408    } else {
1409        last
1410    }
1411}
1412
1413/// Determine if cubic curve needs splitting before offsetting.
1414fn cubic_offset_should_split(cubic: Cubic) -> bool {
1415    let Cubic([p0, p1, p2, p3]) = cubic;
1416    // angle(p3 - p0, p2 - p1) > 90 or < -90
1417    if (p3 - p0).dot(p2 - p1) < 0.0 {
1418        return true;
1419    }
1420    // control points should be on the same side of the baseline.
1421    let a0 = (p3 - p0).cross(p1 - p0);
1422    let a1 = (p3 - p0).cross(p2 - p0);
1423    if a0 * a1 < 0.0 {
1424        return true;
1425    }
1426    // distance between center mass and midpoint of a curve,
1427    // should be bigger than 0.1 of the bounding box diagonal
1428    let c_mass = (p0 + p1 + p2 + p3) / 4.0;
1429    let c_mid = cubic.at(0.5);
1430    let dist = (c_mass - c_mid).length();
1431    let bbox_diag = cubic.bbox(None).diag().length();
1432    bbox_diag * 0.1 < dist
1433}
1434
1435#[cfg(test)]
1436mod tests {
1437    use super::*;
1438    use crate::{assert_approx_eq, assert_approx_eq_iter};
1439
1440    #[test]
1441    fn test_polyline_offset() {
1442        let dist = -(2.0 as Scalar).sqrt();
1443
1444        let p0 = Point::new(1.0, 0.0);
1445        let r0 = Point::new(0.0, 1.0);
1446
1447        let p1 = Point::new(2.0, 1.0);
1448        let r1 = Point::new(2.0, 3.0);
1449
1450        let p2 = Point::new(3.0, 0.0);
1451        let r2 = Point::new(4.0, 1.0);
1452
1453        // basic
1454        let mut ps = [p0, p1, p2];
1455        assert!(polyline_offset(&mut ps, dist));
1456        assert_eq!(&ps, &[r0, r1, r2]);
1457
1458        // repeat at start
1459        let mut ps = [p0, p0, p1, p2];
1460        assert!(polyline_offset(&mut ps, dist));
1461        assert_eq!(&ps, &[r0, r0, r1, r2]);
1462
1463        // repeat at start
1464        let mut ps = [p0, p1, p2, p2];
1465        assert!(polyline_offset(&mut ps, dist));
1466        assert_eq!(&ps, &[r0, r1, r2, r2]);
1467
1468        // repeat in the middle
1469        let mut ps = [p0, p1, p1, p2];
1470        assert!(polyline_offset(&mut ps, dist));
1471        assert_eq!(&ps, &[r0, r1, r1, r2]);
1472
1473        // all points are close to each other, can not offset
1474        let mut ps = [p0, p0, p0];
1475        assert!(!polyline_offset(&mut ps, dist));
1476
1477        // splitted single line
1478        let mut ps = [p0, p1, Point::new(3.0, 2.0)];
1479        assert!(polyline_offset(&mut ps, dist));
1480        assert_eq!(&ps, &[r0, Point::new(1.0, 2.0), r1]);
1481
1482        // four points
1483        let mut ps = [p0, p1, p2, Point::new(2.0, -1.0)];
1484        assert!(polyline_offset(&mut ps, dist));
1485        assert_eq!(&ps, &[r0, r1, Point::new(5.0, 0.0), Point::new(3.0, -2.0)]);
1486    }
1487
1488    #[test]
1489    fn test_roots() {
1490        let l: Line = "M0-1 2 1".parse().unwrap();
1491        assert_eq!(l.roots().collect::<Vec<_>>(), vec![0.5]);
1492
1493        let q: Quad = "M0-2Q7,6 6-4".parse().unwrap();
1494        assert_approx_eq_iter!(q.roots(), [0.73841681234051, 0.15047207654837882]);
1495
1496        let c = Cubic::new((0.0, -2.0), (2.0, 4.0), (4.0, -3.0), (9.0, 1.0));
1497        assert_approx_eq_iter!(
1498            c.roots(),
1499            [0.8812186869024836, 0.1627575589800928, 0.5810237541174236]
1500        );
1501
1502        let c: Cubic = "M8,-1 C1,3 6,-3 9,1".parse().unwrap();
1503        assert_approx_eq_iter!(
1504            c.roots().collect::<Vec<_>>(),
1505            [0.8872983346207419, 0.11270166537925835, 0.4999999999999995]
1506        );
1507    }
1508
1509    #[test]
1510    fn test_curve_matrices() {
1511        #[rustfmt::skip]
1512        let i3 = [
1513            1.0, 0.0, 0.0,
1514            0.0, 1.0, 0.0,
1515            0.0, 0.0, 1.0
1516        ];
1517        let M3x3(q) = Q * QI;
1518        assert_approx_eq_iter!(i3, q);
1519
1520        #[rustfmt::skip]
1521        let i4 = [
1522            1.0, 0.0, 0.0, 0.0,
1523            0.0, 1.0, 0.0, 0.0,
1524            0.0, 0.0, 1.0, 0.0,
1525            0.0, 0.0, 0.0, 1.0,
1526        ];
1527        let M4x4(c) = C * CI;
1528        assert_approx_eq_iter!(i4, c);
1529    }
1530
1531    #[test]
1532    fn test_ends() {
1533        let p0 = Point::new(1.0, 0.0);
1534        let p1 = Point::new(2.0, 1.0);
1535        let p2 = Point::new(3.0, 0.0);
1536        let p3 = Point::new(2.0, 0.0);
1537
1538        let c = Cubic::new(p0, p1, p2, p3);
1539        let (start, end) = c.ends();
1540        assert_eq!(start, Line::new(p0, p1));
1541        assert_eq!(end, Line::new(p2, p3));
1542
1543        let c = Cubic::new(p0, p0, p1, p2);
1544        let (start, end) = c.ends();
1545        assert_eq!(start, Line::new(p0, p1));
1546        assert_eq!(end, Line::new(p1, p2));
1547
1548        let c = Cubic::new(p0, p1, p2, p2);
1549        let (start, end) = c.ends();
1550        assert_eq!(start, Line::new(p0, p1));
1551        assert_eq!(end, Line::new(p1, p2));
1552
1553        let q = Quad::new(p0, p1, p2);
1554        let (start, end) = q.ends();
1555        assert_eq!(start, Line::new(p0, p1));
1556        assert_eq!(end, Line::new(p1, p2));
1557
1558        let q = Quad::new(p0, p0, p1);
1559        let (start, end) = q.ends();
1560        assert_eq!(start, Line::new(p0, p1));
1561        assert_eq!(end, Line::new(p0, p1));
1562
1563        let q = Quad::new(p0, p1, p1);
1564        let (start, end) = q.ends();
1565        assert_eq!(start, Line::new(p0, p1));
1566        assert_eq!(end, Line::new(p0, p1));
1567    }
1568
1569    #[test]
1570    fn test_split() {
1571        let q = Quad::new((0.0, 0.0), (8.0, 5.0), (4.0, 0.0));
1572        let (ql, qr) = q.split();
1573        assert_eq!((ql, qr), q.split_at(0.5));
1574        assert_eq!(ql, q.cut(0.0, 0.5));
1575        assert_eq!(qr, q.cut(0.5, 1.0));
1576
1577        let c = Cubic::new((3.0, 7.0), (2.0, 8.0), (0.0, 3.0), (6.0, 5.0));
1578        let (cl, cr) = c.split();
1579        assert_eq!((cl, cr), c.split_at(0.5));
1580        assert_eq!(cl, c.cut(0.0, 0.5));
1581        assert_eq!(cr, c.cut(0.5, 1.0));
1582    }
1583
1584    #[test]
1585    fn test_bbox() {
1586        let b0 = BBox::new(Point::new(2.0, 2.0), Point::new(4.0, 4.0));
1587        let b1 = b0.extend(Point::new(1.0, 3.0));
1588        assert!(b1.min().is_close_to(Point::new(1.0, 2.0)));
1589        assert!(b1.max().is_close_to(b0.max()));
1590        let b2 = b1.extend(Point::new(5.0, 3.0));
1591        assert!(b2.min().is_close_to(b1.min()));
1592        assert!(b2.max().is_close_to(Point::new(5.0, 4.0)));
1593        let b3 = b2.extend(Point::new(3.0, 1.0));
1594        assert!(b3.min().is_close_to(Point::new(1.0, 1.0)));
1595        assert!(b3.max().is_close_to(b2.max()));
1596        let b4 = b3.extend(Point::new(3.0, 5.0));
1597        assert!(b4.min().is_close_to(b3.min()));
1598        assert!(b4.max().is_close_to(Point::new(5.0, 5.0)));
1599
1600        let cubic = Cubic::new((106.0, 0.0), (0.0, 100.0), (382.0, 216.0), (324.0, 14.0));
1601        let bbox = cubic.bbox(None);
1602        assert_approx_eq!(bbox.x(), 87.308, 0.001);
1603        assert_approx_eq!(bbox.y(), 0.0, 0.001);
1604        assert_approx_eq!(bbox.width(), 242.724, 0.001);
1605        assert_approx_eq!(bbox.height(), 125.140, 0.001);
1606
1607        let quad = Quad::new((30.0, 90.0), (220.0, 200.0), (120.0, 50.0));
1608        let bbox = quad.bbox(None);
1609        assert_approx_eq!(bbox.x(), 30.0, 0.001);
1610        assert_approx_eq!(bbox.y(), 50.0, 0.001);
1611        assert_approx_eq!(bbox.width(), 124.483, 0.001);
1612        assert_approx_eq!(bbox.height(), 86.538, 0.001);
1613
1614        let cubic = Cubic::new((0.0, 0.0), (10.0, -3.0), (-4.0, -3.0), (6.0, 0.0));
1615        let bbox = cubic.bbox(None);
1616        assert_approx_eq!(bbox.x(), 0.0);
1617        assert_approx_eq!(bbox.y(), -2.25);
1618        assert_approx_eq!(bbox.width(), 6.0);
1619        assert_approx_eq!(bbox.height(), 2.25);
1620    }
1621
1622    fn segment_length_naive(segment: Segment) -> Scalar {
1623        let mut length = 0.0;
1624        for line in segment.flatten(Transform::identity(), 1e-5) {
1625            length += line.length();
1626        }
1627        length
1628    }
1629
1630    #[test]
1631    fn test_length() {
1632        let error = 2e-4;
1633
1634        let cubic = Cubic::new((158.0, 70.0), (210.0, 250.0), (25.0, 190.0), (219.0, 89.0));
1635        let length = cubic.length(0.0, 1.0);
1636        let length_naive = segment_length_naive(cubic.into());
1637        assert_approx_eq!(length, length_naive, length_naive * error);
1638
1639        let cubic = Cubic::new((210.0, 30.0), (209.0, 234.0), (54.0, 31.0), (118.0, 158.0));
1640        let length = cubic.length(0.0, 1.0);
1641        let length_naive = segment_length_naive(cubic.into());
1642        assert_approx_eq!(length, length_naive, length_naive * error);
1643
1644        let cubic = Cubic::new((46.0, 41.0), (210.0, 250.0), (25.0, 190.0), (219.0, 89.0));
1645        let length = cubic.length(0.0, 1.0);
1646        let length_naive = segment_length_naive(cubic.into());
1647        assert_approx_eq!(length, length_naive, length_naive * error);
1648
1649        let cubic = Cubic::new((176.0, 73.0), (109.0, 26.0), (190.0, 196.0), (209.0, 23.0));
1650        let length = cubic.length(0.0, 1.0);
1651        let length_naive = segment_length_naive(cubic.into());
1652        assert_approx_eq!(length, length_naive, length_naive * error);
1653
1654        let quad = Quad::new((139.0, 91.0), (20.0, 110.0), (221.0, 154.0));
1655        let length = quad.length(0.0, 1.0);
1656        let length_naive = segment_length_naive(quad.into());
1657        assert_approx_eq!(length, length_naive, length_naive * error);
1658
1659        let cubic = Cubic::new((40.0, 118.0), (45.0, 216.0), (190.0, 196.0), (209.0, 23.0));
1660        let length = cubic.length(0.1, 0.9);
1661        let length_naive = segment_length_naive(cubic.cut(0.1, 0.9).into());
1662        assert_approx_eq!(length, length_naive, length_naive * error);
1663    }
1664
1665    #[test]
1666    fn test_param_at_length() {
1667        let cubic = Cubic::new((40.0, 118.0), (45.0, 216.0), (190.0, 196.0), (209.0, 23.0));
1668        let length = cubic.length(0.0, 1.0);
1669        let error = length / 1000.0;
1670        for index in 0..128 {
1671            let l = length * index as Scalar / 127.0;
1672            let t = cubic.param_at_length(l, Some(error));
1673            assert_approx_eq!(cubic.length(0.0, t), l, error)
1674        }
1675    }
1676
1677    #[test]
1678    fn test_intersect() -> Result<(), SvgParserError> {
1679        let c0: Segment = "M97,177 C120,20 220,95 47,148".parse()?;
1680        let c1: Segment = "M136,82 C59,184 108,228 153,34".parse()?;
1681        // TODO: fix duplicate results
1682        for (i, p) in c0.intersect(c1, 1e-3).into_iter().enumerate() {
1683            println!("{i}: ({}, {})", p.x(), p.y());
1684        }
1685        Ok(())
1686    }
1687}