marching_cubes/marching.rs
1use super::container::*;
2use super::tables::*;
3
4/// A struct representing the Marching Cubes algorithm.
5///
6/// The Marching Cubes algorithm is used to extract a triangular mesh representation of a three-dimensional surface from a scalar field.
7///
8/// # Examples
9///
10/// ```
11/// use iso::{MarchingCubes, GridCell, Triangle};
12///
13/// let grid = GridCell {
14/// positions: [
15/// [0.0, 0.0, 0.0],
16/// [1.0, 0.0, 0.0],
17/// [1.0, 1.0, 0.0],
18/// [0.0, 1.0, 0.0],
19/// [0.0, 0.0, 1.0],
20/// [1.0, 0.0, 1.0],
21/// [1.0, 1.0, 1.0],
22/// [0.0, 1.0, 1.0]
23/// ],
24/// value: [0.0, 0.5, 0.5, 1.0, 0.0, 1.0, 1.0, 0.0],
25/// };
26///
27/// let isolevel = 0.5;
28/// let mc = MarchingCubes::new(isolevel, grid);
29/// let mut triangles = vec![];
30/// let triangle_count = mc.polygonise(&mut triangles);
31///
32/// assert_eq!(triangle_count, 4);
33/// ```
34#[derive(Debug, Clone, Copy)]
35pub struct MarchingCubes {
36 iso_value: f32,
37 grid: GridCell,
38}
39
40#[derive(Debug, Clone, Copy)]
41pub struct Triangle {
42 pub positions: [Vec3; 3],
43}
44
45#[derive(Debug, Clone, Copy)]
46pub struct GridCell {
47 pub positions: [Vec3; 8],
48 pub value: [f32; 8],
49}
50
51impl MarchingCubes {
52 /// Creates a new instance of the `MarchingCubes` struct.
53 ///
54 /// # Parameters
55 ///
56 /// * `iso_value`: A float value representing the isolevel to use when determining which parts of the grid cell are inside and outside the surface.
57 /// * `grid`: A `GridCell` struct representing the cell to be polygonised. It contains the positions of the cell's vertices and their values.
58 ///
59 /// # Returns
60 ///
61 /// A new instance of the `MarchingCubes` struct.
62 pub fn new(iso_value: f32, grid: GridCell) -> Self {
63 Self {
64 iso_value,
65 grid
66 }
67 }
68
69 /// Polygonises a grid cell based on the stored isolevel, populating the given `triangles` slice with the resulting triangles.
70 /// Look at [polygonise] for more information regarding this function.
71 ///
72 /// # Parameters
73 ///
74 /// * `triangles`: A mutable slice of `Triangle` structs to be filled with the resulting triangles.
75 ///
76 /// # Returns
77 ///
78 /// The number of triangles generated, as an `i32`.
79 pub fn polygonise(self, triangles: &mut [Triangle]) -> i32 {
80 polygonise(self.grid, self.iso_value, triangles)
81 }
82}
83
84/// Polygonises a grid cell based on the given isolevel, populating the given
85/// `triangles` slice with the resulting triangles.
86///
87/// # Parameters
88///
89/// * `grid_cell`: A `GridCell` struct representing the cell to be polygonised. It contains the positions of the cell's vertices and their values.
90/// * `isolevel`: A float value representing the isolevel to use when determining
91/// which parts of the cell are inside and outside the surface. The isolevel separates the grid cell into parts that are either inside or outside the surface being represented.
92/// * `triangles`: A mutable slice of `Triangle` structs to be filled with the
93/// resulting triangles. The triangles are generated by triangulating the parts of the grid cell that are inside the surface.
94///
95/// # Returns
96///
97/// The number of triangles generated, as an `i32`.
98///
99/// # Example
100///
101/// ```
102/// use iso::{GridCell, Triangle, polygonise};
103///
104/// let grid_cell = GridCell {
105/// positions: [
106/// [0.0, 0.0, 0.0],
107/// [1.0, 0.0, 0.0],
108/// [1.0, 1.0, 0.0],
109/// [0.0, 1.0, 0.0],
110/// [0.0, 0.0, 1.0],
111/// [1.0, 0.0, 1.0],
112/// [1.0, 1.0, 1.0],
113/// [0.0, 1.0, 1.0]
114/// ],
115/// value: [0.0, 0.5, 0.5, 1.0, 0.0, 1.0, 1.0, 0.0],
116/// };
117/// let mut triangles = vec![];
118///
119/// let isolevel = 0.5;
120/// let result = polygonise(grid_cell, isolevel, &mut triangles);
121///
122/// assert_eq!(result, 4);
123/// ```
124pub fn polygonise(grid_cell: GridCell, isolevel: f32, triangles: &mut [Triangle]) -> i32 {
125 let mut cube_index: usize;
126 let mut vertices_list: [Vec3; 12] = [empty_vec3(); 12];
127
128 cube_index = 0;
129
130 if grid_cell.value[0] < isolevel {
131 cube_index |= 1;
132 }
133
134 if grid_cell.value[1] < isolevel {
135 cube_index |= 2;
136 }
137
138 if grid_cell.value[2] < isolevel {
139 cube_index |= 4;
140 }
141 if grid_cell.value[3] < isolevel {
142 cube_index |= 8;
143 }
144 if grid_cell.value[4] < isolevel {
145 cube_index |= 16;
146 }
147 if grid_cell.value[5] < isolevel {
148 cube_index |= 32;
149 }
150 if grid_cell.value[6] < isolevel {
151 cube_index |= 64;
152 }
153 if grid_cell.value[7] < isolevel {
154 cube_index |= 128;
155 }
156
157 if EDGE_TABLE[cube_index] == 0 {
158 return 0i32;
159 }
160
161 if EDGE_TABLE[cube_index] & 1 == 0 {
162 vertices_list[0] = interpolate_vertex(
163 isolevel,
164 grid_cell.positions[0],
165 grid_cell.positions[1],
166 vector2(grid_cell.value[0], grid_cell.value[1]),
167 )
168 }
169
170 if EDGE_TABLE[cube_index] & 2 == 0 {
171 vertices_list[1] = interpolate_vertex(
172 isolevel,
173 grid_cell.positions[1],
174 grid_cell.positions[2],
175 vector2(grid_cell.value[1], grid_cell.value[2]),
176 )
177 }
178
179 if EDGE_TABLE[cube_index] & 4 == 0 {
180 vertices_list[2] = interpolate_vertex(
181 isolevel,
182 grid_cell.positions[2],
183 grid_cell.positions[3],
184 vector2(grid_cell.value[2], grid_cell.value[3]),
185 )
186 }
187
188 if EDGE_TABLE[cube_index] & 8 == 0 {
189 vertices_list[3] = interpolate_vertex(
190 isolevel,
191 grid_cell.positions[3],
192 grid_cell.positions[0],
193 vector2(grid_cell.value[3], grid_cell.value[0]),
194 )
195 }
196
197 if EDGE_TABLE[cube_index] & 16 == 0 {
198 vertices_list[4] = interpolate_vertex(
199 isolevel,
200 grid_cell.positions[4],
201 grid_cell.positions[5],
202 vector2(grid_cell.value[4], grid_cell.value[5]),
203 )
204 }
205
206 if EDGE_TABLE[cube_index] & 32 == 0 {
207 vertices_list[5] = interpolate_vertex(
208 isolevel,
209 grid_cell.positions[5],
210 grid_cell.positions[6],
211 vector2(grid_cell.value[5], grid_cell.value[6]),
212 )
213 }
214
215 if EDGE_TABLE[cube_index] & 64 == 0 {
216 vertices_list[6] = interpolate_vertex(
217 isolevel,
218 grid_cell.positions[6],
219 grid_cell.positions[7],
220 vector2(grid_cell.value[6], grid_cell.value[7]),
221 )
222 }
223
224 if EDGE_TABLE[cube_index] & 128 == 0 {
225 vertices_list[7] = interpolate_vertex(
226 isolevel,
227 grid_cell.positions[7],
228 grid_cell.positions[4],
229 vector2(grid_cell.value[7], grid_cell.value[4]),
230 )
231 }
232
233 if EDGE_TABLE[cube_index] & 256 == 0 {
234 vertices_list[8] = interpolate_vertex(
235 isolevel,
236 grid_cell.positions[0],
237 grid_cell.positions[4],
238 vector2(grid_cell.value[0], grid_cell.value[4]),
239 )
240 }
241 if EDGE_TABLE[cube_index] & 512 == 0 {
242 vertices_list[9] = interpolate_vertex(
243 isolevel,
244 grid_cell.positions[1],
245 grid_cell.positions[5],
246 vector2(grid_cell.value[1], grid_cell.value[5]),
247 )
248 }
249
250 if EDGE_TABLE[cube_index] & 1024 == 0 {
251 vertices_list[10] = interpolate_vertex(
252 isolevel,
253 grid_cell.positions[2],
254 grid_cell.positions[6],
255 vector2(grid_cell.value[2], grid_cell.value[6]),
256 )
257 }
258
259 if EDGE_TABLE[cube_index] & 2048 == 0 {
260 vertices_list[11] = interpolate_vertex(
261 isolevel,
262 grid_cell.positions[3],
263 grid_cell.positions[7],
264 vector2(grid_cell.value[3], grid_cell.value[7]),
265 )
266 }
267
268 let mut triangle_num = 0;
269
270 for mut i in 0.. {
271 let tri = TRI_TABLE[cube_index][i];
272
273 if tri == -1 {
274 break;
275 }
276
277 i += 3;
278
279 (*triangles)[triangle_num].positions[0] = vertices_list[TRI_TABLE[cube_index][i] as usize];
280 (*triangles)[triangle_num].positions[1] =
281 vertices_list[TRI_TABLE[cube_index][i + 1] as usize];
282 (*triangles)[triangle_num].positions[2] =
283 vertices_list[TRI_TABLE[cube_index][i + 2] as usize];
284
285 triangle_num += 1;
286 }
287
288 0i32
289}
290
291/// Interpolates a vertex using `isolevel`, `point1`, `point2`, and `alpha_points`
292///
293/// The `isolevel` is used to determine the position of the vertex between `point1` and `point2`.
294/// `alpha_points` is used to determine the relative values of `isolevel` and `point1` and `point2`.
295/// `point1` and `point2` represent the two points that the vertex will be interpolated between.
296///
297/// # Parameters
298///
299/// * `isolevel` - The level of iso to interpolate the vertex at
300/// * `point1` - The first point that the vertex will be interpolated between
301/// * `point2` - The second point that the vertex will be interpolated between
302/// * `alpha_points` - A Vec2 that contains the relative values of `isolevel` and `point1` and `point2`
303///
304/// # Returns
305///
306/// A `Vec3` representing the interpolated vertex
307///
308/// # Example
309///
310/// ```rust
311/// let point1 = vec3(0.0, 0.0, 0.0);
312/// let point2 = vec3(1.0, 1.0, 1.0);
313/// let alpha_points = vec2(0.5, 0.5);
314/// let isolevel = 0.5;
315///
316/// let result = interpolate_vertex(isolevel, point1, point2, alpha_points);
317///
318/// assert_eq!(result, vec3(0.5, 0.5, 0.5));
319/// ```
320fn interpolate_vertex(isolevel: f32, point1: Vec3, point2: Vec3, alpha_points: Vec2) -> Vec3 {
321 const ISO_THRESHOLD: f32 = 0.00001;
322
323 let mut point = empty_vec3();
324 let factor = (isolevel - alpha_points[0]) / (alpha_points[1] - alpha_points[0]);
325
326 if (isolevel - alpha_points[0]).abs() < ISO_THRESHOLD
327 || (alpha_points[0] - alpha_points[1]) < ISO_THRESHOLD
328 {
329 return point1;
330 } else if (isolevel - alpha_points[1]).abs() < ISO_THRESHOLD {
331 return point2;
332 }
333
334 point[0] = point1[0] + factor * (point2[0] - point1[0]);
335 point[1] = point1[1] + factor * (point2[1] - point1[1]);
336 point[2] = point1[2] + factor * (point2[2] - point1[2]);
337 point
338}