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