Voronoi Advanced Tutorial

This tutorial consists of two parts. The first one provides two examples of a real world problems that default configuration of Voronoi library is capable to solve. By default configuration we mean the one that accepts signed 32-bit integer and outputs floating-point (64-bit double) coordinates. We provide those examples to convince even the most skeptical users that they probably don't need to configure library for higher-precision input or output coordinate types. However if the posed problem really requires those, fully featured configuration of both input and output coordinate types is provided in the second part of this tutorial.

Red Planet

Problem Statement



Lets imagine that NASA is planning to send a new robot to Mars. Each day the center situated on Earth will send a destination point coordinates the robot needs to reach by the end of the day. Of course we'd like to save as much energy as possible thus choosing the shortest possible path. This would be a straight line in a perfect world (we don't consider curvature of surface), but situation becomes more complicated as there are some rocks and wells on Mars our robot can't go through. Behind of those our robot has some dimensions that might not allow it to pass narrow places.

Application of Voronoi diagram

The problem above could be solved using Voronoi diagram. The first stage would be to discretize obstacles (rocks and wells) with polylines. Afterwards we will compute Voronoi diagram of the input set of segments. As each Voronoi edge is equidistant from the two closest sites we are able to filter edges the robot won't be able to pass due to it's dimensions. The last step would be to run a bit optimized A* algorithm to find the shortest or at least suboptimal path and we are done.

Discretization of input geometries

To show how good is the default input coordinate type provided by the Voronoi library we will discretize the whole area of Mars. That will be approximately 1.44 *  10^8  square kilometres that is equal to 1.44 *  10^18  square centimetres, which could be snapped to the integer grid with a side of 1.2 * 10^9 centimetres.  To make the Voronoi graph precise on the boundaries of that grid we will replicate input map 9 times (3x3), thus Voronoi diagram within a central piece will provide us with a correct connectivity graph. This step will increase the size of our grid to 3.6 * 10^9 centimetres that is less than 2^32. So we are able to discretize the Red Planet's surface within a 1 centimetre precision using the default input coordinate type (signed 32-bit integer). That would imply maximum absolute error to be equal up to 0.5 centimetres per coordinate. Considering the radius of our robot to be 0.3 metres and for security reasons avoiding any large enough obstacles that are within 1 metre distance from it that error would be irrelevant.

Output analysis

Estimates of the resulting Voronoi diagram precision were already explained here. So to avoid those computations again we will simply state that the maximum absolute error of the output geometries will be on the grid boundaries and will be equal to 2^-16 centimetres, which is approximately equal to 150 nanometres and is 100 times larger than a radius of a complex molecule. We would like to notice that the absolute error of the discretization step is much higher than the one produced by the Voronoi diagram construction algorithm.

VLSI Design

Problem Statement



Very-large-scale integration (VLSI) is the process of creating integrated circuits by combining thousands of transistors into a single chip. Considering the fact that it may take weeks or months to get an integrated circuit manufactured, designers often spend large amounts of time analyzing their layouts to avoid costly mistakes. One of the common static analysis checks is minimum distance requirement between the components of an integrated circuit (distance should be greater than specified value).

Application of Voronoi diagram

It appears that the minimum distance between components of the input set of points and segments corresponds to the one of the Voronoi diagram edges. This means that we can iterate through each edge of the Voronoi graph, extract the pair of input geometries that form it and find the distance between those. As the total amount of such edges is O(N) value (N - is the number of input geometries) the minimum distance could be efficiently find in a linear time once we construct the diagram.

Discretization of input geometries

The average size of the modern CPUs is around 2.5 x 2.5 centimetres. Snapping this to the 32-bit integer grid will give discretization precision of 2.5 / 2^33 centimetres or 3 picometres that is 10 times smaller value than radius of an atom. That would be probably good enough precision even for a few next generations of processors.

Output analysis

The maximum absolute error of the output geometries will be 2.5 / 2^47 centimetres or 0.2 femtometres that is 10 times smaller value than the radius of an electron. However in this particular case we are not interested in the precision of the output, rather in its topology. As it was noticed on the Voronoi main page very small edges are removed from the Voronoi diagram. However user should not worry because the edge that correspond to the minimal distance won't be among those. That means that we would be able to 100% correctly identify a pair of closest objects within the discretization precision.

Conclusions

The above two examples show usage of the default Voronoi coordinate types in the macro and micro world. The main point of those was to give the user understanding of a scale that the default coordinate types provide. There are two main points we didn't mention before, but that would be relevant to the most real world problems:
The second statement means that there is actually no point to provide implementation that operates with floating-point input coordinates, because those always could be snapped to the integer grid. In case the user is not satisfied with the precision that the 32-bit integer grid provides or would like to retrieve coordinates of the output geometries within a smaller relative error, follow the next paragraph.

Voronoi Coordinate Types Configuration

