HomeArticlesGeometry

Voronoi Diagrams: Nearest-Neighbour Space, Made Even

How a scatter of points partitions the plane into regions, and how Lloyd's algorithm relaxes those regions into an even, near-hexagonal tessellation.

mysimulator teamUpdated June 2026≈ 8 min read▶ Open the simulation

Every point in the plane, claimed by its nearest site

Scatter a handful of points, called sites, across a plane. For every other point in that plane, ask which site is nearest. The set of points closest to a given site — its Voronoi cell — is a convex polygon, and the collection of all such cells for all sites is the Voronoi diagram. The boundary between two neighbouring cells is always a straight segment of the perpendicular bisector between their two sites, because that bisector is exactly the set of points equidistant from both.

live demo · random sites relaxing toward centroidal cells● LIVE

Named after Georgy Voronoi, who formalised it in 1908, the construction had already appeared in Descartes's sketches of how space partitions around neighbouring stars, and it shows up constantly in nature: the polygonal cracks of dried mud, the domains of competing bacterial colonies on a petri dish, and the territorial ranges of animals that space themselves out from their nearest rivals are all approximately Voronoi diagrams, because "closest to me, farthest from my neighbour" is exactly the rule the mathematics encodes.

Raster computation: the cheat that works

The textbook way to build an exact Voronoi diagram is Fortune's sweepline algorithm, which runs in O(n log n) time and outputs the polygons directly. For an interactive canvas, a simpler brute-force trick is usually good enough: for every pixel, compute the distance to all n sites and colour it by whichever is closest.

for each pixel (x, y):
  best = infinity; owner = -1
  for i in 0..n-1:
    d = (x - sites[i].x)^2 + (y - sites[i].y)^2   // squared distance, no sqrt needed
    if d < best: best = d; owner = i
  pixel(x, y) = colour[owner]

This is O(width · height · n) per frame, which sounds wasteful but is easily fast enough at moderate resolution and site counts on modern hardware, and it sidesteps every edge case of exact polygon computation — degenerate collinear sites, numerical precision at cell boundaries — for free. It is also trivially parallel across pixels, which is why the same approach is a natural fit for a GPU fragment shader (a "jump flooding" variant computes an approximate Voronoi diagram on the GPU in logarithmic passes rather than one pass per site).

Lloyd's algorithm: relax toward the centroid

A Voronoi diagram built from random sites usually looks uneven — some cells cramped together, others sprawling. Stuart Lloyd's 1957 algorithm (originally developed for quantising analog signals, later rediscovered for mesh generation and computer graphics) fixes this with a beautifully simple iteration: compute the Voronoi diagram, then move each site to the centroid — the center of mass — of its own cell, and repeat.

repeat k times:
  cells = voronoi(sites)
  for i in 0..n-1:
    sites[i] = centroid(cells[i])   // area-weighted center of mass

Each iteration nudges the diagram closer to a centroidal Voronoi tessellation, a fixed point where every site sits exactly at its own cell's centroid. Visually, cells even out in size, become more compact, and edges straighten toward the honeycomb-like hexagonal packing that minimises the average squared distance from any point in a cell to its site — the same optimisation problem that shows up in signal quantisation, where Lloyd's algorithm is also known as Lloyd-Max quantisation, and in k-means clustering, which is Lloyd's algorithm applied to arbitrary-dimensional data instead of the 2D plane.

The Delaunay dual

Draw an edge between every pair of sites whose Voronoi cells share a boundary, and the result is the Delaunay triangulation — the geometric dual of the Voronoi diagram. It has the defining property that no site lies inside the circumcircle of any triangle formed by three others, which makes it the triangulation that avoids thin, needle-like triangles as much as possible for a given point set. Because the two structures are dual, most computational-geometry libraries compute them together in a single pass; if you need one, you effectively get the other for the cost of reading off the connectivity differently.

Where this shows up

Centroidal Voronoi tessellations are the standard starting point for finite-element mesh generation, because well-shaped, evenly sized cells produce numerically stable simulations. Stippling algorithms use weighted Lloyd relaxation to place dots that reproduce an image's tone using point density rather than shading. Path-planning and sensor-placement problems on a map often reduce to a weighted Voronoi partition, since it answers "which resource is each location's nearest" directly. And procedural generation in games leans on Voronoi diagrams constantly, for territory maps, cracked-rock textures and cell-shaded terrain, precisely because a few relaxation passes turn a purely random scatter into something that reads as deliberately, organically designed.

Frequently asked questions

What is the difference between a Voronoi diagram and its Delaunay dual?

The Voronoi diagram partitions the plane into regions closest to each site; the Delaunay triangulation connects sites whose Voronoi cells share an edge. They are geometric duals of each other and can be computed together in the same sweep, which is why most Voronoi libraries also expose the triangulation for free.

Why do relaxed Voronoi cells tend toward hexagons?

A regular hexagonal tiling is the unique arrangement that minimises the average squared distance from points in a region to that region's centroid, for a fixed cell area and a fixed number of neighbours per cell (six). Lloyd's relaxation drives every cell toward being its own centroid, so on a large, boundary-free domain the equilibrium naturally converges toward that hexagonal packing, the same reason honeycombs and dried mud cracks look hexagonal.

Does Lloyd's relaxation always converge?

It converges toward a centroidal Voronoi tessellation in the sense that the total variance of each cell strictly decreases every iteration, but it is a gradient-like descent, not an exact solver: it can slow down near the fixed point and can settle into a locally optimal but globally imperfect configuration, particularly near domain boundaries. A handful of iterations is usually enough to visibly even out a random point set.

Try it live

Everything above runs in your browser — open Voronoi Diagrams and Lloyd's Relaxation and change the parameters while it is running. Nothing is installed, nothing is uploaded, the whole model lives in one tab.

▶ Open Voronoi Diagrams and Lloyd's Relaxation simulation

What did you find?

Add reproduction steps (optional)