CGAL 6.1 - 3D Periodic Triangulations

Loading...
Searching...
No Matches
User Manual

Authors
  Manuel Caroli, Aymeric Pellé, Mael Rouxel-Labbé and Monique Teillaud

The periodic 3D-triangulation class of CGAL is designed to represent the triangulations of a set of points in the three-dimensional flat torus. The triangulation forms a partition of the space it is computed in. It is a simplicial complex, i.e. it contains all incident \( j\)-simplices ( \( j<k\)) of any \( k\)-simplex and two \( k\)-simplices either do not intersect or share a common \( j\)-face, \( j<k\). The occurring simplices of dimension up to three are called vertex, edge, facet, and cell, respectively.

# The Flat Torus

The 3D Periodic Triangulation package computes triangulations in the space \( \mathbb T_c^3\), which is defined as follows: Let \( c\in\mathbb R\setminus\{0\}\) and \( G\) be the group \( (c\cdot\mathbb Z^3, +)\), where \( c\cdot\mathbb Z\) denotes the set containing all integer multiples of \( c\). The flat torus is the quotient space: \( \mathbb T_c^3:=\mathbb R^3/G\). The parameter \( c\) defines the period.

The elements of \( \mathbb T_c^3\) are the equivalence classes of sets of points in \( \mathbb R^3\). We call these points representatives of an element of \( \mathbb T_c^3\). The implementation works not directly on elements of \( \mathbb T_c^3\) but on some representatives in \( \mathbb R^3\). So there need to be distinguished representatives to work on. Given \( \alpha\), \( \beta\), and \( \gamma\), the cube \( [\alpha,\alpha+c)\times[\beta,\beta+c)\times[\gamma,\gamma+c)\) contains exactly one representative of each element in \( \mathbb T_c^3\). We call it original domain. From now on, when we talk about points, we generally mean representatives of elements of \( \mathbb T_c^3\) that lie inside the original domain. Note that any input point is required to be an element of the half-open cube representing the original domain as defined above.

There are simplices containing points inside the original domain but also points outside it. The points outside the original domain are periodic copies of points inside the original domain. So, to specify a simplex we need points together with some additional information that determines the respective periodic copy of each point. The set of representatives of an element of \( \mathbb T_c^3\) is a cubic point grid. We address each representative by a three-dimensional integer vector \( (o_x,o_y,o_z)\), called offset. It represents the number of periods a representative in the original domain must be translated in \( x\)-, \( y\)-, and \( z\)-direction. The vector \( (0,0,0)\) corresponds to the representative in the original domain. To specify a \( k\)-simplex we need \( k+1\) point-offset pairs (cf. Figure 51.1).

Figure 51.1 Offsets in a cell.

# Representation

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, their corresponding offsets, 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 orientation of a simplex in \( \mathbb T_c^3\) is defined as the orientation of the corresponding simplex in \( \mathbb R^3\) given by representatives determined by the respective offsets (see Figure 51.2).

Figure 51.2 Orientation of a cell.

As in the underlying combinatorial triangulation (see Chapter 3D Triangulation Data Structure), the neighbors of a cell are 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. Also 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 49.1 of chapter 3D Triangulation Data Structure.

Some point sets do not admit a triangulation in \( \mathbb T_c^3\). In this case we use 27 periodic copies of the point set arranged in a cube of edge length \( 3c\). Any point set constructed in this way has a triangulation in \( \mathbb R^3/G'\) with \( G'=((3c\cdot\mathbb Z)^3,+)\) [1]. So we compute the triangulation in this space, which is a 27-sheeted covering space of \( \mathbb T_c^3\) (see Figure 51.3).

Figure 51.3 The same periodic triangulation in the 1-sheeted covering space and the 27-sheeted covering space.

The machinery that manages the copies is largely hidden from the user. However there are some effects that cannot be ignored. For example if the point set does not permit a triangulation in \( \mathbb T_c^3\) then the combinatorial iterators (`Cell_iterator`, `Facet_iterator`, `Edge_iterator`, and `Vertex_iterator`) return all simplices that are internally stored, which correspond to 27 periodic copies of each geometric primitive (Tetrahedron, Triangle, Segment, and Point). This is necessary to ensure consistency in the adjacency relations. In case it is desired to have only one periodic copy of each primitive, we provide geometric iterators. They return geometric primitives of the triangulation without relations between them. Another effect is that when the algorithm switches from the 27-sheeted covering space to the 1-sheeted covering space, the `Vertex_handle`s and `Cell_handle`s referencing deleted items become invalid.

