Computing the Voronoi diagram of a region inside a box

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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.

Answer 5

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.