User Manual

Authors

Susan Hert and Stefan Schirra

Figure 13.1 The convex hull of a model made of 192135 points.

Introduction

A subset \( S \subseteq \mathbb{R}^3\) is convex if for any two points \( p\) and \( q\) in the set the line segment with endpoints \( p\) and \( q\) is contained in \( S\). The convex hull of a set \( S\) is the smallest convex set containing \( S\). The convex hull of a set of points \( P \in \mathbb{R}^3\) is a convex polytope with vertices in \( P\). A point in \( P\) is an extreme point (with respect to \( P\)) if it is a vertex of the convex hull of \( P\). A set of points is said to be strongly convex if it consists of only extreme points.

This chapter describes the functions provided in CGAL for producing convex hulls in three dimensions as well as functions for checking if sets of points are strongly convex are not. One can compute the convex hull of a set of points in three dimensions in two ways in CGAL: using a static algorithm or using a triangulation to get a fully dynamic computation.

Static Convex Hull Construction

The function convex_hull_3() provides an implementation of the quickhull algorithm [1]. There are two versions of this function available, one that can be used when it is known that the output will be a polyhedron (i.e., there are more than three points and they are not all collinear) and one that handles all degenerate cases and returns an Object, which may be a point, a segment, a triangle, or a polyhedron. Both versions accept a range of input iterators defining the set of points whose convex hull is to be computed and a traits class defining the geometric types and predicates used in computing the hull.

Traits Class

The function convex_hull_3() is parameterized by a traits class, which specifies the types and geometric primitives to be used in the computation. As the function constructs 3D planes from three input points, we cannot simply pass a kernel with inexact constructions as optional argument for the traits class.

If input points from a kernel with exact predicates and non-exact constructions are used, and a certified result is expected, the class Convex_hull_traits_3<R> should be used (R being the input kernel). If the constructions from a kernel are exact this kernel can be used directly as a traits class.

Note that the default traits class takes this into account, that is the above considerations are only important for custom traits classes.

Example

The following program reads points from an input file and computes their convex hull. We assume that the points are not all identical and not all collinear, thus we directly use a polyhedron as output. In the example you see that the convex hull function can write in any model of the concept MutableFaceGraph.

File Convex_hull_3/quickhull_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/convex_hull_3.h>
#include <vector>
#include <fstream>


typedef CGAL::Exact_predicates_inexact_constructions_kernel  K;
typedef CGAL::Polyhedron_3<K>                     Polyhedron_3;
typedef K::Point_3                                Point_3;
typedef CGAL::Surface_mesh<Point_3>               Surface_mesh;


int main(int argc, char* argv[])
{
  std::ifstream in( (argc>1)? argv[1] : "data/cube.xyz");
  std::vector<Point_3> points;
  Point_3 p;
  while(in >> p){
    points.push_back(p);
  }

  // define polyhedron to hold convex hull
  Polyhedron_3 poly;

  // compute convex hull of non-collinear points
  CGAL::convex_hull_3(points.begin(), points.end(), poly);

  std::cout << "The convex hull contains " << poly.size_of_vertices() << " vertices" << std::endl;

  Surface_mesh sm;
  CGAL::convex_hull_3(points.begin(), points.end(), sm);

  std::cout << "The convex hull contains " << num_vertices(sm) << " vertices" << std::endl;

  return 0;
}

Example for Lower Dimensional Results

The following program reads points from an input file and computes their convex hull. Depending on the dimension of the result, we will get a point, a segment, a triangle, or a polyhedral surface. Note that the latter may also be planar polygon with a border.

File Convex_hull_3/quickhull_any_dim_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/convex_hull_3.h>
#include <vector>
#include <fstream>

typedef CGAL::Exact_predicates_inexact_constructions_kernel  K;
typedef CGAL::Polyhedron_3<K>                     Polyhedron_3;
typedef K::Point_3                                Point_3;
typedef K::Segment_3                              Segment_3;
typedef K::Triangle_3                             Triangle_3;


int main(int argc, char* argv[])
{
  std::ifstream in( (argc>1)? argv[1] : "data/cube.xyz");
  std::vector<Point_3> points;
  Point_3 p;

  while(in >> p){
    points.push_back(p);
  }

  CGAL::Object obj;

  // compute convex hull of non-collinear points
  CGAL::convex_hull_3(points.begin(), points.end(), obj);

  if(const Point_3* p = CGAL::object_cast<Point_3>(&obj)){
    std::cout << "Point " << *p << std::endl;
  }
  else if(const Segment_3* s = CGAL::object_cast<Segment_3>(&obj)){
    std::cout << "Segment " << *s << std::endl;
  }
  else if(const Triangle_3* t = CGAL::object_cast<Triangle_3>(&obj)){
    std::cout << "Triangle " << *t << std::endl;
  }
  else  if(const Polyhedron_3* poly = CGAL::object_cast<Polyhedron_3>(&obj)){
    std::cout << "Polyhedron\n " << *poly << std::endl;
  std::cout << "The convex hull contains " << poly->size_of_vertices() << " vertices" << std::endl;

  }
  else {
    std::cout << "something else"<< std::endl;
  }
  return 0;
}

