Geometry Processing
Overview
Geometry processing encompasses algorithms for creating, analyzing, and manipulating 3D meshes and other geometric representations. It bridges computational geometry, numerical methods, and computer graphics to handle real-world surface data.
Mesh Representations
Triangle Mesh
The most common representation. A mesh consists of:
- Vertex list: Positions (and optionally normals, UVs, colors)
- Face list: Triples of vertex indices defining triangles
Simple but lacks topological adjacency information. Finding neighboring faces or edges requires linear search.
Half-Edge Data Structure
The standard for mesh processing. Each edge is split into two directed half-edges, one per adjacent face.
HalfEdge:
vertex // vertex at the head of this half-edge
face // face to the left
next // next half-edge in the face loop
prev // previous half-edge in the face loop (or computed from next)
twin // opposite half-edge (same edge, other face)
Vertex:
position
halfedge // one outgoing half-edge
Face:
halfedge // one half-edge on this face
Traversal operations in O(1):
- Face vertices: Follow
nextpointers around the face - Vertex neighbors (1-ring): Follow
twin->nextaround the vertex - Edge faces:
halfedge->faceandhalfedge->twin->face - Boundary detection: A half-edge with no twin (or a null-face twin) lies on a boundary
Other Representations
- Winged-edge: Stores both faces and both vertices per edge. More complex traversal.
- Corner table: Compact, index-based. Good for GPU processing.
- Polygon soup: Unconnected triangles with no topology. Output of many modeling tools; requires cleanup.
Mesh Simplification
Motivation
Reduce polygon count for LOD (level of detail), transmission, or real-time rendering while preserving visual quality.
Edge Collapse
The fundamental operation: contract an edge (v1, v2) into a single vertex v', removing two triangles (one on each side of the edge).
Before: After:
v1 v'
/ | \ / \
/ | \ / \
a---+---b a-----b
\ | / \ /
\ | / \ /
v2
Quadric Error Metrics (QEM, Garland & Heckbert 1997)
Assigns a cost to each edge collapse based on the distance from the new vertex to the original planes of adjacent triangles.
For each triangle with plane equation n . p + d = 0, define the quadric matrix:
Q = [a^2 ab ac ad]
[ab b^2 bc bd]
[ac bc c^2 cd]
[ad bd cd d^2]
where (a, b, c, d) = (n_x, n_y, n_z, d) // plane coefficients
The error of placing a vertex at position v is: error(v) = v^T * Q * v
Algorithm:
- Compute Q for each vertex (sum of quadrics of adjacent faces)
- For each edge (v1, v2), compute Q_new = Q_v1 + Q_v2
- Find optimal position v' by minimizing
v'^T * Q_new * v'(solve 4x4 linear system, or fall back to midpoint) - Insert all edges into a priority queue ordered by error
- Repeatedly collapse the lowest-cost edge, update neighbors
Produces high-quality simplified meshes. O(N log N) with a priority queue.
View-Dependent LOD
Select different detail levels based on distance, screen-space size, or importance. Continuous LOD (progressive meshes) allows smooth transitions.
Subdivision Surfaces
Motivation
Start with a coarse control mesh and recursively refine to produce a smooth limit surface. Used extensively in film production (Pixar's RenderMan).
Catmull-Clark Subdivision
Works on arbitrary polygon meshes (produces quad meshes after one step).
Rules per subdivision step:
- Face point: Average of face vertices:
F = (1/n) * sum(v_i)for an n-sided face - Edge point: Average of edge endpoints and adjacent face points:
E = (v1 + v2 + F1 + F2) / 4 - Vertex point: Weighted combination:
V' = (F_avg + 2*E_avg + (n-3)*V) / n
where:
F_avg = average of new face points of adjacent faces
E_avg = average of midpoints of adjacent edges
n = valence (number of adjacent edges)
The limit surface is C2 continuous everywhere except at extraordinary vertices (valence != 4), where it is C1.
Loop Subdivision
Works on triangle meshes only. Splits each triangle into four.
Edge vertex (on an interior edge with vertices a, b and opposite vertices c, d):
E = (3/8)(a + b) + (1/8)(c + d)
Vertex update (valence n, with neighbors p_i):
V' = (1 - n*beta) * V + beta * sum(p_i)
beta = (1/n) * (5/8 - (3/8 + (1/4)*cos(2*pi/n))^2)
Limit surface: C2 except at extraordinary vertices (valence != 6).
Comparison
| Property | Catmull-Clark | Loop | |--------------------|-------------------|---------------------| | Input mesh | Any polygons | Triangles only | | Output mesh | Quads | Triangles | | Regular valence | 4 | 6 | | Limit surface | Bicubic B-spline | Quartic box spline | | Continuity | C2 (regular) | C2 (regular) |
Mesh Smoothing
Laplacian Smoothing
Move each vertex toward the centroid of its neighbors:
v_i' = v_i + lambda * L(v_i)
L(v_i) = (1/|N(i)|) * sum_{j in N(i)} (v_j - v_i) // uniform Laplacian
Where lambda is a step size (0 < lambda < 1) and N(i) is the set of neighboring vertices.
Problem: Uniform Laplacian shrinks the mesh. Solutions:
- Taubin smoothing: Alternate between smoothing (lambda) and inflation (-mu, where mu > lambda)
- Cotangent Laplacian: Weight by cotangent of opposite angles for geometric accuracy:
L(v_i) = (1 / 2*A_i) * sum_{j in N(i)} (cot(alpha_ij) + cot(beta_ij)) * (v_j - v_i)
Where alpha_ij and beta_ij are the angles opposite edge (i,j) in the two adjacent triangles, and A_i is the mixed Voronoi area around vertex i.
Bilateral Mesh Filtering
Preserves features (sharp edges) while smoothing flat regions. Weights based on both spatial distance and normal similarity (analogous to bilateral image filtering).
Parameterization
Flatten a 3D surface to 2D (for texture mapping, remeshing, or analysis).
Goals
- Low distortion: Minimize angle distortion (conformal) or area distortion (authalic)
- Bijectivity: No self-overlapping triangles in the UV domain
Methods
- Tutte's embedding: Fix boundary vertices to a convex polygon, solve for interior positions using uniform or harmonic weights. Guarantees bijectivity for disk-topology meshes.
- LSCM (Least Squares Conformal Maps): Minimize angle distortion. Free boundary.
- ABF (Angle-Based Flattening): Optimize triangle angles directly.
- Seamless parameterization: Cut the mesh into charts, parameterize each, align seams.
Boolean Operations (CSG)
Constructive Solid Geometry
Build complex shapes from primitives using union, intersection, and difference:
A union B: Points in A or B
A intersect B: Points in A and B
A difference B: Points in A but not B
Mesh Boolean Algorithm
- Compute intersection curves between meshes A and B
- Split triangles along intersection edges
- Classify each resulting face as inside/outside each operand
- Select faces based on the operation (union: outside of both + intersection boundary)
Robust implementation is notoriously difficult due to floating-point precision. Exact arithmetic (rational numbers, simulation of simplicity) or SNapping techniques ensure correctness.
Signed Distance Fields (SDFs)
Definition
A scalar field where the value at each point is the signed distance to the nearest surface:
SDF(p) > 0 outside
SDF(p) = 0 on the surface
SDF(p) < 0 inside
Analytic SDFs (Primitives)
Sphere: SDF(p) = |p - center| - radius
Box: SDF(p) = |max(|p| - half_extents, 0)| + min(max(|p_x|-h_x, |p_y|-h_y, |p_z|-h_z), 0)
Torus: SDF(p) = |(|p.xz| - R, p.y)| - r
CSG with SDFs
Boolean operations become simple min/max:
Union: min(SDF_A, SDF_B)
Intersection: max(SDF_A, SDF_B)
Difference: max(SDF_A, -SDF_B)
Smooth union: -log(exp(-k*SDF_A) + exp(-k*SDF_B)) / k (smooth min)
SDF Applications
- Ray marching (sphere tracing) for rendering
- Collision detection
- Morphing between shapes
- Font rendering (storing glyphs as SDF textures for resolution-independent text)
- Level set methods for fluid simulation
Marching Cubes
Algorithm (Lorensen & Cline, 1987)
Extract a triangle mesh from an implicit function (SDF, density field, medical scan data).
- Sample the scalar field on a regular 3D grid
- For each cube (8 neighboring grid samples):
- Classify each corner as inside (value < isovalue) or outside
- The 8-bit classification indexes into a lookup table of 256 cases (15 unique by symmetry)
- The table specifies which edges are intersected and how to triangulate
- Interpolate vertex positions along intersected edges:
t = (isovalue - v0) / (v1 - v0)
position = lerp(p0, p1, t)
- Compute normals from the gradient of the scalar field (central differences or analytic)
Ambiguities and Improvements
The original 15-case table has ambiguous configurations (saddle points) that can produce holes. Solutions:
- Marching Cubes 33: Resolves ambiguities with additional cases and face/interior tests
- Dual Contouring: Place vertices inside cells using the Hermite data (position + gradient). Preserves sharp features. Requires intersection points and normals on edges.
- Surface Nets: Simpler dual method; place one vertex per cell, connect adjacent cells. Smooth but loses sharp features.
Adaptive Methods
- Octree-based: Refine the grid only near the surface. Reduces memory and computation. Requires crack patching at resolution boundaries.
- GPU marching cubes: Compute shader classifies cells, generates vertices with atomic counters or prefix sum for compaction.
Applications
- Medical imaging (CT/MRI isosurface extraction)
- Terrain generation from noise functions
- Fluid simulation surface extraction
- 3D scanning and reconstruction
- Procedural geometry from SDFs