In the data structure each vertex stores the input point it corresponds to. If we are computing in the 27-sheeted covering space, each vertex stores the representative inside the original domain it corresponds to. So, the 27 vertices corresponding to the same element of \( \mathbb T_c^3\) all store the same representative in \( \mathbb R^3\), and not different periodic copies.

Validity

A periodic triangulation is said to be `locally valid` iff

(a)-(b) Its underlying combinatorial graph, the triangulation data structure, is `locally valid` (see Section Representation of Chapter 3D Triangulation Data Structure)

(c) Any cell has its vertices ordered according to positive orientation. See Figure 51.2.

# Delaunay Triangulation

The class `Periodic_3_Delaunay_triangulation_3` implements Delaunay triangulations of point sets in \( \mathbb T_c^3\).

Delaunay triangulations have the empty sphere property, that is, the circumscribing sphere of each cell does not contain any other vertex of the triangulation in its interior. These triangulations are uniquely defined except in degenerate cases where five points are co-spherical. Note however that the CGAL implementation computes a uniquely defined triangulation even in these cases [2].

This implementation is fully dynamic: it supports both point insertion and vertex removal.

# Regular Triangulation

The class `Periodic_3_regular_triangulation_3` implements regular triangulations, also known as weighted Delaunay triangulations, of point sets in \( \mathbb T_c^3\).

A regular triangulation is a triangulation in which the power sphere of each simplex is regular. See Section Regular Triangulation for a complete definition. As for Delaunay triangulations, CGAL computes a uniquely defined regular triangulation even in degenerate cases [3].

The implementation is fully dynamic: it supports both point insertion and vertex removal.

# Triangulation Hierarchy

The class `Periodic_3_triangulation_hierarchy_3` is the adaptation of the hierarchical structure described in chapter 3D Triangulations to the periodic case.

# Software Design

We have chosen the prefix "Periodic_3" to emphasize that the triangulation is periodic in all three directions of space.

The main classes `Periodic_3_triangulation_3`, `Periodic_3_Delaunay_triangulation_3`, and `Periodic_3_regular_triangulation_3` provide high-level geometric functionality and are responsible for the geometric validity. They are built as layers on top of a triangulation data structure, which stores their combinatorial structure. This separation between the geometry and the combinatorics is reflected in the software design by the fact that the triangulation classes take two template parameters:

- the geometric traits class, which provides the type of points to use as well as the elementary operations on them (predicates and constructions). Furthermore it contains the offset type. The concept and models for this parameter are described in more detail in Section The Geometric Traits Parameter.
- the triangulation data structure class, which stores the combinatorial structure, described in Section The Triangulation Data Structure Parameter.

The class `Periodic_3_triangulation_3` contains all the functionality that is common to triangulations in general, such as location of a point in the triangulation [4], access functions, geometric queries like the orientation test etc. The class `Periodic_3_Delaunay_triangulation_3` contains all the functionality that is specific to Delaunay triangulations, such as point insertion and vertex removal, the side-of-sphere test, finding the conflicting region of a given point, computation of dual functions, etc. The class `Periodic_3_regular_triangulation_3` does the same for regular triangulations.

## The Geometric Traits Parameter

### Traits for Periodic Triangulations

The first template parameter of the triangulation class `Periodic_3_triangulation_3<Periodic_3TriangulationTraits_3, TriangulationDataStructure_3>` is the geometric traits class, described by the concept `Periodic_3TriangulationTraits_3`. It is different to the `TriangulationTraits_3` (see chapter 3D Triangulations) in that it implements all objects, predicates and constructions using offsets.

The class `CGAL::Periodic_3_triangulation_traits_3<TriangulationTraits_3,Periodic_3Offset_3>` provides the required functionality. It expects two template parameters: a model of the concept `TriangulationTraits_3` and a model of the concept `Periodic_3Offset_3`.

The second parameter `Periodic_3Offset_3` defaults to `CGAL::Periodic_3_offset_3`.

### Traits for Periodic Delaunay Triangulations

The first template parameter of the Delaunay triangulation class `Periodic_3_Delaunay_triangulation_3<Periodic_3DelaunayTriangulationTraits_3, TriangulationDataStructure_3>` is the geometric traits class, described by the concept `Periodic_3DelaunayTriangulationTraits_3`. It is different to the `DelaunayTriangulationTraits_3` (see chapter 3D Triangulations) in that it implements all objects, predicates and constructions using offsets.

The class `Periodic_3_Delaunay_triangulation_traits_3<DelaunayTriangulationTraits_3,Periodic_3Offset_3>` provides the required functionality. It expects two template parameters: a model of the concept `DelaunayTriangulationTraits_3` and a model of the concept `Periodic_3Offset_3`.

The second parameter `Periodic_3Offset_3` defaults to `CGAL::Periodic_3_offset_3`.

