Close to the end of my GSoC project, it is time to present the product of my work. As a quick reminder, I based my project on project proposal 3 and project proposal 4 for Boost.Geometry from the Boost GSoC 2019 Wiki. The goals of my project were to implement

  • robust geometric predicates for Delaunay triangulation
  • a triangulation data structure
  • a triangle data structure
  • a Voronoi diagram data structure
  • a Delaunay triangulation algorithm for 2D points in a cartesian coordinate system
  • uniform point distributions for various types of domains

Work Product

The entirety of the work product can be found in the develop branch of the BooostGSoC19/geometry repository, specifically in the commits of the following Github comparison: https://github.com/BoostGSoC19/geometry/compare/5679ccd08fa5bd09c424a6a7878ad0dafc2c8f93…BoostGSoC19:bad1b9c5a2f4f458284a912a848a25e73c28014b (Updated on 2019-08-25, was: https://github.com/BoostGSoC19/geometry/compare/5679ccd08fa5bd09c424a6a7878ad0dafc2c8f93…BoostGSoC19:4c1cf99ad9f20d85afe30b170501ab625a15b11b). For reference, here is a full list of the commits that make up the work product:

0704d7e660f9faa57f4c68eb1bd95cebb921b704

0c4615575f9b10a4f522cb1096814f20f7a5a8b8

57632e189555f1089fc577d78492dae8ccc9c5d4

a6bfe950c9e9ce4718abe1543316f91c80e1fe07

d906785184d48c84af6625f6960ecfdfc11abacf

b9a3660f47e8ddf818b28add26695a5b84102e12

645f6ed4faf0ce72aba53ff42fcce05d2279c3ef

cff87f3f28ed1088aece0860fa2e5bfc6e347a36

d140ddbb7c7f57f67fa5f7a5080651563bec80c7

74cad8410731a46c41b556563bbac38bd2095d8d

1f55c44bc9363e0c8b1b7c4ff1cbfe695d7dbec7

f150cc1cb84d58458a8e958623d8806d52efb349

1156767e5a8d1243a682270da17b4fae8784dbae

4c1cf99ad9f20d85afe30b170501ab625a15b11b

(Update on 2019-08-25, the following new commits were added)

e981c4656edeb446f09bfb863e87d29499da67f9

361861ca9fb37e37aad36f81bfbee84b686bd3e1

bad1b9c5a2f4f458284a912a848a25e73c28014b

For reference, I will also include a list of the files containing the work product:

extensions/example/Jamfile.v2

extensions/example/random/random_example.cpp

extensions/example/triangulation/Jamfile.v2

extensions/example/triangulation/triangulation_example.cpp

extensions/test/Jamfile.v2

extensions/test/random/random.cpp

extensions/test/triangulation/Jamfile.v2

extensions/test/triangulation/in_circle_robust.cpp

extensions/test/triangulation/side_robust.cpp

extensions/test/triangulation/triangulation.cpp

include/boost/geometry/extensions/random/detail/uniform_point_distribution.hpp

include/boost/geometry/extensions/random/dispatch/uniform_point_distribution.hpp

include/boost/geometry/extensions/random/strategies/agnostic/uniform_envelope_rejection.hpp

include/boost/geometry/extensions/random/strategies/agnostic/uniform_linear.hpp

include/boost/geometry/extensions/random/strategies/agnostic/uniform_point_distribution_discrete.hpp

include/boost/geometry/extensions/random/strategies/cartesian/uniform_point_distribution_box.hpp

include/boost/geometry/extensions/random/strategies/cartesian/uniform_point_distribution_segment.hpp

include/boost/geometry/extensions/random/strategies/spherical/edwilliams_avform_intermediate.hpp

include/boost/geometry/extensions/random/strategies/spherical/uniform_inverse_transform_sampling.hpp

include/boost/geometry/extensions/random/strategies/uniform_point_distribution.hpp

include/boost/geometry/extensions/random/uniform_point_distribution.hpp

include/boost/geometry/extensions/triangulation/algorithms/delaunay_triangulation.hpp

include/boost/geometry/extensions/triangulation/geometries/triangulation.hpp

include/boost/geometry/extensions/triangulation/geometries/voronoi_adaptor.hpp

include/boost/geometry/extensions/triangulation/strategies/cartesian/accelerated_shull.hpp

include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/accelerated_shull.hpp

include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/precise_math.hpp

include/boost/geometry/extensions/triangulation/strategies/cartesian/in_circle_robust.hpp

include/boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp

include/boost/geometry/extensions/triangulation/strategies/delaunay_triangulation.hpp

include/boost/geometry/extensions/triangulation/triangulation.hpp

(Update on 2019-08-25, the following new files were added)

include/boost/geometry/extensions/random/strategies/agnostic/uniform_convex_fan.hpp

include/boost/geometry/extensions/random/strategies/agnostic/uniform_convex_hull_rejection.hpp

