viewport_lib/geometry/
marching_cubes.rs

1//! CPU-side marching cubes isosurface extraction from volumetric scalar data.
2//!
3//! The output is a standard [`MeshData`](crate::resources::MeshData) that can be uploaded via
4//! [`upload_mesh_data()`](crate::ViewportGpuResources::upload_mesh_data) or
5//! [`replace_mesh_data()`](crate::ViewportGpuResources::replace_mesh_data).
6//!
7//! # Example
8//!
9//! ```ignore
10//! let volume = VolumeData {
11//!     data: vec![/* scalar values */],
12//!     dims: [32, 32, 32],
13//!     origin: [0.0, 0.0, 0.0],
14//!     spacing: [0.1, 0.1, 0.1],
15//! };
16//! let mesh = extract_isosurface(&volume, 0.5);
17//! // mesh.positions, mesh.normals, mesh.indices ready for upload.
18//! ```
19
20use crate::resources::MeshData;
21use std::collections::HashMap;
22
23/// A structured 3D scalar field on a regular grid.
24#[derive(Debug, Clone)]
25pub struct VolumeData {
26    /// Flattened scalar values in x-fastest order: `index = x + y*nx + z*nx*ny`.
27    pub data: Vec<f32>,
28    /// Grid dimensions `[nx, ny, nz]`.
29    pub dims: [u32; 3],
30    /// World-space origin of the grid corner `(0, 0, 0)`.
31    pub origin: [f32; 3],
32    /// Cell size in each axis direction.
33    pub spacing: [f32; 3],
34}
35
36impl VolumeData {
37    /// Check whether grid indices are within bounds.
38    pub fn in_bounds(&self, ix: u32, iy: u32, iz: u32) -> bool {
39        ix < self.dims[0] && iy < self.dims[1] && iz < self.dims[2]
40    }
41
42    /// Read the scalar value at grid point `(ix, iy, iz)`.
43    ///
44    /// # Panics
45    ///
46    /// Panics if the index is out of bounds.
47    pub fn sample(&self, ix: u32, iy: u32, iz: u32) -> f32 {
48        let nx = self.dims[0] as usize;
49        let ny = self.dims[1] as usize;
50        self.data[ix as usize + iy as usize * nx + iz as usize * nx * ny]
51    }
52}
53
54/// Trilinear interpolation of the volume at an arbitrary world-space position.
55///
56/// Returns the interpolated scalar value. Points outside the grid are clamped
57/// to the boundary.
58pub fn trilinear_sample(volume: &VolumeData, world_pos: [f32; 3]) -> f32 {
59    let [nx, ny, nz] = volume.dims;
60
61    // Convert world position to continuous grid coordinates.
62    let gx = (world_pos[0] - volume.origin[0]) / volume.spacing[0];
63    let gy = (world_pos[1] - volume.origin[1]) / volume.spacing[1];
64    let gz = (world_pos[2] - volume.origin[2]) / volume.spacing[2];
65
66    // Clamp to valid range.
67    let gx = gx.clamp(0.0, (nx as f32) - 1.001);
68    let gy = gy.clamp(0.0, (ny as f32) - 1.001);
69    let gz = gz.clamp(0.0, (nz as f32) - 1.001);
70
71    let ix = gx.floor() as u32;
72    let iy = gy.floor() as u32;
73    let iz = gz.floor() as u32;
74
75    let fx = gx - ix as f32;
76    let fy = gy - iy as f32;
77    let fz = gz - iz as f32;
78
79    let ix1 = (ix + 1).min(nx - 1);
80    let iy1 = (iy + 1).min(ny - 1);
81    let iz1 = (iz + 1).min(nz - 1);
82
83    // 8-corner trilinear interpolation.
84    let c000 = volume.sample(ix, iy, iz);
85    let c100 = volume.sample(ix1, iy, iz);
86    let c010 = volume.sample(ix, iy1, iz);
87    let c110 = volume.sample(ix1, iy1, iz);
88    let c001 = volume.sample(ix, iy, iz1);
89    let c101 = volume.sample(ix1, iy, iz1);
90    let c011 = volume.sample(ix, iy1, iz1);
91    let c111 = volume.sample(ix1, iy1, iz1);
92
93    let c00 = c000 * (1.0 - fx) + c100 * fx;
94    let c10 = c010 * (1.0 - fx) + c110 * fx;
95    let c01 = c001 * (1.0 - fx) + c101 * fx;
96    let c11 = c011 * (1.0 - fx) + c111 * fx;
97
98    let c0 = c00 * (1.0 - fy) + c10 * fy;
99    let c1 = c01 * (1.0 - fy) + c11 * fy;
100
101    c0 * (1.0 - fz) + c1 * fz
102}
103
104/// Compute the gradient at a world-space position via central differences.
105fn gradient_at(volume: &VolumeData, pos: [f32; 3]) -> [f32; 3] {
106    let hx = volume.spacing[0] * 0.5;
107    let hy = volume.spacing[1] * 0.5;
108    let hz = volume.spacing[2] * 0.5;
109
110    let gx = trilinear_sample(volume, [pos[0] + hx, pos[1], pos[2]])
111        - trilinear_sample(volume, [pos[0] - hx, pos[1], pos[2]]);
112    let gy = trilinear_sample(volume, [pos[0], pos[1] + hy, pos[2]])
113        - trilinear_sample(volume, [pos[0], pos[1] - hy, pos[2]]);
114    let gz = trilinear_sample(volume, [pos[0], pos[1], pos[2] + hz])
115        - trilinear_sample(volume, [pos[0], pos[1], pos[2] - hz]);
116
117    let len = (gx * gx + gy * gy + gz * gz).sqrt();
118    if len > 1e-10 {
119        [gx / len, gy / len, gz / len]
120    } else {
121        [0.0, 1.0, 0.0] // fallback normal
122    }
123}
124
125/// Extract an isosurface from a volume at the given `isovalue` using marching cubes.
126///
127/// Returns a [`MeshData`] with positions, normals (from volume gradient), and triangle
128/// indices. The mesh can be uploaded to the viewport via the standard mesh pipeline.
129///
130/// Vertices along shared edges are deduplicated via an edge-vertex cache, producing
131/// a clean mesh suitable for smooth rendering.
132pub fn extract_isosurface(volume: &VolumeData, isovalue: f32) -> MeshData {
133    let [nx, ny, nz] = volume.dims;
134    if nx < 2 || ny < 2 || nz < 2 {
135        return MeshData {
136            positions: Vec::new(),
137            normals: Vec::new(),
138            indices: Vec::new(),
139            uvs: None,
140            tangents: None,
141            attributes: HashMap::new(),
142        };
143    }
144
145    let mut positions: Vec<[f32; 3]> = Vec::new();
146    let mut normals: Vec<[f32; 3]> = Vec::new();
147    let mut indices: Vec<u32> = Vec::new();
148
149    // Edge-vertex cache: key = (cell_x, cell_y, cell_z, edge_id) -> vertex index.
150    // We encode shared edges using the canonical lower-corner + axis direction.
151    let mut edge_cache: HashMap<(u32, u32, u32, u8), u32> = HashMap::new();
152
153    for iz in 0..(nz - 1) {
154        for iy in 0..(ny - 1) {
155            for ix in 0..(nx - 1) {
156                // 8 corner values in standard marching cubes order (Bourke numbering).
157                let corners = [
158                    volume.sample(ix, iy, iz),             // 0
159                    volume.sample(ix + 1, iy, iz),         // 1
160                    volume.sample(ix + 1, iy + 1, iz),     // 2
161                    volume.sample(ix, iy + 1, iz),         // 3
162                    volume.sample(ix, iy, iz + 1),         // 4
163                    volume.sample(ix + 1, iy, iz + 1),     // 5
164                    volume.sample(ix + 1, iy + 1, iz + 1), // 6
165                    volume.sample(ix, iy + 1, iz + 1),     // 7
166                ];
167
168                // Compute 8-bit cube index.
169                let mut cube_index = 0u8;
170                for (i, &val) in corners.iter().enumerate() {
171                    if val < isovalue {
172                        cube_index |= 1 << i;
173                    }
174                }
175
176                let edge_bits = EDGE_TABLE[cube_index as usize];
177                if edge_bits == 0 {
178                    continue;
179                }
180
181                // Corner world positions.
182                let corner_pos = corner_positions(volume, ix, iy, iz);
183
184                // For each triangle edge in TRI_TABLE, get or create vertex.
185                let tri_row = &TRI_TABLE[cube_index as usize];
186                let mut i = 0;
187                while i < 16 && tri_row[i] >= 0 {
188                    let mut tri_verts = [0u32; 3];
189                    for v in 0..3 {
190                        let edge_id = tri_row[i + v] as u8;
191                        let (a, b) = EDGE_VERTICES[edge_id as usize];
192
193                        // Canonical edge key: lower-corner cell + axis.
194                        let cache_key = canonical_edge_key(ix, iy, iz, edge_id);
195
196                        tri_verts[v] = *edge_cache.entry(cache_key).or_insert_with(|| {
197                            let va = corners[a as usize];
198                            let vb = corners[b as usize];
199                            let t = if (va - vb).abs() > 1e-10 {
200                                (isovalue - va) / (vb - va)
201                            } else {
202                                0.5
203                            };
204                            let t = t.clamp(0.0, 1.0);
205
206                            let pa = corner_pos[a as usize];
207                            let pb = corner_pos[b as usize];
208                            let pos = [
209                                pa[0] + t * (pb[0] - pa[0]),
210                                pa[1] + t * (pb[1] - pa[1]),
211                                pa[2] + t * (pb[2] - pa[2]),
212                            ];
213
214                            let normal = gradient_at(volume, pos);
215
216                            let idx = positions.len() as u32;
217                            positions.push(pos);
218                            normals.push(normal);
219                            idx
220                        });
221                    }
222                    // Emit triangle.
223                    indices.push(tri_verts[0]);
224                    indices.push(tri_verts[1]);
225                    indices.push(tri_verts[2]);
226                    i += 3;
227                }
228            }
229        }
230    }
231
232    MeshData {
233        positions,
234        normals,
235        indices,
236        uvs: None,
237        tangents: None,
238        attributes: HashMap::new(),
239    }
240}
241
242/// Compute world-space positions for the 8 corners of a cell.
243fn corner_positions(volume: &VolumeData, ix: u32, iy: u32, iz: u32) -> [[f32; 3]; 8] {
244    let ox = volume.origin[0] + ix as f32 * volume.spacing[0];
245    let oy = volume.origin[1] + iy as f32 * volume.spacing[1];
246    let oz = volume.origin[2] + iz as f32 * volume.spacing[2];
247    let dx = volume.spacing[0];
248    let dy = volume.spacing[1];
249    let dz = volume.spacing[2];
250
251    [
252        [ox, oy, oz],                // 0
253        [ox + dx, oy, oz],           // 1
254        [ox + dx, oy + dy, oz],      // 2
255        [ox, oy + dy, oz],           // 3
256        [ox, oy, oz + dz],           // 4
257        [ox + dx, oy, oz + dz],      // 5
258        [ox + dx, oy + dy, oz + dz], // 6
259        [ox, oy + dy, oz + dz],      // 7
260    ]
261}
262
263/// Canonical edge key for the vertex deduplication cache.
264///
265/// Shared edges between adjacent cells must produce the same key. We encode each
266/// edge as the lower-index corner cell coordinate + the axis of the edge.
267fn canonical_edge_key(cx: u32, cy: u32, cz: u32, edge_id: u8) -> (u32, u32, u32, u8) {
268    // Each edge is shared by up to 4 cells. We pick the canonical owner as the cell
269    // that has the smallest coordinates among all cells sharing this edge, encoding
270    // the edge as (cell_x, cell_y, cell_z, local_edge_axis).
271    //
272    // Edge axis encoding:
273    //   0-3: edges along X
274    //   4-7: edges along Y
275    //   8-11: edges along Z
276    match edge_id {
277        // X-axis edges
278        0 => (cx, cy, cz, 0),         // edge 0: corner 0-1, bottom-front
279        2 => (cx, cy + 1, cz, 0),     // edge 2: corner 3-2, top-front
280        4 => (cx, cy, cz + 1, 0),     // edge 4: corner 4-5, bottom-back
281        6 => (cx, cy + 1, cz + 1, 0), // edge 6: corner 7-6, top-back
282        // Y-axis edges
283        3 => (cx, cy, cz, 1),         // edge 3: corner 0-3, left-front
284        1 => (cx + 1, cy, cz, 1),     // edge 1: corner 1-2, right-front
285        7 => (cx, cy, cz + 1, 1),     // edge 7: corner 4-7, left-back
286        5 => (cx + 1, cy, cz + 1, 1), // edge 5: corner 5-6, right-back
287        // Z-axis edges
288        8 => (cx, cy, cz, 2),          // edge 8: corner 0-4, bottom-left
289        9 => (cx + 1, cy, cz, 2),      // edge 9: corner 1-5, bottom-right
290        10 => (cx + 1, cy + 1, cz, 2), // edge 10: corner 2-6, top-right
291        11 => (cx, cy + 1, cz, 2),     // edge 11: corner 3-7, top-left
292        _ => (cx, cy, cz, edge_id),    // fallback (should not happen)
293    }
294}
295
296/// Edge vertex pairs: for edge i, EDGE_VERTICES[i] = (corner_a, corner_b).
297const EDGE_VERTICES: [(u8, u8); 12] = [
298    (0, 1), // edge 0
299    (1, 2), // edge 1
300    (2, 3), // edge 2
301    (3, 0), // edge 3
302    (4, 5), // edge 4
303    (5, 6), // edge 5
304    (6, 7), // edge 6
305    (7, 4), // edge 7
306    (0, 4), // edge 8
307    (1, 5), // edge 9
308    (2, 6), // edge 10
309    (3, 7), // edge 11
310];
311
312// ---------------------------------------------------------------------------
313// Standard marching cubes lookup tables (Paul Bourke)
314// ---------------------------------------------------------------------------
315
316/// 12-bit edge intersection bitmask per cube configuration.
317#[rustfmt::skip]
318const EDGE_TABLE: [u16; 256] = [
319    0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
320    0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
321    0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
322    0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
323    0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c,
324    0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
325    0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac,
326    0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
327    0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c,
328    0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
329    0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc,
330    0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
331    0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c,
332    0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
333    0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc,
334    0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
335    0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
336    0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
337    0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
338    0x15c, 0x055, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
339    0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
340    0x2fc, 0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
341    0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
342    0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460,
343    0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
344    0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
345    0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
346    0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230,
347    0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
348    0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190,
349    0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
350    0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000,
351];
352
353/// Triangle table: up to 5 triangles (15 edge indices) per configuration. -1 = sentinel.
354#[rustfmt::skip]
355const TRI_TABLE: [[i8; 16]; 256] = [
356    [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
357    [ 0, 8, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
358    [ 0, 1, 9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
359    [ 1, 8, 3, 9, 8, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
360    [ 1, 2,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
361    [ 0, 8, 3, 1, 2,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
362    [ 9, 2,10, 0, 2, 9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
363    [ 2, 8, 3, 2,10, 8,10, 9, 8,-1,-1,-1,-1,-1,-1,-1],
364    [ 3,11, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
365    [ 0,11, 2, 8,11, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
366    [ 1, 9, 0, 2, 3,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
367    [ 1,11, 2, 1, 9,11, 9, 8,11,-1,-1,-1,-1,-1,-1,-1],
368    [ 3,10, 1,11,10, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
369    [ 0,10, 1, 0, 8,10, 8,11,10,-1,-1,-1,-1,-1,-1,-1],
370    [ 3, 9, 0, 3,11, 9,11,10, 9,-1,-1,-1,-1,-1,-1,-1],
371    [ 9, 8,10,10, 8,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
372    [ 4, 7, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
373    [ 4, 3, 0, 7, 3, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
374    [ 0, 1, 9, 8, 4, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
375    [ 4, 1, 9, 4, 7, 1, 7, 3, 1,-1,-1,-1,-1,-1,-1,-1],
376    [ 1, 2,10, 8, 4, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
377    [ 3, 4, 7, 3, 0, 4, 1, 2,10,-1,-1,-1,-1,-1,-1,-1],
378    [ 9, 2,10, 9, 0, 2, 8, 4, 7,-1,-1,-1,-1,-1,-1,-1],
379    [ 2,10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4,-1,-1,-1,-1],
380    [ 8, 4, 7, 3,11, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
381    [11, 4, 7,11, 2, 4, 2, 0, 4,-1,-1,-1,-1,-1,-1,-1],
382    [ 9, 0, 1, 8, 4, 7, 2, 3,11,-1,-1,-1,-1,-1,-1,-1],
383    [ 4, 7,11, 9, 4,11, 9,11, 2, 9, 2, 1,-1,-1,-1,-1],
384    [ 3,10, 1, 3,11,10, 7, 8, 4,-1,-1,-1,-1,-1,-1,-1],
385    [ 1,11,10, 1, 4,11, 1, 0, 4, 7,11, 4,-1,-1,-1,-1],
386    [ 4, 7, 8, 9, 0,11, 9,11,10,11, 0, 3,-1,-1,-1,-1],
387    [ 4, 7,11, 4,11, 9, 9,11,10,-1,-1,-1,-1,-1,-1,-1],
388    [ 9, 5, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
389    [ 9, 5, 4, 0, 8, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
390    [ 0, 5, 4, 1, 5, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
391    [ 8, 5, 4, 8, 3, 5, 3, 1, 5,-1,-1,-1,-1,-1,-1,-1],
392    [ 1, 2,10, 9, 5, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
393    [ 3, 0, 8, 1, 2,10, 4, 9, 5,-1,-1,-1,-1,-1,-1,-1],
394    [ 5, 2,10, 5, 4, 2, 4, 0, 2,-1,-1,-1,-1,-1,-1,-1],
395    [ 2,10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8,-1,-1,-1,-1],
396    [ 9, 5, 4, 2, 3,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
397    [ 0,11, 2, 0, 8,11, 4, 9, 5,-1,-1,-1,-1,-1,-1,-1],
398    [ 0, 5, 4, 0, 1, 5, 2, 3,11,-1,-1,-1,-1,-1,-1,-1],
399    [ 2, 1, 5, 2, 5, 8, 2, 8,11, 4, 8, 5,-1,-1,-1,-1],
400    [10, 3,11,10, 1, 3, 9, 5, 4,-1,-1,-1,-1,-1,-1,-1],
401    [ 4, 9, 5, 0, 8, 1, 8,10, 1, 8,11,10,-1,-1,-1,-1],
402    [ 5, 4, 0, 5, 0,11, 5,11,10,11, 0, 3,-1,-1,-1,-1],
403    [ 5, 4, 8, 5, 8,10,10, 8,11,-1,-1,-1,-1,-1,-1,-1],
404    [ 9, 7, 8, 5, 7, 9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
405    [ 9, 3, 0, 9, 5, 3, 5, 7, 3,-1,-1,-1,-1,-1,-1,-1],
406    [ 0, 7, 8, 0, 1, 7, 1, 5, 7,-1,-1,-1,-1,-1,-1,-1],
407    [ 1, 5, 3, 3, 5, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
408    [ 9, 7, 8, 9, 5, 7,10, 1, 2,-1,-1,-1,-1,-1,-1,-1],
409    [10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3,-1,-1,-1,-1],
410    [ 8, 0, 2, 8, 2, 5, 8, 5, 7,10, 5, 2,-1,-1,-1,-1],
411    [ 2,10, 5, 2, 5, 3, 3, 5, 7,-1,-1,-1,-1,-1,-1,-1],
412    [ 7, 9, 5, 7, 8, 9, 3,11, 2,-1,-1,-1,-1,-1,-1,-1],
413    [ 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7,11,-1,-1,-1,-1],
414    [ 2, 3,11, 0, 1, 8, 1, 7, 8, 1, 5, 7,-1,-1,-1,-1],
415    [11, 2, 1,11, 1, 7, 7, 1, 5,-1,-1,-1,-1,-1,-1,-1],
416    [ 9, 5, 8, 8, 5, 7,10, 1, 3,10, 3,11,-1,-1,-1,-1],
417    [ 5, 7, 0, 5, 0, 9, 7,11, 0, 1, 0,10,11,10, 0,-1],
418    [11,10, 0,11, 0, 3,10, 5, 0, 8, 0, 7, 5, 7, 0,-1],
419    [11,10, 5, 7,11, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
420    [10, 6, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
421    [ 0, 8, 3, 5,10, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
422    [ 9, 0, 1, 5,10, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
423    [ 1, 8, 3, 1, 9, 8, 5,10, 6,-1,-1,-1,-1,-1,-1,-1],
424    [ 1, 6, 5, 2, 6, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
425    [ 1, 6, 5, 1, 2, 6, 3, 0, 8,-1,-1,-1,-1,-1,-1,-1],
426    [ 9, 6, 5, 9, 0, 6, 0, 2, 6,-1,-1,-1,-1,-1,-1,-1],
427    [ 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8,-1,-1,-1,-1],
428    [ 2, 3,11,10, 6, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
429    [11, 0, 8,11, 2, 0,10, 6, 5,-1,-1,-1,-1,-1,-1,-1],
430    [ 0, 1, 9, 2, 3,11, 5,10, 6,-1,-1,-1,-1,-1,-1,-1],
431    [ 5,10, 6, 1, 9, 2, 9,11, 2, 9, 8,11,-1,-1,-1,-1],
432    [ 6, 3,11, 6, 5, 3, 5, 1, 3,-1,-1,-1,-1,-1,-1,-1],
433    [ 0, 8,11, 0,11, 5, 0, 5, 1, 5,11, 6,-1,-1,-1,-1],
434    [ 3,11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9,-1,-1,-1,-1],
435    [ 6, 5, 9, 6, 9,11,11, 9, 8,-1,-1,-1,-1,-1,-1,-1],
436    [ 5,10, 6, 4, 7, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
437    [ 4, 3, 0, 4, 7, 3, 6, 5,10,-1,-1,-1,-1,-1,-1,-1],
438    [ 1, 9, 0, 5,10, 6, 8, 4, 7,-1,-1,-1,-1,-1,-1,-1],
439    [10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4,-1,-1,-1,-1],
440    [ 6, 1, 2, 6, 5, 1, 4, 7, 8,-1,-1,-1,-1,-1,-1,-1],
441    [ 1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7,-1,-1,-1,-1],
442    [ 8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6,-1,-1,-1,-1],
443    [ 7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9,-1],
444    [ 3,11, 2, 7, 8, 4,10, 6, 5,-1,-1,-1,-1,-1,-1,-1],
445    [ 5,10, 6, 4, 7, 2, 4, 2, 0, 2, 7,11,-1,-1,-1,-1],
446    [ 0, 1, 9, 4, 7, 8, 2, 3,11, 5,10, 6,-1,-1,-1,-1],
447    [ 9, 2, 1, 9,11, 2, 9, 4,11, 7,11, 4, 5,10, 6,-1],
448    [ 8, 4, 7, 3,11, 5, 3, 5, 1, 5,11, 6,-1,-1,-1,-1],
449    [ 5, 1,11, 5,11, 6, 1, 0,11, 7,11, 4, 0, 4,11,-1],
450    [ 0, 5, 9, 0, 6, 5, 0, 3, 6,11, 6, 3, 8, 4, 7,-1],
451    [ 6, 5, 9, 6, 9,11, 4, 7, 9, 7,11, 9,-1,-1,-1,-1],
452    [10, 4, 9, 6, 4,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
453    [ 4,10, 6, 4, 9,10, 0, 8, 3,-1,-1,-1,-1,-1,-1,-1],
454    [10, 0, 1,10, 6, 0, 6, 4, 0,-1,-1,-1,-1,-1,-1,-1],
455    [ 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1,10,-1,-1,-1,-1],
456    [ 1, 4, 9, 1, 2, 4, 2, 6, 4,-1,-1,-1,-1,-1,-1,-1],
457    [ 3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4,-1,-1,-1,-1],
458    [ 0, 2, 4, 4, 2, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
459    [ 8, 3, 2, 8, 2, 4, 4, 2, 6,-1,-1,-1,-1,-1,-1,-1],
460    [10, 4, 9,10, 6, 4,11, 2, 3,-1,-1,-1,-1,-1,-1,-1],
461    [ 0, 8, 2, 2, 8,11, 4, 9,10, 4,10, 6,-1,-1,-1,-1],
462    [ 3,11, 2, 0, 1, 6, 0, 6, 4, 6, 1,10,-1,-1,-1,-1],
463    [ 6, 4, 1, 6, 1,10, 4, 8, 1, 2, 1,11, 8,11, 1,-1],
464    [ 9, 6, 4, 9, 3, 6, 9, 1, 3,11, 6, 3,-1,-1,-1,-1],
465    [ 8,11, 1, 8, 1, 0,11, 6, 1, 9, 1, 4, 6, 4, 1,-1],
466    [ 3,11, 6, 3, 6, 0, 0, 6, 4,-1,-1,-1,-1,-1,-1,-1],
467    [ 6, 4, 8,11, 6, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
468    [ 7,10, 6, 7, 8,10, 8, 9,10,-1,-1,-1,-1,-1,-1,-1],
469    [ 0, 7, 3, 0,10, 7, 0, 9,10, 6, 7,10,-1,-1,-1,-1],
470    [10, 6, 7, 1,10, 7, 1, 7, 8, 1, 8, 0,-1,-1,-1,-1],
471    [10, 6, 7,10, 7, 1, 1, 7, 3,-1,-1,-1,-1,-1,-1,-1],
472    [ 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7,-1,-1,-1,-1],
473    [ 2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9,-1],
474    [ 7, 8, 0, 7, 0, 6, 6, 0, 2,-1,-1,-1,-1,-1,-1,-1],
475    [ 7, 3, 2, 6, 7, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
476    [ 2, 3,11,10, 6, 8,10, 8, 9, 8, 6, 7,-1,-1,-1,-1],
477    [ 2, 0, 7, 2, 7,11, 0, 9, 7, 6, 7,10, 9,10, 7,-1],
478    [ 1, 8, 0, 1, 7, 8, 1,10, 7, 6, 7,10, 2, 3,11,-1],
479    [11, 2, 1,11, 1, 7,10, 6, 1, 6, 7, 1,-1,-1,-1,-1],
480    [ 8, 9, 6, 8, 6, 7, 9, 1, 6,11, 6, 3, 1, 3, 6,-1],
481    [ 0, 9, 1,11, 6, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
482    [ 7, 8, 0, 7, 0, 6, 3,11, 0,11, 6, 0,-1,-1,-1,-1],
483    [ 7,11, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
484    [ 7, 6,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
485    [ 3, 0, 8,11, 7, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
486    [ 0, 1, 9,11, 7, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
487    [ 8, 1, 9, 8, 3, 1,11, 7, 6,-1,-1,-1,-1,-1,-1,-1],
488    [10, 1, 2, 6,11, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
489    [ 1, 2,10, 3, 0, 8, 6,11, 7,-1,-1,-1,-1,-1,-1,-1],
490    [ 2, 9, 0, 2,10, 9, 6,11, 7,-1,-1,-1,-1,-1,-1,-1],
491    [ 6,11, 7, 2,10, 3,10, 8, 3,10, 9, 8,-1,-1,-1,-1],
492    [ 7, 2, 3, 6, 2, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
493    [ 7, 0, 8, 7, 6, 0, 6, 2, 0,-1,-1,-1,-1,-1,-1,-1],
494    [ 2, 7, 6, 2, 3, 7, 0, 1, 9,-1,-1,-1,-1,-1,-1,-1],
495    [ 1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6,-1,-1,-1,-1],
496    [10, 7, 6,10, 1, 7, 1, 3, 7,-1,-1,-1,-1,-1,-1,-1],
497    [10, 7, 6, 1, 7,10, 1, 8, 7, 1, 0, 8,-1,-1,-1,-1],
498    [ 0, 3, 7, 0, 7,10, 0,10, 9, 6,10, 7,-1,-1,-1,-1],
499    [ 7, 6,10, 7,10, 8, 8,10, 9,-1,-1,-1,-1,-1,-1,-1],
500    [ 6, 8, 4,11, 8, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
501    [ 3, 6,11, 3, 0, 6, 0, 4, 6,-1,-1,-1,-1,-1,-1,-1],
502    [ 8, 6,11, 8, 4, 6, 9, 0, 1,-1,-1,-1,-1,-1,-1,-1],
503    [ 9, 4, 6, 9, 6, 3, 9, 3, 1,11, 3, 6,-1,-1,-1,-1],
504    [ 6, 8, 4, 6,11, 8, 2,10, 1,-1,-1,-1,-1,-1,-1,-1],
505    [ 1, 2,10, 3, 0,11, 0, 6,11, 0, 4, 6,-1,-1,-1,-1],
506    [ 4,11, 8, 4, 6,11, 0, 2, 9, 2,10, 9,-1,-1,-1,-1],
507    [10, 9, 3,10, 3, 2, 9, 4, 3,11, 3, 6, 4, 6, 3,-1],
508    [ 8, 2, 3, 8, 4, 2, 4, 6, 2,-1,-1,-1,-1,-1,-1,-1],
509    [ 0, 4, 2, 4, 6, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
510    [ 1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8,-1,-1,-1,-1],
511    [ 1, 9, 4, 1, 4, 2, 2, 4, 6,-1,-1,-1,-1,-1,-1,-1],
512    [ 8, 1, 3, 8, 6, 1, 8, 4, 6, 6,10, 1,-1,-1,-1,-1],
513    [10, 1, 0,10, 0, 6, 6, 0, 4,-1,-1,-1,-1,-1,-1,-1],
514    [ 4, 6, 3, 4, 3, 8, 6,10, 3, 0, 3, 9,10, 9, 3,-1],
515    [10, 9, 4, 6,10, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
516    [ 4, 9, 5, 7, 6,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
517    [ 0, 8, 3, 4, 9, 5,11, 7, 6,-1,-1,-1,-1,-1,-1,-1],
518    [ 5, 0, 1, 5, 4, 0, 7, 6,11,-1,-1,-1,-1,-1,-1,-1],
519    [11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5,-1,-1,-1,-1],
520    [ 9, 5, 4,10, 1, 2, 7, 6,11,-1,-1,-1,-1,-1,-1,-1],
521    [ 6,11, 7, 1, 2,10, 0, 8, 3, 4, 9, 5,-1,-1,-1,-1],
522    [ 7, 6,11, 5, 4,10, 4, 2,10, 4, 0, 2,-1,-1,-1,-1],
523    [ 3, 4, 8, 3, 5, 4, 3, 2, 5,10, 5, 2,11, 7, 6,-1],
524    [ 7, 2, 3, 7, 6, 2, 5, 4, 9,-1,-1,-1,-1,-1,-1,-1],
525    [ 9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7,-1,-1,-1,-1],
526    [ 3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0,-1,-1,-1,-1],
527    [ 6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8,-1],
528    [ 9, 5, 4,10, 1, 6, 1, 7, 6, 1, 3, 7,-1,-1,-1,-1],
529    [ 1, 6,10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4,-1],
530    [ 4, 0,10, 4,10, 5, 0, 3,10, 6,10, 7, 3, 7,10,-1],
531    [ 7, 6,10, 7,10, 8, 5, 4,10, 4, 8,10,-1,-1,-1,-1],
532    [ 6, 9, 5, 6,11, 9,11, 8, 9,-1,-1,-1,-1,-1,-1,-1],
533    [ 3, 6,11, 0, 6, 3, 0, 5, 6, 0, 9, 5,-1,-1,-1,-1],
534    [ 0,11, 8, 0, 5,11, 0, 1, 5, 5, 6,11,-1,-1,-1,-1],
535    [ 6,11, 3, 6, 3, 5, 5, 3, 1,-1,-1,-1,-1,-1,-1,-1],
536    [ 1, 2,10, 9, 5,11, 9,11, 8,11, 5, 6,-1,-1,-1,-1],
537    [ 0,11, 3, 0, 6,11, 0, 9, 6, 5, 6, 9, 1, 2,10,-1],
538    [11, 8, 5,11, 5, 6, 8, 0, 5,10, 5, 2, 0, 2, 5,-1],
539    [ 6,11, 3, 6, 3, 5, 2,10, 3,10, 5, 3,-1,-1,-1,-1],
540    [ 5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2,-1,-1,-1,-1],
541    [ 9, 5, 6, 9, 6, 0, 0, 6, 2,-1,-1,-1,-1,-1,-1,-1],
542    [ 1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8,-1],
543    [ 1, 5, 6, 2, 1, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
544    [ 1, 3, 6, 1, 6,10, 3, 8, 6, 5, 6, 9, 8, 9, 6,-1],
545    [10, 1, 0,10, 0, 6, 9, 5, 0, 5, 6, 0,-1,-1,-1,-1],
546    [ 0, 3, 8, 5, 6,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
547    [10, 5, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
548    [11, 5,10, 7, 5,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
549    [11, 5,10,11, 7, 5, 8, 3, 0,-1,-1,-1,-1,-1,-1,-1],
550    [ 5,11, 7, 5,10,11, 1, 9, 0,-1,-1,-1,-1,-1,-1,-1],
551    [10, 7, 5,10,11, 7, 9, 8, 1, 8, 3, 1,-1,-1,-1,-1],
552    [11, 1, 2,11, 7, 1, 7, 5, 1,-1,-1,-1,-1,-1,-1,-1],
553    [ 0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2,11,-1,-1,-1,-1],
554    [ 9, 7, 5, 9, 2, 7, 9, 0, 2, 2,11, 7,-1,-1,-1,-1],
555    [ 7, 5, 2, 7, 2,11, 5, 9, 2, 3, 2, 8, 9, 8, 2,-1],
556    [ 2, 5,10, 2, 3, 5, 3, 7, 5,-1,-1,-1,-1,-1,-1,-1],
557    [ 8, 2, 0, 8, 5, 2, 8, 7, 5,10, 2, 5,-1,-1,-1,-1],
558    [ 9, 0, 1, 5,10, 3, 5, 3, 7, 3,10, 2,-1,-1,-1,-1],
559    [ 9, 8, 2, 9, 2, 1, 8, 7, 2,10, 2, 5, 7, 5, 2,-1],
560    [ 1, 3, 5, 3, 7, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
561    [ 0, 8, 7, 0, 7, 1, 1, 7, 5,-1,-1,-1,-1,-1,-1,-1],
562    [ 9, 0, 3, 9, 3, 5, 5, 3, 7,-1,-1,-1,-1,-1,-1,-1],
563    [ 9, 8, 7, 5, 9, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
564    [ 5, 8, 4, 5,10, 8,10,11, 8,-1,-1,-1,-1,-1,-1,-1],
565    [ 5, 0, 4, 5,11, 0, 5,10,11,11, 3, 0,-1,-1,-1,-1],
566    [ 0, 1, 9, 8, 4,10, 8,10,11,10, 4, 5,-1,-1,-1,-1],
567    [10,11, 4,10, 4, 5,11, 3, 4, 9, 4, 1, 3, 1, 4,-1],
568    [ 2, 5, 1, 2, 8, 5, 2,11, 8, 4, 5, 8,-1,-1,-1,-1],
569    [ 0, 4,11, 0,11, 3, 4, 5,11, 2,11, 1, 5, 1,11,-1],
570    [ 0, 2, 5, 0, 5, 9, 2,11, 5, 4, 5, 8,11, 8, 5,-1],
571    [ 9, 4, 5, 2,11, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
572    [ 2, 5,10, 3, 5, 2, 3, 4, 5, 3, 8, 4,-1,-1,-1,-1],
573    [ 5,10, 2, 5, 2, 4, 4, 2, 0,-1,-1,-1,-1,-1,-1,-1],
574    [ 3,10, 2, 3, 5,10, 3, 8, 5, 4, 5, 8, 0, 1, 9,-1],
575    [ 5,10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2,-1,-1,-1,-1],
576    [ 8, 4, 5, 8, 5, 3, 3, 5, 1,-1,-1,-1,-1,-1,-1,-1],
577    [ 0, 4, 5, 1, 0, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
578    [ 8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5,-1,-1,-1,-1],
579    [ 9, 4, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
580    [ 4,11, 7, 4, 9,11, 9,10,11,-1,-1,-1,-1,-1,-1,-1],
581    [ 0, 8, 3, 4, 9, 7, 9,11, 7, 9,10,11,-1,-1,-1,-1],
582    [ 1,10,11, 1,11, 4, 1, 4, 0, 7, 4,11,-1,-1,-1,-1],
583    [ 3, 1, 4, 3, 4, 8, 1,10, 4, 7, 4,11,10,11, 4,-1],
584    [ 4,11, 7, 9,11, 4, 9, 2,11, 9, 1, 2,-1,-1,-1,-1],
585    [ 9, 7, 4, 9,11, 7, 9, 1,11, 2,11, 1, 0, 8, 3,-1],
586    [11, 7, 4,11, 4, 2, 2, 4, 0,-1,-1,-1,-1,-1,-1,-1],
587    [11, 7, 4,11, 4, 2, 8, 3, 4, 3, 2, 4,-1,-1,-1,-1],
588    [ 2, 9,10, 2, 7, 9, 2, 3, 7, 7, 4, 9,-1,-1,-1,-1],
589    [ 9,10, 7, 9, 7, 4,10, 2, 7, 8, 7, 0, 2, 0, 7,-1],
590    [ 3, 7,10, 3,10, 2, 7, 4,10, 1,10, 0, 4, 0,10,-1],
591    [ 1,10, 2, 8, 7, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
592    [ 4, 9, 1, 4, 1, 7, 7, 1, 3,-1,-1,-1,-1,-1,-1,-1],
593    [ 4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1,-1,-1,-1,-1],
594    [ 4, 0, 3, 7, 4, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
595    [ 4, 8, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
596    [ 9,10, 8,10,11, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
597    [ 3, 0, 9, 3, 9,11,11, 9,10,-1,-1,-1,-1,-1,-1,-1],
598    [ 0, 1,10, 0,10, 8, 8,10,11,-1,-1,-1,-1,-1,-1,-1],
599    [ 3, 1,10,11, 3,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
600    [ 1, 2,11, 1,11, 9, 9,11, 8,-1,-1,-1,-1,-1,-1,-1],
601    [ 3, 0, 9, 3, 9,11, 1, 2, 9, 2,11, 9,-1,-1,-1,-1],
602    [ 0, 2,11, 8, 0,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
603    [ 3, 2,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
604    [ 2, 3, 8, 2, 8,10,10, 8, 9,-1,-1,-1,-1,-1,-1,-1],
605    [ 9,10, 2, 0, 9, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
606    [ 2, 3, 8, 2, 8,10, 0, 1, 8, 1,10, 8,-1,-1,-1,-1],
607    [ 1,10, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
608    [ 1, 3, 8, 9, 1, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
609    [ 0, 9, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
610    [ 0, 3, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
611    [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
612];
613
614#[cfg(test)]
615mod tests {
616    use super::*;
617
618    #[test]
619    fn test_sphere_isosurface() {
620        let n = 32u32;
621        let mut data = vec![0.0f32; (n * n * n) as usize];
622        let center = n as f32 / 2.0;
623
624        for iz in 0..n {
625            for iy in 0..n {
626                for ix in 0..n {
627                    let dx = ix as f32 - center;
628                    let dy = iy as f32 - center;
629                    let dz = iz as f32 - center;
630                    let dist = (dx * dx + dy * dy + dz * dz).sqrt();
631                    data[(ix + iy * n + iz * n * n) as usize] = dist;
632                }
633            }
634        }
635
636        let volume = VolumeData {
637            data,
638            dims: [n, n, n],
639            origin: [0.0, 0.0, 0.0],
640            spacing: [1.0, 1.0, 1.0],
641        };
642
643        let mesh = extract_isosurface(&volume, 8.0);
644
645        // Mesh was generated.
646        assert!(
647            !mesh.positions.is_empty(),
648            "Sphere isosurface should produce vertices"
649        );
650
651        // Valid triangle list.
652        assert_eq!(mesh.indices.len() % 3, 0, "Indices must be a multiple of 3");
653
654        // Reasonable triangle count for a sphere.
655        assert!(
656            mesh.indices.len() > 100,
657            "Expected > 100 indices for sphere, got {}",
658            mesh.indices.len()
659        );
660
661        // All positions within expected bounding box (center 16 +/- radius 8 + margin).
662        for pos in &mesh.positions {
663            for c in pos {
664                assert!(
665                    *c >= 4.0 && *c <= 28.0,
666                    "Position component {} out of expected range [4, 28]",
667                    c
668                );
669            }
670        }
671
672        // All normals approximately unit length.
673        for n in &mesh.normals {
674            let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
675            assert!(
676                len > 0.95 && len < 1.05,
677                "Normal length {} not approximately 1.0",
678                len
679            );
680        }
681
682        // positions and normals must be same length.
683        assert_eq!(mesh.positions.len(), mesh.normals.len());
684    }
685
686    #[test]
687    fn test_empty_volume() {
688        // All values above isovalue.
689        let data = vec![10.0f32; 8]; // 2x2x2, all 10.0
690        let volume = VolumeData {
691            data,
692            dims: [2, 2, 2],
693            origin: [0.0, 0.0, 0.0],
694            spacing: [1.0, 1.0, 1.0],
695        };
696
697        let mesh = extract_isosurface(&volume, 5.0);
698        assert_eq!(
699            mesh.positions.len(),
700            0,
701            "All-above should produce empty mesh"
702        );
703        assert_eq!(mesh.indices.len(), 0);
704    }
705
706    #[test]
707    fn test_volume_data_sample() {
708        let data = vec![
709            1.0, 2.0, // z=0 row: (0,0,0)=1, (1,0,0)=2
710            3.0, 4.0, // z=0 row: (0,1,0)=3, (1,1,0)=4
711            5.0, 6.0, // z=1 row: (0,0,1)=5, (1,0,1)=6
712            7.0, 8.0, // z=1 row: (0,1,1)=7, (1,1,1)=8
713        ];
714        let vol = VolumeData {
715            data,
716            dims: [2, 2, 2],
717            origin: [0.0, 0.0, 0.0],
718            spacing: [1.0, 1.0, 1.0],
719        };
720
721        assert_eq!(vol.sample(0, 0, 0), 1.0);
722        assert_eq!(vol.sample(1, 0, 0), 2.0);
723        assert_eq!(vol.sample(0, 1, 0), 3.0);
724        assert_eq!(vol.sample(1, 1, 0), 4.0);
725        assert_eq!(vol.sample(0, 0, 1), 5.0);
726        assert_eq!(vol.sample(1, 0, 1), 6.0);
727        assert_eq!(vol.sample(0, 1, 1), 7.0);
728        assert_eq!(vol.sample(1, 1, 1), 8.0);
729    }
730
731    #[test]
732    fn test_trilinear_interpolation() {
733        let data = vec![
734            0.0, 1.0, // (0,0,0)=0, (1,0,0)=1
735            0.0, 1.0, // (0,1,0)=0, (1,1,0)=1
736            0.0, 1.0, // (0,0,1)=0, (1,0,1)=1
737            0.0, 1.0, // (0,1,1)=0, (1,1,1)=1
738        ];
739        let vol = VolumeData {
740            data,
741            dims: [2, 2, 2],
742            origin: [0.0, 0.0, 0.0],
743            spacing: [1.0, 1.0, 1.0],
744        };
745
746        // At grid points.
747        let v00 = trilinear_sample(&vol, [0.0, 0.0, 0.0]);
748        assert!(
749            (v00 - 0.0).abs() < 0.01,
750            "Expected 0.0 at origin, got {}",
751            v00
752        );
753
754        let v10 = trilinear_sample(&vol, [0.999, 0.0, 0.0]);
755        assert!(
756            (v10 - 1.0).abs() < 0.02,
757            "Expected ~1.0 at (1,0,0), got {}",
758            v10
759        );
760
761        // At midpoint along X: linear gradient means value = 0.5.
762        let mid = trilinear_sample(&vol, [0.5, 0.0, 0.0]);
763        assert!(
764            (mid - 0.5).abs() < 0.01,
765            "Expected 0.5 at midpoint, got {}",
766            mid
767        );
768
769        // At midpoint in all axes.
770        let center = trilinear_sample(&vol, [0.5, 0.5, 0.5]);
771        assert!(
772            (center - 0.5).abs() < 0.01,
773            "Expected 0.5 at center, got {}",
774            center
775        );
776    }
777}
viewport_lib/geometry/marching_cubes.rs

viewport_lib/geometry/
marching_cubes.rs
