Pierre Alliez

Christophe Delage

Menelaos I Karavelas
email: mkaravel@tem.uoc.gr

Sylvain Pion

Monique Teillaud

Mariette Yvinec

Delaunay Tessellations and Voronoi Diagrams in Cgal

Introduction

Cgal is a C++ software library of computational geometry algorithms and data structures [CGAa]. Created in 1995, and mostly supported by European research projects, the Cgal library has become a mature Open Source project with approximately half a million lines of code, 10,000 downloads per year, and its availability through the two major Linux distributions Fedora and Debian.

The goal of the Cgal library is to provide easy access to efficient and reliable geometric algorithms to users in industry and academia. The initial motivation for setting up the Cgal library has its roots in the following assessment: While getting a first implementation of a geometric algorithm is a relatively easy task, obtaining a robust implementation is tremendously harder. The main reason is inherent to the dual nature of geometric objects, which combine combinatorial data structures and numerical data. Geometric algorithms have to maintain the consistency between combinatorial and numerical components. This is all the more difficult when geometric primitives are in the neighborhood of a degenerate configuration. Indeed a slight perturbation in numerical data may then yield a drastic change in the combinatorial data structure. Therefore geometric algorithms are highly sensitive to numerical rounding errors and naive implementations, using for instance built-in floating point arithmetic, are bound to fail sooner or later. Robustness is the major concern underlying the development of the Cgal library from the beginning.

More than 50 developers have contributed to Cgal since its creation. The project is managed by an international editorial board whose charge is to review all submitted packages, identify potential contributors, schedule the new releases and ensure the global consistency of the library. In addition to the open source distribution, Cgal is commercialized by the startup company GeometryFactory since 2003. This dual distribution is made available by the Open Source licenses QPL and LGPL, which allow free use of Cgal in Open Source software, in addition to commercial licenses.

The Cgal library provides a rich variety of Voronoi diagrams and Delaunay triangulations. This variety covers several aspects: generators, dimensions and metrics, which we describe in Section 2. One aim of this paper is to present the main paradigms used in Cgal: Generic programming, separation between predicates/constructions and combinatorics, and exact geometric computation (not to be confused with exact arithmetic!). The first two paradigms translate into software design choices, described in Section 4, while the last covers both robustness and efficiency issues, respectively described in Section 6 and 7. Other important aspects of the Cgal library are the interface issues, be they for traversing a tessellation, or for interoperability with other libraries or languages, see Section 5. We present in Section 8 some tessellations at work in the context of surface reconstruction and mesh generation. Section 9 is devoted to some on-going and future work on periodic triangulations (triangulations in periodic spaces), and on high-quality mesh generation with optimized tessellations. Section 10 provides typical numbers in terms of efficiency and scalability for constructing tessellations, and lists the remaining weaknesses. We conclude by listing some of our directions for the future.

Voronoi diagrams

Affine diagrams

A Voronoi diagram is defined by a set of sites and a distance function. The set of sites is a set S = {s 1 , s 2 , . . . , s n } of objects in a space Σ, and the distance function δ(x, s i ) measures the distance from a point x of Σ to a site s i . The Voronoi cell V (T ) of a subset T ⊂ S, is the set of points in Σ which are at equal distance to the sites of T and closer to any site of T than to any site in S \ T . The Voronoi diagram Vor(S, δ), simply noted Vor(S) when there is no ambiguity on the distance function, is the partition of the space Σ formed by the non-empty Voronoi cells.

Standard Voronoi diagrams are obtained in R d when the sites are points and the distance function is the Euclidean distance function, see Figure 1.

Power diagrams (also called Laguerre diagrams) are obtained when the sites are weighted points in R d × R, i.e., s i = (p i , w i ) with p i ∈ R d and w i ∈ R, and the distance function δ(x, s i ) from a point x to the site s i , is the power of x with respect to the sphere with center p i and square radius w i : δ(x, s i ) = (x -p i ) 2 -w i .

Note that the definition of power diagrams accepts negative square radii as well.

The locus of points in space Σ that are equidistant to two sites s i and s j is called hereafter a bisector. Standard Voronoi diagrams and power diagrams are also called affine diagrams (see, e.g., [START_REF] Boissonnat | Curved Voronoi diagrams[END_REF]) because their bisectors are affine subspaces of R d . The nerve of a Voronoi diagram is the set of subsets T ⊂ S with a nonempty Voronoi cell. A degenerate configuration for a Voronoi diagram in R d happens when there are points in space R d equidistant to more than d+1 sites. In the absence of degenerate situations, each subset T in the nerve of Vor(S) has cardinality at most d + 1 and can be naturally embedded as a simplex that is the convex hull conv(T ) of T . Affine diagrams share the property that, under the non-degeneracy assumption, the natural embedding of their nerve is a triangulation. Such a dual triangulation is called a Delaunay triangulation in the case of a standard Voronoi diagram and a regular triangulation in the case of a power diagram.

In the Cgal library, Voronoi diagrams and power diagrams are represented through their dual triangulations. The library provides a rich set of functionalities for Delaunay and regular triangulations in R 2 and R 3 . The triangulations are incrementally computed and can be dynamically maintained through insertion and deletion of sites. They offer localization and nearest site queries, where a localization query looks for the simplex of the triangulation including a given point and a nearest site queries seeks the site closest to a given point. A whole set of functionalities is also available to navigate through the triangulation or to explore a subset like the star or the link of a given vertex. Dual functions provide access to the features of the dual Voronoi or power diagram.

Using the standard lifting map,

(p i , w i ) ∈ R d × R -→ (p i , p 2 i -w i ) ∈ R d+1 ,
Delaunay and regular triangulations in R d can be computed as the projection of the lower part of a convex hull in R d+1 . The Cgal library offers convex hull computations in any dimension and a reduced interface to Delaunay triangulations and Voronoi diagrams in R d , including site insertion, localization and nearest site query.

Quadratic diagrams

Möbius diagrams, introduced by Boissonnat and Karavelas [START_REF] Boissonnat | On the combinatorial complexity of Euclidean Voronoi cells and convex hulls of ddimensional spheres[END_REF], generalize power diagrams and should appear soon in Cgal, at least in their planar version. Möbius diagrams are defined by sites which are doubly weighted points

s i = (p i , λ i , µ i ) ∈ R d × R × R and the distance function δ(x, s i ) from point x ∈ R d to the site s i is: δ(x, s i ) = λ i (x -p i ) 2 -µ i .
The bisectors of Möbius diagrams are hyperspheres (see Figure 2) and in fact the class of Möbius diagrams coincides exactly with the class of diagrams whose bisectors are hyperspheres [START_REF] Boissonnat | Curved Voronoi diagrams[END_REF].

In a Möbius diagram, the cells do not have to be connected and the dual of a Möbius diagram is not a triangulation of the set of sites. Even more general diagrams, whose bisectors are general conics, may be obtained when each site s i is equipped with a positive definite symmetric tensor M i such that the distance δ(x, s i ) is

δ(x, s i ) = (x -s i ) t M i (x -s i ).
These diagrams are called anisotropic Voronoi diagrams [START_REF] Labelle | Anisotropic Voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF].

Planar Euclidean diagrams

In the plane, Cgal also provides Euclidean Voronoi diagrams for extended, i.e., non-punctual, sites, namely segment Voronoi diagrams and Apollonius diagrams. In segment Voronoi diagrams, the sites are straight line segments in R 2 and the distance function δ(x, s i ) from point x to segment s i is just the Euclidean distance function:

Fig. 2. A Möbius diagram. δ(x, s i ) = min y∈si ||x -y||.
The Cgal package for segment Voronoi diagrams accepts as input intersecting segments without restriction, meaning that input segments are allowed to intersect properly, overlap or share some endpoint, see Figure 3. However, to avoid dealing with two-dimensional bisectors and disconnected Voronoi cells the package does not consider the input segments as the sites of the diagram. Instead, the sites of the diagram are related to the subsegments defined by the arrangement of the input segments. Furthermore, the two endpoints and the relative interior of each subsegment are considered as three distinct sites. In such a setting, a segment Voronoi diagram is an instance of an abstract Voronoi diagram as described by Klein [Kle89] and can be computed via a randomized incremental algorithm [START_REF] Klein | Randomized incremental construction of abstract Voronoi diagrams[END_REF]. The Cgal implementation, fully described in [START_REF] Karavelas | A robust and efficient implementation for the segment Voronoi diagram[END_REF], supports on-line insertions and nearest site queries. As an easy alternative to segment Voronoi diagrams, in the case where input segments form the boundary of a polygonal region, Cgal also provides a package to compute the straight skeleton of a polygonal region. This straight skeleton subdivides the polygonal region into subregions which can be described as the intersections of the polygonal region with the Voronoi cells of the Euclidean diagram obtained when using the supporting lines of the input segments as sites.

