Voronoi Diagrams — Fortune's sweep line, Delaunay and GLSL distance fields

1. Definition and properties

Let P = {p₁, p₂, …, pₙ} be a set of n distinct points (seeds) in ℝ². The Voronoi cell of seed pᵢ is:

The boundaries between adjacent cells are bisectors — each edge of the diagram is the perpendicular bisector of the segment between two seeds. Three cells meet at a Voronoi vertex, which is the circumcentre of the triangle formed by those three seeds.

Key structural facts:

• Convexity: as the intersection of half-planes, every Voronoi cell is convex.
• Empty circumcircle: the circumcircle of any Voronoi vertex passes through exactly three seeds and contains no other seed in its interior.
• Dual graph: the dual of the Voronoi diagram is the Delaunay triangulation (connect seeds whose cells share an edge).
• Statistical regularity: for a uniform Poisson point process, the average cell has exactly 6 sides (Euler + density argument).

2. Naïve O(n²) algorithm — rasterisation

The most direct approach: for every pixel, find the nearest seed and colour the pixel with that seed's colour. This is brute-force nearest-neighbour over the image grid.

// Brute-force Voronoi — O(W·H·n), simple but slow for large n
function voronoiBrute(seeds, width, height) {
  const out = new Uint8ClampedArray(width * height * 4);
  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      let minDist = Infinity, nearest = 0;
      for (let i = 0; i < seeds.length; i++) {
        const dx = x - seeds[i].x, dy = y - seeds[i].y;
        const d = dx*dx + dy*dy;  // skip sqrt — only comparing
        if (d < minDist) { minDist = d; nearest = i; }
      }
      const c = seeds[nearest].color;
      const idx = (y * width + x) * 4;
      out[idx] = c[0]; out[idx+1] = c[1]; out[idx+2] = c[2]; out[idx+3] = 255;
    }
  }
  return out;
}

Fine for small n (≤ 100) or low-res outputs. For high-resolution diagrams with thousands of seeds, this is prohibitively slow — the answer is Fortune's algorithm.

3. Fortune's sweep-line algorithm — O(n log n)

Steven Fortune's 1987 algorithm computes the exact Voronoi diagram in O(n log n) time and O(n) space by sweeping a horizontal line downward across the plane. It maintains two data structures:

• Event queue (priority queue sorted by y): two event types — site events (when sweep line hits a new seed) and circle events (when three parabola arcs cause a Voronoi vertex to form).
• Beach line (balanced BST): a sequence of parabolic arcs that form the frontier of the computed diagram. Each arc is the locus of points equidistant from one seed and the sweep line.

When a site event fires (sweep line reaches seed pᵢ):

1. Locate the arc directly above pᵢ on the beach line.
2. Split it and insert a new parabola arc for pᵢ between the two halves.
3. Record two half-edges for the new Voronoi edge being traced.
4. Check for circle events involving the new arc's neighbours.

When a circle event fires (three consecutive arcs converge to a point):

1. Remove the middle arc from the beach line.
2. Record the Voronoi vertex at the circumcentre of the three seeds.
3. Connect the two half-edges from the removed arc into finished edges.
4. Re-check circle events for the new neighbouring arc triplets.

4. Delaunay triangulation — the dual graph

The Delaunay triangulation DT(P) is the straight-line dual of the Voronoi diagram: connect two seeds with an edge whenever their Voronoi cells share an edge. The defining property is the empty-circumcircle criterion: no seed lies strictly inside the circumcircle of any triangle in DT(P).

The Delaunay triangulation maximises the minimum angle over all possible triangulations of a point set — making it ideal for finite-element meshes and numerical simulation grids where sliver triangles degrade accuracy.

Bowyer-Watson incremental insertion

An easy-to-implement alternative to Fortune's algorithm for computing Delaunay triangulations:

// Bowyer-Watson incremental Delaunay triangulation
// (simplified — omits boundary ("super-triangle") removal)
function bowyer(points) {
  const triangles = [superTriangle(points)];  // initially one huge triangle

  for (const p of points) {
    // Find all triangles whose circumcircle contains p
    const bad = triangles.filter(t => inCircumcircle(t, p));

    // Find the boundary of the polygonal hole
    const boundary = uniqueEdges(bad);  // edges not shared by two bad triangles

    // Remove bad triangles
    for (const t of bad) triangles.splice(triangles.indexOf(t), 1);

    // Re-triangulate hole: connect each boundary edge to p
    for (const [a, b] of boundary) {
      triangles.push({ a, b, c: p });
    }
  }
  return triangles.filter(t => !sharesSuperVertex(t));
}

// Empty-circumcircle test
function inCircumcircle({ a, b, c }, p) {
  const ax = a.x - p.x, ay = a.y - p.y;
  const bx = b.x - p.x, by = b.y - p.y;
  const cx = c.x - p.x, cy = c.y - p.y;
  return (
    ax*(by*(cx*cx+cy*cy) - cy*(bx*bx+by*by)) -
    ay*(bx*(cx*cx+cy*cy) - cx*(bx*bx+by*by)) +
    (ax*ax+ay*ay)*(bx*cy - by*cx)
  ) > 0;
}

The Bowyer-Watson algorithm is O(n²) in the worst case (though O(n log n) in practice with randomised insertion). Once you have the Delaunay triangulation, the Voronoi diagram is its geometric dual — connect the circumcentres of adjacent triangles.

5. Lloyd relaxation — centroidal Voronoi tessellation

A centroidal Voronoi tessellation (CVT) is a special Voronoi diagram where each seed coincides with the centroid of its own cell. Lloyd's algorithm iterates toward a CVT by alternating two steps:

1. Compute the Voronoi diagram of the current seeds.
2. Move each seed to the centroid (weighted average of pixel positions in that cell).

