Mikhail Bogdanov

Olivier Devillers

Monique Teillaud

Hyperbolic Delaunay complexes and Voronoi diagrams made practical

Keywords: Hyperbolic geometry, Poincaré ball, Space of spheres Géométrie hyperbolique, boule de Poincaré, espace des sphères

, giving a detailed description of algorithms, and presenting a static and a dynamic variants. All proofs are based on geometric reasoning, they do not resort to any use of the analytic formula of the hyperbolic distance. We also study algebraic and arithmetic issues, observing that only rational computations are needed. This allows for an exact and efficient implementation in 2D. All degenerate cases are handled. The implementation will be submitted to the CGAL editorial board for future integration into the CGAL library.

Introduction

As D. Eppstein states: "Hyperbolic viewpoint may help even for Euclidean problems" [START_REF] Eppstein | Hyperbolic geometry, Möbius transformations, and geometric optimization[END_REF]. 1 He gives two examples: the computation of 3D Delaunay complexes of sets lying in two planes [START_REF] Boissonnat | Outputsensitive construction of the Delaunay triangulation of points lying in two planes[END_REF], and optimal Möbius transformation and conformal mesh generation [START_REF] Bern | Optimal Möbius transformations for information visualization and meshing[END_REF]. Hyperbolic geometry is also used in applications like graph drawing [START_REF] Munzner | Exploring large graphs in 3D hyperbolic space[END_REF][START_REF] Eppstein | Succinct greedy graph drawing in the hyperbolic plane[END_REF].

Several years ago, we showed that the hyperbolic Delaunay complex and Voronoi diagram can easily be deduced from their Euclidean counterparts [START_REF] Devillers | The space of spheres, a geometric tool to unify duality results on Voronoi diagrams[END_REF][START_REF] Boissonnat | Outputsensitive construction of the Delaunay triangulation of points lying in two planes[END_REF]. As far as we know, this was the first time when the computation of hyperbolic Delaunay complexes and Voronoi diagrams was addressed. Since then, the topic appeared again in many publications. Onishi and Takayama write that they "rediscover the algorithm of [START_REF] Devillers | The space of spheres, a geometric tool to unify duality results on Voronoi diagrams[END_REF]", in a way that they consider as "more natural", ie. their proofs rely only on algebraic computations instead of geometric reasoning [START_REF] Onishi | Construction of Voronoi diagrams on the upper half-plane[END_REF]. Nielsen and Nock transform the computation of the Voronoi diagram in the non-conformal Klein model to the computation of an Euclidean power diagram [START_REF] Nielsen | Hyperbolic Voronoi diagrams made easy[END_REF]; however, even when the input sites have rational coordinates, the weighted points on which the power diagram is computed have algebraic coordinates. Many other references can be found in [START_REF] Tanuma | Revisiting hyperbolic Voronoi diagrams from theoretical, applied and generalized viewpoints[END_REF] (which does not mention [START_REF] Devillers | The space of spheres, a geometric tool to unify duality results on Voronoi diagrams[END_REF][START_REF] Boissonnat | Outputsensitive construction of the Delaunay triangulation of points lying in two planes[END_REF], though).

None of the above papers shows interest in practical aspects, especially algebraic and arithmetic aspects, which are well known to be crucial for exactness and efficiency of implementations.

In this paper, we stick to the Poincaré ball model of the hyperbolic space, which is conformal, ie., preserves hyperbolic angles. Due to this property, the model is used in a wide range of applications (see for instance [START_REF] Jin | Discrete surface Ricci flow[END_REF][START_REF] Zeng | Resilient routing for sensor networks using hyperbolic embedding of universal covering space[END_REF][START_REF] Hurdal | Cortical cartography using the discrete conformal approach of circle packings[END_REF]). We elaborate on our preliminary work [START_REF] Bogdanov | Hyperbolic Delaunay triangulations and Voronoi diagrams made practical[END_REF][START_REF] Devillers | The space of spheres, a geometric tool to unify duality results on Voronoi diagrams[END_REF][START_REF] Boissonnat | Outputsensitive construction of the Delaunay triangulation of points lying in two planes[END_REF], giving a detailed description of algorithms allowing to compute the hyperbolic Delaunay complex and Voronoi diagram in any dimension, either in a static or in a dynamic setting. All degenerate cases are handled. All proofs rely on purely geometric proofs, avoiding any computation and any use of the hyperbolic distance formulas. We show that only simple arithmetic computations on rational numbers are needed. The algorithm was implemented in 2D in an exact and efficient way. The implementation will soon be submitted to the CGAL editorial board for future integration into the CGAL library.

We first recall some background on the space of spheres (Section 2), Euclidean Voronoi diagrams and Delaunay triangulations. Section 3 recalls basics on hyperbolic geometry, and their interpretation in the space of spheres. In Section 4, we study hyperbolic Voronoi diagrams and Delaunay complexes, and we present algorithms. Section 5 shows, using geometric reasoning, that the computation and embedding of hyperbolic Delaunay complexes and Voronoi diagrams only use rational computations, but for Voronoi vertices whose coordinates are algebraic numbers of degree two. Section 6 presents the implementation, it gives precisions on algebraic and arithmetic aspects, and presents experimental results, in dimension 2.

2 The space of spheres E d denotes the d-dimensional Euclidean space, < ., . > the scalar product, and . the Euclidean norm.

The space of spheres states a correspondence between spheres of E d and points of E d+1 [5, Chapter 20] [START_REF] Devillers | The space of spheres, a geometric tool to unify duality results on Voronoi diagrams[END_REF]. Let χ denote the last coordinate in the space of spheres E d+1 . The direction of the χ-axis is called vertical.

An Euclidean sphere S centered at c with radius r is denoted as S = (c, r) and has equation S(x) = 0 in E d , where:

S(x) = x 2 -2 < c, x > + c 2 -r 2 . E d χ Π S0 φ(S0) * p ∈ S0 φ(p)
Figure 1: The space of spheres.

The map φ associates S to the point

φ(S) = (c, χ) ∈ E d+1 , χ = c 2 -r 2 .
We can embed E d into E d+1 by identifying it with the hyperplane χ = 0. By the embedding, φ(S) ∈ E d+1 projects vertically on E d to the center c of S. The points of E d , considered as spheres of null radius, map by φ to the unit paraboloid of E d+1

Π : χ = c 2 .
Spheres of E d map to points below Π, whereas a point above Π corresponds to an imaginary sphere, ie., a sphere whose radius is an imaginary complex number.

The pencil of spheres determined by two spheres S 1 and S 2 is the set of spheres whose equations are the affine combinations of the equations of S 1 and S 2 :

S : S(x) = α • S 1 (x) + (1 -α) • S 2 (x), α ∈ E.
This pencil is mapped by φ to the line through φ(S 1 ) and φ(S 2 ) in the space of spheres.

The set of spheres of E d orthogonal to a given sphere S 0 is represented in E d+1 by the polar hyperplane φ(S 0 ) * of point φ(S 0 ) = (c 0 , χ 0 ) with respect to Π. The equation of this polar hyperplane is obtained by polarizing the equation of Π in S 0 (see Figure 1):

φ(S 0 ) * : χ + χ 0 2 =< c 0 , c > .
In particular, a point (seen as a sphere of null radius) and a sphere are orthogonal if and only if the point lies on the sphere. The intersection of φ(S 0 ) * and Π is the image by φ of the set of points of E d lying on S 0 , ie., φ(S 0 ) * ∩ Π vertically projects on E d to S 0 . For a point p ∈ E d , the set of spheres passing through p maps to the hyperplane φ(p) * tangent to Π at point φ(p) ∈ Π.