Apollonius diagrams, also called additively weighted Voronoi diagrams, are closely related to the Euclidean Voronoi diagrams of spheres and have obvious

s i = (p i , w i ) in R d × R and the function measuring the distance from point x ∈ R d to s i is δ(x, s i ) = ||x -s i || -w i .
Because the diagram is invariant under a uniform translation of the weights of all sites, it may be assumed that all sites have non-negative weights. In that case, the weighted point (p i , w i ) may be interpreted as a sphere with center p i and radius w i and the Apollonius distance is just the signed Euclidean distance from a point x to the sphere (p i , w i ), where the signed Euclidean distance from a point x to a sphere (p i , w i ) is the Euclidean distance affected with a negative sign when the point x lies within the ball bounded by the sphere (p i , w i ) and with a positive sign otherwise. As any other type of Voronoi diagrams, Apollonius diagrams can also be described as the projection of the lower envelope of the set of functions δ i (x) = δ(x, s i ). In the case of Apollonius diagrams, the graph of the functions δ i (x) are cones with apexes at the points (p i , w i ) of R d+1 and aperture angle 45 • . In R 2 , Apollonius diagrams have hyperbolic bisectors and they belong to the class of abstract Voronoi diagrams (see Figure 4(right)). Therefore, they can be efficiently constructed using an incremental algorithm [START_REF] Menelaos | The Voronoi diagram of planar convex objects[END_REF]. The Cgal package implements the algorithm described in [START_REF] Menelaos | The Voronoi diagram of planar convex objects[END_REF]. It supports on-line insertions and deletions, as well as nearest-site queries.

Constrained triangulations and Voronoi diagrams

The Cgal library also offers constrained triangulations and constrained Delaunay triangulations. Constrained triangulations are tessellations in which segments of a given set are prescribed to appear as edges. Constrained Delaunay triangulations are constrained triangulations that have the constrained empty circle property, a relaxed version of the empty circle property of Delaunay triangulations. The constrained empty circle property considers the set of constraint segments as obstacles to the visibility and requires that any simplex in the triangulation has a circumsphere enclosing no vertex visible from some point in the interior of the simplex. The Cgal packages for constrained triangulation and constrained Delaunay triangulation may handle intersecting, overlapping and partially overlapping input segments without restriction. When constrained segments intersect, the subsegments resulting from the arrangement of input segments are considered as the prescribed edges of the triangulation.

In some applications, it is also required that Voronoi regions do not cross some input constraint segments. Constrained Voronoi diagrams, also called bounded Voronoi diagrams (BVD) [START_REF] Seidel | Constrained Delaunay triangulations and Voronoi diagrams with obstacles[END_REF] or Voronoi diagrams with barriers [START_REF] Lingas | Voronoi diagrams with barriers and their applications. Inform[END_REF], have been defined for this purpose. A bounded Voronoi diagram is defined for a set of point sites and a set of constrained segments considered as obstacle to the visibility. The distance δ(x, s i ) between a point x of R 2 and a site s i is:

δ(x, s i ) = ||x -s i || if x is visible from s i +∞ otherwise .
Figure 5 illustrates a constrained Delaunay triangulation and its bounded Voronoi diagram.

The notion of face blindness is pivotal for constructing the bounded Voronoi diagram. A triangle is said to be blind if the triangle and its circumcenter c lie on the two different sides of a constrained edge E. Formally, is blind if and only if there exists a constrained edge E such that one can find a point p in (not an endpoint of E), such that the intersection pc ∩ E is non-empty. The BVD construction algorithm initially tags all faces of the triangulation as being blind or not blind. It then constructs each cell of the BVD using these tags.

4 Software design

Geometric traits and triangulation data structure

One of the most important paradigms used in Cgal to resolve robustness issues arising in geometric code, without dropping computational efficiency, is to establish and maintain a clear-cut separation between the combinatorial aspects of data structures/algorithms, and the numerical computations involved in these algorithms. This distinction shows up in the template parameters of Cgal classes: Triangulation classes have two template parameters called respectively the triangulation traits and the triangulation data structure.

All numerical issues are encapsulated in the triangulation traits. A model of triangulation traits provides both the geometric primitives, i.e., points, segments, triangles, tetrahedra, etc., handled by the triangulation classes and the operators required on those primitives. In the case of triangulations, the required operators are only predicates testing the sign of polynomial expressions of input points coordinates. Typically, the predicates required to build a Delaunay triangulation in R 3 are the orientation predicate which, given four points p 0 , p 1 , p 2 , p 3 in R 3 , computes the sign of the determinant orient(p 0 , p 1 , p 2 , p 3 ), and the in sphere predicate which given five points p 0 , p 1 , p 2 , p 3 , p 4 tests if p 4 is enclosed or not by the sphere circumscribed to p 0 , p 1 , p 2 , p 3 (see Section 6.1).

The triangulation data structure parameter provides the data structure to store the triangulation. In the default implementation provided by Cgal for this parameter, the representation of triangulations is mainly based on the vertices and on the full dimensional simplicial cells, i.e., triangles in 2D and tetrahedra in 3D. The triangulation data structure is therefore mainly a container of cells and vertices. The connectivity of the triangulation is encoded in cells and vertices as follows: each cell has a pointer to each of its vertices, and a pointer to each of the adjacent cells; each vertex has a single pointer to one of its incident cells. The faces with intermediate dimensions (i.e. edges in 2D, edges and facets in 3D) are only implicitly encoded through the adjacency relation between cells. Such a representation has the advantage of being compact and to generalize in arbitrary dimension.

Cgal triangulations always encode a triangulation that fills up the whole convex hull of its vertices. Furthermore, they have a fictitious vertex called the infinite vertex. The join of the infinite vertex with any facet of the convex hull is considered as a cell of the triangulation. With this convention, the set of faces encoded in a triangulation is a topological sphere and special cases for boundary facets are avoided.

Any Voronoi diagram in the plane with simply connected cells, as it is the case for abstract Voronoi diagrams and in particular for segment Voronoi diagrams and Apollonius diagrams, admit a dual planar graph called the Delaunay graph, the nodes of which are in bijection with the Voronoi sites. When there are no degenerate configurations, the Delaunay graph consists of triangular elementary cycles. In the presence of degenerate configurations the Delaunay graph may have more complex faces. Handling the degenerate cases via a perturbation scheme yields a triangulated version of the Delaunay graph by adding a few edges. The Cgal concept for two-dimensional triangulation data structure is able to handle triangulated Delaunay graphs although these planar graphs are not triangulations. There are two aspects in which triangulated Delaunay graphs may differ from a standard triangulation. First, Delaunay graphs may have multiple edges joining the same pair of sites, as well as pairs of adjacent triangles sharing two edges. Secondly, Delaunay graphs may not be embeddable on the plane via straight-line segments. See Figure 6.

The flexibility of Cgal triangulations

As any Cgal data structure, the Cgal triangulations gain their flexibility from their template parameters. Each parameter is associated to a concept, describing the set of types and functions expected from any model of this concept. The concepts are documented.

The geometric traits parameters provide the possibility to run Cgal algorithms or to build Cgal data structures, not only on the geometric primitives offered by the Cgal kernels, but also on user-defined geometric primitives provided that those primitives are included in a model, the so-called traits class. For instance, users from GIS may use the flexibility provided by triangulation traits to handle triangulated models of terrains. GIS applications often derive their models of terrains from the Delaunay triangulation of the projections of the terrain vertices on the xy-plane. In such a case, the terrain can be represented as a two-dimensional Delaunay triangulation with a triangulation traits providing three-dimensional points and predicates involving only the x and y coordinates of these points.

A great deal of flexibility also arises from the triangulation data structure parameter. In particular the concept for the triangulation data structure includes subconcepts for the models representing the cells and vertices. Therefore the user can easily change or customize the cells and vertices of the triangulation, for instance by adding to these elements colors or any other attributes needed for his/her applications.

Interfaces

Voronoi interface

The Cgal packages related to 2D Voronoi diagrams actually compute and store the dual Delaunay graph or more precisely its triangulated version. Although a Voronoi diagram and its dual are mathematically equivalent geometric objects, it is typical in applications to prefer one view over the other. To carry the point further, when computing the triangulated version of the Delaunay graph, we loose the degenerate-configuration information existent in the input. From the user's perspective, it might be easier to think of and treat the Voronoi diagram directly rather than through its dual representation. To address these issues we provide in Cgal an adaptor that adapts the 2D Delaunay graphs computed by the various Cgal packages to the corresponding 2D Voronoi diagrams seen as two-dimensional arrangements. The objective is to provide a Voronoi diagram interface that has the look and feel of a planar subdivision represented, e.g., via a Double-Connected Edge List (DCEL) data structure, although the internal representation is a triangulation data structure storing the triangulated Delaunay graph.

