brep_kernel/brep/
pair_classification.rs

1use crate::{project_point_to_surface, KnotVector, NurbsSurface, Vec3};
2use serde::{Deserialize, Serialize};
3
4#[derive(Clone, Copy, Debug, Deserialize, Eq, PartialEq, Serialize)]
5#[serde(rename_all = "snake_case")]
6pub enum SurfacePairRelation {
7    Disjoint,
8    Candidate,
9    Transverse,
10    NearTangent,
11    Cosurface,
12    Singular,
13}
14
15#[derive(Clone, Copy, Debug, Deserialize, Serialize)]
16pub struct SurfacePairClassification {
17    pub relation: SurfacePairRelation,
18    pub minimum_separation: f64,
19    pub minimum_normal_cross: f64,
20    /// Every sample within contact distance of the other surface has
21    /// PARALLEL normals there: the surfaces only touch tangentially, so
22    /// their intersection cannot contain a transverse curve.  Marching
23    /// such a pair walks the tangency band's noise (it neither closes nor
24    /// exits — the distance-0 glue-extrude runaway); the shared topology
25    /// comes from boundary-curve exchange instead.
26    pub tangential_only: bool,
27}
28
29use crate::spatial::Aabb as Bounds;
30
31fn same_surface(first: &NurbsSurface, second: &NurbsSurface, tolerance: f64) -> bool {
32    first.degree_u == second.degree_u
33        && first.degree_v == second.degree_v
34        && first.knots_u.len() == second.knots_u.len()
35        && first.knots_v.len() == second.knots_v.len()
36        && first.control_points.len() == second.control_points.len()
37        && first
38            .knots_u
39            .iter()
40            .zip(&second.knots_u)
41            .all(|(a, b)| (a - b).abs() <= tolerance)
42        && first
43            .knots_v
44            .iter()
45            .zip(&second.knots_v)
46            .all(|(a, b)| (a - b).abs() <= tolerance)
47        && first
48            .control_points
49            .iter()
50            .zip(&second.control_points)
51            .all(|(first_row, second_row)| {
52                first_row.len() == second_row.len()
53                    && first_row.iter().zip(second_row).all(|(a, b)| {
54                        (a.x - b.x).abs() <= tolerance
55                            && (a.y - b.y).abs() <= tolerance
56                            && (a.z - b.z).abs() <= tolerance
57                            && (a.w - b.w).abs() <= tolerance
58                    })
59            })
60}
61
62fn planar_support(surface: &NurbsSurface, tolerance: f64) -> Result<Option<(Vec3, Vec3)>, String> {
63    let u = KnotVector::new(surface.knots_u.clone(), surface.degree_u)?.domain();
64    let v = KnotVector::new(surface.knots_v.clone(), surface.degree_v)?.domain();
65    let origin = surface.evaluate((u[0] + u[1]) / 2.0, (v[0] + v[1]) / 2.0)?;
66    let normal = surface.normal((u[0] + u[1]) / 2.0, (v[0] + v[1]) / 2.0)?;
67    let controls = surface
68        .control_points
69        .iter()
70        .flatten()
71        .map(|control| control.point())
72        .collect::<Result<Vec<_>, _>>()?;
73    let scale = controls
74        .iter()
75        .map(|point| point.sub(origin).length())
76        .fold(1.0f64, f64::max);
77    let plane_tolerance = (tolerance * 100.0).max(1e-7) * scale;
78    Ok(controls
79        .iter()
80        .all(|point| point.sub(origin).dot(normal).abs() <= plane_tolerance)
81        .then_some((origin, normal)))
82}
83
84
85
86/// PLANE PAIR STRADDLE (2026-11-06; `BREP_PLANE_PAIR_STRADDLE`,
87/// default on, `=0` restores the parallel reading). Two planes within the angular tolerance of parallel are
88/// still two planes when one crosses the other inside the patches: the
89/// hair-turned cube (brep-bench `tangent/cube-turned-a-hair`, 0.001 deg =
90/// 1.7e-5 rad) has side pairs whose `|n x n'|` reads 1.7e-5 <= 1e-4 and whose
91/// midpoints stand 7.6e-10 apart, so the parallel branch called them
92/// Cosurface, while their corners stand +-8.7e-5 off each other's plane. A
93/// pair is read as crossing when some control point of either patch stands
94/// more than `band` on each side of the other's plane (a degree-1 patch lies
95/// in its control hull, so the crossing line passes through it). A pair
96/// within the band on one side keeps the parallel reading.
97pub(crate) fn plane_pair_straddle_enabled() -> bool {
98    crate::switch!("BREP_PLANE_PAIR_STRADDLE", true)
99}
100
101/// Whether the planar patches `first` (through `first_origin`, normal
102/// `first_normal`) and `second` (likewise) cross each other by more than
103/// `band` on both sides: see [`plane_pair_straddle_enabled`].
104pub(crate) fn planar_patches_straddle(
105    first: &NurbsSurface,
106    first_origin: Vec3,
107    first_normal: Vec3,
108    second: &NurbsSurface,
109    second_origin: Vec3,
110    second_normal: Vec3,
111    band: f64,
112) -> bool {
113    let straddles = |surface: &NurbsSurface, origin: Vec3, normal: Vec3| {
114        let mut below = false;
115        let mut above = false;
116        for control in surface.control_points.iter().flatten() {
117            let Ok(point) = control.point() else { return false };
118            let distance = point.sub(origin).dot(normal);
119            below |= distance < -band;
120            above |= distance > band;
121        }
122        below && above
123    };
124    straddles(second, first_origin, first_normal) || straddles(first, second_origin, second_normal)
125}
126
127fn surface_samples(surface: &NurbsSurface) -> Result<Vec<(f64, f64, Vec3)>, String> {
128    let u = surface.domain_u()?;
129    let v = surface.domain_v()?;
130    let mut samples = Vec::with_capacity(25);
131    for i in 0..=4 {
132        for j in 0..=4 {
133            let uu = u[0] + (u[1] - u[0]) * i as f64 / 4.0;
134            let vv = v[0] + (v[1] - v[0]) * j as f64 / 4.0;
135            samples.push((uu, vv, surface.evaluate(uu, vv)?));
136        }
137    }
138    Ok(samples)
139}
140
141/// Per-surface data the pair classifier needs, precomputable once per face
142/// and reused across every pair that face participates in.
143pub struct SurfaceClassifyData {
144    bounds: Bounds,
145    planar: Option<(Vec3, Vec3)>,
146    /// (u, v, point, normal); normal is None where the surface is singular
147    /// (poles) — the classifier treats close singular samples specially.
148    samples: Vec<(f64, f64, Vec3, Option<Vec3>)>,
149}
150
151thread_local! {
152    static CLASSIFY_MEMO: std::cell::RefCell<rustc_hash::FxHashMap<[u64; 2], std::rc::Rc<SurfaceClassifyData>>> = Default::default();
153    static CLASSIFY_READS: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
154}
155
156pub(crate) fn classify_reads() -> usize { CLASSIFY_READS.with(std::cell::Cell::get) }
157
158pub(crate) fn clear_classify_memo() {
159    CLASSIFY_MEMO.with(|memo| memo.borrow_mut().clear());
160}
161
162impl SurfaceClassifyData {
163    /// Only geometry and tolerance enter this read. Renumbering a face or
164    /// moving it to the next boolean operand does not require new samples.
165    pub(crate) fn memoized(surface: &NurbsSurface, tolerance: f64) -> Result<std::rc::Rc<Self>, String> {
166        let enabled = std::env::var("BREP_CLASSIFY_MEMO").as_deref() != Ok("0");
167        let mut hash = crate::mass_properties::TwoLane::new();
168        crate::mass_properties::hash_surface(&mut hash, surface);
169        hash.f64(tolerance);
170        let key = hash.finish();
171        if enabled {
172            if let Some(data) = CLASSIFY_MEMO.with(|memo| memo.borrow().get(&key).cloned()) {
173                return Ok(data);
174            }
175        }
176        CLASSIFY_READS.with(|reads| reads.set(reads.get() + 1));
177        let data = std::rc::Rc::new(Self::build(surface, tolerance)?);
178        if enabled {
179            CLASSIFY_MEMO.with(|memo| {
180                let mut memo = memo.borrow_mut();
181                if memo.len() >= 4096 { memo.clear(); }
182                memo.insert(key, data.clone());
183            });
184        }
185        Ok(data)
186    }
187
188    pub fn build(surface: &NurbsSurface, model_tolerance: f64) -> Result<Self, String> {
189        let samples = surface_samples(surface)?
190            .into_iter()
191            .map(|(u, v, point)| (u, v, point, surface.normal(u, v).ok()))
192            .collect();
193        Ok(Self {
194            bounds: Bounds::from_surface_controls(surface)?,
195            planar: planar_support(surface, model_tolerance)?,
196            samples,
197        })
198    }
199}
200
201/// Classify a pair before generic SSI runs.  The classifier is conservative:
202/// only separated control hulls are declared disjoint; uncertain pairs stay
203/// `Candidate` rather than being incorrectly culled.
204pub fn classify_surface_pair(
205    first: &NurbsSurface,
206    second: &NurbsSurface,
207    model_tolerance: f64,
208    angular_tolerance: f64,
209) -> Result<SurfacePairClassification, String> {
210    let _wall = crate::wall_profile::Scope::enter(crate::wall_profile::CLASSIFY);
211    let first_data = SurfaceClassifyData::build(first, model_tolerance)?;
212    let second_data = SurfaceClassifyData::build(second, model_tolerance)?;
213    classify_surface_pair_cached(
214        first,
215        &first_data,
216        second,
217        &second_data,
218        model_tolerance,
219        angular_tolerance,
220    )
221}
222
223/// `classify_surface_pair` with the per-surface work hoisted out — callers
224/// with many pairs per face build each face's `SurfaceClassifyData` once.
225pub fn classify_surface_pair_cached(
226    first: &NurbsSurface,
227    first_data: &SurfaceClassifyData,
228    second: &NurbsSurface,
229    second_data: &SurfaceClassifyData,
230    model_tolerance: f64,
231    angular_tolerance: f64,
232) -> Result<SurfacePairClassification, String> {
233    let _wall = crate::wall_profile::Scope::enter(crate::wall_profile::CLASSIFY);
234    let first_bounds = first_data.bounds;
235    let second_bounds = second_data.bounds;
236    let hull_separation = first_bounds.separation(second_bounds);
237    let scale = first_bounds.pair_scale(second_bounds);
238    let contact = (model_tolerance * 100.0).max(scale * 1e-7);
239    if hull_separation > contact {
240        return Ok(SurfacePairClassification {
241            relation: SurfacePairRelation::Disjoint,
242            minimum_separation: hull_separation,
243            minimum_normal_cross: 1.0,
244            tangential_only: false,
245        });
246    }
247    if same_surface(first, second, model_tolerance) {
248        return Ok(SurfacePairClassification {
249            relation: SurfacePairRelation::Cosurface,
250            minimum_separation: 0.0,
251            minimum_normal_cross: 0.0,
252            tangential_only: true,
253        });
254    }
255    if let (Some((first_origin, first_normal)), Some((second_origin, second_normal))) =
256        (first_data.planar, second_data.planar)
257    {
258        let normal_cross = first_normal.cross(second_normal).length();
259        let crossing = normal_cross <= angular_tolerance
260            && plane_pair_straddle_enabled()
261            && planar_patches_straddle(first, first_origin, first_normal, second, second_origin, second_normal, contact);
262        if normal_cross <= angular_tolerance && !crossing {
263            let separation = second_origin.sub(first_origin).dot(first_normal).abs();
264            return Ok(SurfacePairClassification {
265                relation: if separation <= contact {
266                    SurfacePairRelation::Cosurface
267                } else {
268                    SurfacePairRelation::Disjoint
269                },
270                minimum_separation: separation,
271                minimum_normal_cross: normal_cross,
272                tangential_only: separation <= contact,
273            });
274        }
275        return Ok(SurfacePairClassification {
276            relation: SurfacePairRelation::Transverse,
277            minimum_separation: 0.0,
278            minimum_normal_cross: normal_cross,
279            tangential_only: false,
280        });
281    }
282
283    let mut minimum_separation = f64::INFINITY;
284    let mut minimum_normal_cross = f64::INFINITY;
285    let mut close_samples = 0usize;
286    let mut parallel_close_samples = 0usize;
287    let mut singular_close_sample = false;
288    for (source_data, target) in [(first_data, second), (second_data, first)] {
289        for &(_, _, point, source_normal) in &source_data.samples {
290            let projection = project_point_to_surface(target, point)?;
291            minimum_separation = minimum_separation.min(projection.distance);
292            if projection.distance > contact {
293                continue;
294            }
295            close_samples += 1;
296            match (source_normal, target.normal(projection.u, projection.v)) {
297                (Some(first_normal), Ok(second_normal)) => {
298                    let cross = first_normal.cross(second_normal).length();
299                    minimum_normal_cross = minimum_normal_cross.min(cross);
300                    if cross <= angular_tolerance {
301                        parallel_close_samples += 1;
302                    }
303                }
304                _ => singular_close_sample = true,
305            }
306        }
307    }
308    if minimum_separation == f64::INFINITY {
309        minimum_separation = hull_separation;
310    }
311    if minimum_normal_cross == f64::INFINITY {
312        minimum_normal_cross = 1.0;
313    }
314    let relation = if close_samples == 0 {
315        SurfacePairRelation::Candidate
316    } else if singular_close_sample {
317        SurfacePairRelation::Singular
318    } else if close_samples >= 18 && parallel_close_samples == close_samples {
319        SurfacePairRelation::Cosurface
320    } else if minimum_normal_cross <= angular_tolerance {
321        SurfacePairRelation::NearTangent
322    } else {
323        SurfacePairRelation::Transverse
324    };
325    Ok(SurfacePairClassification {
326        relation,
327        minimum_separation,
328        minimum_normal_cross,
329        tangential_only: close_samples > 0 && parallel_close_samples == close_samples,
330    })
331}
332
333
brep_kernel/brep/pair_classification.rs

brep_kernel/brep/
pair_classification.rs
