linebender · jneem · Nov 27, 2025 · Aug 29, 2025
diff --git a/Cargo.lock b/Cargo.lock
diff --git a/Cargo.toml b/Cargo.toml
@@ -18,6 +18,7 @@ repository = "https://github.com/linebender/kurbo"
 
 [workspace.dependencies]
 kurbo = { version = "0.13.0", path = "kurbo", default-features = false }
+polycool = { version = "0.4.0", path = "polycool", default-features = false }
 
 [workspace.lints]
 rust.unsafe_code = "forbid"

diff --git a/kurbo/Cargo.toml b/kurbo/Cargo.toml
@@ -21,9 +21,9 @@ workspace = true
 [features]
 default = ["std"]
 # Get floating point functions from the standard library (likely using your target's libc)
-std = ["euclid?/std"]
+std = ["polycool/std", "euclid?/std"]
 # Use floating point implementations from libm
-libm = ["dep:libm", "euclid?/libm"]
+libm = ["dep:libm", "polycool/libm", "euclid?/libm"]
 # Enable `From`/`Into` conversion of Kurbo and [mint] types, enabling interoperability
 # with other graphics libraries
 mint = ["dep:mint"]
@@ -37,6 +37,7 @@ schemars = ["schemars/smallvec", "dep:schemars"]
 [dependencies]
 smallvec = "1.15.1"
 euclid = { version = "0.22", optional = true, default-features = false }
+polycool.workspace = true
 
 [dependencies.arrayvec]
 version = "0.7.6"
@@ -91,6 +92,10 @@ harness = false
 name = "quartic"
 harness = false
 
+[[bench]]
+name = "nearest"
+harness = false
+
 [[bench]]
 name = "rect_expand"
 harness = false
