Merge pull request #289 from linebender/fit_robust

Improve robustness of curve fitting

src/cubicbez.rs

@@ -44,6 +44,15 @@ struct ToQuads {
     n: usize,
 }
 
+/// Classification result for cusp detection.
+#[derive(Debug)]
+pub enum CuspType {
+    /// Cusp is a loop.
+    Loop,
+    /// Cusp has two closely spaced inflection points.
+    DoubleInflection,
+}
+
 impl CubicBez {
     /// Create a new cubic Bézier segment.
     #[inline]
@@ -348,6 +357,100 @@ impl CubicBez {
             .filter(|t| *t >= 0.0 && *t <= 1.0)
             .collect()
     }
+
+    /// Preprocess a cubic Bézier to ease numerical robustness.
+    ///
+    /// If the cubic Bézier segment has zero or near-zero derivatives, perturb
+    /// the control points to make it easier to process (especially offset and
+    /// stroke), avoiding numerical robustness problems.
+    pub(crate) fn regularize(&self, dimension: f64) -> CubicBez {
+        let mut c = *self;
+        // First step: if control point is too near the endpoint, nudge it away
+        // along the tangent.
+        let dim2 = dimension * dimension;
+        if c.p0.distance_squared(c.p1) < dim2 {
+            let d02 = c.p0.distance_squared(c.p2);
+            if d02 >= dim2 {
+                // TODO: moderate if this would move closer to p3
+                c.p1 = c.p0.lerp(c.p2, (dim2 / d02).sqrt());
+            } else {
+                c.p1 = c.p0.lerp(c.p3, 1.0 / 3.0);
+                c.p2 = c.p3.lerp(c.p0, 1.0 / 3.0);
+                return c;
+            }
+        }
+        if c.p3.distance_squared(c.p2) < dim2 {
+            let d13 = c.p1.distance_squared(c.p2);
+            if d13 >= dim2 {
+                // TODO: moderate if this would move closer to p0
+                c.p2 = c.p3.lerp(c.p1, (dim2 / d13).sqrt());
+            } else {
+                c.p1 = c.p0.lerp(c.p3, 1.0 / 3.0);
+                c.p2 = c.p3.lerp(c.p0, 1.0 / 3.0);
+                return c;
+            }
+        }
+        if let Some(cusp_type) = self.detect_cusp(dimension) {
+            let d01 = c.p1 - c.p0;
+            let d01h = d01.hypot();
+            let d23 = c.p3 - c.p2;
+            let d23h = d23.hypot();
+            match cusp_type {
+                CuspType::Loop => {
+                    c.p1 += (dimension / d01h) * d01;
+                    c.p2 -= (dimension / d23h) * d23;
+                }
+                CuspType::DoubleInflection => {
+                    // Avoid making control distance smaller than dimension
+                    if d01h > 2.0 * dimension {
+                        c.p1 -= (dimension / d01h) * d01;
+                    }
+                    if d23h > 2.0 * dimension {
+                        c.p2 += (dimension / d23h) * d23;
+                    }
+                }
+            }
+        }
+        c
+    }
+
+    /// Detect whether there is a cusp.
+    ///
+    /// Return a cusp classification if there is a cusp with curvature greater than
+    /// the reciprocal of the given dimension.
+    fn detect_cusp(&self, dimension: f64) -> Option<CuspType> {
+        let d01 = self.p1 - self.p0;
+        let d02 = self.p2 - self.p0;
+        let d03 = self.p3 - self.p0;
+        let d12 = self.p2 - self.p1;
+        let d23 = self.p3 - self.p2;
+        let det_012 = d01.cross(d02);
+        let det_123 = d12.cross(d23);
+        let det_013 = d01.cross(d03);
+        let det_023 = d02.cross(d03);
+        if det_012 * det_123 > 0.0 && det_012 * det_013 < 0.0 && det_012 * det_023 < 0.0 {
+            let q = self.deriv();
+            // accuracy isn't used for quadratic nearest
+            let nearest = q.nearest(Point::ORIGIN, 1e-9);
+            // detect whether curvature at minimum derivative exceeds 1/dimension,
+            // without division.
+            let d = q.eval(nearest.t);
+            let d2 = q.deriv().eval(nearest.t);
+            let cross = d.to_vec2().cross(d2.to_vec2());
+            if nearest.distance_sq.powi(3) < (cross * dimension).powi(2) {
+                let a = 3. * det_012 + det_023 - 2. * det_013;
+                let b = -3. * det_012 + det_013;
+                let c = det_012;
+                let d = b * b - 4. * a * c;
+                if d > 0.0 {
+                    return Some(CuspType::DoubleInflection);
+                } else {
+                    return Some(CuspType::Loop);
+                }
+            }
+        }
+        None
+    }
 }
 
 /// An iterator for cubic beziers.