The lower half space of E d+1 limited by φ(p) * (i.e. the half space which does not contain Π) represents the spheres of E d that enclose p in their interior open ball. In a symmetric way, the upper halfspace, denoted as Φ(p) represents the spheres that do not enclose p in their interior.

In general, for a flat A in E d+1 , we denote as A * the polar of A with respect to Π:

A * = (x, χ) ∈ E d+1 , ∀a = (x a , χ a ) ∈ A, χ + χ a 2 =< x a , x > ,
If A ⊕ B denotes the affine sum of two flats A and B, standard algebra shows that

(A ∩ B) * = A * ⊕ B * . (1)

Inria

Euclidean Voronoi diagram and Delaunay triangulation

Let P be a finite set of points in E d . For the sake of simplicity, we assume points to be in nondegenerate position. This is in fact not a restriction of our method, since we can always come down to this situation by using a symbolic perturbation scheme as in [START_REF] Devillers | Perturbations for Delaunay and weighted Delaunay 3D triangulations[END_REF] (the paper describes the 3D case but the scheme is general).

The Euclidean Voronoi diagram VD E (P) of P is the partition of

E d into Voronoi cells V E (p i ) = {x ∈ E d | ∀p j ∈ P, x -p i ≤ x -p j }.
The Euclidean Delaunay triangulation DT E (P) is the geometric dual of the Voronoi diagram. For further reading on these extensively studied data structures, see for instance [START_REF] Mark De Berg | Computational Geometry: Algorithms and Applications[END_REF][START_REF] Boissonnat | Algorithmic Geometry[END_REF][START_REF] Aurenhammer | Voronoi diagrams: A survey of a fundamental geometric data structure[END_REF][START_REF] Aurenhammer | Voronoi diagrams[END_REF].

Each cell V E (p i ) of the Voronoi diagram can be interpreted as the set of centers of spheres passing through p i and enclosing no point of P. The set of empty spheres, ie. spheres that do not enclose any point of P, is mapped by φ in the space of spheres to the intersection of the upper half spaces Φ(p) of E d+1 , p ∈ P, as defined above. The boundary of this intersection is a convex polyhedron U P , whose facets are tangent to Π. The correspondence between VD E (P) and U P is a different way of seeing the well-known duality between Voronoi diagrams and arrangements [START_REF] Edelsbrunner | Voronoi diagrams and arrangements[END_REF]:

Proposition 1 ([15]) The Voronoi diagram VD E (P) is the cell complex of dimension d in E d obtained by vertically projecting the polyhedron U P onto E d .
If σ is a k-simplex of DT E (P), we denote as P σ the set of its vertices. The dual of σ is a (dk)-face of VD E (P), which is the vertical projection of a (dk)-face u σ of U P . Any point in u σ is the image by φ of the center of a sphere passing through the vertices of σ. Thus, u σ is a convex polyhedron included in the (dk)-flat p∈Pσ φ(p) * that is the intersection of the hyperplanes dual to φ(p) for all vertices p of σ.

The k-simplex σ is incident to m (k + 1)-simplices τ 0 , τ 1 , . . . , τ m-1 of DT E (P). The set of vertices of a simplex τ i of this family is P τi = P σ ∪ {p i }, for some point p i ∈ P \ P σ . The duals of these simplices τ i determine the boundary of u σ on the (dk)-flat p∈Pσ φ(p) * . To summarize:

u σ = p∈Pσ φ(p) * 0≤i<m Φ(p i ).
(

) 2
3 The Poincaré ball model of the hyperbolic space