### Traits for Periodic Regular Triangulations

The first template parameter of the regular triangulation class `Periodic_3_regular_triangulation_3<Periodic_3RegularTriangulationTraits_3, TriangulationDataStructure_3>` is the geometric traits class, described by the concept `Periodic_3RegularTriangulationTraits_3`. It is different to the RegularTriangulationTraits_3 (see chapter 3D Triangulations) in that it implements all objects, predicates and constructions using offsets.

The class `Periodic_3_regular_triangulation_traits_3<RegularTriangulationTraits_3,Periodic_3Offset_3>` provides the required functionality. It expects two template parameters: a model of the concept `RegularTriangulationTraits_3` and a model of the concept `Periodic_3Offset_3`.

The second parameter `Periodic_3Offset_3` defaults to `Periodic_3_offset_3`.

### Compatible Kernels

The kernels `Cartesian`, `Homogeneous`, `Simple_cartesian`, `Simple_homogeneous` and `Filtered_kernel` can all be used as models of the concepts `TriangulationTraits_3`, `DelaunayTriangulationTraits_3`, and `RegularTriangulationTraits_3`. The periodic triangulation classes provide exact predicates and exact constructions if these respective template parameters do. They provide exact predicates but not exact constructions if `Filtered_kernel<CK>` with `CK` an inexact kernel is used as template parameter. Using `Exact_predicates_inexact_constructions_kernel` provides fast and exact predicates and not exact constructions, using `Exact_predicates_exact_constructions_kernel` provides fast and exact predicates and exact constructions. The latter is recommended if the dual constructions and constructions of points, segments, triangles, and tetrahedra are used.

## The Triangulation Data Structure Parameter

The second template parameter of the periodic triangulation classes is a triangulation data structure class. This parameter must meet the requirements described in the concept `TriangulationDataStructure_3`, for which CGAL offers the model `CGAL::Triangulation_data_structure_3`. To represent periodic triangulations, the face and vertex of the triangulation data structure must be also be models of extended concepts. The model `CGAL::Triangulation_data_structure_3` is parameterized by a vertex base class and a face base class, which gives the possibility to customize the vertices and faces used by the triangulation data structure.

A default value for the triangulation data structure parameter is provided in all the triangulation classes, so it does not need to be specified by the user unless he wants to use a different triangulation data structure or a different vertex or cell base class.

## Flexibility of the Design

The class `Periodic_3_triangulation_3` uses the `TriangulationDataStructure_3` in essentially the same way as the class `Triangulation_3` and the flexibility described in Software Design is therefore applicable in exactly the same way. Furthermore, the classes `Triangulation_vertex_base_with_info_3` and `Triangulation_cell_base_with_info_3` can be reused directly, see Example Adding a Color.

# Examples

## Basic Example

This example shows the incremental construction of a 3D Delaunay triangulation, the location of a point and how to perform elementary operations on indices in a cell. It uses the default parameter of the class `Periodic_3_Delaunay_triangulation_3` for the triangulation data structure.

File Periodic_3_triangulation_3/simple_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/periodic_3_triangulation_3_io.h>
#include <iostream>
#include <fstream>
#include <cassert>
#include <list>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> Gt;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<Gt> P3DT3;
typedef P3DT3::Point Point;
typedef P3DT3::Iso_cuboid Iso_cuboid;
typedef P3DT3::Vertex_handle Vertex_handle;
typedef P3DT3::Cell_handle Cell_handle;
typedef P3DT3::Locate_type Locate_type;
int main(int, char**)
{
Iso_cuboid domain(-1,-1,-1,2,2,2); // the fundamental domain
// construction from a list of points :
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));
P3DT3 T(L.begin(), L.end(), domain); // put the domain with the constructor
P3DT3::size_type n = T.number_of_vertices();
// insertion from a vector :
std::vector<Point> V(3);
V[0] = Point(0,0,1);
V[1] = Point(1,1,1);
V[2] = Point(-1,-1,-1);
n = n + T.insert(V.begin(), V.end());
assert( n == 6 ); // 6 points have been inserted
assert( T.is_valid() ); // checking validity of T
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 == P3DT3::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
// writing file output
std::ofstream oFileT("output.tri", std::ios::out); // as a .tri file
oFileT << T;
std::ofstream to_off("output_regular.off"); // as a .off file
CGAL::write_triangulation_to_off(to_off, T);
std::ofstream d_to_off("output_dual.off");
draw_dual_to_off(d_to_off, T);
// reading file output
P3DT3 T1;
std::ifstream iFileT("output.tri",std::ios::in);
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;
}
CGAL::Exact_predicates_inexact_constructions_kernel
CGAL::Periodic_3_Delaunay_triangulation_3
The class Periodic_3_Delaunay_triangulation_3 represents a Delaunay triangulation in three-dimensiona...
Definition: Periodic_3_Delaunay_triangulation_3.h:26
CGAL::Periodic_3_Delaunay_triangulation_traits_3
The class Periodic_3_Delaunay_triangulation_traits_3 is designed as a default traits class for the cl...
Definition: Periodic_3_Delaunay_triangulation_traits_3.h:23

