3D Triangulation

The basic 3D-triangulation class of CGAL is primarily designed to represent the triangulations of a set of points $$A in $$ ³. It is a partition of the convex hull of $$A into tetrahedra whose vertices are the points of $$A. Together with the unbounded cell having the convex hull boundary as its frontier, the triangulation forms a partition of $$ ³. Its cells ($$3-faces) are such that two cells either do not intersect or share a common facet ($$2-face), edge ($$1-face) or vertex ($$0-face).

Representation

In order to deal only with tetrahedra, which is convenient for many applications, the unbounded cell can be subdivided into tetrahedra by considering that each convex hull facet is incident to an infinite cell having as fourth vertex an auxiliary vertex called the infinite vertex. In that way, each facet is incident to exactly two cells and special cases at the boundary of the convex hull are simple to deal with.

The class Triangulation_3<TriangulationTraits_3,TriangulationDataStructure_3> of CGAL implements this point of view and therefore considers the triangulation of the set of points as a set of finite and infinite tetrahedra. Notice that the infinite vertex has no significant coordinates and that no geometric predicate can be applied on it.

A triangulation is a collection of vertices and cells that are linked together through incidence and adjacency relations. Each cell gives access to its four incident vertices and to its four adjacent cells. Each vertex gives access to one of its incident cells.

The four vertices of a cell are indexed with 0, 1, 2 and 3 in positive orientation, the positive orientation being defined by the orientation of the underlying Euclidean space $$ ³. The neighbors of a cell are also indexed with 0, 1, 2, 3 in such a way that the neighbor indexed by $$i is opposite to the vertex with the same index. See Figure .

Figure:  Orientation of a cell (3-dimensional case).

As in the underlying combinatorial triangulation (see Chapter ), edges ($$1-faces) and facets ($$2-faces) are not explicitly represented: a facet is given by a cell and an index (the facet i of a cell c is the facet of c that is opposite to the vertex with index i) and an edge is given by a cell and two indices (the edge (i,j) of a cell c is the edge whose endpoints are the vertices of c with indices i and j). See Figure .

Degenerate Dimensions  The class Triangulation_3<TriangulationTraits_3,TriangulationDataStructure_3> can also deal with triangulations whose dimension is less than 3. A triangulation of a set of points in $$ ^d covers the whole space $$ ^d and consists of cells having $$d+1 vertices: some of them are infinite, they are obtained by linking the additional infinite vertex to each facet of the convex hull of the points.

• dimension 2: when a triangulation only contains coplanar points (which is the case when there are only three points), it consists of triangular faces.
• dimension 1: the triangulation contains only collinear points (which is the case when there are only two points), it consists of edges.
• dimension 0: the triangulation contains only one finite point.
• dimension -1: this is a convention to handle the case when the only vertex of the triangulation is the infinite one.

The same cell class is used in all cases: triangular faces in 2D can be considered as degenerate cells, having only three vertices (resp. neighbors) numbered $$(0,1,2), and one $$NULL vertex (resp. neighbor); edges in 1D have only two vertices (resp. neighbors) numbered $$0 and $$1.

The implicit representation of facets (resp. edges) still holds for degenerate dimensions (i.e.dimensions $$<3) : in dimension 2, each cell has only one facet of index 3, and 3 edges $$(0,1), $$(1,2) and $$(2,0); in dimension 1, each cell has one edge $$(0,1).

Validity  A triangulation of $$ ³ is said to be locally valid iff

(a)-(b) Its underlying combinatorial graph, the triangulation data structure, is locally valid (see Section  of Chapter )
(c) Any cell has its vertices ordered according to positive orientation. See Figure .

When the triangulation is degenerated into a triangulation of dimension 2, the geometric validity reduces to:

(c-2D) For any two adjacent triangles $$(u,v,w₁) and $$(u,v,w₂) with common edge $$(u,v), $$w₁ and $$w₂ lie on opposite sides of $$(u,v) in the plane.

When all the points are collinear, this condition becomes:

(c-1D) For any two adjacent edges $$(u,v) and $$(v,w), $$u and $$w lie on opposite sides of the common vertex $$v on the line.

The is_valid() method provided in triangulation_3 checks the local validity of a given triangulation. This does not always ensure global validity [MNS$$⁺96, DLPT98] but it is sufficient for practical cases.

Software Design