The hyperbolic space H d [5, Chapter 19] [33,[START_REF] Wilson | Curved Spaces[END_REF] can be represented by several widely used models. There are transformations between these models, as recalled in [START_REF] Nielsen | Hyperbolic Voronoi diagrams made easy[END_REF][START_REF] Thurston | The Geometry and Topology of Three-Manifolds[END_REF].

In In the space of spheres, the set of points at infinity H ∞ (ie. the unit sphere in E d ) is mapped to the point φ(H ∞ ) = (0, . . . , 0, -1).

Its polar hyperplane is

π ∞ = φ(H ∞ ) * : χ = 1. O H ∞ l 1 l 2 c H q p ∈ S c H S Figure 2: The Poincaré disk model of H 2 : hyperbolic lines l 1 , l 2 , hyperbolic circle S with cen- ter c H .
H d is mapped to the part of the hyperboloid Π that lies below π ∞ . For a point x ∈ H ∞ , the hyperplane φ(x) * is tangent to Π and passes through φ(H ∞ ).

Hyperbolic spheres are Euclidean spheres contained in B, but the center of a hyperbolic sphere usually does not coincide with its Euclidean center. Let us consider the pencil of spheres determined by H ∞ and a point c H ∈ B considered as a sphere of radius 0. Those spheres of the pencil that are contained in B are the hyperbolic spheres centered at c H . We denote them by S c H . They represent a collection of nested spheres growing from c H to B. In the space of spheres

E d+1 , φ(S c H ) is the half-open line segment [φ(c H ), φ(H ∞ )).
A point p is closer than point q to point c H for the hyperbolic distance if the sphere of S c H that passes through p is inside the sphere of S c H that passes through q (Figure 2). Note that we do not need to consider the explicit expression of the hyperbolic distance. See [START_REF] Beardon | The Geometry of Discrete Groups[END_REF] for its more complete description.

The intersection of the upper half spaces of E d+1 limited by the hyperplanes φ(x) * , x ∈ H ∞ forms the cone C with apex φ(H ∞ ) tangent to Π (see Figure 3). C represents the set of Euclidean spheres that do not intersect H ∞ . The set of hyperbolic spheres is the set of Euclidean spheres inside B. In the space of spheres, it is mapped to the open subset C Π of C that lies below π ∞ and below Π. C Π can also be seen as ∪ x∈B φ(S x ).

Let S ⊂ B a hyperbolic sphere mapped to φ(S) ∈ C Π in the space of spheres. We denote as ψ Π the central projection onto Π centered at φ(H ∞ ): ψ Π (φ(S)) is the intersection of line (φ(H ∞ )φ(S)) with Π. This intersection point is in fact the projection by φ on Π of the hyperbolic center c H of S. Note that hyperbolic spheres are the only spheres on which this central projection is defined: for a Euclidean sphere S that is not contained in B, the line (φ(H ∞ )φ(S)) does not intersect Π.

We denote as ψ π∞ the central projection centered at φ(H ∞ ) onto π ∞ .

Computing the hyperbolic Delaunay complex

Let P be a finite set of points in H d , represented in the Poincaré ball model. The hyperbolic Voronoi diagram VD H (P) is defined as its Euclidean counterparts, replacing the Euclidean dis-

Inria χ O Π c H φ(H ∞ ) B φ(q) q C Π π ∞ = φ(H ∞ ) * φ(S c H ) φ(q) * E d φ(S) C c E φ(c H ) = ψ Π (φ(S)) ψ π∞ (φ(S)) O π∞ B π∞ Φ(q)
Figure 3: The Poincaré ball model of H d in the space of spheres. q ∈ S, c E is the Euclidean center of S, c H is its hyperbolic center.

tance by the hyperbolic distance.

Let us take a point x in the cell V H (p i ) of p in VD H (P). x is closer to p i than to any other point in P for the hyperbolic distance. As already noted, we can avoid considering the hyperbolic distance explicitly, and express this proximity property in other words: a sphere of the pencil generated by the sphere of null radius x and H ∞ , and growing from x to the sphere H ∞ , meets p i before meeting any other point of P. This illustrates the fact that for every point x in V H (p i ) there exists a unique sphere S hyperbolically centered at x passing through p i and enclosing no point of P.

A simplex σ is an element of the hyperbolic Delaunay complex DT H (P) iff there exists a point c ∈ H d that is the hyperbolic center of a hyperbolic sphere (ie. a Euclidean sphere contained in B) passing through the vertices of σ and enclosing no point of P in its open interior ball.

An example, drawn with our implementation (see Section 6), is shown on Figure 4. Note that some k-simplices in the hyperbolic Delaunay complex may have no (k +1)-simplices in the complex incident to them. This is the reason why we don't call DT H (P) a triangulation. DT H (P) is still a simplicial complex.

Let us prove a property of DT H (P), which will not be used in the algorithm, but which is still interesting in itself. The proof is a direct consequence of the following Lemma. Proof. Let us first construct a sphere S 0 centered at the origin O by starting from the unit ball B and reducing its radius until it contains points of only one color. Without loss of generality, S 0 contains only blue points. Let p r be the red point though which S 0 passes.

Then we construct another sphere S pq , by starting from S 0 and reducing the radius, while keeping the sphere in the pencil of spheres through p r and tangent to S 0 , until it is empty. Let p b be the blue point though which S pq passes.

We have constructed a sphere S pq whose open interior ball is empty, that passes though p r and p b , and that is contained in B.

The existence of the sphere shows that for any partition of P into two sets of points, there is an edge of the hyperbolic Delaunay complex between these two sets, which proves Proposition 2.

The correspondence between VD H (P) and VD E (P) can be seen as follows, using the central projection ψ Π (see Section 3). Let f E be a k-face of VD E (P), ie. the set of Euclidean centers of empty spheres passing through a subset

P f of d + 1 -k points of P. By Proposition 1, f E is the vertical projection of a k-face f U of U P . Let c E be a point of f E , i.e.

the Euclidean center

Inria of an empty sphere S ⊂ E d passing through P f . If S is contained in B, the point obtained by centrally projecting φ(S) ∈ f U to Π from φ(H ∞ ) is ψ Π (φ(S)), which projects vertically onto the hyperbolic center c H of S (see Figure 3). This can be summarized as: Proposition 4 ( [START_REF] Devillers | The space of spheres, a geometric tool to unify duality results on Voronoi diagrams[END_REF]) The hyperbolic Voronoi diagram VD H (P) can be obtained by centrally projecting from φ(H ∞ ) the part of the polyhedron U P lying in C Π to the paraboloid Π, and projecting the result vertically onto E d .

As a consequence, DT H (P) is a subcomplex of DT E (P). It consists of faces of DT E (P) that admit at least one empty ball passing though its vertices and included in B.

The computation of DT H (P) thus consists of the following two steps:

• Compute the Euclidean Delaunay triangulation DT E (P),

• Extract from DT E (P) the simplices that also belong to DT H (P).

The rest of this section is devoted to showing how this extraction step can be performed. We first describe variants of this general scheme. The predicate that tests whether a given simplex of DT E (P) is also in DT H (P) is denoted as is hyperbolic and will be detailed in Section 4.2.

Extracting DT H (P) from DT E (P): algorithms

We give several variants of the extraction scheme. The basic one closely follows what was just presented. The second is an improvement that allows to test a smaller number of simplices. Both are static, ie., they first compute the whole Euclidean triangulation before performing the extraction. The third is dynamic: it allows to add a point and to update the hyperbolic Delaunay complex while the Euclidean Delaunay triangulation is updated.

Basic algorithm

Let us first remark that, if a simplex of DT E (P) also belongs to the hyperbolic Delaunay complex DT H (P), then all its faces also belong to DT H (P). Extracting DT H (P) from DT E (P) can be done by examining simplices by decreasing dimensions, starting from d-simplices. The extraction consists in simply marking each k-simplex of DT E (P), k = 1, . . . , d as "hyperbolic" or "nonhyperbolic".

For each dimension k, we maintain a dictionary D k of simplices to be examined. The dictionary D d initially contains all d-simplices of DT E (P), other dictionaries are empty. The algorithm proceeds by decreasing dimensions, starting at k = d.

For each k, d ≥ k ≥ 1; while D k = ∅ we pop the first k-simplex σ of DT E (P) out of D k ;
If it is marked already, don't do anything; If it is not marked yet, we test whether it is a simplex of DT H (P); If yes, we mark σ as "hyperbolic", as well as all its faces of all dimensions i from k -1 down to 1; If not, we mark σ as "non-hyperbolic". Its (k -1)-faces are not examined. We insert them into D k-1 ; end while; end for.

Anticipating on Section 4.2, in which it will be clear that predicate is hyperbolic can be evaluated with constant complexity, we get:

S σ h(σ) B H ∞ C σ h(τ ) O S τ p Sτ w σ τ v p Sσ

Improved extraction scheme

This improved scheme does not reduce the theoretical complexity given in Proposition 5, but it makes the extraction "output-sentitive" in some way: it is not sensitive to the number of simplices in the result DT H (P), but to the number of simplices that are not in the result, ie., to the size of DT E (P) \ DT H (P).

We consider that the Euclidean space E d is compactified to a topological sphere by the addition of a point at infinity. This point at infinity can be linked to all simplices of the convex hull of P. After adding these infinite simplices to it, the Euclidean triangulation DT E (P) becomes a triangulation of a combinatorial sphere.

Proposition 6

The graph G whose nodes are the d-simplices of DT E (P) \ DT H (P) and the infinite d-simplices of DT E (P), and whose arcs are adjacency relations through facets in DT E (P), is connected.

Proof. We first remark that the infinite simplices of DT E (P) form the set of all simplices that are adjacent to the infinite vertex, so, their graph is connected.

Let σ be a finite d-simplex of DT E (P) that is not in DT H (P). The sphere S σ circumscribing σ intersects H ∞ . More precisely, among the d spherical caps on S σ that are limited by the supporting (Euclidean) hyperplanes of the facets of σ and that do not contain any vertex of σ, at least one cap intersects H ∞ . If there are several such caps, we choose one that contains the point p Sσ of S σ that is the farthest to O.

Let us call C σ such a cap and h(σ) the corresponding facet of σ (See Figure 5). Any sphere S passing though the vertices of h(σ) either encloses the vertex v of σ opposite to h(σ), or intersects H ∞ . Thus, the (d -1)-simplex h(σ) does not belong to DT H (P). Moreover, when S does not enclose v, it encloses p Sσ , and its point p S farthest to O is such that Op S > Op Sσ .

Let τ be the neighbor of σ through h(σ) in DT E (P) and w the vertex of τ that is not in h(σ). Observe that τ ∈ DT H (P); otherwise h(σ) would have been in DT H (P) since h(σ) is a face of τ . Then, observe that Op Sτ > Op Sσ . A path of d-simplices can be constructed in graph G, using adjacency relations, starting at σ and ending at an infinite simplex, by choosing at each step the adjacent simplex through h(σ). This path traverses only simplices of DT E (P) \ DT H (P). Proof. In dimension 2, the definition of h(σ) appearing in the proof of Proposition 6 does not involve any choice of a good spherical cap, since in two dimensions the three caps are disjoint and only one intersects H ∞ . The mapping h has been proved to be injective. It remains to prove that any edge e ∈ DT E (P) \ DT H (P) is the image by h of a triangle σ. Let S be a circle circumscribing e. S intersects H ∞ , otherwise e ∈ DT H (P). Then we consider circles of the pencil of circles passing through the vertices of e, starting from from S, and going in the direction that decreases the intersection with H ∞ , until we find a vertex v of P. Vertex v and edge e form the Euclidean Delaunay triangle σ such that h(σ) = f . The improved version of the extraction algorithm consists in "digging" DT E (P): its starts from the infinite d-simplices of DT E (P), and recursively traverses G using adjacency relations. The recursive traversal stops digging as soon as it can only reach d-simplices on which is hyperbolic is true. During the traversal, faces of all d-simplices of G are tested against is hyperbolic and marked accordingly as in the basic algorithm.

Inria v σ i B σ if S i

Dynamic variant

Let us quickly recall the following definitions: The star St X (v) of a vertex v in a simplicial complex X is defined to be the subcomplex consisting of all cofaces of v, ie., all the simplices of X that contain v and their faces. The link Lk X (v) of v is the subcomplex of X consisting of all faces of St(v) that do not contain v.

Lemma 8 If all the (d -1)-simplices of Lk DT E (P) (v) belong to DT H (P), then St DT E (P) (v) is a subcomplex of DT H (P).
Proof. Let σ 0 , σ 1 , . . . , σ m-1 be the d-simplices in St DT E (P) (v), and σ if , 0 ≤ i < m be the facet of σ i opposite to v. The hypothesis means that there is an empty ball B i included in B whose boundary S i passes through the vertices of σ if (See Figure 6). The circumscribing ball of each

σ j , 0 ≤ j < m, is included in the union ∪ m-1 i=0 (B i ∪ σ i ), so, it is also included in B. B S H σ i-1 p σ S E σ i S Figure 7: Proof of Lemma 9.
This variant consists in using a dynamic algorithm, allowing the Euclidean Delaunay triangulation to be updated at each insertion of a new point. At the same time, we can update the hyperbolic complex in an efficient way, by updating the marks "hyperbolic" or "non-hyperbolic". Let p be the i th point, and P i-1 be the set of i -1 points inserted before. The insertion of p removes a set of d-simplices from DT E (P i-1 ), whose union forms a topological ball, and whose boundary is the link of p in DT E (P i-1 ∪ {p}). These simplices are replaced in by simplices of the star of p in DT E (P i-1 ∪ {p}). Let us consider the set Σ of d-simplices of DT E (P i-1 ∪ {p}) that lie outside this topological ball and that are adjacent through a facet to some d-simplex in this star.

• If, for any d-simplex σ ∈ Σ, σ is an element of in DT H (P i-1 ), ie., if it is marked as hyperbolic in DT E (P i-1 ), then all simplices in the link of p are hyperbolic, since they are facets of hyperbolic simplices. So, by Lemma 8, all new simplices of DT E (P i-1 ∪ {p}) created by the insertion of p are in DT H (P i-1 ∪ {p}). In such a good case, there is no need to run is hyperbolic on any simplex.

• Otherwise, we check new simplices and their faces with is hyperbolic, starting from dsimplices and going down in dimensions, as in the basic algorithm. Note that faces of new simplices must be checked, even if they are not new, since their mark may need to be updated.

To show the correctness of this algorithm, it remains to show the following:

Lemma 9 A k-simplex σ that was in DT H (P i-1 ) and stays in DT E (P i-1 ∪ {p}) can become non-hyperbolic after the insertion of p only if σ is in the link of p in DT E (P i-1 ∪ {p}).

Proof. Let σ be a k-simplex of DT H (P i-1 ) that stays in DT E (P i-1 ∪ {p}). Let us assume that σ does not belong to DT H (P i-1 ∪ {p}). Then there was at least one empty sphere S H σ i-1 passing through the vertices of σ and included in B before p was inserted, and there is no such sphere any more afterwards. So, S H σ i-1 encloses p in its open interior ball (see Figure 7). The simplex σ stays in DT E (P i-1 ∪ {p}), so, there is at least one empty (not containing any point of P i-1 ∪ {p}) sphere S E σ i passing though the vertices of σ. Let us consider the pencil of spheres generated by S H σ i-1 and S E σ i . There is a sphere in this pencil that passes through p. This sphere S is included in the union of the two balls interior to S H σ i-1 and S E σ i , so, it is is empty. S passes though all the vertices of σ ∪ {p}. So, σ ∪ {p} is a (k + 1)-simplex in DT E (P i-1 ∪ {p}), and σ is an element of the link of p in this triangulation.

Inria H ∞ φ(H ∞ ) H d E d π ∞ Π C Π u σ n σ O π∞ B π∞
So, simplices outside the link of p in DT E (P i-1 ∪ {p}) don't need to be tested.

Altogether, this discussion shows that the dynamic variant to compute the hyperbolic Delaunay triangulation has the same complexity as the dynamic algorithm to compute the Euclidean triangulation. Using an appropriate point location data structure [START_REF] Devillers | The Delaunay hierarchy[END_REF], we get:

Extracting DT H (P) from DT E (P): predicate

Let us now explain the test is hyperbolic, which checks whether a simplex of DT E (P) is in DT H (P). Let observe right now that this predicate has no degenerate case: as mentioned in Section 3 (page 6), hyperbolic spheres are Euclidean spheres contained in the open ball B. The only candidate for a degenerate case would be the limit case when the only empty sphere passing through the vertices of a given simplex is tangent to ∂B = H ∞ . Then the simplex is just not hyperbolic, and the case is in fact not degenerate.

We first look at d-simplices. Let σ be a d-simplex of DT E (P), dual to a Voronoi vertex, projected from the vertex u σ of U P in E d+1 . u σ is the image by φ of the (empty) sphere circumscribing σ. From Section 3, σ is a d-simplex of DT H (P) iff u σ lies in C Π . Since P ⊂ B, the circumscribing sphere of σ cannot completely lie outside B, so, this equivalence can be rewritten as: σ is a d-simplex of DT H (P) iff u σ lies in C.

For a general dimension k, the discussion at the end of Section 3 straightforwardly shows that the following conditions are equivalent:

1. the k-simplex σ, k ≤ d of DT E (P) is a simplex of DT H (P) RR n°8146 B π∞ ψ π∞ (u σ ) n σ O π∞ (d -k)-flat supporting ψ π∞ (u σ ) in π ∞ Figure 9: The central projection of u σ on π ∞ intersects the open unit d-ball B π∞ in π ∞ . 2. u σ intersects C 3. ψ Π (u σ ) = ∅ 4. ψ π∞ (u σ ) intersects the open unit d-ball B π∞ of π ∞ .
See Figure 8 for an illustration.

If σ is a k-simplex of DT H (P) whose incident (k + 1)-simplices do not belong to DT H (P), the central projection ψ π∞ (u σ ) of u σ onto π ∞ is a convex (dk)-polyhedron whose facets lie outside the unit d-ball B π∞ in π ∞ (See Figure 9). The polyhedron ψ π∞ (u σ ) intersects B π∞ , but its boundary does not intersect B π∞ , so, the whole intersection of the supporting flat of ψ π∞ (u σ ) with B π∞ is contained in ψ π∞ (u σ ). So, it is sufficient to test any well chosen point in the intersection of the supporting flat and the ball B π∞ . In hyperplane π ∞ , let n σ be the point of the supporting (dk)-flat of ψ π∞ (u σ ) that is the nearest to the origin O π∞ of B π∞ . Predicate is hyperbolic will check whether ψ π∞ (u σ ) intersects B π∞ , ie., is hyperbolic(σ) n σ is contained in both the convex polyhedron ψ π∞ (u σ ) and the ball B π∞ .

Geometric proofs for algebraic aspects

In this section we first recall how the Euclidean Delaunay triangulation DT E (P) can be computed only using rational computations. Then we show in Section 5.2 that extracting the hyperbolic Delaunay complex DT H (P) from DT E (P) can also be performed only using rational computations. In Section 5.3, we consider the geometric embedding of DT H (P) in the Poincaré ball model. The algebraic aspects of the geometric embedding of VD H (P) are studied in Section 5.4.

We assume that the coordinates (x 0 , . . . , x d-1 ) of each point x of P are rational.

Computing DT E (P)

This section quickly recalls basic facts, we refer the reader to the mentioned literature for more details. Many standard algorithms have been proposed to compute DT E (P). The CGAL library Inria [1] offers an efficient implementation in 2D and 3D [START_REF] Yvinec | 2D triangulations[END_REF][START_REF] Pion | 3D triangulations[END_REF], based on the incremental construction first proposed by Bowyer [START_REF] Bowyer | Computing Dirichlet tesselations[END_REF] and Watson [START_REF] Watson | Computing the n-dimensional Delaunay tesselation with applications to Voronoi polytopes[END_REF]. Several options are proposed for point location, using either some walking strategy [START_REF] Devillers | Walking in a triangulation[END_REF] or a hierarchical data structure [START_REF] Devillers | The Delaunay hierarchy[END_REF].

The robustness of the implementation against arithmetic issues is ensured by following the exact geometric paradigm [START_REF] Yap | The exact computation paradigm[END_REF]. The computation of the combinatorial structure underlying DT E (P) only relies on the evaluation of two predicates.

The predicate orientation(p 0 , p 1 , . . . , p d ) decides the orientation of d + 1 points, and boils down to the sign of the determinant:

p0 0 p0 1 . . . p0 d-1 1 p1 0 p1 1 . . . p1 d-1 1 . . . . . . . . . . . . . . . p d 0 p d 1 . . . p d d-1 1
The predicate in sphere(p 0 , p 1 , .., p d , q) decides whether q lies in the open interior ball of the sphere passing through the d + 1 first points. If we assume that these d + 1 first points are positively oriented, then the predicate is given by the sign of the determinant:

p0 0 p0 1 . . . p0 d-1 p0 2 1 p1 0 p1 1 . . . p1 d-1 p1 2 1 . . . . . . . . . . . . . . . . . . p d 0 p d 1 . . . p d d-1 p d 2 1 q0 q1 . . . q d-1 q 2 1
, which is in fact exactly the orientation predicate of (φ(p 0 ), φ(p 1 ), . . . , φ(p d ), φ(q)) in the space of spheres.

Only rational operations are needed to evaluate these signs of polynomial expressions. The computation is made both exact and fast by using filtered exact computations [START_REF] Devillers | Efficient exact geometric predicates for Delaunay triangulations[END_REF].

Extracting DT H (P) from DT E (P)

Lemma 11 Each rational point of E d+1 is projected by ψ π∞ on π ∞ to a rational point.

Proof. A point (c, χ) ∈ E d+1 is projected to π ∞ following the line through φ(H ∞ ) to point (x, 1), x ∈ E d , such that x = 2 1 + χ • c.
Note that 1 + χ is the power product of the unit sphere H ∞ and the sphere S such that φ(S) = (c, χ).

Lemma 12

The sphere S circumscribing a d-simplex σ whose vertices are rational points is mapped by φ to a rational point.

Proof. φ(S) = p∈Pσ φ(p) * . Each φ(p) * is a rational hyperplane of E d+1 , so, the intersection is rational.

Let us remark that a sphere S ⊂ E d has rational Euclidean circumcenter and squared radius iff its associated point φ(S) in the space of spheres has rational coordinates. The equation of cone C in E d+1 is given as follows

(x, χ) ∈ C ⇐⇒ x 2 - 1 + χ 2 2 < 0. ( 3
)

RR n°8146

Proposition 13 The evaluation of the predicate is hyperbolic, which tests whether a k-simplex of DT E (P) belongs to DT H (P), can be performed using rational computations only.

Proof. Let σ be a k-simplex in DT E (P) and {τ i , 0 ≤ i < m} be the collection of its incident (k+1)-simplices in DT E (P). Then the corresponding (d-k)-face u σ of U P is given by Equation (2) (see Section 2). Thus, the construction of u σ involves rational computations only. As shown by Lemma 11, the central projection ψ π∞ (u σ ) of u σ onto π ∞ is rational as well.

The projection n σ of the origin O π∞ onto the (dk)-flat ψ π∞ p∈Pσ φ(p) * supporting ψ π∞ (u σ ), as defined in Section 4.2, is given by the intersection of the normal to this (dk)-flat passing through O π∞ and itself, thus it is rational.

It remains to say that the test whether n σ lies in the convex polyhedron given by Equation ( 2) boils down to orientation tests, and inclusion in B π∞ reduces to comparing its square distance to O π∞ with 1.

Altogether, Sections 5.1 and 5.2 show that the combinatorial structure of DT H (P) can be computed using only rational computations.

Computing the embedding of DT H (P) in the Poincaré ball

Let us now focus on the geometric embedding of DT H (P).

As already mentioned, a hyperplane (ie., a (d -1)-flat) in H d is a portion of an Euclidean d-sphere. We say that the hyperplane is rational if the corresponding Euclidean sphere as a rational equation, ie., if its Euclidean center has rational coordinates and its squared Euclidean radius is rational, or equivalently, if it is mapped by φ to a rational point.

A k-flat in H d , for k > 0, is given by the intersection of dk hyperplanes. We say that the k-flat is rational if all these hyperplanes are rational. Then we inductively define a k-simplex to be rational if its supporting k-flat is rational and if its faces of dimensions 0, . . . , k -1 are rational.

Proposition 14

In the Poincaré ball model, the geometric embedding of the hyperplane supporting any facet of DT H (P) is rational.

Proof. A facet of DT H (P) is a (d-1)-simplex supported by the hyperplane containing its vertices p 0 , p 1 , . . . , p d-1 . This hyperplane is embedded in the Poincaré ball model as the Euclidean sphere S that passes through the p i 's and that is orthogonal to H ∞ , ie.,

φ(S) = 0≤i<d-1 φ(p i ) * π ∞ .
It is the intersection of hyperplanes of E d+1 , which are all rational.

Corollary 15 The embedding of DT H (P) is rational.

Proof. By hypothesis on P, each vertex of DT H (P) has rational coordinates. Each edge (1simplex) of DT H (P) is supported by the intersection of d -1 hyperplanes, and its endpoints are rational vertices. The edge will be rational iff all hyperplanes are rational. In the same way, a k-simplex DT H (P), for any k > 0, will be rational iff all d-k hyperplanes defining its supporting flat are rational. All considered hyperplanes are rational from Proposition 14, which concludes the proof.

Inria H ∞ φ(H ∞ ) p q φ(p) φ(q) H 2 E 2 π ∞ Π bisector of p and q H * O φ(p) ⊕ φ(q)
Figure 10: H * is the image by φ of the bisector of p and q in H 2 .

Computing the embedding of VD H (P) in the Poincaré ball

Let us move on to the geometrical embedding of the hyperbolic Voronoi diagram VD H (P).

Proposition 16

The bisector of two points of P is a hyperbolic hyperplane whose equation in E d is rational.

Proof. Let p and q be two points of P. In the space of spheres, We are going to construct their hyperbolic bisector as the locus of hyperbolic centers of spheres passing through both of them.

The intersection φ(p) * ∩ φ(q) * is the (d -1)-flat in the space of spheres E d+1 that represent all spheres of E d passing through p and q. By the construction explained earlier for hyperbolic centers, ψ Π (φ(p) * ∩φ(q) * ) is the image by φ of the set of all centers of these spheres. By definition of ψ Π , it is the intersection of the hyperplane H = (φ(p) * ∩ φ(q) * ) ⊕ φ(H ∞ ) with Π.

The polar of this hyperplane is the point H * = (φ(p) ⊕ φ(q)) ∩ π ∞ , using Equation 1. Consequently, H ∩ Π is the image by φ of a sphere in E d , which is the set of centers of spheres through p and q. All steps in the construction involve only rational computations. Figure 10 illustrates it for d = 2.

Since a k-face in VD H (P) is the bisector of k points of P, it is the intersection of k -1 rational hyperplanes, we deduce that Corollary 17 The bisector of k points, for 2 ≤ k ≤ d, is a hyperbolic (dk + 1)-flat whose equation is rational.

For k = d + 1, this means that the equation of a Voronoi vertex, seen as a Euclidean 0-sphere, is rational. However: Proposition 18 The coordinates of a hyperbolic Voronoi vertex are algebraic numbers of degree 2.

Proof. There are at least two ways of seeing this result. One is to consider a hyperbolic Voronoi vertex as the intersection of d hyperplanes in H d , ie., d Euclidean spheres S 0 , . . . , S d-1 . This is (0, -1)

(0, 0) (c, c 2 -r 2 ) (x, x 2 ) c 2 -r 2 + 1 1+ x 2 E d χ (x, 0) (c, 0)
Figure 11: Computation of the hyperbolic center x of a sphere (c, r). also the intersection of S 0 with the d -1 radical hyperplanes of S 0 and S i , 1 ≤ i < d. Each radical hyperplane is rational, since its equation is obtained as the difference of the equations of the two corresponding spheres. The Voronoi vertex is the intersection point between a sphere and a line that lies in B, which shows the result.

A direct construction of the hyperbolic center x of a sphere of Euclidean center c and radius r, using ψ Π (see Figures 3 and11), shows that x is the smallest solution of

x = 1+ x 2 1+c 2 -r 2 • c , then, x = 1+ x 2 1 + c 2 -r 2 • c.
The intersection of an infinite k-face of VD H (P) and H ∞ is a rational sphere. For k = 1, the intersection is a Euclidean 0-sphere on H ∞ . Although, the equation of the 0-sphere is rational, the two points of the sphere have coordinates which are algebraic numbers of degree 2. To show this, it suffices to repeat the proof of Proposition 18 considering a point at infinity as the intersection of the sphere H ∞ and d -1 Euclidean spheres.

Implementation in H 2

The construction of the Delaunay complex and the Voronoi diagram in H 2 was implemented using CGAL [1]. The implementation will soon be submitted for future integration in CGAL. In this section, we first present the implementation, then we detail the computations used. We finally give some benchmarks.

Software design

CGAL library provides a package to compute the Delaunay triangulation and the Voronoi diagram in E 2 [START_REF] Yvinec | 2D triangulations[END_REF]. First we give an introduction to the design of the package (a general presentation Inria of the CGAL design can be found in [START_REF] Fogel | Generic programming and the CGAL library[END_REF]). Then we discuss the design of our implementation for the hyperbolic case. We consider our implementation in terms of changes relative to the CGAL package.

CGAL 2D triangulations

The main class of the algorithm computing the Delaunay triangulation of E 2 is CGAL::Delaunay triangulation 2<Gt, Tds>. This class has two template parameters Gt, providing basic geometric computations, and Tds, handling combinatorical operations.

The geometric traits class Gt must provide a set of geometric objects, predicates and constructions that are required by the algorithm. Computing a Delaunay triangulation requires only predicates, namely Orientation 2 and In circle 2 are the two main predicates. The computation of the Voronoi diagram requires geometric constructions such as Construct circumcenter 2. The CGAL kernels provide all necessary predicates and constructions to compute Delaunay triangulation and Voronoi diagrams in E 2 and can be used as traits class. The number type used for the computations can be chosen. A number type ensuring that predicates are evaluated exactly is required for the robust computation of the Delaunay triangulation. Constructions are only required for the Voronoi diagram and can be left approximate.

The second template parameter of the triangulation class is the triangulation data structure Tds, which maintains the combinatorial structure of the triangulation. It stores the vertices and the faces. A face stores three pointers to its vertices indexed from 0 to 2 and three pointers to its adjacent faces. A vertex stores a pointer to one of its incident faces. An edge is implicitly represented as a pair of an incident face and the index of the vertex that is opposite to the face. The class CGAL::Triangulation data structure 2 provides a default implementation of Tds. It has two template parameters: a vertex class and a face class.

The 2D hyperbolic complex

Now we describe our implementation of the algorithm for computing the Delaunay complex in H 2 . We discuss the design as a whole. Then we turn to some specific aspects of the implementation.

Design

The class Delaunay hyperbolic triangulation 2<Gt, Tds> is derived from CGAL::Delaunay triangulation 2<Gt, Tds>. It mainly consists of an implementation of extracting schemes of the hyperbolic complex from the Euclidean triangulation. Modifications in the classes used as template parameters are quickly mentioned below. The geometric traits Gt. We introduce class Triangulation hyperbolic traits 2<Kernel> that provides the new predicate is hyperbolic, as well as all geometric constructions (segments, circumcenters, bisectors, etc) in the hyperbolic plane. A CGAL kernel can be used as template parameter.

The triangulation data structure Tds. The same class CGAL::Triangulation data structure 2 can be reused, enriching it with marks "hyperbolic" and "non-hyperbolic". Vertices of DT E (P) are always in DT H (P), so, they do not need to be marked. To be able to mark faces, we replace the default face class CGAL::Triangulation face base 2 by a new class Triangulation hyperbolic face base 2. The new class inherits from CGAL::Triangulation face base 2 and maintains additional information. A Boolean is used to mark a face of DT E (P) as "hyperbolic" or "nonhyperbolic". To mark edges, we use the one-to-one correspondence between edges and faces of DT E (P)\DT H (P) defined in Corollary 7. Since each edge of a face is accessed through the index of the opposite vertex in the face, then it suffices to store also the index of the vertex opposite to the "non-hyperbolic" edge.

Extracting schemes

We implemented the two extracting schemes of the hyperbolic complex from the Euclidean triangulation described in Section 4. if each face of Σ is already marked ''hyperbolic'', then mark faces in star(v) as ''hyperbolic'' else run is hyperbolic on each face and edge of star(v); return v; }

