Geometry: Curves, Surfaces & Meshes

A polygon mesh represents a surface as a set of vertices, edges, and polygonal faces (usually triangles or quads). The naive representation — a flat array of triangles, each storing three vertex positions — is called an indexed face list (or simply an OBJ-style mesh). It is compact and cache-friendly for rendering, but terrible for most geometric queries.

Consider: "given vertex v, list all its neighboring vertices." In an indexed face list, you must scan every triangle to find ones that reference v — O(F) time, where F is the face count. For a mesh with 100,000 faces, that is 100,000 comparisons per vertex per query. Geometry processing algorithms (smoothing, remeshing, simulation) do millions of such queries; the naive representation becomes a bottleneck.

Half-edge (doubly-connected edge list, DCEL) solves this by augmenting the mesh with explicit adjacency pointers. Every undirected edge is split into two half-edges, each pointing in the opposite direction along the same geometric edge. Each half-edge stores:

• vertex: the vertex it points to
• face: the face to its left
• next: the next half-edge around the same face (CCW)
• prev: the previous half-edge around the same face
• twin: the opposite half-edge (across the edge)

With these six pointers, any local adjacency query runs in O(degree) time — constant for typical meshes where vertex degree is bounded (usually 6-8 for triangle meshes).

Winged-edge is an older alternative that stores 4 face-adjacent half-edges per edge. Corner tables (or 'star' representations) are a compact cache-friendly alternative that stores only two integer arrays: the vertex array and the triangle array with opposite-vertex pointers.

For triangle meshes specifically, an important invariant is the Euler formula for manifolds: V - E + F = 2 - 2g, where V = vertices, E = edges, F = faces, and g = genus (number of holes). For a sphere (g=0): V - E + F = 2. For a torus (g=1): V - E + F = 0. This formula is a powerful sanity check and also governs average vertex degree: for a closed triangle mesh, average degree ≈ 6.

# Half-edge mesh sketch in Python
from dataclasses import dataclass, field
from typing import Optional

@dataclass
class HalfEdge:
    vertex: 'Vertex'           # vertex this half-edge points TO
    face: Optional['Face'] = None
    next: Optional['HalfEdge'] = None
    prev: Optional['HalfEdge'] = None
    twin: Optional['HalfEdge'] = None

@dataclass
class Vertex:
    position: tuple            # (x, y, z)
    half_edge: Optional[HalfEdge] = None  # any outgoing half-edge

@dataclass
class Face:
    half_edge: Optional[HalfEdge] = None  # any half-edge on this face

def vertex_neighbors(v: Vertex):
    """Iterate all vertices adjacent to v in O(degree) time."""
    start = v.half_edge        # outgoing half-edge from v
    he = start
    while True:
        yield he.vertex        # neighbor
        he = he.twin.next      # rotate around v
        if he is start:
            break

A Bezier curve of degree n is defined by n+1 control points P₀, P₁, ..., Pₙ. The curve itself is a polynomial map from parameter t ∈ [0,1] into 2D or 3D space:

B(t) = Σᵢ₌₀ⁿ C(n,i) · tⁱ · (1−t)^(n−i) · Pᵢ

where C(n,i) are binomial coefficients. The scalar functions Cᵢⁿ(t) = C(n,i) tⁱ(1-t)^(n-i) are called Bernstein basis polynomials. They have beautiful properties: they are all non-negative on [0,1], they sum to 1 for all t, and each one has exactly one maximum in the interval — which means the curve lies inside the convex hull of its control points (the convex hull property). This is crucial for intersection tests and culling.

The de Casteljau algorithm computes B(t) geometrically via repeated linear interpolation, and is the numerically stable way to evaluate Bezier curves:

Given P₀...Pₙ and parameter t:
For r = 1 to n:
  For i = 0 to n-r:
    Pᵢ = (1-t)·Pᵢ + t·Pᵢ₊₁
Return P₀

This is a triangular recursion — at each level you have one fewer point. The beauty: the two sub-arrays generated midway through the recursion are the control points for the left and right halves of the curve. This is the basis of recursive subdivision: you can split a Bezier curve at any t into two Bezier curves, allowing adaptive refinement.

