Implement simplify for bez paths

Implementation of SimplifyBezPath including ParamCurveFit.

This also fixes some issues with the trait and fit implementation.

Also a bit of cleanup, comment fixes, and micro-optimization.

examples/offset.rs

@@ -1,4 +1,4 @@
-use kurbo::{fit_to_bezpath, offset::CubicOffset, CubicBez, Shape};
+use kurbo::{offset::CubicOffset, CubicBez, Shape};
 
 fn main() {
     println!("<svg width='800' height='600' xmlns='http://www.w3.org/2000/svg'>");

examples/simplify.rs

@@ -0,0 +1,45 @@
+use kurbo::{BezPath, Point};
+
+fn plot_fn(f: &dyn Fn(f64) -> f64, d: &dyn Fn(f64) -> f64, xa: f64, xb: f64, n: usize) -> BezPath {
+    let width = 800.;
+    let dx = (xb - xa) / n as f64;
+    let xs = width / (xb - xa);
+    let ys = 250.;
+    let y_origin = 300.;
+    let plot = |x, y| Point::new((x - xa) * xs, y_origin - y * ys);
+    let mut x0 = xa;
+    let mut y0 = f(xa);
+    let mut d0 = d(xa);
+    let mut path = BezPath::new();
+    path.move_to(plot(x0, y0));
+    for i in 0..n {
+        let x3 = xa + dx * (i + 1) as f64;
+        let y3 = f(x3);
+        let d3 = d(x3);
+        let x1 = x0 + (1. / 3.) * dx;
+        let x2 = x3 - (1. / 3.) * dx;
+        let y1 = y0 + d0 * (1. / 3.) * dx;
+        let y2 = y3 - d3 * (1. / 3.) * dx;
+        path.curve_to(plot(x1, y1), plot(x2, y2), plot(x3, y3));
+        x0 = x3;
+        y0 = y3;
+        d0 = d3;
+    }
+    path
+}
+
+fn main() {
+    println!("<svg width='800' height='600' xmlns='http://www.w3.org/2000/svg'>");
+    let path = plot_fn(&|x| x.sin(), &|x| x.cos(), -8., 8., 20);
+    println!(
+        "  <path d='{}' stroke='#8c8' fill='none' stroke-width='5'/>",
+        path.to_svg()
+    );
+    let simpl = kurbo::simplify::SimplifyBezPath::new(path);
+    let simplified_path = kurbo::fit_to_bezpath_opt(&simpl, 0.1);
+    println!(
+        "  <path d='{}' stroke='#000' fill='none' stroke-width='1'/>",
+        simplified_path.to_svg()
+    );
+    println!("</svg>");
+}

src/fit.rs

@@ -52,19 +52,20 @@ pub trait ParamCurveFit {
     fn moment_integrals(&self, range: Range<f64>) -> (f64, f64, f64) {
         let t0 = 0.5 * (range.start + range.end);
         let dt = 0.5 * (range.end - range.start);
-        GAUSS_LEGENDRE_COEFFS_16
+        let (a, x, y) = GAUSS_LEGENDRE_COEFFS_16
             .iter()
             .map(|(wi, xi)| {
                 let t = t0 + xi * dt;
                 let (p, d) = self.sample_pt_deriv(t);
-                let a = (0.5 * wi) * d.x * p.y;
+                let a = wi * d.x * p.y;
                 let x = p.x * a;
                 let y = p.y * a;
                 (a, x, y)
             })
             .fold((0.0, 0.0, 0.0), |(a0, x0, y0), (a1, x1, y1)| {
                 (a0 + a1, x0 + x1, y0 + y1)
-            })
+            });
+        (a * dt, x * dt, y * dt)
     }
 
     /// Find a cusp or corner within the given range.
@@ -107,7 +108,7 @@ impl CurveFitSample {
         let c3 = p3.dot(self.tangent);
         solve_cubic(c0, c1, c2, c3)
             .into_iter()
-            .filter(|&t| t >= 0. && t <= 1.)
+            .filter(|t| (0.0..=1.0).contains(t))
             .collect()
     }
 }
