DGtal
1.2.0

Part of the Geometry package.
This part of the manual describes the DGtal implementation of the famous "QuickHull" algorithm by Barber et al. [9] , and how to use it to compute convex hulls and Delaunay convex cell decompositions.
The following programs are related to this documentation: exampleLatticeBallQuickHull3D.cpp , exampleLatticeBallDelaunay2D.cpp , exampleRationalBallQuickHull3D.cpp , exampleRationalBallDelaunay3D.cpp , exampleQuickHull3D.cpp , testQuickHull.cpp , testConvexityHelper.cpp
The objective of the QuickHull algorithm is to compute the convex hull of a set of points V lying in a space of arbitrary dimension d (here d is greater than one). This algorithm has the limitation to only process full dimensional convex hulls, because of the way it is initialized. It maintains and updates a list of facets (which defines the faces of the current polytope) by following these steps:
qhull
implementation, this class uses a kernel that can be chosen in order to provide exact computations. This is the case for lattice points, but also rational points. Computation times are very similar.CGAL 3D convex hull
package, or with the arbitrary dimensional CGAL dD Triangulation
package, this algorithm does not build a simplicial convex hull. Facets may not be triangles or simplices in higher dimensions. This happens frequently for lattice points, where coplanar and cospherical situations are common.CGAL
function convex_hull_3
for the usual CGAL kernels Cartesian and Exact_predicates_inexact_constructions_kernel.In order to get a fine control over convex hull computations, you can use class QuickHull directly. It is required to choose a kernel specific to your objective. For now, the available kernels are:
CoordinateInteger
used for representing lattice points, and the integer type InternalInteger
used internally in facet and normal computations, as well as above predicates. Depending on your input (amplitude of coordinates and number of dimensions), you must choose consistently in order to get correct results. Default parameters are DGtal::int64_t for both.As a rule of the thumb, determinant computations raise integers to the power of the dimension, hence the choice of DGtal::BigInteger (requires WITH_GMP) for InternalInteger
is compulsory when \( N^d \ge 4.10^{18} \), for \( N \) the maximum norm of input points and \( d \) the dimension. When CoordinateInteger
is DGtal::int64_t, then the following table sums up these rules.
\( d \) dimension  max \( N \) for DGtal::int64_t  max \( N \) for DGtal::BigInteger 

