Voronoi Diagrams: Dividing Space by Nearest Point

The nearest-point rule and what a cell actually is

Given a set of sites (seed points) P = {p₁, …, pₙ} in the plane, the Voronoi cell of site pᵢ is the set of all locations closer to pᵢ than to any other site:

V(pᵢ) = { x : dist(x, pᵢ) ≤ dist(x, pⱼ) for all j ≠ i }

The key insight is that each pairwise constraint dist(x, pᵢ) ≤ dist(x, pⱼ) defines a half-plane — the side of the perpendicular bisector of segment pᵢpⱼ that contains pᵢ. A cell is therefore the intersection of n − 1 half-planes, which makes it a convex polygon (possibly unbounded). This is why the boundaries in a Voronoi diagram are always straight edges meeting at vertices: an edge is a piece of a bisector equidistant from exactly two sites, and a vertex is equidistant from three or more sites — the center of a circle through them with no site inside.

Under Euclidean distance the whole plane tiles perfectly: every point lands in some cell, cells overlap only on their shared boundaries, and the union is the entire plane.

Why the size is linear even though there are n² pairs

A naive count suggests trouble: there are n(n−1)/2 ≈ n²/2 perpendicular bisectors, so one might fear a quadratic-sized diagram. The saving grace is planarity. The Voronoi diagram is a planar subdivision, and Euler's formula for planar graphs (V − E + F = 2) bounds its complexity tightly.

• Faces = n cells (plus the unbounded outer relation).
• Vertices ≤ 2n − 5.
• Edges ≤ 3n − 6.

So the diagram has Θ(n) total complexity, not Θ(n²). Most of those n² bisectors never appear as a full edge — they are pruned away because some third site is closer. This linear size is what makes an optimal Θ(n log n) construction even possible: you cannot beat the Ω(n log n) sorting lower bound, but you also do not have to pay for a quadratic output.

A caveat: this Θ(n) bound is a planar phenomenon. In 3-D a Voronoi diagram can have Θ(n²) complexity, and in d dimensions up to Θ(n^⌈d/2⌉) — the reason algorithms that are cheap in the plane become expensive fast.

Fortune's sweep line: the canonical O(n log n) build

The elegant optimal algorithm is Steven Fortune's sweep line (1986). A naive top-to-bottom sweep fails because a cell's shape depends on sites below the line that haven't been seen yet. Fortune's trick is to track a beach line: the boundary of the region already known to be correct, made of parabolic arcs. Each arc is the locus of points equidistant from one site and the sweep line, so as the line descends the parabolas grow and the arcs' intersections trace out Voronoi edges.

Two event types drive it, held in a priority queue ordered by y:

• Site event: the sweep reaches a new site; a new arc is inserted into the beach line, splitting an existing arc.
• Circle event: three consecutive arcs define a circle whose bottom the sweep is about to reach; the middle arc disappears, and that point becomes a Voronoi vertex (equidistant from three sites).

The invariant: everything above the beach line is finalized and correct. Both queue operations and the balanced-tree updates on the beach line cost O(log n), and there are O(n) events total (each site adds at most one arc; each circle event removes one), giving Θ(n log n) time and Θ(n) space.

A brute-force cell, in code

If you only need one cell and clarity beats asymptotics, clip a bounding box by every bisector half-plane. This is the half-plane intersection definition made literal — O(n²) overall for all cells, but bulletproof and easy to verify:

def cell(i, sites, box):
  poly = box  # start with a big convex bounding polygon
  for j in range(len(sites)):
    if j == i: continue
    # keep the side of the bisector nearest to site i
    mid = midpoint(sites[i], sites[j])
    normal = sites[i] - sites[j]  # points toward i
    poly = clip_halfplane(poly, mid, normal)
  return poly

Here clip_halfplane is Sutherland–Hodgman polygon clipping against one line. Doing this for every i is O(n²), but for a few thousand sites it runs in milliseconds and is the fastest correct thing to write in an interview. Reach for Fortune's algorithm or a library (e.g. scipy.spatial.Voronoi, which lifts to a 3-D convex hull) only when n grows large.

Delaunay: the diagram's inseparable twin

Every Voronoi diagram has a dual graph called the Delaunay triangulation: connect two sites with an edge whenever their Voronoi cells share a boundary. This duality is not a curiosity — it is how most software actually computes things. Many pipelines build the Delaunay triangulation first (it is numerically more robust and directly gives a convex-hull-based construction) and then read off the Voronoi diagram as its dual, since a Voronoi vertex is exactly the circumcenter of a Delaunay triangle.

Delaunay's defining empty-circle property — no site lies strictly inside the circumcircle of any triangle — is the same fact that makes Voronoi vertices equidistant from three sites with no closer fourth. Delaunay triangulations are prized because they maximize the minimum angle over all triangulations of the point set, avoiding thin slivers, which is why they dominate finite-element meshing and terrain modeling.

Practical rule of thumb: if you want a partition of space, think Voronoi; if you want a good triangular mesh or a nearest-neighbor graph, think Delaunay. You are always one dualization step from the other.

Where it's used, and the traps

Voronoi structures answer proximity queries across many fields: the classic post-office problem (which facility is nearest?), John Snow's cholera map, cellular-network coverage, robot motion planning along cell boundaries that stay maximally far from obstacles, natural-neighbor interpolation, and the Lloyd relaxation loop that repeatedly moves each site to its cell's centroid — the geometric heart of k-means clustering and of blue-noise stippling.

Common pitfalls to respect:

• Degeneracies. Four sites that are exactly cocircular create a Voronoi vertex of degree 4; collinear sites yield unbounded, parallel-edged cells. Robust code perturbs inputs or uses exact/adaptive arithmetic — naive floating point produces non-planar garbage here.
• Unbounded cells. Sites on the convex hull have infinite cells. Libraries return "points at infinity" or ridge directions; forgetting to clip them to a bounding box is the #1 rendering bug.
• Wrong metric. Under Manhattan (L₁) distance bisectors bend and cells are no longer the familiar polygons; a weighted (power/Laguerre) diagram can even produce empty cells. "Voronoi" is a family, not one shape.

Voronoi construction and nearest-point strategies compared

Approach	Build time	Query	Best for
Fortune's sweep line	Θ(n log n)	—	Building the full 2-D diagram optimally
Divide & conquer (Shamos–Hoey)	Θ(n log n)	—	Classic optimal build; recursion-friendly
Incremental / half-plane intersection	O(n²) worst	—	Simple code, one cell at a time
Lifting to 3-D convex hull	Θ(n log n)	—	Unifying Voronoi + Delaunay in higher dim
k-d tree (no diagram)	Θ(n log n) build	O(log n) avg	Nearest-neighbor queries without tessellating