After 10–30 iterations the seeds distribute evenly across the domain, producing the characteristic foam-like, equal-area hexagonal pattern found in biological tissues, geographic districting, and halftoning.

// Raster Lloyd relaxation — one iteration
function lloydStep(seeds, width, height) {
  const sums = seeds.map(() => ({ sx: 0, sy: 0, count: 0 }));

  for (let y = 0; y < height; y++) {
    for (let x = 0; x < width; x++) {
      let best = 0, bestD = Infinity;
      for (let i = 0; i < seeds.length; i++) {
        const dx = x - seeds[i].x, dy = y - seeds[i].y;
        const d = dx*dx + dy*dy;
        if (d < bestD) { bestD = d; best = i; }
      }
      sums[best].sx += x;
      sums[best].sy += y;
      sums[best].count++;
    }
  }
  return seeds.map((s, i) => sums[i].count > 0
    ? { ...s, x: sums[i].sx / sums[i].count,
               y: sums[i].sy / sums[i].count }
    : s
  );
}

6. GLSL distance-field Voronoi

The most visually compelling Voronoi renders are done in fragment shaders, where each pixel computes the nearest seed in parallel on the GPU. The inigo quilez compact variant additionally computes the second-nearest distance, enabling smooth cell borders and interior patterns.

// GLSL Voronoi — passes seeds via texture or uniform array
// Here: procedural seeds generated from grid hash
// Based on Inigo Quilez's voronoise technique

vec2 hash22(vec2 p) {
  p = vec2(dot(p, vec2(127.1, 311.7)),
            dot(p, vec2(269.5, 183.3)));
  return -1.0 + 2.0 * fract(sin(p) * 43758.5453123);
}

// Returns (minDist, secondMinDist, cell index)
vec3 voronoi(vec2 uv, float scale) {
  vec2 p = uv * scale;
  vec2 gv = floor(p);
  float md = 8.0, md2 = 8.0;
  vec2 closest;

  for (int y = -2; y <= 2; y++) {
    for (int x = -2; x <= 2; x++) {
      vec2 cell = gv + vec2(x, y);
      // Jitter seed within cell
      vec2 seed = cell + 0.5 + 0.5 * hash22(cell);
      float d = length(p - seed);
      if (d < md) { md2 = md; md = d; closest = seed; }
      else if (d < md2) { md2 = d; }
    }
  }
  return vec3(md, md2, md2 - md);
}

// Fragment shader usage
void main() {
  vec2 uv = vUv;
  vec3 v = voronoi(uv, 8.0);

  float minD  = v.x;   // distance to nearest seed
  float minD2 = v.y;   // distance to second nearest
  float border = v.z;  // ≈ distance to edge (minD2 - minD)

  // Style: cells with sharp borders
  float edge = smoothstep(0.0, 0.04, border);
  vec3 col = mix(vec3(1.0), vec3(minD * 1.5), edge);
  gl_FragColor = vec4(col, 1.0);
}

By animating the hash22 seed jitter over time (add a uTime offset to the seed position) you get smoothly morphing liquid Voronoi cells — a popular effect in shader demos and UI backgrounds.

7. Applications in science, games and visualisation

• Nearest-neighbour search: the Voronoi cell of a query point q immediately identifies the nearest seed without searching all of them. Combined with the Delaunay graph, exact nearest-neighbour queries run in O(log n).
• Procedural terrain and maps: Voronoi + Lloyd relaxation divides a plane into polygon "plates" or biomes. Assign height/moisture to each cell for instant believable geography (as in Amit Patel's polygon map generator).
• Biological cells: cross-sections of epithelial tissue, soap foam, and bee honeycombs all approximate Voronoi diagrams with the seed at the cell nucleus. Simulation of cell division uses CVT to keep cells uniformly sized.
• Stippling and non-photorealistic rendering: density-weighted CVT produces dot patterns that perceptually match an input image — "Voronoi stippling".
• Coverage and facility location: finding n warehouse locations that minimise the average customer travel distance = centroidal Voronoi tessellation.
• Finite-element meshes: Delaunay triangulation (the Voronoi dual) is the preferred mesh generation strategy because it maximises triangle quality (no sliver triangles), improving numerical stability.
• Fracture simulation: Voronoi cells define the breakage shards in real-time physics engines. Place seeds randomly inside a mesh, slice along the Voronoi boundaries, then simulate the shard rigid bodies with Cannon-es.

8. Extensions: 3D, weighted and power diagrams

• 3D Voronoi: extend naturally to ℝ³ — each seed gets a convex polyhedral cell. Used in computational chemistry (Wigner-Seitz cells), materials science (grain separation in polycrystals), and 3D mesh refinement.
• Power diagrams (Laguerre-Voronoi): each seed pᵢ has a weight wᵢ (e.g. radius of a circle). The power distance replaces Euclidean distance:
  
  pow(x, pᵢ) = d(x, pᵢ)² − wᵢ²
  
  Cells are still convex and the dual is a regular triangulation. Crucial for simulating sphere packing, foam dynamics and the John-Löwner ellipsoid.
• Cone-based Voronoi on GPU: render each seed as a cone (height = distance) from directly above, projecting down onto the framebuffer. The visible cone surface at each pixel gives the nearest seed — an elegant O(n) GPU algorithm for continuous tone Voronoi with any metric that varies with position.
• Geodesic Voronoi: replace ℝ² with a triangle mesh surface. Geodesic distance (shortest path along the surface) defines cell membership. Used for texture atlasing, surface remeshing, and character rig binding.
• Animated and time-varying: move seeds over time (e.g. SPH particles or agent positions). Each frame recompute Voronoi — gives the perception of territory ownership or influence zones in strategy game AI.