include/boost/geometry/extensions/random/strategies/cartesian/uniform_point_distribution_triangle.hpp

Except for two Jamfiles, these files were all created from scratch during the GSoC.

The above comparison contains the entire work product. Because of some branch reorganization, it does not reflect the complete history of the GSoC work and the time of the creation of most of the linked code. If the entire code history is of interest, it can be found here. The correct dates of the individual commits can be found by clicking on the commit message.

In the following passages, I will elaborate on the individual project goals and the work products in greater detail.

Robust geometric predicates

Motivation

The construction of a Delaunay triangulation usually (depending on the algorithm used) involves the repeated evaluation of two geometric predicates, namely the orientation check and the in-circle check. The orientation check determines on which side of a segment, given by two points, a third point lies and, consequently, whether the three given points constitute a valid, positively oriented triangle. The in-circle check determines whether the circumcircle of a given triangle contains the fourth given point and consequently, whether a triangulation containing this triangle as a face and the fourth point as a vertex violates the Delaunay property.

A straightforward way to perform the orientation check is to compute the signed triangle area of the given three points. A triangle is valid and positively oriented if and only if its signed area is positive. We can express the signed areas of a triangle in determinant form. Boost.Geometry already includes a strategy that implements this method of checking orientation in the side_by_triangle strategy. Similarly, although not yet implemented as a strategy in Boost.Geometry, the in-circle check can also be performed by evaluating the sign of a determinant (see here).

However, due to the limitations of floating-point arithmetics, especially cancellation, the direct computation of these determinants may produce incorrect results. Such false results may cause the construction of a Delaunay triangulation to break down; therefore, they motivate the implementation of more robust geometric predicates.

Design

There are various ways to work around the limited precision of floating-point types, some of which, such as using exact number types for all calculations or extended precision incur heavy performance penalties. I have instead chosen to implement the adaptive, robust predicates described by Jonathan Richard Shewchuk in his article Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates (paper and more information is freely available at www.cs.cmu.edu/~quake/robust.html at the time of writing this).

The basic idea is to compute non-overlapping expansions of floating-point computation results that contain the bits of the result that would be eliminated by the rounding aspect of direct floating-point arithmetics. The full expansion of a calculation then represents its result at full precision, which makes the computation robust. The implementation is made adaptive by computing the parts of the expansions of interim results in decreasing order of their significance in the final result. At various stages of the computation, the incomplete result is checked against error bounds and, if the sign of the result is already guaranteed, the calculation is terminated early. This means that, for most inputs, the amount of work done is not much higher than with standard floating-point arithmetics.

I have chosen to implement the robust predicates as described by Shewchuk because they are fast and well-described such that they were easy to implement.

Work Product

I have implemented the methods from the linked paper that were used in the orient2d and incircle functions. The actual implementation of the calculations can be found in include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/precise_math.hpp. Comparing them to the original C implementation, allows me to highlight some design choices that were enabled by the use of modern C++.

Instead of computing epsilon and splitter at run-time and storing them in global variables, my implementation does not involve any global variables and instead computes epsilon and splitter at compile-time (depending on constexpr-support) using the numerical_limits header from the STL.

Rather than macros for smaller math functions, I have used templates, which delegates micro-optimization related decisions about inlining to the compiler and allows better scoping of temporary variables. Rather than passing points, I was able to make use of the std::array class, which makes array bounds explicit throughout the whole code.

Because I used templates, it was easy to add a compile-time parameter that allows the user of the implementation to specify the maximum precision.

Having the mathematics implemented in the detail-namespace, I then added strategies in two classes include/boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp

include/boost/geometry/extensions/triangulation/strategies/cartesian/in_circle_robust.hpp

These classes also contain the documentation of my work. Finally, I have implemented two test cases for the orientation test and for the in-circle test with numbers that would produce incorrect results in a non-robust implementation.

Triangulation, Triangle and Voronoi data structures

Motivation

Delaunay Triangulations and Voronoi diagrams have many applications. For example, they are used in Finite Element Methods or Finite Volume Methods. Because my project involved an algorithm for Delaunay triangulation and Boost.Geometry does not have a triangulation data structure, it was necessary to implement one. A triangle data structure was also planned, but during the work, it became clear that it was sufficient for now and much easier to make the triangles of my triangulation conform to the ring concept.

Design

I implemented the triangulation as a list of vertices and faces with references to each other. Because the Voronoi diagram is the dual of the Delaunay triangulation, I found that it was sufficient, to implement adapter classes that allow viewing Triangulation vertices as Voronoid Diagram faces and viewing Triangulation faces as Voronoi Vertices.

Work Product

