- Author(s) of this documentation:
- Jacques-Olivier Lachaud
- Since
- 1.2
Part of the Geometry package.
This part of the manual describes the DGtal implementation of the famous "QuickHull" algorithm by Barber et al. [7] , 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 , testGenericLatticeConvexHull.cpp , exampleGenericLatticeConvexHull3D.cpp , exampleGenericLatticeConvexHull4D.cpp
- Note
- Since version 2.0, QuickHull detects the affine dimension of input sets of points exactly. You may also use class GenericLatticeConvexHull to compute the convex hull of non full dimensional sets. See Generic convex hulls
The QuickHull convex hull algorithm
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:
- Initialization: an initial d-dimensional simplex is built with \( d+1 \) facets. Furthermore, each remaining point of V is associated to one of these facets if it is above it. Hence, each point of V is either associated to no facet and is thus already inside the current hull, or it is associated to only one facet. Points on each facet and neighbors of facets are computed (easy for a simplex) and all the facets are then queued.
- While the queue of facets is not empty, the front facet f is popped.
- ignore it if contains no "above" point(s)
- otherwise pick the furthest above point p
- find by neighbor adjacencies from f all the facets visible from p (i.e. p is above them)
- define the horizon H as the set of ridges R, where a ridge is a pair of two neighbor facets, where the first is visible from p, whilst the second is not.
- define the new facets as the ones formed by the umbrella from p to its ridges.
- merge possible parallel facets and update "on" points
- reassign "above" points of f to these new facets and put them into the queue
- The algorithm stops when all current facets have no "above" points.
- Optionally it extracts vertices of the convex hull polytope by determining the "on" points that belong to at least d facets.
- Note
- In opposition to the usual 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.
- In opposition to 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.
- This version is generally more than twice faster than CGAL function convex_hull_3 for the usual CGAL kernels Cartesian and Exact_predicates_inexact_constructions_kernel.
- However this implementation is not tailored for incremental dynamic convex hull computations.
- Last, the method for computing the Delaunay triangulation by using the convex hull in a higher dimensional space is generally much slower than CGAL method for computing the Delaunay triangulation by point location (5-10 times slower in our experience).
Convex hull and Delaunay cell complex
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:
- ConvexHullIntegralKernel : it allows the computation of the convex hull of a range of lattice points in arbitrary dimension.
- ConvexHullRationalKernel : it allows the computation of the convex hull of a range of rational points in arbitrary dimension. The user fixes the precision at kernel instantiation.
- DelaunayIntegralKernel : it allows the computation of the Delaunay cell complex of a range of lattice points in arbitrary dimension.
- DelaunayRationalKernel : it allows the computation of the Delaunay cell complex of a range of rational points in arbitrary dimension. The user fixes the precision at kernel instantiation.
- Note
- All kernels are parametrized by the dimension of input points, the integer type 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 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 \).
- Note
- With respect to the choice DGtal::int64_t for both CoordinateInteger and InternalInteger, expect a 25 times slow-down factor if you use DGtal::BigInteger for InternalInteger, and a 35 times slow-down factor if you use DGtal::BigInteger for both CoordinateInteger and InternalInteger.
Convex hull of lattice points
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.
Convex hull of lattice ball with radius 12.5 | 
Convex hull of lattice ball with radius 25 | 
Convex hull of lattice ball with radius 50 |
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 Quad-Core 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.
- Note
- Half of the time is spent on removing duplicated input points. If you know that your range of points does not contain duplicates, you can specify it to QuickHull::setInput.
Convex hull of rational points
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.
- Note
- Different real points may be rounded to the same rational point.
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.
Convex hull of rational ball with radius 10.5 with precision 1 | 
Convex hull of rational ball with radius 10.5 with precision 2 | 
Convex hull of rational ball with radius 10.5 with precision 4 |
Convex hull of rational ball with radius 10.5 with precision 16 | 
Convex hull of rational ball with radius 10.5 with precision 128 | 
Convex hull of rational ball with radius 10.5 with precision 1024 |
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 Quad-Core 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.
- Note
- Half of the time is spent on removing duplicated input points. If you know that your range of points does not contain duplicates, you can specify it to QuickHull::setInput.
- If you choose a too big precision, integer computations required by convex hull computations might overflow. You may consider using DGtal::BigInteger at least for the internal integer representation.
Delaunay cell complex of lattice points (2D)
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 d-dimensional 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.
- Note
- If you wish to visit the infinite facets, you may start from facet numbered hull.nbFiniteFacets().
Delaunay cell decomposition of 100 randomly chosen points in a lattice ball with radius 20 | 
Delaunay cell decomposition of 1000 randomly chosen points in a lattice ball with radius 20 |
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 Quad-Core 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.
Delaunay cell complex of rational points (2D)
It is completely similar with the previous example, just replace the kernel by DelaunayRationalKernel and specify a precision at kernel instantiation.
Further information
A number of public data is stored in object QuickHull and is meaningfull after a convex hull computation:
- QuickHull::input2comp the surjective mapping giving for each input point index its index in the point range processed by quick hull
- QuickHull::comp2input the injective mapping giving for each processed point index its index in the input point range.
- QuickHull::p2v the map giving for each processed point index its index in the range of vertices (or UNASSIGNED if it is not an index)
- QuickHull::v2p the map giving for each vertex index its index in the range of processed points.
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).
You may inquire about the current state of QuickHull object with QuickHull::status. Possible states are
- QuickHull::Status::Uninitialized QuickHull is empty and has just been instantiated.
- QuickHull::Status::InputInitialized A range of input points has been given to QuickHull.
- QuickHull::Status::SimplexCompleted An initial full-dimensional simplex has been found. QuickHull core algorithm can start.
- QuickHull::Status::FacetsCompleted All facets of the convex hull are identified.
- QuickHull::Status::VerticesCompleted All vertices of the convex hull are determined.
- QuickHull::Status::AllCompleted Same as VerticesCompleted.
- QuickHull::Status::NotFullDimensional Error: the initial simplex is not full dimensional.
- QuickHull::Status::InvalidRidge Error: some ridge is not consistent (probably integer overflow).
- QuickHull::Status::InvalidConvexHull Error: the convex hull does not contain all its points (probably integer overflow).
Use ConvexityHelper as a shortcut
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).
Build a lattice polytope (nD)
Building a lattice polytope in nD just requires a call to ConvexityHelper::computeLatticePolytope.
- Note
- For dimension up to 3D, you may indicate if you wish your lattice polytope to be Minkowski summable by a cubical grid cell (see ConvexityHelper::computeLatticePolytope reference doc).
Get boundary surface of lattice convex hull (3D)
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.
Get convex hull cell complex of lattice points (nD)
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
Get Delaunay cell complex of lattice points (nD)
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 post-processing that blowns up the 3D cells for visualization.
Delaunay cell decomposition of randomly chosen points in a 3D lattice ball with radius 5 | 
Delaunay cell decomposition of randomly chosen points in a 3D lattice ball with radius 10 |
Build a rational polytope (nD)
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.
- See also
- exampleQuickHull3D.cpp
- Note
- In opposition with ConvexityHelper::computeConvexHullBoundary the precision must be an integer (no smaller than one), but it acts similarly. This is because the built polytope must be a rational polytope (i.e. constraints are rational).
Get boundary surface of rational convex hull (3D)
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.
- See also
- exampleQuickHull3D.cpp
- Note
- In opposition with ConvexityHelper::computeRationalPolytope the precision might be any double value, but it acts similarly. This is because the build object is a V-representation of a polytope, with floating-point point coordinates.
Coarse model of Stanford bunny | 
Its convex hull with a precision 100 |
Get convex hull cell complex of real points (nD)
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.
- See also
- exampleQuickHull3D.cpp
Get Delaunay cell complex of real points (nD)
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 post-processing that blowns up the 3D cells for visualization.
Generic convex hulls
Since 2.0, you can compute the convex hull of any input set of point X in any dimension, whatever the affine dimension of X. The core class is GenericLatticeConvexHull, parameterized by the dimension of input points, the coordinate type for points and a type for internal computations and output.
Its use is quite simple. The following example draws at random 100000 sets of 3-9 points in 3D, and computes the proportions of the affine dimension of their convex hull as their average number of vertices.
It outputs
0.024% are 0-dimensional #V=1 #F=0 (24/24) 11.2% are 1-dimensional #V=2 #F=0 (11208/11208) 11.3% are 2-dimensional #V=3.02 #F=3.02 (11255/11255) 77.5% are 3-dimensional #V=6.41 #F=8.62 (77513/77513)
Essentially, you call method GenericLatticeConvexHull::compute with an input set of points and the instantiated GenericLatticeConvexHull object holds the convex hull in its following public data members:
- GenericLatticeConvexHull::affine_dimension The affine dimension of the input set.
- GenericLatticeConvexHull::positions The positions of the verticesfacets (a subset of the input points).
- GenericLatticeConvexHull::facets The range giving for each facet the indices of its vertices.
- GenericLatticeConvexHull::points The set of input points, indexed as in the input, possibly converted into the chosen coordinate type.
- GenericLatticeConvexHull::vertex2point The indices of the vertices of the convex hull in the input set.
- GenericLatticeConvexHull::projected_points The set of projected input points, indexed as in the input, converted into the chosen coordinate type. If the affine dimension k is equal to the ambient dimension n, then it is equal to GenericLatticeConvexHull::points. Otherwise the last \( n-k \) coordinates are null).
- GenericLatticeConvexHull::affine_basis The affine basis used for projection (identity if the convex hull is full dimensional, otherwise a basis in reduced echelon form).
You may yave a look at exampleGenericLatticeConvexHull3D.cpp and exampleGenericLatticeConvexHull4D.cpp for examples of use of this class.
Convex hull of 10 points of affine dimension 1 | 
Convex hull of 100 points of affine dimension 2 |
Convex hull of 1000 points of affine dimension 3 | 
Convex hull of 1000 points of affine dimension 3 (interior view) |
You may also count the number of lattice points within the convex polytope spanned by the computed convex hull.
- GenericLatticeConvexHull::count() Returns the number of integer points lying within the convex hull polytope.
- GenericLatticeConvexHull::countInterior() Returns the number of integer points lying within the interior of the convex hull polytope.
- GenericLatticeConvexHull::countBoundary() Returns the number of integer points lying on the convex hull boundary of the polytope.
- GenericLatticeConvexHull::countUpTo Returns the number of integer points within the polytope up to some maximum number max. For instance, a d-dimensional simplex that contains no integer points in its interior contains only d+1 points. If there is more, you know that the simplex has a non empty interior.
- Warning
- Counting results are valid only if GenericLatticeConvexHull::proj_dilation is equal to 1 (i.e. the affine basis for the projection was obtained through an unimodular transformation). Otherwise the result is greater or equal to the expected number.
Further notes
Since 2.0, the QuickHull and ConvexityHelper classes have been slightly improved:
- The determination of the full dimensionality of input points is now exact: it keeps a randomized initial approach to try to get an as big as possible initial simplex, but if no simplex was found, it uses AffineGeometry to determine the exact affine dimension of the input.
ConvexityHelper handles a more general family of degenerated convex hulls. ConvexityHelper::computeLatticePolytope and Computelatticepolytope can now compute convex hulls of:
- 0D and 1D affine dimension in dD
- d-1D and dD affine dimension in dD
Before this version, only 2D affine dimension in 3D was possible.
method ConvexHullCommonKernel::compute( const std::vector< CoordinatePoint >& vpoints, const CombinatorialPlaneSimplex& simplex ), whose aim is to compute a vector orthogonal to the affine space spanned by a given d-simplex, uses now AffineGeometry::orthogonalVector instead of cofactor computations with SimpleMatrix.
This approach is quicker when the dimension d increases (with a complexity upper bounded by \(O(d^3)\)), while keeping the same running times as before in 2D and 3D.
For instance, in 6D, quickhull is almost twice faster than before.