src/fit.rs

@@ -15,7 +15,7 @@ use crate::{
         factor_quartic_inner, solve_cubic, solve_itp_fallible, solve_quadratic,
         GAUSS_LEGENDRE_COEFFS_16,
     },
-    Affine, BezPath, CubicBez, ParamCurve, ParamCurveArclen, Point, Vec2,
+    Affine, BezPath, CubicBez, Line, ParamCurve, ParamCurveArclen, ParamCurveNearest, Point, Vec2,
 };
 
 #[cfg(not(feature = "std"))]
@@ -172,6 +172,17 @@ fn fit_to_bezpath_rec(
 ) {
     let start = range.start;
     let end = range.end;
+    let start_p = source.sample_pt_tangent(range.start, 1.0).p;
+    let end_p = source.sample_pt_tangent(range.end, -1.0).p;
+    if start_p.distance_squared(end_p) <= accuracy * accuracy {
+        if let Some((c, _)) = try_fit_line(source, accuracy, range, start_p, end_p) {
+            if path.is_empty() {
+                path.move_to(c.p0);
+            }
+            path.curve_to(c.p1, c.p2, c.p3);
+            return;
+        }
+    }
     if let Some(t) = source.break_cusp(start..end) {
         fit_to_bezpath_rec(source, start..t, accuracy, path);
         fit_to_bezpath_rec(source, t..end, accuracy, path);
@@ -317,6 +328,42 @@ impl CurveDist {
 const D_PENALTY_ELBOW: f64 = 0.65;
 const D_PENALTY_SLOPE: f64 = 2.0;
 
+/// Try fitting a line.
+///
+/// This is especially useful for very short chords, in which the standard
+/// cubic fit is not numerically stable. The tangents are not considered, so
+/// it's useful in cusp and near-cusp situations where the tangents are not
+/// reliable, as well.
+///
+/// Returns the line raised to a cubic and the error, if within tolerance.
+fn try_fit_line(
+    source: &impl ParamCurveFit,
+    accuracy: f64,
+    range: Range<f64>,
+    start: Point,
+    end: Point,
+) -> Option<(CubicBez, f64)> {
+    let acc2 = accuracy * accuracy;
+    let chord_l = Line::new(start, end);
+    const SHORT_N: usize = 7;
+    let mut max_err2 = 0.0;
+    let dt = (range.end - range.start) / (SHORT_N + 1) as f64;
+    for i in 0..SHORT_N {
+        let t = range.start + (i + 1) as f64 * dt;
+        let p = source.sample_pt_deriv(t).0;
+        let err2 = chord_l.nearest(p, accuracy).distance_sq;
+        if err2 > acc2 {
+            // Not in tolerance; likely subdivision will help.
+            return None;
+        }
+        max_err2 = err2.max(max_err2);
+    }
+    let p1 = start.lerp(end, 1.0 / 3.0);
+    let p2 = end.lerp(start, 1.0 / 3.0);
+    let c = CubicBez::new(start, p1, p2, end);
+    Some((c, max_err2))
+}
+
 /// Fit a single cubic to a range of the source curve.
 ///
 /// Returns the cubic segment and the square of the error.
@@ -326,15 +373,16 @@ pub fn fit_to_cubic(
     range: Range<f64>,
     accuracy: f64,
 ) -> Option<(CubicBez, f64)> {
-    // TODO: handle case where chord is short; return `None` if subdivision
-    // will be useful (ie resulting chord is longer), otherwise simple short
-    // cubic (maybe raised line).
-
     let start = source.sample_pt_tangent(range.start, 1.0);
     let end = source.sample_pt_tangent(range.end, -1.0);
     let d = end.p - start.p;
-    let th = d.atan2();
     let chord2 = d.hypot2();
+    let acc2 = accuracy * accuracy;
+    if chord2 <= acc2 {
+        // Special case very short chords; try to fit a line.
+        return try_fit_line(source, accuracy, range, start.p, end.p);
+    }
+    let th = d.atan2();
     fn mod_2pi(th: f64) -> f64 {
         let th_scaled = th * core::f64::consts::FRAC_1_PI * 0.5;
         core::f64::consts::PI * 2.0 * (th_scaled - th_scaled.round())
@@ -371,7 +419,6 @@ pub fn fit_to_cubic(
     let chord = chord2.sqrt();
     let aff = Affine::translate(start.p.to_vec2()) * Affine::rotate(th) * Affine::scale(chord);
     let mut curve_dist = CurveDist::from_curve(source, range);
-    let acc2 = accuracy * accuracy;
     let mut best_c = None;
     let mut best_err2 = None;
     for (cand, d0, d1) in cubic_fit(th0, th1, unit_area, mx) {

src/offset.rs

@@ -65,6 +65,11 @@ pub struct CubicOffset {
 
 impl CubicOffset {
     /// Create a new curve from Bézier segment and offset.
+    ///
+    /// This method should only be used if the Bézier is smooth. Use
+    /// [`new_regularized`]instead to deal with a wider range of inputs.
+    ///
+    /// [`new_regularized`]: Self::new_regularized
     pub fn new(c: CubicBez, d: f64) -> Self {
         let q = c.deriv();
         let d0 = q.p0.to_vec2();
@@ -80,6 +85,14 @@ impl CubicOffset {
         }
     }
 
+    /// Create a new curve from Bézier segment and offset, with numerical robustness tweaks.
+    ///
+    /// The dimension represents a minimum feature size; the regularization is allowed to
+    /// perturb the curve by this amount in order to improve the robustness.
+    pub fn new_regularized(c: CubicBez, d: f64, dimension: f64) -> Self {
+        Self::new(c.regularize(dimension), d)
+    }
+
     fn eval_offset(&self, t: f64) -> Vec2 {
         let dp = self.q.eval(t).to_vec2();
         let norm = Vec2::new(-dp.y, dp.x);
@@ -158,3 +171,27 @@ impl ParamCurveFit for CubicOffset {
         Some(x)
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::CubicOffset;
+    use crate::{fit_to_bezpath, fit_to_bezpath_opt, CubicBez, PathEl};
+
+    // This test tries combinations of parameters that have caused problems in the past.
+    #[test]
+    fn pathological_curves() {
+        let curve = CubicBez {
+            p0: (-1236.3746269978635, 152.17981429574826).into(),
+            p1: (-1175.18662093517, 108.04721798590596).into(),
+            p2: (-1152.142883879584, 105.76260301083356).into(),
+            p3: (-1151.842639804639, 105.73040758939104).into(),
+        };
+        let offset = 3603.7267536453924;
+        let accuracy = 0.1;
+        let offset_path = CubicOffset::new(curve, offset);
+        let path = fit_to_bezpath_opt(&offset_path, accuracy);
+        assert!(matches!(path.iter().next(), Some(PathEl::MoveTo(_))));
+        let path_opt = fit_to_bezpath(&offset_path, accuracy);
+        assert!(matches!(path_opt.iter().next(), Some(PathEl::MoveTo(_))));
+    }
+}