Extreme points

In addition to the convex_hull_3() function, the function extreme_points_3() is also provided in case only the points on the convex hull are required (without the connectivity information). In addition the traits class adapter CGAL::Extreme_points_traits_adapter_3 is also provided in order to get the indices or more generally any given entity that is associated a 3D point that is on the convex hull.

The following program reads a set of points from an OFF file and outputs the indices of the points that are on the convex hull.
File Convex_hull_3/extreme_indices_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/Extreme_points_traits_adapter_3.h>
#include <CGAL/IO/read_off_points.h>

#include <boost/iterator/counting_iterator.hpp>

#include <vector>
#include <fstream>


typedef CGAL::Exact_predicates_inexact_constructions_kernel      K;
typedef K::Point_3                                               Point_3;

int main(int argc, char* argv[])
{
  std::ifstream in( (argc>1)? argv[1] : "data/star.off");

  if(!in)
  {
    std::cerr<< "Cannot open input file." <<std::endl;
    return 1;
  }

  std::vector<Point_3> points;
  CGAL::read_off_points(in, std::back_inserter(points));

  //This will contain the extreme vertices
  std::vector<std::size_t> extreme_point_indices;

  //call the function with the traits adapter for vertices
  CGAL::extreme_points_3(CGAL::make_range(boost::counting_iterator<std::size_t>(0),
                                          boost::counting_iterator<std::size_t>(points.size())),
                         std::back_inserter(extreme_point_indices),
                         CGAL::make_extreme_points_traits_adapter(CGAL::make_property_map(points)));
  //print the number of extreme vertices
  std::cout << "Indices of points on the convex hull are:\n";
  std::copy(extreme_point_indices.begin(), extreme_point_indices.end(), std::ostream_iterator<std::size_t>(std::cout, " "));
  std::cout << "\n";

  return 0;
}

The following program reads and builds a mesh from an OFF file, and then collects the vertices that are on the convex hull of the mesh.
File Convex_hull_3/extreme_points_3_sm.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Extreme_points_traits_adapter_3.h>
#include <CGAL/convex_hull_3.h>
#include <vector>
#include <fstream>


typedef CGAL::Exact_predicates_inexact_constructions_kernel      K;
typedef K::Point_3                                               Point_3;
typedef CGAL::Surface_mesh<Point_3>                              Mesh;

int main(int argc, char* argv[])
{
  std::ifstream in( (argc>1)? argv[1] : "data/star.off");
  Mesh sm;
  if(! in || !(in >> sm) || !sm.is_valid()){
    std::cerr<<"Not a valid off file."<<std::endl;
  }

  //This will contain the extreme vertices
  std::vector<Mesh::Vertex_index> extreme_vertices;

  //call the function with the traits adapter for vertices
  CGAL::extreme_points_3(vertices(sm), std::back_inserter(extreme_vertices),
                         CGAL::make_extreme_points_traits_adapter(sm.points()));
  //print the number of extreme vertices
  std::cout << "There are  " << extreme_vertices.size() << " extreme vertices in this mesh." << std::endl;

  return 0;
}

Halfspace Intersection

The functions halfspace_intersection_3() and halfspace_intersection_with_constructions_3() uses the convex hull algorithm and the duality to compute the intersection of a list of halfspaces. The first version does not explicitly compute the dual points: the traits class handles this issue. The second one constructs these points and hence is less robust but the computation is faster.

In order to compute the intersection an interior point is needed. It can be either given by the user or computed using linear programming. Notice that the second approach is slower due to the resolution of a linear program.

Example

The following example computes the intersection of halfspaces defined by tangent planes to a sphere.

File Convex_hull_3/halfspace_intersection_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Convex_hull_3/dual/halfspace_intersection_3.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Surface_mesh.h>

#include <list>

typedef CGAL::Exact_predicates_inexact_constructions_kernel   K;
typedef K::Plane_3                                            Plane;
typedef K::Point_3                                            Point;
typedef CGAL::Surface_mesh<Point>                             Surface_mesh;