Algebraic and arithmetic aspects

Let us now detail the computations that are needed when implementing the algorithm.

Predicates

As already mentioned in Section 5.1, we rely on the CGAL package [START_REF] Yvinec | 2D triangulations[END_REF] to exactly compute the Euclidean Delaunay triangulation DT E (P). The hyperbolic complex DT H (P) is then extracted from DT E (P) using the predicate is hyperbolic detailed in Sections 4.2 and 5.2. This predicate is called on Euclidean Delaunay triangles. Non-hyperbolic edges are deduced from non-hyperbolic triangles using the bijection introduced in Corollary 7.

Triangles

Let pp p be a triangle of DT E (P) and S pp p its circumscribing circle. As noticed in Lemma 12,

φ(S pp p ) = φ(p) * ∩ φ(p ) * ∩ φ(p ) * . The equation of φ(p) * is φ(p) * : 2p 0 x 0 + 2p 1 x 1 -χ = p 0 2 + p 1 2 (4)
where (x 0 , x 1 , χ) is a point of E 3 .

Inria

φ(S pp p ) = 1 2p0 2p1 -1 2p 0 2p 1 -1 2p 0 2p 1 -1   p0 2 + p1 2 2p1 -1 p 0 2 + p 1 2 2p 1 -1 p 0 2 + p 1 2 2p 1 -1 , 2p0 p0 2 + p1 2 -1 2p 0 p 0 2 + p 1 2 -1 2p 0 p 0 2 + p 1 2 -1 , 2p0 2p1 p0 2 + p1 2 2p 0 2p 1 p 0 2 + p 1 2 2p 0 2p 1 p 0 2 + p 1 2   (5)
Plugging this coordinates in Equation ( 3) describing the cone C we get that S pp p is hyperbolic if