Continuity at joins: A single cubic Bezier can only represent one smooth arc. Real shapes require piecewise Bezier curves. At a join between two curves, you can have:

• C0 (positional): endpoints coincide
• C1 (tangent): the tangent vectors (derivatives) agree in direction and magnitude
• G1 (geometric tangent): directions agree but magnitudes may differ — visually smooth, sufficient for most rendering
• C2 (curvature): second derivatives agree — required for physically-based simulations

For cubic Bezier joins, C1 continuity requires that P₃ of the left curve, the join point, and P₁ of the right curve are collinear and equidistant from the join — the classic "handle" metaphor in Illustrator and Inkscape.

A degree-3 (cubic) Bezier is by far the most common: it has 4 control points, offers an intuitive handle metaphor, and provides enough flexibility for most design work. SVG paths, PostScript fonts, and CSS animations all use cubic Beziers.

def de_casteljau(points, t):
    """Evaluate Bezier curve at t using de Casteljau's algorithm.
    points: list of control points (each a numpy array or tuple)
    Returns a single point on the curve.
    """
    pts = [list(p) for p in points]  # copy
    n = len(pts)
    for r in range(1, n):
        for i in range(n - r):
            pts[i] = [(1 - t) * pts[i][k] + t * pts[i+1][k]
                      for k in range(len(pts[i]))]
    return pts[0]

# Cubic Bezier example
P = [(0,0), (1,2), (2,2), (3,0)]
curve = [de_casteljau(P, t/100) for t in range(101)]

def bezier_split(points, t):
    """Split Bezier at t; returns (left_ctrl_pts, right_ctrl_pts)."""
    pts = [list(p) for p in points]
    n = len(pts)
    left, right = [pts[0]], []
    for r in range(1, n):
        for i in range(n - r):
            pts[i] = [(1-t)*pts[i][k] + t*pts[i+1][k]
                      for k in range(len(pts[i]))]
        left.append(pts[0])
        right.append(pts[n - r - 1])
    right.reverse()
    return left, [left[-1]] + right

The key weakness of Bezier curves is global support: moving any control point changes the entire curve. For design work involving dozens or hundreds of control points (a ship hull, a typeface glyph), this is unacceptable — you want local control.

B-splines (basis splines) achieve local control by defining the curve as a sum of basis functions that each have local support — each basis function is nonzero only over a few parameter intervals. The curve is:

C(t) = Σᵢ Pᵢ · Nᵢ,ₖ(t)

where Pᵢ are control points and Nᵢ,ₖ(t) are B-spline basis functions of order k (degree k-1), defined recursively over a knot vector T = {t₀, t₁, ..., tₘ}:

Nᵢ,₁(t) = 1 if tᵢ ≤ t < tᵢ₊₁, else 0
Nᵢ,ₖ(t) = (t - tᵢ)/(tᵢ₊ₖ₋₁ - tᵢ) · Nᵢ,ₖ₋₁(t)
         + (tᵢ₊ₖ - t)/(tᵢ₊ₖ - tᵢ₊₁) · Nᵢ₊₁,ₖ₋₁(t)

This is the Cox–de Boor recursion. Each Nᵢ,ₖ is nonzero over at most k knot spans — so moving control point Pᵢ only affects k spans. For cubic B-splines (k=4), each control point influences at most 4 parameter spans. This is the local support property.

Knot vectors are the key design tool:

• Uniform: knots equally spaced (e.g., {0,1,2,3,4,5}). Produces uniform parameterization but the curve does not interpolate any control points.
• Clamped (open): knot vector starts and ends with multiplicity-k knots, e.g., {0,0,0,0,1,2,3,3,3,3} for cubic. This forces the curve to pass through the first and last control points — essential for most design applications.
• Repeated interior knots: a knot with multiplicity m reduces continuity at that parameter value by m — from C^(k-1-m). Repeating a knot to multiplicity k-1 produces a C⁰ corner (cusp); to multiplicity k produces interpolation of the adjacent control point.

NURBS (Non-Uniform Rational B-Splines) add a weight wᵢ to each control point, with the curve defined as:

C(t) = Σ wᵢ Pᵢ Nᵢ,ₖ(t) / Σ wᵢ Nᵢ,ₖ(t)