The adaptation is straightforward under the non-degeneracy assumption. When degeneracies are present in the input, the triangulated Delaunay graph has artifacts that might be undesirable in the Voronoi diagram view of the computed geometric structure. Suppose for example that we have a set of sites that contains subsets in degenerate position. The dual of the computed triangulated Delaunay graph is a Voronoi diagram in which all vertices have degree 3, and for that purpose we are going to call it a degree-3 Voronoi diagram, in order to distinguish it from the true Voronoi diagram of the input sites. A degree-3 Voronoi diagram can have degenerate features, namely Voronoi edges of zero length and/or Voronoi faces of zero area, that do not correspond to the true geometry of the Voronoi diagram.

The way that we treat such issues is by defining an adaptation policy. The adaptation policy is responsible for determining which features in the degree-3 Voronoi diagram are to be rejected and which are not. The policy to be used can vary depending on the application or the intended usage of the resulting Voronoi diagram. What we care about is that, firstly, the policy itself is consistent and, secondly, that the adaptation is also done in a consistent manner. The latter is the responsibility of the adaptor we provide, whereas the former is the responsibility of the implementer of a policy. We currently provide two types of adaptation policies, namely the identity policy and the degeneracy removal policy. The former adapts the triangulated Delaunay graph as is; the resulting Voronoi diagram is in fact a degree-3 Voronoi diagram, which naturally inherits the artifacts present in the triangulated version of the Voronoi diagram. This is the simplest possible adaptation policy and is the most efficient one given the triangulated Delaunay graph: the layer that the adaptor puts on top of the triangulated Delaunay graph is simply a syntactic one, since no geometric computations are performed in addition to those performed by the underlying Delaunay graph. The second type of policy is targeted towards providing to the user a view of the Voronoi diagram that corresponds to the true geometric object. This policy removes the artifacts introduced by either the additional edges introduced in order for the Delaunay graph to consist uniformly of triangles, or the splitting of the input sites into subsites (this is the case for the segment Voronoi diagram, where we split an input segment to three distinct sites, namely its two endpoints and its relative interior).

When the degeneracy removal policy is chosen, we need, within the scope of the adaptor, to identify Voronoi edges of zero length and Voronoi faces of zero area. Depending on whether we want the Voronoi diagram view to be mutable (e.g., allow on-line insertions) or non-mutable (e.g., the input is known beforehand), these tests can be quite expensive when applied to degenerate configurations, which are exactly the ones we want to eliminate. For this purpose we have two possible adaptation traits, one that caches the results of the edge/face length/area tests and one that does not. The latter policy is the preferred one under the off-line scenario, whereas the former one performs better in the on-line scenario. The drawback of the former adaptation policy is that the Delaunay graph adapted needs to support a few additional operations which may not be present.

Currently, in Cgal, we provide adaptation traits and adaptation policies for the 2D Euclidean Voronoi diagram of points, the 2D Apollonius diagram, the 2D Euclidean segment Voronoi diagram and the 2D power diagram.

BGL interface

The Boost Graph Library (BGL) provides a generic open interface for accessing a graph's structure, while hiding the implementation details [BGL]. This interface is interoperable with many graph algorithms, such that breadth-first search, depth-first search, Kruskal's and Prim's minimum spanning tree algorithms or Dijkstra's shortest paths algorithm. Although the BGL provides some general purpose graph classes conforming to the interface, they are not meant to be exhaustive or exclusive.

In that spirit, and in order to take advantage of the graph algorithms provided along with the BGL, BGL interfaces have been implemented for Cgal's triangulations, polyhedra and arrangements. For Delaunay triangulations in particular, Cgal provides an interface for the Delaunay triangulation that makes it an instance of a BGL undirected graph. The infinite edges of the Delaunay triangulation are by default present, but the user can use Boost's filtered graph or BGL's property maps in order to make them invisible. Cgal's manual page devoted to BGL provides various examples on how the users can combine the BGL interfaces provided by Cgal with functionality in the BGL.

Interfaces in other languages

Cgal is written in C++. However, to serve communities which may be used to other languages, most of Cgal packages have been interfaced with other languages.

For example, an interface in Python has been realized for many classes in Cgal: 2D and 3D Triangulations and Alpha Shapes, 2D Meshing, Polyhedron, 2D convex hulls, various geometric optimization algorithms, and a large part of Cgal's geometry kernel. It is named cgal-python [CGAb]. Python is a modern interpreted language which binds easily to C++ code.

Another example is an interface with Scilab [SCI], which is a numerical platform similar to Matlab. The interface is named CG-LAB, and is more targeting the engineers' community [CGL]. It provides interfaces to Delaunay triangulations in 2D, 3D and dD space; convex hulls in 2D and 3D; Delaunay mesh generator in 2D space; and many others.

Finally, interfaces in Java based on JNI (Java Native Interface) are also under development.

Graphical interfaces

Although not the heart of the library, Cgal is also interfaced to some visualization tools, in order to be able to interact easily with geometric data. The main toolkit used in Cgal is Qt [QT]. However, some parts of Cgal also use other toolkits for visualization. A notable one is the Ipe drawing editor [IPE], toward which many plug-ins interfacing Cgal algorithms have been implemented.

Robustness issues

The major issue with providing industrial-strength implementations of geometric algorithms is their robustness. There are two main sources of typical non-robustness issues that arise. The first is due to the approximate computations of geometric predicates and constructions using floating-point arithmetic. The second is the proper handling of all particular cases, called degenerate cases, for which we give an illustration below using symbolic perturbations.

We subdivide geometric primitives used by the higher level algorithms, in two categories: the first category is the geometric predicates, which are functions taking some geometric entities as arguments (e.g., points) and returning a result akin to a Boolean or enumerated value. The second category is the geometric constructions, which construct new geometric objects. An example of the former is the classical orientation test for three points in the plane, while an example of the latter is the construction of the circumcenter of three non-collinear points in the plane.

Robustness issues arise when, for example, the approximate computations inside predicates make them return a wrong sign, which later can draw the algorithm to an inconsistent unexpected state. Some examples of failures due to numerical inaccuracies are described in [KMP + 04], showing for example cases where the algorithm can loop in localizing a point in a 3D Delaunay triangulation. One general solution which is used throughout Cgal is to rely on the Exact Geometric Computation paradigm (EGC), formulated by Yap and Dubé in [START_REF] Yap | The exact computation paradigm[END_REF], which states that robustness of the algorithm is ensured if and only if the predicates are evaluated exactly, since the branches in the algorithms are based on the results of predicates.

From the implementation point of view, the exactness of predicates can be achieved using exact multiple precision integer, rational and floatingpoint numbers. Some wide-spread libraries, like GMP, MPFR and LEDA [GMP, MPF, LED], provide this functionality. Cgal also provides such functionality with its MP Float and Quotient classes. In turn, Cgal's predicates and constructions are parameterized by the so-called number type that they use. This way it is possible to instantiate the formulas using double, MP Float, GMP's mpq class, or other number types that are models of appropriate concepts.

An important aspect here is the study of the algebraic degree of the predicates, the goal being to minimize it, a task that is not at all easy to accomplish. There is, of course, a trade-off to be optimized here, between the number of operations required to compute a predicate and the operations accepted (e.g., usage of square roots can decrease the number of operations required to compute a predicates, but this, on the other hand, necessitates more complicated exact number types). The EGC paradigm, along with arithmetic and geometric filtering techniques (see Section 7), provides a framework for easy development of geometric software (i.e., development without the need for arithmetic considerations), and at the same time provides robustness through exactness, at a moderate additional cost (see Section 7.1). In the rest of this section we describe briefly how two predicates are computed for two of the Voronoi diagrams provided in Cgal. Our presentation is targeted toward illustrating the variety of the techniques implemented, rather than providing an exhaustive presentation.

Degenerate cases: symbolic perturbations

Computing the 3D Delaunay triangulation of a set of n points involves the computation of two major predicates: the orientation predicate and the in sphere predicate.

Let p i , p j , p k , p be four non-coplanar points in R 3 . Point p ν has coordinates (x ν , y ν , z ν ) for each ν. The orientation predicate for these four points determines whether the fourth point p lies on the positive or negative halfspace (or on the plane) with respect to the oriented plane defined by the first three points p i , p j and p k . The orientation predicate is equivalent to determining the sign of the determinant

orient(p i , p j , p k , p ) = 1 1 1 1 x i x j x k x y i y j y k y z i z j z k z .
The in sphere predicate involves five points p i , p j , p k , p and p m , and returns a positive, negative or zero value if p m lies outside, inside or on the boundary of the ball circumscribing p i , p j , p k and p . It is well known that the in sphere test can be computed in the following way:

in sphere(p i , p j , p k , p , p m ) = sign(Det(p i , p j , p k , p , p m )) sign(orient(p i , p j , p k , p )) ,
where

Det(p i , p j , p k , p , p m ) = 1 1 1 1 1 x i x j x k x x m y i y j y k y y m z i z j z k z z m t i t j t k t t m
, and t ν = x 2 ν + y 2 ν + z 2 ν for ν = i, j, k, , m. The sign(Det(p i , p j , p k , p , p m )) predicate in R 3 can be seen as an orientation predicate in R 4 , if each point p = (x, y, z) of R 3 is projected onto a point π(p) = (x, y, z, t) on the unit paraboloid Π of R 4 with equation t = x 2 +y 2 +z 2 [ES86, DMT92].

While exact predicates allow to detect degenerate configurations, the algorithms must solve special cases explicitly. For Delaunay triangulations, in 2D the situation is degenerate as soon as at least four points are cocircular. In 3D, it is degenerate as soon five points are cospherical. In the sequel, we will use the term cospherical for these two cases.

In fact, in degenerate cases, the Delaunay triangulation is not uniquely defined: any triangulation of the set of cospherical points is a Delaunay triangulation. The computed triangulation must be defined in a non-ambiguous way, to avoid inconsistencies when alternating insertions of points and vertex removals. This is performed in the Cgal package by the use of a symbolic perturbation, that allows to decide, for each case of cospherical points, which triangulation is chosen [START_REF] Devillers | Perturbations and vertex removal in a 3D Delaunay triangulation[END_REF].

The general idea of a symbolic perturbation [START_REF] Edelsbrunner | Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms[END_REF][START_REF] Seidel | The nature and meaning of perturbations in geometric computing[END_REF] is to replace the original problem by a different one depending on a parameter ε and such that:

• there exists ε 0 > 0 such that the parameterized problem is in general position for ε ∈ (0, ε 0 ] • the solution of the parametrized problem has a limit when ε goes to zero with positive values. • if the original problem is in general position, the solution of the parameterized problem tends to the solution of the original problem when ε goes to zero.

The symbolic perturbation method then considers the limit solution of the parametrized problem as the perturbed solution of the original problem. More precisely, let's assume we are facing the in sphere test of five cospherical points {p i , p j , p k , p , p m }. Our perturbation scheme needs to rank somehow the five input points. Then, assuming point p i has rank n i , the perturbation scheme adds ε ni to the fourth coordinate of the projected point π(p i ) of R 4 . Then, the orient test is not perturbed and the outcome of the in sphere test depends on the sign of the perturbed determinant

Detε(pi, pj, p k , p , pm) = ˛1 1 1 1 1 xi xj x k x xm yi yj y k y ym zi zj z k z zm ti +ε n i tj +ε n j t k +ε n k t +ε n tm +ε nm ˛.
Developing with respect to the last row yields a polynomial in ε:

Det ε (p i , p j , p k , p ,p m ) = Det(p i , p j , p k , p , p m ) + orient(p i , p j , p k , p )ε nm -orient(p i , p j , p k , p m )ε n + orient(p i , p j , p , p m )ε n k -orient(p i , p k , p , p m )ε nj + orient(p j , p k , p , p m )ε ni .
It can be shown that the polynomial Det ε can never be identically zero. Its sign is the sign of the first non-null coefficient, and depending on this sign, the in sphere predicate will tell whether the point p m must be considered as being inside or outside the tetrahedron formed by p i , p j , p k , p .

As mentioned above, the perturbation depends only on the way we rank the points. Any ordering of the points is potentially possible. In recent Cgal releases, the lexicographical ordering of the points has been implemented. This ordering has the advantage of being intrinsic to the set of points. It is not invariant under some transformations of the set of points, like e.g., rotations.

The same perturbation scheme can also be used for the case of regular triangulations [START_REF] Devillers | Perturbations and vertex removal in Delaunay and regular 3D triangulations[END_REF] and has been implemented in Cgal.

Predicates of the Apollonius diagram: Algebraic tools to the aid of geometric computations

The predicates involved in the computation of the 2D Apollonius diagram in Cgal are substantially more complicated than the ones for the 2D Voronoi diagram. It is straightforward to compute the 2D Voronoi diagram using orientation and incircle tests (much like the 3D case in Section 6.1) which are predicates of algebraic degree 3 and 4 respectively, with respect to the input coordinates. In the Apollonius diagram case six predicates are needed, which can be further decomposed to smaller primitives. In the rest of this subsection we briefly describe how one of these primitives is evaluated, namely the Radii Difference primitive (for a detailed description see [START_REF] Menelaos | Root comparison techniques applied to computing the additively weighted Voronoi diagram[END_REF] or [START_REF] Emiris | The predicates of the Apollonius diagram: algorithmic analysis and implementation[END_REF]). In particular, we show how, using algebraic tools, we can get a lower algebraic degree for the quantities the sign of which is to be evaluated, as compared to the obvious solution.

b i b j B i B j B k B ℓ c ijk c iℓj C ijk C iℓj w ijk w iℓj L ij π ij
The most complicated predicate when computing the Apollonius diagram, is the Edge Conflict Type predicate: given an edge α µν ij of the Apollonius diagram, defined by four disks B i , B j , B µ and B ν , and a fifth disk B λ , we want to determine if B λ is in conflict with one of the disks centered on α µν ij (i.e., intersects it) and tangent to both B i and B j . The Edge Conflict Type predicate can be resolved easily if we know how to order two points p and q lying on the oriented bisector π ij of B i and B j . π ij is either a line (if the radii of B i and B j are equal) or a branch of a hyperbola (if the radii of B i and B j differ). In the Apollonius diagram context, p and q are defined via two other disks B k and B , k = , (B k and B can be any of B µ , B ν and B λ ), see Figure 7. In particular, p is the center c ijk of the ccw-oriented circle C ijk tritangent to B i , B j and B k (touching them in that order) and q is the center c i j of the ccw-oriented circle C i j tritangent to B i , B and B j (again, touching them in that order).

Let L ij be the line passing through the centers of B i and B j ; if π ij is a hyperbolic branch, then L ij is its line of symmetry. If c ijk and c i j lie on different sides of L ij , we can determine their order on π ij by means of the orientation tests orient(b i , b j , c ijk ) and orient(b i , b j , c i j ), where b i and b j are the centers of B i and B j , respectively. This orientation test is of algebraic degree 14 [START_REF] Emiris | The predicates of the Apollonius diagram: algorithmic analysis and implementation[END_REF]. If c ijk and c i j lie on the same side of L ij , determining their order on π ij is equivalent to asking for the sign of the difference w ijkw i j , where w ijk and w i j are the radii of the Voronoi circles C ijk and C i j respectively, hence the name for the Radii Difference primitive. The radii w ijk and w i k can be written as:

w ijk = ρ 1 -r i , w i j = ρ 2 -r i ,
where r i is the radius of B i and, for τ = 1, 2, ρ τ > 0, is a real root of an equation of the form

f τ (x) := α τ x 2 τ -2β τ x τ + γ τ = 0.
Clearly, it suffices to compute the sign of the difference t := ρ 1 -ρ 2 . The algebraic degrees of α τ , β τ and γ τ are 4, 5 and 6, respectively. We assume for now that α 1 α 2 > 0; the case α 1 α 2 = 0 is discussed towards the end of this subsection. Let ∆ τ = β 2 τ -α τ γ τ , be the discriminant of the equation

f τ (x) = 0, and let J = α 1 β 2 -α 2 β 1 . Then t = β 1 + √ ∆ 1 α 1 - β 2 + √ ∆ 2 α 2 = -J + α 2 √ ∆ 1 -α 1 √ ∆ 2 α 1 α 2 .
Since α 1 α 2 > 0, in order to determine the sign of t, it suffices to determine the sign of the numerator of t. The numerator of t is a quantity of the form

Q 1 := X 0 + X 1 √ Y 1 + X 2 √ Y 2 , where, in general, Y 1 = Y 2 and Y 1 , Y 2 > 0.
The algebraic degrees of X 0 , X 1 , X 2 , Y 1 and Y 2 are 11, 8, 8, 6 and 6, respectively. In order to determine the sign of Q 1 , we need to determine, in the worst case, the sign of the quantity (X

2 0 + X 2 1 Y 1 -X 2 2 Y 2 ) 2 -4X 2 0 X 2 1 Y 2 ,
which is a degree 36 quantity in the input.

Let us know consider a slightly different approach. Substitute the value of ρ 1 in terms of t and ρ 2 in f 1 (x) = 0. We thus get:

α 1 (t + ρ 2 ) 2 -2β 1 (t + ρ 2 ) + γ 1 = 0,
which can be rewritten as a quadratic polynomial in terms of t: f1 (t) := ᾱ1 t 2 + β1 t + γ1 = 0, where ᾱ1 = α 1 , β1 = f 1 (ρ 2 ) and γ1 = f 1 (ρ 2 ). Determining the sign of t reduces to determining the sign of the appropriate root of f1 (t). This can be done by applying Descartes' rule of signs to f1 (t), which calls for computing the signs of β1 and γ1 . Both β1 and γ1 are expressions of the form (X 0 + X 1 √ Y )/Z, where Z > 0, and their signs can be evaluated using the following relation:

sign(X 0 + X 1 √ Y ) = sign(sign(X 0 )X 2 0 + sign(X 1 )X 2 1 Y ).
The degrees of X 0 , X 1 and Y are 9, 6 and 6 for β1 , and 14, 11 and 6 for γ1 , respectively. Hence, the highest algebraic degree involved in the evaluation of the signs of the roots of f1 (t) is 28, already an improvement over the algebraic expression of degree 36 we found above. We can further decrease the maximum algebraic degree of the expressions to be evaluated by employing Sturm sequences. Let (P i (x)) 0≤i≤M , be the Sturm sequence of f 1 (x) and f 1 (x)f 2 (x). In the absence of (algebraic) degeneracies M = 4, i.e., the Sturm sequence consists of five polynomials. In this case we can then determine the sign of t by means of the signs of the quantities J, P 3 (∞) (this is P 3 (x)/x 2 evaluated at +∞) and P 4 (P 4 (x) is a constant polynomial), where

P 3 (∞) = -(α 1 K + 2α 2 ∆ 1 ), P 4 = (P 3 (∞)) 2 (4JJ -G 2 ),
and

J = β 1 γ 2 -β 2 γ 1 , G = α 1 γ 2 -α 2 γ 1 , K = α 1 γ 2 + α 2 γ 1 -2β 1 β 2 .
The algebraic degrees of J and P 3 (∞) are 9 and 14, respectively. The algebraic degree of the second factor of P 4 is 20. The expression 4JJ -G 2 can be further factorized to two factors, one of degree 8 and one of degree 12, the latter being a polynomial of 305 monomials: the computation of these factors is very expensive. It is plausible, however, that the Sturm sequence of f 1 (x) and f 1 (x)f 2 (x) contains less than five terms. In fact the only possibility is that (P i (x)) 0≤i≤M contains four instead of five terms, i.e., M = 3, and this happens if P 3 (∞) = -(α 1 K+2α 2 ∆ 1 ) = 0. In this case P 3 (x) is the last term of the Sturm sequence, and is a constant polynomial equal to P 3 = P 3 (x) = 2α 1 ∆ 1 J . Hence the maximum algebraic degree of the quantities encountered in this case is 11.

In our discussion above it is assumed that α 1 α 2 > 0. If α 1 = α 2 = 0, we only need to test the sign of the quantity J which is of degree 11. If α 1 = 0 < α 2 (the case α 2 = 0 < α 1 is symmetric), the Sturm sequence f 1 (x) and f 1 (x)f 2 (x) ends with

P 3 = P 3 (x) = 4β 2 1 α 2 f 2 (γ 1 /(2β 1 )). The quantity 4β 2 1 f 2 (γ 1 /(2β 1 )
) is of degree 16. It is interesting to note that this quantity arises in a geometrically degenerate configuration: it corresponds to three disks having a common tangent line, or equivalently to three disks the Apollonius circle of which has infinite radius.

The degree 16 we just mentioned is the maximum algebraic degree encountered in the evaluation of the predicates of the Apollonius diagram using the Sturm sequence methodology. Recall that the degree 12 factor of the quantity 4JJ -G 2 consists of 305 monomials in the input quantities. As a result its computation involves a dramatically larger number of operations, as compared to computing the quantity 4JJ -G 2 directly. For efficiency reasons, in the current Cgal implementation we do not make use of this factorization, that is we evaluate quantities of degree up to 20.

Efficiency issues

In this section we detail several key aspects which need to be considered to obtain truly efficient implementations of the algorithms.

Arithmetic filtering

As previously mentioned, a critical part of the robustness of the implementation is achieved thanks to the exactness of the geometric predicates and constructions. We are now focusing on the details of their efficient implementation.

Exact predicates can be easily implemented using exact multiple precision integer, rational and floating-point numbers. The main issue with multiple precision number types is their efficiency. Indeed, compared to hardware assisted floating-point types, like double, for which arithmetic operations basically take a few cycles, multiple precision numbers require a memory storage which grows with the accumulation of computations (which can be exponential with rational numbers), and the cost of each operation is at least linear in the size of the operands. In practice, it is reported, e.g., in [START_REF] Devillers | Efficient exact geometric predicates for Delaunay triangulations[END_REF], that the impact of naively using these multiple precision number types, while it solves all numerical robustness issues, incurs a penalty of 2 orders of magnitude in terms of efficiency, compared to using doubles.

People then observed that, in the case of predicates, only the sign of an expression is returned. This means that only the exact sign is required to be computed, not the full bit-wise representation of an expression. Remarks have been made that, given that most of the time floating-point arithmetic returns the correct sign, it should be enough to be able to detect when such floating-point computations have a chance to lead to an incorrect sign. Then, only in these cases, rely on the exact multiple precision costly number types. This scheme is known as arithmetic filtering, since it quickly filters out the easy cases using arithmetic tools, and has been developed in a number of ways [FV96, She97, BFS01, DP03, BBP01]. Some of these are implemented in Cgal, such as dynamic filters based on interval arithmetic.

Interval arithmetic is used in Cgal as a mean to bound an exact value inside an interval whose bounds are doubles. This standard tool is relatively fast to compute with, and it is used inside Cgal's Filtered predicate class which allows to easily create filtered predicates following the scheme: first, we evaluate the predicate using interval arithmetic; then, if the sign cannot be concluded because 0 lies in the resulting interval, we recompute the predicate using an exact number type [START_REF] Brönnimann | Interval arithmetic yields efficient dynamic filters for computational geometry[END_REF]. Cgal does this using template parameterization, and instantiating successively with different number types. This ensures optimal efficiency and ease of use. This method of filtering is quite general, since all arithmetic operations are possible with intervals. The resulting code is typically 3 to 10 times slower than using uncertified doubles, depending on the algorithm (the share of combinatorial operations versus numerical operations can vary) and the data sets (more degeneracies means slower).

Besides this generic method based on intervals, some critical predicates such as those used in 2D and 3D Delaunay triangulations, benefit from additional layers of filtering, named static and semi-static filtering. The idea is to study the sequence of arithmetic operations used in a predicate, and to derive an error bound on the resulting expression we need the sign of, from bounds on the input arguments. Such a method can be much faster than interval arithmetic, since, basically, the difference between this and pure floating-point is that some bounds are computed at the beginning, which can be somehow cached, and the comparison with zero is replaced with a comparison with an error bound. In practice, this allows to reach close to optimal efficiency, as the slowdown compared to uncertified doubles tends to be about 25% [START_REF] Devillers | Efficient exact geometric predicates for Delaunay triangulations[END_REF]. These methods are however heavier to set up and error prone, which leads to formal proving of some of the methods implemented in Cgal [START_REF] Melquiond | Formal certification of arithmetic filters for geometric predicates[END_REF].

p 1 p 2 p 3 q 1 q 2 q 3 s 1 s 2
Fig. 8. Segment site representation. The point s1 is represented by the four points p1, q1, p2 and q2. The subsegment p1s1 is represented by the points p1, q1, p2, q2 and a Boolean which is set to true to indicate that the first endpoint in not an intersection point. The subsegment s1s2 is represented by the six points: p1, q1, p2, q2, p3 and q3. The remaining (non-input) points and subsegments in the figure are represented similarly.

Geometric filtering

The segment Voronoi diagram computation poses additional difficulties when we allow the input segments to intersect arbitrarily [START_REF] Karavelas | A robust and efficient implementation for the segment Voronoi diagram[END_REF]. If we represent the subsegments via the Cartesian coordinates of their endpoints, we may end up with an exponential blow up on the bit complexity of the Cartesian coordinates of intersection points. Because of this blow up we are unable to bound the algebraic degree of the predicates evaluated during execution of the algorithm, which essentially invalidates any attempt for arithmetic filtering (the arithmetic filters will almost always fail due to the high bit complexity of the geometric objects involved in the predicate evaluation).

To overcome this blow up in the Cgal package for computing the segment Voronoi diagram, we take special care on how segments are represented. More precisely, intersection points are represented by four points, namely the endpoints of the two input segments that define them, whereas segments that have intersection points as endpoints are represented either by four points and a Boolean or by six points, depending on whether only one or both endpoints of the segment are intersection points (see Figure 8). This representation allows us:

1. to have a closed-form representation of the sites defining the Voronoi diagram, thus avoiding cascading, 2. to avoid an exponential blow up on the bit complexity of the coordinates of the intersection points, 3. to use geometric filtering, discussed below.

