1#[cfg(test)]
2mod test;
3
4use crate::{
5 geometry::{
6 Coordinate, Coordinates,
7 grid::Voxels,
8 mesh::{Connectivity, Mesh, tessellation::Tessellation},
9 },
10 math::{Tensor, TensorVec},
11};
12use std::{
13 array::from_fn,
14 collections::{HashMap, HashSet},
15};
16
17const DELTAS: [[isize; 3]; 6] = [
18 [1, 0, 0],
19 [-1, 0, 0],
20 [0, 1, 0],
21 [0, -1, 0],
22 [0, 0, 1],
23 [0, 0, -1],
24];
25
26const FACES: [[[usize; 3]; 4]; 6] = [
27 [[1, 0, 0], [1, 1, 0], [1, 1, 1], [1, 0, 1]],
28 [[0, 0, 0], [0, 0, 1], [0, 1, 1], [0, 1, 0]],
29 [[0, 1, 0], [0, 1, 1], [1, 1, 1], [1, 1, 0]],
30 [[0, 0, 0], [1, 0, 0], [1, 0, 1], [0, 0, 1]],
31 [[0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]],
32 [[0, 0, 0], [0, 1, 0], [1, 1, 0], [1, 0, 0]],
33];
34
35const SECTORS: [[isize; 2]; 4] = [[0, 0], [-1, 0], [-1, -1], [0, -1]];
36const NEXT: [(usize, bool); 4] = [(0, false), (1, false), (0, true), (1, true)];
37const PREV: [(usize, bool); 4] = [(1, false), (0, true), (1, true), (0, false)];
38
39impl<T> From<Voxels<T>> for Tessellation
40where
41 T: Copy + Default + Ord,
42{
43 fn from(voxels: Voxels<T>) -> Self {
44 let nel = *voxels.nel();
45 let [nx, ny, _] = nel;
46 let void = T::default();
47 let data = voxels.data();
48 let cell = |idx: [usize; 3]| voxels.flat(idx);
49 let label = |idx: [isize; 3]| -> T {
50 if (0..3).any(|ax| idx[ax] < 0 || idx[ax] >= nel[ax] as isize) {
51 void
52 } else {
53 data[voxels.flat([idx[0] as usize, idx[1] as usize, idx[2] as usize])]
54 }
55 };
56 let mut quads = Vec::new();
57 let mut facemap = HashMap::new();
58 let mut edges = HashSet::new();
59 for k in 0..nel[2] {
60 for j in 0..ny {
61 for i in 0..nx {
62 let a = data[cell([i, j, k])];
63 if a == void {
64 continue;
65 }
66 for (face, delta) in DELTAS.iter().enumerate() {
67 let b = label([
68 i as isize + delta[0],
69 j as isize + delta[1],
70 k as isize + delta[2],
71 ]);
72 if a != b {
73 let corners =
74 FACES[face].map(|off| [i + off[0], j + off[1], k + off[2]]);
75 for edge in 0..4 {
76 let (p, q) = (corners[edge], corners[(edge + 1) % 4]);
77 let axis = (0..3).find(|&ax| p[ax] != q[ax]).unwrap();
78 edges.insert((from_fn(|ax| p[ax].min(q[ax])), axis));
79 }
80 facemap.insert((cell([i, j, k]), face), quads.len());
81 quads.push(corners);
82 }
83 }
84 }
85 }
86 }
87 let mut parent: Vec<usize> = (0..4 * quads.len()).collect();
88 for &(p, t) in &edges {
89 let (u, v) = ((t + 1) % 3, (t + 2) % 3);
90 let mut labels = [void; 4];
91 let mut cells = [0usize; 4];
92 for (sector, offset) in SECTORS.iter().enumerate() {
93 let mut idx = [0isize; 3];
94 idx[t] = p[t] as isize;
95 idx[u] = p[u] as isize + offset[0];
96 idx[v] = p[v] as isize + offset[1];
97 labels[sector] = label(idx);
98 if labels[sector] != void {
99 cells[sector] = cell([idx[0] as usize, idx[1] as usize, idx[2] as usize]);
100 }
101 }
102 for start in 0..4 {
103 if labels[start] == void || labels[start] == labels[(start + 3) % 4] {
104 continue;
105 }
106 let mut finish = start;
107 while labels[(finish + 1) % 4] == labels[start] && (finish + 1) % 4 != start {
108 finish = (finish + 1) % 4;
109 }
110 let face =
111 |(local, positive): (usize, bool)| [u, v][local] * 2 + usize::from(!positive);
112 let q1 = facemap[&(cells[start], face(PREV[start]))];
113 let q2 = facemap[&(cells[finish], face(NEXT[finish]))];
114 for offset in [0, 1] {
115 let mut endpoint = p;
116 endpoint[t] += offset;
117 let a = 4 * q1 + local(&quads[q1], endpoint);
118 let b = 4 * q2 + local(&quads[q2], endpoint);
119 union(&mut parent, a, b);
120 }
121 }
122 }
123 let mut coordinates = Coordinates::new();
124 let mut remap = HashMap::new();
125 let mut triangles = Vec::with_capacity(2 * quads.len());
126 for (q, corners) in quads.iter().enumerate() {
127 let ids: [usize; 4] = from_fn(|l| {
128 let root = find(&mut parent, 4 * q + l);
129 *remap.entry(root).or_insert_with(|| {
130 let id = coordinates.len();
131 coordinates.push(Coordinate::const_from(from_fn(|ax| corners[l][ax] as f64)));
132 id
133 })
134 });
135 triangles.push([ids[0], ids[1], ids[2]]);
136 triangles.push([ids[0], ids[2], ids[3]]);
137 }
138 let triangles = resolve_pinches(&triangles, &mut coordinates);
139 let connectivities = vec![Connectivity::Triangular(triangles.into())];
140 Tessellation::from(Mesh::from((connectivities, coordinates)))
141 }
142}
143
144fn resolve_pinches(triangles: &[[usize; 3]], coordinates: &mut Coordinates<3>) -> Vec<[usize; 3]> {
145 let mut vertex_faces: HashMap<usize, Vec<usize>> = HashMap::new();
146 let mut edge_faces: HashMap<[usize; 2], Vec<usize>> = HashMap::new();
147 for (f, tri) in triangles.iter().enumerate() {
148 for i in 0..3 {
149 let (v, w) = (tri[i], tri[(i + 1) % 3]);
150 vertex_faces.entry(v).or_default().push(f);
151 edge_faces.entry(edge(v, w)).or_default().push(f);
152 }
153 }
154 let mut assignment = triangles.to_vec();
155 for (&v, incident) in &vertex_faces {
156 let index: HashMap<usize, usize> =
157 incident.iter().enumerate().map(|(i, &f)| (f, i)).collect();
158 let mut parent: Vec<usize> = (0..incident.len()).collect();
159 for &f in incident {
160 let tri = &triangles[f];
161 let pos = tri.iter().position(|&n| n == v).unwrap();
162 for x in [tri[(pos + 2) % 3], tri[(pos + 1) % 3]] {
163 if let Some(shared) = edge_faces.get(&edge(v, x))
164 && shared.len() == 2
165 {
166 union(&mut parent, index[&shared[0]], index[&shared[1]]);
167 }
168 }
169 }
170 let mut fan: HashMap<usize, usize> = HashMap::new();
171 for i in 0..incident.len() {
172 let root = find(&mut parent, i);
173 let next = fan.len();
174 fan.entry(root).or_insert(next);
175 }
176 if fan.len() < 2 {
177 continue;
178 }
179 let mut copies: HashMap<usize, usize> = HashMap::new();
180 for (i, &f) in incident.iter().enumerate() {
181 let root = find(&mut parent, i);
182 if fan[&root] == 0 {
183 continue;
184 }
185 let id = *copies.entry(root).or_insert_with(|| {
186 let point = Coordinate::const_from(from_fn(|ax| coordinates[v][ax]));
187 let id = coordinates.len();
188 coordinates.push(point);
189 id
190 });
191 let pos = triangles[f].iter().position(|&n| n == v).unwrap();
192 assignment[f][pos] = id;
193 }
194 }
195 assignment
196}
197
198fn edge(a: usize, b: usize) -> [usize; 2] {
199 if a < b { [a, b] } else { [b, a] }
200}
201
202fn local(corners: &[[usize; 3]; 4], corner: [usize; 3]) -> usize {
203 corners.iter().position(|&c| c == corner).unwrap()
204}
205
206fn find(parent: &mut [usize], mut i: usize) -> usize {
207 while parent[i] != i {
208 parent[i] = parent[parent[i]];
209 i = parent[i];
210 }
211 i
212}
213
214fn union(parent: &mut [usize], a: usize, b: usize) {
215 let (ra, rb) = (find(parent, a), find(parent, b));
216 if ra != rb {
217 parent[ra] = rb;
218 }
219}