On this page you will find software that I made for demonstration purposes in connection with the book Algorithmische Geometrie by Klein, Driemel, and myself. You can download a package with C++ source files here (version 0.06). These implement a quad-edge data structure and several efficient algorithms to compute Delaunay triangulations and Voronoi diagrams, about a thousand times faster than Vorosketch. In particular, the package contains the following files.
This file implements the quad-edge data structure (QEDS). The QEDS stores a connected planar graph and its dual. Thus, when we compute a Voronoi diagram and store it as a QEDS, we get the Delaunay triangulation for free, and vice versa (provided we deal with cases of more than three points on a circle appropriately, which we will). The difficulty with maintaining a QEDS is that for each operation that changes the graph, you have to change several pointers. It is typically only a constant number per operation, but enough to worry: how can you be sure that you always connect things up properly and never introduce an inconsistency? The answer is: all operations are reduced to a single operation called splice, which guarantees that the QEDS is kept in a consistent state at all times. This implementation provides functions for everything we need (cutting faces, sewing one graph onto another, edge flips etc.), all based on splice. For more information, see the QEDS page.
This file implements the necessary basic geometric operations with exact integer arithmetic, using only additions, subtractions and multiplications; no divisions, no roots, no trigonometric functions. In particular, this code implements:
This file implements the randomised incremental construction algorithm for Delaunay triangulations with the Delaunay DAG point location data structure, see Section 6.2 of the book.
This file implements a variant of the randomised incremental construction algorithm, which is briefly alluded to in Section 7.4, combining techniques from there and from Section 6.2. The points are randomly distributed over rounds of increasing size, and within each round, the points are inserted in the triangulation in order along a Pólya/Sierpinski plane-filling curve. We do not use a point location data structure: to find the triangle T that contains the point to be inserted, we simply walk through the triangulation from the previous point that was inserted. Intuitively, the random points from the previous rounds usually ensure that the existing triangles are well-shaped, while the plane-filling curve order from the current round usually ensures that we do not need to walk far: thus we need to cross only few triangles to reach T.
This file implements the Divide-and-Conquer algorithm for Voronoi diagrams, see Section 6.4 of the book.
These are figures which are referred to in the comments in agb-dtvd-dc.cpp. The figures illustrate how the algorithm assembles the Voronoi diagram in the conquer step, using various pointers.
This file contains several simple pre- and postprocessing functions for the triangulation and Voronoi diagram algorithms:
Compile this file to get a working program. For example:
g++ -D INCLUDE_AGB_DTVD_DC -o agb-dtvd agb-dtvd-main.cpp
should produce an executable agb-dtvd that uses the divide-and-conquer algorithm. Replace DC by RICDD or RICPC to use one of the randomised incremental construction algorithms instead.
To use the program, create a text file that contains integer numbers separated by spaces; the first number is the number n of points in the file, the next 2n numbers give, for each point, its x-coordinate and its y-coordinate. For example, a very simple input could be:
3 -1 -1 0 1 1 -1
If you call this file points.txt, then this:
agb-dtvd <points.txt >vd.ps
should produce a Postscript file vd.ps with a colourful drawing of the Voronoi diagram and the Delaunay triangulation of the point set. You will find that with each of the algorithms mentioned above, the program can handle thousands of points within a second—try, for example, the file 10000points.txt in the package.
This is just a simple demo. The software does not attempt to deal gracefully with all possible errors or special cases in the input. In particular, the coordinates of the input points should be integers between -17605 and 17605 and the input should contain at least two distinct points; moreover, in the case of the RIC algorithms, they should not all lie on one line. To handle bigger or non-integer coordinates, you should replace the implementation of coordinates in agb-simplepoint.cpp by something more powerful.
I have tried not to use system-specific C++ features, but I have not compiled these sources on a large variety of systems to check. Please let me know if you have trouble compiling these sources.
For sure, this software is not the fastest code around. Its humble purpose is just to show that (and how) these algorithms work, and that it saves a lot of time as compared to naive solutions. Therefore I invested time in writing elaborate comments in the code. But when it comes to performance, there is lots of room for optimisation. For example, one may use a data structure for triangulations instead of a QEDS for general planar graphs, or one may alternate splitting directions in the divide-and-conquer algorithm to decrease the length of the bisectors that need to be traced. Of course, the plane sweep algorithm would also deserve to be included in this package.