Geometric filtering amounts to performing simple geometric tests exploiting the representation and construction history of the geometric data, as well as the geometric structure inherent in the problem, in order to evaluate predicates in seemingly degenerate configurations. Geometric filtering can be seen as a preprocessing step before performing arithmetic filtering. It can help by eliminating situations in which the arithmetic filter will fail, thus decreasing the number of times we need to evaluate a predicate using exact arithmetic.

Let us consider a simple, yet very illustrative, example of geometric filtering. Suppose we want to determine if two intersection points in the segment Voronoi diagram are identical. We assume here that the input points are represented by doubles via their Cartesian coordinates. Input segments are represented by their two endpoints, whereas intersection points and the subsegments are represented as described above. In order to determine if the two non-input points are identical we need to compute their Cartesian coordinates and compare them. If the two points are geometrically identical, the answer to our question using double arithmetic may be uncertain (due to numerical errors), in which case we will have to reside to the more expensive exact computation. Instead, before testing the coordinates for equality, we can use the representation of the points to potentially answer the question. More specifically, and this is the geometric filtering part of the computation, we can first test if the defining segments of the two points are the same. If they are not, then we proceed to comparing their coordinates as usual. Testing the defining segments for equality does not involve any arithmetic operations on the input, but rather only comparisons on doubles.

Geometric filtering techniques, such as the one presented in the paragraph above, is used in all predicate evaluations in the Cgal segment Voronoi diagram package.

Point location

Point location is crucial for the efficiency of the package. It is an important functionality in its own for the user, but it is also a major ingredient of the incremental construction of the triangulation. The Cgal 2D and 3D Delaunay triangulation packages offer two main strategies: a walk though the triangulation, and the use of a hierarchical structure.

Visibility walk.

This simple strategy is in fact quite popular for Delaunay triangulations. Let's explain it in 2D. Assume we want to locate the point p. We start the walk at an arbitrary triangle. Then, at each step, we move from a triangle t not containing p, to a neighbor of t sharing with t an edge e such that the line supporting e separates t from p. If there are two such edges, then we move to any of these two neighbors. This generalizes to higher dimensions, replacing the edges by facets.

In the case of the Delaunay triangulation, this walk has been proved to be acyclic and to reach the correct triangle [START_REF] Edelsbrunner | An acyclicity theorem for cell complexes in d dimensions[END_REF][START_REF] Floriani | On sorting triangles in a Delaunay tessellation[END_REF]. In the case of an arbitrary triangulation, the walk may loop. A variant, however, the stochastic walk, has been proved to terminate [START_REF] Devillers | Walking in a triangulation[END_REF]. This variant consists in choosing at random the neighboring triangle of t to proceed to, whenever there are more than one edge of t whose supporting line separates t from p.

Delaunay hierarchy.

Even if the walk obtains good running times, a more sophisticated strategy is proposed for large data sets. The point location uses an additional hierarchical data structure [START_REF] Devillers | The Delaunay hierarchy[END_REF]: the last level of the structure is the Delaunay triangulation DT(S), and each level i contains the Delaunay triangulation DT i of a subset of sites S i ⊂ S.

The subsets of sites S i form a decreasing sequence of random subsets of S. The number k of levels is not fixed; for each point, random trials decide its level, and the point with the highest level determines k.

A point p ∈ S such that p ∈ S i ⊆ . . . ⊆ S 0 = S and p ∈ S i+1 has a link to a Delaunay triangle of DT j incident to p for all j between 0 and i.

To locate a query point q, we start at a known vertex v k+1 of the highest level k. Using the basic location strategy in DT k we search for v k , the vertex of DT k nearest to q. Since v k is also a vertex of DT k-1 , we then start from v k to search for v k-1 , the nearest neighbor of q in DT k-1 . The search is continued descending the different levels. At each level i, the nearest vertex v i of q in DT i is determined.

It has been shown that the construction of the Delaunay hierarchy of a set of n points can be done in expected time O(n log n) in 2D and O(n 2 ) in 3D. The space is linear in the size of the triangulation, with a constant factor very close to 1. The walk at each level performs a constant expected number of steps and the expected query time is O(log n). The expectation is on the randomized sampling and the order of insertion, with no assumption on the point distribution.

In practice, instead of looking for the true nearest vertex v i of q in DT i , the point q is located in DT i using the basic location strategy, and the vertex of the triangle containing q that is closest to q is used in place of v i . This approximation invalidates the theoretical complexity analysis, but performs very well in practice.

Voronoi hierarchy.

A variant of the Delaunay hierarchy is the Voronoi hierarchy [START_REF] Menelaos | The Voronoi diagram of planar convex objects[END_REF]. Conceptually, instead of thinking of the hierarchy as consisting of a series of Delaunay triangulations, we understand the hierarchy as a series of Voronoi diagrams. The Voronoi hierarchy is a more convenient hierarchical data structure when we want to compute the Voronoi diagram for non-point objects. In fact, it is applicable to quite a few types of abstract Voronoi diagrams, thus making it more general than the Delaunay hierarchy at a moderate additional cost in query time.

In the Delaunay hierarchy the location of the nearest neighbor v i of the query point q at level i is done by walking from v i+1 (the nearest neighbor of q at level i + 1) along the line connecting v i+1 and q. In the Voronoi hierarchy we move from v i+1 to v i as follows. Let V j be the cell of the Voronoi diagram currently visited at level i (initially, V j is the cell V i+1 of v i+1 ). We look at the neighbors of v j : if there exists a neighbor v m such that δ(q, v m ) < δ(q, v j ), we visit V m and continue as before; otherwise, the distance from q to all neighbors of v j is larger than the distance δ(q, v j ), i.e., v j is the nearest neighbor v i of q at level i.

Moving from a cell V j to a cell V m takes linear time in the size of the cell. Finding the next cell to go to can be sped up by using balanced binary trees attached to the Voronoi cells: using these balanced binary trees we can find the next cell to go to in time O(log n). By means of these trees the nearest-neighbor location query can be performed in O(log 2 n) expected time, where n is the number of sites already inserted in the Voronoi diagram. In fact, we can compute the two-dimensional Euclidean Voronoi diagram for a set of n (possibly intersecting) disks, disjoint line segments or disjoint convex objects, in O(n log 2 n) expected time and O(n) expected space [START_REF] Menelaos | The Voronoi diagram of planar convex objects[END_REF]. For a set of n possibly intersecting line segments, their Voronoi diagram can be computed in O((m + n) log n) expected time and O(m + n) expected space, where m = O(n 2 ) is the complexity of the arrangement of the n line segments [START_REF] Karavelas | A robust and efficient implementation for the segment Voronoi diagram[END_REF]. As in the case of the Delaunay hierarchy, the expectation is on the randomized sampling for determining the maximum level of the hierarchy in which a site is inserted, and the order of insertion.

In reality, maintaining a balanced binary tree per cell, in a incremental or dynamic setting can be quite complicated. The current Cgal implementation of the segment Voronoi diagram and the Apollonius diagram does not maintain these trees. Although the theoretical results are no longer valid, the Voronoi hierarchy works very well in practice.

Spatial sorting

Sometimes, knowledge about a triangulation is not needed before all input points (or a large set of points) have been inserted to the triangulation. In this case, a clever ordering of the points may improve efficiency in building the triangulation. If the points are sorted along a space-filling curve, a simplex containing the previously inserted point will be a very good start for the visibility walk to locate the new point. Furthermore, part of the data structure that will be accessed during each insertion is likely to have been recently accessed. In practice, ordering points along a space-filling curve prior to insertion induces constant time visibility walk for the location of almost all new sites and, much fewer cache misses.

Cgal implements an even more subtle ordering. A set of n points is randomly subdivided in Θ(log n) rounds, where the i th round contains Θ(α i ) points, α being a constant. This subdivision is done in such a way that the probability of a point participating in round i is proportional to size of this round. Each round is then sorted along a Hilbert curve, using a quicksort-like algorithm. This random subdivision in rounds does not change the expected time complexity for Delaunay triangulations, while giving very good results in practice. See [START_REF] Amenta | Incremental constructions con BRIO[END_REF] for a comprehensive analysis. We have observed improvements of up to a factor of 2 in 3D (4 in 2D) over a random insertion order.

Spatial sorting can also be seen as a first step toward supporting better streaming representations of very large triangulations as in [START_REF] Isenburg | Streaming computation of Delaunay triangulations[END_REF].

Further scalability of Cgal algorithms