p 0 2 + p 1 2 2p 1 -1 p 0 2 + p 1 2 2p 1 -1 p 0 2 + p 1 2 2p 1 -1 2 + 2p 0 p 0 2 + p 1 2 -1 2p 0 p 0 2 + p 1 2 -1 2p 0 p 0 2 + p 1 2 -1 2 - 1 4   2p 0 2p 1 -1 2p 0 2p 1 -1 2p 0 2p 1 -1 + 2p 0 2p 1 p 0 2 + p 1 2 2p 0 2p 1 p 0 2 + p 1 2 2p 0 2p 1 p 0 2 + p 1 2   2 < 0
that is

p0 2 + p1 2 p1 1 p 0 2 + p 1 2 p 1 1 p 0 2 + p 1 2 p 1 1 2 + p0 p0 2 + p1 2 1 p 0 p 0 2 + p 1 2 1 p 0 p 0 2 + p 1 2 1 2 - p0 p1 p0 2 + p1 2 -1 p 0 p 1 p 0 2 + p 1 2 -1 p 0 p 1 p 0 2 + p 1 2 -1 2 < 0
Thus we get the following lemma:

Lemma 19 In H 2 , the predicate is hyperbolic applied to a triangle can be evaluated as the sign of a polynomial of degree 8 in the coordinates of its vertices.