My implementation of a triangulation data structure can be found in include/boost/geometry/extensions/triangulation/geometries/triangulation.hpp. I will again highlight some design choices:

  • The vertex struct contains a point and a face-iterator to one of its adjacent faces. All necessary traits are declared so that vertex instances can be directly treated as points by other Boost.Geometry algorithms. Holding a face-iterator rather than a face-index allows iterating over neighboring vertices and faces without passing a reference to the triangulation itself.
  • The face struct contains an array of three vertex-iterators. Using boost:indirect_iterator the face can be directly treated as a range of vertices so that it can be made to conform to the ring concept. Other than that, it stores iterators to its three neighboring faces and additionally the indices for half-edge reversal.
  • The triangulation requires a point type as a template parameter. It optionally allows specifying the triangle orientation and the type of containers that are used, as well as their allocators. By default, triangles are oriented counterclockwise, and vectors are used to store faces and vertices. Since by construction, iterator stability is necessary, the number of vertices must be provided in advance if vectors are used. If that is undesirable, a user may opt to use container types with better iterator stability.
  • Several invariants are maintained to simplify navigation over the triangulation.
  • The i-th neighbor of each triangle is the face opposite to its i-th vertex. This simplifies flip operations.
  • The triangulation always holds an iterator to a boundary vertex to simplify iteration over the boundary edges.
  • If a vertex lies on the boundary, the iterator stored to an incident face will point to a face with two boundary vertices. This simplifies the stop-criterion for the iteration over its incident faces or its neighboring vertices.
  • Several free functions that return ranges are provided, such us, for vertices, the range of incident faces and incident vertices, and, for faces, the range of adjacent faces, incident faces, and vertices.

My implementation of the Voronoi diagram data structure can be found in include/boost/geometry/extensions/triangulation/geometries/voronoi_adaptor.hpp

It consists of two classes

  • The voronoi_vertex_view is constructed from a triangulation face and allows to access it as the corresponding Voronoi vertex (which is equal to the circumcircle center of the face.
  • The voronoi_face_view is constructed from a triangulation vertex.

Delaunay Triangulation

Design

I adapted the S-Hull algorithm as described in S-hull: a fast radial sweep-hull routine for Delaunay triangulation by David Sinclair. To accelerate the construction, I changed step 7 to perform a binary search over the convex hull, rather than a linear search. This optimization greatly reduces the number of calls to the orientation predicate.

Work Product

The function that is used to call the construction is implemented in include/boost/geometry/extensions/triangulation/algorithms/delaunay_triangulation.hpp, but the actual implementation has been delegated to a construction strategy because it is coordinate-system and dimension-specific and more construction strategies may be added in the future. The construction has been implemented in include/boost/geometry/extensions/triangulation/strategies/cartesian/accelerated_shull.hpp.

Additionally, I produced a unit test in extensions/test/triangulation/triangulation.cpp that produces a small triangulation and an example in extensions/example/triangulation/triangulation_example.cpp that demonstrates basic usage of the triangulation and the Voronoi adapter. The output file of the example is the following graphic: The blue triangles are the faces of the Delaunay triangulation and the red lines are the Voronoi edges excluding infinite edges.

Random Point Distributions

Design

Since there was no concept of geometry distributions in Boost.Geometry before, I decided to model it after the existing concept of number distributions as found in Boost.Random or the STL. In place of the min() and max(), I have put a domain() function because the former functions are not applicable in a multidimensional space. Other than that, I think I was able to adapt the existing distribution concept to geometries.

Work Product

The general interface of uniform point distributions can be found in include/boost/geometry/extensions/random/uniform_point_distribution.hpp. To stay true to the strategy-paradigm of Boost.Geometry, I implemented the CS- and Geometry-specific code in strategies.

After implementing the strategies, I added a unit test in extensions/test/random/random.cpp. The test checks whether

  • point distributions are correctly recovered from operator<< and operator>>,
  • the generated points lie in the domain,
  • the results are roughly uniformly distributed.

Finally, I’ve added an example for Random Point Geometries. The example shows how to generate points in various cartesian geometries as well as one instance of spherical distributions. Its output file is the following graphic: From top to bottom and left to right, the drawn images represent points from instances of uniform_point_distribution with the following domains:

  1. A box (0° 0°, 90° 90°) with spherical coordinates. The points are drawn by directly using longitude and latitude as x- and y-coordinates. It should be noted that this representation is deceiving with regard to the uniformity of the points because it makes points around the north pole (at the top of the image) look further apart than they actually are.
  2. A multipoint with cartesian coordinates.
  3. A polygon with cartesian coordinates.
  4. A 2D box with cartesian coordinates.
  5. A 2D segment with cartesian coordinates.

(Update 2019-08-26 19:50 CEST) The following two images show two more interesting usages. The first is a Voronoi diagram of German airports. The second is an image of 1,000 uniformy distributed points sampled in the boundaries of Germany.

All SVG files in this blog post are licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license. The airport data in the Voronoi Diagram example is derived from the airports.dat from openflights.org.