Besides all these improvements which are applicable to most Cgal algorithms, it is also possible to introduce some parallelization in the algorithms. This aspect is more and more important today, as even all last generation laptop computers are provided with multi-core processors, and the trend in future hardware is strongly following this direction. Following the practical work described in [START_REF] Daniel | Engineering a compact parallel Delaunay algorithm in 3D[END_REF], we aim at having good support for parallelization in the key Cgal algorithms for triangulations and meshing. Without going into a lot of details, this is performed by carefully splitting the work among all available processors, and managing proper locking of the data structures.

There are other techniques which are important to allow the handling of huge data sets. One such technique is the compact representation of the triangulation in memory. Cgal already provides a relatively good trade-off in this direction thanks to its Compact container class, a generic STL-like container offering the flexibility of a std::list (constant time insertions and deletions), together with a memory overhead almost as low as that of an array. More can be done in terms of compact representation and local compression [START_REF] Castelli | 2D triangulation representation using stable catalogs[END_REF], most probably by offering the user a way to reach the best balance between raw speed and memory usage for his/her application.

Tessellations at work

Alpha-shapes

Assuming a set of points S in R d and a real parameter α, the alpha-complex AC(S, α) is a subcomplex of the Delaunay triangulation DT(S) whose definition is based on the radii of empty circumspheres. A circumsphere of a simplex in DT(S) is a sphere passing through the vertices of this simplex. A circumsphere is said to be empty if it encloses no vertex of the triangulation DT(S). The characteristic property of Delaunay simplexes is to have empty circumspheres. The alpha-complex AC(S, α) includes all Delaunay simplexes which have an empty circumsphere with square radius not bigger than α. The alpha-shape is just the domain covered by simplexes of the alpha-complex. Weighted alpha-complexes and weighted alpha-shapes are defined similarly on top of regular triangulations. Recall that the sites associated to the vertices of a regular triangulation are weighted points which can be considered as spheres as well. The power product of two spheres σ 1 and σ 2 with centers c 1 and c 2 and radii r 1 and r 2 is defined as:

π(σ 1 , σ 2 ) = ||c 1 c 2 || 2 -r 2 1 -r 2 2 .
Two spheres are said to be orthogonal when their power product is zero. Instead of considering the circumspheres of simplexes, the definition of weighted alpha-shapes consider the spheres orthogonal to the sites associated with the simplex vertices. Such an orthogonal sphere is said to be empty when it has a positive power product with any site in the triangulation. The weighted alpha-complex WAC(S, α) is the subcomplex of the regular triangulation of RT(S) including all the simplexes that have an empty orthogonal sphere with square radius not bigger than α.

Alpha-shapes and weighted alpha-shapes have been introduced by Edelsbrunner and Mücke [START_REF] Edelsbrunner | Three-dimensional alpha shapes[END_REF]. Alpha-shapes have also been generalized by Kim et al. [KSK + 06] into beta-shapes. Beta-shapes are similar to weighted alphashapes, but based on the Apollonius diagram of weighted points instead of their power diagram. These notions are widely used in the context of shape reconstruction from point clouds and find natural applications in computational biology and molecular modeling. In the Cgal library, alpha-shapes are implemented as a data structure attached to the Delaunay or regular triangulations. This data structure allows, for any value of the parameter α, (1) to classify any face as included or not in the alpha-complex or (2) to retrieve the whole alpha-complex. It also provides a filtration of the triangulation, where faces are output according to the order they are included in the alpha-complex, as the parameter α is grown from -∞ to +∞.

Meshing and reconstruction

As they encode distance relationships, Voronoi diagrams, Delaunay and regular triangulations are basic tools in the field of meshing and reconstruction. The current tools offered in Cgal for mesh generation are all based upon the Delaunay refinement paradigm pioneered by Chew [START_REF] Chew | Constrained Delaunay triangulations[END_REF] and Ruppert [START_REF] Ruppert | A Delaunay refinement algorithm for quality 2-dimensional mesh generation[END_REF].

Cgal currently offers a triangle mesh generator for two-dimensional domains with constraints defined as a planar straight-line graph (PSLG), see Figure 9. Cgal also offers a module to generate surface triangle meshes and a three-dimensional mesh generator able to handle domains bounded by smooth curved surfaces or piecewise smooth curved surfaces. To approximate curved triangulation of a dense sample also yields a good approximation for normals, area, curvature and other differential properties. The required density of the sampling is relative to the so-called local feature size lfs(x), defined for each point x of the surface Σ as the distance from x to the medial axis of R 3 \ Σ. The medial axis of R 3 \ Σ is the locus of centers of balls included in R 3 \ Σ that are maximal with respect to inclusion.

The surface meshing tool of Cgal uses the notion of restricted Delaunay triangulation both to extract the approximating surface mesh from the Delaunay triangulation of the current sampling and to guide the refinement process, yielding a dense enough sampling.

A noticeable feature of this surface and volume mesh generator is that it only requires to know the surface through a geometric oracle which: (1) tells on which side of the surface lies a query point, (2) detects if a query line segment intersects the surface or not, and (3) computes the intersection points if any. Such an oracle is given to the mesh generator as a template parameter. This mechanism provides a mesh generator flexible enough to be applied in a wide variety of situations, ranging from domains defined by implicit surfaces to domains defined by level-sets in 3D gray-scaled images through point-set surfaces. Figure 11 illustrates the output of the surface mesh generator in a surface reconstruction application. An implicit function approximating the distance to the surface is first derived from the input point cloud, and the zero level set of this function is then tessellated. In Figure 12 the output of the Cgal surface mesh generator is compared with the output of the marching cubes algorithm applied to an octree. 9 On-going and future work

Periodic triangulations and meshes

The scope of computational geometry research has been mostly limited to manipulation of geometric elements in the Euclidean space R d . Other geometric spaces have hardly, if at all, been considered in the computational geometry literature so far. However, they are relevant to needs in applied domains. In particular, periodic spaces such as the torus and the cylinder in 2D or 3D are the natural spaces in a variety of situations, especially for simulations requiring periodic boundary conditions as it is often the case in fluid dynamics or astronomy.

The standard workaround consists in triangulating an elementary domain, and duplicating some fraction of the points outside that domain to simulate periodicity. This is at the same time inelegant and inefficient: it increases the size of the input set of points to be triangulated, which results in a significant overhead both in time and memory.

Some work was performed to solve theoretical problems to adapt the algorithms initially designed for R 3 to the case of the flat torus R 3 /Z 3 [CT08, CT09]. The package recently released in Cgal allows users to triangulate directly the points in the periodic space and avoid any duplication whenever possible. See Figures 13 and 14 for illustrations.

More work is in progress to compute meshes of periodic surfaces.

Optimized meshing

When high quality meshes are sought after, one approach is to resort to an optimization procedure. Two questions now arise: Which criterion should we optimize? Which degrees of freedom should we use? As the degrees of freedom are both continuous and discrete (vertex positions and mesh connectivity), there is a need for narrowing the space of possible triangulations. A class of mesh optimization techniques rely on the observation that evenly distributed points in 2D lead to well-shaped triangles [START_REF] Eppstein | Global optimization of mesh quality[END_REF]. Isotropic 2D meshing can therefore be casted into the problem of isotropic point sampling, which amounts to distributing a set of points on the input domain in as even a manner as possible. One approach to evenly distribute a set of points in 2D is to construct a centroidal Voronoi tessellation [START_REF] Du | Centroidal Voronoi tesselations: Applications and algorithms[END_REF]. Given a density function defined over a bounded domain Ω, a centroidal Voronoi tessellation (denoted CVT) of Ω is a class of Voronoi tessellations, where each site happens to coincide with the centroid (i.e., center of mass) of its Voronoi region. The centroid c i of a Voronoi region V i is calculated as:

c i = Vi x • ρ(x) dx Vi ρ(x) dx , (1)
where ρ(x) is the density function defined over V i . This structure turns out to have a surprisingly broad range of applications for numerical analysis, location optimization, optimal repartition of resources, cell growth, vector quantization, etc. This follows from the mathematical importance of its relationship with the energy function

E(Z, V) = n i=1 Vi ρ(x)||x -z i || 2 dx , (2)
where V = {V i } is a partition of Ω and Z = {z i }. It has been shown [START_REF] Du | Centroidal Voronoi tesselations: Applications and algorithms[END_REF] that (i) for a given partition V, the energy E(Z, V) is minimized when z i is the mass centroid of V i , and (ii) for a given set of centers Z, the energy function E(Z, V) is minimized when V is the Voronoi tessellation of the centers. One way to build a centroidal Voronoi tessellation is to use Lloyd's relaxation method. Lloyd's algorithm is a deterministic, fixed point iteration [START_REF] Lloyd | Least square quantization in PCM[END_REF]. Given a density function and an initial set of n sites, Lloyd's algorithm repeats the following two steps until satisfactory convergence has been achieved (see Figure 15 The Lloyd iteration can be applied to the problem of 2D quality mesh generation. Figure 16 illustrates the Lloyd iteration applied to an initial Voronoi diagram, with a variable density map derived from the local feature size of the domain boundary. The Lloyd iteration can furthermore be applied to isotropic surface remeshing. One approach is to use a global conformal planar parameterization of the input surface triangle mesh and to apply Lloyd relaxation in the parameter space using a density function designed to compensate for the area distortion due to flattening (see Figure 17).

