What's a Voronoi Diagram?

Given a set $$$S$$$ of $$$n$$$ points in the 2D plane, the Voronoi diagram is a partition of the plane into regions. The region associated with a point $$$p\in S$$$ contains precisely the points $$$q$$$ such that $$$q$$$ is closer to $$$p$$$ than any other point in $$$S$$$. In other words, a point $$$q$$$ belongs to the region of its nearest neighbor.

The Voronoi diagram can be represented by the planar graph of boundaries between regions. A vertex in this graph is where three or more segments meet (or a point at infinity), and an edge is a segment connecting two of these vertices. The way the regions were defined, we see that a segment contains the points equidistant to two points in $$$S$$$, so it is part of their perpendicular bisector. Similarly, a Voronoi vertex is equidistant to three or more vertices in $$$S$$$, so it is their circumcenter. Recall that the circumcenter is the intersection of perpendicular bisectors.

What's a Delaunay Triangulation?

The Delaunay triangulation is the dual graph of the Voronoi diagram. That is, a Voronoi vertex is a Delaunay face and a Delaunay vertex is a Voronoi region. A Delaunay edge connects two points if and only if their corresponding Voronoi regions share a border.

In the following image, the set $$$S$$$ corresponds to the black points. The red vertices and edges are the Voronoi diagram. The black points and edges are the Delaunay triangulation.

The Delaunay triangulation only involves edges between existing points, while the Voronoi diagram creates new vertices and needs a way to represent edges going infinitely in one direction. For this reason, it is often more convenient to compute the Delaunay triangulation and only convert to the Voronoi diagram if needed.

Why are these useful?

The definition of the Voronoi diagram immediately shows signs of applications. Given a set $$$S$$$ of $$$n$$$ points and $$$m$$$ query points $$$p_1,\ldots, p_m$$$, we can answer for each query point, its nearest neighbor in $$$S$$$. This can be done in $$$O((n + q)\log(n+q))$$$ offline by sweeping the Voronoi diagram and query points. Or it can be done online with persistent data structures.

Other than this, the Delaunay triangulation has many amazing properties of its own that can make it extremely useful for solving problems.

• For each Delaunay triangle, its circumcircle does not strictly contain any points in $$$S$$$. (In fact, you can also consider this the defining property of Delaunay triangulation)
• The number of Delaunay edges is at most $$$3n-6$$$, so there is hope for an efficient construction.
• Each point $$$p\in S$$$ is adjacent to its nearest neighbor with a Delaunay edge.
• The Delaunay triangulation maximizes the minimum angle in the triangles among all possible triangulations.
• The Euclidean minimum spanning tree is a subset of Delaunay edges.

Degenerate Cases

