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
86fn surface_samples(surface: &NurbsSurface) -> Result<Vec<(f64, f64, Vec3)>, String> {
87    let u = surface.domain_u()?;
88    let v = surface.domain_v()?;
89    let mut samples = Vec::with_capacity(25);
90    for i in 0..=4 {
91        for j in 0..=4 {
92            let uu = u[0] + (u[1] - u[0]) * i as f64 / 4.0;
93            let vv = v[0] + (v[1] - v[0]) * j as f64 / 4.0;
94            samples.push((uu, vv, surface.evaluate(uu, vv)?));
95        }
96    }
97    Ok(samples)
98}
99
100/// Per-surface data the pair classifier needs, precomputable once per face
101/// and reused across every pair that face participates in.
102pub struct SurfaceClassifyData {
103    bounds: Bounds,
104    planar: Option<(Vec3, Vec3)>,
105    /// (u, v, point, normal); normal is None where the surface is singular
106    /// (poles) — the classifier treats close singular samples specially.
107    samples: Vec<(f64, f64, Vec3, Option<Vec3>)>,
108}
109
110thread_local! {
111    static CLASSIFY_MEMO: std::cell::RefCell<rustc_hash::FxHashMap<[u64; 2], std::rc::Rc<SurfaceClassifyData>>> = Default::default();
112    static CLASSIFY_READS: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
113}
114
115pub(crate) fn classify_reads() -> usize { CLASSIFY_READS.with(std::cell::Cell::get) }
116
117pub(crate) fn clear_classify_memo() {
118    CLASSIFY_MEMO.with(|memo| memo.borrow_mut().clear());
119}
120
121impl SurfaceClassifyData {
122    /// Only geometry and tolerance enter this read. Renumbering a face or
123    /// moving it to the next boolean operand does not require new samples.
124    pub(crate) fn memoized(surface: &NurbsSurface, tolerance: f64) -> Result<std::rc::Rc<Self>, String> {
125        let enabled = std::env::var("BREP_CLASSIFY_MEMO").as_deref() != Ok("0");
126        let mut hash = crate::mass_properties::TwoLane::new();
127        crate::mass_properties::hash_surface(&mut hash, surface);
128        hash.f64(tolerance);
129        let key = hash.finish();
130        if enabled {
131            if let Some(data) = CLASSIFY_MEMO.with(|memo| memo.borrow().get(&key).cloned()) {
132                return Ok(data);
133            }
134        }
135        CLASSIFY_READS.with(|reads| reads.set(reads.get() + 1));
136        let data = std::rc::Rc::new(Self::build(surface, tolerance)?);
137        if enabled {
138            CLASSIFY_MEMO.with(|memo| {
139                let mut memo = memo.borrow_mut();
140                if memo.len() >= 4096 { memo.clear(); }
141                memo.insert(key, data.clone());
142            });
143        }
144        Ok(data)
145    }
146
147    pub fn build(surface: &NurbsSurface, model_tolerance: f64) -> Result<Self, String> {
148        let samples = surface_samples(surface)?
149            .into_iter()
150            .map(|(u, v, point)| (u, v, point, surface.normal(u, v).ok()))
151            .collect();
152        Ok(Self {
153            bounds: Bounds::from_surface_controls(surface)?,
154            planar: planar_support(surface, model_tolerance)?,
155            samples,
156        })
157    }
158}
159
160/// Classify a pair before generic SSI runs.  The classifier is conservative:
161/// only separated control hulls are declared disjoint; uncertain pairs stay
162/// `Candidate` rather than being incorrectly culled.
163pub fn classify_surface_pair(
164    first: &NurbsSurface,
165    second: &NurbsSurface,
166    model_tolerance: f64,
167    angular_tolerance: f64,
168) -> Result<SurfacePairClassification, String> {
169    let first_data = SurfaceClassifyData::build(first, model_tolerance)?;
170    let second_data = SurfaceClassifyData::build(second, model_tolerance)?;
171    classify_surface_pair_cached(
172        first,
173        &first_data,
174        second,
175        &second_data,
176        model_tolerance,
177        angular_tolerance,
178    )
179}
180
181/// `classify_surface_pair` with the per-surface work hoisted out — callers
182/// with many pairs per face build each face's `SurfaceClassifyData` once.
183pub fn classify_surface_pair_cached(
184    first: &NurbsSurface,
185    first_data: &SurfaceClassifyData,
186    second: &NurbsSurface,
187    second_data: &SurfaceClassifyData,
188    model_tolerance: f64,
189    angular_tolerance: f64,
190) -> Result<SurfacePairClassification, String> {
191    let first_bounds = first_data.bounds;
192    let second_bounds = second_data.bounds;
193    let hull_separation = first_bounds.separation(second_bounds);
194    let scale = first_bounds.pair_scale(second_bounds);
195    let contact = (model_tolerance * 100.0).max(scale * 1e-7);
196    if hull_separation > contact {
197        return Ok(SurfacePairClassification {
198            relation: SurfacePairRelation::Disjoint,
199            minimum_separation: hull_separation,
200            minimum_normal_cross: 1.0,
201            tangential_only: false,
202        });
203    }
204    if same_surface(first, second, model_tolerance) {
205        return Ok(SurfacePairClassification {
206            relation: SurfacePairRelation::Cosurface,
207            minimum_separation: 0.0,
208            minimum_normal_cross: 0.0,
209            tangential_only: true,
210        });
211    }
212    if let (Some((first_origin, first_normal)), Some((second_origin, second_normal))) =
213        (first_data.planar, second_data.planar)
214    {
215        let normal_cross = first_normal.cross(second_normal).length();
216        if normal_cross <= angular_tolerance {
217            let separation = second_origin.sub(first_origin).dot(first_normal).abs();
218            return Ok(SurfacePairClassification {
219                relation: if separation <= contact {
220                    SurfacePairRelation::Cosurface
221                } else {
222                    SurfacePairRelation::Disjoint
223                },
224                minimum_separation: separation,
225                minimum_normal_cross: normal_cross,
226                tangential_only: separation <= contact,
227            });
228        }
229        return Ok(SurfacePairClassification {
230            relation: SurfacePairRelation::Transverse,
231            minimum_separation: 0.0,
232            minimum_normal_cross: normal_cross,
233            tangential_only: false,
234        });
235    }
236
237    let mut minimum_separation = f64::INFINITY;
238    let mut minimum_normal_cross = f64::INFINITY;
239    let mut close_samples = 0usize;
240    let mut parallel_close_samples = 0usize;
241    let mut singular_close_sample = false;
242    for (source_data, target) in [(first_data, second), (second_data, first)] {
243        for &(_, _, point, source_normal) in &source_data.samples {
244            let projection = project_point_to_surface(target, point)?;
245            minimum_separation = minimum_separation.min(projection.distance);
246            if projection.distance > contact {
247                continue;
248            }
249            close_samples += 1;
250            match (source_normal, target.normal(projection.u, projection.v)) {
251                (Some(first_normal), Ok(second_normal)) => {
252                    let cross = first_normal.cross(second_normal).length();
253                    minimum_normal_cross = minimum_normal_cross.min(cross);
254                    if cross <= angular_tolerance {
255                        parallel_close_samples += 1;
256                    }
257                }
258                _ => singular_close_sample = true,
259            }
260        }
261    }
262    if minimum_separation == f64::INFINITY {
263        minimum_separation = hull_separation;
264    }
265    if minimum_normal_cross == f64::INFINITY {
266        minimum_normal_cross = 1.0;
267    }
268    let relation = if close_samples == 0 {
269        SurfacePairRelation::Candidate
270    } else if singular_close_sample {
271        SurfacePairRelation::Singular
272    } else if close_samples >= 18 && parallel_close_samples == close_samples {
273        SurfacePairRelation::Cosurface
274    } else if minimum_normal_cross <= angular_tolerance {
275        SurfacePairRelation::NearTangent
276    } else {
277        SurfacePairRelation::Transverse
278    };
279    Ok(SurfacePairClassification {
280        relation,
281        minimum_separation,
282        minimum_normal_cross,
283        tangential_only: close_samples > 0 && parallel_close_samples == close_samples,
284    })
285}
286
287
brep_kernel/brep/pair_classification.rs

brep_kernel/brep/
pair_classification.rs