The class Triangulation_3<TriangulationTraits_3,TriangulationDataStructure_3> is designed to be used as a layer upon a 3D-triangulation data structure as presented in Section  of Chapter . It provides high-level geometric operations such as location of a point in the triangulation and insertion of a point, and is responsible for the geometric validity. This class is parameterized by two classes:

• the geometric traits class, where the user can specify the type of points to use as well as the elementary operations on them (predicates,...). The concept of such a class is introduced in Section  and the homegeneous and cartesian kernels of CGAL are models for this concept (see also CGAL::Regular_triangulation_euclidean_traits_3<R,Weight> for the Regular triangulation).
• the triangulation data structure class of the middle level, described in Chapter .

Delaunay triangulations as well as triangulation hierarchies [Dev98] are also implemented in the package: Delaunay_triangulation_3<DelaunayTriangulationTraits_3,TriangulationDataStructure_3> inherits from Triangulation_3<DelaunayTriangulationTraits_3,TriangulationDataStructure_3>, and Triangulation_hierarchy_3<Tr> inherits from Tr.

Triangulation_3<TriangulationTraits_3,TriangulationDataStructure_3> derives from Triangulation_utils_3<TriangulationTraits_3,TriangulationDataStructure_3>, which defines a set of tools working on the indices of vertices in cells.

The Geometric Traits

The first template parameter of the triangulation class Triangulation_3<TriangulationTraits_3,TriangulationDataStructure_3> of CGAL is the geometric traits class.

The geometric traits class must define the geometric objects (points, segments, triangles and tetrahedra) forming the triangulation together with a few geometric predicates on these objects: equality, coordinate comparison, orientation in space, orientation in case of coplanar points, and collinearity tests.

In addition to the requirements described before, the geometric traits class of a Delaunay triangulation must define predicates for the empty sphere property.

CGAL provides two predefined geometric traits classes Cartesian<R> and Homogeneous<R> using the kernel geometric objects and predicates. These classes are templated with a representation class. They supply the user with all the functionalities described for the concepts TriangulationTraits_3, DelaunayTriangulationTraits_3 and TriangulationHierarchyTraits_3. So, any of them can be used as a default traits class for Triangulation_3<TriangulationTraits_3,TriangulationDataStructure_3>, Delaunay_triangulation_3<DelaunayTriangulationTraits_3,TriangulationDataStructure_3> and for the parameter Tr of TriangulationHierarchy_3<Tr>.

To be used as the traits class for the regular triangulation provided by CGAL, a class must provide definitions for the power tests. (See Section .) Regular_triangulation_euclidean_traits_3<R,Weight> is a traits class designed to be used by the class Regular_triangulation_3<RegularTriangulationTraits_3,TriangulationDataStructure_3>. It provides Weighted_point, a class for weighted points provided by the class, which derives from the CGAL three dimensional point class. It supplies the user with all the functionalities described for the concept RegularTriangulationTraits_3. It can be used as a default traits class for Regular_triangulation_3<RegularTriangulationTraits_3,TriangulationDataStructure_3>.

The Triangulation Data Structure Parameter

The second template parameter of the basic triangulation class Triangulation_3<TriangulationTraits_3,TriangulationDataStructure_3> is a triangulation data structure class. This class can be seen as a container for the cells and vertices maintaining incidence and adjacency relations (see Chapter ). A model of this triangulation data structure is Triangulation_data_structure_3.

A default parameter is defined in all the triangulation classes, so, it need not be specified by the user.

Basic Triangulation

Triangulation_3<TriangulationTraits_3,TriangulationDataStructure_3> expects a model of a geometric traits class as its first template argument and a model of a triangulation data structure as its second argument.

The requirements and defaults for the traits classes are described in the reference pages for TriangulationTraits_3, DelaunayTriangulationTraits_3, TriangulationHierarchyTraits_3, RegularTriangulationTraits_3, and CGAL::Regular_triangulation_euclidean_traits_3<R,Weight>.

The requirements and default for the triangulation data structure are described in Chapter . However, a default parameter is defined in Triangulation_3, so that the user who has no additional needs can use Triangulation_3<TriangulationTraits_3> without specifying the second argument.

Delaunay Triangulation

The class Delaunay_triangulation_3<DelaunayTriangulationTraits_3,TriangulationDataStructure_3> represents a three-dimensional Delaunay triangulation. This Delaunay triangulation is fully dynamic: it supports both insertions and vertex removal.

