From 071762be7a08f715bf92910ba534d527bc0d04dc Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Mon, 1 Feb 2021 09:25:59 -0800 Subject: [PATCH 01/12] Experiment with Euler spirals Checkpoint of experiment work --- Cargo.toml | 3 + examples/euler.rs | 75 ++++++++++++++++ src/euler.rs | 213 ++++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 5 ++ 4 files changed, 296 insertions(+) create mode 100644 examples/euler.rs create mode 100644 src/euler.rs diff --git a/Cargo.toml b/Cargo.toml index e881e536..30f210c5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -28,3 +28,6 @@ features = ["derive"] # This is used for research but not really needed; maybe refactor. [dev-dependencies] rand = "0.7.2" + +[features] +euler = [] diff --git a/examples/euler.rs b/examples/euler.rs new file mode 100644 index 00000000..4768c1b7 --- /dev/null +++ b/examples/euler.rs @@ -0,0 +1,75 @@ +use kurbo::{Affine, CubicBez, FitEulerResult, Point, Shape, fit_euler}; +use rand::random; + +fn check_euler_int_accuracy() { + for _ in 0..1000000 { + let k0= 8.0 * rand::random::(); + let k1 = 8.0 * rand::random::(); + let accurate = kurbo::integ_euler_12n(k0, k1, 8); + let est = kurbo::integ_euler_12n(k0, k1, 4); + let err = (est.0 - accurate.0).abs().max(est.1 - accurate.1).abs(); + let est_err_raw = 0.2 * k0 * k0 + k1.abs(); + if est_err_raw < 14.2 && err > 1e-9 { + println!("est_err_raw = {:?}, err = {:e}", est_err_raw, err); + } + } +} + +// errors for 1e-9 threshold +// n == 1: 1.01 +// n == 2: 3.96 +// n == 3: 8.23 +// n == 4: 14.2 + +// integ_euler_8 +// est = 0.112: 1e-9 +// 0.62: 1e-6 + +fn rand_cubic() -> CubicBez { + let p0 = Point::new(0., 0.); + let p1 = Point::new(rand::random::() * 0.5, rand::random::() * 2.0 - 1.0); + let p2 = Point::new(rand::random::() * 0.5 + 0.5, rand::random::() * 2.0 - 1.0); + let p3 = Point::new(1., 0.); + CubicBez::new(p0, p1, p2, p3) +} + +fn cubic_err_scatter() { + println!( + r##" + + + + "##); + for _ in 0..100000 { + let c = rand_cubic(); + let es = from_cubic(c); + let err = es.cubic_euler_err(c, 100); + //println!("{:?} {}", c, es.cubic_euler_err(c)); + const ERR_SCALE: f64 = 100.0; + let x = 950. * (err * ERR_SCALE).min(1.0); + let est_err = es.est_cubic_err(c); + let y = 750. * (est_err * ERR_SCALE).min(1.0); + if y < 750. { + let a = Affine::new([20., 0., 0., 20., x, y]); + let path = (a * c).into_path(1e-6); + println!(r##""##, + path.to_svg() + ); + } + } + println!( + r#" + +"# + ); +} + +fn main() { + cubic_err_scatter() +} + +fn from_cubic(c: CubicBez) -> FitEulerResult { + let th0 = (c.p1 - c.p0).atan2(); + let th1 = -(c.p3 - c.p2).atan2(); + fit_euler(th0, th1) +} diff --git a/src/euler.rs b/src/euler.rs new file mode 100644 index 00000000..f82bc65e --- /dev/null +++ b/src/euler.rs @@ -0,0 +1,213 @@ +use crate::{CubicBez, ParamCurve, ParamCurveArclen, Point}; + +fn integ_euler_12(k0: f64, k1: f64) -> (f64, f64) { + let t1_1 = k0; + let t1_2 = 0.5 * k1; + let t2_2 = t1_1 * t1_1; + let t2_3 = 2. * (t1_1 * t1_2); + let t2_4 = t1_2 * t1_2; + let t3_4 = t2_2 * t1_2 + t2_3 * t1_1; + let t3_6 = t2_4 * t1_2; + let t4_4 = t2_2 * t2_2; + let t4_5 = 2. * (t2_2 * t2_3); + let t4_6 = 2. * (t2_2 * t2_4) + t2_3 * t2_3; + let t4_7 = 2. * (t2_3 * t2_4); + let t4_8 = t2_4 * t2_4; + let t5_6 = t4_4 * t1_2 + t4_5 * t1_1; + let t5_8 = t4_6 * t1_2 + t4_7 * t1_1; + let t5_10 = t4_8 * t1_2; + let t6_6 = t4_4 * t2_2; + let t6_7 = t4_4 * t2_3 + t4_5 * t2_2; + let t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2; + let t6_9 = t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2; + let t6_10 = t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2; + let t7_8 = t6_6 * t1_2 + t6_7 * t1_1; + let t7_10 = t6_8 * t1_2 + t6_9 * t1_1; + let t8_8 = t6_6 * t2_2; + let t8_9 = t6_6 * t2_3 + t6_7 * t2_2; + let t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2; + let t9_10 = t8_8 * t1_2 + t8_9 * t1_1; + let t10_10 = t8_8 * t2_2; + let mut u = 1.; + u -= (1./24.) * t2_2 + (1./160.) * t2_4; + u += (1./1920.) * t4_4 + (1./10752.) * t4_6 + (1./55296.) * t4_8; + u -= (1./322560.) * t6_6 + (1./1658880.) * t6_8 + (1./8110080.) * t6_10; + u += (1./92897280.) * t8_8 + (1./454164480.) * t8_10; + u -= 2.4464949595157930e-11 * t10_10; + let mut v = (1./12.) * t1_2; + v -= (1./480.) * t3_4 + (1./2688.) * t3_6; + v += (1./53760.) * t5_6 + (1./276480.) * t5_8 + (1./1351680.) * t5_10; + v -= (1./11612160.) * t7_8 + (1./56770560.) * t7_10; + v += 2.4464949595157932e-10 * t9_10; + (u, v) +} + +/* +pub fn integ_euler_8(k0: f64, k1: f64) -> (f64, f64) { + let t1_1 = k0; + let t1_2 = 0.5 * k1; + let t2_2 = t1_1 * t1_1; + let t2_3 = 2. * (t1_1 * t1_2); + let t2_4 = t1_2 * t1_2; + let t3_4 = t2_2 * t1_2 + t2_3 * t1_1; + let t3_6 = t2_4 * t1_2; + let t4_4 = t2_2 * t2_2; + let t4_5 = 2. * (t2_2 * t2_3); + let t4_6 = 2. * (t2_2 * t2_4) + t2_3 * t2_3; + let t5_6 = t4_4 * t1_2 + t4_5 * t1_1; + let t6_6 = t4_4 * t2_2; + let mut u = 1.; + u -= (1./24.) * t2_2 + (1./160.) * t2_4; + u += (1./1920.) * t4_4 + (1./10752.) * t4_6; + u -= (1./322560.) * t6_6; + let mut v = (1./12.) * t1_2; + v -= (1./480.) * t3_4 + (1./2688.) * t3_6; + v += (1./53760.) * t5_6; + (u, v) +} +*/ + +#[doc(hidden)] +pub fn integ_euler_12n(mut k0: f64, mut k1: f64, n: usize) -> (f64, f64) { + let th1 = k0; + let th2 = 0.5 * k1; + let ds = (n as f64).recip(); + + k0 *= ds; + k1 *= ds; + + let mut x = 0.0; + let mut y = 0.0; + let mut s = 0.5 * ds - 0.5; + + for i in 0..n { + let km0 = k1 * s + k0; + let km1 = k1 * ds; + + let (u, v) = integ_euler_12(km0, km1); + + let th = (th2 * s + th1) * s; + let cth = th.cos(); + let sth = th.sin(); + + x += cth * u - sth * v; + y += cth * v + sth * u; + s += ds; + } + (x * ds, y * ds) +} + +fn integ_euler(k0: f64, k1: f64) -> (f64, f64) { + let est_err_raw = 0.2 * k0 * k0 + k1.abs(); + if est_err_raw < 1.01 { + integ_euler_12(k0, k1) + } else { + let n = if est_err_raw < 3.96 { + 2 + } else if est_err_raw < 8.23 { + 3 + } else if est_err_raw < 14.2 { + 4 + } else { + // Maybe determine the formula? + // Also, for these huge deflections, cephes-style computation + // of the fresnel integrals is probably the winning strategy. + 8 + }; + integ_euler_12n(k0, k1, n) + } +} + +/// Parameters for an Euler spiral segment. Does not include endpoint geometry. +pub struct FitEulerResult { + // These can be recovered from k0, k1, and chth, but it's annoying. + th0: f64, + th1: f64, + k0: f64, + k1: f64, + chord: f64, + chth: f64, +} + +/// Find the Euler spiral parameters for the given deflection. +pub fn fit_euler(th0: f64, th1: f64) -> FitEulerResult { + // Note: we could skip the solving for very small deflection + let mut k1_old = 0.0; + let dth = th1 - th0; + let k0 = th0 + th1; + let mut k1 = (6.0 - (1./70.) * dth * dth - 0.1 * k0 * k0) * dth; + let mut error_old = dth; + for _ in 0..10 { + let (u, v) = integ_euler(k0, k1); + let chth = v.atan2(u); + let error = dth - (0.25 * k1 - 2.0 * chth); + if error.abs() < 1e-9 { + let chord = u.hypot(v); + return FitEulerResult { th0, th1, k0, k1, chord, chth }; + } + let new_k1 = k1 + (k1_old - k1) * error / (error - error_old); + k1_old = k1; + error_old = error; + k1 = new_k1; + } + panic!("fit_euler diverged on {}, {}", th0, th1); +} + +impl FitEulerResult { + fn th(&self, t: f64) -> f64 { + let u = t - 0.5; + // Maybe sign of following is wrong; this is confusing. But it matches spiro code. + (0.5 * self.k1 * u + self.k0) * u - self.chth + } + + /// TODO + // Param t in [0, 1]. Return value assumes chord (0, 0) - (1, 0) + pub fn xy(&self, t: f64) -> Point { + let thm = self.th(t * 0.5); + let k0 = self.k0; + let k1 = self.k1; + let (u, v) = integ_euler((k0 + k1 * 0.5 * (t - 1.0)) * t, k1 * t * t); + let s = t / self.chord * thm.sin(); + let c = t / self.chord * thm.cos(); + let x = u * c - v * s; + let y = -v * c - u * s; + Point::new(x, y) + } + + /// Calculate error from cubic bezier. + /// + /// This is a fairly brute-force technique, finely sampling the bezier and + /// reporting RMS distance error. + pub fn cubic_euler_err(&self, cubic: CubicBez, n: usize) -> f64 { + // One way to improve this would be compute arclengths for each segment + // and cumulative sum them, rather than from 0 each time. + // + // We can also consider Legendre-Gauss quadrature. + // + // It's likely a rough approximation to arclength will be effective + // here, for example LGQ with a low order. + let cubic_len = cubic.arclen(1e-9); + let mut err = 0.0; + for i in 0..n { + let t = (i + 1) as f64 / ((n + 1) as f64); + let norm_len = cubic.subsegment(0.0..t).arclen(1e-9) / cubic_len; + let cubic_xy = cubic.eval(t); + let euler_xy = self.xy(norm_len); + err += (cubic_xy - euler_xy).hypot2(); + } + (err / (n as f64)).sqrt() + } + + /// Estimate the error wrt the cubic. + /// + /// It's an experiment to see how useful arclength is. It correlates fairly well, + /// but can both underestimate and overestimate, so it's better to just make + /// cubic_euler_err better. + pub fn est_cubic_err(&self, cubic: CubicBez) -> f64 { + let cubic_len = cubic.arclen(1e-9); + let err_arclen = cubic_len * self.chord - 1.0; + + let err = self.cubic_euler_err(cubic, 3); + err.max(err_arclen) + } +} diff --git a/src/lib.rs b/src/lib.rs index 82f24f08..83480f3d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -118,3 +118,8 @@ pub use crate::size::*; pub use crate::svg::*; pub use crate::translate_scale::*; pub use crate::vec2::*; + +#[cfg(feature = "euler")] +mod euler; +#[cfg(feature = "euler")] +pub use crate::euler::*; From ec51c788c810d3950c323ff30b653c24eb3f5839 Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Sun, 7 Feb 2021 15:57:50 -0800 Subject: [PATCH 02/12] Continue development of Euler WIP --- examples/euler.rs | 24 +++-- src/euler.rs | 233 +++++++++++++++++++++++++++++++++++++--------- 2 files changed, 203 insertions(+), 54 deletions(-) diff --git a/examples/euler.rs b/examples/euler.rs index 4768c1b7..fb1ce267 100644 --- a/examples/euler.rs +++ b/examples/euler.rs @@ -1,9 +1,9 @@ -use kurbo::{Affine, CubicBez, FitEulerResult, Point, Shape, fit_euler}; -use rand::random; +use kurbo::{Affine, CubicBez, FitEulerResult, Point, Shape}; +#[allow(unused)] fn check_euler_int_accuracy() { for _ in 0..1000000 { - let k0= 8.0 * rand::random::(); + let k0 = 8.0 * rand::random::(); let k1 = 8.0 * rand::random::(); let accurate = kurbo::integ_euler_12n(k0, k1, 8); let est = kurbo::integ_euler_12n(k0, k1, 4); @@ -27,8 +27,14 @@ fn check_euler_int_accuracy() { fn rand_cubic() -> CubicBez { let p0 = Point::new(0., 0.); - let p1 = Point::new(rand::random::() * 0.5, rand::random::() * 2.0 - 1.0); - let p2 = Point::new(rand::random::() * 0.5 + 0.5, rand::random::() * 2.0 - 1.0); + let p1 = Point::new( + rand::random::() * 0.5, + rand::random::() * 2.0 - 1.0, + ); + let p2 = Point::new( + rand::random::() * 0.5 + 0.5, + rand::random::() * 2.0 - 1.0, + ); let p3 = Point::new(1., 0.); CubicBez::new(p0, p1, p2, p3) } @@ -39,7 +45,8 @@ fn cubic_err_scatter() { - "##); + "## + ); for _ in 0..100000 { let c = rand_cubic(); let es = from_cubic(c); @@ -52,7 +59,8 @@ fn cubic_err_scatter() { if y < 750. { let a = Affine::new([20., 0., 0., 20., x, y]); let path = (a * c).into_path(1e-6); - println!(r##""##, + println!( + r##""##, path.to_svg() ); } @@ -71,5 +79,5 @@ fn main() { fn from_cubic(c: CubicBez) -> FitEulerResult { let th0 = (c.p1 - c.p0).atan2(); let th1 = -(c.p3 - c.p2).atan2(); - fit_euler(th0, th1) + FitEulerResult::fit_euler(th0, th1) } diff --git a/src/euler.rs b/src/euler.rs index f82bc65e..e690bf85 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -1,4 +1,36 @@ -use crate::{CubicBez, ParamCurve, ParamCurveArclen, Point}; +use crate::{CubicBez, Line, ParamCurve, ParamCurveArclen, Point}; + +/// An Euler spiral segment. +pub struct EulerSeg { + p0: Point, + p1: Point, + es: FitEulerResult, +} +/// Parameters for an Euler spiral segment. Does not include endpoint geometry. +/// +/// This is something of an internal detail for `EulerSeg` and might not make +/// it to the public interface. It's public here for experimentation. +pub struct FitEulerResult { + k0: f64, + k1: f64, + chord: f64, + chth: f64, +} + +/// A path consisting of piecewise Euler spiral segments. +pub struct EulerPath(Vec); + +/// An element of a piecewise Euler spiral path. +pub enum EulerPathEl { + /// Start a new subpath at the given point. + MoveTo(Point), + /// A line segment to the given point. + LineTo(Point), + /// An Euler spiral segment to the given point. + EulerTo(FitEulerResult, Point), + /// Close the subpath. + ClosePath, +} fn integ_euler_12(k0: f64, k1: f64) -> (f64, f64) { let t1_1 = k0; @@ -29,19 +61,22 @@ fn integ_euler_12(k0: f64, k1: f64) -> (f64, f64) { let t9_10 = t8_8 * t1_2 + t8_9 * t1_1; let t10_10 = t8_8 * t2_2; let mut u = 1.; - u -= (1./24.) * t2_2 + (1./160.) * t2_4; - u += (1./1920.) * t4_4 + (1./10752.) * t4_6 + (1./55296.) * t4_8; - u -= (1./322560.) * t6_6 + (1./1658880.) * t6_8 + (1./8110080.) * t6_10; - u += (1./92897280.) * t8_8 + (1./454164480.) * t8_10; + u -= (1. / 24.) * t2_2 + (1. / 160.) * t2_4; + u += (1. / 1920.) * t4_4 + (1. / 10752.) * t4_6 + (1. / 55296.) * t4_8; + u -= (1. / 322560.) * t6_6 + (1. / 1658880.) * t6_8 + (1. / 8110080.) * t6_10; + u += (1. / 92897280.) * t8_8 + (1. / 454164480.) * t8_10; u -= 2.4464949595157930e-11 * t10_10; - let mut v = (1./12.) * t1_2; - v -= (1./480.) * t3_4 + (1./2688.) * t3_6; - v += (1./53760.) * t5_6 + (1./276480.) * t5_8 + (1./1351680.) * t5_10; - v -= (1./11612160.) * t7_8 + (1./56770560.) * t7_10; + let mut v = (1. / 12.) * t1_2; + v -= (1. / 480.) * t3_4 + (1. / 2688.) * t3_6; + v += (1. / 53760.) * t5_6 + (1. / 276480.) * t5_8 + (1. / 1351680.) * t5_10; + v -= (1. / 11612160.) * t7_8 + (1. / 56770560.) * t7_10; v += 2.4464949595157932e-10 * t9_10; (u, v) } +/// A cheaper approximation to the Euler spiral integral. This might be useful +/// if we have lots of low-deflection segments, but it's not clear it'll be a +/// whole lot faster. /* pub fn integ_euler_8(k0: f64, k1: f64) -> (f64, f64) { let t1_1 = k0; @@ -68,19 +103,21 @@ pub fn integ_euler_8(k0: f64, k1: f64) -> (f64, f64) { */ #[doc(hidden)] +/// Computation of the Euler spiral integral using subdivision. pub fn integ_euler_12n(mut k0: f64, mut k1: f64, n: usize) -> (f64, f64) { let th1 = k0; let th2 = 0.5 * k1; - let ds = (n as f64).recip(); + let ds = (n as f64).recip(); k0 *= ds; k1 *= ds; let mut x = 0.0; let mut y = 0.0; - let mut s = 0.5 * ds - 0.5; + let s0 = 0.5 * ds - 0.5; for i in 0..n { + let s = s0 + ds * (i as f64); let km0 = k1 * s + k0; let km1 = k1 * ds; @@ -92,11 +129,19 @@ pub fn integ_euler_12n(mut k0: f64, mut k1: f64, n: usize) -> (f64, f64) { x += cth * u - sth * v; y += cth * v + sth * u; - s += ds; } (x * ds, y * ds) } +/// The Euler spiral integral. +/// +/// This approximates Equation 8.1 of the thesis. +/// +/// Discussion question: should this take an accuracy parameter? +/// This version basically hardcodes 1e-9. The cost of accuracy is +/// expected to be modest; the polynomial approximation converges +/// very rapidly and we expect small deflections, especially for +/// piecewise Euler spiral representations of other curves. fn integ_euler(k0: f64, k1: f64) -> (f64, f64) { let est_err_raw = 0.2 * k0 * k0 + k1.abs(); if est_err_raw < 1.01 { @@ -111,57 +156,73 @@ fn integ_euler(k0: f64, k1: f64) -> (f64, f64) { } else { // Maybe determine the formula? // Also, for these huge deflections, cephes-style computation - // of the fresnel integrals is probably the winning strategy. + // of the Fresnel integrals is probably the winning strategy. 8 }; integ_euler_12n(k0, k1, n) } } -/// Parameters for an Euler spiral segment. Does not include endpoint geometry. -pub struct FitEulerResult { - // These can be recovered from k0, k1, and chth, but it's annoying. - th0: f64, - th1: f64, - k0: f64, - k1: f64, - chord: f64, - chth: f64, -} +impl FitEulerResult { + /// Find the Euler spiral parameters for the given deflection. + /// + /// TODO: use research for direct solution. + /// + /// Discussion question: should this take an accuracy parameter? + /// This version basically hardcodes 1e-9. + pub fn fit_euler(th0: f64, th1: f64) -> FitEulerResult { + // Note: we could skip the solving for very small deflection + let mut k1_old = 0.0; + let dth = th1 - th0; + let k0 = th0 + th1; + let mut k1 = (6.0 - (1. / 70.) * dth * dth - 0.1 * k0 * k0) * dth; + let mut error_old = dth; + for _ in 0..10 { + let (u, v) = integ_euler(k0, k1); + let chth = v.atan2(u); + let error = dth - (0.25 * k1 - 2.0 * chth); + if error.abs() < 1e-9 { + let chord = u.hypot(v); + return FitEulerResult { + k0, + k1, + chord, + chth, + }; + } + let new_k1 = k1 + (k1_old - k1) * error / (error - error_old); + k1_old = k1; + error_old = error; + k1 = new_k1; + } + panic!("fit_euler diverged on {}, {}", th0, th1); + } -/// Find the Euler spiral parameters for the given deflection. -pub fn fit_euler(th0: f64, th1: f64) -> FitEulerResult { - // Note: we could skip the solving for very small deflection - let mut k1_old = 0.0; - let dth = th1 - th0; - let k0 = th0 + th1; - let mut k1 = (6.0 - (1./70.) * dth * dth - 0.1 * k0 * k0) * dth; - let mut error_old = dth; - for _ in 0..10 { + /// Create a `FitEulerResult` from k0 and k1 parameters. + pub fn from_k0_k1(k0: f64, k1: f64) -> FitEulerResult { let (u, v) = integ_euler(k0, k1); let chth = v.atan2(u); - let error = dth - (0.25 * k1 - 2.0 * chth); - if error.abs() < 1e-9 { - let chord = u.hypot(v); - return FitEulerResult { th0, th1, k0, k1, chord, chth }; + let chord = u.hypot(v); + FitEulerResult { + k0, + k1, + chord, + chth, } - let new_k1 = k1 + (k1_old - k1) * error / (error - error_old); - k1_old = k1; - error_old = error; - k1 = new_k1; } - panic!("fit_euler diverged on {}, {}", th0, th1); -} -impl FitEulerResult { + /// Determine tangent angle at the given parameter. + /// + /// The sign may be confusing, but it matches the spiro code. When `t = 0`, + /// the result is `-th0`, and when `t = 1`, the result is `th1`. fn th(&self, t: f64) -> f64 { let u = t - 0.5; - // Maybe sign of following is wrong; this is confusing. But it matches spiro code. (0.5 * self.k1 * u + self.k0) * u - self.chth } - /// TODO - // Param t in [0, 1]. Return value assumes chord (0, 0) - (1, 0) + /// Evaluate the curve at the given parameter. + /// + /// The parameter is in the range 0..1, and the result goes from (0, 0) to (1, 0). pub fn xy(&self, t: f64) -> Point { let thm = self.th(t * 0.5); let k0 = self.k0; @@ -203,6 +264,8 @@ impl FitEulerResult { /// It's an experiment to see how useful arclength is. It correlates fairly well, /// but can both underestimate and overestimate, so it's better to just make /// cubic_euler_err better. + /// + /// This is an experiment and will probably be removed. pub fn est_cubic_err(&self, cubic: CubicBez) -> f64 { let cubic_len = cubic.arclen(1e-9); let err_arclen = cubic_len * self.chord - 1.0; @@ -211,3 +274,81 @@ impl FitEulerResult { err.max(err_arclen) } } + +impl EulerSeg { + /// Create a new Euler segment. + /// + /// TODO: document the conventions. An SVG would be especially nice. + pub fn new(p0: Point, p1: Point, th0: f64, th1: f64) -> EulerSeg { + let es = FitEulerResult::fit_euler(th0, th1); + EulerSeg { p0, p1, es } + } + + /// Report whether the segment is a straight line. + pub fn is_line(&self) -> bool { + self.es.k0 == 0.0 && self.es.k1 == 0.0 + } +} + +impl ParamCurve for EulerSeg { + fn eval(&self, t: f64) -> Point { + let Point { x, y } = self.es.xy(t); + let chord = self.p1 - self.p0; + Point::new( + self.p0.x + chord.x * x - chord.y * y, + self.p0.y + chord.x * y + chord.y * x, + ) + } + + fn subsegment(&self, range: std::ops::Range) -> Self { + let p0 = self.eval(range.start); + let p1 = self.eval(range.end); + let dt = range.end - range.start; + let k0 = dt * (self.es.k0 + 0.5 * (range.start + range.end - 1.0) * self.es.k1); + let k1 = dt * dt * self.es.k1; + let es = FitEulerResult::from_k0_k1(k0, k1); + EulerSeg { p0, p1, es } + } + + fn start(&self) -> Point { + self.p0 + } + + fn end(&self) -> Point { + self.p1 + } +} + +impl ParamCurveArclen for EulerSeg { + /// The arc length of the curve. + /// + /// Note that this implementation is fast and accurate. + fn arclen(&self, _accuracy: f64) -> f64 { + (self.p1 - self.p0).hypot() / self.es.chord + } + + /// The parameter that results in the given arc length. + /// + /// This implementation is also fast and accurate. + fn inv_arclen(&self, arclen: f64, _accuracy: f64) -> f64 { + arclen * self.es.chord / (self.p1 - self.p0).hypot() + } +} + +// TODO: other ParamCurve traits. The derivative is pretty straightforward but will +// need a struct. + +impl From for EulerSeg { + fn from(l: Line) -> EulerSeg { + EulerSeg { + p0: l.p0, + p1: l.p1, + es: FitEulerResult { + k0: 0., + k1: 0., + chord: 1., + chth: 0., + }, + } + } +} From 1c31c748c8a5861a3401d71d4d94165c446a8637 Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Tue, 9 Feb 2021 11:12:00 -0800 Subject: [PATCH 03/12] Add ParamCurve traits for Euler Adds derivatives and curvature. --- src/euler.rs | 167 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 166 insertions(+), 1 deletion(-) diff --git a/src/euler.rs b/src/euler.rs index e690bf85..18739ce4 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -1,15 +1,35 @@ -use crate::{CubicBez, Line, ParamCurve, ParamCurveArclen, Point}; +use crate::{ + CubicBez, Line, ParamCurve, ParamCurveArclen, ParamCurveCurvature, ParamCurveDeriv, Point, Vec2, +}; /// An Euler spiral segment. +#[derive(Clone, Copy, Debug)] pub struct EulerSeg { p0: Point, p1: Point, es: FitEulerResult, } + +/// The derivative of an Euler spiral segment. +#[derive(Clone, Copy)] +pub struct EulerSegDeriv { + c0: f64, + c1: f64, + c2: f64, + scale: f64, +} + +/// The second derivative of an Euler spiral segment. +pub struct EulerSegDeriv2(EulerSegDeriv); + /// Parameters for an Euler spiral segment. Does not include endpoint geometry. /// /// This is something of an internal detail for `EulerSeg` and might not make /// it to the public interface. It's public here for experimentation. +/// +/// It's entirely possible the disposition of this is to be inlined into `EulerSeg`. +/// I'm not sure it's useful by itself, or isn +#[derive(Clone, Copy, Debug)] pub struct FitEulerResult { k0: f64, k1: f64, @@ -32,6 +52,16 @@ pub enum EulerPathEl { ClosePath, } +/// An iterator producing euler segments from a cubic bezier. +pub struct CubicToEulerIter { + c: CubicBez, + tolerance: f64, + // [t0 * dt .. (t0 + 1) * dt] is the range we're + // currently considering. + t0: u64, + dt: f64, +} + fn integ_euler_12(k0: f64, k1: f64) -> (f64, f64) { let t1_1 = k0; let t1_2 = 0.5 * k1; @@ -284,6 +314,14 @@ impl EulerSeg { EulerSeg { p0, p1, es } } + /// Create an Euler segment from a cubic Bézier. + /// + /// The curve is fit according to G1 geometric Hermite interpolation, in + /// other words the endpoints and tangents match the given curve. + pub fn from_cubic(c: CubicBez) -> EulerSeg { + todo!() + } + /// Report whether the segment is a straight line. pub fn is_line(&self) -> bool { self.es.k0 == 0.0 && self.es.k1 == 0.0 @@ -335,6 +373,65 @@ impl ParamCurveArclen for EulerSeg { } } +impl ParamCurveDeriv for EulerSeg { + type DerivResult = EulerSegDeriv; + + fn deriv(&self) -> Self::DerivResult { + let FitEulerResult { k0, k1, chth, .. } = self.es; + EulerSegDeriv { + c0: 0.5 * k0 - 0.125 * k1 + chth + (self.p1 - self.p0).atan2(), + c1: -k0 + 0.5 * k1, + c2: -0.5 * k1, + scale: self.arclen(0.0), + } + } +} + +impl ParamCurveCurvature for EulerSeg { + fn curvature(&self, t: f64) -> f64 { + (self.es.k0 + (t - 0.5) * self.es.k1) * self.es.chord / (self.p1 - self.p0).hypot() + } +} + +impl ParamCurve for EulerSegDeriv { + fn eval(&self, t: f64) -> Point { + let theta = self.c0 + t * self.c1 + t * t * self.c2; + (self.scale * Vec2::from_angle(theta)).to_point() + } + + fn subsegment(&self, range: std::ops::Range) -> Self { + let t0 = range.start; + let t1 = range.end; + let dt = t1 - t0; + EulerSegDeriv { + c0: self.c0 + t0 * self.c1 + t0 * t0 * self.c2, + c1: dt * (self.c1 + t0 * self.c2), + c2: dt * dt * self.c2, + scale: dt * self.scale, + } + } +} + +impl ParamCurveDeriv for EulerSegDeriv { + type DerivResult = EulerSegDeriv2; + + fn deriv(&self) -> Self::DerivResult { + EulerSegDeriv2(self.clone()) + } +} + +impl ParamCurve for EulerSegDeriv2 { + fn eval(&self, t: f64) -> Point { + let p = self.0.eval(t); + let scale = self.0.c1 + 2.0 * t * self.0.c2; + Point::new(-p.y * scale, p.x * scale) + } + + fn subsegment(&self, range: std::ops::Range) -> Self { + EulerSegDeriv2(self.0.subsegment(range)) + } +} + // TODO: other ParamCurve traits. The derivative is pretty straightforward but will // need a struct. @@ -352,3 +449,71 @@ impl From for EulerSeg { } } } + +impl Iterator for CubicToEulerIter { + type Item = EulerPathEl; + + fn next(&mut self) -> Option { + let t0 = (self.t0 as f64) * self.dt; + if t0 == 1.0 { + return None; + } + loop { + let t1 = t0 + self.dt; + let cubic = self.c.subsegment(t0..t1); + let es: EulerSeg = todo!("from cubic"); + let err: f64 = todo!("compute error"); + if err <= self.tolerance { + self.t0 += 1; + let shift = self.t0.trailing_zeros(); + self.t0 >>= shift; + self.dt *= (1 << shift) as f64; + return Some(EulerPathEl::EulerTo(es.es, es.p1)) + } + self.dt *= 0.5; + // TODO: should probably limit recursion here. + } + } +} + +#[cfg(test)] +mod tests { + use crate::{EulerSeg, ParamCurve, ParamCurveCurvature, ParamCurveDeriv, Point}; + + #[test] + fn euler_deriv() { + let es = EulerSeg::new(Point::ORIGIN, Point::new(1., 1.), 0.2, 0.3); + let esd = es.deriv(); + const DELTA: f64 = 0.001; + for i in 0..11 { + let t = (i as f64) * 0.1; + let est_deriv = (es.eval(t + 0.5 * DELTA) - es.eval(t - 0.5 * DELTA)) * DELTA.recip(); + let computed_deriv = esd.eval(t); + assert!((est_deriv.to_point() - computed_deriv).hypot() < 1e-7); + } + } + + #[test] + fn euler_deriv_2() { + let es = EulerSeg::new(Point::ORIGIN, Point::new(1., 1.), 0.2, 0.3); + let esd = es.deriv(); + let esd2 = esd.deriv(); + const DELTA: f64 = 0.001; + for i in 0..11 { + let t = (i as f64) * 0.1; + let est_deriv = (esd.eval(t + 0.5 * DELTA) - esd.eval(t - 0.5 * DELTA)) * DELTA.recip(); + let computed_deriv = esd2.eval(t); + assert!((est_deriv.to_point() - computed_deriv).hypot() < 1e-7); + } + } + + #[test] + fn euler_curvature() { + let th = std::f64::consts::FRAC_PI_2; + let es = EulerSeg::new(Point::ORIGIN, Point::new(0., 1.), th, th); + for i in 0..11 { + let t = (i as f64) * 0.1; + assert!((es.curvature(t) - 2.0).abs() < 1e-9); + } + } +} From 38f791a22e402fff32620c7fa75ca65faba4e0c0 Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Thu, 11 Feb 2021 09:42:08 -0800 Subject: [PATCH 04/12] Start curve fitting Generate Euler spirals from cubics. --- examples/euler.rs | 14 ++--- src/euler.rs | 133 ++++++++++++++++++++++++++++++---------------- 2 files changed, 92 insertions(+), 55 deletions(-) diff --git a/examples/euler.rs b/examples/euler.rs index fb1ce267..bdfa11da 100644 --- a/examples/euler.rs +++ b/examples/euler.rs @@ -1,4 +1,4 @@ -use kurbo::{Affine, CubicBez, FitEulerResult, Point, Shape}; +use kurbo::{Affine, CubicBez, EulerSeg, Point, Shape}; #[allow(unused)] fn check_euler_int_accuracy() { @@ -49,12 +49,12 @@ fn cubic_err_scatter() { ); for _ in 0..100000 { let c = rand_cubic(); - let es = from_cubic(c); - let err = es.cubic_euler_err(c, 100); + let es = EulerSeg::from_cubic(c); + let err = es.cubic_euler_err(c, 20); //println!("{:?} {}", c, es.cubic_euler_err(c)); const ERR_SCALE: f64 = 100.0; let x = 950. * (err * ERR_SCALE).min(1.0); - let est_err = es.est_cubic_err(c); + let est_err = es.cubic_euler_err(c, 4); let y = 750. * (est_err * ERR_SCALE).min(1.0); if y < 750. { let a = Affine::new([20., 0., 0., 20., x, y]); @@ -75,9 +75,3 @@ fn cubic_err_scatter() { fn main() { cubic_err_scatter() } - -fn from_cubic(c: CubicBez) -> FitEulerResult { - let th0 = (c.p1 - c.p0).atan2(); - let th1 = -(c.p3 - c.p2).atan2(); - FitEulerResult::fit_euler(th0, th1) -} diff --git a/src/euler.rs b/src/euler.rs index 18739ce4..6bccd4e5 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -28,7 +28,7 @@ pub struct EulerSegDeriv2(EulerSegDeriv); /// it to the public interface. It's public here for experimentation. /// /// It's entirely possible the disposition of this is to be inlined into `EulerSeg`. -/// I'm not sure it's useful by itself, or isn +/// I'm not sure it's useful by itself. #[derive(Clone, Copy, Debug)] pub struct FitEulerResult { k0: f64, @@ -41,6 +41,7 @@ pub struct FitEulerResult { pub struct EulerPath(Vec); /// An element of a piecewise Euler spiral path. +#[derive(Clone, Copy, Debug)] pub enum EulerPathEl { /// Start a new subpath at the given point. MoveTo(Point), @@ -53,6 +54,8 @@ pub enum EulerPathEl { } /// An iterator producing euler segments from a cubic bezier. +/// +/// Discussion: should this be an anonymous (`from_fn`) type? pub struct CubicToEulerIter { c: CubicBez, tolerance: f64, @@ -253,7 +256,7 @@ impl FitEulerResult { /// Evaluate the curve at the given parameter. /// /// The parameter is in the range 0..1, and the result goes from (0, 0) to (1, 0). - pub fn xy(&self, t: f64) -> Point { + pub fn eval(&self, t: f64) -> Point { let thm = self.th(t * 0.5); let k0 = self.k0; let k1 = self.k1; @@ -264,11 +267,43 @@ impl FitEulerResult { let y = -v * c - u * s; Point::new(x, y) } +} + +impl EulerSeg { + /// Create a new Euler segment. + /// + /// TODO: document the conventions. An SVG would be especially nice. + pub fn new(p0: Point, p1: Point, th0: f64, th1: f64) -> EulerSeg { + let es = FitEulerResult::fit_euler(th0, th1); + EulerSeg { p0, p1, es } + } + + /// Create an Euler segment from a cubic Bézier. + /// + /// The curve is fit according to G1 geometric Hermite interpolation, in + /// other words the endpoints and tangents match the given curve. + pub fn from_cubic(c: CubicBez) -> EulerSeg { + let d01 = c.p1 - c.p0; + let d23 = c.p3 - c.p2; + let d03 = c.p3 - c.p0; + let th0 = d03.cross(d01).atan2(d03.dot(d01)); + let th1 = d23.cross(d03).atan2(d23.dot(d03)); + let es = FitEulerResult::fit_euler(th0, th1); + EulerSeg { + p0: c.p0, + p1: c.p3, + es, + } + } /// Calculate error from cubic bezier. /// - /// This is a fairly brute-force technique, finely sampling the bezier and + /// This is a fairly brute-force technique, sampling the bezier and /// reporting RMS distance error. + /// + /// Note: experimentation suggests that n = 4 is enough to estimate the + /// error fairly accurately. A future evolution may get rid of that + /// parameter, and possibly also a tolerance parameter. pub fn cubic_euler_err(&self, cubic: CubicBez, n: usize) -> f64 { // One way to improve this would be compute arclengths for each segment // and cumulative sum them, rather than from 0 each time. @@ -283,45 +318,12 @@ impl FitEulerResult { let t = (i + 1) as f64 / ((n + 1) as f64); let norm_len = cubic.subsegment(0.0..t).arclen(1e-9) / cubic_len; let cubic_xy = cubic.eval(t); - let euler_xy = self.xy(norm_len); + let euler_xy = self.eval(norm_len); err += (cubic_xy - euler_xy).hypot2(); } (err / (n as f64)).sqrt() } - /// Estimate the error wrt the cubic. - /// - /// It's an experiment to see how useful arclength is. It correlates fairly well, - /// but can both underestimate and overestimate, so it's better to just make - /// cubic_euler_err better. - /// - /// This is an experiment and will probably be removed. - pub fn est_cubic_err(&self, cubic: CubicBez) -> f64 { - let cubic_len = cubic.arclen(1e-9); - let err_arclen = cubic_len * self.chord - 1.0; - - let err = self.cubic_euler_err(cubic, 3); - err.max(err_arclen) - } -} - -impl EulerSeg { - /// Create a new Euler segment. - /// - /// TODO: document the conventions. An SVG would be especially nice. - pub fn new(p0: Point, p1: Point, th0: f64, th1: f64) -> EulerSeg { - let es = FitEulerResult::fit_euler(th0, th1); - EulerSeg { p0, p1, es } - } - - /// Create an Euler segment from a cubic Bézier. - /// - /// The curve is fit according to G1 geometric Hermite interpolation, in - /// other words the endpoints and tangents match the given curve. - pub fn from_cubic(c: CubicBez) -> EulerSeg { - todo!() - } - /// Report whether the segment is a straight line. pub fn is_line(&self) -> bool { self.es.k0 == 0.0 && self.es.k1 == 0.0 @@ -330,7 +332,7 @@ impl EulerSeg { impl ParamCurve for EulerSeg { fn eval(&self, t: f64) -> Point { - let Point { x, y } = self.es.xy(t); + let Point { x, y } = self.es.eval(t); let chord = self.p1 - self.p0; Point::new( self.p0.x + chord.x * x - chord.y * y, @@ -432,8 +434,7 @@ impl ParamCurve for EulerSegDeriv2 { } } -// TODO: other ParamCurve traits. The derivative is pretty straightforward but will -// need a struct. +// TODO: other ParamCurve traits. impl From for EulerSeg { fn from(l: Line) -> EulerSeg { @@ -450,8 +451,20 @@ impl From for EulerSeg { } } +impl CubicToEulerIter { + /// Subdivide a cubic Bézier into piecewise Euler spirals. + pub fn new(c: CubicBez, tolerance: f64) -> CubicToEulerIter { + CubicToEulerIter { + c, + tolerance, + t0: 0, + dt: 1.0, + } + } +} + impl Iterator for CubicToEulerIter { - type Item = EulerPathEl; + type Item = EulerSeg; fn next(&mut self) -> Option { let t0 = (self.t0 as f64) * self.dt; @@ -461,15 +474,16 @@ impl Iterator for CubicToEulerIter { loop { let t1 = t0 + self.dt; let cubic = self.c.subsegment(t0..t1); - let es: EulerSeg = todo!("from cubic"); - let err: f64 = todo!("compute error"); + let es = EulerSeg::from_cubic(cubic); + let err: f64 = es.cubic_euler_err(cubic, 4); if err <= self.tolerance { self.t0 += 1; let shift = self.t0.trailing_zeros(); self.t0 >>= shift; self.dt *= (1 << shift) as f64; - return Some(EulerPathEl::EulerTo(es.es, es.p1)) + return Some(es); } + self.t0 *= 2; self.dt *= 0.5; // TODO: should probably limit recursion here. } @@ -478,7 +492,10 @@ impl Iterator for CubicToEulerIter { #[cfg(test)] mod tests { - use crate::{EulerSeg, ParamCurve, ParamCurveCurvature, ParamCurveDeriv, Point}; + use crate::{ + CubicBez, CubicToEulerIter, EulerSeg, ParamCurve, ParamCurveCurvature, ParamCurveDeriv, + ParamCurveNearest, Point, + }; #[test] fn euler_deriv() { @@ -516,4 +533,30 @@ mod tests { assert!((es.curvature(t) - 2.0).abs() < 1e-9); } } + + #[test] + fn euler_from_cubic() { + let c = CubicBez::new((0., 0.), (1., 1.1), (2., 2.2), (3., 3.)); + let es = EulerSeg::from_cubic(c); + // Check endpoints match. These should actually be bit-identical. + assert!((es.p0 - c.p0).hypot() < 1e-15); + assert!((es.p1 - c.p3).hypot() < 1e-15); + // Check tangents match. This can have some error but should be small. + let des = es.deriv(); + let dc = c.deriv(); + assert!(dc.eval(0.0).to_vec2().atan2() - des.eval(0.0).to_vec2().atan2() < 1e-9); + assert!(dc.eval(1.0).to_vec2().atan2() - des.eval(1.0).to_vec2().atan2() < 1e-9); + } + + #[test] + fn cubic_to_euler() { + const TOLERANCE: f64 = 1e-4; + let c = CubicBez::new((0., 0.), (1., 1.1), (2., 2.2), (3., 3.)); + for es in CubicToEulerIter::new(c, TOLERANCE) { + for i in 0..11 { + let t = (i as f64) * 0.1; + assert!(c.nearest(es.eval(t), 1e-9).distance_sq < TOLERANCE.powi(2)); + } + } + } } From b62c929f2c2bf927576e8bc7d826f2f5697f71d0 Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Thu, 11 Feb 2021 13:47:18 -0800 Subject: [PATCH 05/12] Euler to cubic conversion Iterate from Euler spirals to cubics, with an error tolerance. --- src/euler.rs | 84 ++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 79 insertions(+), 5 deletions(-) diff --git a/src/euler.rs b/src/euler.rs index 6bccd4e5..8a1e6607 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -1,5 +1,6 @@ use crate::{ - CubicBez, Line, ParamCurve, ParamCurveArclen, ParamCurveCurvature, ParamCurveDeriv, Point, Vec2, + Affine, CubicBez, Line, ParamCurve, ParamCurveArclen, ParamCurveCurvature, ParamCurveDeriv, + PathEl, Point, Vec2, }; /// An Euler spiral segment. @@ -328,6 +329,65 @@ impl EulerSeg { pub fn is_line(&self) -> bool { self.es.k0 == 0.0 && self.es.k1 == 0.0 } + + /// Convert to cubic beziers. + pub fn to_cubics(&self, accuracy: f64) -> impl Iterator { + let this = *self; + let mut t0_int = 0usize; + let mut dt = 1.0; + let mut p0 = self.p0; + let chord_atan = (self.p1 - self.p0).atan2(); + let thresh = accuracy * self.es.chord / (self.p1 - self.p0).hypot(); + std::iter::from_fn(move || { + let t0 = (t0_int as f64) * dt; + if t0 == 1.0 { + return None; + } + loop { + let t1 = t0 + dt; + let k0 = dt * (this.es.k0 + 0.5 * (t0 + t1 - 1.0) * this.es.k1); + let k1 = dt * dt * this.es.k1; + let a0 = k0.abs(); + let a1 = k1.abs(); + // Error metric empirically determined. I can make a Python notebook available. + let err = 1.5e-5 * a0.powi(5) + + 6e-4 * a0 * a0 * a1 + + 1e-4 * a0 * a1 * a1 + + 3e-6 * a1.powi(3); + // TODO: scale error by arc length + if err * dt <= thresh { + let p1 = if t1 == 1.0 { this.p1 } else { this.eval(t1) }; + + let dp = p1 - p0; + // Transform to take (0, 0) - (1, 0) chord to p0 - p1. + let a = Affine::new([dp.x, dp.y, -dp.y, dp.x, p0.x, p0.y]); + + // Note: it's possible to this with rotation and normalization, + // avoiding the trig. + let d_atan = chord_atan - dp.atan2(); + let th0 = d_atan - this.es.th(t0); + let th1 = -d_atan + this.es.th(t1); + let v0 = Vec2::from_angle(th0); + let c0 = Point::new(0., 0.); + let c1 = c0 + 2. / 3. / (1. + v0.x) * v0; + let c3 = Point::new(1., 0.); + let v1 = Vec2::from_angle(-th1); + let c2 = c3 - 2. / 3. / (1. + v1.x) * v1; + + // Advance subdivision parameters + t0_int += 1; + let shift = t0_int.trailing_zeros(); + t0_int >>= shift; + dt *= (1 << shift) as f64; + p0 = p1; + + return Some(PathEl::CurveTo(a * c1, a * c2, p1)); + } + t0_int *= 2; + dt *= 0.5; + } + }) + } } impl ParamCurve for EulerSeg { @@ -492,10 +552,7 @@ impl Iterator for CubicToEulerIter { #[cfg(test)] mod tests { - use crate::{ - CubicBez, CubicToEulerIter, EulerSeg, ParamCurve, ParamCurveCurvature, ParamCurveDeriv, - ParamCurveNearest, Point, - }; + use crate::{CubicBez, CubicToEulerIter, EulerSeg, ParamCurve, ParamCurveArclen, ParamCurveCurvature, ParamCurveDeriv, ParamCurveNearest, PathEl, Point}; #[test] fn euler_deriv() { @@ -559,4 +616,21 @@ mod tests { } } } + + #[test] + fn euler_to_cubic() { + const TOLERANCE: f64 = 1e-4; + let es = EulerSeg::new(Point::new(0., 0.), Point::new(10., 0.), 0.5, -0.5); + let arclen = es.arclen(1e-9); + let mut approx_arclen = 0.0; + let mut p0 = es.p0; + for seg in es.to_cubics(TOLERANCE) { + if let PathEl::CurveTo(p1, p2, p3) = seg { + let seg = CubicBez::new(p0, p1, p2, p3); + approx_arclen += seg.arclen(1e-9); + p0 = p3; + } + } + assert!((arclen - approx_arclen).abs() < TOLERANCE); + } } From e1b18059ed21e2e4a127bf30c2d3f1fd7b81fc88 Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Sun, 14 Feb 2021 15:35:24 -0800 Subject: [PATCH 06/12] Refine Euler integration Provide accuracy (with verified error bounds). --- examples/euler.rs | 14 +++++++++ src/cubicbez.rs | 1 - src/euler.rs | 79 ++++++++++++++++++++++++++++++----------------- src/lib.rs | 3 +- src/point.rs | 1 - src/svg.rs | 1 - 6 files changed, 66 insertions(+), 33 deletions(-) diff --git a/examples/euler.rs b/examples/euler.rs index bdfa11da..e60ed29c 100644 --- a/examples/euler.rs +++ b/examples/euler.rs @@ -1,3 +1,17 @@ +// Copyright 2021 The kurbo Authors. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + use kurbo::{Affine, CubicBez, EulerSeg, Point, Shape}; #[allow(unused)] diff --git a/src/cubicbez.rs b/src/cubicbez.rs index 21937516..fc5234e2 100644 --- a/src/cubicbez.rs +++ b/src/cubicbez.rs @@ -431,7 +431,6 @@ mod tests { } #[test] - #[allow(clippy::float_cmp)] fn cubicbez_signed_area_linear() { // y = 1 - x let c = CubicBez::new( diff --git a/src/euler.rs b/src/euler.rs index 8a1e6607..4066de1c 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -1,3 +1,17 @@ +// Copyright 2021 The kurbo Authors. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + use crate::{ Affine, CubicBez, Line, ParamCurve, ParamCurveArclen, ParamCurveCurvature, ParamCurveDeriv, PathEl, Point, Vec2, @@ -167,32 +181,32 @@ pub fn integ_euler_12n(mut k0: f64, mut k1: f64, n: usize) -> (f64, f64) { (x * ds, y * ds) } -/// The Euler spiral integral. +/// Evaulate the Euler spiral integral. +/// +/// Compute the following integral to the desired accuracy. +/// +/// $$ +/// \int_{-0.5}^{0.5} \exp(i(k_0 s + 1/2 k_1 s^2)) ds +/// $$ /// -/// This approximates Equation 8.1 of the thesis. +/// This is discussed in section 8.1 of [Raph's thesis], and the error bounds +/// are validated in the notebook attached to the parallel curve blog post. /// -/// Discussion question: should this take an accuracy parameter? -/// This version basically hardcodes 1e-9. The cost of accuracy is -/// expected to be modest; the polynomial approximation converges -/// very rapidly and we expect small deflections, especially for -/// piecewise Euler spiral representations of other curves. -fn integ_euler(k0: f64, k1: f64) -> (f64, f64) { - let est_err_raw = 0.2 * k0 * k0 + k1.abs(); - if est_err_raw < 1.01 { +/// [Raph's thesis]: https://www.levien.com/phd/thesis.pdf +#[inline] +pub fn integ_euler(k0: f64, k1: f64, accuracy: f64) -> (f64, f64) { + let thresh = accuracy.powf(1.0 / 6.0); + integ_euler_thresh(k0, k1, thresh) +} + +fn integ_euler_thresh(k0: f64, k1: f64, thresh: f64) -> (f64, f64) { + let c1 = k1.abs(); + let c0 = k0.abs() + 0.5 * c1; + let est_err_raw = 0.006 * c0 * c0 + 0.029 * c1; + if est_err_raw < thresh { integ_euler_12(k0, k1) } else { - let n = if est_err_raw < 3.96 { - 2 - } else if est_err_raw < 8.23 { - 3 - } else if est_err_raw < 14.2 { - 4 - } else { - // Maybe determine the formula? - // Also, for these huge deflections, cephes-style computation - // of the Fresnel integrals is probably the winning strategy. - 8 - }; + let n = (est_err_raw / thresh).ceil() as usize; integ_euler_12n(k0, k1, n) } } @@ -212,7 +226,7 @@ impl FitEulerResult { let mut k1 = (6.0 - (1. / 70.) * dth * dth - 0.1 * k0 * k0) * dth; let mut error_old = dth; for _ in 0..10 { - let (u, v) = integ_euler(k0, k1); + let (u, v) = integ_euler(k0, k1, 1e-12); let chth = v.atan2(u); let error = dth - (0.25 * k1 - 2.0 * chth); if error.abs() < 1e-9 { @@ -234,7 +248,7 @@ impl FitEulerResult { /// Create a `FitEulerResult` from k0 and k1 parameters. pub fn from_k0_k1(k0: f64, k1: f64) -> FitEulerResult { - let (u, v) = integ_euler(k0, k1); + let (u, v) = integ_euler(k0, k1, 1e-12); let chth = v.atan2(u); let chord = u.hypot(v); FitEulerResult { @@ -257,11 +271,11 @@ impl FitEulerResult { /// Evaluate the curve at the given parameter. /// /// The parameter is in the range 0..1, and the result goes from (0, 0) to (1, 0). - pub fn eval(&self, t: f64) -> Point { + pub fn eval(&self, t: f64, accuracy: f64) -> Point { let thm = self.th(t * 0.5); let k0 = self.k0; let k1 = self.k1; - let (u, v) = integ_euler((k0 + k1 * 0.5 * (t - 1.0)) * t, k1 * t * t); + let (u, v) = integ_euler((k0 + k1 * 0.5 * (t - 1.0)) * t, k1 * t * t, accuracy); let s = t / self.chord * thm.sin(); let c = t / self.chord * thm.cos(); let x = u * c - v * s; @@ -392,7 +406,9 @@ impl EulerSeg { impl ParamCurve for EulerSeg { fn eval(&self, t: f64) -> Point { - let Point { x, y } = self.es.eval(t); + // The accuracy here is somewhat arbitrary, but should be adequate + // for most work, and not entail loss of efficiency. + let Point { x, y } = self.es.eval(t, 1e-9); let chord = self.p1 - self.p0; Point::new( self.p0.x + chord.x * x - chord.y * y, @@ -400,6 +416,8 @@ impl ParamCurve for EulerSeg { ) } + // TODO: I'm not sure this is right, possibly chord needs to be taken into + // account too? Need to test. fn subsegment(&self, range: std::ops::Range) -> Self { let p0 = self.eval(range.start); let p1 = self.eval(range.end); @@ -478,7 +496,7 @@ impl ParamCurveDeriv for EulerSegDeriv { type DerivResult = EulerSegDeriv2; fn deriv(&self) -> Self::DerivResult { - EulerSegDeriv2(self.clone()) + EulerSegDeriv2(*self) } } @@ -552,7 +570,10 @@ impl Iterator for CubicToEulerIter { #[cfg(test)] mod tests { - use crate::{CubicBez, CubicToEulerIter, EulerSeg, ParamCurve, ParamCurveArclen, ParamCurveCurvature, ParamCurveDeriv, ParamCurveNearest, PathEl, Point}; + use crate::{ + CubicBez, CubicToEulerIter, EulerSeg, ParamCurve, ParamCurveArclen, ParamCurveCurvature, + ParamCurveDeriv, ParamCurveNearest, PathEl, Point, + }; #[test] fn euler_deriv() { diff --git a/src/lib.rs b/src/lib.rs index 83480f3d..eca32a20 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -75,7 +75,8 @@ #![allow( clippy::unreadable_literal, clippy::many_single_char_names, - clippy::excessive_precision + clippy::excessive_precision, + clippy::float_cmp )] mod affine; diff --git a/src/point.rs b/src/point.rs index c4549cc5..ae4cf16d 100644 --- a/src/point.rs +++ b/src/point.rs @@ -284,7 +284,6 @@ mod tests { } #[test] - #[allow(clippy::float_cmp)] fn distance() { let p1 = Point::new(0., 10.); let p2 = Point::new(0., 5.); diff --git a/src/svg.rs b/src/svg.rs index 15b3205a..832f7cd3 100644 --- a/src/svg.rs +++ b/src/svg.rs @@ -503,7 +503,6 @@ mod tests { // Regression test for #51 #[test] - #[allow(clippy::float_cmp)] fn test_parse_svg_arc_pie() { let path = BezPath::from_svg("M 100 100 h 25 a 25 25 0 1 0 -25 25 z").unwrap(); // Approximate figures, but useful for regression testing From a5dc6f379be513bcb469e81ec46a05e3ae36091d Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Wed, 17 Feb 2021 08:13:11 -0800 Subject: [PATCH 07/12] Performance tuning on euler integral Add benchmarks and refine error estimation (the powf was potentially expensive). --- benches/euler.rs | 156 +++++++++++++++++++++++++++++++++++++++++++++++ src/euler.rs | 43 +++---------- 2 files changed, 163 insertions(+), 36 deletions(-) create mode 100644 benches/euler.rs diff --git a/benches/euler.rs b/benches/euler.rs new file mode 100644 index 00000000..99c707aa --- /dev/null +++ b/benches/euler.rs @@ -0,0 +1,156 @@ +// Copyright 2021 The kurbo Authors. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// Note: this `cfg` may not work, TODO either move to criterion or sort it out. +// I've been locally removing it to run benches. +#![cfg(nightly)] +#![feature(test)] +extern crate test; +use test::{Bencher, black_box}; + +use kurbo::*; + +// About 15ns +#[bench] +fn bench_euler_integ(b: &mut Bencher) { + b.iter(|| integ_euler(black_box(0.1), black_box(0.2), 1e-9)) +} + +// About 21ns. This is evidence that computing this root directly as part of the +// error estimation would be very expensive, unless it could be constant folded. +#[bench] +fn bench_sixth_root(b: &mut Bencher) { + b.iter(|| black_box(0.1f64).powf(1.0 / 6.0)) +} + +// A lower order approximation to the main integral. The code is here +// if we want to use it, and might be a speedup if the workload consists +// of very low deflection segments (and/or the error threshold is high), +// but for general use any loss of branch prediction coherence would +// probably be a lose. +fn integ_euler_8(k0: f64, k1: f64) -> (f64, f64) { + let t1_1 = k0; + let t1_2 = 0.5 * k1; + let t2_2 = t1_1 * t1_1; + let t2_3 = 2. * (t1_1 * t1_2); + let t2_4 = t1_2 * t1_2; + let t3_4 = t2_2 * t1_2 + t2_3 * t1_1; + let t3_6 = t2_4 * t1_2; + let t4_4 = t2_2 * t2_2; + let t4_5 = 2. * (t2_2 * t2_3); + let t4_6 = 2. * (t2_2 * t2_4) + t2_3 * t2_3; + let t5_6 = t4_4 * t1_2 + t4_5 * t1_1; + let t6_6 = t4_4 * t2_2; + let mut u = 1.; + u -= (1./24.) * t2_2 + (1./160.) * t2_4; + u += (1./1920.) * t4_4 + (1./10752.) * t4_6; + u -= (1./322560.) * t6_6; + let mut v = (1./12.) * t1_2; + v -= (1./480.) * t3_4 + (1./2688.) * t3_6; + v += (1./53760.) * t5_6; + (u, v) +} + +// About 4ns +#[bench] +fn bench_euler_integ_8(b: &mut Bencher) { + b.iter(|| integ_euler_8(black_box(0.1), black_box(0.2))) +} + +// About 240ns +#[bench] +fn bench_fit_euler(b: &mut Bencher) { + b.iter(|| FitEulerResult::fit_euler(black_box(0.1), black_box(0.2))) +} + +// Pretty much straight out of the notebook. Looking at the generated +// asm, the compiler does a good job +fn fast_fit_euler(th0: f64, th1: f64) -> f64 { + let k0 = th0 + th1; + let dth = th1 - th0; + let mut est = dth * 6.; + est += dth.powi(3) * (1. / -70.); + est += dth.powi(5) * (1. / -10780.); + est += dth.powi(7) * 2.769178184818219e-07; + est += dth * k0.powi(2) * (1. / -10.); + est += dth.powi(3) * k0.powi(2) * (1. / 4200.); + est += dth.powi(5) * k0.powi(2) * 1.6959677820260655e-05; + est += dth * k0.powi(4) * (1. / -1400.); + est += dth.powi(3) * k0.powi(4) * 6.84915970574303e-05; + est += dth * k0.powi(6) * -7.936475029053326e-06; + est +} + +// Same as above but using mostly even powers. A small improvement. +fn fast_fit_euler2(th0: f64, th1: f64) -> f64 { + let k0 = th0 + th1; + let dth = th1 - th0; + let mut est = 6.; + est += dth.powi(2) * (1. / -70.); + est += dth.powi(4) * (1. / -10780.); + est += dth.powi(6) * 2.769178184818219e-07; + est += k0.powi(2) * (1. / -10.); + est += dth.powi(2) * k0.powi(2) * (1. / 4200.); + est += dth.powi(4) * k0.powi(2) * 1.6959677820260655e-05; + est += k0.powi(4) * (1. / -1400.); + est += dth.powi(2) * k0.powi(4) * 6.84915970574303e-05; + est += k0.powi(6) * -7.936475029053326e-06; + est * dth +} + +// Compute both k1 and chord at the same time, which is what an +// actual implementation would use. +fn fast_fit_euler3(th0: f64, th1: f64) -> (f64, f64) { + let k0 = th0 + th1; + let dth = th1 - th0; + let mut k1 = 6.; + k1 += dth.powi(2) * (1. / -70.); + k1 += dth.powi(4) * (1. / -10780.); + k1 += dth.powi(6) * 2.769178184818219e-07; + k1 += k0.powi(2) * (1. / -10.); + k1 += dth.powi(2) * k0.powi(2) * (1. / 4200.); + k1 += dth.powi(4) * k0.powi(2) * 1.6959677820260655e-05; + k1 += k0.powi(4) * (1. / -1400.); + k1 += dth.powi(2) * k0.powi(4) * 6.84915970574303e-05; + k1 += k0.powi(6) * -7.936475029053326e-06; + let mut ch = 1.; + ch += dth.powi(2) * (1. / -40.); + ch += dth.powi(4) * 0.00034226190482569864; + ch += dth.powi(6) * -1.9349474568904524e-06; + ch += k0.powi(2) * (1. / -24.); + ch += dth.powi(2) * k0.powi(2) * 0.0024702380951963226; + ch += dth.powi(4) * k0.powi(2) * -3.7297408997537985e-05; + ch += k0.powi(4) * (1. / 1920.); + ch += dth.powi(2) * k0.powi(4) * -4.87350869747975e-05; + ch += k0.powi(6) * -3.1001936068463107e-06; + (k1 * dth, ch) +} + +// About 3ns +#[bench] +fn bench_fast_fit_euler(b: &mut Bencher) { + b.iter(|| fast_fit_euler(black_box(0.1), black_box(0.2))) +} + +// About 2ns +#[bench] +fn bench_fast_fit_euler2(b: &mut Bencher) { + b.iter(|| fast_fit_euler2(black_box(0.1), black_box(0.2))) +} + +// About 7ns +#[bench] +fn bench_fast_fit_euler3(b: &mut Bencher) { + b.iter(|| fast_fit_euler3(black_box(0.1), black_box(0.2))) +} diff --git a/src/euler.rs b/src/euler.rs index 4066de1c..f8d7616b 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -122,34 +122,6 @@ fn integ_euler_12(k0: f64, k1: f64) -> (f64, f64) { (u, v) } -/// A cheaper approximation to the Euler spiral integral. This might be useful -/// if we have lots of low-deflection segments, but it's not clear it'll be a -/// whole lot faster. -/* -pub fn integ_euler_8(k0: f64, k1: f64) -> (f64, f64) { - let t1_1 = k0; - let t1_2 = 0.5 * k1; - let t2_2 = t1_1 * t1_1; - let t2_3 = 2. * (t1_1 * t1_2); - let t2_4 = t1_2 * t1_2; - let t3_4 = t2_2 * t1_2 + t2_3 * t1_1; - let t3_6 = t2_4 * t1_2; - let t4_4 = t2_2 * t2_2; - let t4_5 = 2. * (t2_2 * t2_3); - let t4_6 = 2. * (t2_2 * t2_4) + t2_3 * t2_3; - let t5_6 = t4_4 * t1_2 + t4_5 * t1_1; - let t6_6 = t4_4 * t2_2; - let mut u = 1.; - u -= (1./24.) * t2_2 + (1./160.) * t2_4; - u += (1./1920.) * t4_4 + (1./10752.) * t4_6; - u -= (1./322560.) * t6_6; - let mut v = (1./12.) * t1_2; - v -= (1./480.) * t3_4 + (1./2688.) * t3_6; - v += (1./53760.) * t5_6; - (u, v) -} -*/ - #[doc(hidden)] /// Computation of the Euler spiral integral using subdivision. pub fn integ_euler_12n(mut k0: f64, mut k1: f64, n: usize) -> (f64, f64) { @@ -193,20 +165,19 @@ pub fn integ_euler_12n(mut k0: f64, mut k1: f64, n: usize) -> (f64, f64) { /// are validated in the notebook attached to the parallel curve blog post. /// /// [Raph's thesis]: https://www.levien.com/phd/thesis.pdf -#[inline] pub fn integ_euler(k0: f64, k1: f64, accuracy: f64) -> (f64, f64) { - let thresh = accuracy.powf(1.0 / 6.0); - integ_euler_thresh(k0, k1, thresh) -} - -fn integ_euler_thresh(k0: f64, k1: f64, thresh: f64) -> (f64, f64) { let c1 = k1.abs(); let c0 = k0.abs() + 0.5 * c1; let est_err_raw = 0.006 * c0 * c0 + 0.029 * c1; - if est_err_raw < thresh { + // Fun performance note: if the accuracy were always known at compile time, + // it would be theoretically cheaper to compare against accuracy^(1/6), which + // is computed anyway in the subdivision case. But the cost of the powi(6) is + // basically not measurable, and the cost of the ^(1/6) is ballpark double + // the integration itself. + if est_err_raw.powi(6) < accuracy { integ_euler_12(k0, k1) } else { - let n = (est_err_raw / thresh).ceil() as usize; + let n = (est_err_raw / accuracy.powf(1.0 / 6.0)).ceil() as usize; integ_euler_12n(k0, k1, n) } } From 58ef380586cc477bfab045cb84295bff7874e1ac Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Wed, 17 Feb 2021 13:02:55 -0800 Subject: [PATCH 08/12] Slight cleanup of types and signatures Rename "FitEulerResult" to "EulerParams". --- examples/euler.rs | 82 +++++++++++++++++++++++++++++++++++++++-------- src/euler.rs | 72 +++++++++++++++++++++++------------------ 2 files changed, 110 insertions(+), 44 deletions(-) diff --git a/examples/euler.rs b/examples/euler.rs index e60ed29c..c89d420d 100644 --- a/examples/euler.rs +++ b/examples/euler.rs @@ -11,8 +11,7 @@ // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License. - -use kurbo::{Affine, CubicBez, EulerSeg, Point, Shape}; +use kurbo::{Affine, CubicBez, EulerParams, EulerSeg, Point, Shape, Vec2}; #[allow(unused)] fn check_euler_int_accuracy() { @@ -29,16 +28,6 @@ fn check_euler_int_accuracy() { } } -// errors for 1e-9 threshold -// n == 1: 1.01 -// n == 2: 3.96 -// n == 3: 8.23 -// n == 4: 14.2 - -// integ_euler_8 -// est = 0.112: 1e-9 -// 0.62: 1e-6 - fn rand_cubic() -> CubicBez { let p0 = Point::new(0., 0.); let p1 = Point::new( @@ -86,6 +75,73 @@ fn cubic_err_scatter() { ); } +fn fit_cubic_plot() { + const N: usize = 513; + println!("P3"); + println!("{} {}", N, N); + println!("255"); + const KMAX: f64 = 0.1; + for y in 0..N { + let k1 = KMAX * ((y as f64) * (2.0 / ((N - 1) as f64)) - 1.0); + let k1 = 4. * k1; + for x in 0..N { + let k0 = KMAX * ((x as f64) * (2.0 / ((N - 1) as f64)) - 1.0); + let params = EulerParams::from_k0_k1(k0, k1); + let es = EulerSeg::from_params(Point::new(0., 0.), Point::new(1., 0.), params); + let th0 = -params.th(0.0); + let th1 = params.th(1.0); + let v0 = Vec2::from_angle(th0); + let p0 = Point::new(0., 0.); + let p1 = p0 + 2. / 3. / (1. + v0.x) * v0; + let p3 = Point::new(1., 0.); + let v1 = Vec2::from_angle(-th1); + let p2 = p3 - 2. / 3. / (1. + v1.x) * v1; + let c = CubicBez::new(p0, p1, p2, p3); + let err = es.cubic_euler_err(c, 10); + let mut est = 0.0; + est += 1.5e-5 * k0.powi(5).abs(); + est += 3e-6 * k1.powi(3).abs(); + est += 6e-4 * (k0 * k0 * k1).abs(); + est += 10e-5 * (k0 * k1 * k1).abs(); + //est += 5e-3 * (k0 * k0 * k1 * k1).abs(); + let err = err.powf(0.33); + let est = est.powf(0.33); + let escale = 1.5e4; + let r = escale * err; + let g = escale * est; + let b = g; + println!( + "{} {} {}", + (r.round() as u32).min(255), + (g.round() as u32).min(255), + (b.round() as u32).min(255) + ); + //println!("{} {}: {} {}", k0, k1, err, est); + } + } +} + +fn arc_toy() { + let k0 = 0.2; + let k1 = 0.0; + let params = EulerParams::from_k0_k1(k0, k1); + let es = EulerSeg::from_params(Point::new(0., 0.), Point::new(1., 0.), params); + let th0 = -params.th(0.0); + let th1 = params.th(1.0); + let v0 = Vec2::from_angle(th0); + let p0 = Point::new(0., 0.); + let p1 = p0 + 2. / 3. / (1. + v0.x) * v0; + let p3 = Point::new(1., 0.); + let v1 = Vec2::from_angle(-th1); + let p2 = p3 - 2. / 3. / (1. + v1.x) * v1; + let c = CubicBez::new(p0, p1, p2, p3); + println!("{:?}", c); + let err = es.cubic_euler_err(c, 10); + println!("err = {:e}", err); +} + fn main() { - cubic_err_scatter() + //cubic_err_scatter() + fit_cubic_plot(); + //arc_toy(); } diff --git a/src/euler.rs b/src/euler.rs index f8d7616b..65d76674 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -22,7 +22,7 @@ use crate::{ pub struct EulerSeg { p0: Point, p1: Point, - es: FitEulerResult, + params: EulerParams, } /// The derivative of an Euler spiral segment. @@ -45,7 +45,7 @@ pub struct EulerSegDeriv2(EulerSegDeriv); /// It's entirely possible the disposition of this is to be inlined into `EulerSeg`. /// I'm not sure it's useful by itself. #[derive(Clone, Copy, Debug)] -pub struct FitEulerResult { +pub struct EulerParams { k0: f64, k1: f64, chord: f64, @@ -63,7 +63,7 @@ pub enum EulerPathEl { /// A line segment to the given point. LineTo(Point), /// An Euler spiral segment to the given point. - EulerTo(FitEulerResult, Point), + EulerTo(EulerParams, Point), /// Close the subpath. ClosePath, } @@ -182,14 +182,14 @@ pub fn integ_euler(k0: f64, k1: f64, accuracy: f64) -> (f64, f64) { } } -impl FitEulerResult { +impl EulerParams { /// Find the Euler spiral parameters for the given deflection. /// /// TODO: use research for direct solution. /// /// Discussion question: should this take an accuracy parameter? /// This version basically hardcodes 1e-9. - pub fn fit_euler(th0: f64, th1: f64) -> FitEulerResult { + pub fn fit_euler(th0: f64, th1: f64) -> EulerParams { // Note: we could skip the solving for very small deflection let mut k1_old = 0.0; let dth = th1 - th0; @@ -202,7 +202,7 @@ impl FitEulerResult { let error = dth - (0.25 * k1 - 2.0 * chth); if error.abs() < 1e-9 { let chord = u.hypot(v); - return FitEulerResult { + return EulerParams { k0, k1, chord, @@ -217,12 +217,12 @@ impl FitEulerResult { panic!("fit_euler diverged on {}, {}", th0, th1); } - /// Create a `FitEulerResult` from k0 and k1 parameters. - pub fn from_k0_k1(k0: f64, k1: f64) -> FitEulerResult { + /// Create a `EulerParams` from k0 and k1 parameters. + pub fn from_k0_k1(k0: f64, k1: f64) -> EulerParams { let (u, v) = integ_euler(k0, k1, 1e-12); let chth = v.atan2(u); let chord = u.hypot(v); - FitEulerResult { + EulerParams { k0, k1, chord, @@ -234,7 +234,7 @@ impl FitEulerResult { /// /// The sign may be confusing, but it matches the spiro code. When `t = 0`, /// the result is `-th0`, and when `t = 1`, the result is `th1`. - fn th(&self, t: f64) -> f64 { + pub fn th(&self, t: f64) -> f64 { let u = t - 0.5; (0.5 * self.k1 * u + self.k0) * u - self.chth } @@ -260,8 +260,8 @@ impl EulerSeg { /// /// TODO: document the conventions. An SVG would be especially nice. pub fn new(p0: Point, p1: Point, th0: f64, th1: f64) -> EulerSeg { - let es = FitEulerResult::fit_euler(th0, th1); - EulerSeg { p0, p1, es } + let params = EulerParams::fit_euler(th0, th1); + EulerSeg { p0, p1, params } } /// Create an Euler segment from a cubic Bézier. @@ -274,11 +274,21 @@ impl EulerSeg { let d03 = c.p3 - c.p0; let th0 = d03.cross(d01).atan2(d03.dot(d01)); let th1 = d23.cross(d03).atan2(d23.dot(d03)); - let es = FitEulerResult::fit_euler(th0, th1); + let params = EulerParams::fit_euler(th0, th1); EulerSeg { p0: c.p0, p1: c.p3, - es, + params, + } + } + + /// Create a segment from params and endpoints. + /// + /// Mostly used for experimentation. + #[doc(hidden)] + pub fn from_params(p0: Point, p1: Point, params: EulerParams) -> EulerSeg { + EulerSeg { + p0, p1, params } } @@ -312,7 +322,7 @@ impl EulerSeg { /// Report whether the segment is a straight line. pub fn is_line(&self) -> bool { - self.es.k0 == 0.0 && self.es.k1 == 0.0 + self.params.k0 == 0.0 && self.params.k1 == 0.0 } /// Convert to cubic beziers. @@ -322,7 +332,7 @@ impl EulerSeg { let mut dt = 1.0; let mut p0 = self.p0; let chord_atan = (self.p1 - self.p0).atan2(); - let thresh = accuracy * self.es.chord / (self.p1 - self.p0).hypot(); + let thresh = accuracy * self.params.chord / (self.p1 - self.p0).hypot(); std::iter::from_fn(move || { let t0 = (t0_int as f64) * dt; if t0 == 1.0 { @@ -330,11 +340,11 @@ impl EulerSeg { } loop { let t1 = t0 + dt; - let k0 = dt * (this.es.k0 + 0.5 * (t0 + t1 - 1.0) * this.es.k1); - let k1 = dt * dt * this.es.k1; + let k0 = dt * (this.params.k0 + 0.5 * (t0 + t1 - 1.0) * this.params.k1); + let k1 = dt * dt * this.params.k1; let a0 = k0.abs(); let a1 = k1.abs(); - // Error metric empirically determined. I can make a Python notebook available. + // Error metric empirically determined, using `fit_cubic_plot` in example. let err = 1.5e-5 * a0.powi(5) + 6e-4 * a0 * a0 * a1 + 1e-4 * a0 * a1 * a1 @@ -350,8 +360,8 @@ impl EulerSeg { // Note: it's possible to this with rotation and normalization, // avoiding the trig. let d_atan = chord_atan - dp.atan2(); - let th0 = d_atan - this.es.th(t0); - let th1 = -d_atan + this.es.th(t1); + let th0 = d_atan - this.params.th(t0); + let th1 = -d_atan + this.params.th(t1); let v0 = Vec2::from_angle(th0); let c0 = Point::new(0., 0.); let c1 = c0 + 2. / 3. / (1. + v0.x) * v0; @@ -379,7 +389,7 @@ impl ParamCurve for EulerSeg { fn eval(&self, t: f64) -> Point { // The accuracy here is somewhat arbitrary, but should be adequate // for most work, and not entail loss of efficiency. - let Point { x, y } = self.es.eval(t, 1e-9); + let Point { x, y } = self.params.eval(t, 1e-9); let chord = self.p1 - self.p0; Point::new( self.p0.x + chord.x * x - chord.y * y, @@ -393,10 +403,10 @@ impl ParamCurve for EulerSeg { let p0 = self.eval(range.start); let p1 = self.eval(range.end); let dt = range.end - range.start; - let k0 = dt * (self.es.k0 + 0.5 * (range.start + range.end - 1.0) * self.es.k1); - let k1 = dt * dt * self.es.k1; - let es = FitEulerResult::from_k0_k1(k0, k1); - EulerSeg { p0, p1, es } + let k0 = dt * (self.params.k0 + 0.5 * (range.start + range.end - 1.0) * self.params.k1); + let k1 = dt * dt * self.params.k1; + let params = EulerParams::from_k0_k1(k0, k1); + EulerSeg { p0, p1, params } } fn start(&self) -> Point { @@ -413,14 +423,14 @@ impl ParamCurveArclen for EulerSeg { /// /// Note that this implementation is fast and accurate. fn arclen(&self, _accuracy: f64) -> f64 { - (self.p1 - self.p0).hypot() / self.es.chord + (self.p1 - self.p0).hypot() / self.params.chord } /// The parameter that results in the given arc length. /// /// This implementation is also fast and accurate. fn inv_arclen(&self, arclen: f64, _accuracy: f64) -> f64 { - arclen * self.es.chord / (self.p1 - self.p0).hypot() + arclen * self.params.chord / (self.p1 - self.p0).hypot() } } @@ -428,7 +438,7 @@ impl ParamCurveDeriv for EulerSeg { type DerivResult = EulerSegDeriv; fn deriv(&self) -> Self::DerivResult { - let FitEulerResult { k0, k1, chth, .. } = self.es; + let EulerParams { k0, k1, chth, .. } = self.params; EulerSegDeriv { c0: 0.5 * k0 - 0.125 * k1 + chth + (self.p1 - self.p0).atan2(), c1: -k0 + 0.5 * k1, @@ -440,7 +450,7 @@ impl ParamCurveDeriv for EulerSeg { impl ParamCurveCurvature for EulerSeg { fn curvature(&self, t: f64) -> f64 { - (self.es.k0 + (t - 0.5) * self.es.k1) * self.es.chord / (self.p1 - self.p0).hypot() + (self.params.k0 + (t - 0.5) * self.params.k1) * self.params.chord / (self.p1 - self.p0).hypot() } } @@ -490,7 +500,7 @@ impl From for EulerSeg { EulerSeg { p0: l.p0, p1: l.p1, - es: FitEulerResult { + params: EulerParams { k0: 0., k1: 0., chord: 1., From 462098cf2f274ea48db1b2516b7c7dea37dcc077 Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Wed, 17 Feb 2021 15:58:44 -0800 Subject: [PATCH 09/12] Make Euler examples accessible Also make them easier to find and improve docs a bit. --- Cargo.toml | 2 +- README.md | 4 ++++ examples/euler.rs | 20 +++++++++++++++++--- src/euler.rs | 10 ++++++++-- 4 files changed, 30 insertions(+), 6 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 30f210c5..21ee4650 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,7 +11,7 @@ readme = "README.md" categories = ["graphics"] [package.metadata.docs.rs] -features = ["mint", "serde"] +features = ["mint", "serde", "euler"] [dependencies] arrayvec = "0.5.1" diff --git a/README.md b/README.md index 6bf3868f..ca6fe4bf 100644 --- a/README.md +++ b/README.md @@ -23,6 +23,10 @@ Here we mention a few other curves libraries and touch on some of the decisions Some code has been copied from lyon_geom with adaptation, thus the author of lyon_geom, Nicolas Silva, is credited in the [AUTHORS] file. +## Euler spirals + +In addition to the Bézier functionality, kurbo also contains support for Euler spirals. Enable the `euler` feature to get these. + ## More info To learn more about Bézier curves, [A Primer on Bézier Curves] by Pomax is indispensable. diff --git a/examples/euler.rs b/examples/euler.rs index c89d420d..57ae60a7 100644 --- a/examples/euler.rs +++ b/examples/euler.rs @@ -11,6 +11,9 @@ // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License. + +//! A few test programs and experiments to exercise Euler segments. + use kurbo::{Affine, CubicBez, EulerParams, EulerSeg, Point, Shape, Vec2}; #[allow(unused)] @@ -141,7 +144,18 @@ fn arc_toy() { } fn main() { - //cubic_err_scatter() - fit_cubic_plot(); - //arc_toy(); + if let Some(cmd) = std::env::args().skip(1).next() { + match cmd.as_str() { + "cubic_err_scatter" => cubic_err_scatter(), + "fit_cubic_plot" => fit_cubic_plot(), + "arc_toy" => arc_toy(), + _ => println!("unknown cmd"), + } + println!("arg = {}", cmd); + } else { + println!("usage: euler "); + println!(" cubic_err_scatter: scatter plot of cubic->ES error estimates"); + println!(" fit_cubic_plot: 2d image of ES->cubic error"); + println!(" arc_toy: something that prints single ES->cubic fit"); + } } diff --git a/src/euler.rs b/src/euler.rs index 65d76674..a0e46187 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -18,6 +18,8 @@ use crate::{ }; /// An Euler spiral segment. +/// +/// This is only enabled when the `euler` feature is selected. #[derive(Clone, Copy, Debug)] pub struct EulerSeg { p0: Point, @@ -39,7 +41,7 @@ pub struct EulerSegDeriv2(EulerSegDeriv); /// Parameters for an Euler spiral segment. Does not include endpoint geometry. /// -/// This is something of an internal detail for `EulerSeg` and might not make +/// This is something of an internal detail for [`EulerSeg`] and might not make /// it to the public interface. It's public here for experimentation. /// /// It's entirely possible the disposition of this is to be inlined into `EulerSeg`. @@ -53,6 +55,10 @@ pub struct EulerParams { } /// A path consisting of piecewise Euler spiral segments. +/// +/// TODO: develop this further, including implementing the [`Shape`][crate::Shape] trait. +/// +/// This is only enabled when the `euler` feature is selected. pub struct EulerPath(Vec); /// An element of a piecewise Euler spiral path. @@ -217,7 +223,7 @@ impl EulerParams { panic!("fit_euler diverged on {}, {}", th0, th1); } - /// Create a `EulerParams` from k0 and k1 parameters. + /// Create `EulerParams` from k0 and k1 parameters. pub fn from_k0_k1(k0: f64, k1: f64) -> EulerParams { let (u, v) = integ_euler(k0, k1, 1e-12); let chth = v.atan2(u); From d7a85f28dc2e14524ba44ba92a3a3ae6e2490d58 Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Thu, 18 Feb 2021 16:35:08 -0800 Subject: [PATCH 10/12] Parallel curves A pretty good implementation of parallel curves, including handling the cusp, and applying the closed-form error metric for subdivision. --- benches/euler.rs | 16 +++--- examples/euler.rs | 87 ++++++++++++++++++++++++++++--- src/euler.rs | 129 +++++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 211 insertions(+), 21 deletions(-) diff --git a/benches/euler.rs b/benches/euler.rs index 99c707aa..03cdade0 100644 --- a/benches/euler.rs +++ b/benches/euler.rs @@ -17,7 +17,7 @@ #![cfg(nightly)] #![feature(test)] extern crate test; -use test::{Bencher, black_box}; +use test::{black_box, Bencher}; use kurbo::*; @@ -53,12 +53,12 @@ fn integ_euler_8(k0: f64, k1: f64) -> (f64, f64) { let t5_6 = t4_4 * t1_2 + t4_5 * t1_1; let t6_6 = t4_4 * t2_2; let mut u = 1.; - u -= (1./24.) * t2_2 + (1./160.) * t2_4; - u += (1./1920.) * t4_4 + (1./10752.) * t4_6; - u -= (1./322560.) * t6_6; - let mut v = (1./12.) * t1_2; - v -= (1./480.) * t3_4 + (1./2688.) * t3_6; - v += (1./53760.) * t5_6; + u -= (1. / 24.) * t2_2 + (1. / 160.) * t2_4; + u += (1. / 1920.) * t4_4 + (1. / 10752.) * t4_6; + u -= (1. / 322560.) * t6_6; + let mut v = (1. / 12.) * t1_2; + v -= (1. / 480.) * t3_4 + (1. / 2688.) * t3_6; + v += (1. / 53760.) * t5_6; (u, v) } @@ -75,7 +75,7 @@ fn bench_fit_euler(b: &mut Bencher) { } // Pretty much straight out of the notebook. Looking at the generated -// asm, the compiler does a good job +// asm, the compiler does a good job fn fast_fit_euler(th0: f64, th1: f64) -> f64 { let k0 = th0 + th1; let dth = th1 - th0; diff --git a/examples/euler.rs b/examples/euler.rs index 57ae60a7..b5fdf43a 100644 --- a/examples/euler.rs +++ b/examples/euler.rs @@ -14,7 +14,7 @@ //! A few test programs and experiments to exercise Euler segments. -use kurbo::{Affine, CubicBez, EulerParams, EulerSeg, Point, Shape, Vec2}; +use kurbo::{Affine, BezPath, CubicBez, EulerParams, EulerSeg, ParamCurve, Point, Shape, Vec2}; #[allow(unused)] fn check_euler_int_accuracy() { @@ -50,7 +50,7 @@ fn cubic_err_scatter() { r##" - + "## ); for _ in 0..100000 { @@ -65,10 +65,7 @@ fn cubic_err_scatter() { if y < 750. { let a = Affine::new([20., 0., 0., 20., x, y]); let path = (a * c).into_path(1e-6); - println!( - r##""##, - path.to_svg() - ); + println!("", path.to_svg()); } } println!( @@ -143,19 +140,95 @@ fn arc_toy() { println!("err = {:e}", err); } +fn print_svg_path(bp: &BezPath, color: &str) { + println!( + r" ", + bp.to_svg(), + color + ); +} + +fn print_euler_seg(es: &EulerSeg, color: &str) { + let mut bp = BezPath::new(); + bp.move_to(es.start()); + bp.extend(es.to_cubics(0.1)); + print_svg_path(&bp, color); +} + +fn parallel() { + println!( + r##" + + + "## + ); + let es = EulerSeg::new(Point::new(100., 300.), Point::new(500., 300.), 0.9, 1.1); + print_euler_seg(&es, "black"); + + let offset = -250.0; + let es_par = es.parallel_approx(offset, -1.); + print_euler_seg(&es_par, "#c00"); + + let mut bp = BezPath::new(); + for seg in es.parallel_curve(offset, 1e-1) { + if bp.is_empty() { + bp.move_to(seg.start()); + } + bp.extend(seg.to_cubics(0.1)); + } + print_svg_path(&bp, "#00c"); + + println!( + r#" + +"# + ); +} + +const COLORS: &[&str] = &["#00f", "#aaf"]; + +fn parallel_multi() { + println!( + r##" + + + "## + ); + let es = EulerSeg::new(Point::new(300., 600.), Point::new(700., 600.), 0.7, 1.1); + //let es = EulerSeg::from_params(Point::new(300., 400.), Point::new(700., 400.), EulerParams::from_k0_k1(5.0, 10.0)); + print_euler_seg(&es, "black"); + + for i in 1..25 { + let offset = -25.0 * i as f64; + + for (j, seg) in es.parallel_curve(offset, 1e-1).enumerate() { + print_euler_seg(&seg, COLORS[j % COLORS.len()]); + } + } + + println!( + r#" + +"# + ); +} + fn main() { if let Some(cmd) = std::env::args().skip(1).next() { match cmd.as_str() { "cubic_err_scatter" => cubic_err_scatter(), "fit_cubic_plot" => fit_cubic_plot(), "arc_toy" => arc_toy(), + "parallel" => parallel(), + "parallel_multi" => parallel_multi(), _ => println!("unknown cmd"), } - println!("arg = {}", cmd); } else { println!("usage: euler "); println!(" cubic_err_scatter: scatter plot of cubic->ES error estimates"); println!(" fit_cubic_plot: 2d image of ES->cubic error"); println!(" arc_toy: something that prints single ES->cubic fit"); + println!(" parallel: experiment with parallel curves"); + println!(" parallel_multi: multiple offsets from one ES"); } } diff --git a/src/euler.rs b/src/euler.rs index a0e46187..59d706c7 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -293,9 +293,7 @@ impl EulerSeg { /// Mostly used for experimentation. #[doc(hidden)] pub fn from_params(p0: Point, p1: Point, params: EulerParams) -> EulerSeg { - EulerSeg { - p0, p1, params - } + EulerSeg { p0, p1, params } } /// Calculate error from cubic bezier. @@ -389,6 +387,112 @@ impl EulerSeg { } }) } + + /// Generate an approximate parallel curve as a single segment. + /// + /// The `dir` argument should be -1.0 when the direction is reversed + /// due to the curvature exceeding the offset; in other words it should + /// have opposite sign on either side of a cusp. + pub fn parallel_approx(&self, offset: f64, dir: f64) -> EulerSeg { + let th0 = self.params.th(0.0); + let th1 = self.params.th(1.0); + let v0 = Vec2::new(offset * th0.sin(), offset * th0.cos()); + let v1 = Vec2::new(offset * th1.sin(), offset * th1.cos()); + let chord = (self.p1 - self.p0).hypot(); + let dth = dir * (v1.y - v0.y).atan2(dir * (chord + v1.x - v0.x)); + let c = (self.p1.x - self.p0.x) / chord; + let s = (self.p1.y - self.p0.y) / chord; + let p0 = self.p0 + Vec2::new(v0.x * c - v0.y * s, v0.x * s + v0.y * c); + let p1 = self.p1 + Vec2::new(v1.x * c - v1.y * s, v1.x * s + v1.y * c); + EulerSeg::new(p0, p1, -th0 - dth, th1 + dth) + } + + /// Generate a parallel curve as a sequence of segments. + pub fn parallel_curve(&self, offset: f64, accuracy: f64) -> impl Iterator { + let chord = (self.p1 - self.p0).hypot(); + let a = self.params.k0 * self.params.chord + chord / offset; + let b = 0.5 * self.params.k1 * self.params.chord; + let u0 = a - b; + let u1 = a + b; + let (es0, mut es1) = if u0.signum() * u1.signum() < 0.0 { + let t_split = u0 / (u0 - u1); + let es0 = self.subsegment(0.0..t_split); + let es1 = self.subsegment(t_split..1.0); + (es0, Some(es1)) + } else { + (*self, None) + }; + let mut inner = es0.parallel_curve_raw(offset, accuracy); + std::iter::from_fn(move || { + if let Some(seg) = inner.next() { + Some(seg) + } else { + let es1 = es1.take()?; + inner = es1.parallel_curve_raw(offset, accuracy); + inner.next() + } + }) + } + + /// Generate a parallel curve assuming no cusp. + fn parallel_curve_raw(&self, offset: f64, accuracy: f64) -> impl Iterator { + let this = *self; + + // Estimate error for a single-segment approximation + let chord = (this.p1 - this.p0).hypot(); + // TODO: this predicts the L2 norm error, should rescale to accurately predict + // Hausdorff norm. From quick experimentation, that's probably multiplying by + // about 1.5, but it should be measured more carefully. + let arc = chord / this.params.chord; + let a = this.params.k0 * this.params.chord + chord / offset; + let dir = a.signum() * offset.signum(); + let est_err = 0.005 * (this.params.k1.powi(2) / a).abs() * arc; + let b = 0.5 * this.params.k1 * this.params.chord; + let u0 = a - b; + let u1 = a + b; + let u0 = u0.abs(); + let u1 = u1.abs(); + let (n, c0, dc, du_recip); + if est_err < accuracy { + n = 1; + c0 = 0.0; + dc = 0.0; + du_recip = 0.0; + } else { + n = (est_err * predict_rel(u0, u1) / accuracy).powf(0.25).ceil() as usize; + c0 = u0.powf(0.75); + dc = (u1.powf(0.75) - c0) / (n as f64); + du_recip = (u1 - u0).recip(); + }; + println!( + "", + est_err, u0, u1, n + ); + let mut t0 = 0.0; + let mut i = 0; + std::iter::from_fn(move || { + if i == n { + None + } else { + i += 1; + let t1 = if i == n { + 1.0 + } else { + ((c0 + (i as f64) * dc).powf(4. / 3.) - u0) * du_recip + }; + let seg_par = this.subsegment(t0..t1).parallel_approx(offset, dir); + t0 = t1; + Some(seg_par) + } + }) + } +} + +/// The error of the worst subdivision compared to the error of the +/// whole segment, relative to that predicted by n^4 scaling. +fn predict_rel(t0: f64, t1: f64) -> f64 { + let a = t0.min(t1) / t0.max(t1); + (128_f64 / 81.).powf(0.25) * ((1.0 - a.powf(0.75)) / (1.0 - a)).powi(4) * (1.0 + a) } impl ParamCurve for EulerSeg { @@ -403,8 +507,6 @@ impl ParamCurve for EulerSeg { ) } - // TODO: I'm not sure this is right, possibly chord needs to be taken into - // account too? Need to test. fn subsegment(&self, range: std::ops::Range) -> Self { let p0 = self.eval(range.start); let p1 = self.eval(range.end); @@ -456,7 +558,8 @@ impl ParamCurveDeriv for EulerSeg { impl ParamCurveCurvature for EulerSeg { fn curvature(&self, t: f64) -> f64 { - (self.params.k0 + (t - 0.5) * self.params.k1) * self.params.chord / (self.p1 - self.p0).hypot() + (self.params.k0 + (t - 0.5) * self.params.k1) * self.params.chord + / (self.p1 - self.p0).hypot() } } @@ -562,6 +665,20 @@ mod tests { ParamCurveDeriv, ParamCurveNearest, PathEl, Point, }; + #[test] + fn euler_subsegment() { + let es = EulerSeg::new(Point::ORIGIN, Point::new(1., 1.), 0.2, 0.3); + let t0 = 0.3; + let t1 = 0.8; + let subseg = es.subsegment(t0..t1); + for i in 0..11 { + let t = (i as f64) * 0.1; + let orig_p = es.eval(t0 + (t1 - t0) * t); + let subseg_p = subseg.eval(t); + assert!((orig_p - subseg_p).hypot() < 1e-11); + } + } + #[test] fn euler_deriv() { let es = EulerSeg::new(Point::ORIGIN, Point::new(1., 1.), 0.2, 0.3); From 91a6c9b7c32a2d7c3b048e30493b2249338c5ad0 Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Thu, 18 Feb 2021 23:03:50 -0800 Subject: [PATCH 11/12] Make clippy happy --- examples/euler.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/euler.rs b/examples/euler.rs index b5fdf43a..9574d406 100644 --- a/examples/euler.rs +++ b/examples/euler.rs @@ -214,7 +214,7 @@ fn parallel_multi() { } fn main() { - if let Some(cmd) = std::env::args().skip(1).next() { + if let Some(cmd) = std::env::args().nth(1) { match cmd.as_str() { "cubic_err_scatter" => cubic_err_scatter(), "fit_cubic_plot" => fit_cubic_plot(), From cbbc40c2155c8ad963899205d15b8116fad3511c Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Fri, 19 Feb 2021 09:49:57 -0800 Subject: [PATCH 12/12] More pictures in the examples This commit captures the source code for the illustrations in the parallel curve blog post. --- examples/euler.rs | 99 +++++++++++++++++++++++++++++++++-------------- src/euler.rs | 25 +++++------- 2 files changed, 78 insertions(+), 46 deletions(-) diff --git a/examples/euler.rs b/examples/euler.rs index 9574d406..183f5b16 100644 --- a/examples/euler.rs +++ b/examples/euler.rs @@ -148,26 +148,26 @@ fn print_svg_path(bp: &BezPath, color: &str) { ); } -fn print_euler_seg(es: &EulerSeg, color: &str) { +fn print_euler_seg(es: &EulerSeg, color: &str, accuracy: f64) { let mut bp = BezPath::new(); bp.move_to(es.start()); - bp.extend(es.to_cubics(0.1)); + bp.extend(es.to_cubics(accuracy)); print_svg_path(&bp, color); } -fn parallel() { +// Note: this is the code that procedurally draws the figures in the +// "Cleaner parallel curves with Euler spirals" blog post. + +fn parallel_1() { println!( - r##" - - - "## + "" ); - let es = EulerSeg::new(Point::new(100., 300.), Point::new(500., 300.), 0.9, 1.1); - print_euler_seg(&es, "black"); + let es = EulerSeg::new(Point::new(25., 140.), Point::new(425., 140.), 0.7, 1.1); + print_euler_seg(&es, "black", 0.1); - let offset = -250.0; - let es_par = es.parallel_approx(offset, -1.); - print_euler_seg(&es_par, "#c00"); + let offset = -150.0; + let es_par = es.parallel_approx(offset, 1.); + print_euler_seg(&es_par, "#c00", 0.1); let mut bp = BezPath::new(); for seg in es.parallel_curve(offset, 1e-1) { @@ -178,39 +178,76 @@ fn parallel() { } print_svg_path(&bp, "#00c"); + println!(""); +} + +fn parallel_2() { println!( - r#" - -"# + "" + ); + let es = EulerSeg::from_params( + Point::new(180., 270.), + Point::new(330., 70.), + EulerParams::from_k0_k1(0.0, 40.0), ); + print_euler_seg(&es, "black", 0.1); + + let offset = -60.0; + + let mut bp = BezPath::new(); + for seg in es.parallel_curve(offset, 1e-1) { + if bp.is_empty() { + bp.move_to(seg.start()); + } + bp.extend(seg.to_cubics(0.1)); + } + print_svg_path(&bp, "#00c"); + + println!(""); } const COLORS: &[&str] = &["#00f", "#aaf"]; -fn parallel_multi() { +fn parallel_multi_1() { println!( - r##" - - - "## + "" ); - let es = EulerSeg::new(Point::new(300., 600.), Point::new(700., 600.), 0.7, 1.1); + let es = EulerSeg::new(Point::new(200., 515.), Point::new(600., 515.), 0.7, 1.1); //let es = EulerSeg::from_params(Point::new(300., 400.), Point::new(700., 400.), EulerParams::from_k0_k1(5.0, 10.0)); - print_euler_seg(&es, "black"); + print_euler_seg(&es, "black", 0.1); for i in 1..25 { let offset = -25.0 * i as f64; for (j, seg) in es.parallel_curve(offset, 1e-1).enumerate() { - print_euler_seg(&seg, COLORS[j % COLORS.len()]); + print_euler_seg(&seg, COLORS[j % COLORS.len()], 0.1); } } + println!(""); +} + +fn parallel_multi_2() { println!( - r#" - -"# + "" + ); + let es = EulerSeg::from_params( + Point::new(10., 370.), + Point::new(410., 370.), + EulerParams::from_k0_k1(5.0, 10.0), ); + let accuracy = 1e-5; + print_euler_seg(&es, "black", accuracy); + + for i in 1..25 { + let offset = -25.0 * i as f64; + + for (j, seg) in es.parallel_curve(offset, accuracy).enumerate() { + print_euler_seg(&seg, COLORS[j % COLORS.len()], accuracy); + } + } + + println!(""); } fn main() { @@ -219,8 +256,10 @@ fn main() { "cubic_err_scatter" => cubic_err_scatter(), "fit_cubic_plot" => fit_cubic_plot(), "arc_toy" => arc_toy(), - "parallel" => parallel(), - "parallel_multi" => parallel_multi(), + "parallel_1" => parallel_1(), + "parallel_2" => parallel_2(), + "parallel_multi_1" => parallel_multi_1(), + "parallel_multi_2" => parallel_multi_2(), _ => println!("unknown cmd"), } } else { @@ -228,7 +267,7 @@ fn main() { println!(" cubic_err_scatter: scatter plot of cubic->ES error estimates"); println!(" fit_cubic_plot: 2d image of ES->cubic error"); println!(" arc_toy: something that prints single ES->cubic fit"); - println!(" parallel: experiment with parallel curves"); - println!(" parallel_multi: multiple offsets from one ES"); + println!(" parallel_[12]: experiment with parallel curves"); + println!(" parallel_multi_[12]: multiple offsets from one ES"); } } diff --git a/src/euler.rs b/src/euler.rs index 59d706c7..38cb56e3 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -447,27 +447,20 @@ impl EulerSeg { let a = this.params.k0 * this.params.chord + chord / offset; let dir = a.signum() * offset.signum(); let est_err = 0.005 * (this.params.k1.powi(2) / a).abs() * arc; - let b = 0.5 * this.params.k1 * this.params.chord; - let u0 = a - b; - let u1 = a + b; - let u0 = u0.abs(); - let u1 = u1.abs(); - let (n, c0, dc, du_recip); - if est_err < accuracy { - n = 1; - c0 = 0.0; - dc = 0.0; - du_recip = 0.0; - } else { + let mut n = 1; + let mut c0 = 0.0; + let mut dc = 0.0; + let mut u0 = 0.0; + let mut du_recip = 0.0; + if est_err > accuracy { + let b = 0.5 * this.params.k1 * this.params.chord; + u0 = (a - b).abs(); + let u1 = (a + b).abs(); n = (est_err * predict_rel(u0, u1) / accuracy).powf(0.25).ceil() as usize; c0 = u0.powf(0.75); dc = (u1.powf(0.75) - c0) / (n as f64); du_recip = (u1 - u0).recip(); }; - println!( - "", - est_err, u0, u1, n - ); let mut t0 = 0.0; let mut i = 0; std::iter::from_fn(move || {