Fortune's Sweep and the Beach Line of Parabolic Arcs

The Voronoi diagram of a set of sites partitions the plane by nearest site: every location belongs to the cell of the site closest to it. Cells are convex polygons, because each is the intersection of the half-planes defined by the perpendicular bisectors between its site and every other site. For n sites the diagram has at most 2n − 5 vertices and 3n − 6 edges, so its size is linear, which means a construction that is quadratic in n is doing far more work than the output requires.

The straightforward construction computes each cell by itself, as the intersection of the n − 1 half-planes cut by the perpendicular bisectors to the other sites; a single cell costs O(n log n), so the whole diagram costs Θ(n² log n) and the result is a structure of linear size. Shamos and Hoey gave an O(n log n) divide-and-conquer construction in the mid-1970s, matching the sorting lower bound, and Fortune’s sweep-line algorithm, published in 1987, achieves the same bound with a structure that is simple enough to implement directly. It also differs from the segment sweeps in its status structure: not a chain of straight segments but a chain of parabolic arcs, which accounts for most of the added implementation difficulty.

Why a straight sweep line does not work

In the Bentley–Ottmann sweep the status holds segments that the line actually crosses, ordered by their intersection with it. For the Voronoi diagram, the objects being built are cells that extend to infinity, and the sweep line crosses each cell over a range of y rather than at a point. Ordering cells along the line is possible, but a cell’s extent along the line changes as the line moves, and the changes are not driven by a finite set of events.

Fortune’s insight was to sweep a line downward and maintain not the diagram but the part of it that is already determined. A site above the line has a region of points that are closer to it than to any site below the line, and that region is bounded below by a parabola, the locus of points equidistant from the site and from the sweep line itself. The sweep line stands in for the nearest site below, so the parabola is exactly the boundary between the site’s territory and the not-yet-processed territory. The union of these parabolas, one per processed site, is the beach line, and it is a monotone curve in x.

The beach line is where the diagram’s edges are born. Each arc of the beach line belongs to one site, and where two arcs meet, the meeting point is equidistant from both sites and from the sweep line, so it lies on the Voronoi edge between them. As the line advances, the meeting points trace out that edge.

import math

def breakpoint_x(a, b, ly):
    """x where the parabolas of sites a and b meet, given directrix y = ly."""
    d1, d2 = 2.0 * (a[1] - ly), 2.0 * (b[1] - ly)
    # Parabola of site s: y = ((x - sx)^2 + sy^2 - ly^2) / (2*(sy - ly)).
    # Equating the two and clearing denominators gives A*x^2 + B*x + C = 0.
    A = d2 - d1
    B = -2.0 * d2 * a[0] + 2.0 * d1 * b[0]
    C = d2 * (a[0] ** 2 + a[1] ** 2 - ly ** 2) - d1 * (b[0] ** 2 + b[1] ** 2 - ly ** 2)
    if A == 0.0:                       # sites at equal height: perpendicular bisector
        return -C / B
    disc = B * B - 4.0 * A * C
    r = math.sqrt(max(disc, 0.0))
    x1, x2 = (-B - r) / (2.0 * A), (-B + r) / (2.0 * A)
    lo, hi = (a[0], b[0]) if a[0] <= b[0] else (b[0], a[0])
    return x1 if lo <= x1 <= hi else x2   # the root between the two sites

The function is written for exposition: it returns a coordinate, which is exactly the thing an implementation should avoid depending on. A production version compares a query’s x against the breakpoint by evaluating both sides of the equality in a form that keeps the comparison exact, rather than by computing the breakpoint and comparing floats. The formula above also shows why the sweep’s inner loop is more delicate than the outer event loop: the breakpoints move continuously, so a cached value is valid only until the sweep advances.

The binary search over the beach line compares a query’s x against a sequence of breakpoints, and each comparison is a rational expression in site coordinates and the current sweep position. Cross-multiplying turns the comparison into a sign test on an integer or exactly representable expression, which is the only way the search stays consistent as the sweep advances and the breakpoints drift.

Two kinds of event

A site event occurs when the sweep line reaches a new site. The new site is below every processed site, so its parabola is flatter than the arcs above it and must be inserted into the beach line. It intersects exactly one existing arc, splitting it into two pieces; the new site’s arc appears between them. Two new breakpoints appear, and each starts tracing a new Voronoi edge downward.

