isomesh 0.0.12

Engine-agnostic isosurface extraction: signed distance field in, triangles out. no_std, f32 and f64, no math library in the public API.
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
//! A-014f's tests.
//!
//! The load-bearing ones are `gyroids_flipped_edges_go_to_zero` — the acceptance
//! — and `the_two_sheets_of_a_thin_plate_are_not_merged`, which is the case
//! A-014e's per-triangle vote was written to protect and the one P-7 was
//! registered to de-risk.

use alloc::vec;
use alloc::vec::Vec;

use super::{OrientReport, orient};
use crate::mesh::MeshBuffer;

/// A quad as two triangles sharing edge 1-2, coherently wound and facing the
/// way its normals say is out.
///
/// Coherent means the two traverse their shared edge in *opposite* directions:
/// triangle `0 1 2` goes 1 → 2, so its neighbour must go 2 → 1, which is
/// `1 3 2` and not `1 2 3`.
fn quad() -> MeshBuffer<f64> {
    let mut mesh = MeshBuffer::<f64>::new();
    mesh.positions = vec![
        [0.0, 0.0, 0.0],
        [1.0, 0.0, 0.0],
        [0.0, 1.0, 0.0],
        [1.0, 1.0, 0.0],
    ];
    mesh.normals = vec![[0.0, 0.0, 1.0]; 4];
    mesh.indices = vec![0, 1, 2, 1, 3, 2];
    mesh
}

/// The same quad with its second triangle wound backwards.
fn quad_with_one_reversed() -> MeshBuffer<f64> {
    let mut mesh = quad();
    mesh.indices[4..6].swap(0, 1);
    mesh
}

/// How many interior edges are traversed the *same* way by both their faces —
/// the defect this whole pass exists to remove, counted independently of
/// `validate` so the two are not proving each other.
fn same_direction_edges(mesh: &MeshBuffer<f64>) -> usize {
    use alloc::collections::BTreeMap;
    let mut seen: BTreeMap<(u32, u32), i32> = BTreeMap::new();
    for tri in mesh.indices.as_chunks::<3>().0 {
        for k in 0..3 {
            let (a, b) = (tri[k], tri[(k + 1) % 3]);
            let key = if a <= b { (a, b) } else { (b, a) };
            *seen.entry(key).or_default() += if a <= b { 1 } else { -1 };
        }
    }
    // +-2 means both faces went the same way round; 0 means they opposed.
    seen.values().filter(|v| v.abs() == 2).count()
}

#[test]
fn an_empty_mesh_is_accepted_and_does_nothing() {
    let mut mesh = MeshBuffer::<f64>::new();
    assert_eq!(orient(&mut mesh).expect("empty"), OrientReport::default());
}

#[test]
fn an_index_past_the_end_is_refused_before_anything_moves() {
    let mut mesh = quad_with_one_reversed();
    mesh.indices[4] = 99;
    let before = mesh.indices.clone();
    assert!(orient(&mut mesh).is_err());
    assert_eq!(mesh.indices, before, "a refused call must not have edited");
}

#[test]
fn a_reversed_neighbour_is_flipped_back() {
    let mut mesh = quad_with_one_reversed();
    assert_eq!(same_direction_edges(&mesh), 1);

    let report = orient(&mut mesh).expect("orientable");
    assert_eq!(report.components, 1);
    assert_eq!(report.triangles_flipped, 1);
    assert!(report.is_orientable());
    assert_eq!(same_direction_edges(&mesh), 0);
}

#[test]
fn an_already_coherent_mesh_is_left_alone() {
    let mut mesh = quad();
    let before = mesh.indices.clone();
    let report = orient(&mut mesh).expect("orientable");
    assert!(report.is_noop());
    assert_eq!(mesh.indices, before);
}

/// Winding is not merely *consistent*, it is consistent **the right way round**.
/// A pass that made everything agree with each other and disagreed with the
/// field would satisfy every edge check and render inside out.
#[test]
fn a_wholly_inverted_mesh_is_turned_back_outward() {
    let mut mesh = quad();
    // Reverse both triangles: internally coherent, globally inside out.
    for tri in mesh.indices.as_chunks_mut::<3>().0 {
        tri.swap(1, 2);
    }
    assert_eq!(same_direction_edges(&mesh), 0, "already self-consistent");

    let report = orient(&mut mesh).expect("orientable");
    assert_eq!(report.triangles_flipped, 2, "both, not neither");
    assert_eq!(mesh.indices, quad().indices);
}

