Marching Cubes: From Voxels to Smooth Meshes

When someone says "voxel terrain", the mental image is usually Minecraft: grid-aligned cubes, flat faces, that blocky look. But the cubes are not the important part of a voxel system. The important part is that the world is described by a 3D field of data, not by a mesh someone sculpted by hand. A cube is simply the laziest way to draw that field.

The moment you want rounded caves, smooth hills or a deformable planet like the one in No Man's Sky, an intermediate step is needed: turning that field of numbers into a triangle surface the GPU can paint. That step is called isosurface extraction, and the reference algorithm since 1987 is still Marching Cubes.

During my degree, I took on the challenge of studying and building a system to generate this kind of mesh, and working out how to use it in game development. This article covers how it actually works (not just "there are some magic tables"), why interpolation is the only thing separating a smooth mesh from a stair-stepped one, how to implement it in Unity, what breaks the moment you go from a test cube to a world with chunks (splitting the three-dimensional world into cubes) and levels of detail, plus some of the lessons from that same study. So, fair warning: there will be plenty of maths.

1. The Starting Point: a Density Field

Before generating geometry you need a function that, given a point in space, returns a number. That number is usually called density, and its sign decides whether the point is inside or outside the object.

The most common convention: negative inside the rock, positive in the air, zero exactly on the surface. A sphere of radius r centered at the origin is literally length(p) - r. An infinite plane is p.y. Terrain with relief is p.y - noise(p.x, p.z). Subtracting a sphere from the terrain to carve a cave is a max between the terrain and the negated sphere.

1// A whole world's field fits inside a pure function.
2// That is the big appeal of working with densities.
3public static float Sample(float3 p)
4{
5    float terrain = p.y - Noise.FBM(p.x * 0.02f, p.z * 0.02f) * 24f;
6    float cave     = 0.6f - Noise.Simplex3D(p * 0.05f); // tunnels
7    return math.max(terrain, -cave);
8}

That continuous field is sampled on a regular grid. For a chunk of 32x32x32 cells you store 33x33x33 values, because every cell needs the eight values at its corners and cells share corners with their neighbours.

2. The Idea: Marching Cell by Cell

This algorithm walks every cell of the grid independently. It looks at the eight values at its corners and answers a single question: where does the surface cross this cell?

Each corner is either inside (density below the isovalue) or outside (density above). Eight corners with two states give 2^8 = 256 possible configurations. And here is the trick that makes the algorithm viable: those 256 configurations reduce, through rotation and symmetry, to 15 base cases. All the rest are transformations of those 15.

Since the number of configurations is finite and small, there is no need to reason geometrically at runtime. Two tables are precomputed and the inner loop becomes a handful of memory accesses:

• edgeTable[256]: a 12-bit mask indicating which of the cube's 12 edges are crossed by the surface.
• triTable[256][16]: for each configuration, the list of edge indices grouped in threes, each triplet a triangle, terminated by -1.

1// Building the configuration index: one bit per corner.
2int cubeIndex = 0;
3if (density[0] < isoLevel) cubeIndex |= 1;
4if (density[1] < isoLevel) cubeIndex |= 2;
5if (density[2] < isoLevel) cubeIndex |= 4;
6if (density[3] < isoLevel) cubeIndex |= 8;
7if (density[4] < isoLevel) cubeIndex |= 16;
8if (density[5] < isoLevel) cubeIndex |= 32;
9if (density[6] < isoLevel) cubeIndex |= 64;
10if (density[7] < isoLevel) cubeIndex |= 128;
11
12// 0 = whole cell outside, 255 = whole cell inside. Nothing to emit.
13if (cubeIndex == 0 || cubeIndex == 255) continue;

The corner ordering and the edge ordering are completely arbitrary, but they have to match the ordering used to generate the tables. If you copy triTable from somewhere and number your corners differently, you will get a mesh that looks correct from a distance and is full of holes up close. It is the number one mistake in any implementation written from scratch.

3. Interpolation is the Algorithm

Here is the detail most tutorials skim over, and the one that decides whether the result looks professional or amateur.

Once you know an edge is crossed, you have to place a vertex on it. The naive option is the midpoint. It works, it produces a closed and correct mesh... and it yields a surface with a sort of square-bubble look, because every vertex is anchored to fixed grid positions. It is the equivalent of rounding everything to half a voxel.

The correct option is to interpolate linearly using the density values at the two endpoints:

1static float3 VertexOnEdge(float3 pa, float3 pb, float da, float db, float iso)
2{
3    float denom = db - da;
4    // Endpoints with nearly identical density: any t is equally valid,
5    // so it is pinned to 0.5f to avoid dividing by something tiny.
6    if (math.abs(denom) < 1e-6f) return (pa + pb) * 0.5f;
7
8    float t = (iso - da) / denom;
9    return pa + t * (pb - pa);
10}

With this line, the surface stops being glued to the grid and starts following the real field. The visual difference is enormous for a three-line change: a sphere of radius 8 on a 1-unit grid goes from looking like a rounded die to looking like a sphere.

4. Normals: the Gradient, not the Face

With the mesh already generated, the temptation is to compute normals from each triangle's cross product and average them per vertex. That is what Unity's mesh.RecalculateNormals() does, and for static geometry it is fine. For an isosurface it is the worst option available.

The problem is twofold. First, averaged face normals depend on how the triangles were tessellated, so shading changes that follow the grid pattern appear across the surface. Second, if you generate the mesh in separate chunks, each chunk averages only its own triangles and the normals do not match at the borders, which produces perfectly visible lighting seams.