The ratio of polynomials (a rational function) allows exact representation of conic sections — circles, ellipses, parabolas — which ordinary polynomial splines can only approximate. This is why NURBS became the CAD industry standard (IGES, STEP formats): you can represent both a circle and a freeform surface exactly in the same language. NURBS are supported in every major CAD system (CATIA, Solidworks, Rhino, Maya).

import numpy as np

def cox_de_boor(t, i, k, knots):
    """Evaluate B-spline basis N_{i,k}(t) recursively."""
    if k == 1:
        if knots[i] <= t < knots[i+1]:
            return 1.0
        return 0.0
    denom1 = knots[i+k-1] - knots[i]
    denom2 = knots[i+k]   - knots[i+1]
    left  = ((t - knots[i]) / denom1 * cox_de_boor(t, i,   k-1, knots)
             if denom1 != 0 else 0.0)
    right = ((knots[i+k] - t) / denom2 * cox_de_boor(t, i+1, k-1, knots)
             if denom2 != 0 else 0.0)
    return left + right

def bspline_point(t, control_pts, knots, degree):
    """Evaluate B-spline at parameter t."""
    k = degree + 1  # order
    n = len(control_pts)
    result = np.zeros(len(control_pts[0]))
    for i in range(n):
        b = cox_de_boor(t, i, k, knots)
        result += b * np.array(control_pts[i])
    return result

# Clamped cubic B-spline with 5 control points
P = [(0,0),(1,2),(2,-1),(3,2),(4,0)]
# Clamped knot vector: repeat k=4 times at ends
knots = [0,0,0,0,1,2,3,3,3,3]  # m = n+k = 5+4 = 9 intervals
curve = [bspline_point(t, P, knots, degree=3) for t in np.linspace(0,3-1e-9,200)]

There are two fundamental mathematical paradigms for describing surfaces, each with distinct strengths:

Parametric surfaces map a 2D parameter domain (u,v) ∈ [0,1]² into 3D:

S(u,v) = (x(u,v), y(u,v), z(u,v))

Examples:

• Sphere: x = r·sin(v)cos(u), y = r·sin(v)sin(u), z = r·cos(v)
• Bezier patch: S(u,v) = Σᵢ Σⱼ Pᵢⱼ · Bᵢⁿ(u) · Bⱼⁿ(v) — a tensor product of two Bezier curves
• NURBS surface: same idea with B-spline bases and weights

Parametric surfaces are easy to sample (just plug in (u,v) values), easy to render (tessellate by sampling a grid), and easy to animate (deform control points). However, they are hard to intersect with rays (ray-surface intersection requires solving a nonlinear system) and hard to represent closed surfaces of arbitrary topology (a sphere requires a singularity at the poles, or two patches).

Implicit surfaces define the surface as the zero level set of a scalar function F: ℝ³ → ℝ:

{x ∈ ℝ³ : F(x) = 0}

Examples:

• Sphere: F(x,y,z) = x² + y² + z² − r²
• Torus: F(x,y,z) = (√(x²+y²) − R)² + z² − r²
• Signed Distance Function (SDF): F(x) = signed distance to nearest surface point. SDFs are the modern workhorse: F > 0 outside, F < 0 inside, |∇F| = 1 everywhere (the Eikonal equation).

Implicit surfaces excel at inside/outside queries (just evaluate F), CSG (constructive solid geometry) operations — union is min(F₁,F₂), intersection is max(F₁,F₂), difference is max(F₁,−F₂) — and ray-surface intersection via sphere-tracing / ray marching. They are terrible at parameterization and UV mapping.

The fundamental duality: every parametric surface can be converted to implicit via elimination (hard in general, but tractable for polynomials). Going the other way — from implicit to parametric — is the contouring problem, solved approximately by the Marching Cubes algorithm (1987). Marching Cubes samples F on a 3D grid, then triangulates each cell whose corners straddle the zero level set using a 256-case lookup table (reducible to 15 by symmetry). The result is a triangle mesh — an approximate parametric representation of the implicit surface.

Modern neural approaches (NeRF, occupancy networks, SDF networks) use implicit representations because gradient-based optimization is natural in ℝ³, while rendering uses mesh approximations obtained by marching cubes or differentiable variants.

