From 6f8af5d6420686832dacb5ff7c064f68ed4bc50d Mon Sep 17 00:00:00 2001
From: Tom Churchman <thomas@churchman.nl>
Date: Sun, 14 Sep 2025 23:51:11 +0200
Subject: [PATCH] Tighten bounds for direct curve-to-line-segment flattening

---
 .../vello_common/src/flatten_simd.rs          | 26 ++++++++++++++++---
 1 file changed, 23 insertions(+), 3 deletions(-)

diff --git a/sparse_strips/vello_common/src/flatten_simd.rs b/sparse_strips/vello_common/src/flatten_simd.rs
index 4fb5ef323..3d5ac83e0 100644
--- a/sparse_strips/vello_common/src/flatten_simd.rs
+++ b/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);