brepkit_math/nurbs/intersection/
curve_surface.rs

1//! Curve-surface intersection.
2
3use crate::MathError;
4use crate::nurbs::curve::NurbsCurve;
5use crate::nurbs::decompose::{curve_to_bezier_segments, surface_to_bezier_patches};
6use crate::nurbs::surface::NurbsSurface;
7use crate::vec::{Point3, Vec3};
8
9use super::MAX_NEWTON_ITER;
10
11/// A curve-surface intersection hit with parameters on both entities.
12#[derive(Debug, Clone, Copy)]
13pub struct CurveSurfaceHit {
14    /// 3D position of the intersection.
15    pub point: Point3,
16    /// Parameter on the curve.
17    pub t: f64,
18    /// Parameters on the surface (u, v).
19    pub uv: (f64, f64),
20}
21
22/// Find all intersection points between a NURBS curve and a NURBS surface.
23///
24/// Decomposes both into Bezier segments/patches, performs AABB filtering,
25/// then refines candidates with Newton iteration on the coupled 3x3 system
26/// `C(t) - S(u,v) = 0`.
27///
28/// # Parameters
29///
30/// - `curve`: The NURBS curve
31/// - `surface`: The NURBS surface
32/// - `tolerance`: Distance tolerance for convergence (e.g. 1e-7)
33///
34/// # Errors
35///
36/// Returns an error if Bezier decomposition fails.
37#[allow(clippy::too_many_lines)]
38pub fn intersect_curve_surface(
39    curve: &NurbsCurve,
40    surface: &NurbsSurface,
41    tolerance: f64,
42) -> Result<Vec<CurveSurfaceHit>, MathError> {
43    let segments = curve_to_bezier_segments(curve)?;
44    let patches = surface_to_bezier_patches(surface)?;
45
46    let mut hits = Vec::new();
47    // Use a generous dedup threshold: many seeds converge to the same root,
48    // and we want to collapse them. 1e-6 in 3D space is well below any
49    // meaningful geometric distinction at CAD tolerances.
50    let dedup_dist = tolerance.max(1e-6) * 10.0;
51    let dedup_sq = dedup_dist * dedup_dist;
52
53    for seg in &segments {
54        let seg_aabb = seg.aabb();
55        let (t_lo, t_hi) = seg.domain();
56
57        for patch in &patches {
58            let patch_aabb = patch.aabb();
59            if !seg_aabb.intersects(patch_aabb) {
60                continue;
61            }
62
63            // Sample the curve segment and find seed points where C(t) is
64            // close to the surface patch.
65            let n_samples = 8;
66            let (u0, u1) = patch.u_range;
67            let (v0, v1) = patch.v_range;
68
69            for i in 0..=n_samples {
70                let frac = i as f64 / n_samples as f64;
71                let t = t_lo + (t_hi - t_lo) * frac;
72                let pt = seg.evaluate(t);
73
74                // Find closest (u, v) on this patch to pt via a quick
75                // grid search over a 3x3 sub-grid of the patch.
76                let mut best_u = (u0 + u1) * 0.5;
77                let mut best_v = (v0 + v1) * 0.5;
78                let mut best_dist = f64::MAX;
79                for iu in 0..=2_usize {
80                    let uf = u0 + (u1 - u0) * (iu as f64 / 2.0);
81                    for iv in 0..=2_usize {
82                        let vf = v0 + (v1 - v0) * (iv as f64 / 2.0);
83                        let sp = surface.evaluate(uf, vf);
84                        let d = (sp - pt).length();
85                        if d < best_dist {
86                            best_dist = d;
87                            best_u = uf;
88                            best_v = vf;
89                        }
90                    }
91                }
92                let u_guess = best_u;
93                let v_guess = best_v;
94
95                if let Some(hit) =
96                    refine_curve_surface_point(curve, surface, t, u_guess, v_guess, tolerance)
97                {
98                    // Deduplicate against existing hits.
99                    let dup = hits.iter().any(|h: &CurveSurfaceHit| {
100                        let d = h.point - hit.point;
101                        d.x().mul_add(d.x(), d.y().mul_add(d.y(), d.z() * d.z())) < dedup_sq
102                    });
103                    if !dup {
104                        hits.push(hit);
105                    }
106                }
107            }
108        }
109    }
110
111    // Sort by curve parameter for deterministic output.
112    hits.sort_by(|a, b| a.t.partial_cmp(&b.t).unwrap_or(std::cmp::Ordering::Equal));
113    Ok(hits)
114}
115
116/// Refine a curve-surface intersection seed `(t, u, v)` via Newton iteration
117/// on the 3x3 system `F(t,u,v) = C(t) - S(u,v) = 0`.
118fn refine_curve_surface_point(
119    curve: &NurbsCurve,
120    surface: &NurbsSurface,
121    t_guess: f64,
122    u_guess: f64,
123    v_guess: f64,
124    tolerance: f64,
125) -> Option<CurveSurfaceHit> {
126    let mut t = t_guess;
127    let mut u = u_guess;
128    let mut v = v_guess;
129
130    let (t_min, t_max) = curve.domain();
131    let (u_min, u_max) = surface.domain_u();
132    let (v_min, v_max) = surface.domain_v();
133
134    for _ in 0..MAX_NEWTON_ITER {
135        let c_pt = curve.evaluate(t);
136        let s_pt = surface.evaluate(u, v);
137        let residual = c_pt - s_pt;
138        let dist_sq = residual.x().mul_add(
139            residual.x(),
140            residual
141                .y()
142                .mul_add(residual.y(), residual.z() * residual.z()),
143        );
144
145        if dist_sq < tolerance * tolerance {
146            return Some(CurveSurfaceHit {
147                point: c_pt,
148                t,
149                uv: (u, v),
150            });
151        }
152
153        // Jacobian columns: dC/dt, -dS/du, -dS/dv
154        let c_derivs = curve.derivatives(t, 1);
155        let ct = if c_derivs.len() > 1 {
156            c_derivs[1]
157        } else {
158            return None;
159        };
160
161        let s_derivs = surface.derivatives(u, v, 1);
162        let su = s_derivs[1][0];
163        let sv = s_derivs[0][1];
164
165        // J = [ct | -su | -sv], F = residual (as Vec3)
166        // Solve J * [dt, du, dv]^T = -F
167        let r = Vec3::new(residual.x(), residual.y(), residual.z());
168
169        // 3x3 matrix: columns are ct, -su, -sv
170        // Use Cramer's rule.
171        let col0 = ct;
172        let col1 = Vec3::new(-su.x(), -su.y(), -su.z());
173        let col2 = Vec3::new(-sv.x(), -sv.y(), -sv.z());
174
175        let det = col0.x() * (col1.y() * col2.z() - col1.z() * col2.y())
176            - col1.x() * (col0.y() * col2.z() - col0.z() * col2.y())
177            + col2.x() * (col0.y() * col1.z() - col0.z() * col1.y());
178
179        if det.abs() < 1e-30 {
180            // Singular -- try Tikhonov regularization.
181            let lambda = 1e-6;
182            let jtj = [
183                [col0.dot(col0) + lambda, col0.dot(col1), col0.dot(col2)],
184                [col1.dot(col0), col1.dot(col1) + lambda, col1.dot(col2)],
185                [col2.dot(col0), col2.dot(col1), col2.dot(col2) + lambda],
186            ];
187            let jtr = [col0.dot(r), col1.dot(r), col2.dot(r)];
188            if let Some((dt, du, dv)) = solve_3x3_cramer(&jtj, &jtr) {
189                t -= dt;
190                u -= du;
191                v -= dv;
192            } else {
193                break;
194            }
195        } else {
196            let inv_det = 1.0 / det;
197
198            // Replace each column with -r to get Cramer numerators.
199            let neg_r = Vec3::new(-r.x(), -r.y(), -r.z());
200
201            let dt = inv_det
202                * (neg_r.x() * (col1.y() * col2.z() - col1.z() * col2.y())
203                    - col1.x() * (neg_r.y() * col2.z() - neg_r.z() * col2.y())
204                    + col2.x() * (neg_r.y() * col1.z() - neg_r.z() * col1.y()));
205            let du = inv_det
206                * (col0.x() * (neg_r.y() * col2.z() - neg_r.z() * col2.y())
207                    - neg_r.x() * (col0.y() * col2.z() - col0.z() * col2.y())
208                    + col2.x() * (col0.y() * neg_r.z() - col0.z() * neg_r.y()));
209            let dv = inv_det
210                * (col0.x() * (col1.y() * neg_r.z() - col1.z() * neg_r.y())
211                    - col1.x() * (col0.y() * neg_r.z() - col0.z() * neg_r.y())
212                    + neg_r.x() * (col0.y() * col1.z() - col0.z() * col1.y()));
213
214            t += dt;
215            u += du;
216            v += dv;
217        }
218
219        t = t.clamp(t_min, t_max);
220        u = u.clamp(u_min, u_max);
221        v = v.clamp(v_min, v_max);
222    }
223
224    // Final check.
225    let c_pt = curve.evaluate(t);
226    let s_pt = surface.evaluate(u, v);
227    let residual = c_pt - s_pt;
228    let dist_sq = residual.x().mul_add(
229        residual.x(),
230        residual
231            .y()
232            .mul_add(residual.y(), residual.z() * residual.z()),
233    );
234
235    // 10x relaxation: seed points may be imprecise
236    if dist_sq < (tolerance * 10.0).powi(2) {
237        Some(CurveSurfaceHit {
238            point: c_pt,
239            t,
240            uv: (u, v),
241        })
242    } else {
243        None
244    }
245}
246
247/// Solve a 3x3 system `A * x = b` using Cramer's rule.
248fn solve_3x3_cramer(a: &[[f64; 3]; 3], b: &[f64; 3]) -> Option<(f64, f64, f64)> {
249    let det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
250        - a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0])
251        + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
252
253    if det.abs() < 1e-30 {
254        return None;
255    }
256    let inv = 1.0 / det;
257
258    let x0 = inv
259        * (b[0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
260            - a[0][1] * (b[1] * a[2][2] - a[1][2] * b[2])
261            + a[0][2] * (b[1] * a[2][1] - a[1][1] * b[2]));
262    let x1 = inv
263        * (a[0][0] * (b[1] * a[2][2] - a[1][2] * b[2])
264            - b[0] * (a[1][0] * a[2][2] - a[1][2] * a[2][0])
265            + a[0][2] * (a[1][0] * b[2] - b[1] * a[2][0]));
266    let x2 = inv
267        * (a[0][0] * (a[1][1] * b[2] - b[1] * a[2][1])
268            - a[0][1] * (a[1][0] * b[2] - b[1] * a[2][0])
269            + b[0] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]));
270
271    Some((x0, x1, x2))
272}
brepkit_math/nurbs/intersection/curve_surface.rs

brepkit_math/nurbs/intersection/
curve_surface.rs
