Marching Cubes

Marching cubes converts a scalar field sampled on a regular 3D grid into a triangle mesh of one of its level sets. It visits each cube of eight neighbouring samples, decides from the signs at the corners which cube edges the surface crosses, looks up how to join those crossings into triangles, and places each vertex on its edge by linear interpolation. Lorensen and Cline (1987) introduced it for CT, MR and SPECT slices. In 3D vision it is the usual last step that turns a fused signed distance function, a Poisson indicator function or a neural field into a polygon mesh.

Problem

Given samples of $f : \mathbb{R}^3 \to \mathbb{R}$ on a grid with spacing $h$, find a triangle mesh approximating the level set $\{\mathbf{x} : f(\mathbf{x}) = \alpha\}$, the zero set for a signed distance. Between samples $f$ is unknown; marching cubes assumes it is linear along each grid edge.

Inputs and Outputs

The input is a 3D array, its spacing and $\alpha$; the output is one vertex per grid edge whose ends lie on opposite sides of $\alpha$, triangles indexing them, and optionally vertex normals. Lorensen and Cline count a corner as inside when its value is at least $\alpha$; a signed distance is negative inside, which flips only the direction the triangles face.

Algorithm

For each cube (Lorensen and Cline, 1987, §4):

1. Classify the eight corners against $\alpha$ and pack the results into an 8-bit index, giving $2^8 = 256$ cases.
2. Look up the case in a precomputed table listing the intersected edges and their triangles. Cases 0 and 255, all outside or all inside, produce nothing.
3. Interpolate a vertex on each intersected edge.
4. Compute normals at the corners by central differences and interpolate them to the vertices.
5. Output the triangles and march to the next cube.

Neighbouring cubes share edges: an interior cube needs only three new edge intersections, and reusing those of the previous pixel and line sped up Lorensen and Cline’s implementation threefold (§5.1). Keying each vertex by its grid edge, as the example does, lets the four cubes around an edge share one vertex.

From 256 cases to 15 patterns

Lorensen and Cline triangulated the cases by hand using two symmetries: swapping inside and outside leaves the topology of the surface unchanged, so a case and its complement share a triangulation, and rotations of the cube map cases onto each other. In their words, these “reduce the problem from 256 cases to 14 patterns”, but their Figure 3 draws 15 cubes, numbered 0 to 14, of which pattern 0 is empty; later papers count 15 (Nielson & Hamann, 1991; Chernyaev, 1995). A non-empty pattern produces one to four triangles.

Mathematical Formulation

On an edge from $\mathbf{p}_a$ to $\mathbf{p}_b$ with values $f_a$ and $f_b$ on opposite sides of $\alpha$, the linear interpolant crosses $\alpha$ at

$$\mathbf{x} = \mathbf{p}_a + t\,(\mathbf{p}_b - \mathbf{p}_a), \qquad t = \frac{\alpha - f_a}{f_b - f_a} \in [0, 1].$$

Lorensen and Cline found higher-degree interpolation of little benefit with at most four triangles per cube. For a smooth field, linear interpolation over an edge of length $h$ errs by $O(h^2)$, so vertices lie $O(h^2)$ from the surface, while edge midpoints give $O(h)$ (reasoning; the example measures both).

On a level set the gradient has no tangential component, so it is normal to the surface. Lorensen and Cline (Eqs. 2–4) estimate it at grid point $(i, j, k)$ by central differences, printed as

$$G_x(i, j, k) = \frac{D(i+1, j, k) - D(i-1, j, k)}{\Delta x}$$

and likewise for $y$ and $z$, with $D$ the sampled density; normalizing removes the missing factor $\tfrac12$. The unit normals are interpolated linearly to each vertex. For a signed distance function, the gradient is the outward normal.

Ambiguities and Holes

Ambiguous faces

A face whose inside corners lie on one diagonal and outside corners on the other has four intersected edges, and the signs do not say how to pair them: the two segments either cut off the inside corners, separating them, or cut off the outside corners, joining the inside corners across the face. Every ambiguous face drawn in Lorensen and Cline’s Figure 3 separates the marked corners, so a table extended to all cases by complement symmetry separates whichever sign is in the minority in each cube (our reading of the figure), and two cubes sharing a face can decide differently and leave a hole. Nielson and Hamann (1991, Figure 3) show one between configuration 6 and the complement of configuration 3, crediting Dürst with the problem; Lewiner et al. (preprint, Figure 3) show a crack between cases 12 and 3.

