Marching Cubes Variants¶

Marching cubes turns each cell’s corner signs into triangles by table lookup. Some sign patterns have more than one valid triangulation, and the four variants differ in how they pick one. The consequences range from cracked meshes to wrong topology, shown below on fields of a few cells. The example meshes are open surfaces, so they are rendered double-sided; the overlays show the grid, with red dots for corners inside the surface and blue dots for corners outside.

import torch

import isoext
from isoext import viewer
from isoext.utils import gaussian_smooth

Watertightness¶

A face of a cell is ambiguous when two diagonally opposite corners are inside and the other two outside: the surface can cut the face along either diagonal. Neighboring cells share faces, so both cells must pick the same diagonal, or the mesh gets a hole. The lorensen tables (Lorensen and Cline [1987]) were built with reflections, which can make the two cells disagree. The nagae tables (Nagae et al. [1993]) avoid reflections and always agree; so do lewiner and vega. Counting broken meshes over 50 random closed fields:

open_meshes = {"lorensen": 0, "nagae": 0, "lewiner": 0, "vega": 0}
for seed in range(50):
    torch.manual_seed(seed)
    values = gaussian_smooth(torch.randn(12, 12, 12, device="cuda"), sigma=1.2)
    values[0], values[-1] = 1.0, 1.0
    values[:, 0], values[:, -1] = 1.0, 1.0
    values[:, :, 0], values[:, :, -1] = 1.0, 1.0
    grid = isoext.UniformGrid([12, 12, 12])
    grid.set_values(values)
    for method in open_meshes:
        v, f = isoext.marching_cubes(grid, method=method)
        edges = torch.cat([f[:, [0, 1]], f[:, [1, 2]], f[:, [2, 0]]]).sort(dim=-1).values
        _, counts = torch.unique(edges, dim=0, return_counts=True)
        open_meshes[method] += int((counts != 2).any())
print("fields with broken meshes:", open_meshes)

fields with broken meshes: {'lorensen': 49, 'nagae': 0, 'lewiner': 0, 'vega': 0}

A Crack, Up Close¶

Here are two cells sharing an ambiguous face. lorensen cuts the face along a different diagonal on each side, so the two surface pieces do not meet. The zigzag lines crossing the middle are the two unmatched borders of the crack, one from each cell:

values = torch.tensor(
    [0.1, 0.5, 0.4, 0.7, -0.6, 0.7, 1.0, -0.4, -1.0, -0.3, -0.2, -0.5],
    device="cuda",
).reshape(3, 2, 2)
pair = isoext.UniformGrid([3, 2, 2], aabb_min=[-1, -0.5, -0.5], aabb_max=[1, 0.5, 0.5])
pair.set_values(values)

v, f = isoext.marching_cubes(pair, method="lorensen")
viewer.embed(v, f, color="steelblue", flat_shading=True, side="double", height=300, grid=pair)

The same field with nagae seals the crack:

v, f = isoext.marching_cubes(pair, method="nagae")
viewer.embed(v, f, color="seagreen", flat_shading=True, side="double", height=300, grid=pair)

Topology of Ambiguous Cells¶

Sealing the cracks still leaves a choice: are the surface pieces in an ambiguous cell connected or separate? Fixed tables like nagae always answer the same way. lewiner and vega instead evaluate the trilinear interpolant – the smooth field you get by blending the eight corner values across the cell – and follow its topology, implementing the MC33 method of Chernyaev [1995]. In the cell below, the fixed table gives three separate pieces, while the interpolant connects two of them:

values = torch.tensor([0.5, -0.5, -0.2, 0.9, 0.2, 0.1, 0.3, -0.2], device="cuda").reshape(2, 2, 2)
cell = isoext.UniformGrid([2, 2, 2])
cell.set_values(values)

v, f = isoext.marching_cubes(cell, method="nagae")
print(f"nagae:   {f.shape[0]} triangles")
viewer.embed(v, f, color="seagreen", flat_shading=True, side="double", height=300, grid=cell)

nagae:   3 triangles

v, f = isoext.marching_cubes(cell, method="vega")
print(f"vega:    {f.shape[0]} triangles")
viewer.embed(v, f, color="steelblue", flat_shading=True, side="double", height=300, grid=cell)

vega:    5 triangles

Lewiner Limitations¶

Deciding whether two pieces connect through the inside of a cell takes an interior test. The lewiner variant ports the implementation of Lewiner et al. [2003], whose interior test checks the interpolant on a single cross-section instead of analyzing it fully. That has two consequences, documented by Custodio et al. [2013]: the test can misjudge tunnels in rare configurations, and mirroring the input field can change the answer. The field below triggers both. lewiner connects the surface into a tunnel, and extracting the mirrored field gives two separate sheets instead. The meshes stay watertight and crack free either way.