To alleviate the numerical issues for high isoperimetric distortion, as well as the artificial cuts required for closed or models with non-trivial topology, another approach is to apply the Lloyd relaxation procedure over a set of local overlapping parameterizations [START_REF] Surazhsky | Isotropic Remeshing of Surfaces: a Local Parameterization Approach[END_REF], see Figure 18.

It has been shown that the energy minimized by the Lloyd iteration corresponds to the volume between a paraboloid and an underlaid piecewise linear approximant, which is formed by planar patches tangent to the paraboloid (the lifted Voronoi diagram) [START_REF] Chen | Robust and accurate algorithms for solving anisotropic singularities[END_REF]. In other words, isotropic point sampling can be casted as a function approximation problem in a higher dimension. In 2D, this approach leads to isotropically sampled meshes since it has been shown [START_REF] Shewchuk | What is a good linear element? interpolation, conditioning, and quality measure[END_REF] that any L p optimal approximation of a smooth function asymptotically tends to align and shape its elements according to the eigenvectors and eigenvalues of its Hessian: since the Hessian of the isotropic paraboloid function is isotropic, the resulting meshes must have nearly hexagonal Voronoi cells, i.e., nearly equilateral triangles in the dual Delaunay mesh. Unfortunately, such a property does not hold in 3D due to the presence of slivers. We can attribute the slivers to the fact that the energy tends to optimize the compactness of the dual Voronoi cells, but not the compactness of simplexes in the primal Delaunay triangulation: therefore, the presence of a sliver is not penalized by this energy. In other words, minimizing this energy ensures that the vertices in the domain are well-spaced (i.e., isotropic point Fig. 17. Isotropic remeshing of the Michelangelo David head using global planar conformal parameterization (top). An initial set of generators is randomly generated in parameter space with a variable density so as to compensate for area distortion due to the parameterization (bottom right). The Lloyd iteration is applied in parameter space (bottom middle). The output optimized mesh is uniform once lifted back in 3D (bottom left).

sampling), but having well-spaced vertices guarantees nothing in terms of the quality of the 3D mesh [START_REF] Eppstein | Global optimization of mesh quality[END_REF].

Chen [START_REF] Chen | Mesh smoothing schemes based on optimal Delaunay triangulations[END_REF] proposed an approach dual to the above in the context of mesh optimization. The main idea is to minimize a dual energy, the volume between a paraboloid and an overlaid piecewise linear approximant formed by a linear interpolation of points on the paraboloid, see Figure 19.

Casting the meshing problem as a function approximation problem in a higher dimension, we derive from this approach a variational tetrahedral meshing algorithm [START_REF] Alliez | Variational tetrahedral meshing[END_REF], which alternates between updates of connectivity and vertex positions (see Figures 20 and21). Connectivity optimization is easily achieved: for a given set of vertex positions, its Delaunay triangulation has the optimal connectivity which minimizes the dual energy. Vertex relocations are computed so as to satisfy a necessary condition for the local optimal L 1 approximation of the paraboloid.

Conclusion

We have presented the current offerings of the Cgal library in terms of Voronoi diagrams and Delaunay tessellations, as well as the main paradigms used, and the design choices made to obtain a unique combination of genericity, robustness and efficiency. In the long run, the exact geometric computation paradigm along with symbolic perturbations for handling degeneracies turns out to be a valid choice for constructing ever complex tessellations. Although the loss of performance is small when parameterizing the algorithms with exact predicates, more work remains to be done in order to elaborate upon efficient cascaded constructions for complex, scaffolded geometric data structures.

The current 3D Delaunay triangulation from Cgal inserts 80K points per second on a Intel Xeon 2.33GHz, when spatial sorting is activated. It consumes 300 MBytes per million points, which limits us to 10M points on a 32-bit computer with 3 GBytes of memory. The maximum number of points increases "only" to 30M points on a 64-bit computer with 16 GBytes, as all pointers used in the data structure consumes twice as much memory (64 vs. 32 bits). We plan to improve both scalability using size-adapted pointers and compact data structures, and efficiency by elaborating upon parallel implementations specialized to modern computer architectures with multi-core CPUs. Furthermore, our preliminary work on spatial sorting can be seen as a first step toward out-of-core implementation.

Although the variety of possible tessellations and Voronoi diagrams is huge, some are particularly relevant for applications. After variations in generators, dimension, metric and space, it is planned to generate ever nicer tessellations using variational techniques aimed at optimizing not only the connectivity and vertex locations as already presented here, but also the metric or even the generators.

Fig. 1 .

 1 Fig. 1. 2D Euclidean Voronoi diagram (left) and its dual Delaunay triangulation (right).

Fig. 3 .

 3 Fig. 3. Segment Voronoi diagrams. Left: input segments are sharing endpoints. Right: input segments are intersecting (intersection points are shown in green).

Fig. 4 .

 4 Fig. 4. Power diagram (left) and Apollonius diagram (right) for sites corresponding to the same set of circles.

Fig. 5 .

 5 Fig. 5. Constrained Delaunay triangulation (left) and its bounded Voronoi diagram (right).

Fig. 6 .

 6 Fig. 6. An Apollonius diagram (left) and its triangulated Delaunay graph shown in green (right).

Fig. 7 .

 7 Fig.7. The Radii Difference primitive: if c ijk and c i j lie on the same side of Lij, comparing their order on the oriented bisector πij reduces to comparing the radii w ijk and w i j .

Fig. 9 .

 9 Fig. 9. Triangle mesh generation by Delaunay refinement. Left: Map of Estonia provided as input PSLG (7K segments). Right: Mesh generated by adding 7.5K Steiner points (no sizing constraints were provided as input).

Fig. 10 .

 10 Fig. 10. Delaunay triangulation restricted to the input domain boundary (blue solid edges shown on the left) or the input domain (gray triangles shown on the right).

Fig. 11 .

 11 Fig. 11. Surface meshing for surface reconstruction from point clouds. Left: Input point set. Middle and Right: Implicit functions derived from the point set, and surface triangle meshes generated by Delaunay refinement and filtering for two different mesh sizing parameters.

Fig. 12 .

 12 Fig.12. Delaunay-based surface meshing vs. marching cubes in an octree. All triangles of the mesh shown on the left are better shaped compared to the triangles of the mesh to the right, and the mesh contains fewer triangles for the same approximation accuracy.

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Periodic triangulation in R 2 .

  ): 1. Construct the Voronoi tessellation corresponding to the n sites. 2. Compute the centroids of the n Voronoi regions with respect to the density function, and move the n sites to their respective centroids.

Fig. 15 .

 15 Fig. 15. Left: ordinary Voronoi tessellation. Middle left: Voronoi tessellation after one Lloyd iteration. Middle right: Voronoi tessellation after three Lloyd iterations. Right: Centroidal Voronoi tessellation obtained after convergence of the Lloyd iteration. Each generator coincides with the center of mass of its Voronoi cell.

Fig. 16 .

 16 Fig. 16. Optimized meshing. Left: Initial Voronoi diagram of a set of points uniformly sampled over the domain. Middle: Voronoi diagram after convergence of the Lloyd iteration, using a variable density function derived from the shape of the domain boundary. The distribution of dual triangle angles is shown. Right: Optimized triangle mesh.

Fig. 18 .

 18 Fig. 18. Isotropic remeshing of the Michelangelo David using overlapping parameterizations. Left: Initial triangle surface mesh (350K vertices). Right: Optimized surface mesh (250K vertices).

Fig. 19 .

 19 Fig. 19. Piecewise linear approximations: The paraboloid can be approximated by an underlaid piecewise linear function (left), or by an overlaid one (right).

Fig. 20 .

 20 Fig. 20. Tetrahedral mesh optimization. The initial mesh is obtained by uniformly sampling the interior of the torus. The mesh and distribution of radius ratios are shown at iterations 1, 2 and 50.

Fig. 21 .

 21 Fig. 21. Optimized tetrahedral meshing. Left: Variable sizing function defined inside the domain. The element size is specified in accordance with the local feature size of the domain boundary, and increases with the distance to the boundary. Middle: Cut view of the optimized tetrahedral mesh (80K vertices). All vertices are shown to highlight the increased density near the ear where the local feature size is small. Right: Another cut view of the optimized mesh.

Acknowledgments

This work has been partially supported by the IST Programme of the EU (FET Open) Project under Contract No IST-006413 (ACS -Algorithms for Complex Shapes with Certified Numerics and Topology).