From 62fb3a7d024d66423329c885905379864a33c621 Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Wed, 7 Jun 2023 16:54:48 -0700 Subject: [PATCH 1/3] Improve robustness of curve fitting The main improvement is to try to fit a line when the chord is very short. Very short segments can occur when the cusp finding reports imprecise results, and also the main cubic fit logic is not robust when the chord is short. I believe the simple subdivision case is now fairly robust in that it won't keep recursing. This is in spite of not having an explicit bound on recursion depth; the theory is that each subdivision reduces the parameter space in half, and at some point you'll reach the point where the chord is shorter than the given tolerance. This is true even for an adversarial break_cusp, as it's bypassed in the short chord case. The logic for the optimized case is not as careful, in particular break_cusp is not bypassed. I'm curious whether there are still failure cases. This is likely not the end of the numerical robustness work. In particular, break_cusp can be better tuned to not over-report cusps, and the tangent detection can be improved for cusp and near-cusp source curves. But hopefully this commit gets us to where we at least return a valid fit. Also adds a test lightly adapted from #269. Progress toward #269 and #279 --- src/fit.rs | 61 +++++++++++++++++++++++++++++++++++++++++++++------ src/offset.rs | 24 ++++++++++++++++++++ 2 files changed, 78 insertions(+), 7 deletions(-) diff --git a/src/fit.rs b/src/fit.rs index ed7dda74..e090e6ce 100644 --- a/src/fit.rs +++ b/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, + 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, 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) { diff --git a/src/offset.rs b/src/offset.rs index 1e1df5b3..f8eff7e2 100644 --- a/src/offset.rs +++ b/src/offset.rs @@ -158,3 +158,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(_)))); + } +} From 6c4907d3b77059ccc752ee9253969e307ede567e Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Tue, 20 Jun 2023 12:51:48 -0700 Subject: [PATCH 2/3] Draft of regularization method This is one possible technique for improving the robustness of parallel curve generation - a method to detect cusps and "regularize" by perturbing the control points. The implementation handles cases where the control points are coincident with the endpoints, and also simple cusps, but hasn't been explicitly tuned to handle cases where the control points are nearly co-linear. It is likely that additional tuning of the "dimension" parameter would be worthwhile. As a first approximation, it can just be the accuracy parameter. --- src/cubicbez.rs | 103 ++++++++++++++++++++++++++++++++++++++++++++++++ src/offset.rs | 13 ++++++ 2 files changed, 116 insertions(+) diff --git a/src/cubicbez.rs b/src/cubicbez.rs index 68be82e1..286149fd 100644 --- a/src/cubicbez.rs +++ b/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 { + 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. diff --git a/src/offset.rs b/src/offset.rs index f8eff7e2..e6728142 100644 --- a/src/offset.rs +++ b/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); From 08f0905565b1ba208cc4a91dedfbe82d5edd244e Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Tue, 20 Jun 2023 14:15:57 -0700 Subject: [PATCH 3/3] Rustfmt --- src/cubicbez.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cubicbez.rs b/src/cubicbez.rs index 286149fd..2100db6c 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -443,9 +443,9 @@ impl CubicBez { let c = det_012; let d = b * b - 4. * a * c; if d > 0.0 { - return Some(CuspType::DoubleInflection) + return Some(CuspType::DoubleInflection); } else { - return Some(CuspType::Loop) + return Some(CuspType::Loop); } } }