@@ -171,17 +172,19 @@ pub fn fit_to_cubic(
     let th0 = start.tangent.atan2() - th;
     let th1 = th - end.tangent.atan2();
 
-    let (a, x, y) = source.moment_integrals(range.clone());
+    let (mut area, mut x, mut y) = source.moment_integrals(range.clone());
     let (x0, y0) = (start.p.x, start.p.y);
     let (dx, dy) = (d.x, d.y);
-    let dt = range.end - range.start;
-    let area = a * dt - dx * (y0 + 0.5 * dy);
-    // `area` is signed area of closed curve segment (ie with chord subtracted).
+    // Subtract off area of chord
+    area -= dx * (y0 + 0.5 * dy);
+    // `area` is signed area of closed curve segment.
     // This quantity is invariant to translation and rotation.
+
     // Subtract off moment of chord
-    let mut x = x * dt - dx * (x0 * y0 + 0.5 * (x0 * dy + y0 * dx) + (1. / 3.) * dy * dx);
-    let mut y = y * dt - dx * (y0 * y0 + y0 * dy + (1. / 3.) * dy * dy);
-    // Translate start point to origin; convert to moments.
+    let dy_3 = dy * (1. / 3.);
+    x -= dx * (x0 * y0 + 0.5 * (x0 * dy + y0 * dx) + dy_3 * dx);
+    y -= dx * (y0 * y0 + y0 * dy + dy_3 * dy);
+    // Translate start point to origin; convert raw integrals to moments.
     x -= x0 * area;
     y = 0.5 * y - y0 * area;
     // Rotate into place (this also scales up by chordlength for efficiency).
@@ -193,7 +196,7 @@ pub fn fit_to_cubic(
     let chord2_inv = chord2.recip();
     let unit_area = area * chord2_inv;
     let mx = moment * chord2_inv.powi(2);
-    // `area` is signed area scaled to unit chord; `moment` is scaled x moment
+    // `unit_area` is signed area scaled to unit chord; `mx` is scaled x moment
 
     let chord = chord2.sqrt();
     let aff = Affine::translate(start.p.to_vec2()) * Affine::rotate(th) * Affine::scale(chord);
@@ -379,17 +382,21 @@ pub fn fit_to_bezpath_opt(source: &impl ParamCurveFit, accuracy: f64) -> BezPath
     path
 }
 
+fn measure_one_seg(source: &impl ParamCurveFit, range: Range<f64>) -> Option<f64> {
+    fit_to_cubic(source, range, HUGE).map(|(_, err2)| err2.sqrt())
+}
+
 /// An Ok return is the t1 that hits the desired accuracy.
 /// An Err return is the error of t0..1.
 fn fit_opt_segment(source: &impl ParamCurveFit, accuracy: f64, t0: f64) -> Result<f64, f64> {
-    let (_, err2) = fit_to_cubic(source, t0..1.0, HUGE).unwrap();
-    let err = err2.sqrt();
+    let missing_err = accuracy * 2.0;
+    let err = measure_one_seg(source, t0..1.0).unwrap_or(missing_err);
     if err <= accuracy {
         return Err(err);
     }
     let f = |x| {
-        let (_, err2) = fit_to_cubic(source, t0..x, HUGE).unwrap();
-        err2.sqrt() - accuracy
+        let err = measure_one_seg(source, t0..x).unwrap_or(missing_err);
+        err - accuracy
     };
     const EPS: f64 = 1e-9;
     let k1 = 2.0 / (1.0 - t0);
@@ -405,6 +412,6 @@ fn fit_opt_err_delta(source: &impl ParamCurveFit, accuracy: f64, n: usize) -> f6
         // error is highly non-monotonic. We should probably harvest that solution.
         t0 = fit_opt_segment(source, accuracy, t0).unwrap();
     }
-    let (_, err2) = fit_to_cubic(source, t0..1.0, HUGE).unwrap();
-    accuracy - err2.sqrt()
+    let err = measure_one_seg(source, t0..1.0).unwrap_or(accuracy * 2.0);
+    accuracy - err
 }

src/simplify.rs

@@ -1,124 +1,22 @@
 //! Simplification of a Bézier path
 
-use crate::CubicBez;
+use std::ops::Range;
 
