# Computing the Voronoi diagram of a region inside a box

I am facing a problem as follows: I have a box full of points with a certain unknown distribution and I would like to calculate its Voronoï Diagram. The problem is that the number of points is so huge that this may be impossible to do for the full distribution.

Therefore, I have planned to do that for just a region inside the box, where the number of points was not that big. In order to do so I need to know how to calculate the minimum region that may affect to the Voronoi diagram of a certain smaller region inside that box.

In other words, I would like to compute the Voronoï diagram of the points inside the small cube of the figure below that fits with the Voronoï diagram that would have the points of the full box storing the smallest Voronoï diagram possible on the memory. • I have seen people, especially in astrophysics, computing voronoi tessalation of huge point sets. See the work of Volker Springel. E.g., there is even an open source code here github.com/regonzar/paravt which may be useful to you. Also see arxiv.org/abs/1601.06429 – cfdlab Jul 26 '16 at 9:28
• I would still need some way to do that using the strategy of the small cube, but thank you very much. I will take a look at it. – ccorbella Jul 26 '16 at 9:44
• Sorry, I have problems trying to understand this sentence "I would like to compute the Voronoï diagram of the points inside the small cube of the figure below that fits with the Voronoï diagram that would have the points of the full box storing the smallest Voronoï diagram possible on the memory." – nicoguaro Jul 27 '16 at 23:16
• Sorry, I just meant that I would like to compute the Voronoï diagram of the small cube, taking into account that this should be the same in that region that the diagram I would obtain if I calculated using all the points of the box. In order to do so, I expect to need more points than those that are inside the box (otherwise, I do not believe it may fit with another cube if I follow the same strategy there) but I would like to store as less as possible points. – ccorbella Jul 28 '16 at 7:32
• @ccorbella Not an answer but with what tool you did provide this nice figure please? Maybe add in a caption the tool. – Jan Hackenberg Aug 9 '16 at 9:14

## 5 Answers

To compute the Voronoi diagram of huge (>100 millions) sets of points, you can use the following algorithm:

1) create a kd-tree with all the points
2) for each point p [in parallel optionally]
N = 10
while not finished
compute the N nearest neighbors of the point p
compute the intersection of the N half-spaces defined by p
and the neighbors
if there is a neighbor further away than
twice the radius of the ball centered on p and
bounding the intersection, finished = true
N = N * 1.5
// when exiting the loop, the computed intersection
// corresponds to the Voronoi cell of p, because no other bisector
// can contribute to the Voronoi cell.


The algorithm is explained with further details in my article. It can be trivially parallel-ized (just add "#pragma omp parallel for" before the main loop), since there is no data dependency. It is implemented in my GEOGRAM C++ programming library (together with a memory efficient Kd-Tree that scales up to more than 100 million points). Note that in GEOGRAM there is also a parallel standard Delaunay/Voronoi implementation that works well with up to 100M sites.

Concerning parallel implementation of the classical (Boywer-Watson) algorithm, the GEOGRAM implementation is documented here (see also the associated c++ source file that has extensive comments). I have no published article about it, I will write one if time permits. The main idea is to use spinlocks associated with the tetrahedra to ensure that only one individual thread can modify a tetrahedron.

• First of all, thank you for your answer. I am sorry to tell you that your article seems to be not longer uploaded (at least the full version from the web you link). Anyway, could you explain me how would you implement that Voronoi diagram algorithm in parallel? – ccorbella Aug 8 '16 at 15:26
• Thank you for notifying that, I have fixed the link (the new link has the PDF, click on the PDF icon to get it). I have also added short explanation / link to parallel Delaunay implementation. – BrunoLevy Aug 8 '16 at 19:40
• Note: it will work well provided that the points are uniformly distributed (performance may drop if you have a high variation of points density). – BrunoLevy Sep 2 '16 at 11:40

