conspire/geometry/mesh/quality/surface/triangles/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    ABS_TOL,
6    geometry::{Coordinates, bvh::BoundingVolumeHierarchy, mesh::Mesh},
7    math::{CrossProduct, Scalar},
8};
9
10const D: usize = 3;
11const N: usize = 3;
12
13impl Mesh<D> {
14    pub fn self_intersections(&self) -> Vec<[usize; 2]> {
15        let coordinates = self.coordinates();
16        let faces: Vec<[usize; N]> = self
17            .iter()
18            .flat_map(|block| block.iter())
19            .map(|element| [element[0], element[1], element[2]])
20            .collect();
21        let boxes = self.bounding_boxes();
22        let bvh = BoundingVolumeHierarchy::from(self);
23        let mut hits = Vec::new();
24        for (i, face) in faces.iter().enumerate() {
25            for j in bvh.overlapping(&boxes[i]) {
26                if j > i
27                    && !face.iter().any(|node| faces[j].contains(node))
28                    && triangles_intersect(*face, faces[j], coordinates)
29                {
30                    hits.push([i, j]);
31                }
32            }
33        }
34        hits
35    }
36}
37
38fn triangles_intersect(t1: [usize; N], t2: [usize; N], coordinates: &Coordinates<D>) -> bool {
39    let v = t1.map(|i| &coordinates[i]);
40    let u = t2.map(|i| &coordinates[i]);
41    let n1 = (v[1] - v[0]).cross(v[2] - v[0]);
42    let du = u.map(|p| &n1 * &(p - v[0]));
43    if du[0] * du[1] > 0.0 && du[0] * du[2] > 0.0 {
44        return false;
45    }
46    let n2 = (u[1] - u[0]).cross(u[2] - u[0]);
47    let dv = v.map(|p| &n2 * &(p - u[0]));
48    if dv[0] * dv[1] > 0.0 && dv[0] * dv[2] > 0.0 {
49        return false;
50    }
51    let direction = n1.cross(&n2);
52    if &direction * &direction < ABS_TOL * (&n1 * &n1) * (&n2 * &n2) {
53        return false; // coplanar (or degenerate)
54    }
55    let axis = {
56        let d = [direction[0].abs(), direction[1].abs(), direction[2].abs()];
57        if d[0] >= d[1] && d[0] >= d[2] {
58            0
59        } else if d[1] >= d[2] {
60            1
61        } else {
62            2
63        }
64    };
65    let interval1 = interval([v[0][axis], v[1][axis], v[2][axis]], dv);
66    let interval2 = interval([u[0][axis], u[1][axis], u[2][axis]], du);
67    interval1[0] <= interval2[1] && interval2[0] <= interval1[1]
68}
69
70fn interval(projection: [Scalar; N], distance: [Scalar; N]) -> [Scalar; 2] {
71    let (pivot, a, b) = if distance[0] * distance[1] > 0.0 {
72        (2, 0, 1)
73    } else if distance[0] * distance[2] > 0.0 {
74        (1, 0, 2)
75    } else if distance[1] * distance[2] > 0.0 || distance[0] != 0.0 {
76        (0, 1, 2)
77    } else if distance[1] != 0.0 {
78        (1, 0, 2)
79    } else {
80        (2, 0, 1)
81    };
82    let crossing = |q| {
83        projection[pivot]
84            + (projection[q] - projection[pivot]) * distance[pivot]
85                / (distance[pivot] - distance[q])
86    };
87    let (ta, tb) = (crossing(a), crossing(b));
88    if ta <= tb { [ta, tb] } else { [tb, ta] }
89}
conspire/geometry/mesh/quality/surface/triangles/mod.rs

conspire/geometry/mesh/quality/surface/triangles/
mod.rs