import numpy as np

# --- Implicit sphere SDF ---
def sphere_sdf(p, center, radius):
    """Signed distance from point p to sphere."""
    return np.linalg.norm(p - center) - radius

# CSG union of two spheres
def union_sdf(p, c1, r1, c2, r2):
    return min(sphere_sdf(p, c1, r1), sphere_sdf(p, c2, r2))

# --- Sphere-tracing (ray marching) ---
def ray_march(origin, direction, sdf_fn, max_steps=64, epsilon=1e-4):
    """Find ray-surface intersection via sphere tracing."""
    t = 0.0
    for _ in range(max_steps):
        p = origin + t * direction
        d = sdf_fn(p)
        if abs(d) < epsilon:
            return t, p          # hit!
        t += d                   # safe to step by SDF value
        if t > 100.0:
            break
    return None, None            # miss

# --- Bezier (parametric) surface patch ---
def bezier_patch(u, v, ctrl_pts_16):
    """Evaluate bicubic Bezier patch. ctrl_pts_16: 4x4 grid of 3D points."""
    from functools import reduce
    def bern(n, i, t):
        from math import comb
        return comb(n,i) * t**i * (1-t)**(n-i)
    result = np.zeros(3)
    for i in range(4):
        for j in range(4):
            result += bern(3,i,u) * bern(3,j,v) * np.array(ctrl_pts_16[i][j])
    return result

The dilemma: splines are smooth but topologically rigid (rectangular parameter domains); meshes are topologically flexible but faceted. Subdivision surfaces resolve this tension. The idea: define a surface not by a formula, but by a refinement process. Start with a coarse control mesh. Repeatedly apply subdivision rules — split each face, then reposition every vertex using a weighted average of its neighbors. In the limit, this process converges to a smooth surface.

Catmull-Clark subdivision (1978) works on quad meshes (and generalizes to arbitrary polygons on the first step):

For each subdivision step, compute:

1. Face point Fᵢ = average of face's vertices
2. Edge point Eᵢⱼ = (Vᵢ + Vⱼ + Fᵢ + Fⱼ) / 4, where Fᵢ,Fⱼ are the two adjacent face points
3. Updated vertex for interior vertex V of valence n (number of edges): V' = (F̄ + 2Ē + (n−3)V) / n where F̄ = average of adjacent face points, Ē = average of adjacent edge midpoints

After one step, all faces become quads. Each subsequent step creates a 4× finer quad mesh. The limit surface is C² everywhere except at extraordinary vertices (vertices with valence ≠ 4), where it is C¹. In practice, this is imperceptible — Pixar's production renderer uses Catmull-Clark extensively.

Loop subdivision (1987) works on triangle meshes. For an interior edge vertex:

E' = 3/8 · (V₁ + V₂) + 1/8 · (V₃ + V₄)

where V₁,V₂ are the edge endpoints and V₃,V₄ are the opposite vertices of the two adjacent triangles. For an interior vertex V of valence n:

V' = (1 − nβ)V + β · Σ(neighbors)

where β = (1/n)(5/8 − (3/8 + cos(2π/n)/4)²). The limit surface is C² at regular vertices (valence 6) and C¹ at extraordinary vertices.

Why this is profound: the subdivision rules are linear, so the limit surface is a piecewise polynomial — it can be analyzed exactly. The limit surface for Catmull-Clark coincides with a bicubic B-spline surface at regular regions. This means the mathematical machinery of B-splines applies exactly where the mesh is regular, and the subdivision rules seamlessly handle the irregular extraordinary vertices that B-splines cannot represent.

GPU tessellation: Modern graphics APIs (DirectX 11+, Vulkan, Metal) expose a tessellation stage that can apply subdivision rules on-chip. The application sends a coarse control mesh; the GPU refines it to the desired level of detail at render time — adaptive tessellation based on screen-space edge length. This is how AAA games render smooth characters without storing millions of triangles.

import numpy as np
from collections import defaultdict

