conspire/geometry/mesh/quality/manifold/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, Tensor},
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])).value());
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])).value());
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.full_contraction(&direction)
53        < ABS_TOL * n1.full_contraction(&n1) * n2.full_contraction(&n2)
54    {
55        return false;
56    }
57    let axis = {
58        let d = [direction[0].abs(), direction[1].abs(), direction[2].abs()];
59        if d[0] >= d[1] && d[0] >= d[2] {
60            0
61        } else if d[1] >= d[2] {
62            1
63        } else {
64            2
65        }
66    };
67    let interval1 = interval(
68        [v[0][axis].value(), v[1][axis].value(), v[2][axis].value()],
69        dv,
70    );
71    let interval2 = interval(
72        [u[0][axis].value(), u[1][axis].value(), u[2][axis].value()],
73        du,
74    );
75    interval1[0] <= interval2[1] && interval2[0] <= interval1[1]
76}
77
78fn interval(projection: [Scalar; N], distance: [Scalar; N]) -> [Scalar; 2] {
79    let (pivot, a, b) = if distance[0] * distance[1] > 0.0 {
80        (2, 0, 1)
81    } else if distance[0] * distance[2] > 0.0 {
82        (1, 0, 2)
83    } else if distance[1] * distance[2] > 0.0 || distance[0] != 0.0 {
84        (0, 1, 2)
85    } else if distance[1] != 0.0 {
86        (1, 0, 2)
87    } else {
88        (2, 0, 1)
89    };
90    let crossing = |q| {
91        projection[pivot]
92            + (projection[q] - projection[pivot]) * distance[pivot]
93                / (distance[pivot] - distance[q])
94    };
95    let (ta, tb) = (crossing(a), crossing(b));
96    if ta <= tb { [ta, tb] } else { [tb, ta] }
97}
conspire/geometry/mesh/quality/manifold/triangles/mod.rs

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