## Changing the Vertex Base

The following two examples show how the user can plug his own vertex base in a triangulation. Changing the cell base is similar.

### Adding a Color

If the user does not need to add a type in a vertex that depends on the `TriangulationDataStructure_3` (e.g. a `Vertex_handle` or `Cell_handle`), then he can use the `Triangulation_vertex_base_with_info_3` class to add his own information easily in the vertices. The example below shows how to add a `Color` this way.

File Periodic_3_triangulation_3/colored_vertices.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/IO/Color.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> Gt;
typedef CGAL::Periodic_3_triangulation_ds_vertex_base_3<> VbDS;
typedef CGAL::Triangulation_vertex_base_3<Gt, VbDS> Vb;
typedef CGAL::Periodic_3_triangulation_ds_cell_base_3<> CbDS;
typedef CGAL::Triangulation_cell_base_3<Gt, CbDS> Cb;
typedef CGAL::Triangulation_vertex_base_with_info_3<CGAL::IO::Color, Gt, Vb> VbInfo;
typedef CGAL::Triangulation_data_structure_3<VbInfo, Cb> TDS;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<Gt, TDS> P3DT3;
typedef P3DT3::Point Point;
int main(int, char**)
{
P3DT3 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(.9,0,.1));
// Set the color of finite vertices of degree 6 to red.
P3DT3::Vertex_iterator vit;
for (vit = T.vertices_begin(); vit != T.vertices_end(); ++vit) {
if (T.degree(vit) == 16) {
vit->info() = CGAL::IO::red();
}
}
return 0;
}
CGAL::Periodic_3_Delaunay_triangulation_3::insert
Vertex_handle insert(const Point &p, Cell_handle start=Cell_handle())
Inserts point p in the triangulation and returns the corresponding vertex.
CGAL::Periodic_3_triangulation_ds_cell_base_3
The class Periodic_3_triangulation_ds_cell_base_3 is a model of the concept Periodic_3TriangulationDS...
Definition: Periodic_3_triangulation_ds_cell_base_3.h:20
CGAL::Periodic_3_triangulation_ds_vertex_base_3
The class Periodic_3_triangulation_ds_vertex_base_3 is a model of the concept Periodic_3Triangulation...
Definition: Periodic_3_triangulation_ds_vertex_base_3.h:20
CGAL::Triangulation_cell_base_3
CGAL::Triangulation_data_structure_3
CGAL::Triangulation_vertex_base_3
CGAL::Triangulation_vertex_base_with_info_3

### Adding Handles

If the user needs to add a type in a vertex that depends on the `TriangulationDataStructure_3` (e.g. a `Vertex_handle` or `Cell_handle`), then he has to derive his own vertex base class, as the following example shows.

File Periodic_3_triangulation_3/periodic_adding_handles.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Periodic_3_triangulation_ds_vertex_base_3.h>
#include <CGAL/Triangulation_vertex_base_3.h>
template < class Gt, class VbDS,
class Vb = CGAL::Triangulation_vertex_base_3<Gt,VbDS> >
class My_vertex_base
: public Vb
{
public:
typedef typename Vb::Vertex_handle Vertex_handle;
typedef typename Vb::Cell_handle Cell_handle;
typedef typename Vb::Point Point;
template < class TDS2 >
struct Rebind_TDS {
typedef typename Vb::template Rebind_TDS<TDS2>::Other Vb2;
typedef My_vertex_base<Gt, Vb2> Other;
};
My_vertex_base() {}
My_vertex_base(const Point& p)
: Vb(p) {}
My_vertex_base(const Point& p, Cell_handle c)
: Vb(p, c) {}
Vertex_handle vh;
Cell_handle ch;
};
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> Gt;
typedef CGAL::Periodic_3_triangulation_ds_vertex_base_3<> VbDS;
typedef CGAL::Periodic_3_triangulation_ds_cell_base_3<> CbDS;
typedef CGAL::Triangulation_cell_base_3<Gt,CbDS> Cb;
typedef CGAL::Triangulation_data_structure_3<My_vertex_base<Gt, VbDS>, Cb> TDS;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<Gt, TDS> P3DT3;
typedef P3DT3::Vertex_handle Vertex_handle;
typedef P3DT3::Point Point;
int main(int, char**)
{
P3DT3 T;
Vertex_handle v0 = T.insert(Point(0,0,0));
Vertex_handle v1 = T.insert(Point(.1,0,0));
Vertex_handle v2 = T.insert(Point(0,.1,0));
Vertex_handle v3 = T.insert(Point(0,0,.1));
Vertex_handle v4 = T.insert(Point(.2,.2,.2));
Vertex_handle v5 = T.insert(Point(.9,0,.1));
// Now we can link the vertices as we like.
v0->vh = v1;
v1->vh = v2;
v2->vh = v3;
v3->vh = v4;
v4->vh = v5;
v5->vh = v0;
return 0;
}