diff --git a/kurbo/benches/nearest.rs b/kurbo/benches/nearest.rs
@@ -0,0 +1,53 @@
+// Copyright 2022 the Kurbo Authors
+// SPDX-License-Identifier: Apache-2.0 OR MIT
+
+//! Benchmarks of the quartic equation solver.
+
+#![allow(missing_docs, reason = "criterion emits undocumented functions")]
+
+use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion};
+use std::hint::black_box;
+
+use kurbo::{CubicBez, ParamCurveNearest as _, Point, QuadBez};
+
+fn bench_nearest_quadratic(cc: &mut Criterion) {
+    let q = QuadBez {
+        p0: (-1.0, -1.0).into(),
+        p1: (0.0, 2.0).into(),
+        p2: (1.0, -1.0).into(),
+    };
+    let p = Point::new(0.0, 0.0);
+
+    for acc in [1e-3, 1e-6, 1e-12] {
+        cc.bench_with_input(
+            BenchmarkId::new("quadratic nearest point", acc),
+            &acc,
+            |bb, acc| {
+                bb.iter(|| black_box(q).nearest(black_box(p), *acc));
+            },
+        );
+    }
+}
+
+fn bench_nearest_cubic(cc: &mut Criterion) {
+    let c = CubicBez {
+        p0: (-1.0, -1.0).into(),
+        p1: (0.0, 2.0).into(),
+        p2: (1.0, -1.0).into(),
+        p3: (2.0, 2.0).into(),
+    };
+    let p = Point::new(0.0, 0.0);
+
+    for acc in [1e-3, 1e-6, 1e-12] {
+        cc.bench_with_input(
+            BenchmarkId::new("cubic nearest point", acc),
+            &acc,
+            |bb, acc| {
+                bb.iter(|| black_box(c).nearest(black_box(p), *acc));
+            },
+        );
+    }
+}
+
+criterion_group!(benches, bench_nearest_quadratic, bench_nearest_cubic);
+criterion_main!(benches);
diff --git a/kurbo/src/cubicbez.rs b/kurbo/src/cubicbez.rs
@@ -654,23 +654,52 @@ impl ParamCurveArea for CubicBez {
 }
 
 impl ParamCurveNearest for CubicBez {
-    /// Find the nearest point, using subdivision.
+    /// Find the nearest point using a quintic solver.
+    ///
+    /// The polynomial `|self - p|^2` has degree 6, so we find its critical
+    /// points and evaluate them all to find the best one.
     fn nearest(&self, p: Point, accuracy: f64) -> Nearest {
-        let mut best_r = None;
-        let mut best_t = 0.0;
-        for (t0, t1, q) in self.to_quads(accuracy) {
-            let nearest = q.nearest(p, accuracy);
-            if best_r
-                .map(|best_r| nearest.distance_sq < best_r)
-                .unwrap_or(true)
-            {
-                best_t = t0 + nearest.t * (t1 - t0);
-                best_r = Some(nearest.distance_sq);
+        fn eval_t(p: Point, t_best: &mut f64, r_best: &mut Option<f64>, t: f64, p0: Point) {
+            let r = (p0 - p).hypot2();
+            if r_best.map(|r_best| r < r_best).unwrap_or(true) {
+                *r_best = Some(r);
+                *t_best = t;
             }
         }
+
+        let mut r_best = None;
+        let mut t_best = 0.0;
+
+        // Reparameterize `self - p` as q0 + q1 t + q2 t^2 + q3 t^3.
+        let q0 = self.p0 - p;
+        let q1 = 3.0 * (self.p1 - self.p0);
+        let q2 = 3.0 * (self.p0.to_vec2() - 2.0 * self.p1.to_vec2() + self.p2.to_vec2());
+        let q3 = -self.p0.to_vec2() + 3.0 * self.p1.to_vec2() - 3.0 * self.p2.to_vec2()
+            + self.p3.to_vec2();
+
+        // Coefficients of the degree-5 polynomial (self - p) \cdot tangent.
+        let c0 = q0.dot(q1);
+        let c1 = q1.hypot2() + 2.0 * q2.dot(q0);
+        let c2 = 3.0 * (q2.dot(q1) + q3.dot(q0));
+        let c3 = 4.0 * q3.dot(q1) + 2.0 * q2.hypot2();
+        let c4 = 5.0 * q3.dot(q2);
+        let c5 = 3.0 * q3.hypot2();
+
+        let roots = polycool::Poly::new([c0, c1, c2, c3, c4, c5]).roots_between(0.0, 1.0, accuracy);
+
+        for &t in &roots {
+            eval_t(p, &mut t_best, &mut r_best, t, self.eval(t));
+        }
+
+        // If we found all 5 critical points, we can skip evaluating the endpoints.
+        if roots.len() != 5 {
+            eval_t(p, &mut t_best, &mut r_best, 0.0, self.p0);
+            eval_t(p, &mut t_best, &mut r_best, 1.0, self.p3);
+        }
+
         Nearest {
-            t: best_t,
-            distance_sq: best_r.unwrap(),
+            t: t_best,
+            distance_sq: r_best.unwrap(),
         }
     }
 }
@@ -910,6 +939,22 @@ mod tests {
         verify(c.nearest((-0.1, 0.0).into(), 1e-6), 0.0);
         let a = Affine::rotate(0.5);
         verify((a * c).nearest(a * Point::new(0.1, 0.001), 1e-6), 0.1);
+
+        // Here's a case that tripped up the old solver because the start is close to
+        // degenerate and the end is actually degenerate; see #446.
+        let curve = CubicBez::new(
+            (461.0, 123.0),
+            (460.99999999999994, 123.00000000000004),
+            (111.0, 319.0),
+            (111.0, 319.0),
+        );
+        let p = Point::new(282.0379003395483, 223.21877580985594);
+        let eps = 0.0005;
+        let nearest = curve.nearest(p, eps);
+        let q = curve.eval(nearest.t);
+        let r = curve.eval(0.5075474297354187);
+
+        assert!((q - p).hypot() <= (r - p).hypot() + eps);
     }
 
     // ensure to_quads returns something given collinear points