lib3mf_core/validation/
bvh.rs

1use crate::model::{Mesh, Triangle};
2use glam::Vec3;
3
4#[derive(Debug, Clone, Copy)]
5pub struct AABB {
6    pub min: Vec3,
7    pub max: Vec3,
8}
9
10impl AABB {
11    pub fn from_triangle(mesh: &Mesh, tri: &Triangle) -> Self {
12        let v1 = mesh.vertices[tri.v1 as usize];
13        let v2 = mesh.vertices[tri.v2 as usize];
14        let v3 = mesh.vertices[tri.v3 as usize];
15
16        let min = Vec3::new(
17            v1.x.min(v2.x).min(v3.x),
18            v1.y.min(v2.y).min(v3.y),
19            v1.z.min(v2.z).min(v3.z),
20        );
21        let max = Vec3::new(
22            v1.x.max(v2.x).max(v3.x),
23            v1.y.max(v2.y).max(v3.y),
24            v1.z.max(v2.z).max(v3.z),
25        );
26        Self { min, max }
27    }
28
29    pub fn intersects(&self, other: &Self) -> bool {
30        self.min.x <= other.max.x
31            && self.max.x >= other.min.x
32            && self.min.y <= other.max.y
33            && self.max.y >= other.min.y
34            && self.min.z <= other.max.z
35            && self.max.z >= other.min.z
36    }
37}
38
39pub struct BvhNode {
40    pub aabb: AABB,
41    pub content: BvhContent,
42}
43
44pub enum BvhContent {
45    Leaf(Vec<usize>), // Indices of triangles
46    Branch(Box<BvhNode>, Box<BvhNode>),
47}
48
49impl BvhNode {
50    pub fn build(mesh: &Mesh, tri_indices: Vec<usize>) -> Self {
51        let aabbs: Vec<AABB> = tri_indices
52            .iter()
53            .map(|&i| AABB::from_triangle(mesh, &mesh.triangles[i]))
54            .collect();
55
56        let mut min = Vec3::splat(f32::INFINITY);
57        let mut max = Vec3::splat(f32::NEG_INFINITY);
58        for aabb in &aabbs {
59            min = min.min(aabb.min);
60            max = max.max(aabb.max);
61        }
62
63        let node_aabb = AABB { min, max };
64
65        if tri_indices.len() <= 8 {
66            return BvhNode {
67                aabb: node_aabb,
68                content: BvhContent::Leaf(tri_indices),
69            };
70        }
71
72        // Split along largest axis
73        let size = max - min;
74        let axis = if size.x > size.y && size.x > size.z {
75            0
76        } else if size.y > size.z {
77            1
78        } else {
79            2
80        };
81
82        let mid = (min[axis] + max[axis]) / 2.0;
83
84        let mut left_indices = Vec::new();
85        let mut right_indices = Vec::new();
86
87        for (i, &tri_idx) in tri_indices.iter().enumerate() {
88            let center = (aabbs[i].min[axis] + aabbs[i].max[axis]) / 2.0;
89            if center < mid {
90                left_indices.push(tri_idx);
91            } else {
92                right_indices.push(tri_idx);
93            }
94        }
95
96        // Fallback if split failed to partition
97        if left_indices.is_empty() || right_indices.is_empty() {
98            return BvhNode {
99                aabb: node_aabb,
100                content: BvhContent::Leaf(tri_indices),
101            };
102        }
103
104        BvhNode {
105            aabb: node_aabb,
106            content: BvhContent::Branch(
107                Box::new(BvhNode::build(mesh, left_indices)),
108                Box::new(BvhNode::build(mesh, right_indices)),
109            ),
110        }
111    }
112
113    pub fn find_intersections(
114        &self,
115        mesh: &Mesh,
116        tri_idx: usize,
117        tri_aabb: &AABB,
118        results: &mut Vec<usize>,
119    ) {
120        if !self.aabb.intersects(tri_aabb) {
121            return;
122        }
123
124        match &self.content {
125            BvhContent::Leaf(indices) => {
126                for &idx in indices {
127                    if idx > tri_idx {
128                        // Avoid double-counting and self-check
129                        if tri_aabb.intersects(&AABB::from_triangle(mesh, &mesh.triangles[idx])) {
130                            // Precise check:
131                            if intersect_triangles(mesh, tri_idx, idx) {
132                                results.push(idx);
133                            }
134                        }
135                    }
136                }
137            }
138            BvhContent::Branch(left, right) => {
139                left.find_intersections(mesh, tri_idx, tri_aabb, results);
140                right.find_intersections(mesh, tri_idx, tri_aabb, results);
141            }
142        }
143    }
144}
145
146/// Robust triangle-triangle intersection (simplified Moller-Trumbore)
147fn intersect_triangles(mesh: &Mesh, i1: usize, i2: usize) -> bool {
148    let t1 = &mesh.triangles[i1];
149    let t2 = &mesh.triangles[i2];
150
151    // Shared vertices: check if they share 1, 2 or 3 vertices.
152    // If they share 2 vertices, they share an edge.
153    // If they share 3, they are identical.
154    // In many contexts, edge/vert sharing is NOT considered intersection for "self-intersection" reporting.
155    // 3MF spec: triangles MUST NOT intersect each other EXCEPT on common edges.
156    let shared = count_shared_vertices(t1, t2);
157    if shared >= 2 {
158        return false;
159    }
160
161    let p1 = to_vec3(mesh.vertices[t1.v1 as usize]);
162    let p2 = to_vec3(mesh.vertices[t1.v2 as usize]);
163    let p3 = to_vec3(mesh.vertices[t1.v3 as usize]);
164
165    let q1 = to_vec3(mesh.vertices[t2.v1 as usize]);
166    let q2 = to_vec3(mesh.vertices[t2.v2 as usize]);
167    let q3 = to_vec3(mesh.vertices[t2.v3 as usize]);
168
169    tri_tri_intersect(p1, p2, p3, q1, q2, q3)
170}
171
172fn to_vec3(v: crate::model::Vertex) -> Vec3 {
173    Vec3::new(v.x, v.y, v.z)
174}
175
176fn count_shared_vertices(t1: &Triangle, t2: &Triangle) -> usize {
177    let mut count = 0;
178    let v1 = [t1.v1, t1.v2, t1.v3];
179    let v2 = [t2.v1, t2.v2, t2.v3];
180    for &va in &v1 {
181        for &vb in &v2 {
182            if va == vb {
183                count += 1;
184            }
185        }
186    }
187    count
188}
189
190// ----------------------------------------------------------------------------
191// Simplified Triangle-Triangle Intersection logic (adapted)
192// ----------------------------------------------------------------------------
193
194fn tri_tri_intersect(p1: Vec3, p2: Vec3, p3: Vec3, q1: Vec3, q2: Vec3, q3: Vec3) -> bool {
195    // 1. Plane of tri 2
196    let n2 = (q2 - q1).cross(q3 - q1);
197    if n2.length_squared() < 1e-12 {
198        return false;
199    } // Degenerate tri 2
200    let d2 = -n2.dot(q1);
201
202    // Distances of tri 1 vertices to plane 2
203    let du0 = n2.dot(p1) + d2;
204    let du1 = n2.dot(p2) + d2;
205    let du2 = n2.dot(p3) + d2;
206
207    if (du0.abs() > 1e-6 && du1.abs() > 1e-6 && du2.abs() > 1e-6)
208        && ((du0 > 0.0 && du1 > 0.0 && du2 > 0.0) || (du0 < 0.0 && du1 < 0.0 && du2 < 0.0))
209    {
210        return false; // Tri 1 entirely on one side of plane 2
211    }
212
213    // 2. Plane of tri 1
214    let n1 = (p2 - p1).cross(p3 - p1);
215    if n1.length_squared() < 1e-12 {
216        return false;
217    } // Degenerate tri 1
218    let d1 = -n1.dot(p1);
219
220    // Distances of tri 2 vertices to plane 1
221    let dv0 = n1.dot(q1) + d1;
222    let dv1 = n1.dot(q2) + d1;
223    let dv2 = n1.dot(q3) + d1;
224
225    if (dv0.abs() > 1e-6 && dv1.abs() > 1e-6 && dv2.abs() > 1e-6)
226        && ((dv0 > 0.0 && dv1 > 0.0 && dv2 > 0.0) || (dv0 < 0.0 && dv1 < 0.0 && dv2 < 0.0))
227    {
228        return false; // Tri 2 entirely on one side of plane 1
229    }
230
231    // 3. Line of intersection L
232    let ld = n1.cross(n2);
233    let index = if ld.x.abs() > ld.y.abs() && ld.x.abs() > ld.z.abs() {
234        0
235    } else if ld.y.abs() > ld.z.abs() {
236        1
237    } else {
238        2
239    };
240
241    // Projection onto L
242    let get_interval =
243        |v1: Vec3, v2: Vec3, v3: Vec3, d1: f32, d2: f32, d3: f32| -> Option<(f32, f32)> {
244            // Find which vertices are on which side
245            if (d1 > 0.0 && d2 > 0.0 && d3 > 0.0) || (d1 < 0.0 && d2 < 0.0 && d3 < 0.0) {
246                return None;
247            }
248
249            let mut pts = Vec::new();
250            let tris = [(v1, v2, d1, d2), (v2, v3, d2, d3), (v3, v1, d3, d1)];
251            for (a, b, da, db) in tris {
252                if (da >= 0.0) != (db >= 0.0) {
253                    let t = da / (da - db);
254                    let p = a + (b - a) * t;
255                    pts.push(p[index]);
256                } else if da.abs() < 1e-7 {
257                    pts.push(a[index]);
258                }
259            }
260            if pts.len() < 2 {
261                return None;
262            }
263            let mut min = pts[0];
264            let mut max = pts[0];
265            for &p in &pts {
266                min = min.min(p);
267                max = max.max(p);
268            }
269            Some((min, max))
270        };
271
272    let i1 = get_interval(p1, p2, p3, du0, du1, du2);
273    let i2 = get_interval(q1, q2, q3, dv0, dv1, dv2);
274
275    match (i1, i2) {
276        (Some((t1_min, t1_max)), Some((t2_min, t2_max))) => {
277            // Check overlap of intervals on line L
278            // Use small epsilon to avoid reporting touching as intersection
279            t1_min + 1e-6 < t2_max && t2_min + 1e-6 < t1_max
280        }
281        _ => false,
282    }
283}