The asymptotic decider

Nielson and Hamann (1991) assume that the field varies bilinearly over each face. With corner values $B_{00}, B_{01}, B_{10}, B_{11}$ on the unit square, the level curves are hyperbolas, and the pairing follows from comparing $\alpha$ with the interpolant’s value where the asymptotes cross, its saddle point:

$$B(S_\alpha, T_\alpha) = \frac{B_{00} B_{11} - B_{01} B_{10}}{B_{00} + B_{11} - B_{01} - B_{10}}.$$

The paper prints a plus sign in the numerator; expanding the bilinear form at $(S_\alpha, T_\alpha)$ gives the minus sign, as in Chernyaev’s face test (our derivation). The diagonal pair on the saddle’s side of $\alpha$ is joined across the face. The test uses only the face’s four values, so both cubes sharing it decide alike, and Nielson and Hamann added triangulations for the resulting subcases.

Marching cubes 33

Consistent faces remove cracks, but with the same corner signs and face decisions the trilinear interpolant can still differ inside the cube, for example by a tunnel joining opposite corners. Chernyaev (1995) enumerated its topologies, 33 configurations instead of 15, resolving faces with the Nielson–Hamann test and the interior with a further test. Lewiner et al. (2003) completed his interior resolution, computed and tested the 730 subcases of the extended table, and state that their output is a crack-free manifold with the topology of the trilinear interpolant. scikit-image uses their implementation by default.

Marching tetrahedra

Splitting each cube into tetrahedra removes the ambiguity: a linear function on a tetrahedron has a planar level set, giving one triangle or a quadrilateral, provided neighbouring cubes split their shared faces alike (reasoning from linearity). Lewiner et al. (2003) summarize the trade-off: topological consistency and a small table, but many more triangles and weaker geometric accuracy.

Parameters

• Iso-value $\alpha$: zero for a signed distance or TSDF; the mean indicator value at the input samples in Poisson reconstruction (Kazhdan et al., 2006, §4.4); a cross-validated threshold for Occupancy Networks (Mescheder et al., 2019).
• Grid spacing $h$: sets accuracy and output size; features thinner than a cell can vanish or merge (reasoning from sampling).

Complexity

Each cube is visited once at constant cost, $O(N^3)$ for an $N^3$ grid, while the output grows with the surface area: Lorensen and Cline (§6) report triangle counts proportional to area, over 500,000 in some models. For a surface most cubes are empty, and implementations can skip them; Occupancy Networks evaluate a coarse grid, $32^3$ in practice, and subdivide only cubes whose corners disagree before running marching cubes at the final resolution (Mescheder et al., 2019).

Implementation

The NumPy code generates the case table rather than listing 256 entries: it contours each face like 2D marching squares, directs each segment with the inside on its left (seen from outside the cube), chains the segments into loops and fan-triangulates them. On a face with four crossings, a bit chooses whether the inside corners are separated or joined, so one generator builds three tables: complement, separating the minority sign as described above (at four against four it separates the inside corners, a tie the paper’s figure leaves open); separate, always separating the inside corners; and the asymptotic decider. There is no interior test, so loops never become tunnels.

The test surface is a torus with $R = 0.6$, $r = 0.25$ and exact signed distance $f(\mathbf{p}) = \sqrt{(\sqrt{x^2 + y^2} - R)^2 + z^2} - r$, whose area is $4\pi^2 R r$; since $f$ is exact, $|f(\mathbf{v})|$ is a vertex’s distance to the surface. The code also tries edge midpoints, compares with scikit-image 0.25.2, and checks watertightness by counting boundary and non-manifold edges, here and on a random field full of ambiguous faces.

import numpy as np
from functools import lru_cache
from skimage import measure, __version__ as skv

# Cube geometry. Corner c has offset (c & 1, c >> 1 & 1, c >> 2 & 1) along (x, y, z).
CORNER = np.array([[c & 1, c >> 1 & 1, c >> 2 & 1] for c in range(8)])
EDGES = [(a, b) for a in range(8) for b in range(a + 1, 8)
         if np.abs(CORNER[a] - CORNER[b]).sum() == 1]               # the 12 cube edges
