Marching Cubes

2 min read•1 page

Voxels are great for physics and storage, but they look terrible. To make a game look realistic, we need to convert the blocky grid back into a smooth, continuous polygon mesh. The industry standard algorithm for this is called Marching Cubes.

The Isosurface:

1. Every corner of our grid has a number (a scalar value). E.g. Density, Temperature, or Distance. 2. We pick a magic number called the Isosurface Level (e.g., 0.5). 3. We want to draw a 3D skin exactly where the grid's value equals 0.5. 4. The algorithm "marches" through every single cube in the grid, one by one. 5. If it finds a cube where the skin should pass through, it generates 1 to 5 tiny triangles inside that cube.

python

1# Scalar Fields to Meshes
2def extract_isosurface(voxel_grid, threshold):
3    """The magic behind Medical MRI scans and Minecraft mods."""
4    
5    mesh = Mesh()
6    
7    # 1. Loop through every single box in the grid
8    for voxel in voxel_grid:
9        # 2. Does the surface pass through this specific box?
10        if voxel_contains_surface(voxel, threshold):
11            
12            # 3. If yes, generate triangles for just this box!
13            triangles = marching_cubes(voxel)
14            mesh.add(triangles)
15            
16    return mesh

View (Voxels / Wireframe / Isosurface)

0.00

2 min read•1 page

To figure out what the triangles should look like inside a single cube, we only need to look at the cube's8 Corners. The algorithm uses clever Binary Math (Bitwise Operations) to classify the cube lightning-fast.

The 8-Bit Index:

1. A cube has exactly 8 corners. 2. We check each corner. If its value is less than the Threshold, it is "Inside" the solid object (Bit = 1). If it is greater, it is "Outside" in the empty air (Bit = 0). 3. We pack these 8 boolean true/false values into a single 8-bit binary number (e.g., 10010110). 4. An 8-bit number has exactly 256 possible values (from 0 to 255). 5. This means there are exactly 256 possible ways a surface can intersect a cube!

python

1# Classifying the 8 Corners
2def get_cube_index(cube_corners, threshold):
3    """Generates an 8-bit binary number (0 to 255)."""
4    
5    cube_index = 0
6    
7    # Check all 8 corners of the cube
8    for i in range(8):
9        corner_value = cube_corners[i].value
10        
11        # If the corner is INSIDE the surface
12        if corner_value < threshold:
13            # Flip the i-th bit to 1!
14            # Example: 1 << 0 = 1, 1 << 1 = 2, 1 << 2 = 4...
15            cube_index |= (1 << i)
16            
17    return cube_index

Corner Scenarios

0.00

2 min read•1 page

Calculating exactly where a plane intersects a cube requires heavy 3D math. If we do that for 1 million cubes 60 times a second, the computer will crash. Instead, we cheat.

Pre-computed Answers:

1. Because there are exactly 256 possible scenarios, we can sit down with a piece of paper and figure out all 256 answers ahead of time. 2. We type these answers into a giant array (a Lookup Table) in our code. 3. The computer generates the 8-bit index (e.g., 142). 4. Instead of doing math, the computer just reads Array[142] and instantly gets the list of cube edges that need to be connected to form the triangles!

python

1# The Lookup Table (LUT)
2def lookup_triangles(cube_index):
3    """Avoids calculating geometry in real-time."""
4    
5    # We pre-calculate all 256 possible triangle configurations
6    # and store them in a massive Array/Dictionary when the game loads.
7    
8    # TRIANGLE_TABLE[0] = [] # Completely outside, no triangles
9    # TRIANGLE_TABLE[255] = [] # Completely inside, no triangles
10    # TRIANGLE_TABLE[1] = [ [0, 8, 3] ] # One corner cut off
11    
12    # We simply use our 8-bit index to look up the answer instantly!
13    intersected_edges = TRIANGLE_TABLE[cube_index]
14    
15    return intersected_edges

Table Index (0-255)

1.00

3 min read•1 page

The Lookup Table tells us which edges the skin passes through. But it doesn't tell us whereon the edge it passes through. If we just place vertices at the exact center of every edge, our mesh will look jagged (like Minecraft).

Linear Interpolation (Lerp):

1. We look at the two corners that make up an intersected edge. 2. Corner A has a value of 0.1. Corner B has a value of 0.9. 3. Our Isosurface Threshold is 0.5. 4. Since 0.5 is exactly halfway between 0.1 and 0.9, we place the vertex exactly in the middle of the physical edge! 5. If Corner A was 0.4 instead, the threshold (0.5) would be much closer to Corner A. We slide the vertex closer to Corner A. This creates perfectly smooth curves!

python