Edges

v = p -c E , v = p -c E , v = p -c E , which verify v = v = v and v, v , v are in counterclockwise order.
Let us define v c = --→ Oc E . The vectors v, v , v split the circle S pp p into three (2D) spherical caps. The ray of origin c E and direction v c intersects S pp p in p S pp p , which lies in one of these three caps, and outside B . The cap that contains p S pp p is supported by h(pp p ) and can be uniquely determined by the two vectors from v, v , v that define a wedge in counterclockwise order containing v c (See Figure12).

The computation of h(pp p ) boils down to computing Counterclockwise( v, v c , v ), Counterclockwise( v , v c , v ), and Counterclockwise( v , v c , v), where Counterclockwise( w, v c , w ) tests whether the three vectors w, v c , and w are in counterclockwise order. 2It remains to detail the computation of Counterclockwise( w, v c , w ), where w = w (notice that v c has no reason to be equal to w ). Let v c ⊥ be the vector obtained by rotating v c by

π 2 . If sign( w • v c ⊥ ) = -sign( w • v c ⊥ ) = -1 (
which holds e.g. for v and v in Figure 12-left) then

Counterclockwise( w, v c , w ) is true. If sign( w • v c ⊥ ) = sign( w • v c ⊥ ) = -1 (see Figure 12-right)
then Counterclockwise( w, v c , w ) is equal to Orientation( w, v c , w ). Other cases are symmetric.

