Monogon's blog

By Monogon, history, 3 years ago, In English

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.

Implementation

Example Problems

Panda Preserve solution
Caravans solution
Good Manners solution

Resources

  • Vote: I like it
  • +499
  • Vote: I do not like it

»
3 years ago, # |
  Vote: I like it +77 Vote: I do not like it

Many thanks to tfg for helping me understand the background on Voronoi diagrams better, and revising my code!

»
3 years ago, # |
  Vote: I like it +10 Vote: I do not like it

Please add the "[Tutorial]" prefix to the blog's name.

»
3 years ago, # |
  Vote: I like it +74 Vote: I do not like it

Fun fact: This started in discord with the question "how to solve nearest point fast without point location + voronoi? (not using kd-tree)", the topic shifted to voronoi/delaunay then I realized "Well, you can use fortune's algo with some mutant set probably (to make the implementation easier/shorter). Someone needs to code it tho. Mutant set as something like that CHT set hack".

Thank you for coding it! I'd take way longer adjusting that set hack.

  • »
    »
    3 years ago, # ^ |
      Vote: I like it +1 Vote: I do not like it

    You should really start a youtube channel or twitch or blogs of yours for non-trivial ds and algorithms,you have a great knowledge and understanding of ds algo
    P.S. great fan of yours.

»
3 years ago, # |
  Vote: I like it +140 Vote: I do not like it

I'm going to pretend like I understood this

»
3 years ago, # |
  Vote: I like it +128 Vote: I do not like it

Other people: Write useful blogs and gain contribution

Monogon: First becomes top contributor, then writes useful blogs.

  • »
    »
    3 years ago, # ^ |
    Rev. 3   Vote: I like it +2 Vote: I do not like it

    I mean, he has already written previous useful blogs on 3D Convex Hull and also was authors of one Global Rounds(12) this year and several other Div 1 and Div 2 rounds, which were equally useful. So he deserves all the contribution or am I missing something.

    • »
      »
      »
      3 years ago, # ^ |
        Vote: I like it +39 Vote: I do not like it

      Yeah I never disagreed. All my respects to Monogon.

  • »
    »
    3 years ago, # ^ |
      Vote: I like it +10 Vote: I do not like it

    Meanwhile mesanu writes meme blog and gains contribution.

»
3 years ago, # |
Rev. 2   Vote: I like it +10 Vote: I do not like it

Wow! That implementation looks really nice and clean (and short as well). I have tried to implement this once as well using the "hacked set" but I got too anxious after about 4 hours and gave up.

I'm just wondering, what are the problems that you encountered with points with same x coordinate? Why wouldn't it work? Also, how did you test the correctness of implementation vs degenerate cases (and especially precision errors)?

  • »
    »
    3 years ago, # ^ |
      Vote: I like it 0 Vote: I do not like it

    When points have same x coordinate, the gety function divides by 0. I could add a simple check and return med.y in that case, but I still see issues. When I test on a grid of points, the leftmost points don't look right, and I'm still not entirely sure why.

    I tested the implementation by generating large cases and degenerate cases, and seeing the output with graphics. I also got the problem Caravans accepted. There is definitely room for improvement and further testing.

  • »
    »
    3 years ago, # ^ |
    Rev. 4   Vote: I like it 0 Vote: I do not like it

    To add on what Monogon said, I think the case that broke for without rotation is actually when we add a point and some beach line intersection is directly to the left of the point, so that means that there's a point and a point in its voronoi cell that have the same x coordinate. That makes lower_bound with doubles not work as well. Good thing is that as long as the points are not retarded (as in we don't input floating number coordinates, so min distance is at least 1) rotating by a random angle makes that case not exist with probability near 1 because there are finite bad angles and theoretically infinite angles to be chosen (not really because floating point is finite but close enough).

    You could also manipulate the getY formula by passing things around and squaring both sides of the inequality to get an exact algorithm but you'd have to handle some other things like circle events not being exactly in integer points and the things after squaring that inequality having values up to O(MAX_COORDINATE^4).

    Overall I trust in the rotation more than I don't trust this and you could theoretically run incremental delaunay (the O(flips + N) algorithm that you can run from a starting triangulation) from the triangulation you get from this because in case of any misbehavior due to precision, as long as it's not in the lower_bound thing, the result should be a few flips away from a delaunay triangulation (aka close enough to delaunay if we don't use int128 or something with absurd precision to check for flips which is, if you ask me, good enough).

    • »
      »
      »
      3 years ago, # ^ |
        Vote: I like it 0 Vote: I do not like it

      Ok, so as I understand, it’s not only that points have the same x coordinate, but more that vertices or edges of voronoi cells may have the same y coordinate as the points. So I guess this means that we should rotate in all cases?

      • »
        »
        »
        »
        3 years ago, # ^ |
          Vote: I like it 0 Vote: I do not like it

        Yes, this is also an issue. It should be possible to fix by using exact gety comparisons, but rotating all the input points seems like an easier fix to me.

»
3 years ago, # |
  Vote: I like it 0 Vote: I do not like it

Great Post!

P.S.) The implementation code isn't in the code block.

  • »
    »
    3 years ago, # ^ |
      Vote: I like it +28 Vote: I do not like it

    It was in the spoiler box before. It's broken now because codeforces is feeling 2020 and things are randomly breaking.