1# Finding the Exact Intersection
2def interpolate_edge(corner_A, corner_B, threshold):
3    """Slides the vertex along the edge for a perfect fit."""
4    
5    # 1. Get the scalar values at both ends of the edge
6    val_A = corner_A.value
7    val_B = corner_B.value
8    
9    # 2. How far along the edge is our threshold?
10    # e.g., if A=0.2, B=0.8, and threshold=0.5, then ratio = 0.5 (exactly halfway)
11    ratio = (threshold - val_A) / (val_B - val_A)
12    
13    # 3. Use that ratio to mix the XYZ coordinates
14    x = corner_A.x + ratio * (corner_B.x - corner_A.x)
15    y = corner_A.y + ratio * (corner_B.y - corner_A.y)
16    z = corner_A.z + ratio * (corner_B.z - corner_A.z)
17    
18    return Point(x, y, z)

Value at Corner B

50.00

2 min read•1 page

We now have the exact XYZ coordinates of all the intersection points floating on the edges of our cube. The final step inside this cube is to connect those dots to form physical polygon Triangles.

Winding Order:

1. The Lookup Table doesn't just throw random edges at us; it provides them in a very specific order. 2. It groups the edges into sets of three. 3. We take the first 3 interpolated points and connect them to form Triangle 1. 4. We take the next 3 points and form Triangle 2. 5. The order of the 3 points ensures the Triangle's "Normal" points outward (Counter-Clockwise winding).

python

1# Connecting the Dots
2def generate_cube_geometry(intersected_edges, corner_values):
3    """Turns the edge dots into solid triangles."""
4    
5    triangles = []
6    
7    # The Lookup Table returns edges in sets of 3.
8    # e.g., [Edge_0, Edge_8, Edge_3, Edge_1, Edge_9, Edge_2]
9    
10    # We step through the list 3 edges at a time
11    for i in range(0, len(intersected_edges), 3):
12        
13        # Interpolate exact XYZ for the 3 edges
14        pt1 = interpolate_edge(intersected_edges[i])
15        pt2 = interpolate_edge(intersected_edges[i+1])
16        pt3 = interpolate_edge(intersected_edges[i+2])
17        
18        # Connect them to form a triangle!
19        triangles.append(Triangle(pt1, pt2, pt3))
20        
21    return triangles

Scenarios (1 / 2 / 4 Triangles)

0.00

3 min read•1 page

If we just stop at generating triangles, the mesh will look faceted and low-poly, because the lighting will be perfectly flat across each triangle. We want it to look buttery smooth. We do this by calculatingVertex Normals based on the grid data.

Central Differences:

1. A Normal is an arrow pointing straight away from a surface. 2. Instead of calculating the Normal of the flat triangle, we calculate it using the 3D Grid data (the Scalar Field). 3. We look at the grid values immediately to the Left and Right of our vertex, and find the difference. We repeat for Up/Down and Forward/Back. 4. This is called calculating the Gradient. The gradient always points in the exact direction the scalar field is changing the fastest. 5. We assign this smooth gradient vector to the vertex, and the graphics card uses it to fake perfectly smooth lighting!

python

1# Calculating Smooth Normals
2def calculate_normal(x, y, z):
3    """How to make a jagged mesh look perfectly smooth."""
4    
5    # 1. We look at the grid values immediately around our point
6    # We check +X and -X, +Y and -Y, +Z and -Z
7    
8    # 2. We calculate the difference (Gradient)
9    nx = value_at(x - 1, y, z) - value_at(x + 1, y, z)
10    ny = value_at(x, y - 1, z) - value_at(x, y + 1, z)
11    nz = value_at(x, y, z - 1) - value_at(x, y, z + 1)
12    
13    # 3. Combine them into a Vector and Normalize it (length = 1.0)
14    normal = Vector3(nx, ny, nz).normalize()
15    
16    return normal

Shading (Flat / Smooth)

0.00

2 min read•1 page

We have reached the end of the pipeline. We started with a messy Point Cloud, discretized it intoVoxels, mapped it to a solid Scalar Field, and finally used Marching Cubesto generate a brand new, perfectly manifold, smooth 3D polygon mesh.

The Geometry Pipeline:

1. This exact pipeline is used in Medical MRI visualization to generate 3D bones from 2D slices. 2. It is used in Fluid Dynamics to render splashy water in video games. 3. It is used in 3D Scanning to repair holes and artifacts in raw LiDAR data. 4. By tweaking the Grid Resolution and the Isosurface Threshold, you have complete control over the detail and "thickness" of the final geometry.

python

1# Full Marching Cubes Pipeline
2def generate_bunny_mesh(voxel_grid, threshold):
3    """From blocky Voxels back to a smooth Mesh."""
4    
5    mesh = Mesh()
6    
7    # Run the algorithm on the grid
8    for voxel in voxel_grid:
9        # Get the 8 corner values
10        corners = get_corners(voxel)
11        
12        # Get the 8-bit index
13        idx = get_cube_index(corners, threshold)
14        
15        # Look up which edges are cut
16        edges = lookup_triangles(idx)
17        
18        # Interpolate exact points and form triangles
19        triangles = generate_cube_geometry(edges, corners)
20        
21        mesh.add(triangles)
22        
23    # Calculate smooth lighting normals
24    mesh.calculate_smooth_normals()
25    
26    return mesh

Reconstruction Detail

0.00