def catmull_clark_step(vertices, faces):
    """
    One Catmull-Clark subdivision step.
    vertices: list of np.array([x,y,z])
    faces: list of lists of vertex indices (quads or polys)
    Returns new_vertices, new_faces (all quads).
    """
    new_verts = list(vertices)
    face_pts = []
    edge_to_faces = defaultdict(list)
    
    # 1. Face points
    for fi, face in enumerate(faces):
        fp = np.mean([vertices[v] for v in face], axis=0)
        face_pts.append(fp)
        for i in range(len(face)):
            e = tuple(sorted([face[i], face[(i+1)%len(face)]]))
            edge_to_faces[e].append(fi)
    
    # 2. Edge points
    edge_pt = {}
    for (a,b), adj_faces in edge_to_faces.items():
        mid = (vertices[a] + vertices[b]) / 2
        face_avg = np.mean([face_pts[fi] for fi in adj_faces], axis=0)
        edge_pt[(a,b)] = (mid + face_avg) / 2
    
    # 3. Vertex points (simplified: interior only, valence n)
    new_v = []
    vertex_faces = defaultdict(list)
    vertex_edges = defaultdict(set)
    for fi, face in enumerate(faces):
        for i, v in enumerate(face):
            vertex_faces[v].append(fi)
            nb = face[(i+1)%len(face)]
            vertex_edges[v].add(nb)
    
    for vi, v in enumerate(vertices):
        n = len(vertex_faces[vi])
        F_avg = np.mean([face_pts[fi] for fi in vertex_faces[vi]], axis=0)
        E_avg = np.mean([vertices[nb] for nb in vertex_edges[vi]], axis=0)
        new_v.append((F_avg + 2*E_avg + (n-3)*v) / n)
    
    # 4. Build new quad faces (abbreviated structure)
    # ... (full topology wiring omitted for brevity)
    return new_v, face_pts, edge_pt  # caller assembles new quads

Not all triangle soups are valid meshes. For most geometric operations — subdivision, simulation, boolean operations, 3D printing — you need a manifold mesh: one that locally looks like a 2D disk at every interior point. Concretely:

• Every edge is shared by exactly two faces (no boundary edges in a closed mesh; one allowed for open meshes)
• Every vertex has a single connected fan of faces around it (no bowtie vertices — where two fans share only a vertex)

Violating manifoldness causes half-edge construction to fail, subdivision to produce artifacts, and 3D printing slicers to produce incorrect toolpaths. Tools like MeshLab and libigl can detect and repair non-manifold geometry.

Face orientation determines which side of the surface is the "outside." Consistently oriented meshes have all faces wound in the same rotational sense relative to the outward normal. The face normal for triangle (A,B,C) is:

n̂ = (B−A) × (C−A) / ||(B−A) × (C−A)||

If two adjacent triangles share edge (A,B), one must traverse it as A→B and the other as B→A for the mesh to be consistently oriented. You can check and fix orientation with a flood-fill: pick a seed triangle, mark it consistently, propagate to neighbors via edge parity. If the mesh is non-orientable (like a Möbius strip), this flood-fill will encounter a contradiction.

Vertex normals for smooth shading are computed as a weighted average of adjacent face normals. The weighting scheme matters:

• Uniform: each adjacent face contributes equally — fast but can produce artifacts on irregular meshes
• Area-weighted: weight by triangle area — better for smoothly varying normals
• Angle-weighted: weight by the angle subtended at the vertex — shown by Thürmer & Wüthrich (1998) to be the most accurate approximation to the limit normal of a smooth surface

UV (texture) parameterization is another mesh quality issue: assigning (u,v) coordinates to vertices so texture maps can be applied. A perfectly flat triangle has a trivially correct UV mapping, but on a curved mesh you must cut seams and flatten patches. Minimizing distortion (angle distortion = conformal, area distortion = authalic) is an active research area, with methods ranging from LSCM (Least Squares Conformal Maps) to the spectral method ABF++.

import numpy as np

def face_normal(A, B, C):
    """Outward unit normal of triangle ABC (CCW winding = outward)."""
    n = np.cross(B - A, C - A)
    length = np.linalg.norm(n)
    return n / length if length > 1e-12 else np.zeros(3)