-/// Compute moment integrals.
-pub fn moment_integrals(c: CubicBez) -> (f64, f64, f64) {
-    let (x0, y0) = (c.p0.x, c.p0.y);
-    let (x1, y1) = (c.p1.x, c.p1.y);
-    let (x2, y2) = (c.p2.x, c.p2.y);
-    let (x3, y3) = (c.p3.x, c.p3.y);
-    let r0 = x0 * y3;
-    let r1 = x0 * y0;
-    let r2 = x0 * y1;
-    let r3 = 3. * y2;
-    let r4 = r3 * x0;
-    let r5 = x1 * y3;
-    let r6 = 3. * r5;
-    let r7 = x2 * y3;
-    let r8 = 45. * x1;
-    let r9 = 45. * x0;
-    let r10 = 15. * r0;
-    let r11 = 18. * x0;
-    let r12 = 12. * r0;
-    let r13 = 27. * y2;
-    let r14 = x0.powi(2);
-    let r15 = 105. * y1;
-    let r16 = 30. * y2;
-    let r17 = x1.powi(2);
-    let r18 = x2.powi(2);
-    let r19 = x3.powi(2);
-    let r20 = 45. * y2;
-    let r21 = y0.powi(2);
-    let r22 = y1.powi(2);
-    let r23 = y2.powi(2);
-    let r24 = y3.powi(2);
-
-    let a = -r0 - 10. * r1 - 6. * r2 - r3 * x1 - r4 - r6 - 6. * r7
-        + 6. * x1 * y0
-        + 3. * x2 * y0
-        + 3. * x2 * y1
-        + x3 * y0
-        + 3. * x3 * y1
-        + 6. * x3 * y2
-        + 10. * x3 * y3;
-    let x = -5. * r0 * x3
-        - r10 * x1
-        - r11 * x2 * y2
-        - r12 * x2
-        - r13 * r17
-        - r13 * x1 * x2
-        - r14 * r15
-        - r14 * r16
-        - 280. * r14 * y0
-        - 5. * r14 * y3
-        + 45. * r17 * y0
-        - 18. * r17 * y3
-        + 18. * r18 * y0
-        + 27. * r18 * y1
-        - 45. * r18 * y3
-        + 5. * r19 * y0
-        + 30. * r19 * y1
-        + 105. * r19 * y2
-        + 280. * r19 * y3
-        - r2 * r8
-        - r4 * x3
-        - 30. * r5 * x3
-        - r7 * r8
-        - 105. * r7 * x3
-        - r9 * x1 * y2
-        + 105. * x0 * x1 * y0
-        + 30. * x0 * x2 * y0
-        + 5. * x0 * x3 * y0
-        + 3. * x0 * x3 * y1
-        + 45. * x1 * x2 * y0
-        + 27. * x1 * x2 * y1
-        + 12. * x1 * x3 * y0
-        + 15. * x2 * x3 * y0
-        + 18. * x1 * x3 * y1
-        + 45. * x2 * x3 * y1
-        + 45. * x2 * x3 * y2;
-    let y = -5. * r0 * y0
-        - r1 * r15
-        - r1 * r16
-        - r10 * y2
-        - r11 * r23
-        - r12 * y1
-        - r13 * x1 * y1
-        - r2 * r20
-        - r20 * r5
-        - r20 * r7
-        - 140. * r21 * x0
-        + 105. * r21 * x1
-        + 30. * r21 * x2
-        + 5. * r21 * x3
-        - r22 * r9
-        + 27. * r22 * x2
-        + 18. * r22 * x3
-        - 27. * r23 * x1
-        + 45. * r23 * x3
-        - 5. * r24 * x0
-        - 30. * r24 * x1
-        - 105. * r24 * x2
-        + 140. * r24 * x3
-        - 18. * r5 * y1
-        - r6 * y0
-        + 45. * x1 * y0 * y1
-        + 45. * x2 * y0 * y1
-        + 18. * x2 * y0 * y2
-        + 3. * x2 * y0 * y3
-        + 27. * x2 * y1 * y2
-        + 15. * x3 * y0 * y1
-        + 12. * x3 * y0 * y2
-        + 5. * x3 * y0 * y3
-        + 45. * x3 * y1 * y2
-        + 30. * x3 * y1 * y3
-        + 105. * x3 * y2 * y3;
-    (a * (1. / 20.), x * (1. / 840.), y * (1. / 420.))
+use crate::{
+    CubicBez, CurveFitSample, ParamCurve, ParamCurveDeriv, ParamCurveFit, PathEl, Point, Vec2,
+};
+
+/// A Bézier path which has been prepared for simplification.
+pub struct SimplifyBezPath(Vec<SimplifyCubic>);
+
+struct SimplifyCubic {
+    c: CubicBez,
+    // The inclusive prefix sum of the moment integrals
+    moments: (f64, f64, f64),
 }
 
 /// Compute moment integrals.
-pub fn moment_integrals2(c: CubicBez) -> (f64, f64, f64) {
+pub fn moment_integrals(c: CubicBez) -> (f64, f64, f64) {
     let (x0, y0) = (c.p0.x, c.p0.y);
     let (x1, y1) = (c.p1.x - x0, c.p1.y - y0);
     let (x2, y2) = (c.p2.x - x0, c.p2.y - y0);
@@ -173,7 +71,100 @@ pub fn moment_integrals2(c: CubicBez) -> (f64, f64, f64) {
         - r8 * y2;
 
     let mx = x * (1. / 840.) + x0 * area + 0.5 * x3 * lift;
-    let my = y * (1. / 420.) + y0 * a * 0.1 + x0 * lift;
+    let my = y * (1. / 420.) + y0 * a * 0.1 + y0 * lift;
 
     (area, mx, my)
 }
+
+impl SimplifyBezPath {
+    /// Set up a new Bézier path for simplification.
+    ///
+    /// Currently this is not dealing with discontinuities at all, but it
+    /// could be extended to do so.
+    pub fn new(path: impl IntoIterator<Item = PathEl>) -> Self {
+        let (mut a, mut x, mut y) = (0.0, 0.0, 0.0);
+        let els = crate::segments(path)
+            .map(|seg| {
+                let c = seg.to_cubic();
+                let (ai, xi, yi) = moment_integrals(c);
+                a += ai;
+                x += xi;
+                y += yi;
+                SimplifyCubic {
+                    c,
+                    moments: (a, x, y),
+                }
+            })
+            .collect();
+        SimplifyBezPath(els)
+    }
+
+    /// Resolve a `t` value to a cubic.
+    ///
+    /// Also return the resulting `t` value for the selected cubic.
+    fn scale(&self, t: f64) -> (usize, f64) {
+        let t_scale = t * self.0.len() as f64;
+        let t_floor = t_scale.floor();
+        (t_floor as usize, t_scale - t_floor)
+    }
+
+    fn moment_integrals(&self, i: usize, range: Range<f64>) -> (f64, f64, f64) {
+        if range.end == range.start {
+            (0.0, 0.0, 0.0)
+        } else {
+            moment_integrals(self.0[i].c.subsegment(range))
+        }
+    }
+}
+
+impl ParamCurveFit for SimplifyBezPath {
+    fn sample_pt_deriv(&self, t: f64) -> (Point, Vec2) {
+        let (mut i, mut t0) = self.scale(t);
+        let n = self.0.len();
+        if i == n {
+            i -= 1;
+            t0 = 1.0;
+        }
+        let c = self.0[i].c;
+        (c.eval(t0), c.deriv().eval(t0).to_vec2() * n as f64)
+    }
+
+    fn sample_pt_tangent(&self, t: f64, _: f64) -> CurveFitSample {
+        let (mut i, mut t0) = self.scale(t);
+        if i == self.0.len() {
+            i -= 1;
+            t0 = 1.0;
+        }
+        let c = self.0[i].c;
+        let p = c.eval(t0);
+        let tangent = c.deriv().eval(t0).to_vec2();
+        CurveFitSample { p, tangent }
+    }
+
+    // We could use the default implementation, but provide our own, mostly
+    // because it is possible to efficiently provide an analytically accurate
+    // answer.
+    fn moment_integrals(&self, range: Range<f64>) -> (f64, f64, f64) {
+        let (i0, t0) = self.scale(range.start);
+        let (i1, t1) = self.scale(range.end);
+        if i0 == i1 {
+            self.moment_integrals(i0, t0..t1)
+        } else {
+            let (a0, x0, y0) = self.moment_integrals(i0, t0..1.0);
+            let (a1, x1, y1) = self.moment_integrals(i1, 0.0..t1);
+            let (mut a, mut x, mut y) = (a0 + a1, x0 + x1, y0 + y1);
+            if i1 > i0 + 1 {
+                let (a2, x2, y2) = self.0[i0].moments;
+                let (a3, x3, y3) = self.0[i1 - 1].moments;
+                a += a3 - a2;
+                x += x3 - x2;
+                y += y3 - y2;
+            }
+            (a, x, y)
+        }
+    }
+
+    fn break_cusp(&self, _: Range<f64>) -> Option<f64> {
+        None
+    }
+}