15

I need to compute the Voronoi diagram of a set of points in $R^n$. I'm quite unschooled on the topic, could someone point me to the right references so that I can

a) understand the theory behind it;

b) implement an actual algorithm to compute it.

P.S. Googling I found out I need to know about euclidean graphs, but I couldn't find any decent introduction to it. Pointers are appreciated!

YCor
  • 60,149
Alessandro
  • 253

3 Answers3

19

Often Voronoi diagrams are constructed from their duals, Delaunay triangulations. The Delaunay triangulation of a point set in $d$ dimensions can be obtained from the convex hull of a lifting of the points into $d+1$ dimensions. This relationship is explained in many sources, including The Handbook of Discrete and Computational Geometry (Chapters 22, 23).

Implementing this on your own is quite a project, so you might first investigate what is available. Qhull performs the computations; see www.qhull.org. $d$-dimensional Delaunay triangulation computations are part of CGAL, specifically the module on dD Convex Hulls and Delaunay Triangulations. You could use polymake to convert from the Delaunay triangulation to its corresponding Voronoi diagram. Finally, there are papers written on specific versions of your problem, for example: "An Explicit Solution for Computing the Euclidean $d$-dimensional Voronoi Diagram of Spheres in a Floating-Point Arithmetic," by Marina Gavrilova, 2003.

(Incidentally, "Euclidean graphs" play a relatively minor role in this topic.)

The Amplitwist
  • 1,446
Joseph O'Rourke
  • 149,182
  • 34
  • 342
  • 933
  • Thank you very much for your answer! I read some more about the topic and indeed it seems quite complex. Furthermore, the run-time complexity of the algorithms is probably too high for my needs -I would like not to go over the n^2 worst case complexity. Do you know of any fast approximate algorithms that can help me? – Alessandro Dec 15 '10 at 14:05
  • 1
    @Alessandro: It may not be as bad as you fear. See, e.g., "Higher-dimensional Voronoi diagrams in linear expected time," by Rex Dwyer, and track papers that cite this into the future. http://portal.acm.org/citation.cfm?id=73869 – Joseph O'Rourke Dec 15 '10 at 14:10
  • 1
    I'm afraid it is as bad as Allessandro fears. The main problem is the number of the hyperplanes that these algorithm generates. For a very simple example, I used qhull for 40 points in 15-d and the voronoi decomposition has more than 5 million hyperplanes

    Voronoi diagram by the convex hull of 40 points in 16-d:

    Number of Voronoi regions: 40 Number of Voronoi vertices: 1361237

    Statistics for: rbox 40 s D15 | qvoronoi s p TO test

    Number of points processed: 40 Number of hyperplanes created: 5539466 Number of facets in hull: 2774645 Number of distance tests for qhull: 4729150

     – Adel Ahmadyan Apr 05 '13 at 05:12
5

I think the "lifting of the points into $d+1$ dimensions" referred to in Joseph O'Rourke's answer means something like this: $(x_1,\dots,x_d) \mapsto (x_1,\dots,x_d,x_1^2+\cdots+x_d^2)$. Then the edges of the convex hull of those points in $d+1$ space connect the images of points that are connected in the Delaunay triangulation.

Michael Hardy
  • 11,922
  • 11
  • 81
  • 119
  • @Michael: Yes, that's it; thanks for clarifying. But it is only the lower hull that projects to the Delaunay triangulation. – Joseph O'Rourke Dec 15 '10 at 16:00
  • Doesn't the whole convex hull lie above its boundary in this case? Thus the whole boundary is the lower boundary; the convex hull is not bounded above. – Michael Hardy Dec 15 '10 at 17:42
  • Oh! You meant when the set of points whose Delaunay triangulation you're trying to find it bounded. I was picturing an unbounded set of points. – Michael Hardy Dec 15 '10 at 17:45
  • I like the "sphere" version of the geometric proposition involved, since it lacks these messy asymmetries. Take a finite set of points on a sphere. Form their convex hull. Then the edges of that polygon are the Delaunay triangulation and the dual is the Voronoi tesselation. No need to talk about "lower" and "upper". – Michael Hardy Dec 15 '10 at 17:47
  • Is the sphere version in the literature somewhere? I can imagine it having a nice proof, but I'd have to think about the specifics. – Michael Hardy Dec 15 '10 at 17:48
  • 1
    @Michael: Whether the asymmetries are "messy" may be a matter of taste :-). The top of the paraboloid-lifted hull records furthest-point diagram adjacencies, while the bottom of the hull records closest-point adjacencies. For the sphere viewpoint, see Kevin Q. Brown: Voronoi Diagrams from Convex Hulls. Inf. Process. Lett. 9(5): 223-228 (1979). – Joseph O'Rourke Dec 15 '10 at 18:33
  • Thank you. I've put in an interlibrary loan request for Brown's paper. – Michael Hardy Dec 15 '10 at 20:07
4

R (a favorite open source, object oriented interactive statistical package) has what appears to be a suitably fast algorithm, but I've only used it in a problem with 20 some odd points or so.

As the others have mentioned above, the Voronoi tessellation or mosaic is computed using the (dual graph which is the) Delaunay triangulation. The R function 'veronoi.mosaic' in the package 'tripack' calls a FORTRAN routine. If its convenient enough to do it in R then just call the veronoi.mosaic function (details on getting R and getting the tripack package follow). Otherwise, if you'd rather just use the internal FORTRAN subroutine, you can figure out the meaning of the arguments by looking at the veronoi.mosaic function in R (at the point at which the FORTRAN subroutine 'veronoi' is called) and work backwards.

In any case, all of its open source and works like a charm.

getting R:

http://cran.r-project.org/

the tripack package

http://cran.r-project.org/web/packages/tripack/index.html

Enjoy!

Best Regards, Grant Izmirlian

izmirlig
  • 61
  • I would just caution that robustness of the computation can be an issue. I have no idea if R uses arbitrary precision arithmetic computations when necessary (as does, e.g., CGAL), but it warrants attention. And that second link above is to "constrained two-dimensional Delaunay triangulation package," not $d$-dimensions for abitrary $d$. – Joseph O'Rourke Dec 15 '10 at 19:32