Created
November 14, 2025 05:04
-
-
Save raphlinus/0c024e27289de3052749c40543e72f38 to your computer and use it in GitHub Desktop.
Compute curvature bounds for cubic Bézier using interval arithmetic
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| use kurbo::{CubicBez, ParamCurveDeriv}; | |
| /// Compute bounds on curvature | |
| pub fn kbound(c: CubicBez) -> (f64, f64) { | |
| let q = c.deriv(); | |
| let p1xp0 = q.p1.to_vec2().cross(q.p0.to_vec2()); | |
| let p2xp0 = q.p2.to_vec2().cross(q.p0.to_vec2()); | |
| let p2xp1 = q.p2.to_vec2().cross(q.p1.to_vec2()); | |
| let c0 = 2. * p1xp0; | |
| let c1 = 2. * (p2xp0 - 2.0 * p1xp0); | |
| let c2 = 2. * (p2xp1 - p2xp0 + p1xp0); | |
| let (num_min, num_max) = quadratic_min_max(c0, c1, c2); | |
| let norm_chord = (c.p3 - c.p0).normalize(); | |
| let r0 = q.p0.to_vec2().dot(norm_chord); | |
| let r1 = q.p1.to_vec2().dot(norm_chord); | |
| let r2 = q.p2.to_vec2().dot(norm_chord); | |
| let d0 = r0; | |
| let d1 = 2. * (r1 - r0); | |
| let d2 = (r2 - r1) - (r1 - r0); | |
| let (deriv_min, _deriv_max) = quadratic_min_max(d0, d1, d2); | |
| if deriv_min <= 0.0 { | |
| return (-f64::INFINITY, f64::INFINITY); | |
| } | |
| let scale = 1.0 / (deriv_min * deriv_min * deriv_min); | |
| (num_min.min(0.) * scale, num_max.max(0.) * scale) | |
| } | |
| fn quadratic_min_max(c0: f64, c1: f64, c2: f64) -> (f64, f64) { | |
| let y1 = c0 + c1 + c2; | |
| let mut min = c0.min(y1); | |
| let mut max = c0.max(y1); | |
| if c2 * c1 < 0. && c2 * (c1 + 2.0 * c2) > 0. { | |
| let x = -c1 / (2.0 * c2); | |
| let y = c0 + x * (c1 + x * c2); | |
| min = min.min(y); | |
| max = max.max(y); | |
| } | |
| (min, max) | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment