dotloom_geometry/
intersect.rs

1//! Curve/curve intersections.
2//!
3//! * Segment/segment decisions use exact orientation predicates; the reported point
4//!   is computed in floating point.
5//! * Line/circle and circle/circle use numerically stable closed forms; near-tangent
6//!   cases within the model tolerance report a single tangent point.
7//! * Arcs reuse the circle results filtered by their angular span.
8//! * Bézier intersections flatten to a fine polyline, intersect chords, then refine
9//!   the parameters with Newton iterations on the exact curves.
10
11use core::f64::consts::TAU;
12
13use serde::{Deserialize, Serialize};
14
15use crate::{Arc, Circle, CubicBez, Curve, ModelTolerance, Orientation, Point, Segment, Vector, orientation};
16
17/// One intersection point with the parameters on both curves.
18#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
19pub struct Intersection {
20    /// Intersection point.
21    pub point: Point,
22    /// Parameter on the first curve.
23    pub t_a: f64,
24    /// Parameter on the second curve.
25    pub t_b: f64,
26}
27
28/// Result of intersecting two curves.
29#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
30pub struct Intersections {
31    /// Isolated intersection points (for overlaps: the overlap endpoints).
32    pub points: Vec<Intersection>,
33    /// The curves share a stretch of positive length (collinear or co-circular overlap).
34    pub overlap: bool,
35}
36
37impl Intersections {
38    fn none() -> Self {
39        Self::default()
40    }
41
42    fn swapped(mut self) -> Self {
43        for p in &mut self.points {
44            core::mem::swap(&mut p.t_a, &mut p.t_b);
45        }
46        self
47    }
48}
49
50/// Intersect two curves.
51#[must_use]
52pub fn intersect(a: &Curve, b: &Curve, tol: ModelTolerance) -> Intersections {
53    if !a.bbox().inflate(tol.abs).intersects(b.bbox().inflate(tol.abs)) {
54        return Intersections::none();
55    }
56    let mut r = match (*a, *b) {
57        (Curve::Line(s1), Curve::Line(s2)) => segment_segment(s1, s2, tol),
58        (Curve::Line(s), Curve::Circle(c)) => line_circle(s, c, tol, Span::Full),
59        (Curve::Circle(c), Curve::Line(s)) => line_circle(s, c, tol, Span::Full).swapped(),
60        (Curve::Line(s), Curve::Arc(a2)) => line_circle(s, circle_of(a2), tol, Span::Arc(a2)),
61        (Curve::Arc(a1), Curve::Line(s)) => line_circle(s, circle_of(a1), tol, Span::Arc(a1)).swapped(),
62        (Curve::Circle(c1), Curve::Circle(c2)) => circle_circle(c1, Span::Full, c2, Span::Full, tol),
63        (Curve::Circle(c1), Curve::Arc(a2)) => circle_circle(c1, Span::Full, circle_of(a2), Span::Arc(a2), tol),
64        (Curve::Arc(a1), Curve::Circle(c2)) => circle_circle(circle_of(a1), Span::Arc(a1), c2, Span::Full, tol),
65        (Curve::Arc(a1), Curve::Arc(a2)) => {
66            circle_circle(circle_of(a1), Span::Arc(a1), circle_of(a2), Span::Arc(a2), tol)
67        }
68        (Curve::Cubic(c), other) => cubic_curve(c, &other, tol),
69        (other, Curve::Cubic(c)) => cubic_curve(c, &other, tol).swapped(),
70    };
71    dedup(&mut r.points, tol);
72    r
73}
74
75#[derive(Clone, Copy)]
76enum Span {
77    Full,
78    Arc(Arc),
79}
80
81impl Span {
82    /// Parameter for a point at absolute angle `ang` or `None` when outside.
83    fn param(self, ang: f64, radius: f64, tol: ModelTolerance) -> Option<f64> {
84        match self {
85            Self::Full => Some(crate::normalize_angle(ang) / TAU),
86            Self::Arc(a) => a.param_of_angle(ang, tol.at_scale(radius) / radius.max(f64::MIN_POSITIVE)),
87        }
88    }
89}
90
91fn circle_of(a: Arc) -> Circle {
92    Circle { center: a.center, radius: a.radius }
93}
94
95fn segment_segment(s1: Segment, s2: Segment, tol: ModelTolerance) -> Intersections {
96    let o1 = orientation(s1.a, s1.b, s2.a);
97    let o2 = orientation(s1.a, s1.b, s2.b);
98    let o3 = orientation(s2.a, s2.b, s1.a);
99    let o4 = orientation(s2.a, s2.b, s1.b);
100
101    let d1 = s1.vector();
102    let d2 = s2.vector();
103    let l1 = d1.length_sq();
104    let l2 = d2.length_sq();
105
106    if o1 == Orientation::Collinear && o2 == Orientation::Collinear {
107        // Collinear (or degenerate) segments: compute the overlap on s1's parameter.
108        if l1 == 0.0 && l2 == 0.0 {
109            return if s1.a.distance(s2.a) <= tol.abs { single(s1.a, 0.0, 0.0) } else { Intersections::none() };
110        }
111        if l1 == 0.0 {
112            let t = s2.line_param(s1.a);
113            return if (0.0..=1.0).contains(&t) { single(s1.a, 0.0, t) } else { Intersections::none() };
114        }
115        let ta = s1.line_param(s2.a);
116        let tb = s1.line_param(s2.b);
117        let lo = ta.min(tb).max(0.0);
118        let hi = ta.max(tb).min(1.0);
119        let eps = tol.at_scale(l1.sqrt()) / l1.sqrt();
120        if hi < lo - eps {
121            return Intersections::none();
122        }
123        let param_on_2 = |p: Point| s2.line_param(p).clamp(0.0, 1.0);
124        if hi - lo <= eps {
125            let p = s1.point_at(lo.clamp(0.0, 1.0));
126            return single(p, lo.clamp(0.0, 1.0), param_on_2(p));
127        }
128        let p_lo = s1.point_at(lo);
129        let p_hi = s1.point_at(hi);
130        return Intersections {
131            points: vec![
132                Intersection { point: p_lo, t_a: lo, t_b: param_on_2(p_lo) },
133                Intersection { point: p_hi, t_a: hi, t_b: param_on_2(p_hi) },
134            ],
135            overlap: true,
136        };
137    }
138
139    let crosses = o1 != o2 && o3 != o4;
140    if !crosses {
141        // Near misses within tolerance (endpoint touching a segment).
142        return near_touch(s1, s2, tol);
143    }
144    let denom = d1.cross(d2);
145    if denom == 0.0 {
146        return near_touch(s1, s2, tol);
147    }
148    let w = s2.a - s1.a;
149    let t = (w.cross(d2) / denom).clamp(0.0, 1.0);
150    let u = (w.cross(d1) / denom).clamp(0.0, 1.0);
151    // Average both evaluations for a symmetric result.
152    let p = s1.point_at(t).midpoint(s2.point_at(u));
153    single(p, t, u)
154}
155
156fn near_touch(s1: Segment, s2: Segment, tol: ModelTolerance) -> Intersections {
157    let mut out = Intersections::none();
158    for (p, on_first) in [(s1.a, true), (s1.b, true), (s2.a, false), (s2.b, false)] {
159        let (other, own) = if on_first { (s2, s1) } else { (s1, s2) };
160        let (t_other, q) = other.closest(p);
161        if q.distance(p) <= tol.at_scale(p.to_vector().length()) {
162            let t_own = own.line_param(p).clamp(0.0, 1.0);
163            let (t_a, t_b) = if on_first { (t_own, t_other) } else { (t_other, t_own) };
164            out.points.push(Intersection { point: p, t_a, t_b });
165        }
166    }
167    out
168}
169
170fn single(point: Point, t_a: f64, t_b: f64) -> Intersections {
171    Intersections { points: vec![Intersection { point, t_a, t_b }], overlap: false }
172}
173
174/// Intersections of the *infinite* line through `s` with a circle: `(line_t, angle)`.
175pub(crate) fn line_circle_raw(s: Segment, c: Circle, tol: ModelTolerance) -> Vec<(f64, f64)> {
176    let d = s.vector();
177    let dl = d.length();
178    if dl == 0.0 || !dl.is_finite() {
179        return Vec::new();
180    }
181    let u = d / dl;
182    // Foot of the perpendicular from the center.
183    let along = (c.center - s.a).dot(u);
184    let foot = s.a + u * along;
185    let dist = foot.distance(c.center);
186    let r = c.radius;
187    let ttol = tol.at_scale(r.max(dist));
188    if dist > r + ttol {
189        return Vec::new();
190    }
191    let h2 = r * r - dist * dist;
192    let mk = |off: f64| {
193        let p = foot + u * off;
194        ((along + off) / dl, (p - c.center).angle())
195    };
196    if h2 <= 0.0 || (r - dist).abs() <= ttol {
197        return vec![mk(0.0)];
198    }
199    let h = h2.sqrt();
200    vec![mk(-h), mk(h)]
201}
202
203fn line_circle(s: Segment, c: Circle, tol: ModelTolerance, span: Span) -> Intersections {
204    let len = s.length();
205    let eps = if len > 0.0 { tol.at_scale(len) / len } else { 0.0 };
206    let mut out = Intersections::none();
207    for (t, ang) in line_circle_raw(s, c, tol) {
208        if t < -eps || t > 1.0 + eps {
209            continue;
210        }
211        if let Some(tc) = span.param(ang, c.radius, tol) {
212            let t = t.clamp(0.0, 1.0);
213            out.points.push(Intersection { point: s.point_at(t), t_a: t, t_b: tc });
214        }
215    }
216    out
217}
218
219/// Intersection angles of two full circles: `(angle_on_1, angle_on_2)`.
220/// Returns `None` for coincident circles.
221pub(crate) fn circle_circle_raw(c1: Circle, c2: Circle, tol: ModelTolerance) -> Option<Vec<(f64, f64)>> {
222    let dv = c2.center - c1.center;
223    let d = dv.length();
224    let ttol = tol.at_scale(c1.radius.max(c2.radius).max(d));
225    if d <= ttol {
226        return if (c1.radius - c2.radius).abs() <= ttol { None } else { Some(Vec::new()) };
227    }
228    if d > c1.radius + c2.radius + ttol || d < (c1.radius - c2.radius).abs() - ttol {
229        return Some(Vec::new());
230    }
231    let u = dv / d;
232    let a = (c1.radius * c1.radius - c2.radius * c2.radius + d * d) / (2.0 * d);
233    let h2 = c1.radius * c1.radius - a * a;
234    let base = c1.center + u * a;
235    let ang = |p: Point| ((p - c1.center).angle(), (p - c2.center).angle());
236    let tangent = (d - (c1.radius + c2.radius)).abs() <= ttol || (d - (c1.radius - c2.radius).abs()).abs() <= ttol;
237    if h2 <= 0.0 || tangent {
238        return Some(vec![ang(base)]);
239    }
240    let h = h2.sqrt();
241    let n: Vector = u.perp();
242    Some(vec![ang(base + n * h), ang(base - n * h)])
243}
244
245fn circle_circle(c1: Circle, s1: Span, c2: Circle, s2: Span, tol: ModelTolerance) -> Intersections {
246    let Some(raw) = circle_circle_raw(c1, c2, tol) else {
247        return coincident_circles(c1, s1, s2, tol);
248    };
249    let mut out = Intersections::none();
250    for (a1, a2) in raw {
251        if let (Some(t1), Some(t2)) = (s1.param(a1, c1.radius, tol), s2.param(a2, c2.radius, tol)) {
252            out.points.push(Intersection { point: c1.point_at_angle(a1), t_a: t1, t_b: t2 });
253        }
254    }
255    out
256}
257
258fn coincident_circles(c: Circle, s1: Span, s2: Span, tol: ModelTolerance) -> Intersections {
259    // Co-circular curves: report overlap plus the span endpoints that lie on the other.
260    let mut out = Intersections { points: Vec::new(), overlap: false };
261    let ends = |s: Span| match s {
262        Span::Full => Vec::new(),
263        Span::Arc(a) => vec![a.start, a.end_angle()],
264    };
265    for ang in ends(s1) {
266        if let (Some(t1), Some(t2)) = (s1.param(ang, c.radius, tol), s2.param(ang, c.radius, tol)) {
267            out.points.push(Intersection { point: c.point_at_angle(ang), t_a: t1, t_b: t2 });
268        }
269    }
270    for ang in ends(s2) {
271        if let (Some(t1), Some(t2)) = (s1.param(ang, c.radius, tol), s2.param(ang, c.radius, tol)) {
272            out.points.push(Intersection { point: c.point_at_angle(ang), t_a: t1, t_b: t2 });
273        }
274    }
275    out.overlap = match (s1, s2) {
276        (Span::Full, _) | (_, Span::Full) => true,
277        (Span::Arc(a1), Span::Arc(a2)) => {
278            // Overlap if a midpoint of either lies on the other.
279            s2.param((a1.mid_point() - c.center).angle(), c.radius, tol).is_some()
280                || s1.param((a2.mid_point() - c.center).angle(), c.radius, tol).is_some()
281        }
282    };
283    if !out.overlap {
284        out.points.truncate(out.points.len().min(2));
285    }
286    out
287}
288
289/// Flattened cubic with parameter values, fine enough for intersection seeding.
290fn cubic_samples(c: CubicBez, tol: f64) -> Vec<(f64, Point)> {
291    let mut out = vec![(0.0, c.p0)];
292    sample_rec(c, 0.0, 1.0, tol, 0, &mut out);
293    out
294}
295
296fn sample_rec(c: CubicBez, t0: f64, t1: f64, tol: f64, depth: u32, out: &mut Vec<(f64, Point)>) {
297    if depth >= 14 || c.flatness() <= tol {
298        out.push((t1, c.p3));
299        return;
300    }
301    let (l, r) = c.split(0.5);
302    let tm = 0.5 * (t0 + t1);
303    sample_rec(l, t0, tm, tol, depth + 1, out);
304    sample_rec(r, tm, t1, tol, depth + 1, out);
305}
306
307fn cubic_curve(c: CubicBez, other: &Curve, tol: ModelTolerance) -> Intersections {
308    let size = c.hull_bbox().union(other.bbox()).size().length().max(1e-300);
309    let flat_tol = size * 1e-4;
310    let a_samples = cubic_samples(c, flat_tol);
311    let mut out = Intersections::none();
312    match *other {
313        Curve::Cubic(d) => {
314            let b_samples = cubic_samples(d, flat_tol);
315            for wa in a_samples.windows(2) {
316                let sa = Segment::new(wa[0].1, wa[1].1);
317                for wb in b_samples.windows(2) {
318                    let sb = Segment::new(wb[0].1, wb[1].1);
319                    let r = segment_segment(sa, sb, ModelTolerance { abs: flat_tol, rel: 0.0 });
320                    for p in r.points {
321                        let s0 = wa[0].0 + (wa[1].0 - wa[0].0) * p.t_a;
322                        let t0 = wb[0].0 + (wb[1].0 - wb[0].0) * p.t_b;
323                        if let Some((s, t)) = refine_cubic_cubic(c, d, s0, t0, tol) {
324                            out.points.push(Intersection { point: c.eval(s).midpoint(d.eval(t)), t_a: s, t_b: t });
325                        }
326                    }
327                }
328            }
329        }
330        _ => {
331            for wa in a_samples.windows(2) {
332                let sa = Curve::Line(Segment::new(wa[0].1, wa[1].1));
333                let r = intersect(&sa, other, ModelTolerance { abs: flat_tol, rel: 0.0 });
334                for p in r.points {
335                    let s0 = wa[0].0 + (wa[1].0 - wa[0].0) * p.t_a;
336                    if let Some((s, q, t)) = refine_cubic_curve(c, other, s0, tol) {
337                        out.points.push(Intersection { point: q, t_a: s, t_b: t });
338                    }
339                }
340            }
341        }
342    }
343    out
344}
345
346/// Newton on `dist(B(s), other)` using the closest-point projection onto `other`.
347fn refine_cubic_curve(c: CubicBez, other: &Curve, s0: f64, tol: ModelTolerance) -> Option<(f64, Point, f64)> {
348    let mut s = s0;
349    for _ in 0..30 {
350        let p = c.eval(s);
351        let (_, q) = other.closest(p);
352        let (t_on, _) = other.closest(q);
353        let n = match other.tangent_at(t_on) {
354            Some(tan) => tan.perp(),
355            None => break,
356        };
357        let f = (p - q).dot(n);
358        let fp = c.deriv(s).dot(n);
359        if fp.abs() < 1e-300 {
360            break;
361        }
362        let ns = (s - f / fp).clamp(0.0, 1.0);
363        if (ns - s).abs() < 1e-15 {
364            s = ns;
365            break;
366        }
367        s = ns;
368    }
369    let p = c.eval(s);
370    let (t, q) = other.closest(p);
371    (p.distance(q) <= tol.at_scale(p.to_vector().length()) * 10.0 + 1e-9).then_some((s, p.midpoint(q), t))
372}
373
374fn refine_cubic_cubic(a: CubicBez, b: CubicBez, s0: f64, t0: f64, tol: ModelTolerance) -> Option<(f64, f64)> {
375    let (mut s, mut t) = (s0, t0);
376    for _ in 0..30 {
377        let f = a.eval(s) - b.eval(t);
378        let da = a.deriv(s);
379        let db = b.deriv(t);
380        // Solve [da, -db] [ds, dt]^T = -f
381        let det = da.x * (-db.y) - (-db.x) * da.y;
382        if det.abs() < 1e-300 {
383            break;
384        }
385        let ds = (-f.x * (-db.y) - (-db.x) * (-f.y)) / det;
386        let dt = (da.x * (-f.y) - (-f.x) * da.y) / det;
387        s = (s + ds).clamp(0.0, 1.0);
388        t = (t + dt).clamp(0.0, 1.0);
389        if ds.abs() < 1e-15 && dt.abs() < 1e-15 {
390            break;
391        }
392    }
393    let p = a.eval(s);
394    (p.distance(b.eval(t)) <= tol.at_scale(p.to_vector().length()) * 10.0 + 1e-9).then_some((s, t))
395}
396
397fn dedup(points: &mut Vec<Intersection>, tol: ModelTolerance) {
398    let mut out: Vec<Intersection> = Vec::with_capacity(points.len());
399    for p in points.drain(..) {
400        let dup = out.iter().any(|q| q.point.distance(p.point) <= tol.at_scale(p.point.to_vector().length()) * 10.0);
401        if !dup {
402            out.push(p);
403        }
404    }
405    out.sort_by(|a, b| a.t_a.total_cmp(&b.t_a));
406    *points = out;
407}
408
409#[cfg(test)]
410mod tests {
411    use super::*;
412    use crate::Point;
413
414    fn tol() -> ModelTolerance {
415        ModelTolerance::DEFAULT
416    }
417
418    fn line(ax: f64, ay: f64, bx: f64, by: f64) -> Curve {
419        Curve::Line(Segment::new(Point::new(ax, ay), Point::new(bx, by)))
420    }
421
422    #[test]
423    fn crossing_segments() {
424        let r = intersect(&line(0.0, 0.0, 2.0, 2.0), &line(0.0, 2.0, 2.0, 0.0), tol());
425        assert_eq!(r.points.len(), 1);
426        assert!(r.points[0].point.distance(Point::new(1.0, 1.0)) < 1e-12);
427        assert!((r.points[0].t_a - 0.5).abs() < 1e-12);
428    }
429
430    #[test]
431    fn parallel_segments_do_not_intersect() {
432        let r = intersect(&line(0.0, 0.0, 2.0, 0.0), &line(0.0, 1.0, 2.0, 1.0), tol());
433        assert!(r.points.is_empty() && !r.overlap);
434    }
435
436    #[test]
437    fn collinear_overlap() {
438        let r = intersect(&line(0.0, 0.0, 4.0, 0.0), &line(2.0, 0.0, 6.0, 0.0), tol());
439        assert!(r.overlap);
440        assert_eq!(r.points.len(), 2);
441        assert!(r.points[0].point.distance(Point::new(2.0, 0.0)) < 1e-12);
442        assert!(r.points[1].point.distance(Point::new(4.0, 0.0)) < 1e-12);
443    }
444
445    #[test]
446    fn touching_endpoint() {
447        let r = intersect(&line(0.0, 0.0, 2.0, 0.0), &line(2.0, 0.0, 2.0, 5.0), tol());
448        assert_eq!(r.points.len(), 1);
449        assert!(r.points[0].point.distance(Point::new(2.0, 0.0)) < 1e-12);
450    }
451
452    #[test]
453    fn line_circle_two_and_tangent() {
454        let c = Curve::Circle(Circle::new(Point::ORIGIN, 1.0).unwrap());
455        let r = intersect(&line(-2.0, 0.0, 2.0, 0.0), &c, tol());
456        assert_eq!(r.points.len(), 2);
457        let t = intersect(&line(-2.0, 1.0, 2.0, 1.0), &c, tol());
458        assert_eq!(t.points.len(), 1);
459        assert!(t.points[0].point.distance(Point::new(0.0, 1.0)) < 1e-9);
460    }
461
462    #[test]
463    fn circle_circle_cases() {
464        let a = Curve::Circle(Circle::new(Point::ORIGIN, 1.0).unwrap());
465        let b = Curve::Circle(Circle::new(Point::new(1.0, 0.0), 1.0).unwrap());
466        let r = intersect(&a, &b, tol());
467        assert_eq!(r.points.len(), 2);
468        for p in &r.points {
469            assert!((p.point.distance(Point::ORIGIN) - 1.0).abs() < 1e-12);
470            assert!((p.point.distance(Point::new(1.0, 0.0)) - 1.0).abs() < 1e-12);
471        }
472        let ext = Curve::Circle(Circle::new(Point::new(2.0, 0.0), 1.0).unwrap());
473        assert_eq!(intersect(&a, &ext, tol()).points.len(), 1);
474        let same = intersect(&a, &a, tol());
475        assert!(same.overlap);
476    }
477
478    #[test]
479    fn arc_filters_by_span() {
480        let upper = Curve::Arc(Arc::new(Point::ORIGIN, 1.0, 0.0, core::f64::consts::PI).unwrap());
481        let r = intersect(&line(-2.0, 0.5, 2.0, 0.5), &upper, tol());
482        assert_eq!(r.points.len(), 2);
483        let r2 = intersect(&line(-2.0, -0.5, 2.0, -0.5), &upper, tol());
484        assert!(r2.points.is_empty());
485    }
486
487    #[test]
488    fn cubic_line_refined() {
489        let c = Curve::Cubic(CubicBez {
490            p0: Point::new(0.0, 0.0),
491            p1: Point::new(1.0, 2.0),
492            p2: Point::new(2.0, -2.0),
493            p3: Point::new(3.0, 0.0),
494        });
495        let l = line(-1.0, 0.0, 4.0, 0.0);
496        let r = intersect(&c, &l, tol());
497        assert_eq!(r.points.len(), 3, "{r:?}");
498        for p in &r.points {
499            assert!(p.point.y.abs() < 1e-9);
500            assert!(c.point_at(p.t_a).distance(p.point) < 1e-9);
501        }
502    }
503
504    #[test]
505    fn cubic_cubic() {
506        let a = Curve::Cubic(CubicBez {
507            p0: Point::new(0.0, 0.0),
508            p1: Point::new(1.0, 1.0),
509            p2: Point::new(2.0, 1.0),
510            p3: Point::new(3.0, 0.0),
511        });
512        let b = Curve::Cubic(CubicBez {
513            p0: Point::new(0.0, 0.6),
514            p1: Point::new(1.0, 0.0),
515            p2: Point::new(2.0, 0.0),
516            p3: Point::new(3.0, 0.6),
517        });
518        let r = intersect(&a, &b, tol());
519        assert_eq!(r.points.len(), 2, "{r:?}");
520        for p in &r.points {
521            assert!(a.point_at(p.t_a).distance(b.point_at(p.t_b)) < 1e-8);
522        }
523    }
524}
dotloom_geometry/intersect.rs

dotloom_geometry/
intersect.rs