// compute the tangent plane of a point
template <typename K>
typename K::Plane_3 tangent_plane (typename K::Point_3 const& p) {
    typename K::Vector_3 v(p.x(), p.y(), p.z());
    v = v / sqrt(v.squared_length());
    typename K::Plane_3 plane(v.x(), v.y(), v.z(), -(p - CGAL::ORIGIN) * v);

    return plane;
}

int main (void) {
    // number of generated planes
    int N = 200;

    // generates random planes on a sphere
    std::list<Plane> planes;
    CGAL::Random_points_on_sphere_3<Point> g;
    for (int i = 0; i < N; i++) {
        planes.push_back(tangent_plane<K>(*g++));
    }

    // define polyhedron to hold the intersection
    Surface_mesh chull;

    // compute the intersection
    // if no point inside the intersection is provided, one
    // will be automatically found using linear programming
    CGAL::halfspace_intersection_3(planes.begin(),
                                   planes.end(),
                                   chull );

    std::cout << "The convex hull contains " << num_vertices(chull) << " vertices" << std::endl;

    return 0;
}

Convexity Checking

The function is_strongly_convex_3() implements the algorithm of Mehlhorn et al. [2] to determine if the vertices of a given polytope constitute a strongly convex point set or not. This function is used in postcondition testing for convex_hull_3().

Dynamic Convex Hull Construction

Fully dynamic maintenance of a convex hull can be achieved by using the class Delaunay_triangulation_3. This class supports insertion and removal of points (i.e., vertices of the triangulation) and the convex hull edges are simply the finite edges of infinite faces. The following example illustrates the dynamic construction of a convex hull. First, random points from a sphere of a certain radius are generated and are inserted into a triangulation. Then the number of points of the convex hull are obtained by counting the number of triangulation vertices incident to the infinite vertex. Some of the points are removed and then the number of points remaining on the hull are determined. Notice that the vertices incident to the infinite vertex of the triangulation are on the convex hull but it may be that not all of them are vertices of the hull.

Example

The following example shows how to compute a convex hull with a triangulation. The vertices incident to the infinite vertex are on the convex hull.

The function convex_hull_3_to_face_graph() can be used to obtain a polyhedral surface that is model of the concept MutableFaceGraph, e.g. Polyhedron_3 and Surface_mesh.

File Convex_hull_3/dynamic_hull_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/algorithm.h>
#include <CGAL/convex_hull_3_to_face_graph.h>

#include <list>

typedef CGAL::Exact_predicates_inexact_constructions_kernel     K;
typedef K::Point_3                                              Point_3;
typedef CGAL::Delaunay_triangulation_3<K>                       Delaunay;
typedef Delaunay::Vertex_handle                                 Vertex_handle;
typedef CGAL::Surface_mesh<Point_3>                             Surface_mesh;

int main()
{
  CGAL::Random_points_in_sphere_3<Point_3> gen(100.0);
  std::list<Point_3>   points;

  // generate 250 points randomly in a sphere of radius 100.0
  // and insert them into the triangulation
  std::copy_n(gen, 250, std::back_inserter(points) );
  Delaunay T;
  T.insert(points.begin(), points.end());

  std::list<Vertex_handle>  vertices;
  T.incident_vertices(T.infinite_vertex(), std::back_inserter(vertices));
  std::cout << "This convex hull of the 250 points has "
            << vertices.size() << " points on it." << std::endl;

  // remove 25 of the input points
  std::list<Vertex_handle>::iterator v_set_it = vertices.begin();
  for (int i = 0; i < 25; i++)
  {
     T.remove(*v_set_it);
     v_set_it++;
  }

  //copy the convex hull of points into a polyhedron and use it
  //to get the number of points on the convex hull
  Surface_mesh chull;
  CGAL::convex_hull_3_to_face_graph(T, chull);

  std::cout << "After removal of 25 points, there are "
            << num_vertices(chull) << " points on the convex hull." << std::endl;

  return 0;
}

Performance

In the following, we compare the running times of the two approaches to compute 3D convex hulls. For the static version (using convex_hull_3()) and the dynamic version (using Delaunay_triangulation_3 and convex_hull_3_to_face_graph()), the kernel used was Exact_predicates_inexact_constructions_kernel.

To compute the convex hull of a million of random points in a unit ball the static approach needed 1.63s, while the dynamic approach needed 9.50s. To compute the convex hull of the model of Figure 13.1 featuring 192135 points, the static approach needed 0.18s, while the dynamic approach needed 1.90s.

The measurements have been performed using CGAL 3.9, using the Gnu C++ compiler version 4.3.5, under Linux (Debian distribution), with the compilation options -O3 -DCGAL_NDEBUG. The computer used was equipped with a 64bit Intel Xeon 2.27GHz processor and 12GB of RAM.