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

brep_kernel/intersect/
curve_curve_intersection.rs
