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
155impl SurfaceClassifyData {
156    pub fn build(surface: &NurbsSurface, model_tolerance: f64) -> Result<Self, String> {
157        let samples = surface_samples(surface)?
158            .into_iter()
159            .map(|(u, v, point)| (u, v, point, surface.normal(u, v).ok()))
160            .collect();
161        Ok(Self {
162            bounds: Bounds::of(surface)?,
163            planar: planar_support(surface, model_tolerance)?,
164            samples,
165        })
166    }
167}
168
169/// Classify a pair before generic SSI runs.  The classifier is conservative:
170/// only separated control hulls are declared disjoint; uncertain pairs stay
171/// `Candidate` rather than being incorrectly culled.
172pub fn classify_surface_pair(
173    first: &NurbsSurface,
174    second: &NurbsSurface,
175    model_tolerance: f64,
176    angular_tolerance: f64,
177) -> Result<SurfacePairClassification, String> {
178    let first_data = SurfaceClassifyData::build(first, model_tolerance)?;
179    let second_data = SurfaceClassifyData::build(second, model_tolerance)?;
180    classify_surface_pair_cached(
181        first,
182        &first_data,
183        second,
184        &second_data,
185        model_tolerance,
186        angular_tolerance,
187    )
188}
189
190/// `classify_surface_pair` with the per-surface work hoisted out — callers
191/// with many pairs per face build each face's `SurfaceClassifyData` once.
192pub fn classify_surface_pair_cached(
193    first: &NurbsSurface,
194    first_data: &SurfaceClassifyData,
195    second: &NurbsSurface,
196    second_data: &SurfaceClassifyData,
197    model_tolerance: f64,
198    angular_tolerance: f64,
199) -> Result<SurfacePairClassification, String> {
200    let first_bounds = first_data.bounds;
201    let second_bounds = second_data.bounds;
202    let hull_separation = first_bounds.separation(second_bounds);
203    let scale = first_bounds.scale(second_bounds);
204    let contact = (model_tolerance * 100.0).max(scale * 1e-7);
205    if hull_separation > contact {
206        return Ok(SurfacePairClassification {
207            relation: SurfacePairRelation::Disjoint,
208            minimum_separation: hull_separation,
209            minimum_normal_cross: 1.0,
210            tangential_only: false,
211        });
212    }
213    if same_surface(first, second, model_tolerance) {
214        return Ok(SurfacePairClassification {
215            relation: SurfacePairRelation::Cosurface,
216            minimum_separation: 0.0,
217            minimum_normal_cross: 0.0,
218            tangential_only: true,
219        });
220    }
221    if let (Some((first_origin, first_normal)), Some((second_origin, second_normal))) =
222        (first_data.planar, second_data.planar)
223    {
224        let normal_cross = first_normal.cross(second_normal).length();
225        if normal_cross <= angular_tolerance {
226            let separation = second_origin.sub(first_origin).dot(first_normal).abs();
227            return Ok(SurfacePairClassification {
228                relation: if separation <= contact {
229                    SurfacePairRelation::Cosurface
230                } else {
231                    SurfacePairRelation::Disjoint
232                },
233                minimum_separation: separation,
234                minimum_normal_cross: normal_cross,
235                tangential_only: separation <= contact,
236            });
237        }
238        return Ok(SurfacePairClassification {
239            relation: SurfacePairRelation::Transverse,
240            minimum_separation: 0.0,
241            minimum_normal_cross: normal_cross,
242            tangential_only: false,
243        });
244    }
245
246    let mut minimum_separation = f64::INFINITY;
247    let mut minimum_normal_cross = f64::INFINITY;
248    let mut close_samples = 0usize;
249    let mut parallel_close_samples = 0usize;
250    let mut singular_close_sample = false;
251    for (source_data, target) in [(first_data, second), (second_data, first)] {
252        for &(_, _, point, source_normal) in &source_data.samples {
253            let projection = project_point_to_surface(target, point)?;
254            minimum_separation = minimum_separation.min(projection.distance);
255            if projection.distance > contact {
256                continue;
257            }
258            close_samples += 1;
259            match (source_normal, target.normal(projection.u, projection.v)) {
260                (Some(first_normal), Ok(second_normal)) => {
261                    let cross = first_normal.cross(second_normal).length();
262                    minimum_normal_cross = minimum_normal_cross.min(cross);
263                    if cross <= angular_tolerance {
264                        parallel_close_samples += 1;
265                    }
266                }
267                _ => singular_close_sample = true,
268            }
269        }
270    }
271    if minimum_separation == f64::INFINITY {
272        minimum_separation = hull_separation;
273    }
274    if minimum_normal_cross == f64::INFINITY {
275        minimum_normal_cross = 1.0;
276    }
277    let relation = if close_samples == 0 {
278        SurfacePairRelation::Candidate
279    } else if singular_close_sample {
280        SurfacePairRelation::Singular
281    } else if close_samples >= 18 && parallel_close_samples == close_samples {
282        SurfacePairRelation::Cosurface
283    } else if minimum_normal_cross <= angular_tolerance {
284        SurfacePairRelation::NearTangent
285    } else {
286        SurfacePairRelation::Transverse
287    };
288    Ok(SurfacePairClassification {
289        relation,
290        minimum_separation,
291        minimum_normal_cross,
292        tangential_only: close_samples > 0 && parallel_close_samples == close_samples,
293    })
294}
295
296#[cfg(test)]
297mod tests {
298    use super::*;
299    use crate::{make_plane, Vec3};
300
301    #[test]
302    fn parallel_separated_planes_are_disjoint() {
303        let first = make_plane(
304            Vec3::default(),
305            Vec3::new(1.0, 0.0, 0.0),
306            Vec3::new(0.0, 1.0, 0.0),
307            2.0,
308            2.0,
309        )
310        .unwrap();
311        let second = make_plane(
312            Vec3::new(0.0, 0.0, 1.0),
313            Vec3::new(1.0, 0.0, 0.0),
314            Vec3::new(0.0, 1.0, 0.0),
315            2.0,
316            2.0,
317        )
318        .unwrap();
319        let class = classify_surface_pair(&first, &second, 1e-7, 1e-4).unwrap();
320        assert_eq!(class.relation, SurfacePairRelation::Disjoint);
321    }
322
323    #[test]
324    fn crossing_planes_are_transverse() {
325        let first = make_plane(
326            Vec3::default(),
327            Vec3::new(1.0, 0.0, 0.0),
328            Vec3::new(0.0, 1.0, 0.0),
329            2.0,
330            2.0,
331        )
332        .unwrap();
333        let second = make_plane(
334            Vec3::default(),
335            Vec3::new(1.0, 0.0, 0.0),
336            Vec3::new(0.0, 0.0, 1.0),
337            2.0,
338            2.0,
339        )
340        .unwrap();
341        let class = classify_surface_pair(&first, &second, 1e-7, 1e-4).unwrap();
342        assert_eq!(class.relation, SurfacePairRelation::Transverse);
343    }
344}
brep_kernel/brep/pair_classification.rs

brep_kernel/brep/
pair_classification.rs
