THE BOOST.POLYGON VORONOI LIBRARY


The Boost.Polygon Voronoi library provides functionality to construct the Voronoi diagram of a set of points and linear segments in 2D space with the following set of limitations:
  • coordinates of the input points and endpoints of the segments should be of the integer type;
  • input segments should not overlap except their endpoints.
While the first restriction is permanent (it allows to give the exact warranties about the output precision and algorithm execution), the second one may be resolved using the Boost.Polygon functionality. The strong sides of the library and main benefits comparing to the other implementations are discussed in the following paragraphs.

Robustness and Efficiency

Let's explain a bit those terms. The efficiency is simply measured by the time it takes the algorithm to execute. The robustness is a bit more harder to explain. But those of you who had the experience with the following situations would understand what it doesn't mean: application segfaults randomly, algorithm output contains degeneracies or is completely invalid (e.g. a point is considered to be outside of the polygon, while should be inside). In other words robust implementation doesn't fail and produces expected output in 100% of cases, thus user can rely on it. Robustness is a weak place of the most non-commercial implementations of any complex geometric algorithm. The main issues could be divided onto two main categories: memory management issues, numeric stability issues. Our implementation avoids the first type of issues using pure STL data structures, thus you won't find any operator new in the code. The second category of problems is resolved using multiprecision geometric predicates. Even for commercial implementations usage of such predicates usually results in a huge performance slowdown. Here is another strong side of the Boost.Polygon Voronoi library: we avoid multiprecision computations in 95% of cases using extremely fast floating-point predicates. Yes, those are not always exact, but we developed the relative error arithmetic apparatus to identify them and switch to the higher precision predicates when required.

Precision of the Output Structures

One of the extremely important results of using two types of predicates is that library efficiently computes relatively precise coordinates of the output geometries. Here we will explain a bit what exactly "relatively precise" means and how the received output may differ from the theoretically correct one (here and after we assume that output coordinate type is the IEEE-754 floating-point type).

Voronoi implementation guaranties that the relative error of the coordinates of the output geometries is always not higher than 64 machine epsilons (6 bits of mantissa), while in many cases it is slightly less. That also means that using floating-point type with the larger mantissa will produce more precise output. Let's consider the following example: the output Voronoi vertex has double (53-bit mantissa) x-coordinate equal to 1.0, then the absolute error is at most 2^-53 * 2^6 = 2^-47 and the exact value of x-coordinate lies in the range [1.0 - 2^-47, 1.0 + 2^-47]. For x-coordinate equal to 2^31, the absolute error will be at most 2^-53 * 2^31 * 2^6 = 2^-16 and the exact value of x-coordinate lies in the range [2^31 - 2^-16, 2^31 + 2^16]. For the output Voronoi vertex with long double (64-bit mantissa) x-coordinate equal to 2^31, the absolute error will be at most 2^-64 * 2 ^31 * 2^6 = 2^-27 and the exact value of x-coordinate lies in the range [2^31-2^-27, 2^31+2^-27]. If you'd like to become master of the absolute and relative error try this article.

During the finalization step the implementation unites Voronoi vertices whose both coordinates are situated within the relative error range equal to 128 machine epsilons and removes any Voronoi edges between them. That is the only case that might cause differences between the algorithm output topology and theoretically precise one.  Now let's see what is the practical impact of this. Consider the following example: we are going to construct the Voronoi diagram of our Solar System. The radius of our Solar System is approximately 2^42 metres, and we are going to snap it to the integer grid of [-2^42; 2^42] x [-2^42; 2^42].  Let's choose the long double (64 bit mantissa) output coordinate type, then the maximum absolute error for the output geometries within our Solar System will be on its boundaries and equal to 2^-64 * 2^42 * 2^6 = 2^-18 metres. In the output we are going to consider vertices with both coordinates that are within 2^-17 metres (8 micrometres) distance to be equal. That distance is equal to the size of a bacteria and is relative to the Solar System size.

Fully Functional with Segments

There are not many implementations of the Voronoi diagram construction algorithm that could handle linear segment inputs, even considering the commercial libraries. Support of the segments allows to discretize any input geometry (circle, ellipse, parabola). Of course as the result those might have floating-point coordinates, but that is resolved using scaling and snapping to the integer grid. This functionality is very handy as it allows to compute the medial axis transform of the arbitrary set of input geometries. So one may start using it for the next generation pattern recognition or computer vision project.

Basic and Advanced Usage Cases

The main library header voronoi.hpp defines the following static functions to integrate the Voronoi library functionality with the Boost.Polygon interfaces:

template <typename Point, typename VB>
size_t insert(const Point &point, VB *vb)
Inserts a point into the Voronoi builder data structure.
Point type should model the point concept.
Returns index number of the inserted site.
template <typename PointIterator, typename VB>
void insert(PointIterator first,
            PointIterator last,
            VB *vb)
