Voronoi Diagram (Fortune's Algorithm)

What a Voronoi diagram answers

Drop a handful of points — call them sites — onto a plane: cell towers, coffee shops, weather stations, pixels in a stippling. A Voronoi diagram answers a single question for every other point in the plane: which site is closest? It carves the plane into one convex cell per site, where every point inside a cell is nearer to that site than to any other. The cell boundaries are pieces of the perpendicular bisectors between pairs of sites; three bisectors meeting at a point form a Voronoi vertex equidistant from three sites.

The brute-force construction is honest but slow: for each site, intersect the n − 1 half-planes defined by its bisectors with every other site. That is O(n) work per site, O(n²) overall. Steven Fortune's 1986 sweepline algorithm gets it down to O(n log n), which is optimal — the problem has an Ω(n log n) lower bound because it can sort. The trick is to never look at all pairs at once.

How the sweepline and beach line work

Imagine a horizontal line sweeping from the top of the plane downward. Above the line, every site has already been seen; below it, sites are still unknown. The problem is that a Voronoi edge can be influenced by a site the sweep hasn't reached yet — so you cannot finalize anything just because the sweep passed it.

Fortune's insight is the beach line: a lower envelope of parabolas. For a site already passed by the sweep, the set of points equidistant from that site (the focus) and the sweepline (the directrix) is a parabola. The beach line is the front-most boundary of all these parabolas — the locus of points that are as close to some known site as they are to the still-unexplored region below the sweep. Any point above the beach line already knows its nearest site for good; the sweep can never reveal anything closer.

As the sweepline descends, each parabola widens, and the breakpoints where two parabolas cross slide sideways. Those breakpoints trace exactly the Voronoi edges. The whole algorithm reduces to maintaining the beach line through two kinds of discrete events:

• Site event. The sweep reaches a new site. A new parabola appears as a degenerate vertical ray and immediately splits the arc directly above it into two, inserting itself in the middle. This adds one arc (net +2 arcs, counting the split) and starts two new edges.
• Circle event. An arc gets squeezed between its two neighbors until it shrinks to zero width. The point where it vanishes is equidistant from three sites — the centers of the three arcs — so it is a Voronoi vertex. The disappearing arc is removed and a single edge continues from the vertex.

Site events are all known in advance (sort sites by y). Circle events are computed on the fly: whenever the beach line changes, you check each triple of consecutive arcs to see whether they will converge, and if so, schedule a circle event at the bottom of that circle. Both event types go into a single priority queue keyed by y, and you process them top to bottom.

The two data structures that make it fast

Fortune's algorithm needs two carefully chosen structures, and getting either wrong drags it back to O(n²):

• Event queue — a priority queue (binary heap). Holds site events (known up front) and circle events (added as discovered), ordered by descending y. Each pop and push is O(log n). There are O(n) events, so O(n log n) queue work.
• Beach line — a balanced binary search tree. The arcs are stored left-to-right in a self-balancing tree (red-black or similar), keyed implicitly by x position at the current sweep height. Inserting an arc on a site event, deleting an arc on a circle event, and locating "which arc is directly above this site" are all O(log n). A naive linked list makes the locate step O(n) and kills the bound.

The output is usually a doubly connected edge list (DCEL) — also called a half-edge structure — so that afterward you can walk around any cell or hop between adjacent cells in O(1) per step.

When to reach for Fortune's algorithm

• Static 2D point sets where you need the full diagram once: stippling and dithering, terrain meshing, procedural texture generation, biological cell modeling.
• Nearest-site queries at scale. Build the diagram once in O(n log n), then point-locate in O(log n) per query — far better than scanning all sites each time.
• Anywhere you also want Delaunay. Fortune's algorithm yields the Voronoi diagram, and its dual is the Delaunay triangulation for free.
• Memory-resident workloads where the O(n) DCEL fits in RAM and pointer chasing is fine.

If your points arrive incrementally or move, prefer randomized incremental Delaunay (easier to make dynamic). If you only need an approximate diagram for rendering, the GPU "jump flooding" trick computes a pixel-accurate Voronoi image in a few passes without any of this machinery. And in higher dimensions, switch to building the dual (Delaunay) via convex hull lifting — Fortune's sweepline is inherently 2D.

Fortune's vs other Voronoi constructions

	Fortune's sweepline	Bowyer–Watson (incremental Delaunay)	Divide & conquer	Lifting to 3D hull	Jump flooding (GPU)
Time	O(n log n)	O(n log n) expected, O(n²) worst	O(n log n)	O(n log n)	O(p log p) per pixel grid
What it builds	Voronoi (Delaunay dual)	Delaunay (Voronoi dual)	Either	Delaunay via 3D convex hull	Pixel-approx Voronoi
Incremental / dynamic	No (batch only)	Yes (add points one at a time)	No	No	Recompute each frame
Implementation difficulty	High (events + beach line)	Moderate	Hard (the merge step)	Moderate (needs a hull)	Low (a shader)
Exactness	Exact (with robust predicates)	Exact	Exact	Exact	Approximate
Best for	Static 2D, optimal worst case	Streaming / moving points	Textbook optimality proofs	Higher dimensions	Real-time rendering

In practice, most production libraries (CGAL, Triangle, scipy's Voronoi) build Delaunay incrementally or via Qhull's hull-lifting and read off the Voronoi dual, because incremental code is easier to make robust and dynamic. Fortune's algorithm wins on textbook optimality and on the elegance of producing the Voronoi diagram directly — it is the canonical answer when an interviewer asks for an optimal 2D construction.

What the numbers actually say

• Exactly n site events. One per input site, all known before you start.
• At most 2n − 5 circle events for n ≥ 3, because a planar Voronoi diagram has at most 2n − 5 vertices (Euler's formula on the planar graph). So the queue processes fewer than 3n events total.
• At most 3n − 6 edges and 2n − 5 vertices. The diagram is a planar graph of linear size, so the entire output fits in O(n) memory regardless of how the sites are scattered.
• The beach line never holds more than 2n − 1 arcs — each site contributes arcs, but a single site can own several disconnected arcs at once, which is the most counter-intuitive fact in the whole algorithm.
• For n = 1,000,000 sites, that is roughly 3 million events, each an O(log n) ≈ 20-step heap and tree operation — about 60 million pointer-chasing operations, milliseconds of work on modern hardware versus the trillion-operation O(n²) brute force.

JavaScript implementation

A complete robust Fortune's implementation is several hundred lines (the circle-event geometry and beach-line surgery are fiddly). The structure below shows the event loop and the two handlers — the heart of the algorithm — with the beach line as an ordered arc list and lazy circle-event invalidation.

// Sites: [{x, y}]. Sweep moves down (decreasing y). Returns Voronoi edges.
function fortune(sites) {
  const edges = [];
  const beach = [];                // arcs left-to-right: {site, leftEdge, rightEdge, event}
  // Min-heap on y (highest y = processed first since we sweep downward)
  const queue = new Heap((a, b) => b.y - a.y);

  for (const s of sites) queue.push({ type: 'site', y: s.y, site: s });

  while (!queue.isEmpty()) {
    const ev = queue.pop();
    if (ev.type === 'site') handleSite(ev.site);
    else if (ev.valid) handleCircle(ev);   // lazy deletion: skip stale events
  }
  finishEdges(beach, edges);               // clip unbounded edges to a box
  return edges;

  function handleSite(site) {
    if (beach.length === 0) { beach.push({ site }); return; }
    const i = arcAbove(beach, site.x, site.y);   // O(log n) with a balanced tree
    invalidate(beach[i].event);                  // its circle event is now stale
    const left  = { site: beach[i].site };
    const mid    = { site };
    const right = { site: beach[i].site };
    beach.splice(i, 1, left, mid, right);        // split the arc in three
    // two new edge records start growing from the breakpoints here...
    checkCircle(beach, i, site.y);               // the left copy may now converge
    checkCircle(beach, i + 2, site.y);           // the right copy may now converge
  }

  function handleCircle(ev) {
    const i = ev.arcIndex;
    const vertex = ev.center;                    // equidistant from 3 sites
    edges.push({ from: vertex });                // a Voronoi vertex is born
    invalidate(beach[i - 1].event);
    invalidate(beach[i + 1].event);
    beach.splice(i, 1);                          // the squeezed arc vanishes
    checkCircle(beach, i - 1, ev.y);
    checkCircle(beach, i, ev.y);
  }

  function checkCircle(beach, i, sweepY) {
    if (i <= 0 || i >= beach.length - 1) return; // need a left and right neighbor
    const a = beach[i - 1].site, b = beach[i].site, c = beach[i + 1].site;
    const o = circumcenter(a, b, c);             // null if the three sites are collinear
    if (!o) return;
    const r = Math.hypot(o.x - b.x, o.y - b.y);
    const bottom = o.y - r;                      // circle event fires here
    if (bottom > sweepY) return;                 // already passed — not a future event
    const ev = { type: 'circle', y: bottom, center: o, arcIndex: i, valid: true };
    beach[i].event = ev;
    queue.push(ev);
  }
}
function invalidate(ev) { if (ev) ev.valid = false; }

Two details to flag. First, arcAbove is the step that must be O(log n) — done over a linear array with splice as shown it is really O(n) per event, fine for a demo but quadratic at scale; a production version swaps in a balanced BST. Second, circle events are never removed from the heap; they are flagged valid = false and skipped when popped (lazy deletion), because pulling an arbitrary element out of a binary heap is more trouble than the wasted pops.

Python implementation

The same skeleton in Python, leaning on heapq for the event queue. The circumcenter helper is the one piece of real geometry; everything else is bookkeeping.

import heapq, math

def circumcenter(a, b, c):
    ax, ay = a; bx, by = b; cx, cy = c
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    if abs(d) < 1e-12:                       # collinear: no circle
        return None
    ux = ((ax*ax + ay*ay) * (by - cy)
        + (bx*bx + by*by) * (cy - ay)
        + (cx*cx + cy*cy) * (ay - by)) / d
    uy = ((ax*ax + ay*ay) * (cx - bx)
        + (bx*bx + by*by) * (ax - cx)
        + (cx*cx + cy*cy) * (bx - ax)) / d
    return (ux, uy)

def fortune(sites):
    edges, beach = [], []                    # beach: list of [site, event]
    heap = []                                # entries: (-y, kind, payload)
    counter = 0                              # tie-breaker so tuples never compare payloads
    for s in sites:
        heapq.heappush(heap, (-s[1], counter, 'site', s)); counter += 1

    def check_circle(i, sweep_y):
        nonlocal counter
        if i <= 0 or i >= len(beach) - 1:
            return
        a, b, c = beach[i-1][0], beach[i][0], beach[i+1][0]
        o = circumcenter(a, b, c)
        if o is None:
            return
        r = math.dist(o, b)
        bottom = o[1] - r
        if bottom > sweep_y:                  # event lies above sweep — already gone
            return
        ev = {'center': o, 'arc': i, 'valid': True}
        beach[i][1] = ev
        heapq.heappush(heap, (-bottom, counter, 'circle', ev)); counter += 1

    while heap:
        negy, _, kind, payload = heapq.heappop(heap)
        if kind == 'site':
            site = payload
            if not beach:
                beach.append([site, None]); continue
            i = arc_above(beach, site)        # O(log n) with a BST in production
            if beach[i][1]: beach[i][1]['valid'] = False
            beach[i:i+1] = [[beach[i][0], None], [site, None], [beach[i][0], None]]
            check_circle(i, site[1]); check_circle(i + 2, site[1])
        else:                                 # circle event
            if not payload['valid']:          # stale: skip (lazy deletion)
                continue
            i = payload['arc']
            edges.append(payload['center'])   # Voronoi vertex
            for n in (i-1, i+1):
                if 0 <= n < len(beach) and beach[n][1]:
                    beach[n][1]['valid'] = False
            beach.pop(i)
            check_circle(i - 1, -negy); check_circle(i, -negy)
    return edges

Note the counter tie-breaker pushed into every heap tuple. Without it, two events at the same y force Python to compare the next tuple field — eventually the dict payload — and raise TypeError: '<' not supported between dicts. A monotonic counter makes every tuple unique before the payload is ever reached.

Variants worth knowing

Weighted (power / Laguerre) Voronoi. Give each site a weight and measure distance as |p − s|² − w. Cells become power cells, and the structure models things like disks of different radii (molecular surfaces, influence zones). The beach line generalizes; the bisectors are still straight lines.

Manhattan / L∞ Voronoi. Swap Euclidean distance for L1 or L∞. Cell boundaries become axis-aligned and 45° segments instead of bisectors, which is exactly what VLSI routing and grid pathfinding want.

Centroidal Voronoi tessellation (CVT). Iterate Lloyd's algorithm — move each site to the centroid of its cell, rebuild, repeat. The diagram relaxes into evenly spaced, blue-noise cells used for stippling, remeshing, and de-clumping samples.

Additively / power-weighted diagrams in 3D. The sweepline does not survive in 3D; instead lift sites to a paraboloid, take the lower convex hull (that is the Delaunay tetrahedralization), and dualize. This is what Qhull and CGAL do for any dimension above two.

Approximate GPU Voronoi (jump flooding). Seed each site's pixel, then flood with halving step sizes (n/2, n/4, …, 1). After ⌈log₂ n⌉ passes every pixel holds its nearest seed. Pixel-accurate, embarrassingly parallel, and the standard real-time technique for soft shadows, distance fields, and live cell rendering.

Common bugs and edge cases

• Collinear or co-circular sites. Three collinear sites have no circumcenter (the bisectors are parallel); four co-circular sites create a degenerate vertex. Robust predicates (exact arithmetic or careful epsilon handling) are what separate a toy from a library.
• Two sites with identical y. The first site event has no arc to split. Handle the empty/first-arc case explicitly, and break y-ties with x so the order is total.
• Scheduling circle events that point upward. Only schedule a circle event whose bottom lies at or below the current sweep — a "left turn" triple converges upward and must be rejected, or you create vertices the sweep already passed.
• Forgetting lazy invalidation. When an arc is removed, its neighbors' previously scheduled circle events are now stale. If you process a stale event you delete the wrong arc and corrupt the beach line. Mark and skip.
• Leaving edges unbounded. Hull cells have rays to infinity. Clip every dangling breakpoint to a bounding box at the end, or rendering and point-location both break.
• Using a linked-list beach line at scale. The arcAbove locate step is O(n) on a list, silently turning the whole algorithm O(n²). Use a balanced BST keyed by breakpoint x for the advertised O(n log n).