EDGE_ID = {frozenset(e): i for i, e in enumerate(EDGES)}
FACES = []                                                          # 6 faces: corners in cyclic order, outward normal
for ax in range(3):
    u, v = [a for a in range(3) if a != ax]
    for side in (0, 1):
        cyc = []
        for du, dv in [(0, 0), (1, 0), (1, 1), (0, 1)]:
            off = np.zeros(3, int); off[ax], off[u], off[v] = side, du, dv
            cyc.append(int(off @ [1, 2, 4]))
        n = np.zeros(3); n[ax] = 1 if side else -1
        FACES.append((cyc, n))
MID = np.array([(CORNER[a] + CORNER[b]) / 2 for a, b in EDGES])   # edge midpoints, local coords


@lru_cache(maxsize=None)
def cell_triangles(config, face_bits):
    """Triangles (as local edge triples) for one cube, generated from its corner signs.

    config: 8-bit mask, bit c set when corner c is inside (f < level).
    face_bits: 6-bit mask; on an ambiguous face (inside corners on a diagonal),
    bit 1 cuts off the inside corners (keeps them separated), bit 0 joins them.
    Each face is contoured like marching squares; the directed segments chain into
    closed loops on the cube's surface, which are fan-triangulated.
    """
    inside = [(config >> c) & 1 for c in range(8)]
    nxt = {}
    for f, (cyc, n) in enumerate(FACES):
        ins = [inside[c] for c in cyc]
        cut = [EDGE_ID[frozenset((cyc[i], cyc[(i + 1) % 4]))] for i in range(4)]
        crossed = [i for i in range(4) if ins[i] != ins[(i + 1) % 4]]
        if len(crossed) == 2:
            segs = [(cut[crossed[0]], cut[crossed[1]])]
            corner_in = cyc[ins.index(1)]
        elif len(crossed) == 4:
            sep_in = (face_bits >> f) & 1
            k0 = [i for i in range(4) if ins[i] == sep_in]          # corners to cut off
            segs = [(cut[(i - 1) % 4], cut[i]) for i in k0]          # the two edges at corner i
            corner_in = None
        else:
            continue
        for a, b in segs:
            # Direct each segment so that the inside lies on its left, seen from outside the cube.
            if corner_in is None:
                i = k0[segs.index((a, b))]
                c, sgn = cyc[i], (1 if sep_in else -1)
            else:
                c, sgn = corner_in, 1
            d, w = MID[b] - MID[a], CORNER[c] - MID[a]
            if sgn * np.dot(n, np.cross(d, w)) < 0:
                a, b = b, a
            nxt[a] = b
    tris = []
    while nxt:                                   # chain segments into loops, fan-triangulate
        start = next(iter(nxt)); loop = [start]
        while nxt[loop[-1]] != start:
            loop.append(nxt.pop(loop[-1]))
        nxt.pop(loop[-1])
        tris += [(loop[0], loop[i + 1], loop[i]) for i in range(1, len(loop) - 1)]
    return tuple(tris)


