Tighten bounds for direct curve-to-line-segment flattening

sparse_strips/vello_common/src/flatten_simd.rs

@@ -8,7 +8,7 @@
 use crate::flatten::TOL_2;
 #[cfg(not(feature = "std"))]
 use crate::kurbo::common::FloatFuncs as _;
-use crate::kurbo::{CubicBez, ParamCurve, PathEl, Point, QuadBez};
+use crate::kurbo::{CubicBez, Line, ParamCurve, PathEl, Point, QuadBez};
 use alloc::vec;
 use alloc::vec::Vec;
 use bytemuck::{Pod, Zeroable};
@@ -37,6 +37,24 @@ pub(crate) fn flatten<S: Simd>(
     let sqrt_tol = tolerance.sqrt();
     let mut last_pt = None;
 
+    fn dist2(line: crate::kurbo::Line, point: crate::kurbo::Point) -> f64 {
+        let d = line.p1 - line.p0;
+        let v = point - line.p0;
+
+        // Calculate projection parameter `t` of the point onto s(t), with s(t) the line segment
+        // such that s(t) = (1-t) * p0 + t * p1.
+        //
+        // Note this will be inf when the segment has 0 length; see the clamping below.
+        let t = d.dot(v) / d.hypot2();
+
+        // Clamp the parameter to be on the line segment. This results in `t==0` if `t==inf` above.
+        let t = t.max(0.).min(1.);
+
+        // Calculate ||p - s(t)||^2.
+        let distance_sq = (v - t * d).hypot2();
+        distance_sq
+    }
+
     for el in path {
         match el {
             PathEl::MoveTo(p) => {
@@ -67,7 +85,8 @@ pub(crate) fn flatten<S: Simd>(
                     //
                     // The following takes the square to elide the square root of the Euclidean
                     // distance.
-                    if f64::max((p1 - p0).hypot2(), (p1 - p2).hypot2()) <= 4. * TOL_2 {
+                    let line = Line::new(p0, p2);
+                    if dist2(line, p1) <= 4. * TOL_2 {
                         callback.callback(PathEl::LineTo(p2));
                     } else {
                         let q = QuadBez::new(p0, p1, p2);
@@ -107,7 +126,8 @@ pub(crate) fn flatten<S: Simd>(
                     //
                     // The following takes the square to elide the square root of the Euclidean
                     // distance.
-                    if f64::max((p0 - p1).hypot2(), (p3 - p2).hypot2()) <= 16. / 9. * TOL_2 {
+                    let line = Line::new(p0, p3);
+                    if f64::max(dist2(line, p1), dist2(line, p2)) <= 16. / 9. * TOL_2 {
                         callback.callback(PathEl::LineTo(p3));
                     } else {
                         let c = CubicBez::new(p0, p1, p2, p3);