## Voronoi diagrams on the sphere

Lately I wanted to be able to compute the Voronoi diagram of a family of points on the unit sphere in an efficient way. I started by looking all over the place for implemented results and I found a few: a Python implementation, a Matlab implementation or an online interactive Java implementation. Each of these have their inconveniences: I’m not familiar with Python or Java and the Matlab code is not vectorized, which means it is really slow. So the challenge was on. Luckily, the idea for a fast, efficient algorithm was found here. I present below the main idea, along with some justifications. For the sake of simplicity, suppose we work on the unit sphere.

In order to compute the Voronoi diagram of a set of points on the sphere we need to

- compute the convex hull of this set of points;
- the normals to each of the faces give the Voronoi points;
- the face connectivities provide the connectivities between Voronoi points.

Now, each of these steps has well optimized algorithms which are implemented in numerical software like Matlab. Just to give you an idea, it is possible to compute the Voronoi diagram of a set of points in about 5 seconds.

Now let’s give some mathematical justifications for the above steps. First, let’s begin with the definition of the Voronoi diagram. A *Voronoi* diagram associated to a set of points on the sphere is a partition of the sphere into regions associated to the points, such that the points in region are closer to the -th point than to any of the other points. To any Voronoi diagram we associate a set of Voronoi points which are the vertices of the Voronoi cells. This Voronoi diagram is the dual of the Delaunay triangulation of the set of points. Recall that the Delaunay triangulation has the property that the circumcircle of each triangle does not contain any of the other remaining points in its interior. Now, if we have the Delaunay triangulation, the Voronoi points are exactly the projection on the sphere of the circumcenters of the triangles.

So a first step of any algorithm would be to compute the Delaunay triangulation. Now, what’s up with the convex hull step? It turns out that in the case of the family of points on the sphere, taking the convex hull, i.e. finding the faces of the convex polyhedron which is determined by these points on the sphere, gives the Delaunay triangulation. Let’s justify this. Suppose is a triangular face in the convex hull of the points. In fact, all faces in the convex hull can be given as triangles, since a polygonal face can be divided into triangles. Suppose that the circumcircle of contains another point in its interior. Since the circumcircle is the intersection of the plane with the sphere, it follows that is strictly above the plane (opposite to the center of the sphere), which means that is outside the convex hull, which is not possible. Thus, the triangulation given by the convex hull is Delaunay. (of course, in the above, we left out the case where all the points are in a hemi-sphere, but this can be treated in the same way)

Thus the convex hull gives the Delaunay triangulation. The circumcenters of the corresponding triangles projected on the sphere are the Voronoi centers. To validate the second step we only need to notice that the the line passing through the center of the sphere, orthogonal to the face, passes through the circumcenter of the triangle. The use of the face normals is obvious once this fact is observed. In order to prove this, denote with the center of the sphere and with the projection of on . We have as are on the sphere, so the right triangles are congruent, because they also have as a common side. Thus and is the circumcenter of .

The connectivity property comes naturally from the definition of the Voronoi cells. The only way that the arc between the circumcenters of and can be separated is to have another Voronoi cell which intersects . But that would imply that one of the circmcircles of or should contain another point in their interior, which contradicts the fact that we have a Delaunay triangulation.

After these arguments justifying the validity of the method, here are some numerical results obtained with Matlab.