Marching Cubes¶
Marching cubes is the workhorse of iso-surface extraction: it visits every grid cell and cuts triangles out of it based on which corners are inside the surface. It is the default choice for most fields.
import torch
import isoext
from isoext import viewer
# Setup: create a grid with a sphere
grid = isoext.UniformGrid([64, 64, 64])
points = grid.get_points()
grid.set_values(points.norm(dim=-1) - 0.7)
Basic Usage¶
vertices, faces = isoext.marching_cubes(grid, level=0.0, method="vega")
grid: AUniformGridorSparseGridwith values setlevel: The iso-value to extract (default: 0.0)method: Algorithm variant ("vega","lewiner","nagae", or"lorensen")
vertices, faces = isoext.marching_cubes(grid)
print(f"Vertices: {vertices.shape}") # (N, 3) float32
print(f"Faces: {faces.shape}") # (M, 3) uint32
viewer.embed(vertices, faces)
Vertices: torch.Size([9168, 3])
Faces: torch.Size([18332, 3])
How It Works¶
Marching cubes visits every cell independently. Each of a cell’s 8 corners is either inside or outside the surface, giving 256 possible sign patterns, and a lookup table maps each pattern to a set of triangles. The triangle vertices sit on the cell’s edges. If one end of an edge is inside the surface and the other is outside, the surface has to cross that edge somewhere in between, and interpolating the two values estimates where: with level-relative endpoint values \(v_a\) and \(v_b\), the crossing sits at the fraction \(t = v_a / (v_a - v_b)\) of the way along the edge. An endpoint whose value is close to zero is close to the surface, so it pulls the crossing toward itself.
The simplest pattern has one corner inside, the red dot below. The surface crosses the three edges that connect it to outside corners; the gold dots mark those crossings. For this pattern the lookup table produces one triangle, and the three crossings are its vertices:
cell = isoext.UniformGrid([2, 2, 2])
cell.set_values(torch.tensor(
[-0.5, 0.3, 0.4, 0.5, 0.6, 0.4, 0.3, 0.2], device="cuda"
).reshape(2, 2, 2))
crossings = isoext.get_intersection(cell).get_points()
v, f = isoext.marching_cubes(cell)
print(f"{f.shape[0]} triangle(s)")
viewer.embed(v, f, color="steelblue", flat_shading=True, side="double", height=300,
grid=cell, draw=lambda s: viewer.add_points(s, crossings, point_size=0.09))
1 triangle(s)
Algorithm Variants¶
The method argument selects between four variants that share the
same interface:
vega(default) resolves topological ambiguities by following the field’s interpolant, using the corrected interior test.lewinerfollows the interpolant with the older interior test; kept for comparison with other implementations.nagaeuses fixed tables; always closed and slightly faster.lorensenis the original 1987 algorithm; ambiguous cells can leave small cracks.
See Marching Cubes Variants for the differences in detail.
for method in ["nagae", "lorensen", "lewiner", "vega"]:
v, f = isoext.marching_cubes(grid, method=method)
print(f"{method:9s} {f.shape[0]:,} triangles")
nagae 18,332 triangles
lorensen 18,332 triangles
lewiner 18,332 triangles
vega 18,332 triangles
Iso-Level¶
The level parameter controls which iso-surface to extract. For signed distance fields, level=0 gives the surface. Other values give offset surfaces:
# Extract at different iso-levels
for level in [-0.2, 0.0, 0.2]:
v, f = isoext.marching_cubes(grid, level=level)
print(f"level={level:+.1f}: {v.shape[0]:,} vertices (radius ≈ {0.7 - level:.1f})")
level=-0.2: 4,728 vertices (radius ≈ 0.9)
level=+0.0: 9,168 vertices (radius ≈ 0.7)
level=+0.2: 15,072 vertices (radius ≈ 0.5)
Saving Meshes¶
Use write_obj to save the mesh to an OBJ file:
vertices, faces = isoext.marching_cubes(grid)
isoext.write_obj("mesh.obj", vertices, faces)
print("Saved to mesh.obj")
Saved to mesh.obj