From Equation ( 5)

c E = 1 2p0 2p1 -1 2p 0 2p 1 -1 2p 0 2p 1 -1   p0 2 + p1 2 2p1 -1 p 0 2 + p 1 2 2p 1 -1 p 0 2 + p 1 2 2p 1 -1 , 2p0 p0 2 + p1 2 -1 2p 0 p 0 2 + p 1 2 -1 2p 0 p 0 2 + p 1 2 -1   .
So, the coordinates of v, v , v , and v c are rational fractions with numerators of degree 3 and a common denominator of degree 2. Thus the signs of the above scalar products and Orientation tests boil down to signs of polynomials of degree 6, and we get the following lemma: Lemma 20 The non hyperbolic edge associated to a non hyperbolic triangle by the map h can be determined by the evaluation of the signs of polynomials of degree 6 in the coordinates of its vertices.

Exact evaluation

We have seen that all predicate evaluations boil down to computing signs of polynomials. As for the CGAL Euclidean Delaunay triangulations (section 5.1), this can be done in a fast and exact way using filtered exact computations, providing an efficient and fully robust implementation.

Constructions

To draw the Delaunay triangulation we need to construct the hyperbolic line through two points, and for and the Voronoi diagram we need to construct the hyperbolic bisector of two points and the hyperbolic center of the circle through three points.

Hyperbolic line

The hyperbolic line through p and p is supported by the Euclidean circle S through p, p , and orthogonal to B:

φ(S) = φ(p) * ∩ φ(p ) * ∩ π∞ = 1 p 0 p 1 -1 p 0 p 1 -1 0 0 1   p 0 2 + p 1 2 p 1 -1 p 0 2 + p 1 2 p 1 -1 1 0 1 , p 0 p 0 2 + p 1 2 -1 p 0 p 0 2 + p 1 2 -1 0 1 1 , p 0 p 1 p 0 2 + p 1 2 p 0 p 1 p 0 2 + p 1 2 0 0 1   =     - p 1 p 0 2 + p 1 2 -1 p 1 p 0 2 + p 1 2 -1 p 0 p 1 -p 0 p 1 , p 0 p 0 2 + p 1 2 -1 p 0 p 0 2 + p 1 2 -1 p 0 p 1 -p 0 p 1 , 1     .

Hyperbolic bisector

As noticed in the proof of Lemma 16, the hyperbolic bisector of p and p is suported by the Euclidean circle S such that φ(S) = line(φ(p), φ(p

)) ∩ π ∞ . φ(S) = p 0 2 + p 1 2 2(p 0 2 + p 1 2 )(p 0 2 + p 1 2 ) (p 0 , p 1 , p 0 2 + p 1 2 ) + p 0 2 + p 1 2 2(p 0 2 + p 1 2 )(p 0 2 + p 1 2 ) (p 0 , p 1 , p 0 2 + p 1 2 )

Hyperbolic circumcenter

As observed in the proof of Proposition 18 the hyperbolic center x of the triangle pp p can be computed as

ψ Π (φ(p) * ∩ φ(p ) * ∩ φ(p ) * ). Let 1 δ (α, β, γ) be φ(p) * ∩ φ(p ) * ∩ φ(p ) * ), we have already computed in Equation 5 δ = 2p0 2p1 -1 2p 0 2p 1 -1 2p 0 2p 1 -1 , α = p0 2 + p1 2 2p1 -1 p 0 2 + p 1 2 2p 1 -1 p 0 2 + p 1 2 2p 1 -1 , β = 2p0 p0 2 + p1 2 -1 2p 0 p 0 2 + p 1 2 -1 2p 0 p 0 2 + p 1 2 -1 , γ = 2p0 2p1 p0 2 + p1 2 2p 0 2p 1 p 0 2 + p 1 2 2p 0 2p 1 p 0 2 + p 1 2 then let λ be the smallest solution of λ 2 + γ-δ √ α 2 +β 2 λ -1 = 0 that is λ = - γ -δ 4α 2 + 4β 2 - δ 2 -2δγ + γ 2 4α 2 + 4β 2 -1 = δ -γ -(γ -δ) 2 -4α 2 -4β 2 4α 2 + 4β 2
and the hyperbolic center is