def marching_cubes(F, h, origin, rule="decider", midpoint=False):
    """Mesh the zero level set of grid F (spacing h); inside is F < 0.

    rule: 'decider'    asymptotic decider on ambiguous faces (bilinear saddle sign),
          'separate'   always keep inside corners apart on ambiguous faces,
          'complement' separate the minority sign (the inside corners at four against four).
    """
    n = np.array(F.shape)
    ins = (F < 0).astype(np.int64)
    sl = lambda c: tuple(slice(o, o + m - 1) for o, m in zip(CORNER[c], n))
    config = sum(ins[sl(c)] << c for c in range(8))
    if rule == "complement":
        nin = sum(ins[sl(c)] for c in range(8))
        bits = np.where(nin <= 4, 63, 0)
    elif rule == "separate":
        bits = np.full(config.shape, 63)
    else:                                         # saddle of the bilinear interpolant on each face
        bits = np.zeros(config.shape, np.int64)
        for f, (cyc, _) in enumerate(FACES):
            p0, p1, p2, p3 = (F[sl(c)] for c in cyc)
            den = p0 + p2 - p1 - p3
            with np.errstate(divide="ignore", invalid="ignore"):
                saddle = (p0 * p2 - p1 * p3) / den
            bits |= (saddle >= 0).astype(np.int64) << f   # saddle outside: inside corners apart
    key = config * 64 + bits
    active = np.flatnonzero((config > 0) & (config < 255))
    keys, inv = np.unique(key.ravel()[active], return_inverse=True)
    cells = np.stack(np.unravel_index(active, config.shape), 1)
    tri_edges = []                                # (cell index, local edge) triples
    for k, kk in enumerate(keys):
        t = np.array(cell_triangles(int(kk // 64), int(kk % 64)), int).reshape(-1, 3)
        idx = np.flatnonzero(inv.ravel() == k)
        if len(t):
            tri_edges.append((np.repeat(idx, len(t)), np.tile(t, (len(idx), 1))))
    ci = np.concatenate([a for a, _ in tri_edges]); le = np.concatenate([b for _, b in tri_edges])
    # Global id of each grid edge: lower grid point + axis, so neighbouring cells share vertices.
    ea = np.array([CORNER[a] for a, b in EDGES]); axis = np.array([np.argmax(CORNER[b] - CORNER[a]) for a, b in EDGES])
    lo = cells[ci][:, None, :] + ea[le]                    # (T, 3, 3) lower grid point of each edge
    gid = (axis[le] * n.prod() + np.ravel_multi_index(np.moveaxis(lo, -1, 0), n))
    uid, faces = np.unique(gid.ravel(), return_inverse=True)
    faces = faces.reshape(-1, 3)
    ax = uid // n.prod(); p = np.stack(np.unravel_index(uid % n.prod(), n), 1)
    q = p.copy(); q[np.arange(len(q)), ax] += 1
    fa, fb = F[tuple(p.T)], F[tuple(q.T)]
    t = np.full(len(p), 0.5) if midpoint else fa / (fa - fb)                  # linear interpolation
    verts = origin + h * (p + t[:, None] * (q - p))
    # Normals as in Lorensen & Cline: central-difference gradient at grid points, interpolated.
    G = np.stack(np.gradient(F, h), -1)
    nrm = (1 - t)[:, None] * G[tuple(p.T)] + t[:, None] * G[tuple(q.T)]
    nrm /= np.linalg.norm(nrm, axis=1, keepdims=True)
    return verts, faces, nrm

# ---------------------------------------------------------------- experiments
def mesh_stats(V, F):
    """Area, signed volume and edge counts of a triangle mesh (vertices merged by position)."""
    V, inv = np.unique(np.round(V, 9), axis=0, return_inverse=True)
    F = inv.ravel()[F]
    F = F[(F[:, 0] != F[:, 1]) & (F[:, 1] != F[:, 2]) & (F[:, 2] != F[:, 0])]  # drop zero-area
    a, b, c = V[F[:, 0]], V[F[:, 1]], V[F[:, 2]]
    cr = np.cross(b - a, c - a)
    he = np.concatenate([F[:, [0, 1]], F[:, [1, 2]], F[:, [2, 0]]])
    und, cnt = np.unique(np.sort(he, axis=1), axis=0, return_counts=True)
    nv = len(np.unique(F))
    return dict(area=0.5 * np.linalg.norm(cr, axis=1).sum(),
                vol=np.einsum("ij,ij->i", a, np.cross(b, c)).sum() / 6,
                tris=len(F), boundary=int(np.sum(cnt == 1)), nonmanifold=int(np.sum(cnt > 2)),
                chi=nv - len(und) + len(F),
                min_angle=np.degrees(min_angles(a, b, c)))

def min_angles(a, b, c):
    def ang(p, q, r):
        u, v = q - p, r - p
        return np.arccos(np.clip(np.einsum("ij,ij->i", u, v) /
                                 (np.linalg.norm(u, axis=1) * np.linalg.norm(v, axis=1) + 1e-300), -1, 1))
    return np.minimum(np.minimum(ang(a, b, c), ang(b, c, a)), ang(c, a, b))

R_T, r_T = 0.6, 0.25                                   # torus radii; area 4 pi^2 R r
def torus(p):
    return np.sqrt((np.hypot(p[..., 0], p[..., 1]) - R_T) ** 2 + p[..., 2] ** 2) - r_T
def torus_normal(p):
    rho = np.hypot(p[..., 0], p[..., 1])
    q = np.stack([p[..., 0] * (1 - R_T / rho), p[..., 1] * (1 - R_T / rho), p[..., 2]], -1)
    return q / np.linalg.norm(q, axis=-1, keepdims=True)
AREA = 4 * np.pi ** 2 * R_T * r_T

print(f"torus R = {R_T}, r = {r_T}: area {AREA:.4f}, genus 1 (chi = 0)")
print("grid  h      variant    tris    |f(v)| mean   max      area err  normal err  min angle<5deg  bnd  chi")
for n in (17, 33, 65, 129):
    g = np.linspace(-1, 1, n); h = g[1] - g[0]
    F = torus(np.stack(np.meshgrid(g, g, g, indexing="ij"), -1))
    runs = [("ours", *marching_cubes(F, h, -1.0)),
            ("midpoint", *marching_cubes(F, h, -1.0, midpoint=True))]
    for m in ("lorensen", "lewiner"):
        V, T, N, _ = measure.marching_cubes(F, 0.0, spacing=(h, h, h), method=m)
        runs.append((m, V - 1.0, T, -N))                # skimage normals follow the descending field
    for name, V, T, N in runs:
        s = mesh_stats(V, T)
        d = np.abs(torus(V))
        ne = np.degrees(np.arccos(np.clip((N * torus_normal(V)).sum(1), -1, 1)))
        print(f"{n:4d}  {h:.4f} {name:9s} {s['tris']:6d}  {d.mean():.1e}  {d.max():.1e}  "
              f"{(s['area'] - AREA) / AREA:+8.2%}  {ne.max():8.2f}  {np.mean(s['min_angle'] < 5):12.1%}"
              f"  {s['boundary']:4d}  {s['chi']:3d}")

# Ambiguous faces: random values (fixed seeds), closed by a positive border.
print("\nrandom 20^3 field   cracks (boundary edges) / non-manifold edges, seeds 0, 1, 2")
rows = {k: [] for k in ("complement", "separate", "decider", "skimage lorensen", "skimage lewiner")}
for seed in range(3):
    rng = np.random.default_rng(seed)
    F = np.ones((20, 20, 20)); F[1:-1, 1:-1, 1:-1] = rng.uniform(-1, 1, (18, 18, 18))
    for rule in ("complement", "separate", "decider"):
        s = mesh_stats(*marching_cubes(F, 1.0, 0.0, rule=rule)[:2])
        rows[rule].append(f"{s['boundary']:5d}/{s['nonmanifold']}")
    for m in ("lorensen", "lewiner"):
        s = mesh_stats(*measure.marching_cubes(F, 0.0, method=m)[:2])
        rows["skimage " + m].append(f"{s['boundary']:5d}/{s['nonmanifold']}")
for k, v in rows.items():
    print(f"{k:17s}", "  ".join(v))

g = np.linspace(-1, 1, 33)
F = np.linalg.norm(np.stack(np.meshgrid(g, g, g, indexing="ij"), -1), axis=-1) - 0.7
V, T, N, _ = measure.marching_cubes(F, 0.0, spacing=(g[1] - g[0],) * 3)
print(f"\nscikit-image {skv}, sphere SDF, defaults: volume {mesh_stats(V - 1, T)['vol']:+.3f} "
      f"(winding outward), normals pointing outward: {np.mean((N * (V - 1)).sum(1) > 0):.0%}")

Output (about 2 s):

torus R = 0.6, r = 0.25: area 5.9218, genus 1 (chi = 0)
grid  h      variant    tris    |f(v)| mean   max      area err  normal err  min angle<5deg  bnd  chi
  17  0.1250 ours         960  1.6e-03  7.0e-03    -1.48%      2.58          0.8%     0    0
  17  0.1250 midpoint     960  2.6e-02  6.1e-02    +0.69%      3.07          0.0%     0    0
  17  0.1250 lorensen     960  1.6e-03  7.0e-03    -1.48%     12.56          0.8%     0    0
  17  0.1250 lewiner      960  1.6e-03  7.0e-03    -1.48%     12.56          0.8%     0    0
  33  0.0625 ours        4080  3.4e-04  1.9e-03    -0.37%      0.84          5.5%     0    0
  33  0.0625 midpoint    4080  1.3e-02  3.1e-02    +5.38%      0.83          0.0%     0    0
  33  0.0625 lorensen    4080  3.4e-04  1.9e-03    -0.37%      6.38          5.5%     0    0
  33  0.0625 lewiner     4080  3.4e-04  1.9e-03    -0.37%      6.38          5.5%     0    0
  65  0.0312 ours       16912  9.9e-05  4.8e-04    -0.09%      0.22          5.9%     0    0
  65  0.0312 midpoint   16912  6.0e-03  1.6e-02    +5.91%      0.23          0.0%     0    0
  65  0.0312 lorensen   16912  9.9e-05  4.8e-04    -0.09%      3.38          6.1%     0    0
  65  0.0312 lewiner    16912  9.9e-05  4.8e-04    -0.09%      3.38          6.1%     0    0
 129  0.0156 ours       68896  2.5e-05  1.2e-04    -0.02%      0.06          5.8%     0    0
 129  0.0156 midpoint   68896  2.9e-03  7.8e-03    +7.02%      0.06          0.0%     0    0
 129  0.0156 lorensen   68896  2.5e-05  1.2e-04    -0.02%      1.80          6.1%     0    0
 129  0.0156 lewiner    68896  2.5e-05  1.2e-04    -0.02%      1.80          6.1%     0    0

random 20^3 field   cracks (boundary edges) / non-manifold edges, seeds 0, 1, 2
complement         3364/0   3320/0   3336/0
separate              0/0      0/0      0/0
decider               0/0      0/0      0/0
skimage lorensen      0/0      0/0      0/0
skimage lewiner       0/0      0/0      0/0

scikit-image 0.25.2, sphere SDF, defaults: volume +1.430 (winding outward), normals pointing outward: 0%

Every torus mesh is closed, with no boundary edges and $\chi = 0$, as genus 1 requires. With interpolation, each halving of $h$ divides the mean vertex distance by about four, to $2.5 \times 10^{-5}$, and the area error falls from $-1.48\%$ to $-0.02\%$; with midpoints the distance only halves and the area error grows to $+7\%$, as the vertices zigzag around the surface. scikit-image matches our counts, distances and areas. Interpolated central-difference normals err by at most $2.6^\circ$, falling to $0.06^\circ$; scikit-image’s returned normals, documented only as calculated from the data, err more, $12.6^\circ$ to $1.8^\circ$. The triangle count grows fourfold per halving, to 68,896, and on the three finer grids about 6% of triangles have an angle below $5^\circ$; midpoint meshes have none.

On the random field, the complement-symmetric table leaves 3,320–3,364 cracked edges over the three seeds; the other rules and both scikit-image methods leave none. Given a signed distance with default settings, scikit-image winds triangles outward but returns inward normals, along the descending field.

In skimage.measure.marching_cubes, level defaults to the mean of the volume’s minimum and maximum, so pass 0 for a signed distance; spacing, step_size, allow_degenerate, mask and gradient_direction are the other options. The default method='lewiner' is documented as an improved marching cubes 33 ported from Lewiner’s C++ code; method='lorensen' selects the original. As of October 2026 the newest release, 0.26.0, documents the same signature.

Properties and Behavior

• Watertightness. If the level set is closed and inside the grid, and shared faces are decided consistently, every edge belongs to exactly two triangles. A surface that leaves the grid stays open at the border; padding the volume with outside values closes it.
• Mesh quality. Where the surface passes near a grid point, the crossings on its edges crowd together into slivers (reasoning; 6% in the example).

Limitations

• Sharp features. Vertices restricted to grid edges cut across corners and creases, which come out rounded at any resolution. Dual contouring (Ju et al., 2002) stores the exact intersection $\mathbf{p}_i$ and normal $\mathbf{n}_i$ on each crossed edge, places one vertex per cube at the minimizer of $E(\mathbf{x}) = \sum_i (\mathbf{n}_i \cdot (\mathbf{x} - \mathbf{p}_i))^2$, and joins the four vertices around each crossed edge into a quad, which can reproduce sharp edges and corners without explicitly detecting features.
• Thin triangles. Surface nets (Gibson, 1998), made for binary segmented data, place one node per surface cube, link neighbouring nodes and relax them while keeping each in its cube. Ju et al. note that such dual meshes typically have better aspect ratios, since vertices move freely inside cubes.

Use in 3D Vision

• Depth fusion. Curless and Levoy (1996, §5.2.3) meshed their fused range images with marching cubes and a table that resolves ambiguous cases; TSDF fusion covers the fusion and the raycasting alternative.
• Poisson surface reconstruction. Kazhdan et al. (2006, §4.4) adapt marching cubes to their octree, avoiding cracks between cells of different sizes by projecting the finer faces’ iso-curves onto the coarser face.
• Neural fields. DeepSDF’s zero level set can be raycast or meshed by marching cubes (Park et al., 2019); Occupancy Networks refine their marching cubes mesh with the network’s gradients (Mescheder et al., 2019).