2  \( \approx 2e9 \)  \( \approx 4e18 \) 
3  \( \approx 1.5e6 \)  \( \approx 4e18 \) 
4  \( \approx 4e4 \)  \( \approx 4e18 \) 
5  \( \approx 5e3 \)  \( \approx 4e18 \) 
6  \( \approx 1.2e3 \)  \( \approx 4e18 \) 
If you use rational kernels with a given precision \( p \), you should divide the above values by \( p \).
CoordinateInteger
and InternalInteger
, expect a 25 times slowdown factor if you use DGtal::BigInteger for InternalInteger
, and a 35 times slowdown factor if you use DGtal::BigInteger for both CoordinateInteger
and InternalInteger
.To compute the convex hull of lattice points, you need to include QuickHull.h
.
You should then define some typedefs:
Then we assume that V contains a range of lattice 3D points (with either 32 bits or 64 bits precision). You may compute their convex hull (facets and vertices) with (QuickHull::setInput and QuickHull::computeConvexHull):
You may check computation times with (QuickHull::timings)
You can easily retrieve the vertex positions with QuickHull::getVertexPositions (you can put them in a vector of real or lattice points), and get the consistently oriented vertices for each face with QuickHull::getFacetVertices. Below we show how to build a SurfaceMesh that represents the convex hull boundary and save it as OBJ file.
To give an idea of computation times, for 1e7 lattice points randomly chosen in a ball of radius 1000, computations times are as follows on a Macbook pro (processor 2,7 GHz QuadCore Intel Core i7, memory 16 GB 2133 MHz LPDDR3):
examples/geometry/tools/exampleLatticeBallQuickHull3D 10000000 1000 #points=9990515 #vertices=14183 #facets=28036 purge duplicates= 4221 ms. init simplex = 223 ms. quickhull core = 3214 ms. compute vertices= 255 ms. total time = 7913 ms.
QuickHull cannot handle exactly "real" points, but can handle rational approximations of such points. The user must specify a precision, the rational denominator, at kernel creation. To compute the convex hull of rational points, you need to include QuickHull.h
.
You should then define some typedefs:
Then we assume that V contains a range of real 3D points (with double precision). You may compute their approximate convex hull (facets and vertices) with the following lines (QuickHull::setInput and QuickHull::computeConvexHull):
The precision is the denominator used for representing rational coordinates. Real points with double
precision are rounded to these nearest rational coordinates.
You may check computation times with (QuickHull::timings)
You can easily retrieve the rational vertex positions with QuickHull::getVertexPositions (you can put them in a vector of real points), and get the consistently oriented vertices for each face with QuickHull::getFacetVertices. Below we show how to build a SurfaceMesh that represents the convex hull boundary and save it as OBJ file.
To give an idea of computation times, for 1e7 real points randomly chosen in a ball of radius 1000, computations times are as follows on a Macbook pro (processor 2,7 GHz QuadCore Intel Core i7, memory 16 GB 2133 MHz LPDDR3), for a precision of 1024:
examples/geometry/tools/exampleRationalBallQuickHull3D 10000000 1000 1024 #points=10000000 #vertices=14267 #facets=28503 purge duplicates= 4049 ms. init simplex = 322 ms. quickhull core = 3271 ms. compute vertices= 239 ms. total time = 7882 ms.
Computing the Delaunay complex of lattice points is very similar to computing its convex hull. You should just use the appropriate kernel as follows:
Then we assume that V contains a range of lattice 3D points (with either 32 bits or 64 bits precision). You may compute their Delaunay complex hull (facets and vertices) with (QuickHull::setInput and QuickHull::computeConvexHull):
You can easily retrieve the vertex positions with QuickHull::getVertexPositions (you can put them in a vector of real or lattice points), and get the consistently oriented vertices for each cell with QuickHull::getFacetVertices.
In opposition with convex hull computation, facets here designate the ddimensional Delaunay cells (here d=2). You also have finite facets (the ones within the convex hull) and infinite facets (the ones of the furthest site Delaunay decomposition). Below we show how to build a SurfaceMesh that represents the finite Delaunay cells and save it as OBJ file.
hull.nbFiniteFacets()
.To give an idea of computation times, for 1e6 lattice points randomly chosen in a ball of radius 1000, computations times are as follows on a Macbook pro (processor 2,7 GHz QuadCore Intel Core i7, memory 16 GB 2133 MHz LPDDR3):
examples/geometry/tools/exampleLatticeBallDelaunay2D 1000000 1000 #points=856202 #vertices=856202 #facets=1459906 purge duplicates= 284 ms. init simplex = 15 ms. quickhull core = 19316 ms. compute vertices= 1109 ms. total time = 20723 ms.
It is completely similar with the previous example, just replace the kernel by DelaunayRationalKernel and specify a precision at kernel instantiation.
A number of public data is stored in object QuickHull and is meaningfull after a convex hull computation:
To sum up, the maps are as follows and give indices in respective arrays:
input2comp p2v input points > processed points > convex hull vertices (user given) (no duplicates) (output) comp2input v2p input points < processed points < convex hull vertices (user given) (no duplicates) (output)
QuickHull:timings stores also the respective times taken by each step of the computation (see examples).
Class ConvexityHelper offers several functions that makes easier the computation of convex hull or Delaunay complex (of course with a little bit less flexibility than using QuickHull directly).
Building a lattice polytope in nD just requires a call to ConvexityHelper::computeLatticePolytope.
You may simply use one of the two static methods ConvexityHelper::computeConvexHullBoundary to build either a SurfaceMesh or a PolygonalSurface that represents the boundary convex hull of a set of 3D lattice points.
You can use ConvexityHelper::computeConvexHullCellComplex to build a cell complex that represents the convex hull of the given lattice points. This complex has exactly 1 full dimension cell, as many codimension 1 facets as the number of polytope faces, and as many vertices than the points lying on the convex hull boundary. w
You may use ConvexityHelper::computeDelaunayCellComplex to compute the Delaunay cell complex of a range of lattice points.
You can have a look at exampleLatticeBallDelaunay3D.cpp for a complete example, with a postprocessing that blowns up the 3D cells for visualization.
Building a rational polytope in nD that approximately encloses the given range of real points just requires a call to ConvexityHelper::computeRationalPolytope. You specify the precision as the denominator used in all rational numbers approximating the real point coordinates.
You may simply use one of the two static methods ConvexityHelper::computeConvexHullBoundary to build either a SurfaceMesh or a PolygonalSurface that represents the boundary convex hull of a set of 3D real points. The precision is given as a scale factor that transforms real coordinates before rounding them at the nearest integer.
You can use ConvexityHelper::computeConvexHullCellComplex to build a cell complex that represents the convex hull of the given real points. This complex has exactly 1 full dimension cell, as many codimension 1 facets as the number of polytope faces, and as many vertices than the points lying on the convex hull boundary.
You may use ConvexityHelper::computeDelaunayCellComplex to approximate the Delaunay cell complex of a range of real points, by specifying a precision.
You can have a look at exampleRationalBallDelaunay3D.cpp for a complete example, with a postprocessing that blowns up the 3D cells for visualization.