x = 1 + λ 2 δ(1 + γ) (α, β).
Evaluation CGAL provides us with number type Sqrt extension for exact computations on algebraic numbers of degree 2 [START_REF] Hemmer | CGAL User and Reference Manual[END_REF], allowing us to construct in an exact way the hyperbolic center of a circle as well as hyperbolic bisector between two input points. These constructions are rounded only when displaying the hyperbolic Voronoi diagram.

Experimental results

Experiments are run on point sets P that are uniformly distributed (up to rounding errors when generating them), both according to the Euclidean metric and according to the hyperbolic metric in open balls in B (See Figure 13 for the hyperbolic metric). Each open ball is centered at the origin of B with Euclidean radius 1-. The point sets in the open balls are denoted by small sph for = 0.001 and by big sph for = 0.0000001. The tables for the results are organized as follows. The first column gives the number of input points. The second and the third column give running times in seconds for the computations of DT E (P) and DT H (P), respectively. The fourth column shows the overhead factor of the computation of DT H (P) compared to DT E (P). Experiments are conducted on a MacBookPro 2.6 GHz running CGAL 4.0 in release mode, using GCC 4.2. Running times are averaged on 10 trials. For both DT E (P) and DT H (P), we use CGAL::Exact predicates inexact constructions kernel, which provides filtered exact geometric predicates.

Static variant

We insert all points in P at once, and we measure the running times of the computations of DT E (P) with the CGAL implementation [START_REF] Yvinec | 2D triangulations[END_REF] and of DT H (P) with our implementation. This allows us to measure the overhead due the static extraction of the Delaunay complex in H 2 .

Results on points uniformly distributed according to the hyperbolic metric: In this distribution in big sph with 10 7 points, the ratio (# edges in DT E (P)) / (# edges in DT H (P)) is 1.0005. We observe that the overhead of the extraction decreases with the size of the input point set, and becomes very small for large point sets (1% or 2%). The reason is that the graph of faces that is traversed by the extraction scheme grows slower than the whole graph of the Euclidean Inria triangulation. Faces examined during the extraction are in some sense "close" to the convex hull of the point set.

small
As could be expected, the overhead for small sph is better than the overhead obtained on big sph for larger point sets, since a small fraction of Euclidean Delaunay circles intersect H ∞ .

Dynamic variant

We insert points one by one and update DT H (P) at each insertion. As can be seen in the following table, the overhead is better with the static scheme, which can be explained by the additional is hyperbolic tests done on triangles that appear during the construction but are not present in the final triangulation.

Results on points uniformly distributed according to the hyperbolic metric: The overhead is already very small when points are uniformly distributed according to the Euclidean metric, since most points are internal to the triangulation, and a small fraction of insertions actually require performing is hyperbolic tests (using Lemma 8). This is not the case for the distribution according to the hyperbolic metric in big sph, since most points actually lie close to the Euclidean convex hull.

small

  the Poincaré ball model, the d-dimensional hyperbolic space H d is represented as the open unit ball B = {x ∈ E d : x < 1}. The points on the boundary of B are the points at infinity. The set of such points is H ∞ = {x ∈ E d : x = 1}. Hyperbolic lines, or geodesics, are represented either as arcs of Euclidean circles orthogonal to H ∞ or as diameters of B. See Figure 2 for an illustration in 2D. The Poincaré ball model of H d can be embedded into E d+1 by identifying it with the open unit d-disk χ = 0, x < 1.

Figure 4 :

 4 Figure 4: Hyperbolic Delaunay complex.

Lemma 3

 3 Let P be partitioned into a set R of red points and a set B of blue points. Then there exist two points p r ∈ R and p b ∈ B, and an empty sphere S pq passing through p r and p b that is contained in B.

Figure 5 :Proposition 5

 55 Figure 5: Proof of Prop 6.

Figure 6 :

 6 Figure 6: Proof of Lemma 8.

Figure 8 :

 8 Figure 8: Condition for a simplex σ of DT E (P) to be a simplex of DT H (P).

Proposition 10 d 2 + 2 .

 1022 The dynamic variant has optimal randomized worst-case time complexity O n n log n and space complexity O n d

1 .

 1 Static scheme. The depth-first search of the graph G defined in Proposition 6 and marking its nodes are implemented as a private member function of the class Delaunay hyperbolic triang- ulation 2<Gt, Tds>. Since the static scheme first computes the whole Euclidean triangulation before performing the extraction, then the member function insert(InputIterator first, InputIterator last) of Delaunay hyperbolic triangulation 2<Gt, Tds> inherited from the base class Delaunay triangulation 2<Gt, Tds> (denoted by Base below) is overridden as follows in the pseudo-code listing: void insert(InputIterator first, InputIterator last) { Base::insert(InputIterator first, InputIterator last); mark hyperbolic(); // depth-first search Dynamic variant. At each insertion of a new point p, we first insert the point in the underlying Euclidean triangulation and then examine the star of p. The member function insert(const Point& p, Face handle start = Face handle() ) is overridden as shown in the following pseudocode (Σ is defined in Section 4.1 page 12): Vertex handle insert(const Point& p, Face handle start = Face handle() ) { Vertex handle v = Base::insert(p, start);

Figure 12 :

 12 Figure 12: Determining the edge h(pp p

Figure 13 :

 13 Figure 13: 100 uniformly distributed points according to the hyperbolic metric in an open ball in H 2 .

  Let pp p be a triangle of DT E (P) \ DT H (P) with circumscribing circle S pp p . The vertices p, p , p are given in counterclockwise order. Corollary 7 and Proposition 2 applied to DT E ({p, p , p }), yield that h(pp p ) is the unique edge in DT E ({p, p , p }) \ DT H ({p, p , p }).Let c E be the Euclidean center of S pp p . We define the vectors

  Let us additionaly mention that we observed that, in this distribution, the ratio (# edges in DT E (P)) / (# edges in DT H (P)) quickly decreases with the number of vertices. For instance, on big sph with 10 7 points, this ratio is 1.006.Results on points uniformly distributed according to the Euclidean metric:

	10 4	sph E 2 0.00800 0.0108 1.34 H 2 factor	big sph E 2 10 4 0.00745 0.0181 2.43 factor H 2
	10 5	0.0867	0.0982 1.13	10 5	0.0833	0.177	2.12
	10 6	0.890	0.935	1.05	10 6	0.898	1.068	1.18
	10 7	9.61	9.80	1.02	10 7	10.1	10.4	1.02
	small sph E 2 10 4 0.00782 0.00854 1.09 H 2 factor	big sph E 2 10 4 0.00776 0.00837 1.07 factor H 2
	10 5	0.0852	0.0875	1.02	10 5	0.0853	0.0870	1.02
	10 6	0.882	0.903	1.02	10 6	0.892	0.912	1.02
	10 7	9.50	9.60	1.01	10 7	9.58	9.74	1.01

see also http://www.ics.uci.edu/ ~eppstein/pubs/geom-hyperbolic.html RR n°8146

RR n°8146

Degeneracies in the Counterclockwise test cannot occur since they would correspond either to equality of two points of {p, p , p } ⊂ B or equality of one of these points to p S pp p ∈ B.

Inria

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Conclusion and Open Problems

We proposed a algorithm to compute hyperbolic Voronoi diagrams and Delaunay complexes, as well as a complete and efficient implementation, to be submitted to CGAL.

We are pursuing research on periodic hyperbolic Delaunay triangulations, in the flavor of what we had proposed for the Euclidean case [11,12]. This is motivated by applications in various fields such as geometry processing [31] and neuro mathematics [22].