The alternative is to sample the gradient of the density field via central differences. The surface normal is, by definition, the direction in which the field grows fastest:

1static float3 FieldNormal(float3 p, float h = 0.05f)
2{
3    float dx = Sample(p + new float3(h, 0, 0)) - Sample(p - new float3(h, 0, 0));
4    float dy = Sample(p + new float3(0, h, 0)) - Sample(p - new float3(0, h, 0));
5    float dz = Sample(p + new float3(0, 0, h)) - Sample(p - new float3(0, 0, h));
6    return math.normalize(new float3(dx, dy, dz));
7}

It costs six field samples per vertex, which is not free, but it gives continuous normals independent of tessellation and chunking. If the field is cached in a density buffer, it can be approximated with the differences of the already sampled values and then the normal can be interpolated along the edge exactly like the position is.

5. Where the Real Problems Begin

A world is not generated all at once. It is split into chunks so they can be generated in parallel, unloaded by distance and regenerated only where something changes when the player fires the terraforming tool.

The problem is that every chunk has to emit the cells that reach right up to its border, and those cells need the density values of the neighbouring chunk. If each chunk samples only its own range, a one-cell gap is left between adjacent chunks: the classic terrain with a grid of cracks.

The usual solution is the overlapping border. A chunk of 32 cells per side samples a grid of 34x34x34 values (32 cells plus one extra corner per side, plus a guard ring), and only emits triangles for the 32 interior cells. The ring values are recomputed from the field, not copied from the neighbour, so a chunk can be generated without knowing anyone.

1// A chunk is generated without depending on its neighbours: it resamples
2// the field at the border instead of asking anyone for data.
3const int Cells  = 32;
4const int Points = Cells + 3;        // 32 cells + final corner + guard
5float3 origin = chunkCoord * Cells * voxelSize - voxelSize; // shifted by one cell
6
7for (int z = 0; z < Points; z++)
8for (int y = 0; y < Points; y++)
9for (int x = 0; x < Points; x++)
10    density[Index(x, y, z)] = Sample(origin + new float3(x, y, z) * voxelSize);

And here a performance warning is in order. A 32^3 chunk is 32,768 cells, each with eight samples of a field with several octaves of noise. On the main thread that is tens of milliseconds per chunk, enough for visible hitches as the player moves. Chunk generation belongs in the Job System with Burst, or in a compute shader, not in Update.

6. LOD and the Seam Problem

With chunks solved comes the next step up: distant chunks do not need the same grid density as close ones. Doubling the voxel size cuts the cell count to an eighth, so the temptation is obvious.

The problem is that Marching Cubes is not compatible with itself across resolutions. A chunk with 1-unit voxels and its neighbour with 2-unit voxels do not place their vertices in the same spots on the shared face, so holes appear that let you see the sky through them. These are not subtle artifacts, they are cracks.

There are three ways to live with this:

• Skirts: emit a vertical geometry skirt at the edges of every chunk, dropping low enough to cover any possible crack. It is a workaround, it is cheap and it is what many commercial games use. You notice it if you move the camera level with the ground.
• Transvoxel: Eric Lengyel's extension to the algorithm, with an additional set of tables for the transition cells between two resolutions. It is the correct solution and it really does close the cracks, in exchange for considerably more complexity in the generator.
• A single level of detail with large chunks: perfectly valid if the world is not a planet. Many cave games do not need LOD at all.

7. Ambiguous Cases and the Alternatives

Marching Cubes has a flaw known since its publication: some configurations are ambiguous. When two diagonally opposite corners of a face are inside and the other two outside, the surface can connect them in two topologically distinct ways, and the original tables pick one arbitrarily. If two adjacent cells pick incompatible interpretations on their shared face, a hole is left behind.

In practice, with smooth fields and continuous noise, these cases appear rarely and the holes are tiny. If the project needs meshes that are guaranteed closed (for physics simulation, volume computation or 3D printing) there are two paths:

• Marching Tetrahedra: split each cube into 5 or 6 tetrahedra and apply the same reasoning. A tetrahedron only has 16 configurations and none of them is ambiguous, so the mesh is always closed. The price is considerably more triangles and a surface with a slight directional bias inherited from how the cube was cut.
• Dual Contouring: instead of placing vertices on the edges, it places one vertex per cell and positions it by solving a least-squares system with the normals of the intersections. It preserves sharp edges, which is exactly what Marching Cubes rounds off with no way around it. It needs Hermite data (position and normal per intersection) and is noticeably harder to implement well.

8. When Not to Use Marching Cubes

The algorithm is a specific tool, not a mandatory step in every voxel system.

If the game's aesthetic is blocky, generating cube faces with interior-face culling is faster, simpler, easier to texture and easier to map to UVs. Marching Cubes only pays off if the surface has to be smooth.

If the world has intentionally sharp edges (architecture, artificial structures, crystals) Marching Cubes is going to round them off and no parameter prevents it: it is a consequence of placing vertices only on grid edges. That is Dual Contouring's territory.

And if the problem is 2D, the two-dimensional version is called Marching Squares, has 16 cases instead of 256 and fits in an afternoon. It is what sits underneath any contour generator, 2D destructible level maps and most heatmap systems.

Marching Cubes is one of those algorithms you can explain in twenty minutes and master in several weeks. The main loop is short and the tables are published. What consumes the time is everything else: the corner ordering, the chunk guard ring, the gradient normals, the LOD cracks and moving generation off the main thread. None of those things are in the paper.