values = torch.tensor(
    [0.5625, -0.4375, 0.0625, 0.5625, 0.1875, 0.125, -0.4375, 0.125],
    device="cuda",
).reshape(2, 2, 2)
cell.set_values(values)
v, f = isoext.marching_cubes(cell, method="lewiner")
print(f"lewiner, original: {f.shape[0]} triangles")
viewer.embed(v, f, color="tomato", flat_shading=True, side="double", height=300, grid=cell)

lewiner, original: 6 triangles

cell.set_values(values.permute(2, 1, 0).contiguous())  # mirrored field
v, f = isoext.marching_cubes(cell, method="lewiner")
print(f"lewiner, mirrored: {f.shape[0]} triangles")
viewer.embed(v, f, color="gold", flat_shading=True, side="double", height=300, grid=cell)

lewiner, mirrored: 2 triangles

The vega variant ports the implementation of Vega et al. [2019], which replaces the interior test with one that analyzes the interpolant fully. Here it extracts the two sheets the interpolant actually contains, and mirroring the field does not change the answer:

cell.set_values(values)
v, f = isoext.marching_cubes(cell, method="vega")
print(f"vega, original: {f.shape[0]} triangles")
viewer.embed(v, f, color="steelblue", flat_shading=True, side="double", height=300, grid=cell)

vega, original: 2 triangles

cell.set_values(values.permute(2, 1, 0).contiguous())  # mirrored field
v, f = isoext.marching_cubes(cell, method="vega")
print(f"vega, mirrored: {f.shape[0]} triangles")
viewer.embed(v, f, color="gold", flat_shading=True, side="double", height=300, grid=cell)

vega, mirrored: 2 triangles

Choosing a Variant¶

vega is the default: watertight, faithful to the trilinear interpolant, and about 10% slower than the fixed tables. lewiner follows the interpolant too but has the limitations shown above; it is kept for comparison with other ports of the same implementation, such as scikit-image’s. Use nagae when that margin matters and the topology of ambiguous cells does not. lorensen is included for reference and comparison with other implementations.

Lookup Tables¶

The lorensen and nagae tables are generated by luts/gen_mc_lut.py, which reads base cases from a JSON file and expands them to all 256 configurations via rotations and, optionally, reflections. New table-driven variants can be added the same way.

The lewiner tables are converted from the reference implementation distributed with scikit-image (BSD) by luts/convert_mc33_luts.py, rather than generated: MC33 tables couple triangulations with face and interior test descriptors, and re-deriving that machinery is where historical implementations accumulated bugs.

The vega table is converted from the MC33_c_library of Vega et al. [2019] (MIT) by luts/convert_vega_luts.py. The converter re-derives every offset the case dispatch can reach and decodes the triangle strip at each one, so a conversion error fails at generation time instead of on the GPU.

References¶

[1]

William E. Lorensen and Harvey E. Cline. Marching cubes: a high resolution 3D surface construction algorithm. In Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH '87), 163–169. 1987. doi:10.1145/37401.37422.

[2]

Takanori Nagae, Takeshi Agui, and Hiroshi Nagahashi. Surface construction and contour generation from volume data. In Medical Imaging 1993: Image Processing, Proc. SPIE 1898, 74–84. 1993. doi:10.1117/12.154567.

[3]

Evgeni Chernyaev. Marching cubes 33: construction of topologically correct isosurfaces. Technical Report CERN-CN-95-17, CERN, Geneva, 1995. URL: https://repository.cern/records/7zfxg-q0t96.

[4]

Thomas Lewiner, Hélio Lopes, Antônio Wilson Vieira, and Geovan Tavares. Efficient implementation of marching cubes' cases with topological guarantees. Journal of Graphics Tools, 8(2):1–15, 2003. doi:10.1080/10867651.2003.10487582.

[5]

Lis Custodio, Tiago Etiene, Sinesio Pesco, and Claudio Silva. Practical considerations on marching cubes 33 topological correctness. Computers & Graphics, 37(7):840–850, 2013. doi:10.1016/j.cag.2013.04.004.

[6] (1,2)

David Vega, Javier Abache, and David Coll. A fast and memory-saving marching cubes 33 implementation with the correct interior test. Journal of Computer Graphics Techniques (JCGT), 8(3):1–18, 2019. URL: https://jcgt.org/published/0008/03/01/.