Voronoi Diagram, in depth: bisectors and convex cells, a hand-computed example, Fortune's sweep, working clipping code and robustness traps

Definition and structure

For two sites p and q, the points closer to p form a half-plane bounded by the perpendicular bisector of the segment pq. The cell of p is the set of points closer to p than to every other site, so it is the intersection of the n - 1 half-planes defined by p and each other site. That single sentence gives most of the diagram's properties.

• Cells are convex, because an intersection of half-planes is convex. A cell may be unbounded; exactly the sites on the convex hull of the input have unbounded cells.
• Edges are pieces of bisectors. A point on the edge between cells of p and q is equidistant from both and farther from every other site.
• Vertices are circumcentres. A vertex is equidistant from at least three sites, so it is the centre of a circle through them, and that circle contains no other site in its interior. This empty-circle property is exactly the Delaunay condition.
• Delaunay duality. Connect two sites whenever their cells share an edge and you get the Delaunay triangulation. Either structure can be built from the other in linear time.

The size is linear. For n >= 3 sites that are not all on one line, Euler's formula on the planar graph (with one extra vertex at infinity to close the unbounded edges) gives at most 2n - 5 vertices and 3n - 6 edges. So the average cell has fewer than six edges, no matter how the sites are arranged. That is why storing the full diagram for a million points is practical.

Worked example: three sites by hand

Take three sites: A = (0, 0), B = (4, 0), C = (0, 4). The bisector of A and B is the vertical line x = 2. The bisector of A and C is the horizontal line y = 2. The bisector of B and C passes through their midpoint (2, 2) and is perpendicular to C - B = (-4, 4), so it is the line y = x. All three bisectors meet at V = (2, 2), the circumcentre of the triangle, at distance sqrt(8), about 2.83 from each site.

From V three rays leave. The A-B edge runs straight down along x = 2, away from C. The A-C edge runs left along y = 2, away from B. The B-C edge runs up and to the right along y = x, away from A. Check a point: (1, 1) is at distance 1.41 from A, 3.16 from B and 3.16 from C, so it is in A's cell, and indeed it is below y = 2 and left of x = 2. All three cells are unbounded, as they must be, because all three sites are on the hull. The counts match the bound: one vertex is at most 2(3) - 5 = 1, three edges are at most 3(3) - 6 = 3.

Now add D = (4, 4). D lies exactly on the circle through A, B and C, so four sites are cocircular and the diagram has a single vertex of degree four at (2, 2). The dual Delaunay triangulation of the square is then ambiguous: either diagonal is valid. That degeneracy is the one that breaks naive implementations, and it appears constantly in real data on grids.

How to construct the diagram

There are three practical ways to construct the diagram.

Half-plane intersection per cell. For each site, start with a bounding box and clip it by the bisector half-plane for every other site. That is O(n) clips per cell and O(n^2) overall with simple clipping, or O(n^2 log n) with a proper half-plane intersection. It is easy to get right, handles bounded output naturally, and is fine for a few thousand sites. A spatial index such as a k-d tree cuts the work to the near neighbours that can actually bound the cell.

Fortune's sweep (1986) runs in O(n log n), which is optimal because sorting reduces to it. A horizontal sweep line moves across the plane. Behind it, the boundary between finished and unfinished territory is a beach line made of parabolic arcs, one arc per site that is still influencing the diagram: the points equidistant from a site and the sweep line. Two kinds of event change the beach line. A site event, when the line reaches a new site, inserts a new arc and starts two edges. A circle event, when three consecutive arcs have a circumcircle whose bottom the line reaches, removes the middle arc and emits a Voronoi vertex at the circle's centre. The beach line lives in a balanced tree, events in a priority queue, and each event costs O(log n). The tricky part is invalidating circle events that a later site event makes false.

Through Delaunay. Build the Delaunay triangulation with Bowyer-Watson incremental insertion or a lifting map to a 3D convex hull, then take circumcentres of triangles as Voronoi vertices and connect the centres of adjacent triangles. This is what Qhull, and therefore SciPy, does. It generalizes to higher dimensions, which Fortune's sweep does not, although diagram size grows quickly with dimension.

On GPUs, the jump flooding algorithm computes a discrete Voronoi diagram on a pixel grid in log2(N) passes for an N by N grid, each pixel adopting the nearest seed seen among neighbours at shrinking offsets. It is approximate and raster-only, but it is the right tool for distance fields in graphics.

Code: library call and a bounded-cell clipper

For most work, call a library. SciPy wraps Qhull and returns vertices plus index lists.

import numpy as np
from scipy.spatial import Voronoi

sites = np.array([[0, 0], [4, 0], [0, 4], [5, 5], [2, 7]], dtype=float)
vor = Voronoi(sites)
print(vor.vertices)        # Voronoi vertex coordinates
print(vor.ridge_points)    # pairs of site indices whose cells share an edge
print(vor.ridge_vertices)  # vertex indices of each edge; -1 means "goes to infinity"
print(vor.point_region)    # site index -> region index
print(vor.regions)         # region -> list of vertex indices (may contain -1)

