Change Line::nearest implementation to be branchless (#505)

This is a branchless calculation of a point-onto-line-segment
projection. Also marks the method `#[inline]` and adds some testing.

I started looking at this because of the math we introduced in
linebender/vello#1214. This might allow us to
introduce some tighter bounds there without introducing additional
branching.

This should be useful here regardless of whether the above happens.

---------

Co-authored-by: jneem <[email protected]>

Author: tomcur

CHANGELOG.md

@@ -22,6 +22,7 @@ This release has an [MSRV][] of 1.82.
 ## Changed
 
 - Speed up methods like `Ellipse::radii` by reworking the singular value decomposition expression. ([#499][] by [@tomcur][])
+- The `Line::nearest` method to calculate the projection of a point onto a line segment has been made more performant. ([#505][] by [@tomcur][])
 
 ## [0.12.0] (2025-09-04)
 

kurbo/src/line.rs

@@ -158,19 +158,29 @@ impl ParamCurveArea for Line {
 }
 
 impl ParamCurveNearest for Line {
+    #[inline]
     fn nearest(&self, p: Point, _accuracy: f64) -> Nearest {
         let d = self.p1 - self.p0;
-        let dotp = d.dot(p - self.p0);
-        let d_squared = d.dot(d);
-        let (t, distance_sq) = if dotp <= 0.0 {
-            (0.0, (p - self.p0).hypot2())
-        } else if dotp >= d_squared {
-            (1.0, (p - self.p1).hypot2())
-        } else {
-            let t = dotp / d_squared;
-            let dist = (p - self.eval(t)).hypot2();
-            (t, dist)
-        };
+        let v = p - self.p0;
+
+        // Calculate projection parameter `t` of the point onto the line segment s(t), with
+        // s(t) = (1-t) * p0 + t * p1.
+        //
+        // Note when the segment has 0 length, this will be positive or negative infinity or NaN;
+        // see the clamping below.
+        let t = d.dot(v) / d.hypot2();
+
+        // Clamp the parameter to be on the line segment. This clamps negative infinity and NaN to
+        // `0.`, and positive infinity to `1.`.
+        #[expect(
+            clippy::manual_clamp,
+            reason = "`clamp` uses slightly more instructions than chained `max` and `min` on x86 and aarch64"
+        )]
+        let t = { t.max(0.).min(1.) };
+
+        // Calculate ||p - s(t)||^2.
+        let distance_sq = (v - t * d).hypot2();
+
         Nearest { distance_sq, t }
     }
 }
@@ -394,4 +404,39 @@ mod tests {
         })
         .is_finite());
     }
+
+    #[test]
+    fn line_nearest() {
+        use crate::{ParamCurve, ParamCurveNearest};
+
+        const EPSILON: f64 = 1e-9;
+
+        let line = Line::new((-4., 0.), (2., 1.));
+
+        // Projects onto the line segment end point.
+        let point = Point::new(4., 0.);
+        let nearest = line.nearest(point, 0.);
+        assert_eq!(nearest.t, 1.);
+        assert!((nearest.distance_sq - line.p1.distance_squared(point)).abs() < EPSILON);
+
+        // Projects onto the line segment start point.
+        let point = Point::new(0., -50.);
+        let nearest = line.nearest(point, 0.);
+        assert_eq!(nearest.t, 0.);
+        assert!((nearest.distance_sq - line.p0.distance_squared(point)).abs() < EPSILON);
+
+        // Projects onto the line segment proper (not just onto one of its extrema).
+        let point = Point::new(-1., 0.5);
+        let nearest = line.nearest(point, 0.);
+        assert!(nearest.t > 0. && nearest.t < 1.);
+        // Ensure evaluating and calculating distance manually has the same result.
+        assert!(
+            (line.eval(nearest.t).distance_squared(point) - nearest.distance_sq).abs() < EPSILON
+        );
+
+        // Test minimality while avoiding reimplementing projection in this test by checking that
+        // moving to a slightly different point on the segment increases the distance.
+        assert!(line.eval(nearest.t * 0.95).distance_squared(point) > nearest.distance_sq);
+        assert!(line.eval(nearest.t * 1.05).distance_squared(point) > nearest.distance_sq);
+    }
 }