Spotlight #60 – Computational Geometry — Sphere Packing, Space-Filling Curves & Triangulation

I. Sphere Packing — Kepler’s Conjecture and the FCC Lattice

In 1611, Johannes Kepler asserted that the most efficient way to pack equal spheres is the arrangement used by greengrocers stacking oranges: hexagonally-close-packed layers, each offset to sit in the hollows of the layer below. This gives a packing density of:

η = π / (3√2)  ≈  0.7405

That is, about 74% of space is occupied by spheres, leaving 26% void. The identical density is achieved by the face-centred cubic (FCC) lattice and the hexagonally close-packed (HCP) arrangement; the difference is how successive layers are stacked (ABCABC vs ABABAB). Kepler's conjecture sat unproved for nearly four centuries. The challenge was not just geometric intuition but rigorously ruling out every possible irregular arrangement that might secretly pack more densely.

The Hales-Ferguson Proof and Flyspeck

Thomas Hales and Samuel Ferguson announced a proof in 1998 after reducing the conjecture to a finite but enormous linear programming problem, then solving it computationally. The submitted manuscript, accompanied by gigabytes of computer code, was accepted in 2005 by Annals of Mathematics — but the referees famously reported they were "99% certain" of correctness while unable to verify every computational step by hand. This motivated the Flyspeck project: a decade-long effort to produce a fully computer-verified formal proof, completed in 2014 using the HOL Light and Isabelle proof assistants. Flyspeck was a landmark in formal mathematics, demonstrating that a major open problem could be settled to absolute logical certainty by machine.

Random packing — simply pouring spheres into a container — achieves a maximum density of roughly 0.64 (the "random close-packing" limit). The gap between 0.64 and 0.7405 has enormous practical consequences: it determines the porosity of pharmaceutical tablets, the permeability of soil, and the conductivity of composite materials. The simulation allows you to compare these regimes directly.

II. Hilbert Curve — A Line That Fills a Square

In 1890, Giuseppe Peano shocked the mathematical world by constructing a continuous curve that passed through every point in the unit square — a curve with a two-dimensional image from a one-dimensional parameter. The following year, David Hilbert gave a cleaner, geometrically transparent construction. Each iteration of the Hilbert curve divides the square into four quadrants, connects them with a U-shaped path, and recursively applies the same rule inside each quadrant, rotating and reflecting as needed.

The limiting curve has Hausdorff dimension exactly 2, which is why it fills the plane: it is not a curve in the intuitive sense but a fractal object that blurs the distinction between line and area. The formal statement is that a continuous surjection from [0,1] to [0,1]² exists, disproving the intuition that curves and squares must have different "sizes" of infinity of points.

Data Locality and the Cache

The Hilbert curve's most valuable property is its locality: points that are close together along the 1D parameter tend to be close together in 2D space. The converse is not perfectly true (some nearby 2D points are far apart on the curve), but the correlation is strong enough to be exploited in computing. Geographic information systems use Hilbert curve indices to sort spatial data so that nearby features tend to reside in the same disk block or cache line. When a query retrieves all points within a bounding box, spatially-ordered data causes far fewer cache misses than a naive row-major ordering.