Seems like the experts are not answering your question so I will try to provide an idea. But before I do that I strongly suggest that you look up in the literature for some sophisticated methods that have been already developed. However, without guaranteeing that this is a good or fast or efficient suggestion, I propose the following methodology. Keep in mind, I may have made some mistakes, so I do not guarantee that everything is fully correct, but I hope the idea of the method gives you some approach that will help you solve your problem.

Let $V$ be the set of your points in the whole "big" cube. Fix your "small" cube $C$ somewhere in the big cube and let $V_C$ be the set of points that are contained in $C$, i.e. $V_C = V \cap C.$ Initially set $V'_C=V_C$.

Step 1: Generate the Voronoi diagram $Vor(V'_C)$. For each point $v \in V'_C$ denote by $Vor(v)$ its Voronoi cell, which is a convex polyhedron in three-space. Furthermore, denote by $W(v)$ the vertices of the Voronoi cell centered at $v \in V'_C$ and by $W(V'_C) = \cup_{v \in V'_C} W(v)$ the vertices of all Voronoi cells from the Voronoi diagram $Vor(V'_C)$.

Step 2: Color all points from $V'_C$ and all Voronoi vertices $W(V'_C)$ white.

Step 3: For each Voronoi vertex $w \in W(V'_C)$ draw the Delaunay sphere centered at $w$, that is the sphere with center $w$ and radius the distance between $w$ and one of the points from $V'_C$ whose Voronoi cell has $w$ as a vertex (it doesn't matter which point, there are several but the result is always the same).

Case 3.1. If the Delaunay sphere of $w$ is contained in the cube $C$, color $w$ black.

Case 3.2. If the Delaunay sphere is not contained in the cube $C$ but it doesn't contain any point from $V$ in its (open) interior, color the point $w$ black.

Case 3.3. If the Delaunay sphere of $w$ contain points from $V$ in its (open) interior, (1) add the points from $V$ contained in the interior of the sphere to the set $V'_C$ and (2) keep the color of the point $w$ white.

Step 4: For each point $v \in V'_C$ check if all Voronoi vertices $W(v)$ of its Vornoi cell are black. If not all of them are black, keep the color of $v$ white. If they are black, color $v$ black.

Step 5: Check if all points of the original set $V_C$ are black.

Case 5.1. If they are all black, the Voronoi diagram $Vor(V'_C)$ restricted to the cube $C$ is the local portion of the global Voronoi diagram $Vor(V)$ restricted to $C$. End.

Case 5.2. If there are white vertices in $V_C$, then go back to Step 1. In Step 1, when generating the new Voronoi diagram $Vor(V'_C)$, one keeps the Voronoi cells around black points from $V'_C$ the same, keeps all black Voronoi vertices from $W(V'_C)$ and makes alteration only in relation to the white ones.

I hope this helps.

The simplest way of doing this is to surround your iner box with a bigger box that contains at least all the nearest neighbours of the points within your inner box. Note that there will be an issue arising when the inner box is close to the edge of the encompassing data box: you have no external points.

Calculating a Voronoi/Delaunay tessellation may be more subtle than you might think. One of the issues is how to decide precisely whether a point is on one side or other of a tessellation plane/line.

There is the very complete "CGAL" C++ library for doing this at http://www.cgal.org/. My colleagues and I have used this in several published papers in astrophysics: it appears to be rock-solid in addressing all the potential pitfalls in creating these tessellations.

• Thank you very much for your answer. Then, my question should be rewritten, if you want so, as "how to find the nearest neighbours points of the ones that are inside of the cube doing the least number of calculations". My problem was basically that one. Do you know any way to do so? – ccorbella Aug 2 '16 at 7:03
• How does one decide how big the outer box should be? If it's not big enough, you may not get the restriction of the full diagram to the original small box. I think the right decision weather a vertex of a voronoi cell from the local diagram is a vertex of the global diagram, i.e. a vertex that is not going to be altered by any future recalculation of the local voronoi diagram, is based on wether the interior of the corresponding delaunay sphere contains any points from the total set of points or not. This is exactly the definition of a delaynay cell, which is dual to a voronoi cell. – Futurologist Aug 2 '16 at 10:57
• @ConradCorbellaBagot For calculation of the nearest neighbors in a large n-dim data set, very efficient algorithm exists. Maybe you want to state what you are really interested at. – Bort Aug 3 '16 at 9:34
• The Voronoi/Delaunay tessellations are well defined on both an infinite point set and on a point set that is bounded, but not on a subset of points of a larger set. For such subsets you have to make an arbitrary compromise decision. In cosmology where we have a finite box in an otherwise infinite universe we choose periodic boundary conditions. When doing image analysis on part of an image I "dress" the boundary with the first nearest neighbours of the points defining the boundary (there are complications to consider). I find that going the next nearest neighbour has relatively little gain. – JonesTheAstronomer Aug 4 '16 at 14:21
• There is a detailed presentation of some of this in adsabs.harvard.edu/abs/2011MNRAS.416.2494P which is freely downloadable from both the journal and the arXiv. It also has a discussion of nonlinear Kriging reconstruction of a density field using these tessellations. The data source here is astronomical, but the discussion is quite generic (for 3 dimensional data sets). – JonesTheAstronomer Aug 4 '16 at 14:28

I understand your question as: I want to draw Voronoi diagram for a subset of points such that it is same as the one which is obtained when considering the complete set of points. Voronoi diagrams are drawn by first joining neighboring points and then drawing a plane perpendicular to the line at the midpoint. You do this for all the nearest neighbors and you have a voronoi diagram in the neighborhood of a point. Do this for all the points and you have voronoi diagram for all the points. You see, voronoi diagrams are defined locally. There is no second nearest neighbor or third nearest neighbor effect. Only first nearest neighbor effect. So all you have to do to get a voronoi diagram with a subset of point is to identify the points in the subregion of interest, connect them to all of their respective nearest neighbors, and draw a plane passing through midpoints of these line segment and perpendicular to the line segment. This diagram will be same for a local region whether you consider a sub region or complete region.

• Two questions. First of all, "There is second nearest neighbor or third nearest neighbor effect" should be "There is not...", right? – ccorbella Aug 2 '16 at 6:51
• Yes. Thank you for pointing that out. I will update the answer.What is the second question? – Kaustubh Kaluskar Aug 2 '16 at 6:58
• Sorry, I was editing the comment xDD And the other one, I do know that Voronoi diagrams are defined locally (under other conditions my question would have no answer). My question should be rewritten, if you want so, as "how to find the nearest neighbours points of the ones that are inside of the cube doing the least number of calculations". Do you know any way to do so? Anyway, thanks a lot for your time. – ccorbella Aug 2 '16 at 6:59
• I use knnsearch function in MATLAB. My typical dataset is around 1.5 million points and I do it on my laptop. From Mathworks website: 'IDX = knnsearch(X,Y) finds the nearest neighbor in X for each point in Y. IDX is a column vector with my rows. Each row in IDX contains the index of nearest neighbor in X for the corresponding row in Y.' Here X would be your full dataset and Y are points inside the cube. – Kaustubh Kaluskar Aug 2 '16 at 7:05
• First of all, thank you very much once again. I am going to try it tough I am unsure if I am able to run an algorithm as a kd-tree in such a huge number of points, as it says that "For large dimensions (20 is already large) do not expect this to run significantly faster than brute force. High-dimensional nearest-neighbor queries are a substantial open problem in computer science."... – ccorbella Aug 2 '16 at 7:36

I suggest you to have a visual and intuitive approach using Grasshopper for Rhinoceros3D. Although Rhinoceros is a commercial CAD package and Grasshopper is a plugin for it, you can run plugins for free with no limitations and make your experiments (unlicensed Rhino3D limits only the saving of Rhino files). Grasshopper includes a large number of mathematical functions used in a canvas, and 3D Voronoi diagrams is one of them. • The link is interesting but this does not answer the question. – BrunoLevy Aug 8 '16 at 14:01
• That was not what I was asking for, but thank you for the tool. – ccorbella Aug 8 '16 at 15:30