diff --git a/examples/offset.rs b/examples/offset.rs index 6350b644..3d52022e 100644 --- a/examples/offset.rs +++ b/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!(""); diff --git a/examples/simplify.rs b/examples/simplify.rs new file mode 100644 index 00000000..15b5ab48 --- /dev/null +++ b/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!(""); + let path = plot_fn(&|x| x.sin(), &|x| x.cos(), -8., 8., 20); + println!( + " ", + path.to_svg() + ); + let simpl = kurbo::simplify::SimplifyBezPath::new(path); + let simplified_path = kurbo::fit_to_bezpath_opt(&simpl, 0.1); + println!( + " ", + simplified_path.to_svg() + ); + println!(""); +} diff --git a/src/fit.rs b/src/fit.rs index be31e795..f3643182 100644 --- a/src/fit.rs +++ b/src/fit.rs @@ -52,19 +52,20 @@ pub trait ParamCurveFit { fn moment_integrals(&self, range: Range) -> (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) -> Option { + 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 { - 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 } diff --git a/src/simplify.rs b/src/simplify.rs index 4f31a4bd..02650add 100644 --- a/src/simplify.rs +++ b/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); + +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) -> 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) { + 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) { + 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) -> Option { + None + } +}