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
29#[derive(Clone, Copy)]
30struct Bounds {
31    low: Vec3,
32    high: Vec3,
33}
34
35fn same_surface(first: &NurbsSurface, second: &NurbsSurface, tolerance: f64) -> bool {
36    first.degree_u == second.degree_u
37        && first.degree_v == second.degree_v
38        && first.knots_u.len() == second.knots_u.len()
39        && first.knots_v.len() == second.knots_v.len()
40        && first.control_points.len() == second.control_points.len()
41        && first
42            .knots_u
43            .iter()
44            .zip(&second.knots_u)
45            .all(|(a, b)| (a - b).abs() <= tolerance)
46        && first
47            .knots_v
48            .iter()
49            .zip(&second.knots_v)
50            .all(|(a, b)| (a - b).abs() <= tolerance)
51        && first
52            .control_points
53            .iter()
54            .zip(&second.control_points)
55            .all(|(first_row, second_row)| {
56                first_row.len() == second_row.len()
57                    && first_row.iter().zip(second_row).all(|(a, b)| {
58                        (a.x - b.x).abs() <= tolerance
59                            && (a.y - b.y).abs() <= tolerance
60                            && (a.z - b.z).abs() <= tolerance
61                            && (a.w - b.w).abs() <= tolerance
62                    })
63            })
64}
65
66fn planar_support(surface: &NurbsSurface, tolerance: f64) -> Result<Option<(Vec3, Vec3)>, String> {
67    let u = KnotVector::new(surface.knots_u.clone(), surface.degree_u)?.domain();
68    let v = KnotVector::new(surface.knots_v.clone(), surface.degree_v)?.domain();
69    let origin = surface.evaluate((u[0] + u[1]) / 2.0, (v[0] + v[1]) / 2.0)?;
70    let normal = surface.normal((u[0] + u[1]) / 2.0, (v[0] + v[1]) / 2.0)?;
71    let controls = surface
72        .control_points
73        .iter()
74        .flatten()
75        .map(|control| control.point())
76        .collect::<Result<Vec<_>, _>>()?;
77    let scale = controls
78        .iter()
79        .map(|point| point.sub(origin).length())
80        .fold(1.0f64, f64::max);
81    let plane_tolerance = (tolerance * 100.0).max(1e-7) * scale;
82    Ok(controls
83        .iter()
84        .all(|point| point.sub(origin).dot(normal).abs() <= plane_tolerance)
85        .then_some((origin, normal)))
86}
87
88impl Bounds {
89    fn of(surface: &NurbsSurface) -> Result<Self, String> {
90        let mut low = Vec3::new(f64::INFINITY, f64::INFINITY, f64::INFINITY);
91        let mut high = Vec3::new(f64::NEG_INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY);
92        for control in surface.control_points.iter().flatten() {
93            let point = control.point()?;
94            low.x = low.x.min(point.x);
95            low.y = low.y.min(point.y);
96            low.z = low.z.min(point.z);
97            high.x = high.x.max(point.x);
98            high.y = high.y.max(point.y);
99            high.z = high.z.max(point.z);
100        }
101        Ok(Self { low, high })
102    }
103
104    fn separation(self, other: Self) -> f64 {
105        let axis_gap = |a0: f64, a1: f64, b0: f64, b1: f64| {
106            if a1 < b0 {
107                b0 - a1
108            } else if b1 < a0 {
109                a0 - b1
110            } else {
111                0.0
112            }
113        };
114        Vec3::new(
115            axis_gap(self.low.x, self.high.x, other.low.x, other.high.x),
116            axis_gap(self.low.y, self.high.y, other.low.y, other.high.y),
117            axis_gap(self.low.z, self.high.z, other.low.z, other.high.z),
118        )
119        .length()
120    }
121
122    fn scale(self, other: Self) -> f64 {
123        self.high
124            .sub(self.low)
125            .length()
126            .max(other.high.sub(other.low).length())
127            .max(1.0)
128    }
129}
130
131fn surface_samples(surface: &NurbsSurface) -> Result<Vec<(f64, f64, Vec3)>, String> {
132    let u = surface.domain_u()?;
133    let v = surface.domain_v()?;
134    let mut samples = Vec::with_capacity(25);
135    for i in 0..=4 {
136        for j in 0..=4 {
137            let uu = u[0] + (u[1] - u[0]) * i as f64 / 4.0;
138            let vv = v[0] + (v[1] - v[0]) * j as f64 / 4.0;
139            samples.push((uu, vv, surface.evaluate(uu, vv)?));
140        }
141    }
142    Ok(samples)
143}
144
145/// Per-surface data the pair classifier needs, precomputable once per face
146/// and reused across every pair that face participates in.
147pub struct SurfaceClassifyData {
148    bounds: Bounds,
149    planar: Option<(Vec3, Vec3)>,
150    /// (u, v, point, normal); normal is None where the surface is singular
151    /// (poles) — the classifier treats close singular samples specially.
152    samples: Vec<(f64, f64, Vec3, Option<Vec3>)>,
153}
154
155thread_local! {
156    static CLASSIFY_MEMO: std::cell::RefCell<rustc_hash::FxHashMap<[u64; 2], std::rc::Rc<SurfaceClassifyData>>> = Default::default();
157    static CLASSIFY_READS: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
158}
159
160pub(crate) fn classify_reads() -> usize { CLASSIFY_READS.with(std::cell::Cell::get) }
161
162pub(crate) fn clear_classify_memo() {
163    CLASSIFY_MEMO.with(|memo| memo.borrow_mut().clear());
164}
165
166impl SurfaceClassifyData {
167    /// Only geometry and tolerance enter this read. Renumbering a face or
168    /// moving it to the next boolean operand does not require new samples.
169    pub(crate) fn memoized(surface: &NurbsSurface, tolerance: f64) -> Result<std::rc::Rc<Self>, String> {
170        let enabled = std::env::var("BREP_CLASSIFY_MEMO").as_deref() != Ok("0");
171        let mut hash = crate::mass_properties::TwoLane::new();
172        crate::mass_properties::hash_surface(&mut hash, surface);
173        hash.f64(tolerance);
174        let key = hash.finish();
175        if enabled {
176            if let Some(data) = CLASSIFY_MEMO.with(|memo| memo.borrow().get(&key).cloned()) {
177                return Ok(data);
178            }
179        }
180        CLASSIFY_READS.with(|reads| reads.set(reads.get() + 1));
181        let data = std::rc::Rc::new(Self::build(surface, tolerance)?);
182        if enabled {
183            CLASSIFY_MEMO.with(|memo| {
184                let mut memo = memo.borrow_mut();
185                if memo.len() >= 4096 { memo.clear(); }
186                memo.insert(key, data.clone());
187            });
188        }
189        Ok(data)
190    }
191
192    pub fn build(surface: &NurbsSurface, model_tolerance: f64) -> Result<Self, String> {
193        let samples = surface_samples(surface)?
194            .into_iter()
195            .map(|(u, v, point)| (u, v, point, surface.normal(u, v).ok()))
196            .collect();
197        Ok(Self {
198            bounds: Bounds::of(surface)?,
199            planar: planar_support(surface, model_tolerance)?,
200            samples,
201        })
202    }
203}
204
205/// Classify a pair before generic SSI runs.  The classifier is conservative:
206/// only separated control hulls are declared disjoint; uncertain pairs stay
207/// `Candidate` rather than being incorrectly culled.
208pub fn classify_surface_pair(
209    first: &NurbsSurface,
210    second: &NurbsSurface,
211    model_tolerance: f64,
212    angular_tolerance: f64,
213) -> Result<SurfacePairClassification, String> {
214    let first_data = SurfaceClassifyData::build(first, model_tolerance)?;
215    let second_data = SurfaceClassifyData::build(second, model_tolerance)?;
216    classify_surface_pair_cached(
217        first,
218        &first_data,
219        second,
220        &second_data,
221        model_tolerance,
222        angular_tolerance,
223    )
224}
225
226/// `classify_surface_pair` with the per-surface work hoisted out — callers
227/// with many pairs per face build each face's `SurfaceClassifyData` once.
228pub fn classify_surface_pair_cached(
229    first: &NurbsSurface,
230    first_data: &SurfaceClassifyData,
231    second: &NurbsSurface,
232    second_data: &SurfaceClassifyData,
233    model_tolerance: f64,
234    angular_tolerance: f64,
235) -> Result<SurfacePairClassification, String> {
236    let first_bounds = first_data.bounds;
237    let second_bounds = second_data.bounds;
238    let hull_separation = first_bounds.separation(second_bounds);
239    let scale = first_bounds.scale(second_bounds);
240    let contact = (model_tolerance * 100.0).max(scale * 1e-7);
241    if hull_separation > contact {
242        return Ok(SurfacePairClassification {
243            relation: SurfacePairRelation::Disjoint,
244            minimum_separation: hull_separation,
245            minimum_normal_cross: 1.0,
246            tangential_only: false,
247        });
248    }
249    if same_surface(first, second, model_tolerance) {
250        return Ok(SurfacePairClassification {
251            relation: SurfacePairRelation::Cosurface,
252            minimum_separation: 0.0,
253            minimum_normal_cross: 0.0,
254            tangential_only: true,
255        });
256    }
257    if let (Some((first_origin, first_normal)), Some((second_origin, second_normal))) =
258        (first_data.planar, second_data.planar)
259    {
260        let normal_cross = first_normal.cross(second_normal).length();
261        if normal_cross <= angular_tolerance {
262            let separation = second_origin.sub(first_origin).dot(first_normal).abs();
263            return Ok(SurfacePairClassification {
264                relation: if separation <= contact {
265                    SurfacePairRelation::Cosurface
266                } else {
267                    SurfacePairRelation::Disjoint
268                },
269                minimum_separation: separation,
270                minimum_normal_cross: normal_cross,
271                tangential_only: separation <= contact,
272            });
273        }
274        return Ok(SurfacePairClassification {
275            relation: SurfacePairRelation::Transverse,
276            minimum_separation: 0.0,
277            minimum_normal_cross: normal_cross,
278            tangential_only: false,
279        });
280    }
281
282    let mut minimum_separation = f64::INFINITY;
283    let mut minimum_normal_cross = f64::INFINITY;
284    let mut close_samples = 0usize;
285    let mut parallel_close_samples = 0usize;
286    let mut singular_close_sample = false;
287    for (source_data, target) in [(first_data, second), (second_data, first)] {
288        for &(_, _, point, source_normal) in &source_data.samples {
289            let projection = project_point_to_surface(target, point)?;
290            minimum_separation = minimum_separation.min(projection.distance);
291            if projection.distance > contact {
292                continue;
293            }
294            close_samples += 1;
295            match (source_normal, target.normal(projection.u, projection.v)) {
296                (Some(first_normal), Ok(second_normal)) => {
297                    let cross = first_normal.cross(second_normal).length();
298                    minimum_normal_cross = minimum_normal_cross.min(cross);
299                    if cross <= angular_tolerance {
300                        parallel_close_samples += 1;
301                    }
302                }
303                _ => singular_close_sample = true,
304            }
305        }
306    }
307    if minimum_separation == f64::INFINITY {
308        minimum_separation = hull_separation;
309    }
310    if minimum_normal_cross == f64::INFINITY {
311        minimum_normal_cross = 1.0;
312    }
313    let relation = if close_samples == 0 {
314        SurfacePairRelation::Candidate
315    } else if singular_close_sample {
316        SurfacePairRelation::Singular
317    } else if close_samples >= 18 && parallel_close_samples == close_samples {
318        SurfacePairRelation::Cosurface
319    } else if minimum_normal_cross <= angular_tolerance {
320        SurfacePairRelation::NearTangent
321    } else {
322        SurfacePairRelation::Transverse
323    };
324    Ok(SurfacePairClassification {
325        relation,
326        minimum_separation,
327        minimum_normal_cross,
328        tangential_only: close_samples > 0 && parallel_close_samples == close_samples,
329    })
330}
331
332
brep_kernel/brep/pair_classification.rs

brep_kernel/brep/
pair_classification.rs