#[test]
fn two_disjoint_patches_are_two_components() {
    let mut mesh = quad();
    let offset = mesh.positions.len() as u32;
    let second = quad_with_one_reversed();
    mesh.positions
        .extend(second.positions.iter().map(|p| [p[0] + 10.0, p[1], p[2]]));
    mesh.normals.extend(second.normals.iter().copied());
    mesh.indices
        .extend(second.indices.iter().map(|i| i + offset));

    let report = orient(&mut mesh).expect("orientable");
    assert_eq!(report.components, 2);
    assert_eq!(
        report.triangles_flipped, 1,
        "only the second patch was wrong"
    );
    assert_eq!(same_direction_edges(&mesh), 0);
}

/// The seed must be the most *confident* triangle, not the first one.
///
/// Here triangle 0 has a normal lying in its own plane — no information — while
/// triangle 1's is decisive and says the mesh is inside out. Seeding by index
/// would take triangle 0's non-vote, leave the pair as it found them, and report
/// success on a mesh that is still wrong.
#[test]
fn the_seed_is_the_most_confident_triangle_not_the_first() {
    let mut mesh = quad();
    for tri in mesh.indices.as_chunks_mut::<3>().0 {
        tri.swap(1, 2);
    }
    // Vertex 0 is used only by triangle 0, so this blunts that triangle's vote
    // without touching triangle 1's.
    mesh.normals[0] = [0.0, 0.0, -1.0];

    let report = orient(&mut mesh).expect("orientable");
    assert_eq!(
        report.triangles_flipped, 2,
        "the indecisive triangle was allowed to seed the component"
    );
    assert_eq!(mesh.indices, quad().indices);
}

/// A triangle with no usable vote at all still gets a *defined* answer, from its
/// neighbour rather than from the order the faces happened to be visited in.
/// A-014f's second acceptance clause.
#[test]
fn a_triangle_with_no_vote_of_its_own_inherits_one() {
    let mut mesh = quad_with_one_reversed();
    // Zero normals on the vertices unique to the second triangle: its own vote
    // is now the zero vector, which is neither outward nor inward.
    mesh.normals[3] = [0.0, 0.0, 0.0];

    let report = orient(&mut mesh).expect("orientable");
    assert_eq!(same_direction_edges(&mesh), 0, "it inherited nothing");
    assert!(report.is_orientable());
}

/// Determinism, and specifically that it does not depend on the order the faces
/// arrive in — the failure mode a `HashMap`-backed adjacency would introduce.
#[test]
fn the_result_does_not_depend_on_the_order_the_faces_are_listed_in() {
    let base = {
        let mut m = quad_with_one_reversed();
        orient(&mut m).expect("orientable");
        m.indices
    };
    // The same two triangles, listed the other way round.
    let mut swapped = quad_with_one_reversed();
    let (a, b) = swapped.indices.split_at(3);
    let reordered: Vec<u32> = b.iter().chain(a.iter()).copied().collect();
    swapped.indices = reordered;
    orient(&mut swapped).expect("orientable");
    assert_eq!(same_direction_edges(&swapped), 0);
    assert_eq!(base.len(), swapped.indices.len());
}

/// **A Möbius band is reported, not silently half-fixed.**
///
/// Five quads in a ring with a half-twist: propagation comes back to where it
/// started with the opposite sign, and no assignment of windings can make every
/// edge oppose. The honest answer is to say so.
#[test]
fn a_mobius_band_is_reported_as_non_orientable() {
    let strips = 5;
    let mut mesh = MeshBuffer::<f64>::new();
    for i in 0..strips {
        let t = i as f64 / strips as f64;
        mesh.positions.push([t, 0.0, 0.0]);
        mesh.positions.push([t, 1.0, 0.0]);
        mesh.normals.push([0.0, 0.0, 1.0]);
        mesh.normals.push([0.0, 0.0, 1.0]);
    }
    for i in 0..strips {
        let a = (i * 2) as u32;
        let b = a + 1;
        let next = ((i + 1) % strips) as u32 * 2;
        // The last quad joins back with its two vertices swapped: the twist.
        let (c, d) = if i + 1 == strips {
            (next + 1, next)
        } else {
            (next, next + 1)
        };
        mesh.indices.extend_from_slice(&[a, b, c]);
        mesh.indices.extend_from_slice(&[b, d, c]);
    }

    let report = orient(&mut mesh).expect("indices are in range");
    assert_eq!(report.components, 1);
    assert!(
        !report.is_orientable(),
        "a half-twisted ring was reported orientable"
    );
    assert_eq!(report.non_orientable_components, 1);
}

