1 //! A trait for curves parametrized by a scalar.
2 
3 use std::ops::Range;
4 
5 use arrayvec::ArrayVec;
6 
7 use crate::{Point, Rect};
8 
9 /// A curve parametrized by a scalar.
10 ///
11 /// If the result is interpreted as a point, this represents a curve.
12 /// But the result can be interpreted as a vector as well.
13 pub trait ParamCurve: Sized {
14     /// Evaluate the curve at parameter `t`.
15     ///
16     /// Generally `t` is in the range [0..1].
eval(&self, t: f64) -> Point17     fn eval(&self, t: f64) -> Point;
18 
19     /// Get a subsegment of the curve for the given parameter range.
subsegment(&self, range: Range<f64>) -> Self20     fn subsegment(&self, range: Range<f64>) -> Self;
21 
22     /// Subdivide into (roughly) halves.
23     #[inline]
subdivide(&self) -> (Self, Self)24     fn subdivide(&self) -> (Self, Self) {
25         (self.subsegment(0.0..0.5), self.subsegment(0.5..1.0))
26     }
27 
28     /// The start point.
start(&self) -> Point29     fn start(&self) -> Point {
30         self.eval(0.0)
31     }
32 
33     /// The end point.
end(&self) -> Point34     fn end(&self) -> Point {
35         self.eval(1.0)
36     }
37 }
38 
39 // TODO: I might not want to have separate traits for all these.
40 
41 /// A differentiable parametrized curve.
42 pub trait ParamCurveDeriv {
43     /// The parametric curve obtained by taking the derivative of this one.
44     type DerivResult: ParamCurve;
45 
46     /// The derivative of the curve.
47     ///
48     /// Note that the type of the return value is somewhat inaccurate, as
49     /// the derivative of a curve (mapping of param to point) is a mapping
50     /// of param to vector. We choose to accept this rather than have a
51     /// more complex type scheme.
deriv(&self) -> Self::DerivResult52     fn deriv(&self) -> Self::DerivResult;
53 
54     /// Estimate arclength using Gaussian quadrature.
55     ///
56     /// The coefficients are assumed to cover the range (-1..1), which is
57     /// traditional.
58     #[inline]
gauss_arclen(&self, coeffs: &[(f64, f64)]) -> f6459     fn gauss_arclen(&self, coeffs: &[(f64, f64)]) -> f64 {
60         let d = self.deriv();
61         coeffs
62             .iter()
63             .map(|(wi, xi)| wi * d.eval(0.5 * (xi + 1.0)).to_vec2().hypot())
64             .sum::<f64>()
65             * 0.5
66     }
67 }
68 
69 /// A parametrized curve that can have its arc length measured.
70 pub trait ParamCurveArclen: ParamCurve {
71     /// The arc length of the curve.
72     ///
73     /// The result is accurate to the given accuracy (subject to
74     /// roundoff errors for ridiculously low values). Compute time
75     /// may vary with accuracy, if the curve needs to be subdivided.
arclen(&self, accuracy: f64) -> f6476     fn arclen(&self, accuracy: f64) -> f64;
77 
78     /// Solve for the parameter that has the given arclength from the start.
79     ///
80     /// This implementation is bisection, which is very robust but not
81     /// necessarily the fastest. It does measure increasingly short
82     /// segments, though, which should be good for subdivision algorithms.
inv_arclen(&self, arclen: f64, accuracy: f64) -> f6483     fn inv_arclen(&self, arclen: f64, accuracy: f64) -> f64 {
84         // invariant: the curve's arclen on [0..t_last] + remaining = arclen
85         let mut remaining = arclen;
86         let mut t_last = 0.0;
87         let mut t0 = 0.0;
88         let mut t1 = 1.0;
89         let n = (-accuracy.log2()).ceil();
90         let inner_accuracy = accuracy / n;
91         let n = n as usize;
92         for i in 0..n {
93             let tm = 0.5 * (t0 + t1);
94             let (range, dir) = if tm > t_last {
95                 (t_last..tm, 1.0)
96             } else {
97                 (tm..t_last, -1.0)
98             };
99             let range_size = range.end - range.start;
100             let arc = self.subsegment(range).arclen(inner_accuracy);
101             remaining -= arc * dir;
102             if i == n - 1 || (remaining).abs() < accuracy {
103                 // Allocate remaining arc evenly.
104                 return tm + range_size * remaining / arc;
105             }
106             if remaining > 0.0 {
107                 t0 = tm;
108             } else {
109                 t1 = tm;
110             }
111             t_last = tm;
112         }
113         unreachable!();
114     }
115 }
116 
117 /// A parametrized curve that can have its signed area measured.
118 pub trait ParamCurveArea {
119     /// Compute the signed area under the curve.
120     ///
121     /// For a closed path, the signed area of the path is the sum of signed
122     /// areas of the segments. This is a variant of the "shoelace formula."
123     /// See:
124     /// <https://github.com/Pomax/bezierinfo/issues/44> and
125     /// <http://ich.deanmcnamee.com/graphics/2016/03/30/CurveArea.html>
126     ///
127     /// This can be computed exactly for Béziers thanks to Green's theorem,
128     /// and also for simple curves such as circular arcs. For more exotic
129     /// curves, it's probably best to subdivide to cubics. We leave that
130     /// to the caller, which is why we don't give an accuracy param here.
signed_area(&self) -> f64131     fn signed_area(&self) -> f64;
132 }
133 
134 /// A parametrized curve that reports the nearest point.
135 pub trait ParamCurveNearest {
136     /// Find the point on the curve nearest the given point.
137     ///
138     /// Returns the parameter and the square of the distance.
nearest(&self, p: Point, accuracy: f64) -> (f64, f64)139     fn nearest(&self, p: Point, accuracy: f64) -> (f64, f64);
140 }
141 
142 /// A parametrized curve that reports its curvature.
143 pub trait ParamCurveCurvature: ParamCurveDeriv
144 where
145     Self::DerivResult: ParamCurveDeriv,
146 {
147     /// Compute the signed curvature at parameter `t`.
148     #[inline]
curvature(&self, t: f64) -> f64149     fn curvature(&self, t: f64) -> f64 {
150         let deriv = self.deriv();
151         let deriv2 = deriv.deriv();
152         let d = deriv.eval(t).to_vec2();
153         let d2 = deriv2.eval(t).to_vec2();
154         // TODO: What's the convention for sign? I think it should match signed
155         // area - a positive area curve should have positive curvature.
156         d2.cross(d) * d.hypot2().powf(-1.5)
157     }
158 }
159 
160 /// The maximum number of extrema that can be reported in the `ParamCurveExtrema` trait.
161 ///
162 /// This is 4 to support cubic Béziers. If other curves are used, they should be
163 /// subdivided to limit the number of extrema.
164 pub const MAX_EXTREMA: usize = 4;
165 
166 /// A parametrized curve that reports its extrema.
167 pub trait ParamCurveExtrema: ParamCurve {
168     /// Compute the extrema of the curve.
169     ///
170     /// Only extrema within the interior of the curve count.
171     /// At most four extrema can be reported, which is sufficient for
172     /// cubic Béziers.
173     ///
174     /// The extrema should be reported in increasing parameter order.
extrema(&self) -> ArrayVec<[f64; MAX_EXTREMA]>175     fn extrema(&self) -> ArrayVec<[f64; MAX_EXTREMA]>;
176 
177     /// Return parameter ranges, each of which is monotonic within the range.
extrema_ranges(&self) -> ArrayVec<[Range<f64>; MAX_EXTREMA + 1]>178     fn extrema_ranges(&self) -> ArrayVec<[Range<f64>; MAX_EXTREMA + 1]> {
179         let mut result = ArrayVec::new();
180         let mut t0 = 0.0;
181         for t in self.extrema() {
182             result.push(t0..t);
183             t0 = t;
184         }
185         result.push(t0..1.0);
186         result
187     }
188 
189     /// The smallest rectangle that encloses the curve in the range (0..1).
bounding_box(&self) -> Rect190     fn bounding_box(&self) -> Rect {
191         let mut bbox = Rect::from_points(self.start(), self.end());
192         for t in self.extrema() {
193             bbox = bbox.union_pt(self.eval(t))
194         }
195         bbox
196     }
197 }
198