Apply robustness fixes to CPU shader as well This should make CPU and GPU shaders consistent. Some ideas were borrowed from #520, in particular the divide by zero of the chord length in offset calculation is avoided (though there is now refactoring so the robustness logic is applied once and the effective chord length is now stored in cubic_params). The rustdoc for CubicParams::from_points_derivs also documents a bit of the rationale for the current robustness logic.
diff --git a/src/cpu_shader/euler.rs b/src/cpu_shader/euler.rs index 17c6bc9..575ed5b 100644 --- a/src/cpu_shader/euler.rs +++ b/src/cpu_shader/euler.rs
@@ -9,25 +9,25 @@ use super::util::Vec2; use std::f32::consts::FRAC_PI_4; +// Threshold for tangents to be considered near zero length +pub const TANGENT_THRESH: f32 = 1e-6; + /// This struct contains parameters derived from a cubic Bézier for the /// purpose of fitting a G1 continuous Euler spiral segment and estimating /// the Fréchet distance. /// /// The tangent angles represent deviation from the chord, so that when they /// are equal, the corresponding Euler spiral is a circular arc. -/// -/// Similar, the control point distances are normalized to a chord, so that -/// for small angles values near 1/3 represent a smooth curve. #[derive(Debug)] pub struct CubicParams { /// Tangent angle relative to chord at start. pub th0: f32, /// Tangent angle relative to chord at end. pub th1: f32, - /// Distance of first control point. - pub d0: f32, - /// Distance of second control point. - pub d1: f32, + /// The effective chord length, always a robustly nonzero value. + pub chord_len: f32, + /// The estimated error between the source cubic and the proposed Euler spiral. + pub err: f32, } #[derive(Debug)] @@ -49,13 +49,39 @@ impl CubicParams { /// Compute parameters from endpoints and derivatives. /// - /// Robustness note: this function must be protected from being called when the - /// chord is near zero. + /// This function is designed to be robust across a wide range of inputs. In + /// particular, it splits between near-zero chord length and the happy path. + /// In the former case, the parameters for the Euler spiral would not be valid, + /// so it proposes a straight line and computes a pretty good (conservative) + /// estimate of the Fréchet distance between that line and the source cubic. + /// + /// Computing an accurate estimate here fixes two tricky cases: very short + /// lines, in which the error will be below threshold and the flatten logic will + /// output a single line segment without subdividing, and loop cases with a + /// short chord, in which case the error will exceed the threshold, and the + /// chords of the subdivided pieces will be longer. + /// + /// An additional case is the near-cusp where the proposed Euler spiral has + /// a 180 degree U-turn (or, more generally, one angle exceeds 90 degrees and + /// the other does not). In that case, the resulting Euler spiral is quite + /// well defined (with finite curvature, so that its offset will generate a + /// near-semicircle, preserving G1 continuity), but the analytic error + /// calculation would be a huge overestimate. In that case, we just return + /// a rough estimate of the distance between the chord and the spiral segment. pub fn from_points_derivs(p0: Vec2, p1: Vec2, q0: Vec2, q1: Vec2, dt: f32) -> Self { let chord = p1 - p0; - let length_squared = chord.length_squared(); - assert_ne!(length_squared, 0.0); - let scale = dt / length_squared; + let chord_squared = chord.length_squared(); + let chord_len = chord_squared.sqrt(); + if chord_squared < TANGENT_THRESH.powi(2) { + let chord_err = ((9. / 32.0) * (q0.length_squared() + q1.length_squared())).sqrt() * dt; + return CubicParams { + th0: 0.0, + th1: 0.0, + chord_len: TANGENT_THRESH, + err: chord_err, + }; + } + let scale = dt / chord_squared; let h0 = Vec2::new( q0.x * chord.x + q0.y * chord.y, q0.y * chord.x - q0.x * chord.y, @@ -70,64 +96,61 @@ let d1 = h1.length() * scale; // Robustness note: we may want to clamp the magnitude of the angles to // a bit less than pi. Perhaps here, perhaps downstream. - CubicParams { th0, th1, d0, d1 } - } - // Estimated error of GH to Euler spiral - // - // Return value is normalized to chord - to get actual error, multiply - // by chord. - pub fn est_euler_err(&self) -> f32 { - // Potential optimization: work with unit vector rather than angle - let cth0 = self.th0.cos(); - let cth1 = self.th1.cos(); - if cth0 * cth1 < 0.0 { - // Rationale: this happens when fitting a cusp or near-cusp with - // a near 180 degree u-turn. The actual ES is bounded in that case. - // Further subdivision won't reduce the angles if actually a cusp. - // + // Estimate error of geometric Hermite interpolation to Euler spiral. + let cth0 = th0.cos(); + let cth1 = th1.cos(); + let mut err = if cth0 * cth1 < 0.0 { // A value of 2.0 represents the approximate worst case distance // from an Euler spiral with 0 and pi tangents to the chord. It // is not very critical; doubling the value would result in one more // subdivision in effectively a binary search for the cusp, while too // small a value may result in the actual error exceeding the bound. - return 2.0; + 2.0 + } else { + // Protect against divide-by-zero. This happens with a double cusp, so + // should in the general case cause subdivisions. + let e0 = (2. / 3.) / (1.0 + cth0).max(1e-9); + let e1 = (2. / 3.) / (1.0 + cth1).max(1e-9); + let s0 = th0.sin(); + let s1 = th1.sin(); + // Note: some other versions take sin of s0 + s1 instead. Those are incorrect. + // Strangely, calibration is the same, but more work could be done. + let s01 = cth0 * s1 + cth1 * s0; + let amin = 0.15 * (2. * e0 * s0 + 2. * e1 * s1 - e0 * e1 * s01); + let a = 0.15 * (2. * d0 * s0 + 2. * d1 * s1 - d0 * d1 * s01); + let aerr = (a - amin).abs(); + let symm = (th0 + th1).abs(); + let asymm = (th0 - th1).abs(); + let dist = (d0 - e0).hypot(d1 - e1); + let ctr = 4.625e-6 * symm.powi(5) + 7.5e-3 * asymm * symm.powi(2); + let halo_symm = 5e-3 * symm * dist; + let halo_asymm = 7e-2 * asymm * dist; + /* + println!(" e0: {e0}"); + println!(" e1: {e1}"); + println!(" s0: {s0}"); + println!(" s1: {s1}"); + println!(" s01: {s01}"); + println!(" amin: {amin}"); + println!(" a: {a}"); + println!(" aerr: {aerr}"); + println!(" symm: {symm}"); + println!(" asymm: {asymm}"); + println!(" dist: {dist}"); + println!(" ctr: {ctr}"); + println!(" halo_symm: {halo_symm}"); + println!(" halo_asymm: {halo_asymm}"); + */ + ctr + 1.55 * aerr + halo_symm + halo_asymm + }; + err *= chord_len; + CubicParams { + th0, + th1, + chord_len, + err, } - // Protect against divide-by-zero. This happens with a double cusp, so - // should in the general case cause subdivisions. - let e0 = (2. / 3.) / (1.0 + cth0).max(1e-9); - let e1 = (2. / 3.) / (1.0 + cth1).max(1e-9); - let s0 = self.th0.sin(); - let s1 = self.th1.sin(); - // Note: some other versions take sin of s0 + s1 instead. Those are incorrect. - // Strangely, calibration is the same, but more work could be done. - let s01 = cth0 * s1 + cth1 * s0; - let amin = 0.15 * (2. * e0 * s0 + 2. * e1 * s1 - e0 * e1 * s01); - let a = 0.15 * (2. * self.d0 * s0 + 2. * self.d1 * s1 - self.d0 * self.d1 * s01); - let aerr = (a - amin).abs(); - let symm = (self.th0 + self.th1).abs(); - let asymm = (self.th0 - self.th1).abs(); - let dist = (self.d0 - e0).hypot(self.d1 - e1); - let ctr = 4.625e-6 * symm.powi(5) + 7.5e-3 * asymm * symm.powi(2); - let halo_symm = 5e-3 * symm * dist; - let halo_asymm = 7e-2 * asymm * dist; - /* - println!(" e0: {e0}"); - println!(" e1: {e1}"); - println!(" s0: {s0}"); - println!(" s1: {s1}"); - println!(" s01: {s01}"); - println!(" amin: {amin}"); - println!(" a: {a}"); - println!(" aerr: {aerr}"); - println!(" symm: {symm}"); - println!(" asymm: {asymm}"); - println!(" dist: {dist}"); - println!(" ctr: {ctr}"); - println!(" halo_symm: {halo_symm}"); - println!(" halo_asymm: {halo_asymm}"); - */ - ctr + 1.55 * aerr + halo_symm + halo_asymm } } @@ -204,10 +227,10 @@ ) } + // Note: offset provided is scaled so that 1 = chord length pub fn eval_with_offset(&self, t: f32, offset: f32) -> Vec2 { let chord = self.p1 - self.p0; - let scaled = offset / chord.length(); - let Vec2 { x, y } = self.params.eval_with_offset(t, scaled); + let Vec2 { x, y } = self.params.eval_with_offset(t, offset); Vec2::new( self.p0.x + chord.x * x - chord.y * y, self.p0.y + chord.x * y + chord.y * x,
diff --git a/src/cpu_shader/flatten.rs b/src/cpu_shader/flatten.rs index bfc7d3a..d43837c 100644 --- a/src/cpu_shader/flatten.rs +++ b/src/cpu_shader/flatten.rs
@@ -3,7 +3,9 @@ use std::f32::consts::FRAC_1_SQRT_2; -use super::euler::{espc_int_approx, espc_int_inv_approx, CubicParams, EulerParams, EulerSeg}; +use super::euler::{ + espc_int_approx, espc_int_inv_approx, CubicParams, EulerParams, EulerSeg, TANGENT_THRESH, +}; use super::util::{Transform, Vec2, ROBUST_EPSILON}; use crate::cpu_dispatch::CpuBinding; use vello_encoding::math::f16_to_f32; @@ -258,127 +260,109 @@ break; } log!("@@@ loop1: t0: {t0}, dt: {dt}"); - loop { - let mut t1 = t0 + dt; - let this_p0 = last_p; - let this_q0 = last_q; - let (mut this_p1, mut this_q1) = eval_cubic_and_deriv(p0, p1, p2, p3, t1); - if this_q1.length_squared() < DERIV_THRESH.powi(2) { - let (new_p1, new_q1) = eval_cubic_and_deriv(p0, p1, p2, p3, t1 - DERIV_EPS); - this_q1 = new_q1; - // Change just the derivative at the endpoint, but also move the point so it - // matches the derivative exactly if in the interior. - if t1 < 1. { - this_p1 = new_p1; - t1 -= DERIV_EPS; - } + let mut t1 = t0 + dt; + let this_p0 = last_p; + let this_q0 = last_q; + let (mut this_p1, mut this_q1) = eval_cubic_and_deriv(p0, p1, p2, p3, t1); + if this_q1.length_squared() < DERIV_THRESH.powi(2) { + let (new_p1, new_q1) = eval_cubic_and_deriv(p0, p1, p2, p3, t1 - DERIV_EPS); + this_q1 = new_q1; + // Change just the derivative at the endpoint, but also move the point so it + // matches the derivative exactly if in the interior. + if t1 < 1. { + this_p1 = new_p1; + t1 -= DERIV_EPS; } - let actual_dt = t1 - last_t; - let chord_len = (this_p1 - this_p0).length(); - // Subdivide the loop case when the chord is short, but don't subdivide when it is - // simply a very short segment. - if chord_len >= TANGENT_THRESH - || (this_q0.length_squared() * actual_dt * actual_dt < DERIV_THRESH.powi(2) - && this_q1.length_squared() * actual_dt * actual_dt < DERIV_THRESH.powi(2)) - { - let cubic_params = - CubicParams::from_points_derivs(this_p0, this_p1, this_q0, this_q1, actual_dt); - let est_err = cubic_params.est_euler_err(); - let err = est_err * chord_len; - log!("@@@ loop2: sub:{:?}, {:?} t0: {t0}, t1: {t1}, dt: {dt}, est_err: {est_err}, err: {err}", subcubic, cubic_params); - if err * scale <= tol || dt <= SUBDIV_LIMIT { - log!("@@@ error within tolerance"); - t0_u += 1; - let shift = t0_u.trailing_zeros(); - t0_u >>= shift; - dt *= (1 << shift) as f32; - let euler_params = EulerParams::from_angles(cubic_params.th0, cubic_params.th1); - let es = EulerSeg::from_params(this_p0, this_p1, euler_params); + } + let actual_dt = t1 - last_t; + let cubic_params = + CubicParams::from_points_derivs(this_p0, this_p1, this_q0, this_q1, actual_dt); + log!("@@@ loop2: sub:{:?}, {:?} t0: {t0}, t1: {t1}, dt: {dt}, est_err: {est_err}, err: {err}", subcubic, cubic_params); + if cubic_params.err * scale <= tol || dt <= SUBDIV_LIMIT { + log!("@@@ error within tolerance"); + let euler_params = EulerParams::from_angles(cubic_params.th0, cubic_params.th1); + let es = EulerSeg::from_params(this_p0, this_p1, euler_params); - let (k0, k1) = (es.params.k0 - 0.5 * es.params.k1, es.params.k1); + let (k0, k1) = (es.params.k0 - 0.5 * es.params.k1, es.params.k1); - // compute forward integral to determine number of subdivisions - let dist_scaled = offset * es.params.ch / chord_len; - // The number of subdivisions for curvature = 1 - let scale_multiplier = - 0.5 * FRAC_1_SQRT_2 * (scale * chord_len / (es.params.ch * tol)).sqrt(); - // TODO: tune these thresholds - const K1_THRESH: f32 = 1e-3; - const DIST_THRESH: f32 = 1e-3; - let mut a = 0.0; - let mut b = 0.0; - let mut integral = 0.0; - let mut int0 = 0.0; - let (n_frac, robust) = if k1.abs() < K1_THRESH { - let k = k0 + 0.5 * k1; - let n_frac = (k * (k * dist_scaled + 1.0)).abs().sqrt(); - (n_frac, EspcRobust::LowK1) - } else if dist_scaled.abs() < DIST_THRESH { - let f = |x: f32| x * x.abs().sqrt(); - a = k1; - b = k0; - int0 = f(b); - let int1 = f(a + b); - integral = int1 - int0; - //println!("int0={int0}, int1={int1} a={a} b={b}"); - let n_frac = (2. / 3.) * integral / a; - (n_frac, EspcRobust::LowDist) - } else { - a = -2.0 * dist_scaled * k1; - b = -1.0 - 2.0 * dist_scaled * k0; - int0 = espc_int_approx(b); - let int1 = espc_int_approx(a + b); - integral = int1 - int0; - let k_peak = k0 - k1 * b / a; - let integrand_peak = (k_peak * (k_peak * dist_scaled + 1.0)).abs().sqrt(); - let scaled_int = integral * integrand_peak / a; - let n_frac = scaled_int; - (n_frac, EspcRobust::Normal) - }; - let n = (n_frac * scale_multiplier).ceil().max(1.0); + // compute forward integral to determine number of subdivisions + let offset_chord_normalized = offset / cubic_params.chord_len; + let dist_scaled = offset_chord_normalized * es.params.ch; + // The number of subdivisions for curvature = 1 + let scale_multiplier = 0.5 + * FRAC_1_SQRT_2 + * (scale * cubic_params.chord_len / (es.params.ch * tol)).sqrt(); + // TODO: tune these thresholds + const K1_THRESH: f32 = 1e-3; + const DIST_THRESH: f32 = 1e-3; + let mut a = 0.0; + let mut b = 0.0; + let mut integral = 0.0; + let mut int0 = 0.0; + let (n_frac, robust) = if k1.abs() < K1_THRESH { + let k = k0 + 0.5 * k1; + let n_frac = (k * (k * dist_scaled + 1.0)).abs().sqrt(); + (n_frac, EspcRobust::LowK1) + } else if dist_scaled.abs() < DIST_THRESH { + let f = |x: f32| x * x.abs().sqrt(); + a = k1; + b = k0; + int0 = f(b); + let int1 = f(a + b); + integral = int1 - int0; + //println!("int0={int0}, int1={int1} a={a} b={b}"); + let n_frac = (2. / 3.) * integral / a; + (n_frac, EspcRobust::LowDist) + } else { + a = -2.0 * dist_scaled * k1; + b = -1.0 - 2.0 * dist_scaled * k0; + int0 = espc_int_approx(b); + let int1 = espc_int_approx(a + b); + integral = int1 - int0; + let k_peak = k0 - k1 * b / a; + let integrand_peak = (k_peak * (k_peak * dist_scaled + 1.0)).abs().sqrt(); + let scaled_int = integral * integrand_peak / a; + let n_frac = scaled_int; + (n_frac, EspcRobust::Normal) + }; + let n = (n_frac * scale_multiplier).ceil().max(1.0); - // Flatten line segments - log!("@@@ loop2: lines: {n}"); - // TODO: make all computation above robust and uncomment this assertion - //assert!(!n.is_nan()); - if n.is_nan() { - // Skip the segment if `n` is NaN. This is for debugging purposes only - log!("@@@ NaN: parameters:\n es: {:#?}\n k0: {k0}, k1: {k1}\n dist_scaled: {dist_scaled}\n es_scale: {es_scale}\n a: {a}\n b: {b}\n int0: {int0}, int1: {int1}, integral: {integral}\n k_peak: {k_peak}\n integrand_peak: {integrand_peak}\n scaled_int: {scaled_int}\n n_frac: {n_frac}", es); - } else { - for i in 0..n as usize { - let lp1 = if i == n as usize - 1 && t1 == 1.0 { - t_end - } else { - let t = (i + 1) as f32 / n; - let s = match robust { - EspcRobust::LowK1 => t, - // Note opportunities to minimize divergence - EspcRobust::LowDist => { - let c = (integral * t + int0).cbrt(); - let inv = c * c.abs(); - (inv - b) / a - } - EspcRobust::Normal => { - let inv = espc_int_inv_approx(integral * t + int0); - (inv - b) / a - } - }; - es.eval_with_offset(s, offset) - }; - let l0 = if offset >= 0. { lp0 } else { lp1 }; - let l1 = if offset >= 0. { lp1 } else { lp0 }; - output_line_with_transform( - path_ix, l0, l1, &transform, line_ix, lines, bbox, - ); - lp0 = lp1; + // Flatten line segments + log!("@@@ loop2: lines: {n}"); + assert!(!n.is_nan()); + for i in 0..n as usize { + let lp1 = if i == n as usize - 1 && t1 == 1.0 { + t_end + } else { + let t = (i + 1) as f32 / n; + let s = match robust { + EspcRobust::LowK1 => t, + // Note opportunities to minimize divergence + EspcRobust::LowDist => { + let c = (integral * t + int0).cbrt(); + let inv = c * c.abs(); + (inv - b) / a } - } - last_p = this_p1; - last_q = this_q1; - last_t = t1; - break; - } + EspcRobust::Normal => { + let inv = espc_int_inv_approx(integral * t + int0); + (inv - b) / a + } + }; + es.eval_with_offset(s, offset_chord_normalized) + }; + let l0 = if offset >= 0. { lp0 } else { lp1 }; + let l1 = if offset >= 0. { lp1 } else { lp0 }; + output_line_with_transform(path_ix, l0, l1, &transform, line_ix, lines, bbox); + lp0 = lp1; } + last_p = this_p1; + last_q = this_q1; + last_t = t1; + t0_u += 1; + let shift = t0_u.trailing_zeros(); + t0_u >>= shift; + dt *= (1 << shift) as f32; + } else { t0_u = t0_u.saturating_mul(2); dt *= 0.5; } @@ -661,9 +645,6 @@ const PATH_TAG_CUBICTO: u8 = 3; const PATH_TAG_F32: u8 = 8; -// Threshold for tangents to be considered near zero length -const TANGENT_THRESH: f32 = 1e-6; - fn flatten_main( n_wg: u32, config: &ConfigUniform,