lyon_geom/
cubic_bezier.rs

1use crate::cubic_bezier_intersections::cubic_bezier_intersections_t;
2use crate::scalar::Scalar;
3use crate::segment::Segment;
4use crate::traits::Transformation;
5use crate::utils::{cubic_polynomial_roots, min_max};
6use crate::{point, Box2D, Point, Vector};
7use crate::{Line, LineEquation, LineSegment, QuadraticBezierSegment};
8use arrayvec::ArrayVec;
9
10use core::cmp::Ordering::{Equal, Greater, Less};
11use core::ops::Range;
12
13#[cfg(test)]
14use std::vec::Vec;
15
16/// A 2d curve segment defined by four points: the beginning of the segment, two control
17/// points and the end of the segment.
18///
19/// The curve is defined by equation:²
20/// ```∀ t ∈ [0..1],  P(t) = (1 - t)³ * from + 3 * (1 - t)² * t * ctrl1 + 3 * t² * (1 - t) * ctrl2 + t³ * to```
21#[derive(Copy, Clone, Debug, PartialEq)]
22#[cfg_attr(feature = "serialization", derive(Serialize, Deserialize))]
23pub struct CubicBezierSegment<S> {
24    pub from: Point<S>,
25    pub ctrl1: Point<S>,
26    pub ctrl2: Point<S>,
27    pub to: Point<S>,
28}
29
30impl<S: Scalar> CubicBezierSegment<S> {
31    /// Sample the curve at t (expecting t between 0 and 1).
32    pub fn sample(&self, t: S) -> Point<S> {
33        let t2 = t * t;
34        let t3 = t2 * t;
35        let one_t = S::ONE - t;
36        let one_t2 = one_t * one_t;
37        let one_t3 = one_t2 * one_t;
38
39        self.from * one_t3
40            + self.ctrl1.to_vector() * S::THREE * one_t2 * t
41            + self.ctrl2.to_vector() * S::THREE * one_t * t2
42            + self.to.to_vector() * t3
43    }
44
45    /// Sample the x coordinate of the curve at t (expecting t between 0 and 1).
46    pub fn x(&self, t: S) -> S {
47        let t2 = t * t;
48        let t3 = t2 * t;
49        let one_t = S::ONE - t;
50        let one_t2 = one_t * one_t;
51        let one_t3 = one_t2 * one_t;
52
53        self.from.x * one_t3
54            + self.ctrl1.x * S::THREE * one_t2 * t
55            + self.ctrl2.x * S::THREE * one_t * t2
56            + self.to.x * t3
57    }
58
59    /// Sample the y coordinate of the curve at t (expecting t between 0 and 1).
60    pub fn y(&self, t: S) -> S {
61        let t2 = t * t;
62        let t3 = t2 * t;
63        let one_t = S::ONE - t;
64        let one_t2 = one_t * one_t;
65        let one_t3 = one_t2 * one_t;
66
67        self.from.y * one_t3
68            + self.ctrl1.y * S::THREE * one_t2 * t
69            + self.ctrl2.y * S::THREE * one_t * t2
70            + self.to.y * t3
71    }
72
73    /// Return the parameter values corresponding to a given x coordinate.
74    pub fn solve_t_for_x(&self, x: S) -> ArrayVec<S, 3> {
75        let (min, max) = self.fast_bounding_range_x();
76        if min > x || max < x {
77            return ArrayVec::new();
78        }
79
80        self.parameters_for_xy_value(x, self.from.x, self.ctrl1.x, self.ctrl2.x, self.to.x)
81    }
82
83    /// Return the parameter values corresponding to a given y coordinate.
84    pub fn solve_t_for_y(&self, y: S) -> ArrayVec<S, 3> {
85        let (min, max) = self.fast_bounding_range_y();
86        if min > y || max < y {
87            return ArrayVec::new();
88        }
89
90        self.parameters_for_xy_value(y, self.from.y, self.ctrl1.y, self.ctrl2.y, self.to.y)
91    }
92
93    fn parameters_for_xy_value(
94        &self,
95        value: S,
96        from: S,
97        ctrl1: S,
98        ctrl2: S,
99        to: S,
100    ) -> ArrayVec<S, 3> {
101        let mut result = ArrayVec::new();
102
103        let a = -from + S::THREE * ctrl1 - S::THREE * ctrl2 + to;
104        let b = S::THREE * from - S::SIX * ctrl1 + S::THREE * ctrl2;
105        let c = -S::THREE * from + S::THREE * ctrl1;
106        let d = from - value;
107
108        let roots = cubic_polynomial_roots(a, b, c, d);
109        for root in roots {
110            if root > S::ZERO && root < S::ONE {
111                result.push(root);
112            }
113        }
114
115        result
116    }
117
118    #[inline]
119    fn derivative_coefficients(&self, t: S) -> (S, S, S, S) {
120        let t2 = t * t;
121        (
122            -S::THREE * t2 + S::SIX * t - S::THREE,
123            S::NINE * t2 - S::value(12.0) * t + S::THREE,
124            -S::NINE * t2 + S::SIX * t,
125            S::THREE * t2,
126        )
127    }
128
129    /// Sample the curve's derivative at t (expecting t between 0 and 1).
130    pub fn derivative(&self, t: S) -> Vector<S> {
131        let (c0, c1, c2, c3) = self.derivative_coefficients(t);
132        self.from.to_vector() * c0
133            + self.ctrl1.to_vector() * c1
134            + self.ctrl2.to_vector() * c2
135            + self.to.to_vector() * c3
136    }
137
138    /// Sample the x coordinate of the curve's derivative at t (expecting t between 0 and 1).
139    pub fn dx(&self, t: S) -> S {
140        let (c0, c1, c2, c3) = self.derivative_coefficients(t);
141        self.from.x * c0 + self.ctrl1.x * c1 + self.ctrl2.x * c2 + self.to.x * c3
142    }
143
144    /// Sample the y coordinate of the curve's derivative at t (expecting t between 0 and 1).
145    pub fn dy(&self, t: S) -> S {
146        let (c0, c1, c2, c3) = self.derivative_coefficients(t);
147        self.from.y * c0 + self.ctrl1.y * c1 + self.ctrl2.y * c2 + self.to.y * c3
148    }
149
150    /// Return the sub-curve inside a given range of t.
151    ///
152    /// This is equivalent to splitting at the range's end points.
153    pub fn split_range(&self, t_range: Range<S>) -> Self {
154        let (t0, t1) = (t_range.start, t_range.end);
155        let from = self.sample(t0);
156        let to = self.sample(t1);
157
158        let d = QuadraticBezierSegment {
159            from: (self.ctrl1 - self.from).to_point(),
160            ctrl: (self.ctrl2 - self.ctrl1).to_point(),
161            to: (self.to - self.ctrl2).to_point(),
162        };
163
164        let dt = t1 - t0;
165        let ctrl1 = from + d.sample(t0).to_vector() * dt;
166        let ctrl2 = to - d.sample(t1).to_vector() * dt;
167
168        CubicBezierSegment {
169            from,
170            ctrl1,
171            ctrl2,
172            to,
173        }
174    }
175
176    /// Split this curve into two sub-curves.
177    pub fn split(&self, t: S) -> (CubicBezierSegment<S>, CubicBezierSegment<S>) {
178        let ctrl1a = self.from + (self.ctrl1 - self.from) * t;
179        let ctrl2a = self.ctrl1 + (self.ctrl2 - self.ctrl1) * t;
180        let ctrl1aa = ctrl1a + (ctrl2a - ctrl1a) * t;
181        let ctrl3a = self.ctrl2 + (self.to - self.ctrl2) * t;
182        let ctrl2aa = ctrl2a + (ctrl3a - ctrl2a) * t;
183        let ctrl1aaa = ctrl1aa + (ctrl2aa - ctrl1aa) * t;
184
185        (
186            CubicBezierSegment {
187                from: self.from,
188                ctrl1: ctrl1a,
189                ctrl2: ctrl1aa,
190                to: ctrl1aaa,
191            },
192            CubicBezierSegment {
193                from: ctrl1aaa,
194                ctrl1: ctrl2aa,
195                ctrl2: ctrl3a,
196                to: self.to,
197            },
198        )
199    }
200
201    /// Return the curve before the split point.
202    pub fn before_split(&self, t: S) -> CubicBezierSegment<S> {
203        let ctrl1a = self.from + (self.ctrl1 - self.from) * t;
204        let ctrl2a = self.ctrl1 + (self.ctrl2 - self.ctrl1) * t;
205        let ctrl1aa = ctrl1a + (ctrl2a - ctrl1a) * t;
206        let ctrl3a = self.ctrl2 + (self.to - self.ctrl2) * t;
207        let ctrl2aa = ctrl2a + (ctrl3a - ctrl2a) * t;
208        let ctrl1aaa = ctrl1aa + (ctrl2aa - ctrl1aa) * t;
209
210        CubicBezierSegment {
211            from: self.from,
212            ctrl1: ctrl1a,
213            ctrl2: ctrl1aa,
214            to: ctrl1aaa,
215        }
216    }
217
218    /// Return the curve after the split point.
219    pub fn after_split(&self, t: S) -> CubicBezierSegment<S> {
220        let ctrl1a = self.from + (self.ctrl1 - self.from) * t;
221        let ctrl2a = self.ctrl1 + (self.ctrl2 - self.ctrl1) * t;
222        let ctrl1aa = ctrl1a + (ctrl2a - ctrl1a) * t;
223        let ctrl3a = self.ctrl2 + (self.to - self.ctrl2) * t;
224        let ctrl2aa = ctrl2a + (ctrl3a - ctrl2a) * t;
225
226        CubicBezierSegment {
227            from: ctrl1aa + (ctrl2aa - ctrl1aa) * t,
228            ctrl1: ctrl2a + (ctrl3a - ctrl2a) * t,
229            ctrl2: ctrl3a,
230            to: self.to,
231        }
232    }
233
234    #[inline]
235    pub fn baseline(&self) -> LineSegment<S> {
236        LineSegment {
237            from: self.from,
238            to: self.to,
239        }
240    }
241
242    /// Returns true if the curve can be approximated with a single line segment, given
243    /// a tolerance threshold.
244    pub fn is_linear(&self, tolerance: S) -> bool {
245        // Similar to Line::square_distance_to_point, except we keep
246        // the sign of c1 and c2 to compute tighter upper bounds as we
247        // do in fat_line_min_max.
248        let baseline = self.to - self.from;
249        let v1 = self.ctrl1 - self.from;
250        let v2 = self.ctrl2 - self.from;
251        let c1 = baseline.cross(v1);
252        let c2 = baseline.cross(v2);
253        // TODO: it would be faster to multiply the threshold with baseline_len2
254        // instead of dividing d1 and d2, but it changes the behavior when the
255        // baseline length is zero in ways that breaks some of the cubic intersection
256        // tests.
257        let inv_baseline_len2 = S::ONE / baseline.square_length();
258        let d1 = (c1 * c1) * inv_baseline_len2;
259        let d2 = (c2 * c2) * inv_baseline_len2;
260
261        let factor = if (c1 * c2) > S::ZERO {
262            S::THREE / S::FOUR
263        } else {
264            S::FOUR / S::NINE
265        };
266
267        let f2 = factor * factor;
268        let threshold = tolerance * tolerance;
269
270        d1 * f2 <= threshold && d2 * f2 <= threshold
271    }
272
273    /// Returns whether the curve can be approximated with a single point, given
274    /// a tolerance threshold.
275    pub(crate) fn is_a_point(&self, tolerance: S) -> bool {
276        let tolerance_squared = tolerance * tolerance;
277        // Use <= so that tolerance can be zero.
278        (self.from - self.to).square_length() <= tolerance_squared
279            && (self.from - self.ctrl1).square_length() <= tolerance_squared
280            && (self.to - self.ctrl2).square_length() <= tolerance_squared
281    }
282
283    /// Computes the signed distances (min <= 0 and max >= 0) from the baseline of this
284    /// curve to its two "fat line" boundary lines.
285    ///
286    /// A fat line is two conservative lines between which the segment
287    /// is fully contained.
288    pub(crate) fn fat_line_min_max(&self) -> (S, S) {
289        let baseline = self.baseline().to_line().equation();
290        let (d1, d2) = min_max(
291            baseline.signed_distance_to_point(&self.ctrl1),
292            baseline.signed_distance_to_point(&self.ctrl2),
293        );
294
295        let factor = if (d1 * d2) > S::ZERO {
296            S::THREE / S::FOUR
297        } else {
298            S::FOUR / S::NINE
299        };
300
301        let d_min = factor * S::min(d1, S::ZERO);
302        let d_max = factor * S::max(d2, S::ZERO);
303
304        (d_min, d_max)
305    }
306
307    /// Computes a "fat line" of this segment.
308    ///
309    /// A fat line is two conservative lines between which the segment
310    /// is fully contained.
311    pub fn fat_line(&self) -> (LineEquation<S>, LineEquation<S>) {
312        let baseline = self.baseline().to_line().equation();
313        let (d1, d2) = self.fat_line_min_max();
314
315        (baseline.offset(d1), baseline.offset(d2))
316    }
317
318    /// Applies the transform to this curve and returns the results.
319    #[inline]
320    pub fn transformed<T: Transformation<S>>(&self, transform: &T) -> Self {
321        CubicBezierSegment {
322            from: transform.transform_point(self.from),
323            ctrl1: transform.transform_point(self.ctrl1),
324            ctrl2: transform.transform_point(self.ctrl2),
325            to: transform.transform_point(self.to),
326        }
327    }
328
329    /// Swap the beginning and the end of the segment.
330    pub fn flip(&self) -> Self {
331        CubicBezierSegment {
332            from: self.to,
333            ctrl1: self.ctrl2,
334            ctrl2: self.ctrl1,
335            to: self.from,
336        }
337    }
338
339    /// Approximate the curve with a single quadratic bézier segment.
340    ///
341    /// This is terrible as a general approximation but works if the cubic
342    /// curve does not have inflection points and is "flat" enough. Typically
343    /// usable after subdividing the curve a few times.
344    pub fn to_quadratic(&self) -> QuadraticBezierSegment<S> {
345        let c1 = (self.ctrl1 * S::THREE - self.from) * S::HALF;
346        let c2 = (self.ctrl2 * S::THREE - self.to) * S::HALF;
347        QuadraticBezierSegment {
348            from: self.from,
349            ctrl: ((c1 + c2) * S::HALF).to_point(),
350            to: self.to,
351        }
352    }
353
354    /// Evaluates an upper bound on the maximum distance between the curve
355    /// and its quadratic approximation obtained using `to_quadratic`.
356    pub fn to_quadratic_error(&self) -> S {
357        // See http://caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
358        S::sqrt(S::THREE) / S::value(36.0)
359            * ((self.to - self.ctrl2 * S::THREE) + (self.ctrl1 * S::THREE - self.from)).length()
360    }
361
362    /// Returns true if the curve can be safely approximated with a single quadratic bézier
363    /// segment given the provided tolerance threshold.
364    ///
365    /// Equivalent to comparing `to_quadratic_error` with the tolerance threshold, avoiding
366    /// the cost of two square roots.
367    pub fn is_quadratic(&self, tolerance: S) -> bool {
368        S::THREE / S::value(1296.0)
369            * ((self.to - self.ctrl2 * S::THREE) + (self.ctrl1 * S::THREE - self.from))
370                .square_length()
371            <= tolerance * tolerance
372    }
373
374    /// Computes the number of quadratic bézier segments required to approximate this cubic curve
375    /// given a tolerance threshold.
376    ///
377    /// Derived by Raph Levien from section 10.6 of Sedeberg's CAGD notes
378    /// <https://scholarsarchive.byu.edu/cgi/viewcontent.cgi?article=1000&context=facpub#section.10.6>
379    /// and the error metric from the caffein owl blog post <http://caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html>
380    pub fn num_quadratics(&self, tolerance: S) -> u32 {
381        self.num_quadratics_impl(tolerance)
382    }
383
384    fn num_quadratics_impl(&self, tolerance: S) -> u32 {
385        debug_assert!(tolerance > S::ZERO);
386
387        let x = self.from.x - S::THREE * self.ctrl1.x + S::THREE * self.ctrl2.x - self.to.x;
388        let y = self.from.y - S::THREE * self.ctrl1.y + S::THREE * self.ctrl2.y - self.to.y;
389
390        let err = ((x * x + y * y) / (S::value(432.0) * tolerance * tolerance)).to_f32().unwrap_or(1.0);
391
392        // Avoid computing powf(1/6).ceil() using a lookup table that contains
393        // i^6 for i in 1..25.
394        const MAX_QUADS: usize = 16;
395        const LUT: [f32; MAX_QUADS] = [
396            1.0, 64.0, 729.0, 4096.0, 15625.0, 46656.0, 117649.0, 262144.0, 531441.0, 1000000.0,
397            1771561.0, 2985984.0, 4826809.0, 7529536.0, 11390625.0, 16777216.0,
398        ];
399
400        if err <= 16777216.0 {
401            #[allow(clippy::needless_range_loop)]
402            for i in 0..MAX_QUADS {
403                if err <= LUT[i] {
404                    return (i + 1) as u32;
405                }
406            }
407        }
408
409        // If the number of quads does not fit in the LUT, fall back to the
410        // expensive computation.
411        S::from(err).unwrap().powf(S::ONE / S::SIX).ceil().max(S::ONE).to_u32().unwrap_or(1)
412    }
413
414    /// Returns the flattened representation of the curve as an iterator, starting *after* the
415    /// current point.
416    pub fn flattened(&self, tolerance: S) -> Flattened<S> {
417        Flattened::new(self, tolerance)
418    }
419
420    /// Invokes a callback for each monotonic part of the segment.
421    pub fn for_each_monotonic_range<F>(&self, cb: &mut F)
422    where
423        F: FnMut(Range<S>),
424    {
425        let mut extrema: ArrayVec<S, 4> = ArrayVec::new();
426        self.for_each_local_x_extremum_t(&mut |t| extrema.push(t));
427        self.for_each_local_y_extremum_t(&mut |t| extrema.push(t));
428        extrema.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
429
430        let mut t0 = S::ZERO;
431        for &t in &extrema {
432            if t != t0 {
433                cb(t0..t);
434                t0 = t;
435            }
436        }
437
438        cb(t0..S::ONE);
439    }
440
441    /// Invokes a callback for each monotonic part of the segment.
442    pub fn for_each_monotonic<F>(&self, cb: &mut F)
443    where
444        F: FnMut(&CubicBezierSegment<S>),
445    {
446        self.for_each_monotonic_range(&mut |range| {
447            let mut sub = self.split_range(range);
448            // Due to finite precision the split may actually result in sub-curves
449            // that are almost but not-quite monotonic. Make sure they actually are.
450            let min_x = sub.from.x.min(sub.to.x);
451            let max_x = sub.from.x.max(sub.to.x);
452            let min_y = sub.from.y.min(sub.to.y);
453            let max_y = sub.from.y.max(sub.to.y);
454            sub.ctrl1.x = sub.ctrl1.x.max(min_x).min(max_x);
455            sub.ctrl1.y = sub.ctrl1.y.max(min_y).min(max_y);
456            sub.ctrl2.x = sub.ctrl2.x.max(min_x).min(max_x);
457            sub.ctrl2.y = sub.ctrl2.y.max(min_y).min(max_y);
458            cb(&sub);
459        });
460    }
461
462    /// Invokes a callback for each y-monotonic part of the segment.
463    pub fn for_each_y_monotonic_range<F>(&self, cb: &mut F)
464    where
465        F: FnMut(Range<S>),
466    {
467        let mut t0 = S::ZERO;
468        self.for_each_local_y_extremum_t(&mut |t| {
469            cb(t0..t);
470            t0 = t;
471        });
472
473        cb(t0..S::ONE);
474    }
475
476    /// Invokes a callback for each y-monotonic part of the segment.
477    pub fn for_each_y_monotonic<F>(&self, cb: &mut F)
478    where
479        F: FnMut(&CubicBezierSegment<S>),
480    {
481        self.for_each_y_monotonic_range(&mut |range| {
482            let mut sub = self.split_range(range);
483            // Due to finite precision the split may actually result in sub-curves
484            // that are almost but not-quite monotonic. Make sure they actually are.
485            let min_y = sub.from.y.min(sub.to.y);
486            let max_y = sub.from.y.max(sub.to.y);
487            sub.ctrl1.y = sub.ctrl1.y.max(min_y).min(max_y);
488            sub.ctrl2.y = sub.ctrl2.y.max(min_y).min(max_y);
489            cb(&sub);
490        });
491    }
492
493    /// Invokes a callback for each x-monotonic part of the segment.
494    pub fn for_each_x_monotonic_range<F>(&self, cb: &mut F)
495    where
496        F: FnMut(Range<S>),
497    {
498        let mut t0 = S::ZERO;
499        self.for_each_local_x_extremum_t(&mut |t| {
500            cb(t0..t);
501            t0 = t;
502        });
503
504        cb(t0..S::ONE);
505    }
506
507    /// Invokes a callback for each x-monotonic part of the segment.
508    pub fn for_each_x_monotonic<F>(&self, cb: &mut F)
509    where
510        F: FnMut(&CubicBezierSegment<S>),
511    {
512        self.for_each_x_monotonic_range(&mut |range| {
513            let mut sub = self.split_range(range);
514            // Due to finite precision the split may actually result in sub-curves
515            // that are almost but not-quite monotonic. Make sure they actually are.
516            let min_x = sub.from.x.min(sub.to.x);
517            let max_x = sub.from.x.max(sub.to.x);
518            sub.ctrl1.x = sub.ctrl1.x.max(min_x).min(max_x);
519            sub.ctrl2.x = sub.ctrl2.x.max(min_x).min(max_x);
520            cb(&sub);
521        });
522    }
523
524    /// Approximates the cubic bézier curve with sequence of quadratic ones,
525    /// invoking a callback at each step.
526    pub fn for_each_quadratic_bezier<F>(&self, tolerance: S, cb: &mut F)
527    where
528        F: FnMut(&QuadraticBezierSegment<S>),
529    {
530        self.for_each_quadratic_bezier_with_t(tolerance, &mut |quad, _range| cb(quad));
531    }
532
533    /// Approximates the cubic bézier curve with sequence of quadratic ones,
534    /// invoking a callback at each step.
535    pub fn for_each_quadratic_bezier_with_t<F>(&self, tolerance: S, cb: &mut F)
536    where
537        F: FnMut(&QuadraticBezierSegment<S>, Range<S>),
538    {
539        debug_assert!(tolerance >= S::EPSILON * S::EPSILON);
540
541        let num_quadratics = S::from(self.num_quadratics_impl(tolerance)).unwrap();
542        let step = S::ONE / num_quadratics;
543        let n = num_quadratics.to_u32().unwrap_or(1);
544        let mut t0 = S::ZERO;
545        for _ in 0..(n - 1) {
546            let t1 = t0 + step;
547
548            let quad = self.split_range(t0..t1).to_quadratic();
549            cb(&quad, t0..t1);
550
551            t0 = t1;
552        }
553
554        // Do the last step manually to make sure we finish at t = 1.0 exactly.
555        let quad = self.split_range(t0..S::ONE).to_quadratic();
556        cb(&quad, t0..S::ONE)
557    }
558
559    /// Approximates the curve with sequence of line segments.
560    ///
561    /// The `tolerance` parameter defines the maximum distance between the curve and
562    /// its approximation.
563    pub fn for_each_flattened<F: FnMut(&LineSegment<S>)>(&self, tolerance: S, callback: &mut F) {
564        debug_assert!(tolerance >= S::EPSILON * S::EPSILON);
565
566        let n = flattened_segments_wang(self, tolerance);
567        let step = S::ONE / n;
568        let mut t = step;
569
570        let curve = self.polynomial_form();
571
572        let n = n.to_u32().unwrap_or(1) - 1;
573        let mut from = self.from;
574        for _ in 0..n {
575            let to = curve.sample(t);
576            callback(&LineSegment { from, to });
577            from = to;
578            t += step;
579        }
580
581        callback(&LineSegment { from, to: self.to });
582    }
583
584    /// Approximates the curve with sequence of line segments.
585    ///
586    /// The `tolerance` parameter defines the maximum distance between the curve and
587    /// its approximation.
588    ///
589    /// The end of the t parameter range at the final segment is guaranteed to be equal to `1.0`.
590    pub fn for_each_flattened_with_t<F: FnMut(&LineSegment<S>, Range<S>)>(
591        &self,
592        tolerance: S,
593        callback: &mut F,
594    ) {
595        debug_assert!(tolerance >= S::EPSILON * S::EPSILON);
596
597        let n = flattened_segments_wang(self, tolerance);
598        let step = S::ONE / n;
599        let mut t0 = S::ZERO;
600
601        let curve = self.polynomial_form();
602
603        let n = n.to_u32().unwrap_or(1) - 1;
604        let mut from = self.from;
605        for _ in 0..n {
606            let t1 = t0 + step;
607            let to = curve.sample(t1);
608            callback(&LineSegment { from, to }, t0..t1);
609            from = to;
610            t0 = t1;
611        }
612
613        callback(&LineSegment { from, to: self.to }, t0..S::ONE);
614    }
615
616    /// Compute the length of the segment using a flattened approximation.
617    pub fn approximate_length(&self, tolerance: S) -> S {
618        let mut length = S::ZERO;
619
620        self.for_each_quadratic_bezier(tolerance, &mut |quad| {
621            length += quad.length();
622        });
623
624        length
625    }
626
627    /// Invokes a callback at each inflection point if any.
628    pub fn for_each_inflection_t<F>(&self, cb: &mut F)
629    where
630        F: FnMut(S),
631    {
632        // Find inflection points.
633        // See www.faculty.idc.ac.il/arik/quality/appendixa.html for an explanation
634        // of this approach.
635        let pa = self.ctrl1 - self.from;
636        let pb = self.ctrl2.to_vector() - (self.ctrl1.to_vector() * S::TWO) + self.from.to_vector();
637        let pc = self.to.to_vector() - (self.ctrl2.to_vector() * S::THREE)
638            + (self.ctrl1.to_vector() * S::THREE)
639            - self.from.to_vector();
640
641        let a = pb.cross(pc);
642        let b = pa.cross(pc);
643        let c = pa.cross(pb);
644
645        if S::abs(a) < S::EPSILON {
646            // Not a quadratic equation.
647            if S::abs(b) < S::EPSILON {
648                // Instead of a linear acceleration change we have a constant
649                // acceleration change. This means the equation has no solution
650                // and there are no inflection points, unless the constant is 0.
651                // In that case the curve is a straight line, essentially that means
652                // the easiest way to deal with is is by saying there's an inflection
653                // point at t == 0. The inflection point approximation range found will
654                // automatically extend into infinity.
655                if S::abs(c) < S::EPSILON {
656                    cb(S::ZERO);
657                }
658            } else {
659                let t = -c / b;
660                if in_range(t) {
661                    cb(t);
662                }
663            }
664
665            return;
666        }
667
668        fn in_range<S: Scalar>(t: S) -> bool {
669            t >= S::ZERO && t < S::ONE
670        }
671
672        let discriminant = b * b - S::FOUR * a * c;
673
674        if discriminant < S::ZERO {
675            return;
676        }
677
678        if discriminant < S::EPSILON {
679            let t = -b / (S::TWO * a);
680
681            if in_range(t) {
682                cb(t);
683            }
684
685            return;
686        }
687
688        // This code is derived from https://www2.units.it/ipl/students_area/imm2/files/Numerical_Recipes.pdf page 184.
689        // Computing the roots this way avoids precision issues when a, c or both are small.
690        let discriminant_sqrt = S::sqrt(discriminant);
691        let sign_b = if b >= S::ZERO { S::ONE } else { -S::ONE };
692        let q = -S::HALF * (b + sign_b * discriminant_sqrt);
693        let mut first_inflection = q / a;
694        let mut second_inflection = c / q;
695
696        if first_inflection > second_inflection {
697            core::mem::swap(&mut first_inflection, &mut second_inflection);
698        }
699
700        if in_range(first_inflection) {
701            cb(first_inflection);
702        }
703
704        if in_range(second_inflection) {
705            cb(second_inflection);
706        }
707    }
708
709    /// Return local x extrema or None if this curve is monotonic.
710    ///
711    /// This returns the advancements along the curve, not the actual x position.
712    pub fn for_each_local_x_extremum_t<F>(&self, cb: &mut F)
713    where
714        F: FnMut(S),
715    {
716        Self::for_each_local_extremum(self.from.x, self.ctrl1.x, self.ctrl2.x, self.to.x, cb)
717    }
718
719    /// Return local y extrema or None if this curve is monotonic.
720    ///
721    /// This returns the advancements along the curve, not the actual y position.
722    pub fn for_each_local_y_extremum_t<F>(&self, cb: &mut F)
723    where
724        F: FnMut(S),
725    {
726        Self::for_each_local_extremum(self.from.y, self.ctrl1.y, self.ctrl2.y, self.to.y, cb)
727    }
728
729    fn for_each_local_extremum<F>(p0: S, p1: S, p2: S, p3: S, cb: &mut F)
730    where
731        F: FnMut(S),
732    {
733        // See www.faculty.idc.ac.il/arik/quality/appendixa.html for an explanation
734        // The derivative of a cubic bezier curve is a curve representing a second degree polynomial function
735        // f(x) = a * x² + b * x + c such as :
736
737        let a = S::THREE * (p3 + S::THREE * (p1 - p2) - p0);
738        let b = S::SIX * (p2 - S::TWO * p1 + p0);
739        let c = S::THREE * (p1 - p0);
740
741        fn in_range<S: Scalar>(t: S) -> bool {
742            t > S::ZERO && t < S::ONE
743        }
744
745        // If the derivative is a linear function
746        if a == S::ZERO {
747            if b != S::ZERO {
748                let t = -c / b;
749                if in_range(t) {
750                    cb(t);
751                }
752            }
753            return;
754        }
755
756        let discriminant = b * b - S::FOUR * a * c;
757
758        // There is no Real solution for the equation
759        if discriminant < S::ZERO {
760            return;
761        }
762
763        // There is one Real solution for the equation
764        if discriminant == S::ZERO {
765            let t = -b / (S::TWO * a);
766            if in_range(t) {
767                cb(t);
768            }
769            return;
770        }
771
772        // There are two Real solutions for the equation
773        let discriminant_sqrt = discriminant.sqrt();
774
775        let mut first_extremum = (-b - discriminant_sqrt) / (S::TWO * a);
776        let mut second_extremum = (-b + discriminant_sqrt) / (S::TWO * a);
777        if first_extremum > second_extremum {
778            core::mem::swap(&mut first_extremum, &mut second_extremum);
779        }
780
781        if in_range(first_extremum) {
782            cb(first_extremum);
783        }
784
785        if in_range(second_extremum) {
786            cb(second_extremum);
787        }
788    }
789
790    /// Find the advancement of the y-most position in the curve.
791    ///
792    /// This returns the advancement along the curve, not the actual y position.
793    pub fn y_maximum_t(&self) -> S {
794        let mut max_t = S::ZERO;
795        let mut max_y = self.from.y;
796        if self.to.y > max_y {
797            max_t = S::ONE;
798            max_y = self.to.y;
799        }
800        self.for_each_local_y_extremum_t(&mut |t| {
801            let y = self.y(t);
802            if y > max_y {
803                max_t = t;
804                max_y = y;
805            }
806        });
807
808        max_t
809    }
810
811    /// Find the advancement of the y-least position in the curve.
812    ///
813    /// This returns the advancement along the curve, not the actual y position.
814    pub fn y_minimum_t(&self) -> S {
815        let mut min_t = S::ZERO;
816        let mut min_y = self.from.y;
817        if self.to.y < min_y {
818            min_t = S::ONE;
819            min_y = self.to.y;
820        }
821        self.for_each_local_y_extremum_t(&mut |t| {
822            let y = self.y(t);
823            if y < min_y {
824                min_t = t;
825                min_y = y;
826            }
827        });
828
829        min_t
830    }
831
832    /// Find the advancement of the x-most position in the curve.
833    ///
834    /// This returns the advancement along the curve, not the actual x position.
835    pub fn x_maximum_t(&self) -> S {
836        let mut max_t = S::ZERO;
837        let mut max_x = self.from.x;
838        if self.to.x > max_x {
839            max_t = S::ONE;
840            max_x = self.to.x;
841        }
842        self.for_each_local_x_extremum_t(&mut |t| {
843            let x = self.x(t);
844            if x > max_x {
845                max_t = t;
846                max_x = x;
847            }
848        });
849
850        max_t
851    }
852
853    /// Find the x-least position in the curve.
854    pub fn x_minimum_t(&self) -> S {
855        let mut min_t = S::ZERO;
856        let mut min_x = self.from.x;
857        if self.to.x < min_x {
858            min_t = S::ONE;
859            min_x = self.to.x;
860        }
861        self.for_each_local_x_extremum_t(&mut |t| {
862            let x = self.x(t);
863            if x < min_x {
864                min_t = t;
865                min_x = x;
866            }
867        });
868
869        min_t
870    }
871
872    /// Returns a conservative rectangle the curve is contained in.
873    ///
874    /// This method is faster than `bounding_box` but more conservative.
875    pub fn fast_bounding_box(&self) -> Box2D<S> {
876        let (min_x, max_x) = self.fast_bounding_range_x();
877        let (min_y, max_y) = self.fast_bounding_range_y();
878
879        Box2D {
880            min: point(min_x, min_y),
881            max: point(max_x, max_y),
882        }
883    }
884
885    /// Returns a conservative range of x that contains this curve.
886    #[inline]
887    pub fn fast_bounding_range_x(&self) -> (S, S) {
888        let min_x = self
889            .from
890            .x
891            .min(self.ctrl1.x)
892            .min(self.ctrl2.x)
893            .min(self.to.x);
894        let max_x = self
895            .from
896            .x
897            .max(self.ctrl1.x)
898            .max(self.ctrl2.x)
899            .max(self.to.x);
900
901        (min_x, max_x)
902    }
903
904    /// Returns a conservative range of y that contains this curve.
905    #[inline]
906    pub fn fast_bounding_range_y(&self) -> (S, S) {
907        let min_y = self
908            .from
909            .y
910            .min(self.ctrl1.y)
911            .min(self.ctrl2.y)
912            .min(self.to.y);
913        let max_y = self
914            .from
915            .y
916            .max(self.ctrl1.y)
917            .max(self.ctrl2.y)
918            .max(self.to.y);
919
920        (min_y, max_y)
921    }
922
923    /// Returns a conservative rectangle that contains the curve.
924    #[inline]
925    pub fn bounding_box(&self) -> Box2D<S> {
926        let (min_x, max_x) = self.bounding_range_x();
927        let (min_y, max_y) = self.bounding_range_y();
928
929        Box2D {
930            min: point(min_x, min_y),
931            max: point(max_x, max_y),
932        }
933    }
934
935    /// Returns the smallest range of x that contains this curve.
936    #[inline]
937    pub fn bounding_range_x(&self) -> (S, S) {
938        let min_x = self.x(self.x_minimum_t());
939        let max_x = self.x(self.x_maximum_t());
940
941        (min_x, max_x)
942    }
943
944    /// Returns the smallest range of y that contains this curve.
945    #[inline]
946    pub fn bounding_range_y(&self) -> (S, S) {
947        let min_y = self.y(self.y_minimum_t());
948        let max_y = self.y(self.y_maximum_t());
949
950        (min_y, max_y)
951    }
952
953    /// Returns whether this segment is monotonic on the x axis.
954    pub fn is_x_monotonic(&self) -> bool {
955        let mut found = false;
956        self.for_each_local_x_extremum_t(&mut |_| {
957            found = true;
958        });
959        !found
960    }
961
962    /// Returns whether this segment is monotonic on the y axis.
963    pub fn is_y_monotonic(&self) -> bool {
964        let mut found = false;
965        self.for_each_local_y_extremum_t(&mut |_| {
966            found = true;
967        });
968        !found
969    }
970
971    /// Returns whether this segment is fully monotonic.
972    pub fn is_monotonic(&self) -> bool {
973        self.is_x_monotonic() && self.is_y_monotonic()
974    }
975
976    /// Computes the intersections (if any) between this segment and another one.
977    ///
978    /// The result is provided in the form of the `t` parameters of each point along the curves. To
979    /// get the intersection points, sample the curves at the corresponding values.
980    ///
981    /// Returns endpoint intersections where an endpoint intersects the interior of the other curve,
982    /// but not endpoint/endpoint intersections.
983    ///
984    /// Returns no intersections if either curve is a point.
985    pub fn cubic_intersections_t(&self, curve: &CubicBezierSegment<S>) -> ArrayVec<(S, S), 9> {
986        cubic_bezier_intersections_t(self, curve)
987    }
988
989    /// Computes the intersection points (if any) between this segment and another one.
990    pub fn cubic_intersections(&self, curve: &CubicBezierSegment<S>) -> ArrayVec<Point<S>, 9> {
991        let intersections = self.cubic_intersections_t(curve);
992
993        let mut result_with_repeats = ArrayVec::<_, 9>::new();
994        for (t, _) in intersections {
995            result_with_repeats.push(self.sample(t));
996        }
997
998        // We can have up to nine "repeated" values here (for example: two lines, each of which
999        // overlaps itself 3 times, intersecting in their 3-fold overlaps). We make an effort to
1000        // dedupe the results, but that's hindered by not having predictable control over how far
1001        // the repeated intersections can be from each other (and then by the fact that true
1002        // intersections can be arbitrarily close), so the results will never be perfect.
1003
1004        let pair_cmp = |s: &Point<S>, t: &Point<S>| {
1005            if s.x < t.x || (s.x == t.x && s.y < t.y) {
1006                Less
1007            } else if s.x == t.x && s.y == t.y {
1008                Equal
1009            } else {
1010                Greater
1011            }
1012        };
1013        result_with_repeats.sort_unstable_by(pair_cmp);
1014        if result_with_repeats.len() <= 1 {
1015            return result_with_repeats;
1016        }
1017
1018        #[inline]
1019        fn dist_sq<S: Scalar>(p1: &Point<S>, p2: &Point<S>) -> S {
1020            (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y)
1021        }
1022
1023        let epsilon_squared = S::EPSILON * S::EPSILON;
1024        let mut result = ArrayVec::new();
1025        let mut reference_intersection = &result_with_repeats[0];
1026        result.push(*reference_intersection);
1027        for i in 1..result_with_repeats.len() {
1028            let intersection = &result_with_repeats[i];
1029            if dist_sq(reference_intersection, intersection) < epsilon_squared {
1030                continue;
1031            } else {
1032                result.push(*intersection);
1033                reference_intersection = intersection;
1034            }
1035        }
1036
1037        result
1038    }
1039
1040    /// Computes the intersections (if any) between this segment a quadratic bézier segment.
1041    ///
1042    /// The result is provided in the form of the `t` parameters of each point along the curves. To
1043    /// get the intersection points, sample the curves at the corresponding values.
1044    ///
1045    /// Returns endpoint intersections where an endpoint intersects the interior of the other curve,
1046    /// but not endpoint/endpoint intersections.
1047    ///
1048    /// Returns no intersections if either curve is a point.
1049    pub fn quadratic_intersections_t(
1050        &self,
1051        curve: &QuadraticBezierSegment<S>,
1052    ) -> ArrayVec<(S, S), 9> {
1053        self.cubic_intersections_t(&curve.to_cubic())
1054    }
1055
1056    /// Computes the intersection points (if any) between this segment and a quadratic bézier segment.
1057    pub fn quadratic_intersections(
1058        &self,
1059        curve: &QuadraticBezierSegment<S>,
1060    ) -> ArrayVec<Point<S>, 9> {
1061        self.cubic_intersections(&curve.to_cubic())
1062    }
1063
1064    /// Computes the intersections (if any) between this segment and a line.
1065    ///
1066    /// The result is provided in the form of the `t` parameters of each
1067    /// point along curve. To get the intersection points, sample the curve
1068    /// at the corresponding values.
1069    pub fn line_intersections_t(&self, line: &Line<S>) -> ArrayVec<S, 3> {
1070        if line.vector.square_length() < S::EPSILON {
1071            return ArrayVec::new();
1072        }
1073
1074        let from = self.from.to_vector();
1075        let ctrl1 = self.ctrl1.to_vector();
1076        let ctrl2 = self.ctrl2.to_vector();
1077        let to = self.to.to_vector();
1078
1079        let p1 = to - from + (ctrl1 - ctrl2) * S::THREE;
1080        let p2 = from * S::THREE + (ctrl2 - ctrl1 * S::TWO) * S::THREE;
1081        let p3 = (ctrl1 - from) * S::THREE;
1082        let p4 = from;
1083
1084        let c = line.point.y * line.vector.x - line.point.x * line.vector.y;
1085
1086        let roots = cubic_polynomial_roots(
1087            line.vector.y * p1.x - line.vector.x * p1.y,
1088            line.vector.y * p2.x - line.vector.x * p2.y,
1089            line.vector.y * p3.x - line.vector.x * p3.y,
1090            line.vector.y * p4.x - line.vector.x * p4.y + c,
1091        );
1092
1093        let mut result = ArrayVec::new();
1094
1095        for root in roots {
1096            if root >= S::ZERO && root <= S::ONE {
1097                result.push(root);
1098            }
1099        }
1100
1101        // TODO: sort the intersections?
1102
1103        result
1104    }
1105
1106    /// Computes the intersection points (if any) between this segment and a line.
1107    pub fn line_intersections(&self, line: &Line<S>) -> ArrayVec<Point<S>, 3> {
1108        let intersections = self.line_intersections_t(line);
1109
1110        let mut result = ArrayVec::new();
1111        for t in intersections {
1112            result.push(self.sample(t));
1113        }
1114
1115        result
1116    }
1117
1118    /// Computes the intersections (if any) between this segment and a line segment.
1119    ///
1120    /// The result is provided in the form of the `t` parameters of each
1121    /// point along curve and segment. To get the intersection points, sample
1122    /// the segments at the corresponding values.
1123    pub fn line_segment_intersections_t(&self, segment: &LineSegment<S>) -> ArrayVec<(S, S), 3> {
1124        if !self
1125            .fast_bounding_box()
1126            .inflate(S::EPSILON, S::EPSILON)
1127            .intersects(&segment.bounding_box().inflate(S::EPSILON, S::EPSILON))
1128        {
1129            return ArrayVec::new();
1130        }
1131
1132        let intersections = self.line_intersections_t(&segment.to_line());
1133
1134        let mut result = ArrayVec::new();
1135        if intersections.is_empty() {
1136            return result;
1137        }
1138
1139        let seg_is_mostly_vertical =
1140            S::abs(segment.from.y - segment.to.y) >= S::abs(segment.from.x - segment.to.x);
1141        let (seg_long_axis_min, seg_long_axis_max) = if seg_is_mostly_vertical {
1142            segment.bounding_range_y()
1143        } else {
1144            segment.bounding_range_x()
1145        };
1146
1147        for t in intersections {
1148            let intersection_xy = if seg_is_mostly_vertical {
1149                self.y(t)
1150            } else {
1151                self.x(t)
1152            };
1153            if intersection_xy >= seg_long_axis_min && intersection_xy <= seg_long_axis_max {
1154                let t2 = (self.sample(t) - segment.from).length() / segment.length();
1155                // Don't take intersections that are on endpoints of both curves at the same time.
1156                if (t != S::ZERO && t != S::ONE) || (t2 != S::ZERO && t2 != S::ONE) {
1157                    result.push((t, t2));
1158                }
1159            }
1160        }
1161
1162        result
1163    }
1164
1165    #[inline]
1166    pub fn from(&self) -> Point<S> {
1167        self.from
1168    }
1169
1170    #[inline]
1171    pub fn to(&self) -> Point<S> {
1172        self.to
1173    }
1174
1175    pub fn line_segment_intersections(&self, segment: &LineSegment<S>) -> ArrayVec<Point<S>, 3> {
1176        let intersections = self.line_segment_intersections_t(segment);
1177
1178        let mut result = ArrayVec::new();
1179        for (t, _) in intersections {
1180            result.push(self.sample(t));
1181        }
1182
1183        result
1184    }
1185
1186    fn baseline_projection(&self, t: S) -> S {
1187        // See https://pomax.github.io/bezierinfo/#abc
1188        // We are computing the interpolation factor between
1189        // `from` and `to` to get the position of C.
1190        let one_t = S::ONE - t;
1191        let one_t3 = one_t * one_t * one_t;
1192        let t3 = t * t * t;
1193
1194        t3 / (t3 + one_t3)
1195    }
1196
1197    fn abc_ratio(&self, t: S) -> S {
1198        // See https://pomax.github.io/bezierinfo/#abc
1199        let one_t = S::ONE - t;
1200        let one_t3 = one_t * one_t * one_t;
1201        let t3 = t * t * t;
1202
1203        ((t3 + one_t3 - S::ONE) / (t3 + one_t3)).abs()
1204    }
1205
1206    // Returns a quadratic bézier curve built by dragging this curve's point at `t`
1207    // to a new position, without moving the endpoints.
1208    //
1209    // The relative effect on control points is chosen to give a similar "feel" to
1210    // most vector graphics editors: dragging from near the first endpoint will affect
1211    // the first control point more than the second control point, etc.
1212    pub fn drag(&self, t: S, new_position: Point<S>) -> Self {
1213        // A lot of tweaking could go into making the weight feel as natural as possible.
1214        let min = S::value(0.1);
1215        let max = S::value(0.9);
1216        let weight = if t < min {
1217            S::ZERO
1218        } else if t > max {
1219            S::ONE
1220        } else {
1221            (t - min) / (max - min)
1222        };
1223
1224        self.drag_with_weight(t, new_position, weight)
1225    }
1226
1227    // Returns a quadratic bézier curve built by dragging this curve's point at `t`
1228    // to a new position, without moving the endpoints.
1229    //
1230    // The provided weight specifies the relative effect on control points.
1231    //  - with `weight = 0.5`, `ctrl1` and `ctrl2` are equally affected,
1232    //  - with `weight = 0.0`, only `ctrl1` is affected,
1233    //  - with `weight = 1.0`, only `ctrl2` is affected,
1234    //  - etc.
1235    pub fn drag_with_weight(&self, t: S, new_position: Point<S>, weight: S) -> Self {
1236        // See https://pomax.github.io/bezierinfo/#abc
1237        //
1238        //   From-----------Ctrl1
1239        //    |               \ d1     \
1240        //    C-------P--------A        \  d12
1241        //    |                 \d2      \
1242        //    |                  \        \
1243        //    To-----------------Ctrl2
1244        //
1245        // The ABC relation means we can place the new control points however we like
1246        // as long as the ratios CA/CP, d1/d12 and d2/d12 remain constant.
1247        //
1248        // we use the weight to guide our decisions. A weight of 0.5 would be a uniform
1249        // displacement (d1 and d2 do not change and both control points are moved by the
1250        // same amount).
1251        // The approach is to use the weight interpolate the most constrained control point
1252        // between it's old position and the position it would have with uniform displacement.
1253        // then we determine the position of the least constrained control point such that
1254        // the ratios mentioned earlier remain constant.
1255
1256        let c = self.from.lerp(self.to, self.baseline_projection(t));
1257        let cacp_ratio = self.abc_ratio(t);
1258
1259        let old_pos = self.sample(t);
1260        // Construct A before and after drag using the constance ca/cp ratio
1261        let old_a = old_pos + (old_pos - c) / cacp_ratio;
1262        let new_a = new_position + (new_position - c) / cacp_ratio;
1263
1264        // Sort ctrl1 and ctrl2 such ctrl1 is the least affected (or most constrained).
1265        let mut ctrl1 = self.ctrl1;
1266        let mut ctrl2 = self.ctrl2;
1267        if t < S::HALF {
1268            core::mem::swap(&mut ctrl1, &mut ctrl2);
1269        }
1270
1271        // Move the most constrained control point by a subset of the uniform displacement
1272        // depending on the weight.
1273        let uniform_displacement = new_a - old_a;
1274        let f = if t < S::HALF {
1275            S::TWO * weight
1276        } else {
1277            S::TWO * (S::ONE - weight)
1278        };
1279        let mut new_ctrl1 = ctrl1 + uniform_displacement * f;
1280
1281        // Now that the most constrained control point is placed there is only one position
1282        // for the least constrained control point that satisfies the constant ratios.
1283        let d1_pre = (old_a - ctrl1).length();
1284        let d12_pre = (self.ctrl2 - self.ctrl1).length();
1285
1286        let mut new_ctrl2 = new_ctrl1 + (new_a - new_ctrl1) * (d12_pre / d1_pre);
1287
1288        if t < S::HALF {
1289            core::mem::swap(&mut new_ctrl1, &mut new_ctrl2);
1290        }
1291
1292        CubicBezierSegment {
1293            from: self.from,
1294            ctrl1: new_ctrl1,
1295            ctrl2: new_ctrl2,
1296            to: self.to,
1297        }
1298    }
1299
1300    pub fn to_f32(&self) -> CubicBezierSegment<f32> {
1301        CubicBezierSegment {
1302            from: self.from.to_f32(),
1303            ctrl1: self.ctrl1.to_f32(),
1304            ctrl2: self.ctrl2.to_f32(),
1305            to: self.to.to_f32(),
1306        }
1307    }
1308
1309    pub fn to_f64(&self) -> CubicBezierSegment<f64> {
1310        CubicBezierSegment {
1311            from: self.from.to_f64(),
1312            ctrl1: self.ctrl1.to_f64(),
1313            ctrl2: self.ctrl2.to_f64(),
1314            to: self.to.to_f64(),
1315        }
1316    }
1317
1318    pub fn polynomial_form(&self) -> CubicBezierPolynomial<S> {
1319        CubicBezierPolynomial {
1320            a0: self.from.to_vector(),
1321            a1: (self.ctrl1 - self.from) * S::THREE,
1322            a2: self.from * S::THREE - self.ctrl1 * S::SIX + self.ctrl2.to_vector() * S::THREE,
1323            a3: self.to - self.from + (self.ctrl1 - self.ctrl2) * S::THREE
1324        }
1325    }
1326}
1327
1328impl<S: Scalar> Segment for CubicBezierSegment<S> {
1329    impl_segment!(S);
1330
1331    fn for_each_flattened_with_t(
1332        &self,
1333        tolerance: Self::Scalar,
1334        callback: &mut dyn FnMut(&LineSegment<S>, Range<S>),
1335    ) {
1336        self.for_each_flattened_with_t(tolerance, &mut |s, t| callback(s, t));
1337    }
1338}
1339
1340/// The polynomial form of a cubic bézier segment.
1341///
1342/// The `sample` implementation uses Horner's method and tends to be faster than
1343/// `CubicBezierSegment::sample`.
1344pub struct CubicBezierPolynomial<S> {
1345    pub a0: Vector<S>,
1346    pub a1: Vector<S>,
1347    pub a2: Vector<S>,
1348    pub a3: Vector<S>,
1349}
1350
1351impl<S: Scalar> CubicBezierPolynomial<S> {
1352    #[inline(always)]
1353    pub fn sample(&self, t: S) -> Point<S> {
1354        // Horner's method.
1355        let mut v = self.a0;
1356        let mut t2 = t;
1357        v += self.a1 * t2;
1358        t2 *= t;
1359        v += self.a2 * t2;
1360        t2 *= t;
1361        v += self.a3 * t2;
1362
1363        v.to_point()
1364    }
1365}
1366
1367/// Computes the number of line segments required to build a flattened approximation
1368/// of the curve with segments placed at regular `t` intervals.
1369fn flattened_segments_wang<S: Scalar>(curve: &CubicBezierSegment<S>, tolerance: S) -> S {
1370    let from = curve.from.to_vector();
1371    let ctrl1 = curve.ctrl1.to_vector();
1372    let ctrl2 = curve.ctrl2.to_vector();
1373    let to = curve.to.to_vector();
1374    let v1 = (from - ctrl1 * S::TWO + ctrl2) * S::SIX;
1375    let v2 = (ctrl1 - ctrl2 * S::TWO + to) * S::SIX;
1376    let l = v1.dot(v1).max(v2.dot(v2));
1377    let d = S::ONE / (S::EIGHT * tolerance);
1378    let err4 = l * d * d;
1379    let err4_f32 = err4.to_f32().unwrap_or(1.0);
1380
1381    // Avoid two square roots using a lookup table that contains
1382    // i^4 for  i in 1..25.
1383    const N: usize = 24;
1384    const LUT: [f32; N] = [
1385        1.0, 16.0, 81.0, 256.0, 625.0, 1296.0, 2401.0, 4096.0, 6561.0,
1386        10000.0, 14641.0, 20736.0, 28561.0, 38416.0, 50625.0, 65536.0,
1387        83521.0, 104976.0, 130321.0, 160000.0, 194481.0, 234256.0,
1388        279841.0, 331776.0
1389    ];
1390
1391    // If the value we are looking for is within the LUT, take the fast path
1392    if err4_f32 <= 331776.0 {
1393        #[allow(clippy::needless_range_loop)]
1394        for i in 0..N {
1395            if err4_f32 <= LUT[i] {
1396                return S::from(i + 1).unwrap_or(S::ONE);
1397            }
1398        }
1399    }
1400
1401    // Otherwise fall back to computing via two square roots.
1402    err4.sqrt().sqrt().max(S::ONE)
1403}
1404
1405pub struct Flattened<S: Scalar> {
1406    curve: CubicBezierPolynomial<S>,
1407    to: Point<S>,
1408    segments: u32,
1409    step: S,
1410    t: S,
1411}
1412
1413impl<S: Scalar> Flattened<S> {
1414    pub(crate) fn new(curve: &CubicBezierSegment<S>, tolerance: S) -> Self {
1415        debug_assert!(tolerance >= S::EPSILON * S::EPSILON);
1416
1417        let n = flattened_segments_wang(curve, tolerance);
1418
1419        let step = S::ONE / n;
1420
1421        Flattened {
1422            curve: curve.polynomial_form(),
1423            to: curve.to,
1424            segments: n.to_u32().unwrap_or(1),
1425            step,
1426            t: step
1427        }
1428    }
1429}
1430
1431impl<S: Scalar> Iterator for Flattened<S> {
1432    type Item = Point<S>;
1433
1434    fn next(&mut self) -> Option<Point<S>> {
1435        if self.segments == 0 {
1436            return None;
1437        }
1438
1439        self.segments -= 1;
1440        if self.segments == 0 {
1441            return Some(self.to)
1442        }
1443
1444        let p = self.curve.sample(self.t);
1445        self.t += self.step;
1446
1447        Some(p)
1448    }
1449
1450    fn size_hint(&self) -> (usize, Option<usize>) {
1451        let n = self.segments as usize;
1452        (n, Some(n))
1453    }
1454}
1455
1456#[cfg(test)]
1457fn print_arrays(a: &[Point<f32>], b: &[Point<f32>]) {
1458    std::println!("left:  {a:?}");
1459    std::println!("right: {b:?}");
1460}
1461
1462#[cfg(test)]
1463fn assert_approx_eq(a: &[Point<f32>], b: &[Point<f32>]) {
1464    if a.len() != b.len() {
1465        print_arrays(a, b);
1466        panic!("Lengths differ ({} != {})", a.len(), b.len());
1467    }
1468    for i in 0..a.len() {
1469        let threshold = 0.029;
1470        let dx = f32::abs(a[i].x - b[i].x);
1471        let dy = f32::abs(a[i].y - b[i].y);
1472        if dx > threshold || dy > threshold {
1473            print_arrays(a, b);
1474            std::println!("diff = {dx:?} {dy:?}");
1475            panic!("The arrays are not equal");
1476        }
1477    }
1478}
1479
1480#[test]
1481fn test_iterator_builder_1() {
1482    let tolerance = 0.01;
1483    let c1 = CubicBezierSegment {
1484        from: Point::new(0.0, 0.0),
1485        ctrl1: Point::new(1.0, 0.0),
1486        ctrl2: Point::new(1.0, 1.0),
1487        to: Point::new(0.0, 1.0),
1488    };
1489    let iter_points: Vec<Point<f32>> = c1.flattened(tolerance).collect();
1490    let mut builder_points = Vec::new();
1491    c1.for_each_flattened(tolerance, &mut |s| {
1492        builder_points.push(s.to);
1493    });
1494
1495    assert!(iter_points.len() > 2);
1496    assert_approx_eq(&iter_points[..], &builder_points[..]);
1497}
1498
1499#[test]
1500fn test_iterator_builder_2() {
1501    let tolerance = 0.01;
1502    let c1 = CubicBezierSegment {
1503        from: Point::new(0.0, 0.0),
1504        ctrl1: Point::new(1.0, 0.0),
1505        ctrl2: Point::new(0.0, 1.0),
1506        to: Point::new(1.0, 1.0),
1507    };
1508    let iter_points: Vec<Point<f32>> = c1.flattened(tolerance).collect();
1509    let mut builder_points = Vec::new();
1510    c1.for_each_flattened(tolerance, &mut |s| {
1511        builder_points.push(s.to);
1512    });
1513
1514    assert!(iter_points.len() > 2);
1515    assert_approx_eq(&iter_points[..], &builder_points[..]);
1516}
1517
1518#[test]
1519fn test_iterator_builder_3() {
1520    let tolerance = 0.01;
1521    let c1 = CubicBezierSegment {
1522        from: Point::new(141.0, 135.0),
1523        ctrl1: Point::new(141.0, 130.0),
1524        ctrl2: Point::new(140.0, 130.0),
1525        to: Point::new(131.0, 130.0),
1526    };
1527    let iter_points: Vec<Point<f32>> = c1.flattened(tolerance).collect();
1528    let mut builder_points = Vec::new();
1529    c1.for_each_flattened(tolerance, &mut |s| {
1530        builder_points.push(s.to);
1531    });
1532
1533    assert!(iter_points.len() > 2);
1534    assert_approx_eq(&iter_points[..], &builder_points[..]);
1535}
1536
1537#[test]
1538fn test_issue_19() {
1539    let tolerance = 0.15;
1540    let c1 = CubicBezierSegment {
1541        from: Point::new(11.71726, 9.07143),
1542        ctrl1: Point::new(1.889879, 13.22917),
1543        ctrl2: Point::new(18.142855, 19.27679),
1544        to: Point::new(18.142855, 19.27679),
1545    };
1546    let iter_points: Vec<Point<f32>> = c1.flattened(tolerance).collect();
1547    let mut builder_points = Vec::new();
1548    c1.for_each_flattened(tolerance, &mut |s| {
1549        builder_points.push(s.to);
1550    });
1551
1552    assert_approx_eq(&iter_points[..], &builder_points[..]);
1553
1554    assert!(iter_points.len() > 1);
1555}
1556
1557#[test]
1558fn test_issue_194() {
1559    let segment = CubicBezierSegment {
1560        from: Point::new(0.0, 0.0),
1561        ctrl1: Point::new(0.0, 0.0),
1562        ctrl2: Point::new(50.0, 70.0),
1563        to: Point::new(100.0, 100.0),
1564    };
1565
1566    let mut points = Vec::new();
1567    segment.for_each_flattened(0.1, &mut |s| {
1568        points.push(s.to);
1569    });
1570
1571    assert!(points.len() > 2);
1572}
1573
1574#[test]
1575fn flatten_with_t() {
1576    let segment = CubicBezierSegment {
1577        from: Point::new(0.0f32, 0.0),
1578        ctrl1: Point::new(0.0, 0.0),
1579        ctrl2: Point::new(50.0, 70.0),
1580        to: Point::new(100.0, 100.0),
1581    };
1582
1583    for tolerance in &[0.1, 0.01, 0.001, 0.0001] {
1584        let tolerance = *tolerance;
1585
1586        let mut a = Vec::new();
1587        segment.for_each_flattened(tolerance, &mut |s| {
1588            a.push(*s);
1589        });
1590
1591        let mut b = Vec::new();
1592        let mut ts = Vec::new();
1593        segment.for_each_flattened_with_t(tolerance, &mut |s, t| {
1594            b.push(*s);
1595            ts.push(t);
1596        });
1597
1598        assert_eq!(a, b);
1599
1600        for i in 0..b.len() {
1601            let sampled = segment.sample(ts[i].start);
1602            let point = b[i].from;
1603            let dist = (sampled - point).length();
1604            assert!(dist <= tolerance);
1605
1606            let sampled = segment.sample(ts[i].end);
1607            let point = b[i].to;
1608            let dist = (sampled - point).length();
1609            assert!(dist <= tolerance);
1610        }
1611    }
1612}
1613
1614#[test]
1615fn test_flatten_end() {
1616    let segment = CubicBezierSegment {
1617        from: Point::new(0.0, 0.0),
1618        ctrl1: Point::new(100.0, 0.0),
1619        ctrl2: Point::new(100.0, 100.0),
1620        to: Point::new(100.0, 200.0),
1621    };
1622
1623    let mut last = segment.from;
1624    segment.for_each_flattened(0.0001, &mut |s| {
1625        last = s.to;
1626    });
1627
1628    assert_eq!(last, segment.to);
1629}
1630
1631#[test]
1632fn test_flatten_point() {
1633    let segment = CubicBezierSegment {
1634        from: Point::new(0.0, 0.0),
1635        ctrl1: Point::new(0.0, 0.0),
1636        ctrl2: Point::new(0.0, 0.0),
1637        to: Point::new(0.0, 0.0),
1638    };
1639
1640    let mut last = segment.from;
1641    segment.for_each_flattened(0.0001, &mut |s| {
1642        last = s.to;
1643    });
1644
1645    assert_eq!(last, segment.to);
1646}
1647
1648#[test]
1649fn issue_652() {
1650    use crate::point;
1651
1652    let curve = CubicBezierSegment {
1653        from: point(-1061.0, -3327.0),
1654        ctrl1: point(-1061.0, -3177.0),
1655        ctrl2: point(-1061.0, -3477.0),
1656        to: point(-1061.0, -3327.0),
1657    };
1658
1659    for _ in curve.flattened(1.0) {}
1660    for _ in curve.flattened(0.1) {}
1661    for _ in curve.flattened(0.01) {}
1662
1663    curve.for_each_flattened(1.0, &mut |_| {});
1664    curve.for_each_flattened(0.1, &mut |_| {});
1665    curve.for_each_flattened(0.01, &mut |_| {});
1666}
1667
1668#[test]
1669fn fast_bounding_box_for_cubic_bezier_segment() {
1670    let a = CubicBezierSegment {
1671        from: Point::new(0.0, 0.0),
1672        ctrl1: Point::new(0.5, 1.0),
1673        ctrl2: Point::new(1.5, -1.0),
1674        to: Point::new(2.0, 0.0),
1675    };
1676
1677    let expected_aabb = Box2D {
1678        min: point(0.0, -1.0),
1679        max: point(2.0, 1.0),
1680    };
1681
1682    let actual_aabb = a.fast_bounding_box();
1683
1684    assert_eq!(expected_aabb, actual_aabb)
1685}
1686
1687#[test]
1688fn minimum_bounding_box_for_cubic_bezier_segment() {
1689    let a = CubicBezierSegment {
1690        from: Point::new(0.0, 0.0),
1691        ctrl1: Point::new(0.5, 2.0),
1692        ctrl2: Point::new(1.5, -2.0),
1693        to: Point::new(2.0, 0.0),
1694    };
1695
1696    let expected_bigger_aabb: Box2D<f32> = Box2D {
1697        min: point(0.0, -0.6),
1698        max: point(2.0, 0.6),
1699    };
1700    let expected_smaller_aabb: Box2D<f32> = Box2D {
1701        min: point(0.1, -0.5),
1702        max: point(2.0, 0.5),
1703    };
1704
1705    let actual_minimum_aabb = a.bounding_box();
1706
1707    assert!(expected_bigger_aabb.contains_box(&actual_minimum_aabb));
1708    assert!(actual_minimum_aabb.contains_box(&expected_smaller_aabb));
1709}
1710
1711#[test]
1712fn y_maximum_t_for_simple_cubic_segment() {
1713    let a = CubicBezierSegment {
1714        from: Point::new(0.0, 0.0),
1715        ctrl1: Point::new(0.5, 1.0),
1716        ctrl2: Point::new(1.5, 1.0),
1717        to: Point::new(2.0, 2.0),
1718    };
1719
1720    let expected_y_maximum = 1.0;
1721
1722    let actual_y_maximum = a.y_maximum_t();
1723
1724    assert_eq!(expected_y_maximum, actual_y_maximum)
1725}
1726
1727#[test]
1728fn y_minimum_t_for_simple_cubic_segment() {
1729    let a = CubicBezierSegment {
1730        from: Point::new(0.0, 0.0),
1731        ctrl1: Point::new(0.5, 1.0),
1732        ctrl2: Point::new(1.5, 1.0),
1733        to: Point::new(2.0, 0.0),
1734    };
1735
1736    let expected_y_minimum = 0.0;
1737
1738    let actual_y_minimum = a.y_minimum_t();
1739
1740    assert_eq!(expected_y_minimum, actual_y_minimum)
1741}
1742
1743#[test]
1744fn y_extrema_for_simple_cubic_segment() {
1745    let a = CubicBezierSegment {
1746        from: Point::new(0.0, 0.0),
1747        ctrl1: Point::new(1.0, 2.0),
1748        ctrl2: Point::new(2.0, 2.0),
1749        to: Point::new(3.0, 0.0),
1750    };
1751
1752    let mut n: u32 = 0;
1753    a.for_each_local_y_extremum_t(&mut |t| {
1754        assert_eq!(t, 0.5);
1755        n += 1;
1756    });
1757    assert_eq!(n, 1);
1758}
1759
1760#[test]
1761fn x_extrema_for_simple_cubic_segment() {
1762    let a = CubicBezierSegment {
1763        from: Point::new(0.0, 0.0),
1764        ctrl1: Point::new(1.0, 2.0),
1765        ctrl2: Point::new(1.0, 2.0),
1766        to: Point::new(0.0, 0.0),
1767    };
1768
1769    let mut n: u32 = 0;
1770    a.for_each_local_x_extremum_t(&mut |t| {
1771        assert_eq!(t, 0.5);
1772        n += 1;
1773    });
1774    assert_eq!(n, 1);
1775}
1776
1777#[test]
1778fn x_maximum_t_for_simple_cubic_segment() {
1779    let a = CubicBezierSegment {
1780        from: Point::new(0.0, 0.0),
1781        ctrl1: Point::new(0.5, 1.0),
1782        ctrl2: Point::new(1.5, 1.0),
1783        to: Point::new(2.0, 0.0),
1784    };
1785    let expected_x_maximum = 1.0;
1786
1787    let actual_x_maximum = a.x_maximum_t();
1788
1789    assert_eq!(expected_x_maximum, actual_x_maximum)
1790}
1791
1792#[test]
1793fn x_minimum_t_for_simple_cubic_segment() {
1794    let a = CubicBezierSegment {
1795        from: Point::new(0.0, 0.0),
1796        ctrl1: Point::new(0.5, 1.0),
1797        ctrl2: Point::new(1.5, 1.0),
1798        to: Point::new(2.0, 0.0),
1799    };
1800
1801    let expected_x_minimum = 0.0;
1802
1803    let actual_x_minimum = a.x_minimum_t();
1804
1805    assert_eq!(expected_x_minimum, actual_x_minimum)
1806}
1807
1808#[test]
1809fn derivatives() {
1810    let c1 = CubicBezierSegment {
1811        from: Point::new(1.0, 1.0),
1812        ctrl1: Point::new(1.0, 2.0),
1813        ctrl2: Point::new(2.0, 1.0),
1814        to: Point::new(2.0, 2.0),
1815    };
1816
1817    assert_eq!(c1.dx(0.0), 0.0);
1818    assert_eq!(c1.dx(1.0), 0.0);
1819    assert_eq!(c1.dy(0.5), 0.0);
1820}
1821
1822#[test]
1823fn fat_line() {
1824    use crate::point;
1825
1826    let c1 = CubicBezierSegment {
1827        from: point(1.0f32, 2.0),
1828        ctrl1: point(1.0, 3.0),
1829        ctrl2: point(11.0, 11.0),
1830        to: point(11.0, 12.0),
1831    };
1832
1833    let (l1, l2) = c1.fat_line();
1834
1835    for i in 0..100 {
1836        let t = i as f32 / 99.0;
1837        assert!(l1.signed_distance_to_point(&c1.sample(t)) >= -0.000001);
1838        assert!(l2.signed_distance_to_point(&c1.sample(t)) <= 0.000001);
1839    }
1840
1841    let c2 = CubicBezierSegment {
1842        from: point(1.0f32, 2.0),
1843        ctrl1: point(1.0, 3.0),
1844        ctrl2: point(11.0, 14.0),
1845        to: point(11.0, 12.0),
1846    };
1847
1848    let (l1, l2) = c2.fat_line();
1849
1850    for i in 0..100 {
1851        let t = i as f32 / 99.0;
1852        assert!(l1.signed_distance_to_point(&c2.sample(t)) >= -0.000001);
1853        assert!(l2.signed_distance_to_point(&c2.sample(t)) <= 0.000001);
1854    }
1855
1856    let c3 = CubicBezierSegment {
1857        from: point(0.0f32, 1.0),
1858        ctrl1: point(0.5, 0.0),
1859        ctrl2: point(0.5, 0.0),
1860        to: point(1.0, 1.0),
1861    };
1862
1863    let (l1, l2) = c3.fat_line();
1864
1865    for i in 0..100 {
1866        let t = i as f32 / 99.0;
1867        assert!(l1.signed_distance_to_point(&c3.sample(t)) >= -0.000001);
1868        assert!(l2.signed_distance_to_point(&c3.sample(t)) <= 0.000001);
1869    }
1870}
1871
1872#[test]
1873fn is_linear() {
1874    let mut angle = 0.0;
1875    let center = Point::new(1000.0, -700.0);
1876    for _ in 0..100 {
1877        for i in 0..10 {
1878            for j in 0..10 {
1879                let (sin, cos) = f64::sin_cos(angle);
1880                let endpoint = Vector::new(cos * 100.0, sin * 100.0);
1881                let curve = CubicBezierSegment {
1882                    from: center - endpoint,
1883                    ctrl1: center + endpoint.lerp(-endpoint, i as f64 / 9.0),
1884                    ctrl2: center + endpoint.lerp(-endpoint, j as f64 / 9.0),
1885                    to: center + endpoint,
1886                };
1887                assert!(curve.is_linear(1e-10));
1888            }
1889        }
1890        angle += 0.001;
1891    }
1892}
1893
1894#[test]
1895fn test_monotonic() {
1896    use crate::point;
1897    let curve = CubicBezierSegment {
1898        from: point(1.0, 1.0),
1899        ctrl1: point(10.0, 2.0),
1900        ctrl2: point(1.0, 3.0),
1901        to: point(10.0, 4.0),
1902    };
1903
1904    curve.for_each_monotonic_range(&mut |range| {
1905        let sub_curve = curve.split_range(range);
1906        assert!(sub_curve.is_monotonic());
1907    });
1908}
1909
1910#[test]
1911fn test_line_segment_intersections() {
1912    use crate::point;
1913    fn assert_approx_eq(a: ArrayVec<(f32, f32), 3>, b: &[(f32, f32)], epsilon: f32) {
1914        for i in 0..a.len() {
1915            if f32::abs(a[i].0 - b[i].0) > epsilon || f32::abs(a[i].1 - b[i].1) > epsilon {
1916                std::println!("{a:?} != {b:?}");
1917            }
1918            assert!((a[i].0 - b[i].0).abs() <= epsilon && (a[i].1 - b[i].1).abs() <= epsilon);
1919        }
1920        assert_eq!(a.len(), b.len());
1921    }
1922
1923    let epsilon = 0.0001;
1924
1925    // Make sure we find intersections with horizontal and vertical lines.
1926
1927    let curve1 = CubicBezierSegment {
1928        from: point(-1.0, -1.0),
1929        ctrl1: point(0.0, 4.0),
1930        ctrl2: point(10.0, -4.0),
1931        to: point(11.0, 1.0),
1932    };
1933    let seg1 = LineSegment {
1934        from: point(0.0, 0.0),
1935        to: point(10.0, 0.0),
1936    };
1937    assert_approx_eq(
1938        curve1.line_segment_intersections_t(&seg1),
1939        &[(0.5, 0.5)],
1940        epsilon,
1941    );
1942
1943    let curve2 = CubicBezierSegment {
1944        from: point(-1.0, 0.0),
1945        ctrl1: point(0.0, 5.0),
1946        ctrl2: point(0.0, 5.0),
1947        to: point(1.0, 0.0),
1948    };
1949    let seg2 = LineSegment {
1950        from: point(0.0, 0.0),
1951        to: point(0.0, 5.0),
1952    };
1953    assert_approx_eq(
1954        curve2.line_segment_intersections_t(&seg2),
1955        &[(0.5, 0.75)],
1956        epsilon,
1957    );
1958}
1959
1960#[test]
1961fn test_parameters_for_value() {
1962    use crate::point;
1963    fn assert_approx_eq(a: ArrayVec<f32, 3>, b: &[f32], epsilon: f32) {
1964        for i in 0..a.len() {
1965            if f32::abs(a[i] - b[i]) > epsilon {
1966                std::println!("{a:?} != {b:?}");
1967            }
1968            assert!((a[i] - b[i]).abs() <= epsilon);
1969        }
1970        assert_eq!(a.len(), b.len());
1971    }
1972
1973    {
1974        let curve = CubicBezierSegment {
1975            from: point(0.0, 0.0),
1976            ctrl1: point(0.0, 8.0),
1977            ctrl2: point(10.0, 8.0),
1978            to: point(10.0, 0.0),
1979        };
1980
1981        let epsilon = 1e-4;
1982        assert_approx_eq(curve.solve_t_for_x(5.0), &[0.5], epsilon);
1983        assert_approx_eq(curve.solve_t_for_y(6.0), &[0.5], epsilon);
1984    }
1985    {
1986        let curve = CubicBezierSegment {
1987            from: point(0.0, 10.0),
1988            ctrl1: point(0.0, 10.0),
1989            ctrl2: point(10.0, 10.0),
1990            to: point(10.0, 10.0),
1991        };
1992
1993        assert_approx_eq(curve.solve_t_for_y(10.0), &[], 0.0);
1994    }
1995}
1996
1997#[test]
1998fn test_cubic_intersection_deduping() {
1999    use crate::point;
2000
2001    let epsilon = 0.0001;
2002
2003    // Two "line segments" with 3-fold overlaps, intersecting in their overlaps for a total of nine
2004    // parameter intersections.
2005    let line1 = CubicBezierSegment {
2006        from: point(-1_000_000.0, 0.0),
2007        ctrl1: point(2_000_000.0, 3_000_000.0),
2008        ctrl2: point(-2_000_000.0, -1_000_000.0),
2009        to: point(1_000_000.0, 2_000_000.0),
2010    };
2011    let line2 = CubicBezierSegment {
2012        from: point(-1_000_000.0, 2_000_000.0),
2013        ctrl1: point(2_000_000.0, -1_000_000.0),
2014        ctrl2: point(-2_000_000.0, 3_000_000.0),
2015        to: point(1_000_000.0, 0.0),
2016    };
2017    let intersections = line1.cubic_intersections(&line2);
2018    // (If you increase the coordinates above to 10s of millions, you get two returned intersection
2019    // points; i.e. the test fails.)
2020    assert_eq!(intersections.len(), 1);
2021    assert!(f64::abs(intersections[0].x) < epsilon);
2022    assert!(f64::abs(intersections[0].y - 1_000_000.0) < epsilon);
2023
2024    // Two self-intersecting curves that intersect in their self-intersections, for a total of four
2025    // parameter intersections.
2026    let curve1 = CubicBezierSegment {
2027        from: point(-10.0, -13.636363636363636),
2028        ctrl1: point(15.0, 11.363636363636363),
2029        ctrl2: point(-15.0, 11.363636363636363),
2030        to: point(10.0, -13.636363636363636),
2031    };
2032    let curve2 = CubicBezierSegment {
2033        from: point(13.636363636363636, -10.0),
2034        ctrl1: point(-11.363636363636363, 15.0),
2035        ctrl2: point(-11.363636363636363, -15.0),
2036        to: point(13.636363636363636, 10.0),
2037    };
2038    let intersections = curve1.cubic_intersections(&curve2);
2039    assert_eq!(intersections.len(), 1);
2040    assert!(f64::abs(intersections[0].x) < epsilon);
2041    assert!(f64::abs(intersections[0].y) < epsilon);
2042}
2043
2044#[test]
2045fn cubic_line_intersection_on_endpoint() {
2046    let l1 = LineSegment {
2047        from: Point::new(0.0, -100.0),
2048        to: Point::new(0.0, 100.0),
2049    };
2050
2051    let cubic = CubicBezierSegment {
2052        from: Point::new(0.0, 0.0),
2053        ctrl1: Point::new(20.0, 20.0),
2054        ctrl2: Point::new(20.0, 40.0),
2055        to: Point::new(0.0, 60.0),
2056    };
2057
2058    let intersections = cubic.line_segment_intersections_t(&l1);
2059
2060    assert_eq!(intersections.len(), 2);
2061    assert_eq!(intersections[0], (1.0, 0.8));
2062    assert_eq!(intersections[1], (0.0, 0.5));
2063
2064    let l2 = LineSegment {
2065        from: Point::new(0.0, 0.0),
2066        to: Point::new(0.0, 60.0),
2067    };
2068
2069    let intersections = cubic.line_segment_intersections_t(&l2);
2070
2071    assert!(intersections.is_empty());
2072
2073    let c1 = CubicBezierSegment {
2074        from: Point::new(0.0, 0.0),
2075        ctrl1: Point::new(20.0, 0.0),
2076        ctrl2: Point::new(20.0, 20.0),
2077        to: Point::new(0.0, 60.0),
2078    };
2079
2080    let c2 = CubicBezierSegment {
2081        from: Point::new(0.0, 60.0),
2082        ctrl1: Point::new(-40.0, 4.0),
2083        ctrl2: Point::new(-20.0, 20.0),
2084        to: Point::new(0.0, 00.0),
2085    };
2086
2087    let intersections = c1.cubic_intersections_t(&c2);
2088
2089    assert!(intersections.is_empty());
2090}
2091
2092#[test]
2093fn test_cubic_to_quadratics() {
2094    use euclid::approxeq::ApproxEq;
2095
2096    let quadratic = QuadraticBezierSegment {
2097        from: point(1.0, 2.0),
2098        ctrl: point(10.0, 5.0),
2099        to: point(0.0, 1.0),
2100    };
2101
2102    let mut count = 0;
2103    assert_eq!(quadratic.to_cubic().num_quadratics(0.0001), 1);
2104    quadratic
2105        .to_cubic()
2106        .for_each_quadratic_bezier(0.0001, &mut |c| {
2107            assert_eq!(count, 0);
2108            assert!(c.from.approx_eq(&quadratic.from));
2109            assert!(c.ctrl.approx_eq(&quadratic.ctrl));
2110            assert!(c.to.approx_eq(&quadratic.to));
2111            count += 1;
2112        });
2113
2114    let cubic = CubicBezierSegment {
2115        from: point(1.0f32, 1.0),
2116        ctrl1: point(10.0, 2.0),
2117        ctrl2: point(1.0, 3.0),
2118        to: point(10.0, 4.0),
2119    };
2120
2121    let mut prev = cubic.from;
2122    let mut count = 0;
2123    cubic.for_each_quadratic_bezier(0.01, &mut |c| {
2124        assert!(c.from.approx_eq(&prev));
2125        prev = c.to;
2126        count += 1;
2127    });
2128    assert!(prev.approx_eq(&cubic.to));
2129    assert!(count < 10);
2130    assert!(count > 4);
2131}