A circle event occurs when three consecutive arcs become tangent at a single point: the middle arc shrinks to zero width and disappears. At that moment, the two breakpoints that bounded it meet, which means two Voronoi edges terminate at a Voronoi vertex. The vertex is the centre of the circle through the three sites, and the event fires exactly when the sweep line reaches the lowest point of that circle. After the event, the arcs on either side of the vanished arc become adjacent and a new breakpoint starts a new edge.

The status is the sequence of arcs, stored in a balanced tree ordered left to right along the beach line, with each arc carrying pointers to its two breakpoints. Finding the arc that a new site splits is a search for the x-position of the site, which is O(log n) if the tree supports it; the search compares the site’s x against the parabolas’ intersections. That comparison is itself a predicate on the site coordinates and is subject to the same numerical hazards as any other, as discussed in exact geometric predicates.

Scheduling circle events, and the invalidation problem

Each triple of consecutive arcs potentially defines a circle event, and the algorithm schedules the earliest one by pushing the circle’s lowest point into the event queue. The catch is that triples change constantly: inserting a site destroys one triple and creates two; removing an arc destroys its own triple and creates one new triple from its neighbours’ new adjacency.

A scheduled event is valid only if, at the time it is popped, the three arcs are still consecutive and the middle one still exists. The standard implementation attaches to each arc a pointer to the circle event it is the middle of, and invalidates that event whenever the arc is removed or its neighbours change. When an event is popped, the check is a pointer comparison: if the arc’s stored event pointer does not match the popped event, the event is stale and is discarded.

This is the same stale-event pattern as in Bentley–Ottmann, and it is the reason both algorithms require a status structure with stable references to its elements. An array with shifting indices cannot support it, because the pointer stored with an arc would become meaningless after an insertion elsewhere in the line.

The degenerate cases that dominate the code

Sites with equal y produce parabolas of equal width, so the breakpoint between two such arcs sits at a fixed x and the edge they trace is vertical; the quadratic that locates the breakpoint degenerates to a linear equation, which an implementation has to special-case. Three sites on a common line produce a circle of infinite radius, so no circle event ever fires and the arcs never vanish; the corresponding Voronoi edges are parallel and meet at infinity. Four cocircular sites produce simultaneous circle events, and processing them in the wrong order makes an arc vanish twice or leaves the beach line inconsistent.

The usual treatment mirrors the one for Delaunay triangulation: symbolic perturbation, or an explicit sort key that makes the input a total order with no ties, or a pre-pass that perturbs coordinates by a tiny deterministic amount. All three convert degeneracies into general position, and all three make the output differ from the unperturbed diagram by an amount that is small in coordinates but can be large in combinatorics — a Voronoi vertex that should be a single degree-4 point becomes two degree-3 points joined by a short edge.

Cost, and what the diagram is for

The event queue holds O(n) site events and O(n) circle events, since each arc can be the middle of at most one live event and there are O(n) arcs. Each event costs O(log n) for the queue and O(log n) for the status update, so the total is O(n log n) time and O(n) space. The diagram itself is linear in size, so the space bound is optimal.

The primary use is nearest-neighbour queries. With the diagram built, a query is a point-location problem: find which cell contains the point, and its site is the nearest. That is the same problem solved by the trapezoidal map point location structure, and Fortune’s algorithm can be followed by a conversion to a point-location structure, or the diagram can be stored in a bounding-box hierarchy. The k-d tree range search is the alternative for the same queries when the sites are static and the diagram’s construction cost is not amortised: it answers nearest-neighbour queries without building the diagram, exactly in its standard form and within a bounded error ratio in its approximate variant.

A second use is the dual relationship with Delaunay: the Voronoi vertices correspond to Delaunay triangles, the Voronoi edges to Delaunay edges, and the cells to vertices. Building one gives the other in linear time, which is why libraries that need both build whichever is easier and derive the other, rather than running two constructions. The diagram also underlies the largest-empty-circle problem, medial axis approximations, and the natural-neighbour interpolation scheme, all of which read the diagram’s combinatorial structure rather than its coordinates.