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 (swap v1/v2 to match renderer's CCW front-face convention).
223                    indices.push(tri_verts[0]);
224                    indices.push(tri_verts[2]);
225                    indices.push(tri_verts[1]);
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]
355pub(crate) const 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_tri_table_edge_consistency() {
620        // Every edge index referenced in TRI_TABLE[i] must appear in EDGE_TABLE[i].
621        let mut failures = Vec::new();
622        for cube_index in 0u16..256 {
623            let edge_bits = EDGE_TABLE[cube_index as usize];
624            let tri_row = &TRI_TABLE[cube_index as usize];
625            let mut j = 0;
626            while j < 16 && tri_row[j] >= 0 {
627                let edge_id = tri_row[j] as u8;
628                if edge_bits & (1 << edge_id) == 0 {
629                    failures.push(format!(
630                        "TRI_TABLE[{}]: edge {} not in EDGE_TABLE ({:#014b})",
631                        cube_index, edge_id, edge_bits
632                    ));
633                    break; // one failure per cube_index is enough
634                }
635                j += 1;
636            }
637        }
638        if !failures.is_empty() {
639            panic!(
640                "{} TRI_TABLE entries inconsistent with EDGE_TABLE:\n{}",
641                failures.len(),
642                failures.join("\n")
643            );
644        }
645    }
646
647    #[test]
648    fn test_sphere_isosurface() {
649        let n = 32u32;
650        let mut data = vec![0.0f32; (n * n * n) as usize];
651        let center = n as f32 / 2.0;
652
653        for iz in 0..n {
654            for iy in 0..n {
655                for ix in 0..n {
656                    let dx = ix as f32 - center;
657                    let dy = iy as f32 - center;
658                    let dz = iz as f32 - center;
659                    let dist = (dx * dx + dy * dy + dz * dz).sqrt();
660                    data[(ix + iy * n + iz * n * n) as usize] = dist;
661                }
662            }
663        }
664
665        let volume = VolumeData {
666            data,
667            dims: [n, n, n],
668            origin: [0.0, 0.0, 0.0],
669            spacing: [1.0, 1.0, 1.0],
670        };
671
672        let mesh = extract_isosurface(&volume, 8.0);
673
674        // Mesh was generated.
675        assert!(
676            !mesh.positions.is_empty(),
677            "Sphere isosurface should produce vertices"
678        );
679
680        // Valid triangle list.
681        assert_eq!(mesh.indices.len() % 3, 0, "Indices must be a multiple of 3");
682
683        // Reasonable triangle count for a sphere.
684        assert!(
685            mesh.indices.len() > 100,
686            "Expected > 100 indices for sphere, got {}",
687            mesh.indices.len()
688        );
689
690        // All positions within expected bounding box (center 16 +/- radius 8 + margin).
691        for pos in &mesh.positions {
692            for c in pos {
693                assert!(
694                    *c >= 4.0 && *c <= 28.0,
695                    "Position component {} out of expected range [4, 28]",
696                    c
697                );
698            }
699        }
700
701        // All normals approximately unit length.
702        for n in &mesh.normals {
703            let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
704            assert!(
705                len > 0.95 && len < 1.05,
706                "Normal length {} not approximately 1.0",
707                len
708            );
709        }
710
711        // positions and normals must be same length.
712        assert_eq!(mesh.positions.len(), mesh.normals.len());
713    }
714
715    #[test]
716    fn test_sphere_winding_order() {
717        // Extract a sphere isosurface and verify geometric normals (cross product)
718        // align with the gradient-based vertex normals (which point outward for SDF).
719        // A winding mismatch means some fraction of triangles would be back-face culled.
720        let n = 32u32;
721        let center = n as f32 / 2.0;
722        let mut data = vec![0.0f32; (n * n * n) as usize];
723        for iz in 0..n {
724            for iy in 0..n {
725                for ix in 0..n {
726                    let dx = ix as f32 - center;
727                    let dy = iy as f32 - center;
728                    let dz = iz as f32 - center;
729                    data[(ix + iy * n + iz * n * n) as usize] =
730                        (dx * dx + dy * dy + dz * dz).sqrt();
731                }
732            }
733        }
734        let volume = VolumeData {
735            data,
736            dims: [n, n, n],
737            origin: [0.0, 0.0, 0.0],
738            spacing: [1.0, 1.0, 1.0],
739        };
740        let mesh = extract_isosurface(&volume, 8.0);
741        assert!(!mesh.positions.is_empty(), "expected vertices");
742
743        let mut correct = 0usize;
744        let mut flipped = 0usize;
745        let tri_count = mesh.indices.len() / 3;
746        for t in 0..tri_count {
747            let i0 = mesh.indices[t * 3] as usize;
748            let i1 = mesh.indices[t * 3 + 1] as usize;
749            let i2 = mesh.indices[t * 3 + 2] as usize;
750            let p0 = mesh.positions[i0];
751            let p1 = mesh.positions[i1];
752            let p2 = mesh.positions[i2];
753            // Geometric normal via cross product.
754            let e1 = [p1[0]-p0[0], p1[1]-p0[1], p1[2]-p0[2]];
755            let e2 = [p2[0]-p0[0], p2[1]-p0[1], p2[2]-p0[2]];
756            let gn = [
757                e1[1]*e2[2] - e1[2]*e2[1],
758                e1[2]*e2[0] - e1[0]*e2[2],
759                e1[0]*e2[1] - e1[1]*e2[0],
760            ];
761            // Vertex normal (gradient-based, points outward for this SDF).
762            let vn = mesh.normals[i0];
763            let dot = gn[0]*vn[0] + gn[1]*vn[1] + gn[2]*vn[2];
764            if dot >= 0.0 { correct += 1; } else { flipped += 1; }
765        }
766
767        let total = correct + flipped;
768        let flipped_pct = flipped as f32 / total as f32 * 100.0;
769        assert!(
770            flipped_pct < 5.0,
771            "{}/{} triangles ({:.1}%) have geometric normal opposing the gradient — winding is wrong",
772            flipped, total, flipped_pct
773        );
774    }
775
776    #[test]
777    fn test_empty_volume() {
778        // All values above isovalue.
779        let data = vec![10.0f32; 8]; // 2x2x2, all 10.0
780        let volume = VolumeData {
781            data,
782            dims: [2, 2, 2],
783            origin: [0.0, 0.0, 0.0],
784            spacing: [1.0, 1.0, 1.0],
785        };
786
787        let mesh = extract_isosurface(&volume, 5.0);
788        assert_eq!(
789            mesh.positions.len(),
790            0,
791            "All-above should produce empty mesh"
792        );
793        assert_eq!(mesh.indices.len(), 0);
794    }
795
796    #[test]
797    fn test_volume_data_sample() {
798        let data = vec![
799            1.0, 2.0, // z=0 row: (0,0,0)=1, (1,0,0)=2
800            3.0, 4.0, // z=0 row: (0,1,0)=3, (1,1,0)=4
801            5.0, 6.0, // z=1 row: (0,0,1)=5, (1,0,1)=6
802            7.0, 8.0, // z=1 row: (0,1,1)=7, (1,1,1)=8
803        ];
804        let vol = VolumeData {
805            data,
806            dims: [2, 2, 2],
807            origin: [0.0, 0.0, 0.0],
808            spacing: [1.0, 1.0, 1.0],
809        };
810
811        assert_eq!(vol.sample(0, 0, 0), 1.0);
812        assert_eq!(vol.sample(1, 0, 0), 2.0);
813        assert_eq!(vol.sample(0, 1, 0), 3.0);
814        assert_eq!(vol.sample(1, 1, 0), 4.0);
815        assert_eq!(vol.sample(0, 0, 1), 5.0);
816        assert_eq!(vol.sample(1, 0, 1), 6.0);
817        assert_eq!(vol.sample(0, 1, 1), 7.0);
818        assert_eq!(vol.sample(1, 1, 1), 8.0);
819    }
820
821    #[test]
822    fn test_trilinear_interpolation() {
823        let data = vec![
824            0.0, 1.0, // (0,0,0)=0, (1,0,0)=1
825            0.0, 1.0, // (0,1,0)=0, (1,1,0)=1
826            0.0, 1.0, // (0,0,1)=0, (1,0,1)=1
827            0.0, 1.0, // (0,1,1)=0, (1,1,1)=1
828        ];
829        let vol = VolumeData {
830            data,
831            dims: [2, 2, 2],
832            origin: [0.0, 0.0, 0.0],
833            spacing: [1.0, 1.0, 1.0],
834        };
835
836        // At grid points.
837        let v00 = trilinear_sample(&vol, [0.0, 0.0, 0.0]);
838        assert!(
839            (v00 - 0.0).abs() < 0.01,
840            "Expected 0.0 at origin, got {}",
841            v00
842        );
843
844        let v10 = trilinear_sample(&vol, [0.999, 0.0, 0.0]);
845        assert!(
846            (v10 - 1.0).abs() < 0.02,
847            "Expected ~1.0 at (1,0,0), got {}",
848            v10
849        );
850
851        // At midpoint along X: linear gradient means value = 0.5.
852        let mid = trilinear_sample(&vol, [0.5, 0.0, 0.0]);
853        assert!(
854            (mid - 0.5).abs() < 0.01,
855            "Expected 0.5 at midpoint, got {}",
856            mid
857        );
858
859        // At midpoint in all axes.
860        let center = trilinear_sample(&vol, [0.5, 0.5, 0.5]);
861        assert!(
862            (center - 0.5).abs() < 0.01,
863            "Expected 0.5 at center, got {}",
864            center
865        );
866    }
867}
viewport_lib/geometry/marching_cubes.rs

viewport_lib/geometry/
marching_cubes.rs