Database engines (including Google's Bigtable and its open-source successors) use Hilbert indices to cluster related rows on disk. Image compression algorithms use 2D-to-1D mappings (related to the Hilbert curve) for quantisation codebooks. Even DRAM memory controllers have experimented with Hilbert orderings to improve multi-threaded spatial access patterns. The idea that a curve's topology can optimise hardware performance is a beautiful example of abstract mathematics finding unexpected engineering application.

III. Ulam Spiral — Hidden Geometry of the Primes

In 1963, Stanislaw Ulam was sitting through a tedious meeting and started doodling: he wrote 1 in the centre, then wound the integers outward in a square spiral, and circled the primes. The resulting diagram, when extended to thousands of integers, shows unmistakeable diagonal streaks — lines along which primes cluster far above the density expected from random chance. Ulam's doodle became a minor sensation when Scientific American published it.

The explanation lies in quadratic polynomials. Every diagonal of the spiral corresponds to a polynomial of the form 4n² + bn + c for fixed b and c. Some polynomials generate primes at exceptionally high rates. Euler's famous formula n² + n + 41 produces primes for all n from 0 to 39 — forty consecutive primes. In the Ulam spiral, this polynomial traces a bright diagonal streak. The Hardy-Littlewood prime conjectures (Conjecture B, 1923) provide a heuristic explanation: the density of primes generated by a polynomial is related to how often it avoids divisibility by small primes, quantified by a product of local factors called the singular series.

Why No Proof Exists

Despite the clear visual pattern, no quadratic polynomial is proven to generate infinitely many primes. Dirichlet's theorem guarantees infinitely many primes in linear progressions an + b (when gcd(a,b) = 1), but the analogue for quadratics — Landau's fourth problem — remains open. The Ulam spiral thus visualises one of the deepest open questions in number theory: the irregular distribution of primes is simultaneously locally patterned and globally mysterious.

IV. Delaunay Triangulation — The Best Triangle Mesh

Given a set of points in the plane, many possible triangulations connect them into non-overlapping triangles. The Delaunay triangulation is the unique one that maximises the minimum interior angle across all triangles — equivalently, it avoids "sliver" triangles as much as geometrically possible. The defining criterion is the empty circumcircle property: for every triangle in the triangulation, its circumscribed circle contains no other point of the set in its interior.

The Delaunay triangulation is the geometric dual of the Voronoi diagram: if you connect the centres of adjacent Voronoi cells (cells that share an edge), you obtain the Delaunay triangulation. This duality means that computing one gives you the other, and it explains why the Delaunay triangulation partitions space in the most "natural" way relative to the given point distribution.

Bowyer-Watson Algorithm and Applications

The Bowyer-Watson algorithm builds the triangulation incrementally. Start with a super-triangle containing all points. Insert each new point p: find all triangles whose circumcircle contains p (the "bad" triangles); remove them to form a polygonal hole; re-triangulate the hole by connecting all its boundary edges to p. The algorithm runs in O(n log n) expected time and is straightforward to implement, which explains its prevalence in mesh generation software.

Applications span a remarkable range. In finite element analysis, element quality (and thus numerical accuracy) is directly linked to minimum angle, making Delaunay the default starting mesh. Terrain modelling software triangulates irregular survey point clouds into TIN (triangulated irregular network) models. Geospatial interpolation methods like natural neighbour interpolation use the Delaunay structure to determine weights. Even computational fluid dynamics uses Delaunay meshes as the starting point before further refinement algorithms (Ruppert's algorithm, Chew's algorithms) enforce angle and density constraints.

V. Marching Squares — Extracting Contours from Scalar Fields

The marching squares algorithm extracts isocontours — lines along which a scalar field equals a specified threshold — from a regular grid. For each grid cell, each of the four corners is either above or below the isovalue. There are 2⁴ = 16 possible corner configurations, and each configuration maps to a specific set of line segments crossing the cell's edges. Storing these mappings in a lookup table, the algorithm "marches" through every cell and assembles the appropriate segments into contour lines.

Two of the 16 cases (cases 5 and 10) are ambiguous: both pairs of diagonal corners have the same sign, so two equally valid contour topologies exist. The standard resolution samples the scalar value at the cell centre: if the centre is above the isovalue, the contour forms two separate curves (the "separated" case); if below, the contour forms a saddle-point crossing. Getting ambiguous cases wrong produces topological errors — contours that self-intersect or leave gaps.

From 2D to 3D: Marching Cubes

The 3D generalisation, marching cubes (Lorensen and Cline, 1987), operates on cubic voxels with eight corners, yielding 2⁸ = 256 configurations reduced to 15 unique cases by symmetry. It is the workhorse algorithm of medical imaging: CT and MRI scanners produce 3D scalar fields (Hounsfield units for tissue density), and marching cubes extracts organ surfaces at specific density thresholds. The resulting triangle meshes can be 3D-printed as anatomical models, imported into surgical planning software, or used for radiation therapy planning. Virtually every isosurface you have ever seen in a scientific visualisation was produced by this algorithm or one of its descendants.