brep_kernel/intersect/
curve_curve_intersection.rs

1use crate::curve::KNOT_IDENTITY_TOL as KNOT_TOLERANCE;
2use crate::{NurbsCurve, Vec3};
3use serde::Serialize;
4
5const EPSILON: f64 = 1e-12;
6
7#[derive(Clone, Copy)]
8struct Bounds {
9    minimum: Vec3,
10    maximum: Vec3,
11}
12
13impl Bounds {
14    fn from_curve(curve: &NurbsCurve) -> Result<Self, String> {
15        let points = curve
16            .control_points
17            .iter()
18            .map(|point| point.point())
19            .collect::<Result<Vec<_>, _>>()?;
20        let mut minimum = Vec3::new(f64::INFINITY, f64::INFINITY, f64::INFINITY);
21        let mut maximum = Vec3::new(f64::NEG_INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY);
22        for point in points {
23            minimum.x = minimum.x.min(point.x);
24            minimum.y = minimum.y.min(point.y);
25            minimum.z = minimum.z.min(point.z);
26            maximum.x = maximum.x.max(point.x);
27            maximum.y = maximum.y.max(point.y);
28            maximum.z = maximum.z.max(point.z);
29        }
30        Ok(Self { minimum, maximum })
31    }
32
33    fn intersects(self, other: Self, tolerance: f64) -> bool {
34        self.minimum.x <= other.maximum.x + tolerance
35            && self.maximum.x + tolerance >= other.minimum.x
36            && self.minimum.y <= other.maximum.y + tolerance
37            && self.maximum.y + tolerance >= other.minimum.y
38            && self.minimum.z <= other.maximum.z + tolerance
39            && self.maximum.z + tolerance >= other.minimum.z
40    }
41}
42
43#[derive(Clone, Copy)]
44struct Segment {
45    start: f64,
46    end: f64,
47    bounds: Bounds,
48}
49
50fn interior_knots(curve: &NurbsCurve) -> Vec<f64> {
51    let start = curve.knots[curve.degree];
52    let end = curve.knots[curve.knots.len() - 1 - curve.degree];
53    let mut result = Vec::new();
54    for &knot in &curve.knots {
55        if knot <= start + KNOT_TOLERANCE || knot >= end - KNOT_TOLERANCE {
56            continue;
57        }
58        if result
59            .last()
60            .is_none_or(|previous: &f64| (*previous - knot).abs() > KNOT_TOLERANCE)
61        {
62            result.push(knot);
63        }
64    }
65    result
66}
67
68fn segment_boxes(curve: &NurbsCurve, count: usize) -> Result<Vec<Segment>, String> {
69    let [start, end] = curve.domain()?;
70    let mut parameters: Vec<f64> = (0..=count)
71        .map(|index| start + (end - start) * index as f64 / count as f64)
72        .collect();
73    parameters.extend(interior_knots(curve));
74    parameters.sort_by(f64::total_cmp);
75    parameters.dedup_by(|a, b| (*a - *b).abs() <= KNOT_TOLERANCE);
76
77    let mut segments = Vec::new();
78    let mut rest = curve.clone();
79    let mut rest_start = start;
80    for (index, &parameter) in parameters.iter().enumerate().skip(1) {
81        if parameter - rest_start <= KNOT_TOLERANCE {
82            continue;
83        }
84        if index == parameters.len() - 1 {
85            segments.push(Segment {
86                start: rest_start,
87                end,
88                bounds: Bounds::from_curve(&rest)?,
89            });
90        } else {
91            let (piece, remainder) = rest.split(parameter)?;
92            segments.push(Segment {
93                start: rest_start,
94                end: parameter,
95                bounds: Bounds::from_curve(&piece)?,
96            });
97            rest = remainder;
98            rest_start = parameter;
99        }
100    }
101    Ok(segments)
102}
103
104fn newton_closest_pair(
105    first: &NurbsCurve,
106    second: &NurbsCurve,
107    seed_s: f64,
108    seed_t: f64,
109) -> Result<[f64; 2], String> {
110    let [s0, s1] = first.domain()?;
111    let [t0, t1] = second.domain()?;
112    let mut s = seed_s;
113    let mut t = seed_t;
114    for _ in 0..40 {
115        let first_derivatives = first.derivatives_small(s, 2)?;
116        let second_derivatives = second.derivatives_small(t, 2)?;
117        let residual = first_derivatives[0].sub(second_derivatives[0]);
118        let f = first_derivatives[1].dot(residual);
119        let g = -second_derivatives[1].dot(residual);
120        let j00 = first_derivatives[2].dot(residual) + first_derivatives[1].length_squared();
121        let j01 = -first_derivatives[1].dot(second_derivatives[1]);
122        let j10 = -second_derivatives[1].dot(first_derivatives[1]);
123        let j11 = -second_derivatives[2].dot(residual) + second_derivatives[1].length_squared();
124        let residual_length = residual.length();
125        let converged_first =
126            f.abs() <= EPSILON + 1e-10 * first_derivatives[1].length() * residual_length.max(1e-30);
127        let converged_second = g.abs()
128            <= EPSILON + 1e-10 * second_derivatives[1].length() * residual_length.max(1e-30);
129        if (converged_first && converged_second) || residual_length <= EPSILON {
130            return Ok([s, t]);
131        }
132        let determinant = j00 * j11 - j01 * j10;
133        if determinant.abs() <= EPSILON {
134            return Ok([s, t]);
135        }
136        let mut ds = (-f * j11 + g * j01) / determinant;
137        let mut dt = (-g * j00 + f * j10) / determinant;
138        ds = ds.clamp(-(s1 - s0) / 4.0, (s1 - s0) / 4.0);
139        dt = dt.clamp(-(t1 - t0) / 4.0, (t1 - t0) / 4.0);
140        let next_s = (s + ds).clamp(s0, s1);
141        let next_t = (t + dt).clamp(t0, t1);
142        if (next_s - s).abs() <= 1e-15 * (s1 - s0) && (next_t - t).abs() <= 1e-15 * (t1 - t0) {
143            return Ok([next_s, next_t]);
144        }
145        s = next_s;
146        t = next_t;
147    }
148    Ok([s, t])
149}
150
151fn parameter_tolerance(curve: &NurbsCurve, parameter: f64, tolerance: f64) -> Result<f64, String> {
152    let speed = curve.deriv1(parameter)?.1.length();
153    let [start, end] = curve.domain()?;
154    if speed <= EPSILON {
155        Ok((end - start) * 1e-3)
156    } else {
157        Ok(((tolerance / speed) * 10.0).max((end - start) * 1e-9))
158    }
159}
160
161fn parameter_distance(a: f64, b: f64, start: f64, end: f64, closed: bool) -> f64 {
162    let distance = (a - b).abs();
163    if closed {
164        distance.min(end - start - distance)
165    } else {
166        distance
167    }
168}
169
170#[derive(Clone, Copy, Debug, Serialize)]
171pub struct CurveCurveIntersection {
172    pub s: f64,
173    pub t: f64,
174    pub point: Vec3,
175    pub gap: f64,
176    pub tangential: bool,
177}
178
179/// `BREP_EXACT_LINE_LINE` (lane X row 7): two straight segments that are not
180/// near parallel meet where their exact closest pair says
181/// ([`exact_line_line`]), instead of where Newton runs from box seeds settle.
182/// On by default; `0`/`off` returns to the general search, and `census`
183/// answers by the general search and appends one line per call whose answers
184/// differ to `BREP_EXACT_LINE_LINE_CENSUS`.
185///
186/// Promoted 2026-10-16 on top of the canonical fragment order. 239 documents:
187/// 6,853 calls change, every one with the same hit count; the exact pair
188/// removes residuals up to 2.02e-8, and the per-hit gap goes from a 2.17e-8
189/// maximum to 7.9e-9. The exact pair is exact by construction, and the 1,387
190/// hits whose evaluated gap reads above the general search's, by at most
191/// 1.21e-13 (p50 8.9e-16), are the rounding floor of measuring the gap, not
192/// lost accuracy. 49 documents change bytes; every body is geometrically the
193/// same, with no name moved.
194#[derive(Clone, Copy, PartialEq, Eq)]
195enum LineLineMode {
196    Off,
197    On,
198    Census,
199}
200
201fn line_line_mode() -> LineLineMode {
202    static MODE: std::sync::OnceLock<LineLineMode> = std::sync::OnceLock::new();
203    let mode = *MODE.get_or_init(|| match std::env::var("BREP_EXACT_LINE_LINE").as_deref() {
204        Ok("0") | Ok("off") => LineLineMode::Off,
205        Ok("census") => LineLineMode::Census,
206        _ => LineLineMode::On,
207    });
208    mode
209}
210
211
212/// The segment a degree-1 curve of two equal-weight control points traces
213/// linearly over its clamped knots: (start point, end point, domain).
214fn straight_segment(curve: &NurbsCurve) -> Option<(Vec3, Vec3, [f64; 2])> {
215    if curve.degree != 1 || curve.control_points.len() != 2 || curve.knots.len() != 4 {
216        return None;
217    }
218    let k = &curve.knots;
219    if k[0] != k[1] || k[2] != k[3] || !(k[1] < k[2]) {
220        return None;
221    }
222    if curve.control_points[0].w != curve.control_points[1].w {
223        return None;
224    }
225    Some((curve.control_points[0].point().ok()?, curve.control_points[1].point().ok()?, [k[1], k[2]]))
226}
227
228/// The exact meeting of two straight segments that are not near parallel:
229/// their unique closest pair (each parameter clamped to its segment, the
230/// other then re-solved), one transverse hit when its gap is within
231/// `tolerance`, none otherwise. `None` (the general search answers) for
232/// anything else: a curve that is not a straight segment, a degenerate one,
233/// or two directions within the tangential band of `intersect_curves`.
234pub(crate) fn exact_line_line(first: &NurbsCurve, second: &NurbsCurve, tolerance: f64) -> Option<Vec<CurveCurveIntersection>> {
235    let (p0, p1, [s0, s1]) = straight_segment(first)?;
236    let (q0, q1, [t0, t1]) = straight_segment(second)?;
237    let d1 = p1.sub(p0);
238    let d2 = q1.sub(q0);
239    let a = d1.dot(d1);
240    let e = d2.dot(d2);
241    if a <= EPSILON || e <= EPSILON {
242        return None;
243    }
244    if d1.cross(d2).length() <= 1e-3_f64.sin() * (a * e).sqrt() {
245        return None;
246    }
247    let r = p0.sub(q0);
248    let b = d1.dot(d2);
249    let c = d1.dot(r);
250    let f = d2.dot(r);
251    let denominator = a * e - b * b;
252    let mut s = ((b * f - c * e) / denominator).clamp(0.0, 1.0);
253    let mut t = (b * s + f) / e;
254    if t < 0.0 {
255        t = 0.0;
256        s = (-c / a).clamp(0.0, 1.0);
257    } else if t > 1.0 {
258        t = 1.0;
259        s = ((b - c) / a).clamp(0.0, 1.0);
260    }
261    let first_point = p0.add(d1.scale(s));
262    let second_point = q0.add(d2.scale(t));
263    let gap = first_point.sub(second_point).length();
264    if gap > tolerance {
265        return Some(Vec::new());
266    }
267    Some(vec![CurveCurveIntersection {
268        s: (s0 + (s1 - s0) * s).clamp(s0, s1),
269        t: (t0 + (t1 - t0) * t).clamp(t0, t1),
270        point: first_point.add(second_point).scale(0.5),
271        gap,
272        tangential: false,
273    }])
274}
275
276pub fn intersect_curves(
277    first: &NurbsCurve,
278    second: &NurbsCurve,
279    tolerance: f64,
280) -> Result<Vec<CurveCurveIntersection>, String> {
281    let mode = line_line_mode();
282    if mode != LineLineMode::Off {
283        if let Some(exact) = exact_line_line(first, second, tolerance) {
284            if mode == LineLineMode::On {
285                return Ok(exact);
286            }
287            let general = intersect_curves_general(first, second, tolerance)?;
288            line_line_census(first, second, &general, &exact)?;
289            return Ok(general);
290        }
291    }
292    intersect_curves_general(first, second, tolerance)
293}
294
295/// One census line per call whose general and exact answers differ: both
296/// hit lists, and for each general hit its points' distances from the exact
297/// closest pair (the residual the exact answer removes).
298fn line_line_census(
299    first: &NurbsCurve,
300    second: &NurbsCurve,
301    general: &[CurveCurveIntersection],
302    exact: &[CurveCurveIntersection],
303) -> Result<(), String> {
304    let same = general.len() == exact.len()
305        && general.iter().zip(exact).all(|(g, x)| g.s.to_bits() == x.s.to_bits() && g.t.to_bits() == x.t.to_bits());
306    if same {
307        return Ok(());
308    }
309    let Ok(path) = std::env::var("BREP_EXACT_LINE_LINE_CENSUS") else {
310        return Ok(());
311    };
312    let mut line = format!("line-line census: label {} general {} exact {}", std::env::var("BREP_CENSUS_LABEL").unwrap_or_default(), general.len(), exact.len());
313    if let Some(x) = exact.first() {
314        let (exact_first, exact_second) = (first.evaluate(x.s)?, second.evaluate(x.t)?);
315        line.push_str(&format!(" exact s={:e} t={:e} gap={:e}", x.s, x.t, x.gap));
316        for g in general {
317            let residual = first.evaluate(g.s)?.sub(exact_first).length().max(second.evaluate(g.t)?.sub(exact_second).length());
318            line.push_str(&format!(" | general s={:e} t={:e} gap={:e} ds={:e} dt={:e} residual={:e}", g.s, g.t, g.gap, g.s - x.s, g.t - x.t, residual));
319        }
320    } else {
321        for g in general {
322            line.push_str(&format!(" | general s={:e} t={:e} gap={:e}", g.s, g.t, g.gap));
323        }
324    }
325    use std::io::Write;
326    if let Ok(mut file) = std::fs::OpenOptions::new().create(true).append(true).open(path) {
327        let _ = writeln!(file, "{line}");
328    }
329    Ok(())
330}
331
332/// Every `(i, j)` with `first[i].bounds` and `second[j].bounds` within
333/// `tolerance` of each other, in ascending `(i, j)` order — exactly the pairs
334/// and the order the nested loop produced. The second list is sorted by its
335/// boxes' minimum x once; for each first segment the candidates are the
336/// second segments whose minimum x lies at or below the first's maximum x
337/// plus the tolerance (a prefix of the sorted order), filtered by the full
338/// box test. Each pair is tested by the same `Aabb::intersects` as before.
339fn overlapping_segment_pairs(first: &[Segment], second: &[Segment], tolerance: f64) -> Vec<(usize, usize)> {
340    let mut order: Vec<usize> = (0..second.len()).collect();
341    order.sort_by(|&a, &b| second[a].bounds.minimum.x.total_cmp(&second[b].bounds.minimum.x));
342    let minimum_x: Vec<f64> = order.iter().map(|&j| second[j].bounds.minimum.x).collect();
343    let mut pairs = Vec::new();
344    for (i, segment) in first.iter().enumerate() {
345        let limit = segment.bounds.maximum.x + tolerance;
346        // Candidates: the prefix with minimum x <= limit.
347        let end = minimum_x.partition_point(|&x| x <= limit);
348        let mut row: Vec<usize> = order[..end]
349            .iter()
350            .copied()
351            .filter(|&j| segment.bounds.intersects(second[j].bounds, tolerance))
352            .collect();
353        row.sort_unstable();
354        pairs.extend(row.into_iter().map(|j| (i, j)));
355    }
356    pairs
357}
358
359fn intersect_curves_general(
360    first: &NurbsCurve,
361    second: &NurbsCurve,
362    tolerance: f64,
363) -> Result<Vec<CurveCurveIntersection>, String> {
364    let first_segment_count = 8usize.max((interior_knots(first).len() + 1) * (first.degree + 1));
365    let second_segment_count = 8usize.max((interior_knots(second).len() + 1) * (second.degree + 1));
366    let first_segments = segment_boxes(first, first_segment_count)?;
367    let second_segments = segment_boxes(second, second_segment_count)?;
368    // The overlapping segment pairs, in the order the nested loop over
369    // `first_segments` then `second_segments` visits them: the seeds' order
370    // decides which of two converging seeds is kept as the hit below, so the
371    // pairs are gathered by a sweep over the x extents (sorted once, a binary
372    // search per first segment) and then put back into that order. Two fitted
373    // sections of two thousand segments each cost four million box tests
374    // a call before (2026-10-17, the shelled-box report's shell: 0.57 s of
375    // `lies_on` spans and 0.55 s of section crossings on three face pairs).
376    let seeds = overlapping_segment_pairs(&first_segments, &second_segments, tolerance)
377        .into_iter()
378        .map(|(i, j)| {
379            [
380                (first_segments[i].start + first_segments[i].end) * 0.5,
381                (second_segments[j].start + second_segments[j].end) * 0.5,
382            ]
383        })
384        .collect::<Vec<_>>();
385    let [s0, s1] = first.domain()?;
386    let [t0, t1] = second.domain()?;
387    let closed_first = first.evaluate(s0)?.sub(first.evaluate(s1)?).length() <= tolerance;
388    let closed_second = second.evaluate(t0)?.sub(second.evaluate(t1)?).length() <= tolerance;
389    let mut results: Vec<CurveCurveIntersection> = Vec::new();
390    for [seed_s, seed_t] in seeds {
391        let [s, t] = newton_closest_pair(first, second, seed_s, seed_t)?;
392        let first_point = first.evaluate(s)?;
393        let second_point = second.evaluate(t)?;
394        let gap = first_point.sub(second_point).length();
395        if gap > tolerance {
396            continue;
397        }
398        let first_tangent = first.deriv1(s)?.1;
399        let second_tangent = second.deriv1(t)?.1;
400        let tangential = match (first_tangent.normalized(), second_tangent.normalized()) {
401            (Ok(a), Ok(b)) => a.cross(b).length() <= 1e-3_f64.sin(),
402            _ => false,
403        };
404        let point = first_point.add(second_point).scale(0.5);
405        let scale = 1.0 + point.length();
406        let spatial_tolerance = if tangential {
407            tolerance.sqrt() * scale
408        } else {
409            tolerance * 10.0
410        };
411        let mut duplicate = false;
412        for result in &results {
413            if (parameter_distance(result.s, s, s0, s1, closed_first)
414                <= parameter_tolerance(first, result.s, tolerance)?
415                && parameter_distance(result.t, t, t0, t1, closed_second)
416                    <= parameter_tolerance(second, result.t, tolerance)?)
417                || result.point.sub(point).length()
418                    <= if result.tangential || tangential {
419                        spatial_tolerance
420                    } else {
421                        tolerance * 10.0
422                    }
423            {
424                duplicate = true;
425                break;
426            }
427        }
428        if duplicate {
429            continue;
430        }
431        results.push(CurveCurveIntersection {
432            s: s.clamp(s0, s1),
433            t: t.clamp(t0, t1),
434            point,
435            gap,
436            tangential,
437        });
438    }
439    results.sort_by(|a, b| a.s.total_cmp(&b.s));
440    Ok(results)
441}
442
brep_kernel/intersect/curve_curve_intersection.rs

brep_kernel/intersect/
curve_curve_intersection.rs