## The 27-sheeted Covering Space

The user can check at any time whether a triangulation would be a simplicial complex in \( \mathbb T_c^3\) and force a conversion if so. However, this should be done very carefully in order to ensure that the internal structure always remains a simplicial complex and thus a triangulation.

In this example, we construct a triangulation that can be converted to the 1-sheeted covering space. However, we can insert new points such that the point set does not have a Delaunay triangulation in the 1-sheeted covering space anymore, rendering the triangulation not extensible.

File Periodic_3_triangulation_3/covering.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <iostream>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> Gt;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<Gt> P3DT3;
typedef P3DT3::Point Point;
typedef P3DT3::Covering_sheets Covering_sheets;
int main(int, char**)
{
P3DT3 T;
// Input point grid (27 points)
for (double x=0. ; x < .9 ; x += 0.33) {
for (double y=0. ; y < .9 ; y += 0.33) {
for (double z=0. ; z < .9 ; z += 0.33) {
T.insert(Point(x,y,z));
}
}
}
Covering_sheets cs = T.number_of_sheets();
std::cout<<"Current covering: "<<cs[0]<<' '<<cs[1]<<' '<<cs[2]<<std::endl;
if ( T.is_triangulation_in_1_sheet() ) { // = true
bool is_extensible = T.is_extensible_triangulation_in_1_sheet_h1()
|| T.is_extensible_triangulation_in_1_sheet_h2(); // = false
T.convert_to_1_sheeted_covering();
cs = T.number_of_sheets();
std::cout<<"Current covering: "<<cs[0]<<' '<<cs[1]<<' '<<cs[2]<<std::endl;
if ( is_extensible ) // = false
std::cout<<"It is safe to change the triangulation here."<<std::endl;
else
std::cout<<"It is NOT safe to change the triangulation here!"<<std::endl;
T.convert_to_27_sheeted_covering();
cs = T.number_of_sheets();
std::cout<<"Current covering: "<<cs[0]<<' '<<cs[1]<<' '<<cs[2]<<std::endl;
}
std::cout<<"It is (again) safe to modify the triangulation."<<std::endl;
return 0;
}

## Large Point Set

Two optimizations are available for large point sets. Firstly, spatial sorting can be used to sort the input points according to a Hilbert curve, see chapter Spatial Sorting. Secondly, 36 appropriately chosen dummy points can be inserted to avoid the use of a 27-sheeted covering space in the beginning. These 36 dummy points are deleted in the end. If the point set turns out to not have a Delaunay triangulation in the 1-sheeted covering space, the triangulation is converted to the 27-sheeted covering space during the removal of the 36 dummy points. Note that this might take even longer than computing the triangulation without using this optimization. In general, uniformly distributed random point sets of more than 1000 points have a Delaunay triangulation in the 1-sheeted covering space.

The following example illustrates the use of these optimization techniques.

File Periodic_3_triangulation_3/large_point_set.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Random.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Timer.h>
#include <iostream>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> Gt;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<Gt> P3DT3;
typedef P3DT3::Point Point;
int main(int, char**)
{
CGAL::Timer t;
typedef CGAL::Creator_uniform_3<double, Point> Creator;
CGAL::Random random(7);
CGAL::Random_points_in_cube_3<Point, Creator> in_cube(.5, random);
int n = 10000;
std::vector<Point> pts;
P3DT3 PT1, PT2, PT3;
// Generating n random points
for (int i=0 ; i < n ; i++) {
Point p = *in_cube;
in_cube++;
pts.push_back(Point(p.x()+.5,p.y()+.5,p.z()+.5));
}
// Standard insertion
t.start();
for (int i=0 ; i < n ; i++) {
PT1.insert(pts[i]);
}
t.stop();
std::cout<<" Time: "<<t.time()<<" sec. (Standard insertion)"<<std::endl;
t.reset();
// Iterator range insertion using spatial sorting but no dummy points
t.start();
PT2.insert(pts.begin(), pts.end()); // third parameter defaults to false
t.stop();
std::cout<<" Time: "<<t.time()<<" sec. (with spatial sorting)"<<std::endl;
t.reset();
// Iterator range insertion using spatial sorting and dummy point heuristic
t.start();
PT3.insert(pts.begin(), pts.end(), true);
t.stop();
std::cout<<" Time: "<<t.time()<<" sec. (Dummy point heuristic)"<<std::endl;
return 0;
}
CGAL::Creator_uniform_3