def angle_weighted_vertex_normals(vertices, faces):
    """Compute per-vertex normals weighted by incident angle."""
    normals = np.zeros_like(vertices)
    for face in faces:
        i, j, k = face
        A, B, C = vertices[i], vertices[j], vertices[k]
        fn = np.cross(B-A, C-A)          # unnormalized
        # Angles at each vertex
        def angle(v0, v1, v2):
            a = v1 - v0; b = v2 - v0
            cos_t = np.dot(a,b) / (np.linalg.norm(a)*np.linalg.norm(b) + 1e-12)
            return np.arccos(np.clip(cos_t, -1, 1))
        normals[i] += angle(A, B, C) * fn
        normals[j] += angle(B, C, A) * fn
        normals[k] += angle(C, A, B) * fn
    # Normalize
    norms = np.linalg.norm(normals, axis=1, keepdims=True)
    return normals / (norms + 1e-12)

def is_manifold_edge_check(faces):
    """Return set of non-manifold edges (those not shared by exactly 1-2 faces)."""
    from collections import Counter
    edge_count = Counter()
    for face in faces:
        n = len(face)
        for i in range(n):
            e = tuple(sorted([face[i], face[(i+1)%n]]))
            edge_count[e] += 1
    return {e for e, c in edge_count.items() if c > 2}

Once you have a manifold mesh, you want to process it: denoise a scan, remesh to improve triangle quality, smooth sharp creases, or compute physically meaningful quantities like curvature. The central tool is the discrete Laplacian.

The smooth Laplace–Beltrami operator Δ acts on functions f: S → ℝ on a smooth surface S. Its discrete analog on a triangle mesh is the cotangent Laplacian:

(Lf)ᵢ = (1/Aᵢ) · Σⱼ∈N(i) (cotαᵢⱼ + cotβᵢⱼ)/2 · (fⱼ − fᵢ)

where αᵢⱼ and βᵢⱼ are the two angles opposite edge (i,j) in the adjacent triangles, and Aᵢ is the mixed Voronoi area around vertex i. This discretization converges to the smooth Laplace–Beltrami operator as the mesh is refined.

Laplacian smoothing (naive): move each vertex toward the average of its neighbors: Vᵢ' = Vᵢ + λ · (1/|N(i)|) Σⱼ∈N(i) (Vⱼ − Vᵢ). Simple, but causes shrinkage (the mesh contracts toward its centroid).

Taubin smoothing (1995): alternate between a positive smoothing step (λ > 0) and a negative "inflation" step (μ < 0, |μ| > λ). This eliminates shrinkage while still damping high-frequency noise.

Implicit fairing (Desbrun et al. 1999): solve (I − λL)V' = V. This is unconditionally stable for large λ and produces excellent results at the cost of a sparse linear system solve.

Mean curvature flow: the cotangent Laplacian applied to the vertex positions themselves gives the mean curvature normal vector HN̂. Integrating this flow shrinks the surface toward minimal area (minimal surfaces). This connects mesh processing to differential geometry and has applications in surface reconstruction.

Isotropic remeshing (Botsch & Kobbelt 2004): alternate between edge splits (edges too long), edge collapses (edges too short), edge flips (improve Delaunay condition), and tangential relaxation (Laplacian smoothing in the tangent plane). This produces meshes with near-uniform edge lengths and valence-6 vertices — ideal for subdivision and simulation.

Poisson surface reconstruction (Kazhdan et al. 2006): given an oriented point cloud (positions + normals from a 3D scan), solve a Poisson equation to find an implicit indicator function whose gradient approximates the normal field. Marching Cubes then extracts the surface. This is the standard algorithm in MeshLab and CloudCompare for converting LiDAR scans to meshes.

import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve

def cotangent_laplacian(vertices, faces):
    """Build the cotangent weight matrix L (V x V sparse)."""
    n = len(vertices)
    L = lil_matrix((n, n))
    for tri in faces:
        for k in range(3):
            i, j, l = tri[k], tri[(k+1)%3], tri[(k+2)%3]
            # Angle at l, opposite edge (i,j)
            vi, vj, vl = vertices[i], vertices[j], vertices[l]
            u, v = vi - vl, vj - vl
            cos_a = np.dot(u, v)
            sin_a = np.linalg.norm(np.cross(u, v))
            cot = cos_a / (sin_a + 1e-12)
            L[i, j] -= cot / 2
            L[j, i] -= cot / 2
            L[i, i] += cot / 2
            L[j, j] += cot / 2
    return L.tocsr()