The user is advised to use the class Triangulation_hierarchy_3<Tr> in order to benefit from an increased efficiency.

Triangulation hierarchy

The class Triangulation_hierarchy_3<Tr> implements a triangulation augmented with a data structure that allows fast point location queries. Thus, it allows fast construction of the triangulation. As proved in [Dev98], this structure has an optimal behaviour when it is built for Delaunay triangulations. However it can be used as well for other triangulations and the class is templated by a parameter, which is to be instantiated by one of the CGAL triangulation classes. It offers the same functionalities as the Tr parameter class.

Note that, since the algorithms that are provided are randomized, the running time of constructing a triangulation with a hierarchy may be improved when shuffling the data points.

Regular Triangulation

Regular_triangulation_3<RegularTriangulationTraits_3,TriangulationDataStructure_3> implements incremental regular triangulations.

Let $$S^(w) be a set of weighted points in $$ ³. Let $$p^(w)=(p,w_p), p ³, w_p and $$z^(w)=(z,w_z), z ³, w_z be two weighted points. A weighted point $$p^(w)=(p,w_p) can also be seen as a sphere of center $$p and radius $$w_p. The power product between $$p^(w) and $$z^(w) is defined as

$$(p^(w),z^(w)) = p-z ²-w_p²-w_z²

where $$ p-z is the Euclidean distance between $$p and $$z. $$p^(w) and $$z^(w) are said to be orthogonal if $$(p^(w)-z^(w)) = 0 (see Figure ).

Figure:  Orthogonal weighted points (picture in 2D).

Four weighted points have a unique common orthogonal weighted point called the power sphere. The weighted point orthogonal to three weighted points in the plane defined by these three points is called the power circle. The power segment will denote the weighted point orthogonal to two weighted points on the line defined by these two points.

A sphere $$z^(w) is said to be regular if $$ p^(w) S^(w), (p^(w)-z^(w)) 0.

A triangulation of $$S^(w) is regular if the power spheres of all simplices are regular.

The regular triangulation of $$S^(w) is in fact the projection onto $$ ³ of the convex hull of the four-dimensional points $$(p, p-O ²-w_p), for $$p^(w)=(p,w_p) S^(w). Note that all points of $$S^(w) do not necessarily appear as vertices of the regular triangulation. To know more about regular triangulations, see for example [ES96].

When all weights are 0, power spheres are nothing more than circumscribing spheres, and the regular triangulation is exactly the Delaunay triangulation.

Examples

First example

// Triangulation_3/example_simple.C
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Filtered_kernel.h>

#include <CGAL/Triangulation_3.h>

#include <iostream>
#include <fstream>
#include <cassert>
#include <list>
#include <vector>

typedef CGAL::Filtered_kernel<CGAL::Simple_cartesian<double> > my_K;

// This is just to shorten some symbol names for VC++
struct K : public my_K {};

typedef CGAL::Triangulation_3<K> Triangulation;

typedef Triangulation::Cell_handle Cell_handle;
typedef Triangulation::Vertex_handle Vertex_handle;
typedef Triangulation::Locate_type Locate_type;

typedef K::Point_3 Point;

int main()
{
  Triangulation T;

  // insertion from a list :
  std::list<Point> L;
  L.push_front(Point(0,0,0));
  L.push_front(Point(1,0,0));
  L.push_front(Point(0,1,0));

  int n = T.insert(L.begin(), L.end());

  // insertion from a vector :
  std::vector<Point> V(3);
  V[0] = Point(0,0,1);
  V[1] = Point(1,1,1);
  V[2] = Point(2,2,2);

  n = n + T.insert(V.begin(), V.end());

  // 6 points have been inserted :
  assert( n == 6 );

  // checking validity of T :
  assert( T.is_valid(false) );

  Locate_type lt;
  int li, lj;
  Point p(0,0,0);
  Cell_handle c = T.locate(p, lt, li, lj);
  // p is the vertex of c of index li :
  assert( lt == Triangulation::VERTEX );
  assert(  c->vertex(li)->point() == p );

  Vertex_handle v = c->vertex( (li+1)&3 );
  // v is another vertex of c
  Cell_handle nc = c->neighbor(li);
  // nc = neighbor of c opposite to the vertex associated with p
  // nc must have vertex v :
  int nli;
  assert(  nc->has_vertex( v, nli ) );
  // nli is the index of v in nc

  std::ofstream oFileT("output",std::ios::out);
  // writing file output; 
  oFileT << T; 

  Triangulation T1;
  std::ifstream iFileT("output",std::ios::in);
  // reading file output; 
  iFileT >> T1; 
  assert( T1.is_valid() );
  assert( T1.number_of_vertices() == T.number_of_vertices() );
  assert( T1.number_of_cells() == T.number_of_cells() );

  return 0;
}