Because this is geometry, there are some degenerate cases I have hidden from you up to this point. The first case is when all the points in $$$S$$$ are contained in a line. In this case, the Voronoi diagram has all edges connecting two points at infinity, and the Voronoi diagram is disconnected (in all other cases it's connected). Also in this case, the dual graph is no longer a triangulation. We may still decide to represent it by the path of edges connecting the points.

Another nasty case is co-circular points. Here, there are Voronoi vertices where more than three segments intersect. In this case, the dual graph has a face with more than three sides. However, we want it to be a triangulation, so we can also say that any triangulation of such faces is valid, and we will allow non-unique Delaunay triangulations.

Algorithms

Many algorithms exist to compute the Delaunay triangulation and/or Voronoi diagram of a set of points. Each have their advantages and disadvantages.

Flip Algorithms

In this approach, you start with an arbitrary triangulation of the points. Then you repeatedly fix adjacent triangle pairs until the triangulation has the Delaunay property. The main advantage of this method is its simplicity: an arbitrary triangulation can be constructed with an incremental convex hull algorithm, and the disadvantage is that it could require $$$\Omega(n^2)$$$ flips in the worst case.

Incremental Algorithms

In this approach, you always maintain the Delaunay triangulation of a subset of the points and insert the points one at a time. When a point is inserted, you should delete all triangles whose circumcircle contains the newly added point. After deleting this connected set of triangles, you should repair the face it exposes. For each edge in this exposed face, you should add a triangle with that edge and the new point. This gives the straightforward $$$O(n^2)$$$ Bowyer-Watson algorithm.

If we insert the points in a random order, we can significantly reduce the number of triangle insertions and deletions. But to get expected $$$O(n\log n)$$$ runtime, we would need to avoid processing all the unaffected triangles. And this makes the implementation far more challenging. We need to maintain a search structure of the triangles in order to locate which one contains the new point.

The advantage of this method is that it achieves expected $$$O(n\log n)$$$ runtime with only integer computations, and the construction is online. The disadvantage is that it has a bad constant, isn't deterministic, and the triangle search structure can make for a tricky implementation.

Divide and Conquer

The divide and conquer algorithm first sorts the points by $$$x$$$-coordinate. Then it splits the points into two sets, recursively computes their Delaunay triangulations, and merges them together. As long as you can perform the merge in $$$O(n)$$$ time, the overall algorithm is $$$O(n\log n)$$$. When merging the triangulations, it is actually quite hard to keep track of the important adjacency information. The most popular strategy is perhaps the Quad-Edge data structure.

The advantage is that this gives a deterministic $$$O(n\log n)$$$ algorithm with only integer computations. The disadvantage is the complexity of the quad-edge data structure and the costly memory usage that comes with it.

3D Convex Hull Reduction

An interesting fact is that the problem of Delaunay triangulation can be reduced to 3D convex hull. If for each point $$$(x,y)$$$ we also give it the $$$z$$$-coordinate $$$x^2+y^2$$$ and compute the 3D convex hull, then the downward-facing convex hull faces correspond exactly to the Delaunay triangles. The advantage is that if you happen to have a good 3D convex hull implementation, you get Delaunay triangulation for free. The disadvantages are that $$$O(n\log n)$$$ 3D convex hull is a whole beast of its own, and the large $$$z$$$-coordinates will exacerbate precision and overflow errors in the 3D convex hull implementation. The main strategies for 3D convex hull are incremental and divide-and-conquer algorithms, which are only more complicated than their counterparts in the special case of Delaunay triangulation.

My previous blog on 3D convex hull gave an $$$O(n\log n)$$$ randomized incremental implementation, but it turns out to have quite a bad constant for practical use.

Fortune's Algorithm

In this blog, I will focus on Fortune's Algorithm. It takes a sweep-line approach to the problem, and achieves a deterministic $$$O(n\log n)$$$ time complexity. The main disadvantage is that computations are done with floating points instead of integers, opening up the door to precision errors. However, keep in mind that the algorithms using only integers have a risk of integer overflow even on reasonably small constraints.

The main challenge in designing a sweep-line algorithm for Voronoi diagrams is that a point after the sweep-line can affect the diagram before the sweep-line. That is, if the sweep-line denotes the discovery time of a point, we cannot know the entire Voronoi diagram before the sweep-line. The way Fortune's algorithm overcomes this challenge is by maintaining the Voronoi diagram only for the area we can be certain about, regardless of what points lie after the sweep-line.

Now, which points can be affected by what lies past the sweep-line? Well, such a point must be closer to the sweep-line than to any point in $$$S$$$ to the left of the sweep-line. Recall that the set of points equidistant to a point and a line is given by a parabola, where the point is called the focus and the line is called the directrix. If the focus is a point $$$p\in S$$$ to the left of the sweep-line and the directrix is the sweep-line, then everything to the left of this parabola we can be certain of its region. Now, if we make such a parabola for every point $$$p\in S$$$ to the left of the sweep-line, everything to the right of all parabolas is precisely the uncertain region. The boundary of the certainty and uncertainty regions is called the beach line, consisting of parabolic arcs.

We are ready to understand Fortune's algorithm abstractly. We scan a sweep-line from left to right, maintaining the combinatorial structure of the beach line. The parabolas are only for visual aid, and the ordered list of foci is what is actually stored. As points get added and parabolic arcs shrink to length $$$0$$$, these events tell us precisely the structure of the Voronoi diagram and Delaunay triangulation.

There are two types of events which change the beach line structure: point events and vertex events.

Point Events

In a point event, the sweep-line discovers a point $$$p\in S$$$ and inserts its parabolic arc into the beach line. In order to do this efficiently, we need to search for where it belongs in the beach line. Then it will split the arc it sees into two. An example of a point event is shown below (here the sweep-line goes from top to bottom instead).

Vertex Events

In a vertex event, a parabolic arc disappears from the beach line, corresponding to a Voronoi vertex. To handle a vertex event, we simply remove the point whose arc length shrinks to length $$$0$$$. Such an event is generated dynamically by a triple of points adjacent on the beach line. Every time we modify the beach line, we update the relevant points for their vertex events. To compute the time a vertex event should be handled, we take the circumcircle of the three points generating the event. The rightmost point of the circumcircle is the time we should process this event.

Summary

In summary, the algorithm maintains a priority queue of events. Initially, we push all point events into the priority queue. When processing a point event, we find the correct arc in the beach line, split it up, add the new point, and add new vertex events if needed. When processing a vertex event, we add the generating triple as a Delaunay triangle, remove the arc that should disappear, and add new vertex events if needed.

Fortune's Algorithm Implementation

My code doesn't work when two points have the same $$$x$$$ coordinate. This is handled simply by rotating all input points by $$$1$$$ radian and praying to the geometry gods.

#include <bits/stdc++.h>
 
#define ll long long
#define sz(x) ((int) (x).size())
#define all(x) (x).begin(), (x).end()
#define vi vector<int>
#define pii pair<int, int>
#define rep(i, a, b) for(int i = (a); i < (b); i++)
using namespace std;
template<typename T>
using minpq = priority_queue<T, vector<T>, greater<T>>;
 
using ftype = long double;
const ftype EPS = 1e-12, INF = 1e100;
 
struct pt {
    ftype x, y;
    pt(ftype x = 0, ftype y = 0) : x(x), y(y) {}

    // vector addition, subtraction, scalar multiplication
    pt operator+(const pt &o) const {
        return pt(x + o.x, y + o.y);
    }
    pt operator-(const pt &o) const {
        return pt(x - o.x, y - o.y);
    }
    pt operator*(const ftype &f) const {
        return pt(x * f, y * f);
    }

    // rotate 90 degrees counter-clockwise
    pt rot() const {
        return pt(-y, x);
    }

    // dot and cross products
    ftype dot(const pt &o) const {
        return x * o.x + y * o.y;
    }
    ftype cross(const pt &o) const {
        return x * o.y - y * o.x;
    }

    // length
    ftype len() const {
        return hypotl(x, y);
    }

    // compare points lexicographically
    bool operator<(const pt &o) const {
        return make_pair(x, y) < make_pair(o.x, o.y);
    }
};

// check if two vectors are collinear. It might make sense to use a
// different EPS here, especially if points have integer coordinates
bool collinear(pt a, pt b) {
    return abs(a.cross(b)) < EPS;
}


// intersection point of lines ab and cd. Precondition is that they aren't collinear
pt lineline(pt a, pt b, pt c, pt d) {
    return a + (b - a) * ((c - a).cross(d - c) / (b - a).cross(d - c));
}

// circumcircle of points a, b, c. Precondition is that abc is a non-degenerate triangle.
pt circumcenter(pt a, pt b, pt c) {
    b = (a + b) * 0.5;
    c = (a + c) * 0.5;
    return lineline(b, b + (b - a).rot(), c, c + (c - a).rot());
}

// x coordinate of sweep-line
ftype sweepx;

// an arc on the beacah line is given implicitly by the focus p,
// the focus q of the following arc, and the position of the sweep-line.
struct arc {
    mutable pt p, q;
    mutable int id = 0, i;
    arc(pt p, pt q, int i) : p(p), q(q), i(i) {}

    // get y coordinate of intersection with following arc.
    // don't question my magic formulas
    ftype gety(ftype x) const {
        if(q.y == INF) return INF;
        x += EPS;
        pt med = (p + q) * 0.5;
        pt dir = (p - med).rot();
        ftype D = (x - p.x) * (x - q.x);
        return med.y + ((med.x - x) * dir.x + sqrtl(D) * dir.len()) / dir.y;
    }
    bool operator<(const ftype &y) const {
        return gety(sweepx) < y;
    }
    bool operator<(const arc &o) const {
        return gety(sweepx) < o.gety(sweepx);
    }
};

// the beach line will be stored as a multiset of arc objects
using beach = multiset<arc, less<>>;

// an event is given by
//     x: the time of the event
//     id: If >= 0, it's a point event for index id.
//         If < 0, it's an ID for a vertex event
//     it: if a vertex event, the iterator for the arc to be deleted
struct event {
    ftype x;
    int id;
    beach::iterator it;
    event(ftype x, int id, beach::iterator it) : x(x), id(id), it(it) {}
    bool operator<(const event &e) const {
        return x > e.x;
    }
};

struct fortune {
    beach line; // self explanatory
    vector<pair<pt, int>> v; // (point, original index)
    priority_queue<event> Q; // priority queue of point and vertex events
    vector<pii> edges; // delaunay edges
    vector<bool> valid; // valid[-id] == true if the vertex event with corresponding id is valid
    int n, ti; // number of points, next available vertex ID
    fortune(vector<pt> p) {
        n = sz(p);
        v.resize(n);
        rep(i, 0, n) v[i] = {p[i], i};
        sort(all(v)); // sort points by coordinate, remember original indices for the delaunay edges
    }
    // update the remove event for the arc at position it
    void upd(beach::iterator it) {
        if(it->i == -1) return; // doesn't correspond to a real point
        valid[-it->id] = false; // mark existing remove event as invalid
        auto a = prev(it);
        if(collinear(it->q - it->p, a->p - it->p)) return; // doesn't generate a vertex event
        it->id = --ti; // new vertex event ID
        valid.push_back(true); // label this ID true
        pt c = circumcenter(it->p, it->q, a->p);
        ftype x = c.x + (c - it->p).len();
        // event is generated at time x.
        // make sure it passes the sweep-line, and that the arc truly shrinks to 0
        if(x > sweepx - EPS && a->gety(x) + EPS > it->gety(x)) {
            Q.push(event(x, it->id, it));
        }
    }
    // add Delaunay edge
    void add_edge(int i, int j) {
        if(i == -1 || j == -1) return;
        edges.push_back({v[i].second, v[j].second});
    }
    // handle a point event
    void add(int i) {
        pt p = v[i].first;
        // find arc to split
        auto c = line.lower_bound(p.y);
        // insert new arcs. passing the following iterator gives a slight speed-up
        auto b = line.insert(c, arc(p, c->p, i));
        auto a = line.insert(b, arc(c->p, p, c->i));
        add_edge(i, c->i);
        upd(a); upd(b); upd(c);
    }
    // handle a vertex event
    void remove(beach::iterator it) {
        auto a = prev(it);
        auto b = next(it);
        line.erase(it);
        a->q = b->p;
        add_edge(a->i, b->i);
        upd(a); upd(b);
    }
    // X is a value exceeding all coordinates
    void solve(ftype X = 1e9) {
        // insert two points that will always be in the beach line,
        // to avoid handling edge cases of an arc being first or last
        X *= 3;
        line.insert(arc(pt(-X, -X), pt(-X, X), -1));
        line.insert(arc(pt(-X, X), pt(INF, INF), -1));
        // create all point events
        rep(i, 0, n) {
            Q.push(event(v[i].first.x, i, line.end()));
        }
        ti = 0;
        valid.assign(1, false);
        while(!Q.empty()) {
            event e = Q.top(); Q.pop();
            sweepx = e.x;
            if(e.id >= 0) {
                add(e.id);
            }else if(valid[-e.id]) {
                remove(e.it);
            }
        }
    }
};

Example Problems

• Panda Preserve (ICPC WF 2018)
• Caravans (acmsguru)
• Good Manners (Timus)
• Empire Strikes Back (Timus)

Resources

• Fortune's Algorithm Notes
• Delaunay Triangulation (Wikipedia)
• Voronoi Diagram (Wikipedia)
• CP Algorithms, Divide and Conquer Algorithm