In the following chapter we are going to extend input coordinate type to the 48-bit signed integer and output coordinate type to the 80-bit IEEE floating-point type (long double). The code for this chapter is available in voroni_advanced_tutorial.cpp. While it won't be possible to compile it using the MSVC compiler (it doesn't support 80-bit long double type; ieee754.h header is required), it should give a clear understanding of how the library supports the user provided coordinate types.

So the main step would be to declare the voronoi coordinate type traits that satisfy set of restrictions explained here:

struct my_voronoi_ctype_traits {
  typedef boost::int64_t int_type;
  typedef detail::extended_int<3> int_x2_type;
  typedef detail::extended_int<3> uint_x2_type;
  typedef detail::extended_int<128> big_int_type;
  typedef fpt80 fpt_type;
  typedef fpt80 efpt_type;
  typedef my_ulp_comparison ulp_cmp_type;
  typedef my_fpt_converter to_fpt_converter_type;
  typedef my_fpt_converter to_efpt_converter_type;
};

It is always better to use C++ built-in types wherever it's possible. That's why we use the 64-bit signed integer type to handle our input coordinates. int_x2_type and uint_x2_type is required to handle 96-bit signed/unsigned integers. As there is no such built-in type we use library provided efficient fixed integer type. The big integer type should be capable to handle 48 * 64 bit integers, that is less than 32 * 128, and so far we are good with extended_int<128> type. We use the same floating point type for both fpt_type and efpt_type as it has enough exponent bits to represent both 48 * 32 bit and 48 * 64 bit integers (that is also the reason we don't need two floating-point converter structures). The ulp_cmp_type structure checks weather two IEEE floating-point values are within the given signed integer ulp range (we won't analyze corresponding code here as it requires deep understanding of the floating-point architecture and its usage to compare floating-point values), but just to mention the declaration is following:

struct my_ulp_comparison {
  enum Result {
    LESS = -1,

    EQUAL = 0,
    MORE = 1
  };
  Result operator()(fpt80 a, fpt80 b, unsigned int maxUlps) const;
};

The last step would be to declare the my_fpt_converter structure (converts the integer types to the floating-point type):

struct my_fpt_converter {
  template <typename T>
  fpt80 operator()(const T& that) const {
    return static_cast<fpt80>(that);
  }

  template <size_t N>
  fpt80 operator()(const typename detail::extended_int<N>& that) const {
    fpt80 result = 0.0;
    for (size_t i = 1; i <= (std::min)((size_t)3, that.size()); ++i) {
      if (i != 1)
        result *= static_cast<fpt80>(0x100000000ULL);
      result += that.chunks()[that.size() - i];
    }
    return (that.count() < 0) ? -result : result;
  }
};

At this point we are done with declaration of the Voronoi coordinate type traits. The next step is to declare the Voronoi diagram traits:

struct my_voronoi_diagram_traits {
  typedef fpt80 coordinate_type;
  typedef voronoi_cell<coordinate_type> cell_type;
  typedef voronoi_vertex<coordinate_type> vertex_type;
  typedef voronoi_edge<coordinate_type> edge_type;
  struct vertex_equality_predicate_type {
  public:
    enum { ULPS = 128 };
    bool operator()(const point_type &v1, const point_type &v2) const {
      return (ulp_cmp(v1.x(), v2.x(), ULPS) == my_ulp_comparison::EQUAL &&
              ulp_cmp(v1.y(), v2.y(), ULPS) == my_ulp_comparison::EQUAL);
    }
  private:
    my_ulp_comparison ulp_cmp;
  };
};

Above we simply declared the Voronoi primitive types and vertex equality predicate using the new coordinate type and corresponding ulp comparison structure. As we are done with the declaration of the coordinate type specific structures we are ready to proceed to the construction step itself. The first step would be to initialize voronoi_builder structure with a set of random points:

// Random generator and distribution. MAX is equal to 2^48.
boost::mt19937_64 gen(std::time(0));
boost::random::uniform_int_distribution<boost::int64_t> distr(-MAX, MAX-1);

// Declaring and configuring Voronoi builder with the new coordinate type traits.
voronoi_builder<boost::int64_t, my_voronoi_ctype_traits> vb;


// Voronoi builder initialization step.
for (size_t i = 0; i < GENERATED_POINTS; ++i) {

  boost::int64_t x = distr(gen);
  boost::int64_t y = distr(gen);
  vb.insert_point(x, y);
}

The second step would be to generate the Voronoi diagram and this is done as before with the two lines of code:

// Declaring and configuring Voronoi diagram structure with the new coordinate type traits.
voronoi_diagram<fpt80, my_voronoi_diagram_traits> vd;

vb.construct(&vd);

From this point the user can operate with the Voronoi diagram data structure and in our tutorial we output some simple stats of the resulting Voronoi graph. Probably the hardest part of this tutorial is the declaration of the ulp comparison structure. The library provides efficient well-commented cross-platform implementation for 64-bit floating-point type (double). So the best advice would be to follow that implementation, but before doing that really consider thoughtfully if the default coordinate types are not capable to deal with your problem.

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)