def implicit_fairing(vertices, faces, lam=0.1):
    """One step of implicit Laplacian fairing: solve (I - lam*L) V' = V."""
    L = cotangent_laplacian(vertices, faces)
    n = len(vertices)
    from scipy.sparse import eye
    A = eye(n) - lam * L
    new_verts = np.zeros_like(vertices)
    for d in range(3):    # x, y, z independently
        new_verts[:, d] = spsolve(A, vertices[:, d])
    return new_verts

The Marching Cubes algorithm (Lorensen & Cline, SIGGRAPH 1987) bridges the implicit-to-parametric gap. Given a scalar field F sampled on a regular 3D grid, it extracts a triangle mesh approximating the level set F(x,y,z) = isovalue.

The algorithm processes each grid cube independently. Each of the 8 cube corners is classified as "inside" (F < isovalue) or "outside" (F ≥ isovalue), yielding 2⁸ = 256 configurations. By symmetry and rotation, these reduce to 15 canonical cases. For each case, a lookup table specifies which edges of the cube the surface intersects and how to triangulate them. The intersection point on each edge is computed by linear interpolation:

P = V₁ + (isovalue − F₁)/(F₂ − F₁) · (V₂ − V₁)

This produces a globally consistent mesh because shared edges between adjacent cubes produce the same intersection point.

Limitations:

1. Topology ambiguity: some cases in the original Lorensen-Cline table produce topologically inconsistent meshes ("holes"). Extended tables (33-case, or the asymptotic decider) fix this.
2. Resolution: the output mesh quality depends on grid resolution. Thin features thinner than one voxel are lost.
3. Not watertight: at boundaries of the sampled volume, the mesh may be open.
4. Memory: a 512³ grid requires 128 MB for single-precision floats.

Dual Contouring (Ju et al. 2002) improves on marching cubes by placing one vertex per cube (at the optimal position via QEF — quadric error function minimization) rather than per edge. It better preserves sharp features and handles topology more cleanly.

Marching Tetrahedra subdivides each cube into 5 or 6 tetrahedra, yielding only 2⁴ = 16 cases per tetrahedron — simpler and unambiguous, at the cost of more triangles and a slightly different mesh structure.

Marching Cubes is used in: MRI/CT visualization (isosurfaces of bone density), NeRF and neural SDF reconstruction, fluid simulation surface extraction, and procedural terrain generation (voxel games like Minecraft use a variant).

# Pseudocode for Marching Cubes (one cube)
MC_EDGE_TABLE = [...]  # 256-entry lookup: which edges are cut
MC_TRI_TABLE  = [...]  # 256-entry lookup: edge triplets forming triangles

def march_cube(corners, values, isovalue):
    """
    corners: 8 x 3 array of cube vertex positions
    values:  8-element array of scalar values at corners
    Returns list of triangles (each triangle = 3 xyz positions)
    """
    # Build 8-bit index: bit i set if values[i] < isovalue
    cube_index = 0
    for i in range(8):
        if values[i] < isovalue:
            cube_index |= (1 << i)
    
    if cube_index == 0 or cube_index == 255:
        return []  # fully inside or outside
    
    # Find intersection points on cut edges
    edge_verts = {}  # edge_id -> interpolated position
    cut_edges = MC_EDGE_TABLE[cube_index]
    for edge_id in range(12):
        if cut_edges & (1 << edge_id):
            i, j = EDGE_ENDPOINTS[edge_id]  # vertex indices of this edge
            t = (isovalue - values[i]) / (values[j] - values[i] + 1e-12)
            edge_verts[edge_id] = corners[i] + t * (corners[j] - corners[i])
    
    # Assemble triangles from lookup table
    triangles = []
    tri_list = MC_TRI_TABLE[cube_index]
    for k in range(0, len(tri_list), 3):
        if tri_list[k] == -1: break
        tri = [edge_verts[tri_list[k+d]] for d in range(3)]
        triangles.append(tri)
    return triangles