\( \newcommand{\E}{\mathrm{E}} \) \( \newcommand{\A}{\mathrm{A}} \) \( \newcommand{\R}{\mathrm{R}} \) \( \newcommand{\N}{\mathrm{N}} \) \( \newcommand{\Q}{\mathrm{Q}} \) \( \newcommand{\Z}{\mathrm{Z}} \) \( \def\ccSum #1#2#3{ \sum_{#1}^{#2}{#3} } \def\ccProd #1#2#3{ \sum_{#1}^{#2}{#3} }\)
CGAL 4.4 - High-dimension Triangulation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Groups Pages
User Manual

Authors
Samuel Hornus and Olivier Devillers

This package proposes data structure and algorithms to compute triangulations of points in any dimensions. The Triangulation_data_structure allows to store and manipulate the combinatorial part of a triangulation while the geometric classes Triangulation and Delaunay_triangulation allows to compute a (Delaunay) triangulation of a set of points and to maintain it under insertions (and deletions in the Delaunay case).

Introduction

Some definitions

A finite abstract simplicial complex is built on a finite set of vertices \( V\) and consists of a collection \( S\) of subsets of \( V\) such that

if \( s\) is a set of vertices in \( S\), then all the subsets of \( s\) are also in \( S\).

The sets in \( S\) (which are subsets of \( V\)) are called faces or simplices (the singular of which is simplex). A simplex \( s\in S\) is maximal if it is not a proper subset of some other set in \( S\).

Modifications:

A simplex having \( d+1 \) vertices is said of dimension \( d \). The simplicial complex is pure if all the maximal simplices have the same dimension.

In the sequel, we will call these maximal simplices full cells. A face of a simplex is a subset of it. A proper face of a simplex is a strict subset of it.
Modifications:
A complex has no boundaries if any proper face of a simplex is also a proper face of another simplex. A pure complex is manifold if all faces of dimension \( d-1 \) are proper faces of exactly two simplices.

If the vertices are embedded into Euclidean space \( \mathbb{R}^d\), we deal with finite simplicial complexes which have slightly different simplices and additional requirements:

  • vertices corresponds to points in space.
  • a simplex \( s\in S\) is the convex hull of its vertices.
  • the vertices of a simplex \( s\in S\) are affinely independent.
  • the intersection of any two simplices of \( S\) is a proper face of both simplices (the empty set counts).

See the wikipedia entry for more about simplicial complexes.

What's in this package?

This CGAL package deals with pure manifold simplicial complexes which are connected and have no boundaries, which we will simply call in the sequel triangulations. It provides three main classes for creating and manipulating triangulations.

The class CGAL::Triangulation_data_structure<Dimensionality, TriangulationDSVertex, TriangulationDSFullCell> models an abstract triangulation: vertices in this class are not embedded in Euclidean space but are only of combinatorial nature.

Modifications:
The class CGAL::Triangulation<TriangulationTraits, TriangulationDataStructure> describes an embedded triangulation that has as vertices a given set of points and which fills the convex hull of these points.
Methods are provided for the insertion of points in the triangulation, the traversal of various elements of the triangulation, as well as the localization of a query point inside the triangulation. The convex hull of the points is part of the triangulation, the fact that there is no boundary is ensured by adding an infinite vertex and infinite full cells to triangulate the outside of the convex hull.
Modifications:
See Chapter 3D Triangulations for more details about infinite vertices and cells.

The class CGAL::Delaunay_triangulation<DelaunayTriangulationTraits, TriangulationDataStructure> adds further constraints to a triangulation, in that all its simplices must have the so-called Delaunay or empty-ball property: the interior of a ball circumscribing any simplex (or full cell) must be free from any vertex of the triangulation. The CGAL::Delaunay_triangulation class supports deletion of vertices.

Further definitions

An \( i\)-face denotes an \( i\)-dimensional simplex, or a simplex with \( i+1\) vertices. When these vertices are embedded in Euclidean space, they must be affinely independent.

If the maximal dimension of a simplex in the triangulation is \( d\), we call:

  • an \( i\)-face for some \( i\in[0,d]\) a face;
  • a \( 0\)-face a vertex;
  • a \( 1\)-face an edge;
  • a \( (d-2)\)-face a ridge;
  • a \( (d-1)\)-face a facet; and
  • a \( d\)-face a full cell.

Two faces \( \sigma\) and \( \sigma'\) are incident if and only if \( \sigma'\) is a proper sub-face of \( \sigma\) or vice versa.

Triangulation Data Structure

In this section, we describe the concept TriangulationDataStructure for which CGAL provides one model class: CGAL::Triangulation_data_structure<Dimensionality, TriangulationDSVertex, TriangulationDSFullCell>.

A TriangulationDataStructure can represent an abstract pure complex such that any facet is incident to exactly two full cells.

Modifications:
A TriangulationDataStructure has a maximal dimension which is a
positive integer equal to the maximum dimension a full cell can have. This maximal dimension can be chosen by the user at the creation of a TriangulationDataStructure and can then be queried using the method tds.maximal_dimension(). A TriangulationDataStructure also knows the current dimension of its full cells, which can be queried with tds.current_dimension(). In the sequel, let us denote the maximal dimension with \( D \) and the current dimension with \( d \). The inequalities \( -2 \leq d \leq D\) and \( 0 \le D\) always hold. The special meaning of negative values for \(d\) is explained below.

The data structure triangulates \( \mathcal S^d\)

A TriangulationDataStructure can be viewed as a triangulation of the topological sphere \( \mathcal S^d\), i.e., its faces can be embedded to form a partition of \( \mathcal S^d\) into \(d\)-simplices.

One nice consequence of the above important fact is that a full cell has always exactly \( d+1\) neighbors. Two full cells \( \sigma\) and \( \sigma'\) sharing a facet are called neighbors.

Modifications:
Possible values of \(d\) (the current dimension of the triangulation) include

\(d=-2\)
This corresponds to the non-existence of any object in the TriangulationDataStructure.
\(d=-1\)

This corresponds to a single vertex.

\(d=0\)

This corresponds to two vertices, each corresponding to a full cell; the two full cells being neighbors of each other. This is the unique triangulation of the \( 0\)-sphere.

\( 0< d \le D\)
This corresponds to a standard triangulation of the sphere \( \mathcal S^d\).

The class Triangulation_data_structure

We give here some details about the class Triangulation_data_structure<Dimensionality, TriangulationDSVertex, TriangulationDSFullCell> implementing the concept TriangulationDataStructure.

Storage

A TriangulationDataStructure explicitly stores its vertices and full cells.

Each vertex stores a reference (a handle) to one of its incident full cells.

Each full cell stores references to its \( d+1\) vertices and neighbors. Its vertices and neighbors are indexed from \( 0\) to \( d \). The indices of its neighbors have the following meaning: the \( i\)-th neighbor of \( \sigma\) is the unique neighbor of \( \sigma\) that does not contain the \( i\)-th vertex of \( \sigma\); in other words, it is the neighbor of \( \sigma\) opposite to the \( i\)-th vertex of \( \sigma\) (Figure Figure 40.1).

Modifications:
Faces of dimension between 0 and \( d-1 \) can be accessed as subfaces of a full cell.

simplex-structure.png
Figure 40.1 Indexing the vertices and neighbors of a full cell \( c\) in dimension \( d=2\).

Template parameters

The Triangulation_data_structure<Dimensionality, TriangulationDSVertex, TriangulationDSFullCell> class template is designed in such a way that its user can choose

  • the maximal dimension of the triangulation data structure by specifying the Dimensionality template parameter,
  • the type used to represent vertices by specifying the TriangulationDSVertex template parameter and
  • the type used to represent full cells by specifying the TriangulationDSFullCell template parameter.

The last two parameters have default values and are thus not necessary, unless the user needs custom types (see the reference manual page for this class template). The first template parameter, Dimensionality, must be one of the following:

  • CGAL::Dimension_tag<D> for some integer \( D \). This indicates that the triangulation can store full cells of dimension at most \( D \). The maximum dimension \( D \) is known by the compiler, which triggers some optimizations.
  • CGAL::Dynamic_dimension_tag. In this case, the maximum dimension of the full cells must be passed as an integer argument to an instance constructor (see TriangulationDataStructure).

The TriangulationDSVertex and TriangulationDSFullCell parameters to the class template must be models of the concepts TriangulationDSVertex and TriangulationDSFullCell respectively. CGAL provides models for these concepts: Triangulation_ds_vertex<TriangulationDataStructure> and Triangulation_ds_full_cell<TriangulationDataStructure, TDSFullCellStoragePolicy>, which, as one can see, take the TriangulationDataStructure as a template parameter in order to get access to some nested types in TriangulationDataStructure.

Modifications:
The default values are CGAL::Triangulation_ds_vertex<TDS> and CGAL::Triangulation_ds_full_cell<TDS> where TDS is the current class Triangulation_data_structure<Dimensionality, TriangulationDSVertex, TriangulationDSFullCell> This creates a circular dependency, which we resolve in the same way as in the CGAL Triangulation_2 and Triangulation_3 packages (see Chapters chapterTDS2, chapterTriangulation2, chapterTDS3, and chapterTriangulation3). In particular, models of the concepts TriangulationDSVertex and TriangulationDSFullCell must provide a nested template Rebind_TDS which is documented in those two concept's reference manual pages. This mechanism can be used to provide a custom vertex or full cell class. The user is encouraged to read the documentation of the CGAL Triangulation_2 or Triangulation_3 package.

Examples

Incremental Construction

The following examples shows how to construct a triangulation data structure by inserting vertices. Its main interest is that it demonstrates most of the API to insert new vertices into the triangulation.


File triangulation_data_structure_static.cpp

#include <CGAL/Triangulation_data_structure.h>
#include <iostream>
#include <vector>
int main()
{
TDS S;
assert( 7 == S.maximal_dimension() );
assert( -2 == S.current_dimension() );
assert( S.is_valid() );
std::vector<TDS::Vertex_handle> V(10);
V[0] = S.insert_increase_dimension(); //insert first vertex
assert( -1 == S.current_dimension() );
for( int i = 1; i <= 5; ++i )
V[i] = S.insert_increase_dimension(V[0]);
// the first 6 vertices have created a triangulation
// of the 4-dimensional topological sphere
// (the boundary of a five dimensional simplex).
assert( 4 == S.current_dimension() );
assert( 6 == S.number_of_vertices() );
assert( 6 == S.number_of_full_cells() );
TDS::Full_cell_handle c = V[5]->full_cell();
V[6] = S.insert_in_full_cell(c);
// full cell c is split in 5
assert( 7 == S.number_of_vertices() );
assert( 10 == S.number_of_full_cells() );
c = V[3]->full_cell();
TDS::Facet ft(c, 2); // the Facet opposite to vertex 2 in c
V[7] = S.insert_in_facet(ft);
// facet ft is split in 4 and the two incident cells are split accordingly
assert( 8 == S.number_of_vertices() );
assert( 16 == S.number_of_full_cells() );
c = V[3]->full_cell();
TDS::Face face(c);
// meant to contain the edge joining vertices 2 and 4 of full_cell c
face.set_index(0, 2); // namely vertex 2
face.set_index(1, 4); // and vertex 4
V[8] = S.insert_in_face(face);
// face is split in 2, and all incident full cells also
assert( S.is_valid() );
TDS::Full_cell_handle hole[2];
hole[0] = V[8]->full_cell();
hole[1] = hole[0]->neighbor(0);
// the hole is made of two adjacent full cells
ft = TDS::Facet(hole[0], 1); // a face on the boundary of hole[0]
V[9] = S.insert_in_hole(hole, hole+2, ft);
// the hole is triangulated by linking a new vertex to its boundary
assert( S.is_valid() );
return 0;
}

In the previous example, the maximal dimension is fixed at compile time. It is also possible to fix it at run time, as in the next example.


File triangulation_data_structure_dynamic.cpp

#include <CGAL/Triangulation_data_structure.h>
#include <iostream>
#include <vector>
int main()
{
const int ddim = 5; // dimension of TDS with dynamic dimension
TDS;
typedef TDS::Vertex_handle Vertex_handle;
TDS D(ddim); // the argument is taken into account.
assert( ddim == D.maximal_dimension() );
assert( -2 == D.current_dimension() );
assert( D.is_valid() );
std::vector<Vertex_handle> V(5);
V[0] = D.insert_increase_dimension();
V[1] = D.insert_increase_dimension(V[0]);
V[2] = D.insert_increase_dimension(V[0]);
V[3] = D.insert_increase_dimension(V[0]);
V[4] = D.insert_in_full_cell(V[3]->full_cell());
assert( 6 == D.number_of_full_cells() );
assert( 2 == D.current_dimension() );
assert( D.is_valid() );
return 0;
}

Barycentric subdivision

This example provides a function for computing the barycentric subdivision of a single full cell c in a triangulation data structure. The other full cells adjacent to c are automatically subdivided to match the subdivision of the full cell c. The barycentric subdivision of c is obtained by enumerating all the faces of c in order of decreasing dimension, from the dimension of c to dimension 1, and inserting a new vertex in each face. For the enumeration, we use a combination enumerator, which is not documented, but provided in CGAL.

barycentric-subdivision.png
Figure 40.2 Barycentric subdivision in dimension \( d=2\).


File barycentric_subdivision.cpp

#include <CGAL/Triangulation_data_structure.h>
#include <CGAL/internal/Combination_enumerator.h>
#include <iostream>
#include <vector>
template< typename TDS >
void find_face_from_vertices(const TDS & tds,
const std::vector<typename TDS::Vertex_handle> & face_vertices,
typename TDS::Face & face);
template< typename TDS >
void barycentric_subdivide(TDS & tds, typename TDS::Full_cell_handle fc)
{ /* This function builds the barycentric subdivision of a single
full cell |fc| from a triangulation data structure |tds|. */
typedef typename TDS::Full_cell_handle Full_cell_handle;
typedef typename TDS::Vertex_handle Vertex_handle;
typedef typename TDS::Face Face;
const int dim = tds.current_dimension();
// First, read handles to the cell's vertices
std::vector<Vertex_handle> vertices;
std::vector<Vertex_handle> face_vertices;
for( int i = 0; i <= dim; ++i ) vertices.push_back(fc->vertex(i));
// Then, subdivide the cell |fc| once by inserting one vertex
tds.insert_in_full_cell(fc);
// From now on, we can't use the variable |fc|...
// Then, subdivide faces of |fc| in order of decreasing dimension
for( int d = dim-1; d > 0; --d )
{
face_vertices.resize(d+1);
// The following class
// enumerates all (d+1)-tuple of the set {0, 1, ..., dim}
CGAL::internal::Combination_enumerator combi(d+1, 0, dim);
while( ! combi.end() )
{
for( int i = 0; i <= d; ++i )
face_vertices[i] = vertices[combi[i]];
// we need to find a face with face_vertices
Face face(dim);
find_face_from_vertices(tds, face_vertices, face);
tds.insert_in_face(face);
++combi;
}
}
}
template< typename TDS >
void find_face_from_vertices( const TDS & tds,
const std::vector<typename TDS::Vertex_handle> & face_vertices,
typename TDS::Face & face)
{ /* The main goal of this function is to find a full cell that
contains a given set of vertices |face_vertices|. Then, it
builds a corresponding |face|. */
typedef typename TDS::Face Face;
typedef typename TDS::Vertex_handle Vertex_handle;
typedef typename TDS::Full_cell_handle Full_cell_handle;
typedef typename TDS::Full_cell::Vertex_handle_iterator Vertex_h_iterator;
// get the dimension of the face we want to build
int fdim(face_vertices.size() - 1);
if( fdim <= 0) exit(-1);
// find all full cells incident to the first vertex of |face|
typedef std::vector<Full_cell_handle> Cells;
Cells cells;
std::back_insert_iterator<Cells> out(cells);
tds.incident_full_cells(face_vertices[0], out);
// Iterate over the cells to find one which contains the face_vertices
for( typename Cells::iterator cit = cells.begin(); cit != cells.end(); ++cit){
// find if the cell *cit contains the Face |face|
int i(0);
for( ; i <= fdim; ++i ) {
Vertex_handle face_v(face_vertices[i]);
bool found(false);
Vertex_h_iterator vit = (*cit)->vertices_begin();
for( ; vit != (*cit)->vertices_end(); ++vit ) {
if( *vit == face_v ) {
found = true;
break;
}
}
if( ! found )
break;
}
if( i > fdim ) {// the full cell |*cit| contains |face|
face.set_full_cell(*cit);
for( int i = 0; i <= fdim; ++i )
face.set_index(i, (*cit)->index(face_vertices[i]));
return;
}
}
std::cerr << "Could not build a face from vertices"<<std::endl;
assert(false);
}
int main()
{
const int sdim = 5; // dimension of TDS with compile-time dimension
TDS;
TDS tds(sdim);
TDS::Vertex_handle one_vertex = tds.insert_increase_dimension();
for( int i = 1; i < sdim+2; ++i )
tds.insert_increase_dimension(one_vertex);
// we get a triangulation of space of dim sdim homeomorphic to
// the boundary of simplex of dimension sdim+1 with sdim+2 vertices
assert( sdim == tds.current_dimension() );
assert( 2+sdim == tds.number_of_vertices() );
assert( 2+sdim == tds.number_of_full_cells() );
barycentric_subdivide(tds, tds.full_cells_begin());
// The number of full cells should be twice the factorial of
// |tds.current_dimension()+1|. Eg, 1440 for dimension 5.
std::cout << "Triangulation has "
<< tds.number_of_full_cells() << " full cells";
assert( tds.is_valid() );
std::cout << " and is valid!"<<std::endl;
return 0;
}

Triangulations

Modifications:
The class CGAL::Triangulation<TriangulationTraits, TriangulationDataStructure> maintains a geometric triangulation in Euclidean space. More precisely, it maintains a triangulation (a partition into pairwise interior-disjoint full cells) of the convex hull of the points (the embedded vertices) of the triangulation. A special vertex at infinity is added to the convex hull facets to create infinite full cells and make the triangulation homeomorphic to a sphere of one dimension higher.

Methods are provided for the insertion of points in the triangulation, the contraction of faces, the traversal of various elements of the triangulation as well as the localization of a query point inside the triangulation.

Infinite full cells outside the convex hull are each incident to a finite facet on the convex hull of the triangulation and to a unique vertex at infinity.

Modifications:
The ordering of the vertices of a full cell defines an orientation of that full cell. As long as no advanced class method is called, it is guaranteed that all finite full cells have positive orientation. The infinite full cells are oriented as if the infinite vertex have been on the the other side of the hyperplane supported by the convex hull facets that the other points.

Implementation

The class CGAL::Triangulation<TriangulationTraits, TriangulationDataStructure> stores a model of the concept TriangulationDataStructure which is instantiated with a vertex type that stores a point, and a full cell type that allows the retrieval of the point of its vertices.

The template parameter TriangulationTraits must be a model of the concept TriangulationTraits which provides the geometric Point type as well as various geometric predicates used by the Triangulation class.

Modifications:
The dimension of the predicates provided by TriangulationTraits must match the dimension of the TriangulationDataStructure.
Modifications:
The template parameter TriangulationDataStructure must be a model of the concept TriangulationDataStructure which provides the triangulation data structure as described in the previous section.

Examples

Incremental Construction

The following example shows how to construct a triangulation in which we insert random points. In STEP 1, we generate one hundred random points in \( \mathcal R^5\), which we then insert into a triangulation. In STEP 2, we ask the triangulation to construct the set of edges ( \( 1\) dimensional faces) incident to the vertex at infinity. It is easy to see that these edges are in bijection with the vertices on the convex hull of the points. This gives us a handy way to count the convex hull vertices

Modifications:
(include files triangulation1.cpp and triangulation2.cpp are given and commented below).


File triangulation.cpp

#include <CGAL/Epick_d.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/Triangulation.h>
#include <CGAL/algorithm.h>
typedef CGAL::Triangulation<K> Triangulation;
int main()
{
const int D = 5; // we work in Euclidean 5-space
const int N = 100; // we will insert 100 points
// - - - - - - - - - - - - - - - - - - - - - - - - STEP 1
CGAL::Random_points_in_cube_d<Triangulation::Point> rand_it(D, 1.0);
std::vector<Triangulation::Point> points;
CGAL::cpp11::copy_n(rand_it, N, std::back_inserter(points));
Triangulation t(D); // create triangulation
assert(t.empty());
t.insert(points.begin(), points.end()); // compute triangulation
assert( t.is_valid() );
// - - - - - - - - - - - - - - - - - - - - - - - - STEP 2
typedef Triangulation::Face Face;
typedef std::vector<Face> Faces;
Faces edges;
std::back_insert_iterator<Faces> out(edges);
t.tds().incident_faces(t.infinite_vertex(), 1, out);
// collect faces of dimension 1 (edges) incident to the infinite vertex
std::cout << "There are " << edges.size()
<< " vertices on the convex hull." << std::endl;
#include "triangulation1.cpp" // See below
#include "triangulation2.cpp"
return 0;
}

Traversing the facets of the convex hull

Remember that a triangulation triangulates the convex hull of its vertices. Each facet of the convex hull is incident to one finite full cell and one infinite full cell. In fact there is a bijection between the infinite full cells and the facets of the convex hull. If vertices are not in general position, convex hull faces that are not simplices are triangulated. So, in order to traverse the convex hull facets, there are (at least) two possibilities:

The first is to iterate over the full cells of the triangulation and check if they are infinite or not:


File triangulation1.cpp

{ int i=0;
typedef Triangulation::Full_cell_iterator Full_cell_iterator;
typedef Triangulation::Facet Facet;
for( Full_cell_iterator cit = t.full_cells_begin();
cit != t.full_cells_end(); ++cit )
{
if( ! t.is_infinite(cit) )
continue;
Facet ft(cit, cit->index(t.infinite_vertex()));
++i;// |ft| is a facet of the convex hull
}
std::cout << "There are " << i << " facets on the convex hull."<< std::endl;
}

A second possibility is to ask the triangulation to gather all the full cells incident to the infinite vertex: they form precisely the set of infinite full cells:


File triangulation2.cpp

{ int i=0;
typedef Triangulation::Full_cell_handle Full_cell_handle;
typedef Triangulation::Facet Facet;
typedef std::vector<Full_cell_handle> Full_cells;
Full_cells infinite_full_cells;
std::back_insert_iterator<Full_cells> out(infinite_full_cells);
t.incident_full_cells(t.infinite_vertex(), out);
for( Full_cells::iterator sit = infinite_full_cells.begin();
sit != infinite_full_cells.end(); ++sit )
{
Facet ft(*sit, (*sit)->index(t.infinite_vertex()));
++i; // |ft| is a facet of the convex hull
}
std::cout << "There are " << i << " facets on the convex hull."<< std::endl;
}
Modifications:
One important difference between the two examples above is that the first uses little memory but traverses all the full cells, while the second visits only the infinite full cells but stores handles to them into the potentially big array infinite_full_cells.

Delaunay Triangulations

The class CGAL::Delaunay_triangulation<DelaunayTriangulationTraits, TriangulationDataStructure> derives from CGAL::Triangulation<DelaunayTriangulationTraits, TriangulationDataStructure> and adds further constraints to a triangulation, in that all its full cells must have the so-called Delaunay or empty-ball property: the interior of the ball circumscribing any full cell must be free from any vertex of the triangulation.

The circumscribing ball of a full cell is the ball having all vertices of the full cell on its boundary. In case of degeneracies (co-spherical points) the triangulation is not uniquely defined; Note however that the CGAL implementation computes a unique triangulation even in these cases.

When a new point p is inserted into a Delaunay triangulation, the finite full cells whose circumscribing sphere contain p are said to be in conflict with point p. The set of full cells that are in conflict with p form the conflict zone. That conflict zone is augmented with the infinite full cells whose finite facet does not lie anymore on the convex hull of the triangulation (with p added). The full cells in the conflict zone are removed, leaving a hole that contains p. That hole is ``star shaped'' around p and thus is easily re-triangulated using p as a center vertex.

Delaunay triangulations also support vertex removal.

Implementation

The class CGAL::Delaunay_triangulation<DelaunayTriangulationTraits, TriangulationDataStructure> derives from CGAL::Triangulation<DelaunayTriangulationTraits, TriangulationDataStructure>. It thus stores a model of the concept TriangulationDataStructure which is instantiated with a vertex type that stores a geometric point and allows its retrieval.

The template parameter DelaunayTriangulationTraits must be a model of the concept DelaunayTriangulationTraits which provides the geometric Point type as well as various geometric predicates used by the Delaunay_triangulation class. The concept DelaunayTriangulationTraits refines the concept TriangulationTraits by requiring a few other geometric predicates, necessary for the computation of Delaunay triangulations.

Examples

Access to the conflict zone and creating full cells during point insertion

When using a full cell type containing additional custom information, it may be useful to get an efficient access to the full cells that are going to be erased upon the insertion of a new point in the Delaunay triangulation, and to the newly created full cells. The second part of code example below shows how one can have efficient access to both the conflict zone and the created full cells, while still retaining an efficient update of the Delaunay triangulation.


File delaunay.cpp

#include <CGAL/Epick_d.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/Delaunay_triangulation.h>
#include <CGAL/algorithm.h>
#include <CGAL/Timer.h>
const int D=5;
< void, CGAL::TDS_full_cell_default_storage_policy > TDS_cell;
CGAL::Dimension_tag<D>, Vertex, Cell > TDS;
int main(int argc, char **argv)
{
int N = 100; if( argc > 2 )N = atoi(argv[1]); // number of points
CGAL::Timer cost; // timer
// Instanciate a random point generator
CGAL::Random rng(0);
typedef CGAL::Random_points_in_cube_d<T::Point> Random_points_iterator;
Random_points_iterator rand_it(D, 1.0, rng);
// Generate N random points
std::vector<T::Point> points;
CGAL::cpp11::copy_n(rand_it, N, std::back_inserter(points));
T t(D);
assert(t.empty());
// insert the points in the triangulation
cost.reset();cost.start();
std::cout << " Delaunay triangulation of "<<N<<" points in dim "<<D<< std::flush;
t.insert(points.begin(), points.end());
std::cout << " done in "<<cost.time()<<" seconds." << std::endl;
assert( t.is_valid() );
// insert with special operations in conflict zone and new created cells
cost.reset();
std::cout << " adding "<<N<<" other points "<< std::endl;
for(int i=0; i<N; ++i)
{
T::Vertex_handle v;
T::Face f(t.current_dimension());
T::Facet ft;
T::Full_cell_handle c;
T::Locate_type lt;
typedef std::vector<T::Full_cell_handle> Full_cells;
Full_cells zone, new_full_cells;
std::back_insert_iterator<Full_cells> out(zone);
c = t.locate(*++rand_it, lt, f, ft, v);
// previously inserted vertex v is used as hint for point location (if defined)
T::Facet ftc = t.compute_conflict_zone(*rand_it, c, out);
std::cout<<i<<" conflict zone of size "<<zone.size()<<" -> "<<std::flush;
out = std::back_inserter(new_full_cells);
assert( t.is_valid() );
v = t.insert_in_hole(*rand_it, zone.begin(), zone.end(), ftc, out);
std::cout<<new_full_cells.size()<<" new cells"<<std::endl;
for (Full_cells::iterator it=new_full_cells.begin();
it!=new_full_cells.end(); ++it) (*it)->data() = zone.size();
}
std::cout << " done in "<<cost.time()<<" seconds." << std::endl;
return 0;
}

Complexity and Performances

The current implementation locates points by walking in the triangulation, and sorts the points with spatial sort to insert a set of points. Thus the theoretical complexity are \( O(n\log n)\) for inserting \( n\) random points and \( O(n^{\frac{1}{d}})\) for inserting one point in a triangulation of \( n\) random points. In the worst case, the expected complexity is \( O(n^{\lceil\frac{d}{2}\rceil}+n\log n)\).

The actual timing are the following:

This section will be completed, when the code will be fully ready (and preferably with the new Kernel).

Design and Implementation History

This package is heavily inspired by the works of Monique Teillaud and Sylvain Pion (Triangulation_3) and Mariette Yvinec (Triangulation_2). The first version was written by Samuel Hornus and then pursued by Samuel Hornus and Olivier Devillers.