## Geometric Access

Some application might use the geometric primitives of a triangulation as an input but not require a simplicial complex. We provide for these cases the geometric iterators that return only the geometric primitives fulfilling some properties. In the following example, we use the `Periodic_3_triangulation_3::Periodic_triangle_iterator` with the option `UNIQUE_COVER_DOMAIN`. This means that only the triangles that have a non-empty intersection with the original domain of the 1-sheeted covering space are returned, see Figure P3Triangulation3figgeom_iterators. The `Periodic_3_triangulation_3::Periodic_triangle` is actually a three-dimensional array of point-offset pairs. We check for all three entries of the periodic triangle whether the offset is (0,0,0) using the method `is_null`. If so, we convert the periodic triangle to a `PK::Triangle_3`, which requires exact constructions to be exact.

File Periodic_3_triangulation_3/geometric_access.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> Gt;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<Gt> P3DT3;
typedef Gt::Point_3 Point;
typedef Gt::Triangle_3 Triangle;
typedef P3DT3::Periodic_triangle Periodic_triangle;
typedef P3DT3::Periodic_triangle_iterator Periodic_triangle_iterator;
typedef P3DT3::Iterator_type Iterator_type;
int main(int, char**) {
P3DT3 T;
T.insert(Point(0,0,0));
T.insert(Point(0,0,0.5));
T.insert(Point(0,0.5,0.5));
T.insert(Point(0.5,0,0.5));
Periodic_triangle pt;
Triangle t_bd;
// Extracting the triangles that have a non-empty intersection with
// the original domain of the 1-sheeted covering space
for (Periodic_triangle_iterator ptit = T.periodic_triangles_begin(P3DT3::UNIQUE_COVER_DOMAIN);
ptit != T.periodic_triangles_end(P3DT3::UNIQUE_COVER_DOMAIN); ++ptit) {
pt = *ptit;
if (! (pt[0].second.is_null() && pt[1].second.is_null() && pt[2].second.is_null()) ) {
// Convert the current Periodic_triangle to a Triangle if it is
// not strictly contained inside the original domain.
// Note that this requires EXACT constructions to be exact!
t_bd = T.construct_triangle(pt);
}
}
}

## Periodic Regular Triangulations

The following five examples illustrate various features of 3D periodic regular triangulations.

### Basic example

This example shows the incremental construction of a 3D Delaunay triangulation, the location of a point and how to perform elementary operations on indices in a cell. It uses the default parameter of the class `Periodic_3_regular_triangulation_3` for the triangulation data structure.

