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(_))));
+ }
+}