Example changing the vertex base

// Triangulation_3/example_color.C
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Filtered_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>

#include <cassert>

template < class Traits >
class My_vertex_base
  : public CGAL::Triangulation_vertex_base_3<Traits>
{
public :
  CGAL::Color color;
  typedef typename Traits::Point_3 Point;

  My_vertex_base() 
    : CGAL::Triangulation_vertex_base_3<Traits>(), color(CGAL::WHITE)
    {}
 
  My_vertex_base(const Point & p) 
    : CGAL::Triangulation_vertex_base_3<Traits>(p), color(CGAL::WHITE)
    {}

  My_vertex_base(const Point & p, void* c) 
    : CGAL::Triangulation_vertex_base_3<Traits>(p,c), color(CGAL::WHITE)
    {}

  My_vertex_base(void* c)
    : CGAL::Triangulation_vertex_base_3<Traits>(c), color(CGAL::WHITE)
    {}
};

typedef CGAL::Filtered_kernel<CGAL::Simple_cartesian<double> > my_K;

// This is just to shorten some symbol names for VC++
struct K : public my_K {};

typedef K::Point_3 Point;

typedef CGAL::Triangulation_cell_base_3<K> Cb;
typedef My_vertex_base<K> Vb;
typedef CGAL::Triangulation_data_structure_3<Vb,Cb> Tds;
typedef CGAL::Delaunay_triangulation_3<K, Tds> Delaunay;

typedef Delaunay::Vertex_iterator Vertex_iterator;
typedef Delaunay::Vertex_handle Vertex_handle;

int main()
{
  Delaunay T;

  T.insert(Point(0,0,0));
  T.insert(Point(1,0,0));  
  T.insert(Point(0,1,0));  
  T.insert(Point(0,0,1));  
  T.insert(Point(2,2,2));  
  T.insert(Point(-1,0,1));  

  Vertex_iterator vit;
  std::set<Vertex_handle> adjacent;
  for (vit = T.finite_vertices_begin(); vit != T.vertices_end(); ++vit) {
    T.incident_vertices( &*vit, adjacent);
    if (adjacent.size() == 6) 
      vit->color = CGAL::RED;
  }

  return 0;
}

Use of the Delaunay Hierarchy

// Triangulation_3/example_hierarchy.C
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Filtered_kernel.h>

#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_hierarchy_3.h>

#include <cassert>
#include <vector>

typedef CGAL::Filtered_kernel<CGAL::Simple_cartesian<double> > my_K;

// This is just to shorten some symbol names for VC++
struct K : public my_K {};

typedef CGAL::Triangulation_vertex_base_3<K>             Vb;
typedef CGAL::Triangulation_hierarchy_vertex_base_3<Vb>  Vbh;
typedef CGAL::Triangulation_cell_base_3<void>            Cb;
typedef CGAL::Triangulation_data_structure_3<Vbh,Cb>     Tds;
typedef CGAL::Delaunay_triangulation_3<K,Tds>            Dt;
typedef CGAL::Triangulation_hierarchy_3<Dt>              Dh;

typedef Dh::Vertex_iterator Vertex_iterator;
typedef Dh::Vertex_handle Vertex_handle;
typedef K::Point_3 Point;

int main()
{
  Dh T;

  // insertion of points on a 3D grid
  int x,y,z;
  std::vector<Vertex_handle> V(125);
  int i=0;
  
  for (z=0 ; z<5 ; z++)
    for (y=0 ; y<5 ; y++)
      for (x=0 ; x<5 ; x++) 
	  V[i++] = T.insert(Point(x,y,z));

  assert( T.is_valid() );
  assert( T.number_of_vertices() == 125 );
  assert( T.dimension() == 3 );

  // removal of the vertices in random order
  std::random_shuffle(V.begin(), V.end());

  for (i=0; i<125; ++i)
    assert( T.remove(V[i]) );

  assert( T.is_valid() );
  assert( T.number_of_vertices() == 0 );

  return 0;
}