fix(geom): use scale-aware circle-circle intersections (#384)

crates/geom/src/intersect.rs

@@ -53,19 +53,39 @@ pub fn line_circle(l: &Line, c: &Circle) -> Vec<(Vec2, f64)> {
 
 /// Intersections of two full circles.
 pub fn circle_circle(c1: &Circle, c2: &Circle) -> Vec<Vec2> {
-    let d = c1.center.dist(c2.center);
-    if d < EPS || d > c1.radius + c2.radius + 1e-9 || d < (c1.radius - c2.radius).abs() - 1e-9 {
+    let delta = c2.center - c1.center;
+    let d = delta.len();
+    let (r1, r2) = (c1.radius, c2.radius);
+    if !(d.is_finite() && r1.is_finite() && r2.is_finite()) || d == 0.0 || r1 < 0.0 || r2 < 0.0 {
+        return Vec::new();
+    }
+
+    // Work at the scale of the two circles: fixed absolute tolerances classify
+    // small secants as tangencies, and raw radius squares overflow on large drawings.
+    let scale = d.max(r1).max(r2);
+    let (ds, r1s, r2s) = (d / scale, r1 / scale, r2 / scale);
+    if ds > r1s + r2s || ds < (r1s - r2s).abs() {
+        return Vec::new();
+    }
+    let a = ((r1s - r2s) * (r1s + r2s) + ds * ds) / (2.0 * ds);
+    let h2 = (r1s - a) * (r1s + a);
+    if h2 < -1e-12 {
+        return Vec::new();
+    }
+    let dir = delta / d;
+    let p = c1.center + dir * (a * scale);
+    if !p.is_finite() {
         return Vec::new();
     }
-    let a = (c1.radius * c1.radius - c2.radius * c2.radius + d * d) / (2.0 * d);
-    let h2 = c1.radius * c1.radius - a * a;
-    let dir = (c2.center - c1.center) / d;
-    let p = c1.center + dir * a;
     if h2 <= 1e-12 {
         return vec![p];
     }
-    let h = h2.sqrt();
-    vec![p + dir.perp() * h, p - dir.perp() * h]
+    let step = dir.perp() * (h2.sqrt() * scale);
+    let (p1, p2) = (p + step, p - step);
+    if !p1.is_finite() || !p2.is_finite() {
+        return Vec::new();
+    }
+    vec![p1, p2]
 }
 
 fn on_arc(a: &Arc, p: Vec2) -> bool {

crates/geom/src/tests.rs

@@ -109,6 +109,25 @@ fn circle_circle_hits() {
     assert!(h.iter().all(|p| close(p.x, 4.0) && close(p.y.abs(), 3.0)));
 }
 
+#[test]
+fn tiny_circle_secants_and_separated_circles() {
+    let first = Circle::new(Vec2::ZERO, 1e-9);
+    // Both radii and center separation are below the old absolute tolerance.
+    let hits = circle_circle(&first, &Circle::new(Vec2::new(5e-10, 0.0), 1e-9));
+    assert_eq!(hits.len(), 2);
+    let y = (1e-18_f64 - 2.5e-10_f64.powi(2)).sqrt();
+    assert!(hits.iter().all(|p| (p.x - 2.5e-10).abs() < 1e-20 && (p.y.abs() - y).abs() < 1e-20), "{hits:?}");
+    // A half-nanometre gap is significant relative to these circles.
+    assert!(circle_circle(&first, &Circle::new(Vec2::new(2.5e-9, 0.0), 1e-9)).is_empty());
+}
+
+#[test]
+fn huge_circle_secants_remain_finite() {
+    let hits = circle_circle(&Circle::new(Vec2::ZERO, 1e160), &Circle::new(Vec2::new(1e160, 0.0), 1e160));
+    assert_eq!(hits.len(), 2);
+    assert!(hits.iter().all(|p| p.is_finite() && (p.x / 1e160 - 0.5).abs() < 1e-12 && (p.y.abs() / 1e160 - 3.0_f64.sqrt() / 2.0).abs() < 1e-12), "{hits:?}");
+}
+
 #[test]
 fn spline_interpolates_fit_points() {
     let fit = [Vec2::new(0.0, 0.0), Vec2::new(1.0, 2.0), Vec2::new(3.0, 1.0), Vec2::new(4.0, 4.0), Vec2::new(6.0, 0.0)];