File Periodic_3_triangulation_3/simple_regular_example.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_regular_triangulation_traits_3.h>
#include <CGAL/Periodic_3_regular_triangulation_3.h>
#include <CGAL/periodic_3_triangulation_3_io.h>
#include <iostream>
#include <fstream>
#include <cassert>
#include <list>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_regular_triangulation_traits_3<K> Gt;
typedef CGAL::Periodic_3_regular_triangulation_3<Gt> P3RT3;
typedef P3RT3::Bare_point Point;
typedef P3RT3::Weighted_point Weighted_point;
typedef P3RT3::Iso_cuboid Iso_cuboid;
typedef P3RT3::Vertex_handle Vertex_handle;
typedef P3RT3::Cell_handle Cell_handle;
typedef P3RT3::Locate_type Locate_type;
int main(int, char**)
{
Iso_cuboid domain(-1,-1,-1, 2,2,2); // the fundamental domain
// construction from a list of weighted points :
std::list<Weighted_point> L;
L.push_front(Weighted_point(Point(0,0,0), 0.01));
L.push_front(Weighted_point(Point(1,0,0), 0.02));
L.push_front(Weighted_point(Point(0,1,0), 0.03));
P3RT3 T(L.begin(), L.end(), domain); // put the domain with the constructor
P3RT3::size_type n = T.number_of_vertices();
// insertion from a vector :
std::vector<Weighted_point> V(3);
V[0] = Weighted_point(Point(0,0,1), 0.04);
V[1] = Weighted_point(Point(1,1,1), 0.05);
V[2] = Weighted_point(Point(-1,-1,-1), 0.06);
n = n + T.insert(V.begin(), V.end());
assert( n == 6 ); // 6 points have been inserted
assert( T.is_valid() ); // checking validity of T
Locate_type lt;
int li, lj;
Weighted_point p(Point(0,0,0), 1.);
Cell_handle c = T.locate(p, lt, li, lj);
// p is the vertex of c of index li :
assert( lt == P3RT3::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
// writing file output
std::ofstream oFileT("output_regular.tri", std::ios::out); // as a .tri file
oFileT << T;
std::ofstream to_off("output_regular.off");
CGAL::write_triangulation_to_off(to_off, T);
std::ofstream d_to_off("output_regular_dual.off"); // as a .off file
draw_dual_to_off(d_to_off, T);
// reading file output
P3RT3 T1;
std::ifstream iFileT("output_regular.tri",std::ios::in);
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;
}
CGAL::Periodic_3_regular_triangulation_3
The class Periodic_3_regular_triangulation_3 represents a weighted Delaunay triangulation in three-di...
Definition: Periodic_3_regular_triangulation_3.h:28
CGAL::Periodic_3_regular_triangulation_traits_3
The class Periodic_3_regular_triangulation_traits_3 is designed as a default traits class for the cla...
Definition: Periodic_3_regular_triangulation_traits_3.h:22

### Point Insertion and Vertex Removal

This example shows the incremental construction of a 3D regular triangulation, and the removal of a vertex. It uses the default parameter of the `Periodic_3_regular_triangulation_3` class for the triangulation data structure.

File Periodic_3_triangulation_3/p3rt3_insert_remove.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_regular_triangulation_traits_3.h>
#include <CGAL/Periodic_3_regular_triangulation_3.h>
#include <iostream>
typedef CGAL::Epick K;
typedef K::FT FT;
typedef CGAL::Periodic_3_regular_triangulation_traits_3<K> Gt;
typedef CGAL::Periodic_3_regular_triangulation_3<Gt> P3RT3;
typedef Gt::Weighted_point_3 Weighted_point_3;
typedef Gt::Point_3 Point_3;
typedef P3RT3::Vertex_handle Vertex_handle;
int main (int, char**)
{
P3RT3 p3rt3(P3RT3::Iso_cuboid(0,0,0, 1,1,1));
p3rt3.insert(Weighted_point_3(Point_3(0.1,0.1,0.1), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.9,0.1,0.1), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.1,0.9,0.1), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.1,0.1,0.9), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.4,0.4,0.4), 0.001));
while (p3rt3.number_of_vertices())
p3rt3.remove(p3rt3.vertices_begin());
std::cout << "EXIT SUCCESS" << std::endl;
return 0;
}

### Data Structure

This example shows the incremental construction of a 3D regular triangulation. It uses a more appropriate triangulation data structure, which saves some memory resources.

File Periodic_3_triangulation_3/p3rt3_insert_only.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_regular_triangulation_traits_3.h>
#include <CGAL/Periodic_3_regular_triangulation_3.h>
#include <CGAL/periodic_3_triangulation_3_io.h>
#include <CGAL/Regular_triangulation_vertex_base_3.h>
#include <CGAL/Regular_triangulation_cell_base_3.h>
#include <iostream>
typedef CGAL::Epick K;
typedef CGAL::Periodic_3_regular_triangulation_traits_3<K> Gt;
/* If remove() isn't called in our program, we can use a triangulation data structure
* more appropriate, which saves some memory resources.
*/
typedef CGAL::Regular_triangulation_vertex_base_3<Gt,
CGAL::Periodic_3_triangulation_ds_vertex_base_3<> > Vb;
typedef CGAL::Regular_triangulation_cell_base_3<Gt,
CGAL::Periodic_3_triangulation_ds_cell_base_3<> > Cb;
typedef CGAL::Triangulation_data_structure_3<Vb, Cb> Tds;
typedef CGAL::Periodic_3_regular_triangulation_3<Gt, Tds> P3RT3;
typedef Gt::Iso_cuboid_3 Iso_cuboid;
typedef Gt::Weighted_point_3 Weighted_point_3;
typedef Gt::Point_3 Point_3;
typedef P3RT3::Vertex_handle Vertex_handle;
int main(int, char**)
{
P3RT3 p3rt3(P3RT3::Iso_cuboid(0,0,0, 1,1,1));
p3rt3.insert(Weighted_point_3(Point_3(0.1,0.1,0.1), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.9,0.1,0.1), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.1,0.9,0.1), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.1,0.1,0.9), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.4,0.4,0.4), 0.001));
std::cout << "EXIT SUCCESS" << std::endl;
return 0;
}
CGAL::Regular_triangulation_cell_base_3
CGAL::Regular_triangulation_vertex_base_3

### Hidden Points

