sci_form/eht/
marching_cubes.rs

1//! Marching Cubes isosurface extraction from volumetric grids.
2//!
3//! Given a VolumetricGrid and an isovalue, produces a triangle mesh
4//! (vertices, normals, indices) suitable for WebGL/WebGPU rendering.
5
6use super::volume::VolumetricGrid;
7use serde::{Deserialize, Serialize};
8
9/// A triangle mesh extracted from an isosurface.
10#[derive(Debug, Clone, Serialize, Deserialize)]
11pub struct IsosurfaceMesh {
12    /// Flat array of vertex positions [x0,y0,z0, x1,y1,z1, ...] in Ångström.
13    pub vertices: Vec<f32>,
14    /// Flat array of vertex normals [nx0,ny0,nz0, ...].
15    pub normals: Vec<f32>,
16    /// Triangle indices (every 3 = one triangle).
17    pub indices: Vec<u32>,
18    /// Number of triangles.
19    pub num_triangles: usize,
20    /// The isovalue used.
21    pub isovalue: f32,
22}
23
24impl IsosurfaceMesh {
25    pub fn num_vertices(&self) -> usize {
26        self.vertices.len() / 3
27    }
28}
29
30/// Extract an isosurface from a volumetric grid using Marching Cubes.
31pub fn marching_cubes(grid: &VolumetricGrid, isovalue: f32) -> IsosurfaceMesh {
32    let [nx, ny, nz] = grid.dims;
33    let mut vertices = Vec::new();
34    let mut normals = Vec::new();
35    let mut indices = Vec::new();
36
37    if nx < 2 || ny < 2 || nz < 2 {
38        return IsosurfaceMesh {
39            vertices,
40            normals,
41            indices,
42            num_triangles: 0,
43            isovalue,
44        };
45    }
46
47    for ix in 0..nx - 1 {
48        for iy in 0..ny - 1 {
49            for iz in 0..nz - 1 {
50                process_cube(
51                    grid,
52                    ix,
53                    iy,
54                    iz,
55                    isovalue,
56                    &mut vertices,
57                    &mut normals,
58                    &mut indices,
59                );
60            }
61        }
62    }
63
64    let num_triangles = indices.len() / 3;
65    IsosurfaceMesh {
66        vertices,
67        normals,
68        indices,
69        num_triangles,
70        isovalue,
71    }
72}
73
74fn process_cube(
75    grid: &VolumetricGrid,
76    ix: usize,
77    iy: usize,
78    iz: usize,
79    isovalue: f32,
80    vertices: &mut Vec<f32>,
81    normals: &mut Vec<f32>,
82    indices: &mut Vec<u32>,
83) {
84    // 8 corners of the cube
85    let corners = [
86        (ix, iy, iz),
87        (ix + 1, iy, iz),
88        (ix + 1, iy + 1, iz),
89        (ix, iy + 1, iz),
90        (ix, iy, iz + 1),
91        (ix + 1, iy, iz + 1),
92        (ix + 1, iy + 1, iz + 1),
93        (ix, iy + 1, iz + 1),
94    ];
95
96    let vals: [f32; 8] = [
97        grid.values[grid.index(corners[0].0, corners[0].1, corners[0].2)],
98        grid.values[grid.index(corners[1].0, corners[1].1, corners[1].2)],
99        grid.values[grid.index(corners[2].0, corners[2].1, corners[2].2)],
100        grid.values[grid.index(corners[3].0, corners[3].1, corners[3].2)],
101        grid.values[grid.index(corners[4].0, corners[4].1, corners[4].2)],
102        grid.values[grid.index(corners[5].0, corners[5].1, corners[5].2)],
103        grid.values[grid.index(corners[6].0, corners[6].1, corners[6].2)],
104        grid.values[grid.index(corners[7].0, corners[7].1, corners[7].2)],
105    ];
106
107    // Compute cube index from corner classifications
108    let mut cube_idx = 0u8;
109    for i in 0..8 {
110        if vals[i] >= isovalue {
111            cube_idx |= 1 << i;
112        }
113    }
114
115    if cube_idx == 0 || cube_idx == 255 {
116        return; // Entirely inside or outside
117    }
118
119    let edge_flags = EDGE_TABLE[cube_idx as usize];
120    if edge_flags == 0 {
121        return;
122    }
123
124    let positions: [[f64; 3]; 8] =
125        core::array::from_fn(|i| grid.point_position(corners[i].0, corners[i].1, corners[i].2));
126
127    // Interpolate edge vertices
128    let mut edge_verts = [[0.0f64; 3]; 12];
129    let edges = [
130        (0, 1),
131        (1, 2),
132        (2, 3),
133        (3, 0),
134        (4, 5),
135        (5, 6),
136        (6, 7),
137        (7, 4),
138        (0, 4),
139        (1, 5),
140        (2, 6),
141        (3, 7),
142    ];
143
144    for (ei, &(a, b)) in edges.iter().enumerate() {
145        if edge_flags & (1 << ei) != 0 {
146            edge_verts[ei] =
147                interpolate_vertex(&positions[a], &positions[b], vals[a], vals[b], isovalue);
148        }
149    }
150
151    // Emit triangles
152    let tri = &TRI_TABLE[cube_idx as usize];
153    let mut i = 0;
154    while i < 16 && tri[i] != -1 {
155        let base_idx = (vertices.len() / 3) as u32;
156
157        for k in 0..3 {
158            let edge = tri[i + k] as usize;
159            let v = edge_verts[edge];
160            vertices.push(v[0] as f32);
161            vertices.push(v[1] as f32);
162            vertices.push(v[2] as f32);
163
164            // Approximate normal from gradient at vertex position
165            let n = estimate_normal(grid, &v, isovalue);
166            normals.push(n[0] as f32);
167            normals.push(n[1] as f32);
168            normals.push(n[2] as f32);
169        }
170
171        indices.push(base_idx);
172        indices.push(base_idx + 1);
173        indices.push(base_idx + 2);
174
175        i += 3;
176    }
177}
178
179fn interpolate_vertex(p1: &[f64; 3], p2: &[f64; 3], v1: f32, v2: f32, iso: f32) -> [f64; 3] {
180    let diff = v2 - v1;
181    let t = if diff.abs() < 1e-10 {
182        0.5
183    } else {
184        ((iso - v1) / diff) as f64
185    };
186    [
187        p1[0] + t * (p2[0] - p1[0]),
188        p1[1] + t * (p2[1] - p1[1]),
189        p1[2] + t * (p2[2] - p1[2]),
190    ]
191}
192
193fn estimate_normal(grid: &VolumetricGrid, point: &[f64; 3], _iso: f32) -> [f64; 3] {
194    let h = grid.spacing;
195    let sample = |p: &[f64; 3]| -> f64 {
196        // Nearest-neighbor sampling
197        let ix = ((p[0] - grid.origin[0]) / h).round() as isize;
198        let iy = ((p[1] - grid.origin[1]) / h).round() as isize;
199        let iz = ((p[2] - grid.origin[2]) / h).round() as isize;
200        if ix < 0
201            || iy < 0
202            || iz < 0
203            || ix >= grid.dims[0] as isize
204            || iy >= grid.dims[1] as isize
205            || iz >= grid.dims[2] as isize
206        {
207            return 0.0;
208        }
209        grid.values[grid.index(ix as usize, iy as usize, iz as usize)] as f64
210    };
211
212    let gx =
213        sample(&[point[0] + h, point[1], point[2]]) - sample(&[point[0] - h, point[1], point[2]]);
214    let gy =
215        sample(&[point[0], point[1] + h, point[2]]) - sample(&[point[0], point[1] - h, point[2]]);
216    let gz =
217        sample(&[point[0], point[1], point[2] + h]) - sample(&[point[0], point[1], point[2] - h]);
218
219    let len = (gx * gx + gy * gy + gz * gz).sqrt();
220    if len < 1e-12 {
221        [0.0, 0.0, 1.0]
222    } else {
223        [-gx / len, -gy / len, -gz / len]
224    }
225}
226
227/// Dual-phase isosurface result for orbital visualization.
228///
229/// Extracts both positive and negative lobes of an orbital (or ESP, etc.)
230/// at the given absolute isovalue, producing two separate meshes.
231#[derive(Debug, Clone, Serialize, Deserialize)]
232pub struct DualPhaseMesh {
233    /// Positive lobe (values >= +isovalue).
234    pub positive: IsosurfaceMesh,
235    /// Negative lobe (values <= -isovalue).
236    pub negative: IsosurfaceMesh,
237}
238
239/// Extract dual-phase isosurfaces (positive and negative lobes).
240///
241/// Useful for orbital visualization where positive/negative phases
242/// should be rendered in different colors.
243pub fn marching_cubes_dual(grid: &VolumetricGrid, isovalue: f32) -> DualPhaseMesh {
244    let positive = marching_cubes(grid, isovalue.abs());
245    let negative = marching_cubes(grid, -isovalue.abs());
246    DualPhaseMesh { positive, negative }
247}
248
249/// Mesh simplification: weld duplicate vertices and remove degenerate triangles.
250///
251/// Vertices closer than `weld_distance` are merged. Degenerate triangles
252/// (with zero area or repeated vertex indices) are removed.
253pub fn simplify_mesh(mesh: &IsosurfaceMesh, weld_distance: f32) -> IsosurfaceMesh {
254    let n_verts = mesh.vertices.len() / 3;
255    if n_verts == 0 {
256        return mesh.clone();
257    }
258
259    let weld_dist_sq = weld_distance * weld_distance;
260
261    // Build vertex map: for each vertex, find the canonical (first seen) vertex index
262    let mut vertex_map: Vec<usize> = vec![0; n_verts];
263    let mut new_vertices: Vec<[f32; 3]> = Vec::new();
264    let mut new_normals: Vec<[f32; 3]> = Vec::new();
265    let mut old_to_new: Vec<usize> = vec![0; n_verts];
266
267    for i in 0..n_verts {
268        let vx = mesh.vertices[3 * i];
269        let vy = mesh.vertices[3 * i + 1];
270        let vz = mesh.vertices[3 * i + 2];
271
272        // Check if this vertex matches an existing one
273        let mut found = None;
274        for (j, v) in new_vertices.iter().enumerate() {
275            let dx = vx - v[0];
276            let dy = vy - v[1];
277            let dz = vz - v[2];
278            if dx * dx + dy * dy + dz * dz < weld_dist_sq {
279                found = Some(j);
280                break;
281            }
282        }
283
284        if let Some(j) = found {
285            old_to_new[i] = j;
286        } else {
287            old_to_new[i] = new_vertices.len();
288            new_vertices.push([vx, vy, vz]);
289            new_normals.push([
290                mesh.normals[3 * i],
291                mesh.normals[3 * i + 1],
292                mesh.normals[3 * i + 2],
293            ]);
294        }
295        vertex_map[i] = old_to_new[i];
296    }
297
298    // Rebuild indices with mapped vertices, removing degenerate triangles
299    let mut new_indices = Vec::new();
300    let n_tris = mesh.indices.len() / 3;
301    for t in 0..n_tris {
302        let i0 = vertex_map[mesh.indices[3 * t] as usize] as u32;
303        let i1 = vertex_map[mesh.indices[3 * t + 1] as usize] as u32;
304        let i2 = vertex_map[mesh.indices[3 * t + 2] as usize] as u32;
305
306        // Skip degenerate triangles
307        if i0 == i1 || i1 == i2 || i0 == i2 {
308            continue;
309        }
310
311        // Check for zero-area triangle
312        let v0 = new_vertices[i0 as usize];
313        let v1 = new_vertices[i1 as usize];
314        let v2 = new_vertices[i2 as usize];
315        let e1 = [v1[0] - v0[0], v1[1] - v0[1], v1[2] - v0[2]];
316        let e2 = [v2[0] - v0[0], v2[1] - v0[1], v2[2] - v0[2]];
317        let cx = e1[1] * e2[2] - e1[2] * e2[1];
318        let cy = e1[2] * e2[0] - e1[0] * e2[2];
319        let cz = e1[0] * e2[1] - e1[1] * e2[0];
320        let area_sq = cx * cx + cy * cy + cz * cz;
321        if area_sq < 1e-12 {
322            continue;
323        }
324
325        new_indices.push(i0);
326        new_indices.push(i1);
327        new_indices.push(i2);
328    }
329
330    let flat_verts: Vec<f32> = new_vertices
331        .iter()
332        .flat_map(|v| v.iter().copied())
333        .collect();
334    let flat_norms: Vec<f32> = new_normals.iter().flat_map(|n| n.iter().copied()).collect();
335
336    IsosurfaceMesh {
337        vertices: flat_verts,
338        normals: flat_norms,
339        num_triangles: new_indices.len() / 3,
340        indices: new_indices,
341        isovalue: mesh.isovalue,
342    }
343}
344
345/// Compute angle-weighted vertex normals for smooth shading.
346///
347/// For each vertex, averages the face normals of adjacent triangles,
348/// weighted by the angle at that vertex. This produces smoother shading
349/// than per-face normals.
350pub fn compute_angle_weighted_normals(mesh: &mut IsosurfaceMesh) {
351    let n_verts = mesh.vertices.len() / 3;
352    if n_verts == 0 {
353        return;
354    }
355
356    let mut normals = vec![[0.0f32; 3]; n_verts];
357
358    let n_tris = mesh.indices.len() / 3;
359    for t in 0..n_tris {
360        let idx = [
361            mesh.indices[3 * t] as usize,
362            mesh.indices[3 * t + 1] as usize,
363            mesh.indices[3 * t + 2] as usize,
364        ];
365
366        let v: [[f32; 3]; 3] = [
367            [
368                mesh.vertices[3 * idx[0]],
369                mesh.vertices[3 * idx[0] + 1],
370                mesh.vertices[3 * idx[0] + 2],
371            ],
372            [
373                mesh.vertices[3 * idx[1]],
374                mesh.vertices[3 * idx[1] + 1],
375                mesh.vertices[3 * idx[1] + 2],
376            ],
377            [
378                mesh.vertices[3 * idx[2]],
379                mesh.vertices[3 * idx[2] + 1],
380                mesh.vertices[3 * idx[2] + 2],
381            ],
382        ];
383
384        // Face normal
385        let e1 = [v[1][0] - v[0][0], v[1][1] - v[0][1], v[1][2] - v[0][2]];
386        let e2 = [v[2][0] - v[0][0], v[2][1] - v[0][1], v[2][2] - v[0][2]];
387        let face_normal = [
388            e1[1] * e2[2] - e1[2] * e2[1],
389            e1[2] * e2[0] - e1[0] * e2[2],
390            e1[0] * e2[1] - e1[1] * e2[0],
391        ];
392
393        // For each vertex in the triangle, compute the angle at that vertex
394        for vi in 0..3 {
395            let i0 = vi;
396            let i1 = (vi + 1) % 3;
397            let i2 = (vi + 2) % 3;
398
399            let a = [
400                v[i1][0] - v[i0][0],
401                v[i1][1] - v[i0][1],
402                v[i1][2] - v[i0][2],
403            ];
404            let b = [
405                v[i2][0] - v[i0][0],
406                v[i2][1] - v[i0][1],
407                v[i2][2] - v[i0][2],
408            ];
409
410            let dot = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
411            let la = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt();
412            let lb = (b[0] * b[0] + b[1] * b[1] + b[2] * b[2]).sqrt();
413
414            let cos_angle = if la > 1e-8 && lb > 1e-8 {
415                (dot / (la * lb)).clamp(-1.0, 1.0)
416            } else {
417                0.0
418            };
419            let angle = cos_angle.acos();
420
421            // Weight the face normal by the angle at this vertex
422            normals[idx[vi]][0] += angle * face_normal[0];
423            normals[idx[vi]][1] += angle * face_normal[1];
424            normals[idx[vi]][2] += angle * face_normal[2];
425        }
426    }
427
428    // Normalize all vertex normals
429    for i in 0..n_verts {
430        let n = &mut normals[i];
431        let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
432        if len > 1e-8 {
433            n[0] /= len;
434            n[1] /= len;
435            n[2] /= len;
436        } else {
437            n[0] = 0.0;
438            n[1] = 0.0;
439            n[2] = 1.0;
440        }
441        mesh.normals[3 * i] = n[0];
442        mesh.normals[3 * i + 1] = n[1];
443        mesh.normals[3 * i + 2] = n[2];
444    }
445}
446
447/// Ensure consistent outward-facing normal orientation.
448///
449/// Flips normals that point inward (based on the overall winding of the mesh).
450pub fn flip_normals_outward(mesh: &mut IsosurfaceMesh) {
451    let n_verts = mesh.vertices.len() / 3;
452    if n_verts == 0 {
453        return;
454    }
455
456    // Compute centroid of all vertices
457    let mut cx = 0.0f32;
458    let mut cy = 0.0f32;
459    let mut cz = 0.0f32;
460    for i in 0..n_verts {
461        cx += mesh.vertices[3 * i];
462        cy += mesh.vertices[3 * i + 1];
463        cz += mesh.vertices[3 * i + 2];
464    }
465    cx /= n_verts as f32;
466    cy /= n_verts as f32;
467    cz /= n_verts as f32;
468
469    // For each vertex, check if normal points away from centroid
470    for i in 0..n_verts {
471        let vx = mesh.vertices[3 * i] - cx;
472        let vy = mesh.vertices[3 * i + 1] - cy;
473        let vz = mesh.vertices[3 * i + 2] - cz;
474
475        let dot =
476            vx * mesh.normals[3 * i] + vy * mesh.normals[3 * i + 1] + vz * mesh.normals[3 * i + 2];
477
478        // If normal points inward (dot < 0), flip it
479        if dot < 0.0 {
480            mesh.normals[3 * i] = -mesh.normals[3 * i];
481            mesh.normals[3 * i + 1] = -mesh.normals[3 * i + 1];
482            mesh.normals[3 * i + 2] = -mesh.normals[3 * i + 2];
483        }
484    }
485}
486
487/// Export mesh data in WASM-ready format: interleaved position+normal Float32Arrays.
488///
489/// Returns a flat array of [x, y, z, nx, ny, nz, x, y, z, nx, ny, nz, ...].
490pub fn mesh_to_interleaved(mesh: &IsosurfaceMesh) -> Vec<f32> {
491    let n_verts = mesh.vertices.len() / 3;
492    let mut interleaved = Vec::with_capacity(n_verts * 6);
493    for i in 0..n_verts {
494        interleaved.push(mesh.vertices[3 * i]);
495        interleaved.push(mesh.vertices[3 * i + 1]);
496        interleaved.push(mesh.vertices[3 * i + 2]);
497        interleaved.push(mesh.normals[3 * i]);
498        interleaved.push(mesh.normals[3 * i + 1]);
499        interleaved.push(mesh.normals[3 * i + 2]);
500    }
501    interleaved
502}
503
504// ─── Marching Cubes lookup tables ────────────────────────────────────────────
505// Standard MC edge and triangle tables (Lorensen & Cline, 1987).
506
507#[rustfmt::skip]
508static EDGE_TABLE: [u16; 256] = [
509    0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
510    0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
511    0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
512    0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
513    0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c,
514    0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
515    0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac,
516    0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
517    0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c,
518    0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
519    0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc,
520    0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
521    0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c,
522    0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
523    0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc,
524    0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
525    0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
526    0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
527    0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
528    0x15c, 0x055, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
529    0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
530    0x2fc, 0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
531    0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
532    0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460,
533    0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
534    0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
535    0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
536    0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230,
537    0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
538    0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190,
539    0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
540    0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000,
541];
542
543#[rustfmt::skip]
544static TRI_TABLE: [[i8; 16]; 256] = [
545    [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
546    [ 0, 8, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
547    [ 0, 1, 9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
548    [ 1, 8, 3, 9, 8, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
549    [ 1, 2,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
550    [ 0, 8, 3, 1, 2,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
551    [ 9, 2,10, 0, 2, 9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
552    [ 2, 8, 3, 2,10, 8,10, 9, 8,-1,-1,-1,-1,-1,-1,-1],
553    [ 3,11, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
554    [ 0,11, 2, 8,11, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
555    [ 1, 9, 0, 2, 3,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
556    [ 1,11, 2, 1, 9,11, 9, 8,11,-1,-1,-1,-1,-1,-1,-1],
557    [ 3,10, 1,11,10, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
558    [ 0,10, 1, 0, 8,10, 8,11,10,-1,-1,-1,-1,-1,-1,-1],
559    [ 3, 9, 0, 3,11, 9,11,10, 9,-1,-1,-1,-1,-1,-1,-1],
560    [ 9, 8,10,10, 8,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
561    [ 4, 7, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
562    [ 4, 3, 0, 7, 3, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
563    [ 0, 1, 9, 8, 4, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
564    [ 4, 1, 9, 4, 7, 1, 7, 3, 1,-1,-1,-1,-1,-1,-1,-1],
565    [ 1, 2,10, 8, 4, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
566    [ 3, 4, 7, 3, 0, 4, 1, 2,10,-1,-1,-1,-1,-1,-1,-1],
567    [ 9, 2,10, 9, 0, 2, 8, 4, 7,-1,-1,-1,-1,-1,-1,-1],
568    [ 2,10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4,-1,-1,-1,-1],
569    [ 8, 4, 7, 3,11, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
570    [11, 4, 7,11, 2, 4, 2, 0, 4,-1,-1,-1,-1,-1,-1,-1],
571    [ 9, 0, 1, 8, 4, 7, 2, 3,11,-1,-1,-1,-1,-1,-1,-1],
572    [ 4, 7,11, 9, 4,11, 9,11, 2, 9, 2, 1,-1,-1,-1,-1],
573    [ 3,10, 1, 3,11,10, 7, 8, 4,-1,-1,-1,-1,-1,-1,-1],
574    [ 1,11,10, 1, 4,11, 1, 0, 4, 7,11, 4,-1,-1,-1,-1],
575    [ 4, 7, 8, 9, 0,11, 9,11,10,11, 0, 3,-1,-1,-1,-1],
576    [ 4, 7,11, 4,11, 9, 9,11,10,-1,-1,-1,-1,-1,-1,-1],
577    [ 9, 5, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
578    [ 9, 5, 4, 0, 8, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
579    [ 0, 5, 4, 1, 5, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
580    [ 8, 5, 4, 8, 3, 5, 3, 1, 5,-1,-1,-1,-1,-1,-1,-1],
581    [ 1, 2,10, 9, 5, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
582    [ 3, 0, 8, 1, 2,10, 4, 9, 5,-1,-1,-1,-1,-1,-1,-1],
583    [ 5, 2,10, 5, 4, 2, 4, 0, 2,-1,-1,-1,-1,-1,-1,-1],
584    [ 2,10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8,-1,-1,-1,-1],
585    [ 9, 5, 4, 2, 3,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
586    [ 0,11, 2, 0, 8,11, 4, 9, 5,-1,-1,-1,-1,-1,-1,-1],
587    [ 0, 5, 4, 0, 1, 5, 2, 3,11,-1,-1,-1,-1,-1,-1,-1],
588    [ 2, 1, 5, 2, 5, 8, 2, 8,11, 4, 8, 5,-1,-1,-1,-1],
589    [10, 3,11,10, 1, 3, 9, 5, 4,-1,-1,-1,-1,-1,-1,-1],
590    [ 4, 9, 5, 0, 8, 1, 8,10, 1, 8,11,10,-1,-1,-1,-1],
591    [ 5, 4, 0, 5, 0,11, 5,11,10,11, 0, 3,-1,-1,-1,-1],
592    [ 5, 4, 8, 5, 8,10,10, 8,11,-1,-1,-1,-1,-1,-1,-1],
593    [ 9, 7, 8, 5, 7, 9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
594    [ 9, 3, 0, 9, 5, 3, 5, 7, 3,-1,-1,-1,-1,-1,-1,-1],
595    [ 0, 7, 8, 0, 1, 7, 1, 5, 7,-1,-1,-1,-1,-1,-1,-1],
596    [ 1, 5, 3, 3, 5, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
597    [ 9, 7, 8, 9, 5, 7,10, 1, 2,-1,-1,-1,-1,-1,-1,-1],
598    [10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3,-1,-1,-1,-1],
599    [ 8, 0, 2, 8, 2, 5, 8, 5, 7,10, 5, 2,-1,-1,-1,-1],
600    [ 2,10, 5, 2, 5, 3, 3, 5, 7,-1,-1,-1,-1,-1,-1,-1],
601    [ 7, 9, 5, 7, 8, 9, 3,11, 2,-1,-1,-1,-1,-1,-1,-1],
602    [ 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7,11,-1,-1,-1,-1],
603    [ 2, 3,11, 0, 1, 8, 1, 7, 8, 1, 5, 7,-1,-1,-1,-1],
604    [11, 2, 1,11, 1, 7, 7, 1, 5,-1,-1,-1,-1,-1,-1,-1],
605    [ 9, 5, 8, 8, 5, 7,10, 1, 3,10, 3,11,-1,-1,-1,-1],
606    [ 5, 7, 0, 5, 0, 9, 7,11, 0, 1, 0,10,11,10, 0,-1],
607    [11,10, 0,11, 0, 3,10, 5, 0, 8, 0, 7, 5, 7, 0,-1],
608    [11,10, 5, 7,11, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
609    [10, 6, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
610    [ 0, 8, 3, 5,10, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
611    [ 9, 0, 1, 5,10, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
612    [ 1, 8, 3, 1, 9, 8, 5,10, 6,-1,-1,-1,-1,-1,-1,-1],
613    [ 1, 6, 5, 2, 6, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
614    [ 1, 6, 5, 1, 2, 6, 3, 0, 8,-1,-1,-1,-1,-1,-1,-1],
615    [ 9, 6, 5, 9, 0, 6, 0, 2, 6,-1,-1,-1,-1,-1,-1,-1],
616    [ 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8,-1,-1,-1,-1],
617    [ 2, 3,11,10, 6, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
618    [11, 0, 8,11, 2, 0,10, 6, 5,-1,-1,-1,-1,-1,-1,-1],
619    [ 0, 1, 9, 2, 3,11, 5,10, 6,-1,-1,-1,-1,-1,-1,-1],
620    [ 5,10, 6, 1, 9, 2, 9,11, 2, 9, 8,11,-1,-1,-1,-1],
621    [ 6, 3,11, 6, 5, 3, 5, 1, 3,-1,-1,-1,-1,-1,-1,-1],
622    [ 0, 8,11, 0,11, 5, 0, 5, 1, 5,11, 6,-1,-1,-1,-1],
623    [ 3,11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9,-1,-1,-1,-1],
624    [ 6, 5, 9, 6, 9,11,11, 9, 8,-1,-1,-1,-1,-1,-1,-1],
625    [ 5,10, 6, 4, 7, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
626    [ 4, 3, 0, 4, 7, 3, 6, 5,10,-1,-1,-1,-1,-1,-1,-1],
627    [ 1, 9, 0, 5,10, 6, 8, 4, 7,-1,-1,-1,-1,-1,-1,-1],
628    [10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4,-1,-1,-1,-1],
629    [ 6, 1, 2, 6, 5, 1, 4, 7, 8,-1,-1,-1,-1,-1,-1,-1],
630    [ 1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7,-1,-1,-1,-1],
631    [ 8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6,-1,-1,-1,-1],
632    [ 7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9,-1],
633    [ 3,11, 2, 7, 8, 4,10, 6, 5,-1,-1,-1,-1,-1,-1,-1],
634    [ 5,10, 6, 4, 7, 2, 4, 2, 0, 2, 7,11,-1,-1,-1,-1],
635    [ 0, 1, 9, 4, 7, 8, 2, 3,11, 5,10, 6,-1,-1,-1,-1],
636    [ 9, 2, 1, 9,11, 2, 9, 4,11, 7,11, 4, 5,10, 6,-1],
637    [ 8, 4, 7, 3,11, 5, 3, 5, 1, 5,11, 6,-1,-1,-1,-1],
638    [ 5, 1,11, 5,11, 6, 1, 0,11, 7,11, 4, 0, 4,11,-1],
639    [ 0, 5, 9, 0, 6, 5, 0, 3, 6,11, 6, 3, 8, 4, 7,-1],
640    [ 6, 5, 9, 6, 9,11, 4, 7, 9, 7,11, 9,-1,-1,-1,-1],
641    [10, 4, 9, 6, 4,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
642    [ 4,10, 6, 4, 9,10, 0, 8, 3,-1,-1,-1,-1,-1,-1,-1],
643    [10, 0, 1,10, 6, 0, 6, 4, 0,-1,-1,-1,-1,-1,-1,-1],
644    [ 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1,10,-1,-1,-1,-1],
645    [ 1, 4, 9, 1, 2, 4, 2, 6, 4,-1,-1,-1,-1,-1,-1,-1],
646    [ 3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4,-1,-1,-1,-1],
647    [ 0, 2, 4, 4, 2, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
648    [ 8, 3, 2, 8, 2, 4, 4, 2, 6,-1,-1,-1,-1,-1,-1,-1],
649    [10, 4, 9,10, 6, 4,11, 2, 3,-1,-1,-1,-1,-1,-1,-1],
650    [ 0, 8, 2, 2, 8,11, 4, 9,10, 4,10, 6,-1,-1,-1,-1],
651    [ 3,11, 2, 0, 1, 6, 0, 6, 4, 6, 1,10,-1,-1,-1,-1],
652    [ 6, 4, 1, 6, 1,10, 4, 8, 1, 2, 1,11, 8,11, 1,-1],
653    [ 9, 6, 4, 9, 3, 6, 9, 1, 3,11, 6, 3,-1,-1,-1,-1],
654    [ 8,11, 1, 8, 1, 0,11, 6, 1, 9, 1, 4, 6, 4, 1,-1],
655    [ 3,11, 6, 3, 6, 0, 0, 6, 4,-1,-1,-1,-1,-1,-1,-1],
656    [ 6, 4, 8,11, 6, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
657    [ 7,10, 6, 7, 8,10, 8, 9,10,-1,-1,-1,-1,-1,-1,-1],
658    [ 0, 7, 3, 0,10, 7, 0, 9,10, 6, 7,10,-1,-1,-1,-1],
659    [10, 6, 7, 1,10, 7, 1, 7, 8, 1, 8, 0,-1,-1,-1,-1],
660    [10, 6, 7,10, 7, 1, 1, 7, 3,-1,-1,-1,-1,-1,-1,-1],
661    [ 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7,-1,-1,-1,-1],
662    [ 2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9,-1],
663    [ 7, 8, 0, 7, 0, 6, 6, 0, 2,-1,-1,-1,-1,-1,-1,-1],
664    [ 7, 3, 2, 6, 7, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
665    [ 2, 3,11,10, 6, 8,10, 8, 9, 8, 6, 7,-1,-1,-1,-1],
666    [ 2, 0, 7, 2, 7,11, 0, 9, 7, 6, 7,10, 9,10, 7,-1],
667    [ 1, 8, 0, 1, 7, 8, 1,10, 7, 6, 7,10, 2, 3,11,-1],
668    [11, 2, 1,11, 1, 7,10, 6, 1, 6, 7, 1,-1,-1,-1,-1],
669    [ 8, 9, 6, 8, 6, 7, 9, 1, 6,11, 6, 3, 1, 3, 6,-1],
670    [ 0, 9, 1,11, 6, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
671    [ 7, 8, 0, 7, 0, 6, 3,11, 0,11, 6, 0,-1,-1,-1,-1],
672    [ 7,11, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
673    [ 7, 6,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
674    [ 3, 0, 8,11, 7, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
675    [ 0, 1, 9,11, 7, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
676    [ 8, 1, 9, 8, 3, 1,11, 7, 6,-1,-1,-1,-1,-1,-1,-1],
677    [10, 1, 2, 6,11, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
678    [ 1, 2,10, 3, 0, 8, 6,11, 7,-1,-1,-1,-1,-1,-1,-1],
679    [ 2, 9, 0, 2,10, 9, 6,11, 7,-1,-1,-1,-1,-1,-1,-1],
680    [ 6,11, 7, 2,10, 3,10, 8, 3,10, 9, 8,-1,-1,-1,-1],
681    [ 7, 2, 3, 6, 2, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
682    [ 7, 0, 8, 7, 6, 0, 6, 2, 0,-1,-1,-1,-1,-1,-1,-1],
683    [ 2, 7, 6, 2, 3, 7, 0, 1, 9,-1,-1,-1,-1,-1,-1,-1],
684    [ 1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6,-1,-1,-1,-1],
685    [10, 7, 6,10, 1, 7, 1, 3, 7,-1,-1,-1,-1,-1,-1,-1],
686    [10, 7, 6, 1, 7,10, 1, 8, 7, 1, 0, 8,-1,-1,-1,-1],
687    [ 0, 3, 7, 0, 7,10, 0,10, 9, 6,10, 7,-1,-1,-1,-1],
688    [ 7, 6,10, 7,10, 8, 8,10, 9,-1,-1,-1,-1,-1,-1,-1],
689    [ 6, 8, 4,11, 8, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
690    [ 3, 6,11, 3, 0, 6, 0, 4, 6,-1,-1,-1,-1,-1,-1,-1],
691    [ 8, 6,11, 8, 4, 6, 9, 0, 1,-1,-1,-1,-1,-1,-1,-1],
692    [ 9, 4, 6, 9, 6, 3, 9, 3, 1,11, 3, 6,-1,-1,-1,-1],
693    [ 6, 8, 4, 6,11, 8, 2,10, 1,-1,-1,-1,-1,-1,-1,-1],
694    [ 1, 2,10, 3, 0,11, 0, 6,11, 0, 4, 6,-1,-1,-1,-1],
695    [ 4,11, 8, 4, 6,11, 0, 2, 9, 2,10, 9,-1,-1,-1,-1],
696    [10, 9, 3,10, 3, 2, 9, 4, 3,11, 3, 6, 4, 6, 3,-1],
697    [ 8, 2, 3, 8, 4, 2, 4, 6, 2,-1,-1,-1,-1,-1,-1,-1],
698    [ 0, 4, 2, 4, 6, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
699    [ 1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8,-1,-1,-1,-1],
700    [ 1, 9, 4, 1, 4, 2, 2, 4, 6,-1,-1,-1,-1,-1,-1,-1],
701    [ 8, 1, 3, 8, 6, 1, 8, 4, 6, 6,10, 1,-1,-1,-1,-1],
702    [10, 1, 0,10, 0, 6, 6, 0, 4,-1,-1,-1,-1,-1,-1,-1],
703    [ 4, 6, 3, 4, 3, 8, 6,10, 3, 0, 3, 9,10, 9, 3,-1],
704    [10, 9, 4, 6,10, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
705    [ 4, 9, 5, 7, 6,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
706    [ 0, 8, 3, 4, 9, 5,11, 7, 6,-1,-1,-1,-1,-1,-1,-1],
707    [ 5, 0, 1, 5, 4, 0, 7, 6,11,-1,-1,-1,-1,-1,-1,-1],
708    [11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5,-1,-1,-1,-1],
709    [ 9, 5, 4,10, 1, 2, 7, 6,11,-1,-1,-1,-1,-1,-1,-1],
710    [ 6,11, 7, 1, 2,10, 0, 8, 3, 4, 9, 5,-1,-1,-1,-1],
711    [ 7, 6,11, 5, 4,10, 4, 2,10, 4, 0, 2,-1,-1,-1,-1],
712    [ 3, 4, 8, 3, 5, 4, 3, 2, 5,10, 5, 2,11, 7, 6,-1],
713    [ 7, 2, 3, 7, 6, 2, 5, 4, 9,-1,-1,-1,-1,-1,-1,-1],
714    [ 9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7,-1,-1,-1,-1],
715    [ 3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0,-1,-1,-1,-1],
716    [ 6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8,-1],
717    [ 9, 5, 4,10, 1, 6, 1, 7, 6, 1, 3, 7,-1,-1,-1,-1],
718    [ 1, 6,10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4,-1],
719    [ 4, 0,10, 4,10, 5, 0, 3,10, 6,10, 7, 3, 7,10,-1],
720    [ 7, 6,10, 7,10, 8, 5, 4,10, 4, 8,10,-1,-1,-1,-1],
721    [ 6, 9, 5, 6,11, 9,11, 8, 9,-1,-1,-1,-1,-1,-1,-1],
722    [ 3, 6,11, 0, 6, 3, 0, 5, 6, 0, 9, 5,-1,-1,-1,-1],
723    [ 0,11, 8, 0, 5,11, 0, 1, 5, 5, 6,11,-1,-1,-1,-1],
724    [ 6,11, 3, 6, 3, 5, 5, 3, 1,-1,-1,-1,-1,-1,-1,-1],
725    [ 1, 2,10, 9, 5,11, 9,11, 8,11, 5, 6,-1,-1,-1,-1],
726    [ 0,11, 3, 0, 6,11, 0, 9, 6, 5, 6, 9, 1, 2,10,-1],
727    [11, 8, 5,11, 5, 6, 8, 0, 5,10, 5, 2, 0, 2, 5,-1],
728    [ 6,11, 3, 6, 3, 5, 2,10, 3,10, 5, 3,-1,-1,-1,-1],
729    [ 5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2,-1,-1,-1,-1],
730    [ 9, 5, 6, 9, 6, 0, 0, 6, 2,-1,-1,-1,-1,-1,-1,-1],
731    [ 1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8,-1],
732    [ 1, 5, 6, 2, 1, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
733    [ 1, 3, 6, 1, 6,10, 3, 8, 6, 5, 6, 9, 8, 9, 6,-1],
734    [10, 1, 0,10, 0, 6, 9, 5, 0, 5, 6, 0,-1,-1,-1,-1],
735    [ 0, 3, 8, 5, 6,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
736    [10, 5, 6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
737    [11, 5,10, 7, 5,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
738    [11, 5,10,11, 7, 5, 8, 3, 0,-1,-1,-1,-1,-1,-1,-1],
739    [ 5,11, 7, 5,10,11, 1, 9, 0,-1,-1,-1,-1,-1,-1,-1],
740    [10, 7, 5,10,11, 7, 9, 8, 1, 8, 3, 1,-1,-1,-1,-1],
741    [11, 1, 2,11, 7, 1, 7, 5, 1,-1,-1,-1,-1,-1,-1,-1],
742    [ 0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2,11,-1,-1,-1,-1],
743    [ 9, 7, 5, 9, 2, 7, 9, 0, 2, 2,11, 7,-1,-1,-1,-1],
744    [ 7, 5, 2, 7, 2,11, 5, 9, 2, 3, 2, 8, 9, 8, 2,-1],
745    [ 2, 5,10, 2, 3, 5, 3, 7, 5,-1,-1,-1,-1,-1,-1,-1],
746    [ 8, 2, 0, 8, 5, 2, 8, 7, 5,10, 2, 5,-1,-1,-1,-1],
747    [ 9, 0, 1, 5,10, 3, 5, 3, 7, 3,10, 2,-1,-1,-1,-1],
748    [ 9, 8, 2, 9, 2, 1, 8, 7, 2,10, 2, 5, 7, 5, 2,-1],
749    [ 1, 3, 5, 3, 7, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
750    [ 0, 8, 7, 0, 7, 1, 1, 7, 5,-1,-1,-1,-1,-1,-1,-1],
751    [ 9, 0, 3, 9, 3, 5, 5, 3, 7,-1,-1,-1,-1,-1,-1,-1],
752    [ 9, 8, 7, 5, 9, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
753    [ 5, 8, 4, 5,10, 8,10,11, 8,-1,-1,-1,-1,-1,-1,-1],
754    [ 5, 0, 4, 5,11, 0, 5,10,11,11, 3, 0,-1,-1,-1,-1],
755    [ 0, 1, 9, 8, 4,10, 8,10,11,10, 4, 5,-1,-1,-1,-1],
756    [10,11, 4,10, 4, 5,11, 3, 4, 9, 4, 1, 3, 1, 4,-1],
757    [ 2, 5, 1, 2, 8, 5, 2,11, 8, 4, 5, 8,-1,-1,-1,-1],
758    [ 0, 4,11, 0,11, 3, 4, 5,11, 2,11, 1, 5, 1,11,-1],
759    [ 0, 2, 5, 0, 5, 9, 2,11, 5, 4, 5, 8,11, 8, 5,-1],
760    [ 9, 4, 5, 2,11, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
761    [ 2, 5,10, 3, 5, 2, 3, 4, 5, 3, 8, 4,-1,-1,-1,-1],
762    [ 5,10, 2, 5, 2, 4, 4, 2, 0,-1,-1,-1,-1,-1,-1,-1],
763    [ 3,10, 2, 3, 5,10, 3, 8, 5, 4, 5, 8, 0, 1, 9,-1],
764    [ 5,10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2,-1,-1,-1,-1],
765    [ 8, 4, 5, 8, 5, 3, 3, 5, 1,-1,-1,-1,-1,-1,-1,-1],
766    [ 0, 4, 5, 1, 0, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
767    [ 8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5,-1,-1,-1,-1],
768    [ 9, 4, 5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
769    [ 4,11, 7, 4, 9,11, 9,10,11,-1,-1,-1,-1,-1,-1,-1],
770    [ 0, 8, 3, 4, 9, 7, 9,11, 7, 9,10,11,-1,-1,-1,-1],
771    [ 1,10,11, 1,11, 4, 1, 4, 0, 7, 4,11,-1,-1,-1,-1],
772    [ 3, 1, 4, 3, 4, 8, 1,10, 4, 7, 4,11,10,11, 4,-1],
773    [ 4,11, 7, 9,11, 4, 9, 2,11, 9, 1, 2,-1,-1,-1,-1],
774    [ 9, 7, 4, 9,11, 7, 9, 1,11, 2,11, 1, 0, 8, 3,-1],
775    [11, 7, 4,11, 4, 2, 2, 4, 0,-1,-1,-1,-1,-1,-1,-1],
776    [11, 7, 4,11, 4, 2, 8, 3, 4, 3, 2, 4,-1,-1,-1,-1],
777    [ 2, 9,10, 2, 7, 9, 2, 3, 7, 7, 4, 9,-1,-1,-1,-1],
778    [ 9,10, 7, 9, 7, 4,10, 2, 7, 8, 7, 0, 2, 0, 7,-1],
779    [ 3, 7,10, 3,10, 2, 7, 4,10, 1,10, 0, 4, 0,10,-1],
780    [ 1,10, 2, 8, 7, 4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
781    [ 4, 9, 1, 4, 1, 7, 7, 1, 3,-1,-1,-1,-1,-1,-1,-1],
782    [ 4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1,-1,-1,-1,-1],
783    [ 4, 0, 3, 7, 4, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
784    [ 4, 8, 7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
785    [ 9,10, 8,10,11, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
786    [ 3, 0, 9, 3, 9,11,11, 9,10,-1,-1,-1,-1,-1,-1,-1],
787    [ 0, 1,10, 0,10, 8, 8,10,11,-1,-1,-1,-1,-1,-1,-1],
788    [ 3, 1,10,11, 3,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
789    [ 1, 2,11, 1,11, 9, 9,11, 8,-1,-1,-1,-1,-1,-1,-1],
790    [ 3, 0, 9, 3, 9,11, 1, 2, 9, 2,11, 9,-1,-1,-1,-1],
791    [ 0, 2,11, 8, 0,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
792    [ 3, 2,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
793    [ 2, 3, 8, 2, 8,10,10, 8, 9,-1,-1,-1,-1,-1,-1,-1],
794    [ 9,10, 2, 0, 9, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
795    [ 2, 3, 8, 2, 8,10, 0, 1, 8, 1,10, 8,-1,-1,-1,-1],
796    [ 1,10, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
797    [ 1, 3, 8, 9, 1, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
798    [ 0, 9, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
799    [ 0, 3, 8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
800    [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
801];
802
803#[cfg(test)]
804mod tests {
805    use super::super::volume::VolumetricGrid;
806    use super::*;
807
808    fn sphere_grid(radius: f64, spacing: f64) -> VolumetricGrid {
809        let extent = radius + 1.0;
810        let n = (2.0 * extent / spacing).ceil() as usize + 1;
811        let origin = [-extent, -extent, -extent];
812        let mut values = vec![0.0f32; n * n * n];
813        for ix in 0..n {
814            for iy in 0..n {
815                for iz in 0..n {
816                    let x = origin[0] + ix as f64 * spacing;
817                    let y = origin[1] + iy as f64 * spacing;
818                    let z = origin[2] + iz as f64 * spacing;
819                    let r = (x * x + y * y + z * z).sqrt();
820                    values[ix * n * n + iy * n + iz] = (radius - r) as f32;
821                }
822            }
823        }
824        VolumetricGrid {
825            origin,
826            spacing,
827            dims: [n, n, n],
828            values,
829        }
830    }
831
832    #[test]
833    fn test_marching_cubes_sphere() {
834        let grid = sphere_grid(2.0, 0.3);
835        let mesh = marching_cubes(&grid, 0.0);
836        assert!(
837            mesh.num_triangles > 50,
838            "sphere should produce many triangles, got {}",
839            mesh.num_triangles,
840        );
841        // All indices should be valid
842        let nv = mesh.num_vertices() as u32;
843        for &idx in &mesh.indices {
844            assert!(idx < nv, "index {} out of bound (nv={})", idx, nv);
845        }
846    }
847
848    #[test]
849    fn test_marching_cubes_empty() {
850        let grid = VolumetricGrid {
851            origin: [0.0; 3],
852            spacing: 1.0,
853            dims: [3, 3, 3],
854            values: vec![0.0f32; 27],
855        };
856        let mesh = marching_cubes(&grid, 1.0);
857        assert_eq!(mesh.num_triangles, 0);
858    }
859
860    #[test]
861    fn test_normals_same_count_as_vertices() {
862        let grid = sphere_grid(2.0, 0.3);
863        let mesh = marching_cubes(&grid, 0.0);
864        assert_eq!(mesh.vertices.len(), mesh.normals.len());
865    }
866
867    #[test]
868    fn test_dual_phase_isosurface() {
869        // Scalar field: positive inside sphere -> positive lobe,
870        // negative outside -> negative lobe (with offset)
871        let grid = sphere_grid(2.0, 0.3);
872        let dual = marching_cubes_dual(&grid, 0.5);
873        assert!(
874            dual.positive.num_triangles > 0,
875            "positive lobe should have triangles"
876        );
877        // The negative phase at -0.5 also produces a surface
878        assert!(
879            dual.negative.num_triangles > 0,
880            "negative lobe should have triangles"
881        );
882    }
883
884    #[test]
885    fn test_simplify_mesh_welds_vertices() {
886        let grid = sphere_grid(2.0, 0.3);
887        let mesh = marching_cubes(&grid, 0.0);
888        let before = mesh.num_vertices();
889        let simplified = simplify_mesh(&mesh, 0.01);
890        // Welding should reduce vertex count
891        assert!(
892            simplified.num_vertices() <= before,
893            "simplified should have <= vertices: {} vs {}",
894            simplified.num_vertices(),
895            before
896        );
897        assert!(simplified.num_triangles > 0, "should still have triangles");
898    }
899
900    #[test]
901    fn test_angle_weighted_normals_sphere() {
902        let grid = sphere_grid(2.0, 0.3);
903        let mut mesh = marching_cubes(&grid, 0.0);
904        compute_angle_weighted_normals(&mut mesh);
905
906        // For a sphere centered at origin, normals should point radially outward
907        // Check cosine similarity between vertex position (= radial direction) and normal
908        let n_verts = mesh.num_vertices();
909        let mut cos_sum = 0.0f64;
910        let mut count = 0;
911        for i in 0..n_verts {
912            let vx = mesh.vertices[3 * i] as f64;
913            let vy = mesh.vertices[3 * i + 1] as f64;
914            let vz = mesh.vertices[3 * i + 2] as f64;
915            let vlen = (vx * vx + vy * vy + vz * vz).sqrt();
916            if vlen < 1e-6 {
917                continue;
918            }
919            let nx = mesh.normals[3 * i] as f64;
920            let ny = mesh.normals[3 * i + 1] as f64;
921            let nz = mesh.normals[3 * i + 2] as f64;
922            let nlen = (nx * nx + ny * ny + nz * nz).sqrt();
923            if nlen < 1e-6 {
924                continue;
925            }
926            let cos = (vx * nx + vy * ny + vz * nz) / (vlen * nlen);
927            cos_sum += cos.abs();
928            count += 1;
929        }
930        let avg_cos = cos_sum / count as f64;
931        assert!(
932            avg_cos > 0.90,
933            "angle-weighted sphere normals should align with radial direction, got avg cos = {:.3}",
934            avg_cos
935        );
936    }
937
938    #[test]
939    fn test_flip_normals_outward_sphere() {
940        let grid = sphere_grid(2.0, 0.3);
941        let mut mesh = marching_cubes(&grid, 0.0);
942        flip_normals_outward(&mut mesh);
943
944        // After flipping, all normals should point away from centroid (≈ origin)
945        let n_verts = mesh.num_vertices();
946        let mut inward = 0;
947        for i in 0..n_verts {
948            let dot = mesh.vertices[3 * i] * mesh.normals[3 * i]
949                + mesh.vertices[3 * i + 1] * mesh.normals[3 * i + 1]
950                + mesh.vertices[3 * i + 2] * mesh.normals[3 * i + 2];
951            if dot < 0.0 {
952                inward += 1;
953            }
954        }
955        assert_eq!(inward, 0, "no normals should point inward after flip");
956    }
957
958    #[test]
959    fn test_mesh_to_interleaved() {
960        let grid = sphere_grid(2.0, 0.3);
961        let mesh = marching_cubes(&grid, 0.0);
962        let interleaved = mesh_to_interleaved(&mesh);
963        let n_verts = mesh.num_vertices();
964        assert_eq!(interleaved.len(), n_verts * 6, "6 floats per vertex");
965
966        // Verify interleaving
967        for i in 0..n_verts.min(10) {
968            assert_eq!(interleaved[6 * i], mesh.vertices[3 * i]);
969            assert_eq!(interleaved[6 * i + 1], mesh.vertices[3 * i + 1]);
970            assert_eq!(interleaved[6 * i + 2], mesh.vertices[3 * i + 2]);
971            assert_eq!(interleaved[6 * i + 3], mesh.normals[3 * i]);
972            assert_eq!(interleaved[6 * i + 4], mesh.normals[3 * i + 1]);
973            assert_eq!(interleaved[6 * i + 5], mesh.normals[3 * i + 2]);
974        }
975    }
976}
sci_form/eht/marching_cubes.rs

sci_form/eht/
marching_cubes.rs