Inserts an iterator range of points into the Voronoi builder data structure.
Corresponding point type should model the point concept.
template <typename Segment, typename VB>
size_t insert(const Segment &segment, VB *vb)
Inserts a segment into the Voronoi builder data structure.
Segment type should model the segment concept.
Returns index number of the inserted site.
template <typename SegmentIterator, typename VB>
void insert(SegmentIterator first,
            SegmentIterator last,
            VB *vb)
Inserts an iterator range of segments into the Voronoi builder data structure.
Corresponding segment type should model the segment concept.
template <typename PointIterator, typename VD>
void construct_voronoi(PointIterator first,
                       PointIterator last,
                       VD *vd)
Constructs Voronoi diagram of a set of points.
Corresponding point type should model the point concept.
template <typename SegmentIterator, typename VD>
void construct_voronoi(SegmentIterator first,
                       SegmentIterator last,
                       VD *vd)
Constructs Voronoi diagram of a set of segments.
Corresponding segment type should model the segment concept.
template <typename PointIterator,
          typename SegmentIterator,
          typename VD>
void construct_voronoi(PointIterator p_first,
                       PointIterator p_last,
                       SegmentIterator s_first,
                       SegmentIterator s_last,
                       VD *vd)
Constructs Voronoi diagram of a set of points and segments.
Corresponding point type should model the point concept.
Corresponding segment type should model the segment concept.

This means that it's possible to construct the Voronoi diagram with the following two lines of code (if corresponding input types satisfy the Boost.Polygon concept model):

voronoi_diagram<double> vd;
construct_voronoi(points.begin(), points.end(), &vd);

The library also provides the clear interfaces to associate user data with the output geometries and efficiently traverse Voronoi graph. More details on those are covered in the basic Voronoi tutorial. Advanced usage of the library with the configuration of the coordinate types is explained in the advanced Voronoi tutorial. The library also allows users to implement their own Voronoi diagram / Delaunay triangulation construction routines based on the Voronoi builder API.

No Third Party Dependencies

Yes, the library doesn't depend on any 3rd party code. Even more than that there is only one dependency on the Boost libraries: boost/cstdint.hpp. All the required multiprecision types functionality is implemented as part of the library and is not exposed to the user. Considering the fact that Voronoi implementation consists of just 7 headers (3 public and 4 private) it is easy to compile it within a minute after download. On the other hand voronoi.hpp header provides integration routines with the Boost.Polygon concepts and models with a drawback of additional dependencies.

Extensible for the User Provided Coordinate Types

Our implementation is coordinate type agnostic. That means that as soon as user provided types satisfy the set of restrictions of the Voronoi builder coordinate type traits and implement methods required by the library, no changes are required neither to the algorithm, nor to the implementation of the predicates. So it's possible to construct Voronoi diagram for the 256-bit integer input coordinate type and 512-bit output floating-point type without making any changes to the internal code.

Bright Future

Below one may find the list of the main directions for the future development of the library.
High-priority tasks that already have approximate implementation plan are the following (some of those may be proposed as future GSoC projects):
  • Implementing Delaunay triangulation data structure.
    Note: only data structure needs to be implemented that properly processes events provided by the Voronoi builder.
  • Implementing medial axis transform data structure.
    Note: in general case the Voronoi diagram has completely the same geometry as the medial axis (they are 100% equal), however for many applications user is not interested in the Voronoi edges inside the hole regions. The main point of this data structure is to automatically filter Voronoi edges that belong to those areas.
  • Voronoi diagram data structure could be used to find K nearest neighbors of N sites in O(N*K*log(K) + N*log(N)) time. The return value would be a list of the k nearest neighbors for each site.
  • Using the r-tree data structure built on top of the bounding rectangles around the Voronoi cells to answer the nearest neighbor queries in log(N) time, where N is the number of the Voronoi cells.
    Note: there should be r-tree data structure available soon as part of the Boost libraries.
  • Providing interface to retrieve the convex hull of a set of points and segments from the Voronoi builder once the Voronoi diagram is constructed in O(N) time.
  • Providing serialization utilities for the Voronoi diagram data structure.
High-priority tasks to be considered:
  • Dropping the restriction on the non-intersecting input geometries.
  • Integration of the Voronoi diagram data structure with the BGL (Boost Graph Library).
  • Support of the other types of distance metrics.
  • Construction of the constrained Delaunay triangulation.
  • Support of the circular input geometries.
Based on the community suggestions priorities may be changed.

Theoretical Research

Voronoi was developed as part of the Google Summer of Code 2010. The library was actively maintained for the last two years and involved strong mathematical research in the field of algorithms, data structures, relative error arithmetic and numerical robustness. Nowadays one can often read a scientific article that contains non-practical theoretical results or implementation with benchmarks nobody else can reproduce. The opposite story is with the Boost.Polygon Voronoi library. We provide pure implementation and benchmarks one may run on his PC. In case community finds it useful we will incrementally add more documentation on the theoretical side of our implementation. The authors would like to acknowledge the Steven Fortune's article "A Sweepline algorithm for Voronoi diagrams", that contains the fundamental ideas of the current implementation.
 
Copyright: Copyright © Andrii Sydorchuk 2010-2012.
License: Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)