Voronoi diagram and Delaunay triangulation are dual representations of a set of points to each other. Due to their wide application in science and technology, Voronoi diagram and Delaunay triangulation play important roles in the field of Computational Geometry. In this post, I am going to introduce an implementation of an algorithm to derive both Voronoi Diagram and Delaunay Triangulation of a set of points in the plane.
Voronoi Diagram and Delaunay Triangulation
The Voronoi diagram of a set of points, also known as Thiessen polygons, is a partitioning of a plane into regions by a set of continuous polygons consisting of perpendicular bisectors of the connecting lines of two adjacent points. These regions are called Voronoi cells. And for each point in the set, there is a corresponding Voronoi cell consists of all points closer to that point than to any other.
The Delaunay triangulation of a set of points is dual to its Voronoi diagram. It is a collection of connected but nonoverlapping triangles, and the outer circumcircle of these triangles does not contain any other points in this set.
Design of the Algorithm
There are a lot of ways to generate a Voronoi diagram from a set of points in the plane. In this implementation, I obtained the Voronoi diagram from generating its dual, the Delaunay triangulation. Generally speaking, for the set of $n$ points $P=\{p_1,p_2,\ldots,p_n\}$ in $\mathbb{R}^2$, the algorithm goes in this way: the Delaunay triangulation of the set of points is firstly generated, then we calculate the center of the circumcircle of each triangle and finally we connect these centers with straight lines and form the polygon mesh generated from the vertices of the triangles.
Design of the Data Structure
I implemented this algorithm in object oriented language, so the design of data structure is in the form of class.
Point:


Voronoi edge:


Delaunay edge:


Triangle:


On top of that, I defined a point list and a triangle list under the Collections class to store all the points in the point set $P$ and all the triangles during the procedure of triangulation as global variables:


Thus, in the Point class, DelaunayEdge class and Triangle class, I can use an integer to represent the index of a certain point or triangle from the global lists, and retrieve it directly.
Build up the Vonoroi Diagram
The algorithm of building up a Vonoroi diagram goes in this way:
Begin
Step 1: Obtain the Delaunay triangulation by generating the list of Delaunay edges.
Step 2: Traverse all the Delaunay edges.
Step 3: For each Delaunay edge, traverse the two triangle lists stored with the start point and the end point, and find the two same triangles in the two lists which are the adjacent triangles of this Delaunay edge.
Step 4: Construct a Voronoi edge by connecting the two centers of the circumcles of the two adjacent triangles and add it to the Voronoi edge list.
End
Conduct the Delaunay Triangulation
I planned to apply the Bowyer–Watson algorithm for computing the Delaunay triangulation. A most naïve Bowyer–Watson algorithm goes like this:
Begin
Step 1: Construct a ‘super’ triangle that covers all the points from the point set, add it to the Delaunay triangle list.
Step 2: Insert points from the point set $P=\{p_1,p_2,\ldots,p_n\}$ to the ‘super’ triangle one by one.
Step 3: For each point $p_i$ inserted, traverse the Delaunay triangle list to find all the triangles whose circumcircles cover this point $p_i$ as invalid triangles, delete these triangles from the Delaunay triangle list and delete all the common edges of these triangles, and leave a starshaped polygonal hole.
Step 4: Connect the point $p_i$ to all the vertices of this starshaped polygon, and add the newly formed triangles to the Delaunay triangle list.
Step 5: After all the points are inserted, obtain the Delaunay edge list from the Delaunay triangle list, and delete the edges from the Delaunay edge list that contain a vertex of the original ‘super’ triangle.
End
The following pictures can better illustrate the key steps of this algorithm. As shown in Figure 1, when a new point is inserted, all the triangles whose circumcircles contain this point will be found. The common edges of these triangles, which are highlighted in yellow, will be deleted, leaving the starshaped boundary in red.
And then, the inserted point will be connected to all the vertices of the starshaped polygon, as shown in Figure 2, the new Delaunay triangles will be formed.
However, the aforementioned naïve manner of Delaunay triangulation is clearly an $O(n^2)$ time algorithm that is incapable of handling a massive amount of points.
For improving its efficiency, we can first sort the set of points by xcoordinate, and use an open list and a closed list to store all the Delaunay triangles. In each time a point is inserted, all the triangles with circumcircles to the left of the inserting point are put into the closed list and removed from the open list. All the new triangles generated in an insertion are put into the open list. So in each time of insertion, we just have to traverse the open list to find the invalid triangles, the length of the sequential search of triangles is much reduced.
Besides, in order to derive the Voronoi diagram, when a new point is inserted, all the newly formed triangles that are incident on the point are put into the triangle list that is stored with the point.
Thus, the optimized algorithm goes as follows:
Begin
Step 1: Sort the points in the point set $P=\{p_1,p_2,\ldots,p_n\}$ by xcoordinate.
Step 2: Construct a ‘super’ triangle that covers all the points from the point set, add it to the open list.
Step 3: Insert points from $P$ to the ‘super’ triangle one by one in ascending order of xcoordinates.
Step 4: For each point $p_i$ inserted, traverse the open list to: 1) find all the triangles with circumcircles lying to the left of the point $p_i$, delete these triangles from the open list and add them to the closed list; 2) find all the triangles with circumcircles covering this point $p_i$ as invalid triangles, delete these triangles from the open list and the triangle list stored with $p_i$, delete all the common edges of these triangles, leaving a starshaped polygonal hole.
Step 5: Connect the point $p_i$ to all the vertices of this starshaped polygon, and add the newly formed triangles to the open list and the triangle list stored with $p_i$.
Step 6: After all the points are inserted, merge the open list and the closed list into the Delaunay triangle list, obtain the Delaunay edge list from the Delaunay triangle list, and delete the edges from the Delaunay edge list that contain a vertex of the original ‘super’ triangle.
End
Time Complexity of the Algorithm
The optimized algorithm for Delaunay triangulation takes $O(n\log n)$ time. As Delaunay triangulation is a planar graph, the number of triangles incident on one point is constant, so the procedure of finding adjacent triangles takes constant time and the time of generating Voronoi diagram is $O(n)$. Therefore, the total running time of this algorithm is $O(n\log n)$.
Implementation of the Algorithm
This algorithm is implemented in C#. Figure 3 shows the result of the implementation of this algorithm, it is capable of handling massive input of 10000 points.
The GitHub repository of this implementation is IsaacGuan/VoronoiDelaunay.