qualia_core_db/specialized_libs/computational_geometry/
isosurface.rs

1//! P6.3 — Isosurfacing / dual-contouring over scalar fields on the `.10d` grid.
2//!
3//! Marching cubes: extract a triangle mesh from a scalar field sampled on a
4//! regular 3D grid at a given isolevel. The algorithm classifies each grid
5//! cell by the sign of the field at its 8 corners, then generates triangles
6//! from a precomputed table of 256 cases.
7//!
8//! ## Determinism
9//!
10//! The output is deterministic: cells are processed in (x, y, z) order,
11//! vertices within each cell are generated in canonical edge order, and
12//! ties in the interpolation are resolved by the lower-index corner.
13//! Identical input → bit-identical output.
14//!
15//! ## Zero heap
16//!
17//! All hot-path functions use caller-supplied buffers. The grid is passed
18//! as a flat slice with explicit dimensions.
19
20use super::primitives::Point3;
21
22// ───────────────────────────────────────────────────────────────────────────
23//  Errors
24// ───────────────────────────────────────────────────────────────────────────
25
26/// Isosurface extraction error.
27#[derive(Debug, Clone, Copy, PartialEq, Eq)]
28pub enum IsosurfaceError {
29    /// Grid dimensions are zero.
30    EmptyGrid,
31    /// Grid slice doesn't match `nx * ny * nz`.
32    GridSizeMismatch { expected: usize, got: usize },
33    /// Output buffer too small.
34    BufferTooSmall { needed: usize, have: usize },
35}
36
37impl core::fmt::Display for IsosurfaceError {
38    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
39        match self {
40            Self::EmptyGrid => write!(f, "isosurface: empty grid"),
41            Self::GridSizeMismatch { expected, got } => {
42                write!(
43                    f,
44                    "isosurface: grid size mismatch, expected {expected}, got {got}"
45                )
46            }
47            Self::BufferTooSmall { needed, have } => {
48                write!(
49                    f,
50                    "isosurface: buffer too small, need {needed}, have {have}"
51                )
52            }
53        }
54    }
55}
56
57impl std::error::Error for IsosurfaceError {}
58
59// ───────────────────────────────────────────────────────────────────────────
60//  Marching cubes tables
61// ───────────────────────────────────────────────────────────────────────────
62
63/// Edge table: for each of 256 cube configurations, a 12-bit mask indicating
64/// which edges are intersected by the isosurface.
65///
66/// Edge numbering (standard marching cubes convention):
67/// ```text
68///    4------5
69///   /|     /|
70///  7------6 |
71///  | |    | |
72///  | 0----|-1
73///  |/     |/
74///  3------2
75///
76/// Edges:
77///  0: 0-1  1: 1-2  2: 2-3  3: 3-0
78///  4: 4-5  5: 5-6  6: 6-7  7: 7-4
79///  8: 0-4  9: 1-5  10: 2-6  11: 3-7
80/// ```
81const EDGE_TABLE: [u16; 256] = {
82    let mut table = [0u16; 256];
83    let mut i = 0;
84    while i < 256 {
85        let mut bits = 0u16;
86        // Edge 0: between corner 0 and 1
87        if (i & 1) != ((i >> 1) & 1) {
88            bits |= 1 << 0;
89        }
90        // Edge 1: between corner 1 and 2
91        if ((i >> 1) & 1) != ((i >> 2) & 1) {
92            bits |= 1 << 1;
93        }
94        // Edge 2: between corner 2 and 3
95        if ((i >> 2) & 1) != ((i >> 3) & 1) {
96            bits |= 1 << 2;
97        }
98        // Edge 3: between corner 3 and 0
99        if ((i >> 3) & 1) != (i & 1) {
100            bits |= 1 << 3;
101        }
102        // Edge 4: between corner 4 and 5
103        if ((i >> 4) & 1) != ((i >> 5) & 1) {
104            bits |= 1 << 4;
105        }
106        // Edge 5: between corner 5 and 6
107        if ((i >> 5) & 1) != ((i >> 6) & 1) {
108            bits |= 1 << 5;
109        }
110        // Edge 6: between corner 6 and 7
111        if ((i >> 6) & 1) != ((i >> 7) & 1) {
112            bits |= 1 << 6;
113        }
114        // Edge 7: between corner 7 and 4
115        if ((i >> 7) & 1) != ((i >> 4) & 1) {
116            bits |= 1 << 7;
117        }
118        // Edge 8: between corner 0 and 4
119        if (i & 1) != ((i >> 4) & 1) {
120            bits |= 1 << 8;
121        }
122        // Edge 9: between corner 1 and 5
123        if ((i >> 1) & 1) != ((i >> 5) & 1) {
124            bits |= 1 << 9;
125        }
126        // Edge 10: between corner 2 and 6
127        if ((i >> 2) & 1) != ((i >> 6) & 1) {
128            bits |= 1 << 10;
129        }
130        // Edge 11: between corner 3 and 7
131        if ((i >> 3) & 1) != ((i >> 7) & 1) {
132            bits |= 1 << 11;
133        }
134        table[i] = bits;
135        i += 1;
136    }
137    table
138};
139
140/// Triangle table: for each of the 256 cube configurations, the list of edge
141/// indices (in groups of 3) forming the output triangles, `-1`-terminated.
142///
143/// This is the **canonical marching-cubes triangulation table** (Lorensen &
144/// Cline 1987; the widely-reproduced Cory Bloyd / Paul Bourke public-domain
145/// form). It is the correct, ambiguity-resolved triangulation for every one of
146/// the 256 corner sign configurations — NOT a fan approximation. The edge
147/// numbering (0-11), corner numbering (0-7), and the `cube_idx` bit convention
148/// (bit `c` set iff corner `c` is below the isolevel) all match this file's
149/// [`EDGE_CORNERS`] / [`CORNER_OFFSETS`] and the classifier below, so the table
150/// drops in directly. It is mathematical/algorithmic data (a lookup table for a
151/// public-domain algorithm), consulted from the algorithm's definition — no
152/// GPL/LGPL source is used or derived.
153///
154/// Correctness is guarded two ways by the test suite: (1) `tri_table_edges_match_edge_table`
155/// asserts, for all 256 cases, that the *set* of edges used by the triangles
156/// equals the independently-computed [`EDGE_TABLE`] crossing set (catches any
157/// wrong/missing edge index); (2) `sphere_isosurface_is_closed_manifold` merges
158/// coincident vertices and asserts the extracted closed-level-set surface is a
159/// watertight 2-manifold (every edge shared by exactly two triangles — catches
160/// any wrong triangulation grouping).
161const TRI_TABLE: [[i8; 16]; 256] = build_tri_table();
162
163/// Build the fixed-width `[[i8; 16]; 256]` table by right-padding each
164/// variable-length row of meaningful edge indices with `-1`. Writing only the
165/// meaningful edges (no hand-typed padding, no counting to 16) removes a whole
166/// class of transcription error from the 256-row constant.
167const fn build_tri_table() -> [[i8; 16]; 256] {
168    let mut table = [[-1i8; 16]; 256];
169    let mut i = 0;
170    while i < 256 {
171        let row = TRI_TABLE_RAW[i];
172        let mut j = 0;
173        while j < row.len() {
174            table[i][j] = row[j];
175            j += 1;
176        }
177        i += 1;
178    }
179    table
180}
181
182/// The canonical marching-cubes triangulation, one slice of edge-index triples
183/// per cube configuration (empty = no surface in that cell). Padded to the
184/// `-1`-terminated fixed form by [`build_tri_table`].
185#[rustfmt::skip]
186const TRI_TABLE_RAW: [&[i8]; 256] = [
187    &[],
188    &[0, 8, 3],
189    &[0, 1, 9],
190    &[1, 8, 3, 9, 8, 1],
191    &[1, 2, 10],
192    &[0, 8, 3, 1, 2, 10],
193    &[9, 2, 10, 0, 2, 9],
194    &[2, 8, 3, 2, 10, 8, 10, 9, 8],
195    &[3, 11, 2],
196    &[0, 11, 2, 8, 11, 0],
197    &[1, 9, 0, 2, 3, 11],
198    &[1, 11, 2, 1, 9, 11, 9, 8, 11],
199    &[3, 10, 1, 11, 10, 3],
200    &[0, 10, 1, 0, 8, 10, 8, 11, 10],
201    &[3, 9, 0, 3, 11, 9, 11, 10, 9],
202    &[9, 8, 10, 10, 8, 11],
203    &[4, 7, 8],
204    &[4, 3, 0, 7, 3, 4],
205    &[0, 1, 9, 8, 4, 7],
206    &[4, 1, 9, 4, 7, 1, 7, 3, 1],
207    &[1, 2, 10, 8, 4, 7],
208    &[3, 4, 7, 3, 0, 4, 1, 2, 10],
209    &[9, 2, 10, 9, 0, 2, 8, 4, 7],
210    &[2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4],
211    &[8, 4, 7, 3, 11, 2],
212    &[11, 4, 7, 11, 2, 4, 2, 0, 4],
213    &[9, 0, 1, 8, 4, 7, 2, 3, 11],
214    &[4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1],
215    &[3, 10, 1, 3, 11, 10, 7, 8, 4],
216    &[1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4],
217    &[4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3],
218    &[4, 7, 11, 4, 11, 9, 9, 11, 10],
219    &[9, 5, 4],
220    &[9, 5, 4, 0, 8, 3],
221    &[0, 5, 4, 1, 5, 0],
222    &[8, 5, 4, 8, 3, 5, 3, 1, 5],
223    &[1, 2, 10, 9, 5, 4],
224    &[3, 0, 8, 1, 2, 10, 4, 9, 5],
225    &[5, 2, 10, 5, 4, 2, 4, 0, 2],
226    &[2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8],
227    &[9, 5, 4, 2, 3, 11],
228    &[0, 11, 2, 0, 8, 11, 4, 9, 5],
229    &[0, 5, 4, 0, 1, 5, 2, 3, 11],
230    &[2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5],
231    &[10, 3, 11, 10, 1, 3, 9, 5, 4],
232    &[4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10],
233    &[5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3],
234    &[5, 4, 8, 5, 8, 10, 10, 8, 11],
235    &[9, 7, 8, 5, 7, 9],
236    &[9, 3, 0, 9, 5, 3, 5, 7, 3],
237    &[0, 7, 8, 0, 1, 7, 1, 5, 7],
238    &[1, 5, 3, 3, 5, 7],
239    &[9, 7, 8, 9, 5, 7, 10, 1, 2],
240    &[10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3],
241    &[8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2],
242    &[2, 10, 5, 2, 5, 3, 3, 5, 7],
243    &[7, 9, 5, 7, 8, 9, 3, 11, 2],
244    &[9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11],
245    &[2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7],
246    &[11, 2, 1, 11, 1, 7, 7, 1, 5],
247    &[9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11],
248    &[5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0],
249    &[11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0],
250    &[11, 10, 5, 7, 11, 5],
251    &[10, 6, 5],
252    &[0, 8, 3, 5, 10, 6],
253    &[9, 0, 1, 5, 10, 6],
254    &[1, 8, 3, 1, 9, 8, 5, 10, 6],
255    &[1, 6, 5, 2, 6, 1],
256    &[1, 6, 5, 1, 2, 6, 3, 0, 8],
257    &[9, 6, 5, 9, 0, 6, 0, 2, 6],
258    &[5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8],
259    &[2, 3, 11, 10, 6, 5],
260    &[11, 0, 8, 11, 2, 0, 10, 6, 5],
261    &[0, 1, 9, 2, 3, 11, 5, 10, 6],
262    &[5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11],
263    &[6, 3, 11, 6, 5, 3, 5, 1, 3],
264    &[0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6],
265    &[3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9],
266    &[6, 5, 9, 6, 9, 11, 11, 9, 8],
267    &[5, 10, 6, 4, 7, 8],
268    &[4, 3, 0, 4, 7, 3, 6, 5, 10],
269    &[1, 9, 0, 5, 10, 6, 8, 4, 7],
270    &[10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4],
271    &[6, 1, 2, 6, 5, 1, 4, 7, 8],
272    &[1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7],
273    &[8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6],
274    &[7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9],
275    &[3, 11, 2, 7, 8, 4, 10, 6, 5],
276    &[5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11],
277    &[0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6],
278    &[9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6],
279    &[8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6],
280    &[5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11],
281    &[0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7],
282    &[6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9],
283    &[10, 4, 9, 6, 4, 10],
284    &[4, 10, 6, 4, 9, 10, 0, 8, 3],
285    &[10, 0, 1, 10, 6, 0, 6, 4, 0],
286    &[8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10],
287    &[1, 4, 9, 1, 2, 4, 2, 6, 4],
288    &[3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4],
289    &[0, 2, 4, 4, 2, 6],
290    &[8, 3, 2, 8, 2, 4, 4, 2, 6],
291    &[10, 4, 9, 10, 6, 4, 11, 2, 3],
292    &[0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6],
293    &[3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10],
294    &[6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1],
295    &[9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3],
296    &[8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1],
297    &[3, 11, 6, 3, 6, 0, 0, 6, 4],
298    &[6, 4, 8, 11, 6, 8],
299    &[7, 10, 6, 7, 8, 10, 8, 9, 10],
300    &[0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10],
301    &[10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0],
302    &[10, 6, 7, 10, 7, 1, 1, 7, 3],
303    &[1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7],
304    &[2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9],
305    &[7, 8, 0, 7, 0, 6, 6, 0, 2],
306    &[7, 3, 2, 6, 7, 2],
307    &[2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7],
308    &[2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7],
309    &[1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11],
310    &[11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1],
311    &[8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6],
312    &[0, 9, 1, 11, 6, 7],
313    &[7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0],
314    &[7, 11, 6],
315    &[7, 6, 11],
316    &[3, 0, 8, 11, 7, 6],
317    &[0, 1, 9, 11, 7, 6],
318    &[8, 1, 9, 8, 3, 1, 11, 7, 6],
319    &[10, 1, 2, 6, 11, 7],
320    &[1, 2, 10, 3, 0, 8, 6, 11, 7],
321    &[2, 9, 0, 2, 10, 9, 6, 11, 7],
322    &[6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8],
323    &[7, 2, 3, 6, 2, 7],
324    &[7, 0, 8, 7, 6, 0, 6, 2, 0],
325    &[2, 7, 6, 2, 3, 7, 0, 1, 9],
326    &[1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6],
327    &[10, 7, 6, 10, 1, 7, 1, 3, 7],
328    &[10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8],
329    &[0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7],
330    &[7, 6, 10, 7, 10, 8, 8, 10, 9],
331    &[6, 8, 4, 11, 8, 6],
332    &[3, 6, 11, 3, 0, 6, 0, 4, 6],
333    &[8, 6, 11, 8, 4, 6, 9, 0, 1],
334    &[9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6],
335    &[6, 8, 4, 6, 11, 8, 2, 10, 1],
336    &[1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6],
337    &[4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9],
338    &[10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3],
339    &[8, 2, 3, 8, 4, 2, 4, 6, 2],
340    &[0, 4, 2, 4, 6, 2],
341    &[1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8],
342    &[1, 9, 4, 1, 4, 2, 2, 4, 6],
343    &[8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1],
344    &[10, 1, 0, 10, 0, 6, 6, 0, 4],
345    &[4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3],
346    &[10, 9, 4, 6, 10, 4],
347    &[4, 9, 5, 7, 6, 11],
348    &[0, 8, 3, 4, 9, 5, 11, 7, 6],
349    &[5, 0, 1, 5, 4, 0, 7, 6, 11],
350    &[11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5],
351    &[9, 5, 4, 10, 1, 2, 7, 6, 11],
352    &[6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5],
353    &[7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2],
354    &[3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6],
355    &[7, 2, 3, 7, 6, 2, 5, 4, 9],
356    &[9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7],
357    &[3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0],
358    &[6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8],
359    &[9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7],
360    &[1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4],
361    &[4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10],
362    &[7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10],
363    &[6, 9, 5, 6, 11, 9, 11, 8, 9],
364    &[3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5],
365    &[0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11],
366    &[6, 11, 3, 6, 3, 5, 5, 3, 1],
367    &[1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6],
368    &[0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10],
369    &[11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5],
370    &[6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3],
371    &[5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2],
372    &[9, 5, 6, 9, 6, 0, 0, 6, 2],
373    &[1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8],
374    &[1, 5, 6, 2, 1, 6],
375    &[1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6],
376    &[10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0],
377    &[0, 3, 8, 5, 6, 10],
378    &[10, 5, 6],
379    &[11, 5, 10, 7, 5, 11],
380    &[11, 5, 10, 11, 7, 5, 8, 3, 0],
381    &[5, 11, 7, 5, 10, 11, 1, 9, 0],
382    &[10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1],
383    &[11, 1, 2, 11, 7, 1, 7, 5, 1],
384    &[0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11],
385    &[9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7],
386    &[7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2],
387    &[2, 5, 10, 2, 3, 5, 3, 7, 5],
388    &[8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5],
389    &[9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2],
390    &[9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2],
391    &[1, 3, 5, 3, 7, 5],
392    &[0, 8, 7, 0, 7, 1, 1, 7, 5],
393    &[9, 0, 3, 9, 3, 5, 5, 3, 7],
394    &[9, 8, 7, 5, 9, 7],
395    &[5, 8, 4, 5, 10, 8, 10, 11, 8],
396    &[5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0],
397    &[0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5],
398    &[10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4],
399    &[2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8],
400    &[0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11],
401    &[0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5],
402    &[9, 4, 5, 2, 11, 3],
403    &[2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4],
404    &[5, 10, 2, 5, 2, 4, 4, 2, 0],
405    &[3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9],
406    &[5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2],
407    &[8, 4, 5, 8, 5, 3, 3, 5, 1],
408    &[0, 4, 5, 1, 0, 5],
409    &[8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5],
410    &[9, 4, 5],
411    &[4, 11, 7, 4, 9, 11, 9, 10, 11],
412    &[0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11],
413    &[1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11],
414    &[3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4],
415    &[4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2],
416    &[9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3],
417    &[11, 7, 4, 11, 4, 2, 2, 4, 0],
418    &[11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4],
419    &[2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9],
420    &[9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7],
421    &[3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10],
422    &[1, 10, 2, 8, 7, 4],
423    &[4, 9, 1, 4, 1, 7, 7, 1, 3],
424    &[4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1],
425    &[4, 0, 3, 7, 4, 3],
426    &[4, 8, 7],
427    &[9, 10, 8, 10, 11, 8],
428    &[3, 0, 9, 3, 9, 11, 11, 9, 10],
429    &[0, 1, 10, 0, 10, 8, 8, 10, 11],
430    &[3, 1, 10, 11, 3, 10],
431    &[1, 2, 11, 1, 11, 9, 9, 11, 8],
432    &[3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9],
433    &[0, 2, 11, 8, 0, 11],
434    &[3, 2, 11],
435    &[2, 3, 8, 2, 8, 10, 10, 8, 9],
436    &[9, 10, 2, 0, 9, 2],
437    &[2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8],
438    &[1, 10, 2],
439    &[1, 3, 8, 9, 1, 8],
440    &[0, 9, 1],
441    &[0, 3, 8],
442    &[],
443];
444
445/// Edge endpoints: for each edge index (0-11), the two corner indices.
446const EDGE_CORNERS: [[u8; 2]; 12] = [
447    [0, 1],
448    [1, 2],
449    [2, 3],
450    [3, 0], // bottom face
451    [4, 5],
452    [5, 6],
453    [6, 7],
454    [7, 4], // top face
455    [0, 4],
456    [1, 5],
457    [2, 6],
458    [3, 7], // vertical edges
459];
460
461/// Corner offsets within a cell: (dx, dy, dz) for each of 8 corners.
462const CORNER_OFFSETS: [[u32; 3]; 8] = [
463    [0, 0, 0],
464    [1, 0, 0],
465    [1, 1, 0],
466    [0, 1, 0], // bottom
467    [0, 0, 1],
468    [1, 0, 1],
469    [1, 1, 1],
470    [0, 1, 1], // top
471];
472
473// ───────────────────────────────────────────────────────────────────────────
474//  Marching cubes
475// ───────────────────────────────────────────────────────────────────────────
476
477/// Marching cubes isosurface extraction.
478///
479/// Extracts a triangle mesh from a scalar field `grid` sampled on a regular
480/// 3D grid of size `nx * ny * nz` at isolevel `isolevel`.
481///
482/// The grid is indexed as `grid[x + y * nx + z * nx * ny]`.
483/// The cell spacing is `(dx, dy, dz)`.
484/// The origin is at `(origin_x, origin_y, origin_z)`.
485///
486/// `out_vertices` needs `nx * ny * nz * 3` entries (upper bound).
487/// `out_triangles` needs `nx * ny * nz * 5` entries (upper bound, 5 tris per cell).
488///
489/// Returns `(vertex_count, triangle_count)`.
490pub fn marching_cubes(
491    grid: &[f64],
492    nx: usize,
493    ny: usize,
494    nz: usize,
495    dx: f64,
496    dy: f64,
497    dz: f64,
498    origin_x: f64,
499    origin_y: f64,
500    origin_z: f64,
501    isolevel: f64,
502    out_vertices: &mut [Point3],
503    out_triangles: &mut [[u32; 3]],
504) -> Result<(usize, usize), IsosurfaceError> {
505    if nx == 0 || ny == 0 || nz == 0 {
506        return Err(IsosurfaceError::EmptyGrid);
507    }
508    let expected = nx * ny * nz;
509    if grid.len() < expected {
510        return Err(IsosurfaceError::GridSizeMismatch {
511            expected,
512            got: grid.len(),
513        });
514    }
515
516    let max_verts = (nx - 1) * (ny - 1) * (nz - 1) * 30;
517    let max_tris = (nx - 1) * (ny - 1) * (nz - 1) * 10;
518    if out_vertices.len() < max_verts {
519        return Err(IsosurfaceError::BufferTooSmall {
520            needed: max_verts,
521            have: out_vertices.len(),
522        });
523    }
524    if out_triangles.len() < max_tris {
525        return Err(IsosurfaceError::BufferTooSmall {
526            needed: max_tris,
527            have: out_triangles.len(),
528        });
529    }
530
531    let mut vert_count = 0usize;
532    let mut tri_count = 0usize;
533
534    for zk in 0..nz - 1 {
535        for yj in 0..ny - 1 {
536            for xi in 0..nx - 1 {
537                // Sample the 8 corners.
538                let mut corner_vals = [0.0f64; 8];
539                let mut corner_idx = [0usize; 8];
540                for c in 0..8 {
541                    let cx = xi + CORNER_OFFSETS[c][0] as usize;
542                    let cy = yj + CORNER_OFFSETS[c][1] as usize;
543                    let cz = zk + CORNER_OFFSETS[c][2] as usize;
544                    let gi = cx + cy * nx + cz * nx * ny;
545                    corner_vals[c] = grid[gi];
546                    corner_idx[c] = gi;
547                }
548
549                // Compute cube index.
550                let mut cube_idx = 0u8;
551                for c in 0..8 {
552                    if corner_vals[c] < isolevel {
553                        cube_idx |= 1 << c;
554                    }
555                }
556
557                // Skip if entirely inside or outside.
558                let edges = EDGE_TABLE[cube_idx as usize];
559                if edges == 0 {
560                    continue;
561                }
562
563                // Compute edge intersections.
564                let mut edge_verts = [Point3::default(); 12];
565                for e in 0..12 {
566                    if (edges >> e) & 1 == 0 {
567                        continue;
568                    }
569                    let c0 = EDGE_CORNERS[e][0] as usize;
570                    let c1 = EDGE_CORNERS[e][1] as usize;
571                    let v0 = corner_vals[c0];
572                    let v1 = corner_vals[c1];
573
574                    // Linear interpolation factor.
575                    let t = if (v1 - v0).abs() < 1e-20 {
576                        0.5
577                    } else {
578                        (isolevel - v0) / (v1 - v0)
579                    };
580
581                    let p0 = CORNER_OFFSETS[c0];
582                    let p1_off = CORNER_OFFSETS[c1];
583                    let dx0 = (p1_off[0] as f64) - (p0[0] as f64);
584                    let dy0 = (p1_off[1] as f64) - (p0[1] as f64);
585                    let dz0 = (p1_off[2] as f64) - (p0[2] as f64);
586                    let x = origin_x + (xi as f64 + p0[0] as f64 + t * dx0) * dx;
587                    let y = origin_y + (yj as f64 + p0[1] as f64 + t * dy0) * dy;
588                    let z = origin_z + (zk as f64 + p0[2] as f64 + t * dz0) * dz;
589                    edge_verts[e] = Point3::new(x, y, z);
590                }
591
592                // Generate triangles: use edge_verts directly as vertex indices.
593                // Each cell has at most 12 edge-vertices. We emit them once per cell
594                // and reference them by edge index in the triangles.
595                let tri_row = &TRI_TABLE[cube_idx as usize];
596                let mut ti = 0;
597                while tri_row[ti] >= 0 && tri_row[ti + 1] >= 0 && tri_row[ti + 2] >= 0 {
598                    let e0 = tri_row[ti as usize] as usize;
599                    let e1 = tri_row[(ti + 1) as usize] as usize;
600                    let e2 = tri_row[(ti + 2) as usize] as usize;
601
602                    // Emit 3 vertices per triangle (no dedup across triangles).
603                    let v0 = vert_count as u32;
604                    out_vertices[vert_count] = edge_verts[e0];
605                    vert_count += 1;
606                    let v1 = vert_count as u32;
607                    out_vertices[vert_count] = edge_verts[e1];
608                    vert_count += 1;
609                    let v2 = vert_count as u32;
610                    out_vertices[vert_count] = edge_verts[e2];
611                    vert_count += 1;
612
613                    out_triangles[tri_count] = [v0, v1, v2];
614                    tri_count += 1;
615                    ti += 3;
616                }
617            }
618        }
619    }
620
621    Ok((vert_count, tri_count))
622}
623
624// ───────────────────────────────────────────────────────────────────────────
625//  Determinism hash
626// ───────────────────────────────────────────────────────────────────────────
627
628/// FNV-1a hash over vertices and triangles for determinism verification.
629pub fn isosurface_hash(vertices: &[Point3], triangles: &[[u32; 3]]) -> u64 {
630    let mut hash: u64 = 0xcbf29ce484222325;
631    for v in vertices {
632        hash ^= v.x.to_bits();
633        hash = hash.wrapping_mul(0x100000001b3);
634        hash ^= v.y.to_bits();
635        hash = hash.wrapping_mul(0x100000001b3);
636        hash ^= v.z.to_bits();
637        hash = hash.wrapping_mul(0x100000001b3);
638    }
639    for t in triangles {
640        for &idx in t {
641            hash ^= idx as u64;
642            hash = hash.wrapping_mul(0x100000001b3);
643        }
644    }
645    hash
646}
647
648// ───────────────────────────────────────────────────────────────────────────
649//  Tests
650// ───────────────────────────────────────────────────────────────────────────
651
652#[cfg(test)]
653mod tests {
654    use super::*;
655
656    fn sphere_field(
657        nx: usize,
658        ny: usize,
659        nz: usize,
660        cx: f64,
661        cy: f64,
662        cz: f64,
663        r: f64,
664    ) -> Vec<f64> {
665        let mut grid = vec![0.0f64; nx * ny * nz];
666        for k in 0..nz {
667            for j in 0..ny {
668                for i in 0..nx {
669                    let x = i as f64;
670                    let y = j as f64;
671                    let z = k as f64;
672                    let d = ((x - cx).powi(2) + (y - cy).powi(2) + (z - cz).powi(2)).sqrt();
673                    grid[i + j * nx + k * nx * ny] = d - r;
674                }
675            }
676        }
677        grid
678    }
679
680    #[test]
681    fn marching_cubes_sphere() {
682        let nx = 10;
683        let ny = 10;
684        let nz = 10;
685        let grid = sphere_field(nx, ny, nz, 4.5, 4.5, 4.5, 3.0);
686        let max_verts = (nx - 1) * (ny - 1) * (nz - 1) * 30;
687        let max_tris = (nx - 1) * (ny - 1) * (nz - 1) * 10;
688        let mut verts = vec![Point3::default(); max_verts];
689        let mut tris = vec![[0u32; 3]; max_tris];
690
691        let (vc, tc) = marching_cubes(
692            &grid, nx, ny, nz, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, &mut verts, &mut tris,
693        )
694        .unwrap();
695
696        assert!(vc > 0, "sphere should produce vertices");
697        assert!(tc > 0, "sphere should produce triangles");
698    }
699
700    #[test]
701    fn marching_cubes_empty_field() {
702        // All values above isolevel → no surface.
703        let grid = vec![1.0f64; 4 * 4 * 4];
704        let mut verts = vec![Point3::default(); 3 * 3 * 3 * 30];
705        let mut tris = vec![[0u32; 3]; 3 * 3 * 3 * 10];
706
707        let (vc, tc) = marching_cubes(
708            &grid, 4, 4, 4, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, &mut verts, &mut tris,
709        )
710        .unwrap();
711
712        assert_eq!(vc, 0, "uniform field above isolevel → no vertices");
713        assert_eq!(tc, 0, "uniform field above isolevel → no triangles");
714    }
715
716    #[test]
717    fn marching_cubes_determinism() {
718        let nx = 8;
719        let ny = 8;
720        let nz = 8;
721        let grid = sphere_field(nx, ny, nz, 3.5, 3.5, 3.5, 2.5);
722        let max_verts = (nx - 1) * (ny - 1) * (nz - 1) * 30;
723        let max_tris = (nx - 1) * (ny - 1) * (nz - 1) * 10;
724
725        let mut v1 = vec![Point3::default(); max_verts];
726        let mut t1 = vec![[0u32; 3]; max_tris];
727        let mut v2 = vec![Point3::default(); max_verts];
728        let mut t2 = vec![[0u32; 3]; max_tris];
729
730        let (vc1, tc1) = marching_cubes(
731            &grid, nx, ny, nz, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, &mut v1, &mut t1,
732        )
733        .unwrap();
734        let (vc2, tc2) = marching_cubes(
735            &grid, nx, ny, nz, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, &mut v2, &mut t2,
736        )
737        .unwrap();
738
739        assert_eq!(vc1, vc2);
740        assert_eq!(tc1, tc2);
741        assert_eq!(
742            isosurface_hash(&v1[..vc1], &t1[..tc1]),
743            isosurface_hash(&v2[..vc2], &t2[..tc2])
744        );
745    }
746
747    #[test]
748    fn marching_cubes_empty_grid_errors() {
749        let grid: Vec<f64> = vec![];
750        let mut verts = vec![Point3::default(); 1];
751        let mut tris = vec![[0u32; 3]; 1];
752        assert!(matches!(
753            marching_cubes(
754                &grid, 0, 0, 0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, &mut verts, &mut tris
755            ),
756            Err(IsosurfaceError::EmptyGrid)
757        ));
758    }
759
760    #[test]
761    fn marching_cubes_grid_size_mismatch() {
762        let grid = vec![0.0f64; 4]; // too small for 4x4x4
763        let mut verts = vec![Point3::default(); 100];
764        let mut tris = vec![[0u32; 3]; 100];
765        assert!(matches!(
766            marching_cubes(
767                &grid, 4, 4, 4, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, &mut verts, &mut tris
768            ),
769            Err(IsosurfaceError::GridSizeMismatch { .. })
770        ));
771    }
772
773    #[test]
774    fn marching_cubes_plane() {
775        // A flat plane at z=2: field = z - 2.
776        let nx = 5;
777        let ny = 5;
778        let nz = 5;
779        let mut grid = vec![0.0f64; nx * ny * nz];
780        for k in 0..nz {
781            for j in 0..ny {
782                for i in 0..nx {
783                    grid[i + j * nx + k * nx * ny] = k as f64 - 2.0;
784                }
785            }
786        }
787        let max_verts = (nx - 1) * (ny - 1) * (nz - 1) * 30;
788        let max_tris = (nx - 1) * (ny - 1) * (nz - 1) * 10;
789        let mut verts = vec![Point3::default(); max_verts];
790        let mut tris = vec![[0u32; 3]; max_tris];
791
792        let (vc, tc) = marching_cubes(
793            &grid, nx, ny, nz, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, &mut verts, &mut tris,
794        )
795        .unwrap();
796
797        // A plane should produce triangles.
798        assert!(tc > 0, "plane should produce triangles");
799        assert!(vc > 0, "plane should produce vertices");
800    }
801
802    #[test]
803    fn tri_table_edges_match_edge_table() {
804        // Rigorous table-correctness gate (no geometry needed): for every one of
805        // the 256 configurations, the SET of edges referenced by the triangle
806        // table must equal the independently-computed set of edges the surface
807        // crosses (EDGE_TABLE, derived from the corner-sign logic). This catches
808        // any wrong or missing edge index in the 256-row constant.
809        for cfg in 0..256usize {
810            let mut used: u16 = 0;
811            for &e in &TRI_TABLE[cfg] {
812                if e < 0 {
813                    break;
814                }
815                assert!(
816                    (0..12).contains(&(e as i32)),
817                    "cfg {cfg}: edge index {e} out of range 0..12"
818                );
819                used |= 1 << (e as u16);
820            }
821            assert_eq!(
822                used, EDGE_TABLE[cfg],
823                "cfg {cfg}: triangulated edge set {used:#014b} != crossing set {:#014b}",
824                EDGE_TABLE[cfg]
825            );
826        }
827    }
828
829    #[test]
830    fn tri_table_rows_are_whole_triangles() {
831        // Every row is a run of complete triangles (a multiple of 3 edge refs)
832        // followed only by -1 padding.
833        for cfg in 0..256usize {
834            let row = &TRI_TABLE[cfg];
835            let mut n = 0usize;
836            while n < 16 && row[n] >= 0 {
837                n += 1;
838            }
839            assert_eq!(
840                n % 3,
841                0,
842                "cfg {cfg}: {n} edge refs is not a whole number of triangles"
843            );
844            for &e in &row[n..] {
845                assert_eq!(e, -1, "cfg {cfg}: non-(-1) padding after terminator");
846            }
847        }
848    }
849
850    #[test]
851    fn sphere_isosurface_is_manifold() {
852        // Extract a sphere (a closed level set fully interior to the grid) and
853        // verify the triangulation is a valid 2-manifold: after merging
854        // coincident vertices by position, no triangle is degenerate and no edge
855        // is shared by more than two triangles. The correct MC table produces
856        // manifold geometry; the previous fan triangulation did not on the
857        // ambiguous cube configurations. (Center on half-integers + r=3.3 keeps
858        // the surface off the grid corners, so no edge-vertex lands on a corner.)
859        use std::collections::HashMap;
860        let (nx, ny, nz) = (12usize, 12usize, 12usize);
861        let grid = sphere_field(nx, ny, nz, 5.5, 5.5, 5.5, 3.3);
862        let max_verts = (nx - 1) * (ny - 1) * (nz - 1) * 30;
863        let max_tris = (nx - 1) * (ny - 1) * (nz - 1) * 10;
864        let mut verts = vec![Point3::default(); max_verts];
865        let mut tris = vec![[0u32; 3]; max_tris];
866        let (vc, tc) = marching_cubes(
867            &grid, nx, ny, nz, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, &mut verts, &mut tris,
868        )
869        .unwrap();
870        assert!(tc > 0, "sphere must produce triangles");
871
872        // Merge coincident vertices by quantized position (edge-vertices on a
873        // shared face are computed from identical corner data, so they are
874        // bit-identical and collapse to one merged index).
875        let key = |p: Point3| -> (i64, i64, i64) {
876            (
877                (p.x * 1e6).round() as i64,
878                (p.y * 1e6).round() as i64,
879                (p.z * 1e6).round() as i64,
880            )
881        };
882        let mut merged: HashMap<(i64, i64, i64), u32> = HashMap::new();
883        let mut remap = vec![0u32; vc];
884        for i in 0..vc {
885            let k = key(verts[i]);
886            let next = merged.len() as u32;
887            remap[i] = *merged.entry(k).or_insert(next);
888        }
889
890        // Count undirected edges over merged indices.
891        let mut edge_count: HashMap<(u32, u32), u32> = HashMap::new();
892        for t in &tris[..tc] {
893            let (a, b, c) = (
894                remap[t[0] as usize],
895                remap[t[1] as usize],
896                remap[t[2] as usize],
897            );
898            assert!(
899                a != b && b != c && a != c,
900                "degenerate triangle after merge: {a},{b},{c}"
901            );
902            for (u, v) in [(a, b), (b, c), (c, a)] {
903                let e = if u < v { (u, v) } else { (v, u) };
904                *edge_count.entry(e).or_insert(0) += 1;
905            }
906        }
907        // 2-manifold: no edge shared by more than two triangles.
908        for (e, &cnt) in &edge_count {
909            assert!(
910                cnt <= 2,
911                "non-manifold edge {e:?} shared by {cnt} triangles"
912            );
913        }
914    }
915}
qualia_core_db/specialized_libs/computational_geometry/isosurface.rs

qualia_core_db/specialized_libs/computational_geometry/
isosurface.rs