The -1 entries are where most bugs start: an unbounded region has no finite polygon, so plotting or computing its area needs clipping first. When you need bounded cells, a self-contained clipper is short and transparent:

def clip(poly, a, b, c):
    """Keep the part of convex polygon poly where a*x + b*y <= c (Sutherland-Hodgman, one edge)."""
    out = []
    for i in range(len(poly)):
        P, Q = poly[i], poly[(i + 1) % len(poly)]
        fp, fq = a * P[0] + b * P[1] - c, a * Q[0] + b * Q[1] - c
        if fp <= 0:
            out.append(P)
        if fp * fq < 0:                                  # edge crosses the line
            t = fp / (fp - fq)
            out.append((P[0] + t * (Q[0] - P[0]), P[1] + t * (Q[1] - P[1])))
    return out

def voronoi_cells(sites, box):
    """box = (xmin, ymin, xmax, ymax). Returns one bounded convex polygon per site."""
    x0, y0, x1, y1 = box
    cells = []
    for i, (px, py) in enumerate(sites):
        poly = [(x0, y0), (x1, y0), (x1, y1), (x0, y1)]
        for j, (qx, qy) in enumerate(sites):
            if i == j:
                continue
            # closer to p than q  <=>  2(q - p) . z <= |q|^2 - |p|^2
            a, b = 2 * (qx - px), 2 * (qy - py)
            c = qx * qx + qy * qy - px * px - py * py
            poly = clip(poly, a, b, c)
            if not poly:
                break
        cells.append(poly)
    return cells

cells = voronoi_cells([(0, 0), (4, 0), (0, 4)], (-2, -2, 6, 6))
# cell of A = [(-2,-2), (2,-2), (2,2), (-2,2)]: the worked example, clipped to the box

The inequality comes from expanding |z - p|^2 <= |z - q|^2: the |z|^2 terms cancel and what remains is linear in z. That is the whole reason Voronoi cells are polygons. Duplicate sites make a = b = 0 and c = 0, which keeps the whole polygon; deduplicate inputs first so each location has one owner.

What Voronoi diagrams are used for

Nearest-site queries are the obvious use: build the diagram once and a point-location query finds the owning cell in O(log n). In practice a k-d tree answers the same question with less code, so the diagram earns its keep when you need the regions themselves: service areas, catchment populations, coverage gaps.

Lloyd's algorithm alternates between computing Voronoi cells of a set of centres and moving each centre to the centroid of its cell. On a discrete dataset the cells are just assignment to the nearest centre, which is exactly k-means: every k-means iteration computes an implicit Voronoi partition. On a continuous domain the same loop produces centroidal Voronoi tessellations used for meshing, stippling and sensor placement.

Two more uses are worth knowing. The largest empty circle among the sites, for example the spot farthest from every existing store, has its centre at a Voronoi vertex or on the boundary of the region, so you only need to test those candidates. And the Euclidean minimum spanning tree is a subgraph of the Delaunay triangulation, which turns an O(n^2) problem into O(n log n).

Failure modes

• Cocircular and collinear sites. Four sites on one circle give a degree-four vertex; sites all on one line give parallel edges and no vertices at all. Library calls may raise an error or perturb inputs; Qhull's joggle option perturbs inputs to escape degeneracy.
• Floating-point orientation errors. Deciding which side of a bisector a point lies on, or whether a circle event is real, uses determinants that cancel catastrophically for nearly degenerate inputs. Use robust predicates, as described in 2D geometry algorithms, or integer coordinates.
• Unbounded cells treated as polygons. Taking areas or drawing regions with a -1 vertex index gives garbage. Clip to a bounding box or add far-away dummy sites.
• Wrong metric. Road distance, travel time or geodesic distance on a sphere do not give straight bisectors. For lat-long data over large areas use a spherical Voronoi implementation or project to a local planar coordinate system first.
• Duplicate sites. Two identical sites have no bisector, and some implementations fail or silently drop one. Deduplicate before building.

Trade-offs

Weighted variants change the geometry: multiplicatively weighted diagrams have circular edges, and power diagrams, where each site has a radius, keep straight edges but may leave a site with an empty cell. Pick the variant from the question being asked, not the library at hand.

What to do next

1. Reproduce the three-site example with SciPy and confirm the single vertex at (2, 2) and three infinite ridges.
2. Implement the clipping method above and check its cells against SciPy's on random points.
3. Add D = (4, 4) and see how your code and SciPy handle the cocircular case.
4. Deduplicate and, if needed, project your real data before building the diagram.
5. Decide whether you need regions (build the diagram) or only nearest-site queries (use a k-d tree).
6. Try Lloyd relaxation on random points for ten iterations and watch the cells become regular.