/// **A-014f's acceptance.** `gyroid`'s inconsistently-oriented edges go to zero,
/// and `thin_plate` — the field A-014e's per-triangle vote exists to protect —
/// does not regress. Measured on the real extractor at all three resolutions,
/// because M-182 found `thin_plate` failing at the two the suite never tested.
#[test]
fn every_reference_field_comes_out_coherently_oriented() {
    use crate::shape::RuntimeShape3;
    use crate::subgrid::extract::SubgridMarchingTetrahedra;
    use crate::validate::{ValidateConfig, validate_indexed};

    let mut rows = 0;
    let mut got: Vec<(&str, u32, u64, u64, u64)> = Vec::new();
    for n in [17u32, 25, 33] {
        crate::for_each_reference_field!(f64, |name, field| {
            let (lo, hi) = crate::fields::ReferenceField::domain(&field);
            let cell = (hi[0] - lo[0]) / f64::from(n - 1);
            let shape = RuntimeShape3::new([n; 3]).expect("a cubic grid");
            let cfg = ValidateConfig::from_cell_size(cell).expect("a valid spacing");

            let mut mesh = MeshBuffer::<f64>::default();
            SubgridMarchingTetrahedra::<f64>::new(16)
                .expect("valid")
                .extract(&field, &shape, lo, cell, &mut mesh)
                .unwrap_or_else(|e| panic!("{name} {n}: {e}"));
            crate::weld::Welder::<f64>::new()
                .weld(&mut mesh, crate::weld::epsilon_for(cell))
                .expect("weld");

            let before = validate_indexed(&mesh.positions, &mesh.indices, &cfg);
            let report = orient(&mut mesh).expect("indices in range");
            let after = validate_indexed(&mesh.positions, &mesh.indices, &cfg);

            std::println!(
                "{name:<15} {n}³  flipped edges {:>4} -> {:>3}  \
                 (components {}, triangles flipped {}, non-orientable {})",
                before.inconsistently_oriented_edges,
                after.inconsistently_oriented_edges,
                report.components,
                report.triangles_flipped,
                report.non_orientable_components,
            );

            got.push((
                name,
                n,
                before.inconsistently_oriented_edges,
                after.inconsistently_oriented_edges,
                before.non_manifold_edges,
            ));
            // Orientation must not disturb anything else about the mesh.
            assert_eq!(after.vertices, before.vertices, "{name} {n}³ vertices");
            assert_eq!(after.faces, before.faces, "{name} {n}³ faces");
            assert_eq!(
                (after.non_manifold_edges, after.non_manifold_vertices),
                (before.non_manifold_edges, before.non_manifold_vertices),
                "{name} {n}³ manifoldness moved"
            );
            assert_eq!(
                after.boundary_edges, before.boundary_edges,
                "{name} {n}³ boundary moved"
            );
            rows += 1;
        });
    }
    assert_eq!(rows, 42, "three resolutions over fourteen fields");

    // **The law is now unconditional, and A-019 is what made it so (M-227).**
    // M-187 stated it with a caveat: orientation drives the count to zero on
    // every row whose edges are all manifold, and *"where it does not, the mesh
    // has edges carrying more than two faces, and no assignment of windings can
    // make four faces pairwise oppose across one edge"*. The second half was a
    // property of the **walk**, not of the mesh. A walk that crosses a four-face
    // edge must indeed leave a residue there — but it need not cross it, and once
    // propagation stops instead, every row reaches **zero**: `csg_difference`
    // 6 → 0, `fbm_terrain` 19 → 0, `noise_cavity` 1,629 → 0 and 1,580 → 0 where
    // it used to *grow* to 2,422.
    for &(name, n, before, after, non_manifold) in &got {
        if non_manifold == 0 {
            assert_eq!(
                after, 0,
                "{name} {n}³: {after} flipped edges left on a mesh with no \
                 non-manifold edge, from {before}"
            );
        } else {
            // **`after <= before` is a law again, and A-019 is what restored
            // it (M-227).** It was falsified at A-002e — `noise_cavity` went
            // 1,580 → 2,422 at 25³ — and the cause was propagation *crossing* a
            // four-face edge: with no well-defined neighbour there, the walk took
            // whichever it reached first and carried a consistent winding across
            // a patch that was consistent with the wrong side. Propagation now
            // stops at such an edge, so each sheet is seeded and oriented on its
            // own evidence.
            //
            // **The residue may now reach zero on a non-manifold mesh**, which
            // M-187's law said it could not: that law was about a walk that
            // crossed the edge, where four faces cannot pairwise oppose. A walk
            // that does not cross it has no such obligation, and
            // `csg_difference` at 17³ goes 6 → 0.
            assert!(
                after <= before,
                "{name} {n}³: {before} -> {after} with {non_manifold} \
                 non-manifold edges -- orientation made it worse"
            );
        }
    }

    // The full census, pinned. The law above explains its shape; this is what
    // stops any row of it moving unnoticed.
    let table: Vec<(&str, u32, u64, u64)> =
        got.iter().map(|&(f, n, b, a, _)| (f, n, b, a)).collect();
    assert_eq!(
        table,
        alloc::vec![
            ("sphere", 17, 0, 0),
            ("torus", 17, 0, 0),
            ("box_exact", 17, 0, 0),
            ("csg_difference", 17, 6, 0),
            ("thin_plate", 17, 0, 0),
            ("gyroid", 17, 138, 0),
            ("fbm_terrain", 17, 19, 0),
            ("noise_cavity", 17, 1629, 0),
            // R-177's six prescribed-genus fields, pinned at landing.
            ("ball_drilled_g1", 17, 0, 0),
            ("graph_theta_g2", 17, 17, 0),
            ("ball_drilled_g2", 17, 27, 0),
            ("ball_drilled_g3", 17, 25, 0),
            ("graph_k4_g3", 17, 42, 0),
            ("graph_cube_g5", 17, 0, 0),
            ("sphere", 25, 0, 0),
            ("torus", 25, 6, 0),
            ("box_exact", 25, 0, 0),
            ("csg_difference", 25, 0, 0),
            ("thin_plate", 25, 8, 0),
            ("gyroid", 25, 150, 0),
            ("fbm_terrain", 25, 29, 0),
            ("noise_cavity", 25, 1580, 0),
            ("ball_drilled_g1", 25, 0, 0),
            ("graph_theta_g2", 25, 6, 0),
            ("ball_drilled_g2", 25, 35, 0),
            ("ball_drilled_g3", 25, 44, 0),
            ("graph_k4_g3", 25, 0, 0),
            ("graph_cube_g5", 25, 0, 0),
            ("sphere", 33, 0, 0),
            ("torus", 33, 0, 0),
            ("box_exact", 33, 0, 0),
            ("csg_difference", 33, 36, 0),
            ("thin_plate", 33, 6, 0),
            ("gyroid", 33, 330, 0),
            ("fbm_terrain", 33, 53, 0),
            ("noise_cavity", 33, 1477, 0),
            ("ball_drilled_g1", 33, 56, 0),
            ("graph_theta_g2", 33, 6, 0),
            ("ball_drilled_g2", 33, 60, 0),
            ("ball_drilled_g3", 33, 92, 0),
            ("graph_k4_g3", 33, 0, 0),
            ("graph_cube_g5", 33, 0, 0),
        ],
        "the flipped-edge census moved"
    );

    // A-014f's headline acceptance, named rather than left to the loop.
    let row = |f: &str, n: u32| {
        got.iter()
            .find(|g| g.0 == f && g.1 == n)
            .map(|g| (g.2, g.3))
            .expect("in table")
    };
    assert_eq!(row("gyroid", 17), (138, 0));
    assert_eq!(row("gyroid", 25), (150, 0));
    assert_eq!(row("gyroid", 33), (330, 0));
    // And `thin_plate` -- the field the per-triangle vote exists to protect --
    // is not merely un-regressed, it is fixed at the two resolutions M-182
    // found it failing at.
    assert_eq!(row("thin_plate", 17), (0, 0));
    assert_eq!(row("thin_plate", 25), (8, 0));
    assert_eq!(row("thin_plate", 33), (6, 0));
}