This example shows that points can be hidden during the incremental construction of a 3D regular triangulation.

File Periodic_3_triangulation_3/p3rt3_hidden_points.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_regular_triangulation_traits_3.h>
#include <CGAL/Periodic_3_regular_triangulation_3.h>
#include <iostream>
typedef CGAL::Epick K;
typedef CGAL::Periodic_3_regular_triangulation_traits_3<K> Gt;
typedef CGAL::Periodic_3_regular_triangulation_3<Gt> P3RT3;
typedef Gt::Iso_cuboid_3 Iso_cuboid;
typedef Gt::Point_3 Point_3;
typedef Gt::Weighted_point_3 Weighted_point_3;
typedef P3RT3::Vertex_handle Vertex_handle;
int main ()
{
P3RT3 p3rt3(P3RT3::Iso_cuboid(0,0,0, 1,1,1));
p3rt3.insert(Weighted_point_3(Point_3(0.1,0.1,0.1), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.9,0.1,0.1), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.1,0.91,0.1), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.1,0.1,0.89), 0.01));
p3rt3.insert(Weighted_point_3(Point_3(0.101,0.101,0.101), 0.001));
assert(p3rt3.is_valid());
std::cout << "Number of vertices : " << p3rt3.number_of_vertices() << std::endl;
std::cout << "Number of hidden points : " << p3rt3.number_of_hidden_points() << std::endl;
std::cout << "Removing the first point..." << std::endl;
p3rt3.remove(p3rt3.vertices_begin());
// The first point is removed and the hidden point is revealed.
std::cout << "Number of vertices : " << p3rt3.number_of_vertices() << std::endl;
std::cout << "Number of hidden points : " << p3rt3.number_of_hidden_points() << std::endl;
std::cout << "EXIT SUCCESS" << std::endl;
return 0;
}

### Bad Weights

Finally, this example shows how points whose weight does not satisfy the precondition are handled during the incremental construction of a 3D regular triangulation.

File Periodic_3_triangulation_3/p3rt3_insert_point_with_bad_weight.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_vertex_base_3.h>
#include <CGAL/Regular_triangulation_cell_base_3.h>
#include <CGAL/Periodic_3_regular_triangulation_traits_3.h>
#include <CGAL/Periodic_3_regular_triangulation_3.h>
#include <iostream>
typedef CGAL::Epick K;
typedef CGAL::Periodic_3_regular_triangulation_traits_3<K> Gt;
/* If remove() isn't called in our program, we can use a triangulation data structure more appropriate
* which saves some memory resources.
*/
typedef CGAL::Regular_triangulation_vertex_base_3<Gt,
CGAL::Periodic_3_triangulation_ds_vertex_base_3<> > Vb;
typedef CGAL::Regular_triangulation_cell_base_3<Gt,
CGAL::Periodic_3_triangulation_ds_cell_base_3<> > Cb;
typedef CGAL::Triangulation_data_structure_3<Vb, Cb> Tds;
typedef CGAL::Periodic_3_regular_triangulation_3<Gt, Tds> P3RT3;
typedef Gt::Iso_cuboid_3 Iso_cuboid;
typedef Gt::Weighted_point_3 Weighted_point_3;
typedef Gt::Point_3 Point_3;
int main ()
{
P3RT3 p3rt3(P3RT3::Iso_cuboid(0,0,0, 1,1,1));
// Here, we insert a point with a bad weight.
// p3rt3.insert(Weighted_point_3(Point_3(0.5,0.5,0.5), 1.));
// In debug mode, if we uncomment the previous instruction, the program displays the following error message :
// terminate called after throwing an instance of 'CGAL::Precondition_exception'
// what(): CGAL ERROR: precondition violation!
// Expr: point.weight() < ( FT(0.015625) * (domain().xmax()-domain().xmin()) * (domain().xmax()-domain().xmin()) )
// [...]
// Explanation: point.weight() < 1/64 * domain_size * domain_size
std::cout << "EXIT SUCCESS" << std::endl;
return 0;
}

## Periodic Alpha Shapes

It is possible to use the classes `Periodic_3_Delaunay_triangulation_3` and `Periodic_3_regular_triangulation_3` as underlying triangulations to compute alpha shapes. Examples of usage can be found in Section Example for Periodic Alpha Shapes of the chapter on 3D alpha shapes.

# Design and Implementation History

In 2006, Nico Kruithof started to work with Monique Teillaud on the 3D Periodic Triangulations package.

In 2007, Manuel Caroli continued work on the algorithms [1] and on the package with Monique Teillaud.

In 2015, Aymeric Pellé contributed to adding regular triangulations in the package.

The package follows the design of the 3D Triangulations package (see Chapter 3D Triangulations).