Jean-Daniel Boissonnat

Camille Wormser

Mariette Yvinec

Curved Voronoi Diagrams

Voronoi diagrams are fundamental data structures that have been extensively studied in Computational Geometry. A Voronoi diagram can be defined as the minimization diagram of a finite set of continuous functions. Usually, each of those functions is interpreted as the distance function to an object. The associated Voronoi diagram subdivides the embedding space into regions, each region consisting of the points that are closer to a given object than to the others. We may define many variants of Voronoi diagrams depending on the class of objects, the distance functions and the embedding space. Affine diagrams, i.e. diagrams whose cells are convex polytopes, are well understood. Their properties can be deduced from the properties of polytopes and they can be constructed efficiently. The situation is very different for Voronoi diagrams with curved regions. Curved Voronoi diagrams arise in various contexts where the objects are not punctual or the distance is not the Euclidean distance. We survey the main results on curved Voronoi diagrams. We describe in some detail two general mechanisms to obtain effective algorithms for some classes of curved Voronoi diagrams. The first one consists in linearizing the diagram and applies, in particular, to diagrams whose bisectors are algebraic hypersurfaces. The second one is a randomized incremental paradigm that can construct affine and several planar non-affine diagrams.

We finally introduce the concept of Medial Axis which generalizes the concept of Voronoi diagram to infinite sets. Interestingly, it is possible to efficiently construct a certified approximation of the medial axis of a bounded set from the Voronoi diagram of a sample of points on the boundary of the set.

Introduction

Voronoi diagrams are fundamental data structures that have been extensively studied in Computational Geometry. Given n objects, the associated Voronoi diagram subdivides R d into regions, each region consisting of the points that are closer to a given object than to any other object. We may define many variants of Voronoi diagrams depending on the class of objects, the distance function and the embedding space. Although Voronoi diagrams are most often defined in a metric setting, they can be defined in a more abstract way. In Sect. 1.2, we define them as minimization diagrams of any finite set of continuous functions without referring to a set of objects.

Given a finite set of objects and associated distance functions, we call bisector the locus of the points that are at equal distance from two objects. Voronoi diagrams can be classified according to the nature of the bisectors of the pairs of objects, called the bisectors of the diagram for short. An important class of Voronoi diagrams is the class of affine diagrams, whose bisectors are hyperplanes. Euclidean Voronoi diagrams of finite point sets are affine diagrams. Other examples of affine diagrams are the so-called power (or Laguerre) diagrams, where the objects are no longer points but hyperspheres and the Euclidean distance is replaced by the power of a point to a hypersphere. In Sect. 1.3, we recall well-known facts about affine diagrams. In particular, we characterize affine diagrams and establish a connection between affine diagrams and polytopes. As a consequence, we obtain tight combinatorial bounds and efficient algorithms. We also obtain a dual structure that is a triangulation under a general position assumption.

Non-affine diagrams are by far less well understood. Non-affine diagrams are obtained if one changes the distance function: additively and multiplicatively weighted distances are typical examples. Such diagrams allow to model growing processes and have important applications in biology, ecology, chemistry and other fields (see Sect. 1.9). Euclidean Voronoi diagrams of nonpunctual objects are also non-affine diagrams. They are of particular interest in robotics, CAD and molecular biology. Even for the simplest diagrams, e.g. Euclidean Voronoi diagrams of lines, triangles or spheres in 3-space, obtaining tight combinatorial bounds, efficient algorithms and effective implementations are difficult research questions.

A first class of non-affine diagrams to be discussed in Sect. 1.4 is the case of diagrams whose bisectors are algebraic hypersurfaces. We first consider the case of Möbius diagrams whose bisectors are hyperspheres and the case of anisotropic diagrams whose bisectors are quadratic hypersurfaces (see Sect. 1.4.2). The related case of Apollonius (or Johnson-Mehl) diagrams is also described in Sect. 1. [START_REF] Attali | Stability and computation of medial axes: a state of the art report[END_REF].

The key to obtaining effective algorithms for computing those non-affine diagrams is a linearization procedure that reduces the construction of a nonaffine diagram to intersecting an affine diagram with a manifold in some higher dimensional space. This mechanism is studied in full generality in Sect. 1.5.

In this section, we introduce abstract diagrams, which are diagrams defined in terms of their bisectors. By imposing suitable conditions on these bisectors, any abstract diagram can be built as the minimization diagram of some distance functions, thus showing that the class of abstract diagrams is the same as the class of Voronoi diagrams. Furthermore, the linearization technique introduced in Sect. 1.5 allows to prove that if the bisectors of a diagram belong to a certain class of bisectors, the distance functions defining the diagram can be chosen among a precise class of functions. For instance, affine diagrams are identified with power diagrams, spherical diagrams are identified with Möbius diagrams, and quadratic diagrams with anisotropic diagrams.

In Sect. 1.6, we introduce the incremental paradigm for constructing various diagrams. Under some topological conditions to be satisfied by the diagram, the incremental construction is efficient. The algorithm can be further improved by using a randomized data structure called the Voronoi hierarchy that allows fast localization of new objects. We then obtain fast randomized incremental algorithms for affine diagrams in any dimension and several nonaffine diagrams in the plane. Going beyond those simple cases is difficult. As mentioned above, tight combinatorial bounds and efficient algorithms are lacking even for simple cases. Moreover, the numerical issues are delicate and robust implementations are still far ahead of the state of the art. This motivates the quest for approximate solutions.

In Sect. 1.7, we introduce the concept of Medial Axis of a bounded set Ω, which can be seen as an extension of the notion of Voronoi diagram to infinite sets. Interestingly, it is possible to construct certified approximations of the medial axis of quite general sets efficiently. One approach to be described consists in sampling the boundary of Ω and then computing an appropriate subset of the Voronoi diagram of the sample that approximates the medial axis. Hence the problem of approximating the medial axis of Ω boils down to sampling the boundary of Ω, a problem that is closely related to mesh generation (see Chap. ??).

Sect. 1.8 is devoted to the main Cgal software packages for computing Voronoi diagrams. Sect. 1.9 discusses some applications of curved Voronoi diagrams.

This chapter focuses on curved Voronoi diagrams defined in R d and aims at providing useful background and effective algorithms. Additional material can be found in surveys on Voronoi diagrams [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF][START_REF] Aurenhammer | Voronoi diagrams[END_REF] and in text books on Computational Geometry [17,[START_REF] Boissonnat | Algorithmic Geometry[END_REF]. This chapter does not consider Voronoi diagrams defined in more general spaces. Voronoi diagrams can be defined in hyperbolic geometry without much difficulty [9,[START_REF] Boissonnat | Algorithmic Geometry[END_REF]. In the Poincaré model of hyperbolic geometry, the bisectors are hyperspheres and hyperbolic diagrams of finite point sets are a special case of Möbius diagrams. Computing Voronoi diagrams on Riemannian manifolds is much more involved and very few is known about such diagrams and their construction [38].

Notation: We identify a point x ∈ R d and the vector of its coordinates. We note x • y the dot product of x and y, x 2 = x • x = x 2 the squared Euclidean norm of x, and xy the Euclidean distance between points x and y.

We call hypersurface a manifold of codimension 1. Examples to be used in this chapter are hyperplanes, hyperspheres and quadratic hypersurfaces.

Lower Envelopes and Minimization Diagrams

Let F = {f 1 , . . . , f n } be a set of d-variate continuous functions defined over R d . The lower envelope of F is defined as

F -= min 1≤i≤n f i .
From F and F -, we define a natural partition of R d called the minimization diagram of F. For a point x ∈ R d , we define the index set I(x) of x as the set of all indices i such that F -(x) = f i (x). An equivalence relation noted ≡ can then be defined between two points of R d if they have the same index set:

x ≡ y ⇔ I(x) = I(y). The index set of a face is defined as the largest subset of indices common to all the points of the face. Conversely, the face of index set I is the set of all points x such that I ⊂ I(x).

Observe that the faces of this diagram are not necessarily contractible nor even connected. In particular, a 0-dimensional face may consist of several distinct points.

Lower envelopes and minimization diagrams have been well studied. We recall an important result due to Sharir [45] which provides an almost optimal result when the f i are supposed to be multivariate polynomials of constant maximum degree.

Theorem 1 (Sharir). The number of faces of the minimization diagram of a set F of n multivariate polynomials of constant maximum degree η is O(n d+ε ) for any ε > 0, where the constant of proportionality depends on ε, d and η. The vertices, edges and 2-faces of the diagram can be computed in randomized expected time O(n d+ε ) for any ε > 0.

This general result is close to optimal in the worst-case (see Exercise 2). It has been improved in some special cases. For more information and other related results, one should consult the book by Sharir and Agarwal [46].

Voronoi diagrams, in their general setting, are just minimization diagrams of a finite set of continuous functions. This general definition encompasses the more traditional definition of Voronoi diagrams where the functions are defined as distance functions to a finite set of objects. Theorem 1 provides very general bounds on the complexity of Voronoi diagrams. However, this result calls for improvement. First, in some special cases, much better bounds can be obtained by other approaches to be discussed later in this chapter. In particular, we will see that the most popular Euclidean Voronoi diagram of points has a much smaller combinatorial complexity than the one given in the theorem.

A second issue is the algorithmic complexity. The algorithm mentioned in the theorem fails to provide a complete description of the diagram since only faces of dimensions up to 2 are computed.

Moreover, the implementation of such an algorithm remains a critical issue. As evidenced in Chap. ??, computing lower envelopes of algebraic functions is a formidable task, even in the simplest cases, e.g. quadratic bi-variate functions. We do not know of any implementation for higher degrees and dimensions.

The main goal of the following sections is to present effective algorithms for a variety of Voronoi diagrams for which some additional structure can be exhibited.

Exercise 1. Show that the combinatorial complexity of the lower envelope of n univariate functions whose graphs intersect pairwise in at most two points is O(n). Show that the envelope can be computed in optimal time Θ(n log n).

Exercise 2. Show that the convex hull of n ellipsoids of R d may have Ω(n d-1 ) faces. Since the non-bounded faces of the Euclidean Voronoi diagram of n objects are in 1-1 correspondence with the faces of their convex hull, we get a lower bound on the size of the Voronoi diagram of n ellipsoids of R d . (Hint: consider n ellipsoids inscribed in a (d -1)-sphere S and intersecting S along great n (d -2)-spheres σ 1 , . . . , σ n . The arrangement of the σ i has Θ(n d-1 ) faces.)

Affine Voronoi Diagrams

We first introduce Euclidean Voronoi diagrams of points and establish a correspondence between those diagrams and convex polyhedra in one dimension higher. Polarity allows to associate to a Voronoi diagram its dual cell complex, called a Delaunay triangulation.

Almost identical results can be obtained for power (or Laguerre) diagrams where points are replaced by hyperspheres and the Euclidean distance by the power of a point to a hypersphere. Power diagrams constitute a natural extension of Euclidean Voronoi diagrams and are still affine diagrams. In fact, we will see that any affine diagram is the power diagram of a finite set of hyperspheres.

Euclidean Voronoi Diagrams of Points

Let P = {p 1 , . . . , p n } be a set of points of R d . To each p i , we associate its Voronoi region V (p i )

V (p i ) = {x ∈ R d : x -p i ≤ x -p j , ∀j ≤ n}.
The region V (p i ) is the intersection of n -1 half-spaces. Each such half-space contains p i and is bounded by the bisector of p i and some other point of P. Since the bisectors are hyperplanes, V (p i ) is a convex polyhedron, possibly unbounded.

The Euclidean Voronoi diagram of P, noted Vor(P), is the cell complex whose cells are the Voronoi regions and their faces. Equivalently, the Euclidean Voronoi diagram of P can be defined as the minimization diagram of the distance functions δ i , . . . , δ n , where

δ i (x) = x -p i .
In other words, the Euclidean Voronoi diagram of P is the minimization diagram of a set of functions whose graphs are vertical 1 cones of revolution of 1 By vertical, we mean that the axis of revolution is perpendicular to R d . R d+1 . Since minimizing x-p i over i is the same as minimizing (x-p i ) 2 , the Euclidean Voronoi diagram of P can alternatively be defined as the mimization diagram of the smooth functions (xp i ) 2 whose graphs are translated copies of a vertical paraboloid of revolution of R d+1 . Observing further that, for any x, arg min i (x

-p i ) 2 = arg min i (-2p i • x + p 2
i ), we obtain that the Euclidean Voronoi diagram of P is the minimization diagram of a set of affine functions, namely the functions

d i (x) = -2p i • x + p 2 i
whose graphs are hyperplanes of R d+1 . Let us call h pi , i = 1, . . . , n, those hyperplanes and let h - pi denote the half-space lying below h pi . The minimization diagram of the d i is obtained by projecting the polyhedron

V(P) = h - p1 ∩ • • • ∩ h - pn .
vertically onto R d . See Fig. 1.3. We have therefore proved the following theorem:

Delaunay Triangulation

Two cell complexes V and D are said to be dual if there exists an involutive correspondence between the faces of V and the faces of D that reverses the inclusions, i.e. for any two faces f and g of V , their dual faces f * and g * satisfy: f ⊂ g ⇒ g * ⊂ f * . We introduce now a cell complex that is dual to the Voronoi diagram of a finite set of points P.

We assume for now that the set of points P is in general position, which means that no subset of d + 2 points of P lie on a same hypersphere. Let f be a face of dimension k of the Voronoi diagram of P. All points in the interior of f have the same subset P f of closest points in P. The face dual to f is the convex hull of P f . The Delaunay triangulation of P, noted Del(P), is the cell complex consisting of all the dual faces. Because points of P are assumed to be in general position, |P f | = dk + 1, all the faces of Del(P) are simplices and Del(P) is a simplicial complex. The fact that Del(P) is indeed a triangulation, i.e. a simplicial complex embedded in R d and covering the convex hull of P, will be proved now using a duality between points and hyperplanes in the so-called space of spheres.

Polarity

Let σ be the hypersphere of R d of equation

σ(x) = (x -c) 2 -r 2 = x 2 -2c • x + s = 0,
where c is the center of σ, r its radius and

s = σ(0) = c 2 -r 2 .
We define the following bijective mapping φ :

σ ∈ R d -→ φ(σ) = (c, -s) ∈ R d+1
that maps a hypersphere of R d to a point of R d+1 . We thus consider R d+1 as the images by φ of the hyperspheres of R d and call R d+1 the space of spheres. We note φ(p) = (p, -p 2 ) the image by φ of a point, considered as a hypersphere of radius 0. Observe that φ(p) is a point of the paraboloid Q of R d+1 of equation x 2 + x d+1 = 0. The points of R d+1 that lie above Q are images of imaginary hyperspheres whose squared radii are negative. The points below Q are images of real hyperspheres. We now introduced a mapping between points and hyperplanes of the space of spheres, known as polarity. Polarity associates to the point φ(σ) its polar hyperplane h σ which is the hyperplane of R d+1 of equation 2c•x+x d+1 -s = 0. Observe that the intersection of h σ with Q projects vertically onto σ, and that h σ is the affine hull of the image by φ of the points of σ. If p is a point of R d , the polar hyperplane h p of φ(p) is the hyperplane tangent to Q at φ(p).

We deduce the remarkable following property: x ∈ σ if and only if φ(x) = (x, -x 2 ) ∈ h σ and σ encloses x if and only if φ(x) ∈ h + σ , where h + σ (resp. h - σ ) denotes the closed half-space above (resp. below) h σ . Indeed

σ(x) = 0 ⇐⇒ x 2 -2c • x + s = 0 ⇐⇒ φ(x) ∈ h σ σ(x) < 0 ⇐⇒ x 2 -2c • x + s < 0 ⇐⇒ φ(x) ∈ int h + σ , where int h +
σ denotes the open half-space above h σ . Polarity is an involution that preserves incidences and reverses inclusions. Indeed, if σ and σ are two hyperspheres, we have

φ(σ) ∈ h σ ⇐⇒ 2c • c -s -s = 0 ⇐⇒ φ(σ ) ∈ h σ φ(σ) ∈ h + σ ⇐⇒ 2c • c -s -s > 0 ⇐⇒ φ(σ ) ∈ int h + σ .
Consider now a set P = {p 1 , . . . , p n } of n points and let V(P) denote, as in Sect. 1.3.1, the convex polyhedron defined as the intersection of the n halfspaces below the n polar hyperplanes h p1 , . . . , h pn . Let f be a face of V(P) and assume that f is contained in k + 1 hyperplanes among the h pi . Without loss of generality, we denote those hyperplanes h p1 , . . . , h p k+1 . Let σ denote a hypersphere of R d such that φ(σ) belongs to the relative interior of f . From the above discussion, we have ∀i,

1 ≤ i ≤ k + 1, φ(σ) ∈ h pi ⇐⇒ φ(p i ) ∈ h σ (1.1) ∀i, k + 1 < i ≤ n, φ(σ) ∈ int h - pi ⇐⇒ φ(p i ) ∈ int h - σ (1.2)
Given a convex polyhedron D, we say that a hyperplane h supports D if D ∩ h is non-empty and D is included in one of the two halfspaces, h + or h -, bounded by h. If h is a supporting hyperplane of D, g = D ∩ h is a face of D. If D ⊂ h -, g is called an upper face of D. The collection of all upper faces of D constitutes the upper hull of D, which we denote by ∂ + D.

Let D(P) = conv(φ(P)) be the convex hull of the set φ(P) and consider again the face f of V(P) defined above. Write P f = {p 1 , . . . , p k+1 }. We deduce from (1.1) and (1.2) that, for any φ(σ) in the relative interior of f :

1. The hyperplane h σ is a supporting hyperplane of D(P). 2. h σ supports D(P) along the face f * = h σ ∩ D(P) = conv(φ((P f )).

D(P) ⊂ h -

σ and f * is a face of ∂ + D(P). To each face f of ∂V(P), we associate the face f * of ∂ + D obtained as described above. This correspondence between the faces of ∂V(P) and the faces of ∂ + D(P) is bijective, preserves incidences and reverses inclusions, hence it is a duality.

The upper hull ∂ + D(P) projects vertically onto a cell complex of R d whose vertices are the points of P. Because the projection is 1-1, this projected cell complex is properly embedded in R d and, since the projection preserves convexity, it covers the convex hull of P. Under the general position assumption, the convex polyhedron D(P) is simplicial and the projected complex is a triangulation of P. The duality between the faces of ∂V(P) and the faces of ∂ + D(P) implies that the projection of ∂ + D(P) is the Delaunay triangulation Del(P) of P introduced at the beginning of this section. This concludes the proof that, under the general position assumption, the Delaunay triangulation Del(P) is a triangulation of P. We have the following diagram:

∂V(P) = ∂ h - p1 ∩ • • • ∩ h - pn ←→ ∂ + D(P) = ∂ + (conv(φ(P)))
Voronoi Diagram Vor(P) ←→ Delaunay Triangulation Del(P)

It follows from the above correspondence that the combinatorial complexity of the Delaunay triangulation of n points is the same as the combinatorial complexity of its dual Voronoi diagram. Moreover, the Delaunay triangulation of n points of R d can be deduced from the dual Voronoi diagram or vice versa in time proportional to its size. We also deduce from what precedes that computing the Delaunay triangulation of n points of R d reduces to constructing the convex hull of n points of R d+1 . The following theorem is then a direct consequence of known results on convex hulls [15]. . Both structures can be computed in optimal time Θ n log n + n d+1 2

.

The bounds in this theorem are tight. In particular, the Voronoi diagram of n points of R 3 may be quadratic (see Exercise 4). These bounds are worst-case bounds. Under some assumptions on the point distribution, better bounds can be obtained. In the discussion above, we have assumed that the points of P were in general position. If this is not the case, some faces of D(P) are not simplices, and the complex ∂ + D(P) projects vertically onto a cell complex, dual to the Voronoi diagram and called the Delaunay complex. The faces of the Delaunay complex are convex and any triangulation obtained by triangulating those faces is called a Delaunay triangulation. Since there are several ways of triangulating the faces of the Delaunay complex, the Delaunay triangulation of P is no longer unique. Exercise 4. Show that if we take points on two non coplanar lines of R 3 , say n 1 + 1 on one of the lines and n 2 + 1 on the other, their Voronoi diagram has n 1 n 2 vertices.

Exercise 5. Let S be a hypersphere of R d passing through d + 1 points p 0 , . . . , p d . Show that a point p d+1 of R d lies on S, in the interior of the ball B S bounded by S or outside B S , depending whether the determinant of the (d + 2) × (d + 2) matrix in sphere(p 0 , . . . ,

p d+1 ) = 1 • • • 1 p 0 • • • p d+1 p 2 0 • • • p 2 d+1
is 0, negative or positive. This predicate is the only numerical operation that is required to check if a triangulation is a Delaunay triangulation.

Exercise 6. What are the preimages by φ of the points of R d+1 that lie on a line? (Distinguish the cases where the line intersects Q in 0, 1 or 2 points.)

Exercise 7. Project vertically the faces of the lower hull ∂ -(D(P). Show that we obtain a triangulation of the vertices of conv(P) such that each ball circumscribing a simplex contains all the points of P. Define a dual and make a link with Exercise 3.

Exercise 8 (Empty sphere property). Let s be any k-simplex with vertices in P that can be circumscribed by a a hypersphere that does not enclose any point of P. Show that s is a face of a Delaunay triangulation of P. Moreover, let P be a set of points and T a triangulation of P with the property that any hypersphere circumscribing a d-simplex of T does not enclose any point of P. Show that T is a Delaunay triangulation of P.

Power Diagrams

A construction similar to what we did for the Euclidean Voronoi diagrams of points and their dual Delaunay triangulations can be done for the so-called power or Laguerre diagrams. Here we take as our finite set of objects a set of hyperspheres (instead of points) and consider as distance function of a point x to a hypersphere σ the power of x to σ. As we will see, the class of power diagrams is identical to the class of affine diagrams, i.e. the diagrams whose bisectors are hyperplanes.

Definition of Power Diagrams

We call power of a point x to a hypersphere σ of center c and radius r the real number σ(x) = (xc) 2r 2 .

Let S = {σ 1 , . . . , σ n } be a set of hyperspheres of R d . We denote by c i the center of σ i , r i its radius, σ i (x) = (xc i ) 2r 2 i the power function to σ i , and s i = c 2 ir 2 i the power of the origin. To each σ i , we associate the region L(σ i ) consisting of the points of R d whose power to σ i is not larger than their power to the other hyperspheres of S:

L(σ i ) = {x ∈ R d : σ i (x) ≤ σ j (x), 1 ≤ j ≤ n}.
The set of points that have equal power to two hyperspheres σ i and σ j is a hyperplane, noted π ij , called the radical hyperplane of σ i and σ j . Hyperplane π ij is orthogonal to the line joining the centers of σ i and σ j . We denote by π i ij the half-space bounded by π ij consisting of the points whose power to σ i is smaller than their power to σ j . The region L(σ i ) is the intersection of all halfspaces π i ij , j = i. If this intersection is not empty, it is a convex polyhedron, possibly not bounded. We call power regions the non empty regions L(σ i ).

We define the power diagram of S, noted Pow(S), as the cell complex whose cells are the power regions and their faces. When all hyperspheres have the same radius, their power diagram is identical to the Voronoi diagram of their centers. Equivalently, the power diagram of S can be defined as the minimization diagram of the functions σ i , . . . , σ n . Observing that for any x

arg min i σ i (x) = arg min i (-2c i • x + s i ),
we obtain that the power diagram of S is the minimization diagram of the set of affine functions

d i (x) = -2p i • x + s i
whose graphs are hyperplanes of R d+1 . Let us call h σi , i = 1, . . . , n, those hyperplanes and let h - σi denote the half-space lying below h σi . The minimization diagram of the δ i is obtained by projecting vertically the convex polyhedron

L(S) = h - p1 ∩ • • • ∩ h - pn .
Theorem 4. The faces of the power diagram Pow(S) of S are the vertical projections of the faces of the convex polyhedron L(S).

Power diagrams are very similar to Voronoi diagrams: the only difference is that the hyperplanes supporting the faces of L(S) are not necessarily tangent to the paraboloid Q and that some hyperplane may not contribute a face. In other words, some hypersphere σ i may have an empty power region (see the small circle in the upper left corner of Fig. 1.6).

By proceeding as in Sect. 1.3.2, we can define a convex polyhedron R(S) whose upper hull ∂ + R(S) is dual to ∂L(S). The vertical projection of the faces of ∂ + R(S) constitute the faces of a cell complex which, in general, is a simplicial complex. We call such a complex the regular triangulation of S and denote it by Reg(S). We have the following diagram :

∂L(S) = ∂ h - σ1 ∩ • • • ∩ h - σn ←→ ∂ + R(S) = ∂ + conv(φ(S)) Power diagram Pow(S) ←→ Regular triangulation Reg(S)
We deduce the following theorem that states that computing the power diagram of n hyperspheres of R d (or equivalently its dual regular triangulation) has the same asymptotic complexity as computing the Euclidean Voronoi diagram or the Delaunay triangulation of n points of R d . . Both structures can be computed in optimal time Θ n log n + n d+1 2

.

Affine Voronoi Diagrams

Euclidean Voronoi diagrams of points and power diagrams of hyperspheres are two examples of minimization diagrams whose bisectors are hyperplanes.

It is interesting to classify Voronoi diagrams with respect to their bisectors. A first important class of Voronoi diagrams is the class of affine diagrams which consists of all Voronoi diagrams whose bisectors are hyperplanes.

In Sect. 1.5, we will prove that any affine Voronoi diagram of R d is identical to the power diagram of some set of hyperspheres of R d (Theorem 13), therefore showing that the class of affine Voronoi diagrams is identical to the class of power diagrams.

Exercise 9. Show that the intersection of a power diagram with an affine subspace is still a power diagram and compute the corresponding spheres.

Exercise 10. Show that any power diagram of R d is the intersection of a Voronoi diagram of R d+1 by a hyperplane.

Exercise 11. Show that the only numerical operation that is required to check if a triangulation is the regular triangulation of a set of hyperspheres σ i is the evaluation of the sign of the determinant of the (d + 2) × (d + 2) matrix

power test(σ 0 , . . . , σ d+1 ) = 1 • • • 1 c 0 • • • c d+1 c 2 0 -r 2 0 • • • c 2 d+1 -r 2 d+1
where c i and r i are respectively the center and the radius of σ i .

Voronoi Diagrams with Algebraic Bisectors

In this section, we introduce a first class of non-affine diagrams, namely the class of diagrams whose bisectors are algebraic hypersurfaces. We first consider the case of Möbius diagrams whose bisectors are hyperspheres and the case of anisotropic diagrams whose bisectors are quadratic hypersurfaces. These diagrams can be computed through linearization, a technique to be described in full generality in Sect. 1.5. Apollonius (or Johnson-Mehl) diagrams, although semi-algebraic and not algebraic, are also described in this section since they are closely related to Möbius diagrams and can also be linearized.

Möbius Diagrams

In this section, we introduce a class of non-affine Voronoi diagrams, the socalled Möbius diagrams, introduced by Boissonnat and Karavelas [START_REF] Boissonnat | On the combinatorial complexity of Euclidean Voronoi cells and convex hulls of d-dimensional spheres[END_REF].

The class of Möbius diagrams includes affine diagrams. In fact, as we will see, the class of Möbius diagrams is identical to the class of diagrams whose bisectors are hyperspheres (or hyperplanes).

Definition of Möbius Diagrams

Let ω= {ω 1 , . . . , ω n } be a set of so-called Möbius sites of R d , where ω i is a triple (p i , λ i , µ i ) formed of a point p i of R d , and two real numbers λ i and µ i .

For a point x ∈ R d , the distance δ i (x) from x to the Möbius site ω i is defined as

δ i (x) = λ i (x -p i ) 2 -µ i .
Observe that the graph of δ i is a paraboloid of revolution whose axis is vertical. The Möbius region of the Möbius site ω i , i = 1, . . . , n, is

M (ω i ) = {x ∈ R d : δ i (x) ≤ δ j (x), 1 ≤ j ≤ n}.
Observe that a Möbius region may be non-contractible and even disconnected. Möbius diagrams are generalizations of Euclidean Voronoi and power diagrams. In particular, if all λ i are equal to some positive λ, the Möbius diagram coincides with the power diagram of a set of spheres {σ i , i = 1, . . . , n}, where σ i is the sphere centered at p i of squared radius µ i /λ. If all µ i are equal and all λ i are positive, then the Möbius diagram coincides with the so-called multiplicatively weighted Voronoi diagram of the weighted points (p i , √ λ i ). The following lemma states that the bisector of two Möbius sites is a hypersphere (possibly degenerated in a point or in a hyperplane). Its proof is straightforward.

Lemma 1. Let ω i = {p i , λ i , µ i } and ω j = {p j , λ j , µ j }, ω i = ω j be two Möbius sites. The bisector σ ij of ω i and ω j is the empty set, a single point, a hypersphere or a hyperplane.

Möbius Diagrams and Power Diagrams

We now present an equivalence between Möbius diagrams in R d and power diagrams in R d+1 . This result is a direct generalization of a similar result for multiplicatively weighted diagrams [6]. Given a cell complex C covering a subspace X, we call restriction of C to X the subdivision of X whose faces are the intersections of the faces of C with X. The restriction of C to X is denoted by C X . Note that the restriction C X is not, in general, a cell complex and that its faces may be non-contractible and even non-connected.

We associate to ω= {ω 1 , . . . , ω n } the set of hyperspheres

Σ = {Σ 1 , . . . , Σ n } of R d+1 of equations Σ i (X) = X 2 -2C i • X + s i = 0,
where

C i = (λ i p i , -λi 2 ) and s i = λ i p 2 i -µ i . We denote by Q the paraboloid of R d+1 of equation x d+1 -x 2 = 0 .
Theorem 6 (Linearization). The Möbius diagram Möb(ω) of ω is obtained by projecting vertically the faces of the restriction

Pow Q (Σ) of the power diagram of Σ to Q. Proof. If x ∈ R d is closer to ω i than to ω j with respect to δ M , we have for all j = 1, . . . , n, λ i (x -p i ) 2 -µ i ≤ λ j (x -p j ) 2 -µ j ⇐⇒ λ i x 2 -2λ i p i • x + λ i p 2 i -µ i ≤ λ j x 2 -2λ j p j • x + λ j p 2 j -µ j ⇐⇒ (x 2 + λi 2 ) 2 + (x -λ i p i ) 2 - λ 2 i 4 -λ 2 i p 2 i + λ i p 2 i -µ i ≤ (x 2 + λj 2 ) 2 + (x -λ j p j ) 2 - λ 2 j 4 -λ 2 j p 2 j + λ j p 2 j -µ j ⇐⇒ (X -C i ) 2 -r 2 i ≤ (X -C j ) 2 -r 2 j ⇐⇒ Σ i (X) ≤ Σ j (X)
where

X = (x, x 2 ) ∈ Q ⊂ R d+1 , C i = (λ i p i , -λi 2 ) ∈ R d+1 and r 2 i = λ 2 i p 2 i + λ 2 i 4 -λ i p 2 i + µ i .
The above inequality shows that x is closer to ω i than to ω j if and only if X belongs to the power region of Σ i in the power diagram of the hyperspheres Σ j , j = 1, . . . , n. As X belongs to Q and projects vertically onto x, we have proved the result.

Corollary 1. Let Σ be a finite set of hyperspheres of R d+1 , Pow(Σ) its power diagram and Pow Q (Σ) the restriction of Pow(Σ) to Q. The vertical projection of Pow Q (Σ) is the Möbius diagram Möb(ω) of a set of Möbius sites of R d .
Easy computations give ω.

Combinatorial and Algorithmic Properties

It follows from Theorem 6 that the combinatorial complexity of the Möbius

diagram of n Möbius sites in R d is O(n d 2 +1
). This bound is tight since Aurenhammer [6] has shown that it is tight for multiplicatively weighted Voronoi diagrams.

We easily deduce from the proof of the Linearization Theorem 6 an algorithm for constructing Möbius diagrams. First, we compute the power diagram of the hyperspheres Σ i of R d+1 , intersect each of the faces of this diagram with the paraboloid Q and then project the result on R d .

Theorem 7. Let ω be a set of n Möbius sites in R d , d ≥ 2. The Möbius dia- gram Möb(ω) of ω can be constructed in worst-case optimal time Θ(n log n + n d 2 +1 ).
Another consequence of the linearization theorem is the fact that any Möbius diagram can be represented as a simplicial complex

T Q embedded in R d+1 . T Q is a sub-complex of the regular triangulation T dual to the power diagram Pow(Σ i ) of the hyperspheres Σ i . Since T is embedded in R d+1 , T Q is a simplicial complex of R d+1 .
More precisely, T Q consists of the faces of T that are dual to the faces of Pow Q (Σ), i.e. the faces of the power diagram that intersect Q. We will call

T Q the dual of Pow Q (Σ). Observe that since, in general, no vertex of Pow(Σ) lies on Q, T Q is a d-dimensional simplicial complex (embedded in R d+1 ).
Moreover, if the faces of Pow(Σ) intersect Q transversally and along topological balls, then, by a result of Edelsbrunner and Shah [21], T Q is homeomorphic to Q and therefore to R d . It should be noted that this result states that the simplicial complex T Q has the topology of R d . This result, however, is mainly combinatorial, and does not imply that the embedding of T Q into R d+1 as a sub-complex of the regular triangulation T may be projected in a 1-1 manner onto R d .

Spherical Voronoi Diagrams

Lemma 1 states that the bisectors of two Möbius sites is a hypersphere (possibly degenerated in a hyperplane). More generally, let us consider the Voronoi diagrams such that, for any two objects o i and o j of O, the bisector

σ ij = {x ∈ R d , δ i (x) = δ j (x)} is a hypersphere. Such a diagram is called a spherical Voronoi diagram.
In Sect. 1.5, we will prove that any spherical Voronoi diagram of R d is a Möbius diagram (Theorem 15).

Möbius transformations are the transformations that preserve hyperspheres. An example of a Möbius transformation is the inversion with respect to a hypersphere. If the hypersphere is centered at c and has radius r, the inversion associates to a point x ∈ R d its image

x = c + r(x -c) (x -c) 2 .
Moreover, it is known that any Möbius transformation is the composition of up to four inversions [16]. An immediate consequence of Theorem 15 is that the set of Möbius diagrams in R d is stable under Möbius transformations, hence their name.

Möbius Diagrams on Spheres

Anisotropic Diagrams

The definition of anisotropic Voronoi diagrams presented in this section is a slight extension of the definition proposed by Labelle and Shewchuk [36]. The objects are points and the distance to a point is a quadratic form with an additive weight.

Anisotropic diagrams appear to be a natural generalization of Möbius diagrams and reduce to Möbius diagrams when the matrices are taken to be a scalar times the identity matrix. As will be shown, the class of anisotropic diagrams is identical to the class of diagrams whose bisectors are quadratic hypersurfaces.

Definition and linearization

Consider a finite set of anisotropic sites S = {s 1 , . . . , s n }. Each site

s i , i = 1, . . . , n, is a triple (p i , M i , π i ) formed by a point p i ∈ R d , a d × d symmetric positive definite matrix M i and a scalar weight π i . The distance δ i (x) of point x ∈ R d to site s i is defined by δ i (x) = (x -p i ) t M i (x -p i ) -π i .
The anisotropic Voronoi region of site s is then defined as

AV (s i ) = {x ∈ R d , δ i (x) ≤ δ j (x), ∀1 ≤ j ≤ n},
The anisotropic Voronoi diagram is the minimization diagram of the functions

δ i (x). Let D = d(d+3) 2 .
To each point x = (x 1 , . . . , x d ) ∈ R d , we associate the two points

φ(x) = (x r x i , 1 ≤ r ≤ s ≤ d) ∈ R d(d+1) 2 φ(x) = (x, φ(x)) ∈ R D ,
and we denote by Q the d-manifold of R D defined as

Q = φ(x), x ∈ R d .
To each site s i = (p i , M i , π i ) ∈ S, we associate:

1. the point mi ∈ R d(d+1) 2 defined as mu,u i = - 1 2 M u,u i , for 1 ≤ u ≤ d; mu,v i = -M u,v i , for 1 ≤ u < v ≤ d,
2. the point pi = (M i p i , mi ), 3. the sphere σ i of center pi and radius pi 2p t i M i p iπ i . Let Π be the projection ŷ = (y, ỹ) ∈ R D → y ∈ R d and let Σ be the set of spheres σ i , i = 1, . . . , n.

Theorem 8 (Linearization). The anisotropic diagram of S is the image by Π of the restriction of the power diagram Pow(Σ) to the d-manifold Q.

Proof. We have

δ i (x) = (x -p i ) t M (x -p i ) -π i = x t M i x -2p t i M i x + p t i M i p i -π i = -2p t i φ(x) + p t i M i p i -π i
This implies that δ i (x) < δ j (x) if and only if

( φ(x) -pi ) 2 -(p 2 i -p t i M i p i -π i ) < ( φ(x) -pj ) 2 -(p 2 j -p t j M j p j -π j )
Hence, x is closer to s i than to s j if and only if the power of φ(x) to σ i is smaller than its power to σ j . Equivalently, a point φ(x) ∈ Q belongs to the power cell of σ(s i ) if and only if its projection x = Π( φ(x)) belongs to the anisotropic Voronoi region AV (s i ).

We easily deduce the following theorem. ) where

D = d(d+3) 2 .
This result is to be compared to Theorem 1 which provides a better combinatorial bound. We let as an open question to fill the gap between those two bounds.

Quadratic Voronoi Diagrams

The bisectors of anisotropic diagrams, as defined in the previous section, are quadratic hypersurfaces. A minimization diagram whose bisectors are hyperquadrics is called a quadratic Voronoi diagram. In Sect. 1.5, we will prove that any quadratic Voronoi diagram is the anisotropic Voronoi diagram of a set of anisotropic sites (Theorem 16).

Apollonius Diagrams

In this section, we present diagrams that are closely related to Möbius diagrams: namely, the Euclidean Voronoi diagrams of hyperspheres, also called Apollonius or Johnson-Mehl diagrams. Contrary to Möbius and anisotropic diagrams, the bisectors of Apollonius diagrams are not algebraic hypersurfaces since the bisector between two hyperspheres is only one sheet of a hyperboloid. As a consequence, Apollonius diagrams cannot be linearized in the same way as Möbius and anisotropic diagrams. Nevertheless, another linearization scheme can be applied, leading to interesting combinatorial and algorithmic results.

Definition of Apollonius Diagrams

Let us consider a finite set of weighted points S = {σ 0 , σ 1 , . . . , σ n } where σ i = (p i , r i ), p i ∈ R d and r i ∈ R. We define the distance from x to σ i as

δ i (x) = x -p i -r i .
This distance is also called the additively weighted distance from x to the weighted point σ i . The minimization diagram of the distance functions δ i , i = 1, . . . , n, is called the additively weighted Voronoi diagram, or the Apollonius diagram of S. We denote it by Apo(S) (see Fig. 1.8).

The Apollonius region A(σ i ) of σ i is defined as

A(σ i ) = {x ∈ R d , δ i (x) ≤ δ j (x)}.
It is easy to see that A(σ i ) is either empty or star-shaped from p i . The boundary of A(σ i ) may have a complicated structure. In fact, as we will see, the boundary of A(σ i ) has the same combinatorial structure as a Möbius diagram in R d-1 . Since the diagram is not changed if we replace all r i by r i +r for any r ∈ R, we can assume, without loss of generality, that all r i are non negative. The weighted points are then hyperspheres and the distance to a weighted point is the signed Euclidean distance to the corresponding hypersphere, counted positively outside the hypersphere and negatively inside the hypersphere.

Observe that, in the plane, a vertex of an Apollonius diagram is the center of a circle tangent to three circles of S (assuming all r i non negative). Computing such a point is known as Apollonius' Tenth Problem, hence the name of the diagram.

Apollonius Diagrams and Power Diagrams

The graph of the distance function δ i (x) is the half-cone of revolution C i of equation The bisector of two hyperspheres of S is thus the projection of the intersection of two half-cones. This intersection is a quadratic hypersurface (in fact, a sheet of a two sheet hyperboloid) contained in a hyperplane. Indeed, we have

C i : x d+1 = x -p i -r i , x d+1 + r i ≥ 0
C 1 : (x d+1 + r 1 ) 2 = (x -p 1 ) 2 , x d+1 + r 1 > 0, C 2 : (x d+1 + r 2 ) 2 = (x -p 2 ) 2 , x d+1 + r 2 > 0.
The intersection of the two half-cones is contained in the hyperplane h 12 whose equation is obtained by subtracting the two sides of the above equations:

h 12 : -2(p 1 -p 2 ) • x -2(r 1 -r 2 )x d+1 + p 2 1 -r 2 1 -p 2 2 + r 2 2 = 0.
This shows that there exists a correspondence between the diagram Apo(S) and the power diagram of the hyperspheres Σ i in R d+1 (i = 1, . . . , n), where Σ i is centered at (p i , r i ) and has radius r i √ 2. More precisely, A(σ i ) is the projection of the intersection of the half-cone C i with the power region L(Σ i ). Indeed, x is in A(σ i ) if and only if the projection X i of x onto C i has a smaller x d+1 -coordinate than the projections of x onto the other half-cones C j , j = i. In other words, the coordinates (x, x d+1 ) of X i must obey

(x d+1 + r i ) 2 = (x -p i ) 2 (x d+1 + r j ) 2 ≤ (x -p j ) 2 for any j = i,
and by subtracting both sides, it follows that Σ i (X i ) ≤ Σ j (X i ) for all j.

Algorithm 1 Construction of Apollonius diagrams

Input: a set of hyperspheres S 1. Compute Σi, for i = 1, . . . , n; 2. Compute the power diagram of the Σi's; 3. For all i = 1, . . . , n, project vertically the intersection of the power region L(Σi)

with the half-cone Ci.

Output: the Apollonius diagram of S.

The Apollonius diagram of S can be computed using the following algorithm:

The power diagram of the Σ i can be computed in time O(n This result is optimal in odd dimensions, since the bounds above coincide with the corresponding bounds for the Voronoi diagram of points under the Euclidean distance. It is not optimal in dimension 2 (see Exercise 20). We also conjecture that it is not optimal in any even dimension.

Computing a Single Apollonius Region

We now establish a correspondence, due to Boissonnat and Karavelas [START_REF] Boissonnat | On the combinatorial complexity of Euclidean Voronoi cells and convex hulls of d-dimensional spheres[END_REF], between a single Apollonius region and a Möbius diagram on a hypersphere.

To give the intuition behind the result, we consider first the case where one of the hyperspheres, say σ 0 , is a hyperplane, i.e. a hypersphere of infinite radius. We take for σ 0 the hyperplane x d = 0, and assume that all the other hyperspheres lie the half-space x d > 0. The distance δ 0 (x) from a point x ∈ R d to σ 0 is defined as the Euclidean distance.

The points that are at equal distance from σ 0 and σ i , i > 0, belong to a paraboloid of revolution with vertical axis. Consider such a paraboloid as the graph of a (d -1)-variate function ϑ i defined over R d-1 . If follows from Sect. 1.4.1 that the minimization diagram of the ϑ i , i = 1, . . . , n, is a Möbius diagram (see Fig. 1.9).

Easy computations give the associated weighted points. Write p i = (p i , p i ),

p i ∈ R d-1
, p i ∈ R, i > 0 and let ω= {ω 1 , . . . , ω n } be the set of Möbius sites of R d where ω i = {p i , λ i , µ i }, and We let as an exercise to verify that the vertical projection of the boundary of the Apollonius region A(σ 0 ) of σ 0 onto σ 0 is the Möbius diagram of ω.

λ i = 1 r i + p i , µ i = r i -p i , i > 0.
We have assumed that one of the hyperspheres was a hyperplane. We now consider the case of hyperspheres of finite radii. The crucial observation is that the radial projection of A(σ 0 ) ∩ A(σ i ) ∩ A(σ j ) onto σ 0 , if not empty, is a hypersphere. ). Such a region can be computed in optimal time Θ(n log n+n d+1 2

).

Exercise 17. Show that the cell of hypersphere σ i in the Apollonius diagram of S is empty if and only if σ i is inside another hypersphere σ j .

Exercise 18. The predicates required to construct an Apollonius region are multivariate polynomials of degree at most 8 and 16 when d = 2 and 3 respectively. Detail these predicates [START_REF] Boissonnat | Convex hulls and Voronoi diagrams of additively weighted points[END_REF].

Exercise 19. Show that the convex hull of a finite number of hyperspheres can be deduced from the restriction of a power diagram to a unit hypersphere [START_REF] Boissonnat | Convex hulls and Voronoi diagrams of additively weighted points[END_REF].

Exercise 20. Prove that the combinatorial complexity of the Apollonius diagram of n circles in the plane has linear size.

Exercise 21 (Open problem). Give a tight bound on the combinatorial complexity of the Apollonius diagram of n hyperspheres of R d when d is even.

Linearization

In this section, we introduce abstract diagrams, which are diagrams defined in terms of their bisectors. We impose suitable conditions on these bisectors so that any abstract diagram can be built as the minimization diagram of some distance functions, thus showing that the class of abstract diagrams is the same as the class of Voronoi diagrams. Given a class of bisectors, such as affine or spherical bisectors, we then consider the inverse problem of determining a small class of distance functions that allows to build any diagram having such bisectors. We use a linearization technique to study this question.

Abstract Diagrams

Voronoi diagrams have been defined (see Sect. 1.2) as the minimization diagram of a finite set of continuous functions {δ 1 , . . . , δ n }. It is convenient to interpret each δ i as the distance function to an abstract object o i , i = 1, . . . , n. We define the bisector of two objects o i and

o j of O = {o 1 , . . . , o n } as b ij = {x ∈ R d , δ i (x) = δ j (x)}.
i ij as bi ij = b i ij ∪ b ij .
In a way similar to Klein [START_REF] Klein | Concrete and Abstract Voronoi Diagrams[END_REF], we now define diagrams in terms of bisectors instead of distance functions. Let B = {b ij , i = j} be a set of closed (d -1)manifolds without boundary. We always assume in the following that b ij = b ji for all i = j. We assume further that, for all distinct i, j, k, the following incidence condition (I.C.) holds:

b ij ∩ b jk = b jk ∩ b ki (I.C.)
This incidence condition is obviously needed for B to be the set of bisectors of some distance functions. By Jordan's theorem, each element of B subdivides R d into at least two connected components and crossing a bisector b ij implies moving into another

connected component of R d \ b ij . Hence, once a connected component of R d \ b ij
is declared to belong to i, the assignments of all the other connected components of R d \ b ij to i or j are determined.

Given a set of bisectors B = {b ij , i = j}, an assignment on B associates to each connected component of R d \ b ij a label i or j so that two adjacent connected components have different labels.

Once an assignment on B is defined, the elements of B are called oriented bisectors.

Given B, let us now consider such an assignment and study whether it may derive from some distance functions. In other words, we want to know whether there exists a set ∆ = {δ 1 , . . . , δ n } of distance functions such that 1. the set of bisectors of ∆ is B;

2. for all i = j, a connected component C of R d \ b ij is labeled by i if and only if ∀x ∈ C, δ i (x) ≤ δ j (x).
We define the region of object o i as ∩ j =i bi ij . A necessary condition for the considered assignment to derive from some distance functions is that the regions of any subdiagram cover R d . We call this condition the assignment condition (A.C.):

∀I ⊂ {1, . . . , n}, ∪ i∈I ∩ j∈I\{i} bi ij = R d (A.C.)
Given a set of bisectors B = {b ij , i = j} and an assignment satisfying I.C. and A.C., the abstract diagram of O is the subdivision of R d consisting of the regions of the objects of O and of their faces. The name abstract Voronoi diagram was coined by Klein [START_REF] Klein | Concrete and Abstract Voronoi Diagrams[END_REF], referring to similar objects in the plane.

For any set of distance functions δ i , we can define the corresponding set of oriented bisectors. Obviously, I.C. and A.C. are satisfied and the abstract diagram defined by this set is exactly the minimization diagram for the distance functions δ i . Hence any Voronoi diagram allows us to define a corresponding abstract diagram. Let us now prove the converse: any abstract diagram can be constructed as a Voronoi diagram.

Specifically, we prove that I.C. and A.C. are sufficient conditions for an abstract diagram to be the minimization diagram of some distance functions, thus proving the equivalence between abstract diagrams and Voronoi diagrams. We need the following technical lemmas. Lemma 2. The assignment condition implies that for any distinct i, j, k, we have

b j ij ∩ b k jk ∩ b i ki = ∅. Proof. A.C. implies that R d = ∪ 1≤i≤n ∩ j =i bi ij ⊂ bi ij ∪ bj jk ∪ bk ki . Hence, bi ij ∪ bj jk ∪ bk ki = R d . Taking the complementary sets, we obtain b j ij ∩ b k jk ∩ b i ki = ∅.
Lemma 3. For any distinct i, j, k, we have We can now prove a lemma stating a transitivity relation:

b ij ∩ b k jk ⊂ b k ik and b ij ∩ bk jk ⊂ bk ik (1.3) b ij ∩ b j jk ⊂ b i ik and b ij ∩ bj jk ⊂ bi ik (1.4) Proof. Let us first prove that b ij ∩ b k jk ⊂ b k ik : Consider x ∈ b ij ∩ b k jk . Assume,
Lemma 4. For any distinct i, j, k, we have b i ij ∩ b j jk ⊂ b i ik . Proof. Let x ∈ b i ij ∩ b j jk . Assume, for a contradiction, that x ∈ b i ik . If x ∈ b k ik , we have b i ij ∩ b j jk ∩ b k ik = ∅, contradicting Lemma 2.
Therefore, x has to belong to b ik , which implies that x ∈ b i ij ∩ b ik ⊂ b k kj by Lemma 3. This contradicts x ∈ b j jk . We deduce that x ∈ b i ik , as needed.

The following lemma states that at most two assignments are likely to derive from some Voronoi diagram. All other assignments are determined in a similar way. One can easily see that reversing the sides of b 12 reverses all the assignments. Thus, we have at most two possible global assignments.

Theorem 12. Given a set of bisectors B = {b ij , 1 ≤ i = j ≤ n} that satisfies the incidence condition (I.C.) and an assignment that satisfies the assignment condition (A.C.), there exists a set of distance functions {δ i , 1 ≤ i ≤ n} defining the same bisectors and assignments.

Proof. Let δ 1 be any real continuous function over R d . Let j > 1 and assume the following induction property: for all i < j, the functions δ i have already been constructed so that

∀i, i < j, δ i (x) ≤ δ i (x) ⇔ x ∈ bi ii .
Let us build δ j . We consider the arrangement of all bisectors b ij , for i < j: for each I ⊂ J = {1, . . . , j -1}, we define V I = (∩ i∈I bi ij ) ∩ (∩ k∈J\I bj jk ). The set V I is a non necessarily connected region of the arrangement where we need δ j > δ i if i ∈ I and δ j < δ i if i ∈ J \ I. This leads us to the following construction.

The interior of V I is int

V I = (∩ i∈I b i ij ) ∩ (∩ k∈J\I b j jk )
. Lemma 4 and the induction hypothesis imply that

∀i ∈ I, ∀k ∈ J \ I, ∀x ∈ int V I , δ i (x) < δ k (x).
In particular, if we define ν I = min k∈J\I δ k and µ I = max i∈I δ i on V I , we have µ I < ν I on int V I .

Let us now consider some point x on the boundary of V I . We distinguish two cases. We can first assume that x ∈ b ij for some i ∈ I. Then, by Lemma 3, for any i

∈ I \ {i}, x ∈ b ij ∩ bi i j ⊂ bi i i so that δ i (x) ≤ δ i (x). It follows that µ I (x) = δ i (x).
Consider now the case when x ∈ ∂V I ∩ b jk with k ∈ J \ I, we have

ν I (x) = δ k (x). Finally, if x ∈ ∂V I ∩ b ij ∩ b jk with i ∈ I and k ∈ J \ I, we have µ I (x) = δ i (x) and ν I (x) = δ k (x). By the induction hypothesis, δ i (x) = δ k (x), which implies that µ I (x) = ν I (x).
It follows that we can define a continuous function ρ on ∂V I in the following way:

ρ I (x) = µ I (x) if ∃i ∈ I, x ∈ b ij = ν I (x) if ∃k ∈ J \ I, x ∈ b jk Furthermore, on ∂V I ∩ b ij = ∂V I\{i} ∩ b ij , if i ∈ I, we have ρ I (x) = µ I (x) = ν I\{i} (x) = ρ I\{i} (x).
(1.5)

The definitions of the ρ I are therefore consistent, and we can now use these functions to prove that the following definition of δ j satisfies the induction property. Finally, we require δ j to be any continuous function verifying

µ I < δ j < ν I on each int V I . By continuity of δ j , we deduce from 1.5 that if x ∈ ∂V I ∩ b jk = ∂V I\{i} ∩ b ij with k ∈ J \ I, we have ρ I (x) = µ I (x) = ν I\{i} (x) = ρ I\{i} (x) = δ j (x). It follows that on each V I , for all i < j, δ i (x) < δ j (x) iff x ∈ b i ij and δ i (x) = δ j (x) iff x ∈ b ij . The induction follows.
One can prove that, in the proof of Lemma 5, the assignment we build satisfies the consequences of A.C. stated in Lemmas 2, 3 and 4. The proof of Theorem 12 does not need A.C. but only the consequences of A.C. stated in those three lemmas. It follows that any of the two possible assignments determined in the proof of Lemma 5 allows the construction of distance functions, as in Theorem 12, which implies that A.C. is indeed verified. We thus obtain a stronger version of Lemma 5. Lemma 6. For a given set B satisfying I.C. and assuming that we never have b ij ⊂ b ik for j = k, there are exactly two ways of labeling the connected components of each R d \ b ij as b i ij and b j ij such that A.C. is verified.

Theorem 12 proves the equivalence between Voronoi diagrams and abstract diagrams by constructing a suitable set of distance functions. In the case of affine bisectors, the following result of Aurenhammer [6] allows us to choose the distance functions in a smaller class than the class of continuous functions.

Theorem 13. Any abstract diagram of R d with affine bisectors is identical to the power diagram of some set of spheres of R d .

Proof. In this proof, we first assume that the affine bisectors are in general position, i.e. four of them cannot have a common subspace of co-dimension 2: the general result easily follows by passing to the limit.

Let B = {b ij , 1 ≤ i = j ≤ n} be such a set. We identify R d with the hyperplane x d+1 = 0 of R d+1 . Assume that we can find a set of hyperplanes {H i , 1 ≤ i ≤ n} of R d+1 such that the intersection H i ∩ H j projects onto b ij . Sect. 1.3 then shows that the power diagram of the set of spheres {σ i , 1 ≤ i ≤ n} obtained by projecting the intersection of paraboloid Q with each H i onto R d admits B as its set of bisectors2 (see Fig. 1

.5).

Let us now construct such a set of hyperplanes, before considering the question of the assignment condition.

Let H 1 and H 2 be two non-vertical hyperplanes of R d such that H 1 ∩ H 2 projects vertically onto b 12 . We now define the H i for i > 2: let ∆ 1 i be the maximal subspace of H 1 that projects onto b 1i and let ∆ 2 i be the maximal subspace of H 2 that projects onto b 2i . Both ∆ 1 i and

∆ 2 i have dimension d -1. I.C. implies that b 12 ∩ b 2i ∩ b i1 has co-dimension 2 in R d . Thus ∆ 1 i ∩ ∆ 2 i
, its preimage on H 1 (or H 2 ) by the vertical projection, has the same dimension d-2. This proves that ∆ 1 i and ∆ 2 i span a hyperplane H i of R d+1 . The fact that H i = H 1 and H i = H 2 easily follows from the general position assumption.

We still have to prove that H i ∩ H j projects onto b ij for i = j > 2. I.C. ensures that the projection of H i ∩ H j contains the projection of H i ∩ H j ∩ H 1 and the projection of H i ∩H j ∩H 2 , which are known to be b ij ∩b 1i and b ij ∩b 2i , by construction. The general position assumption implies that there is only

one hyperplane of R d , namely b ij , containing both b ij ∩ b 1i and b ij ∩ b 2i . This is the projection of H i ∩ H j .
As we have seen, building this set of hyperplanes of R d+1 amounts to building a family of spheres whose power diagram admits B as its set of bisectors. At the beginning of the construction, while choosing H 1 and H 2 , we may obtain any of the two possible labellings of the sides of b 12 . Since there is no other degree of freedom, this choice determines all the assignments. Lemma 5 shows that there are at most two possible assignments satisfying A.C., which proves we can build a set of spheres satisfying any of the possible assignments. The result follows.

Exercise 22. Consider the diagram obtained from the Euclidean Voronoi diagram of n points by taking the other assignment. Characterize a region in this diagram in terms of distances to the points and make a link with Exercise 3.

Inverse Problem

We now assume that each bisector is defined as the zero-set of some real-valued function over R d , called a bisector-function in the following. Let us denote by B the set of bisector-functions. By convention, for any bisector-function β ij , we assume that

b i ij = {x ∈ R d : β ij (x) < 0} and b j ij = {x ∈ R d : β ij (x) > 0}.
We now define an algebraic equivalent of the incidence relation in terms of pencil of functions: we say that B satisfies the linear combination condition (L.C.C.) if, for any distinct i, j, k, β ki belongs to the pencil defined by β ij and β jk , i.

i ij = φ -1 (c i ij )
where b i ij denotes the set of points closer to i than to j in D and c i ij denotes the set of points closer to i than to j in D .

Theorem 14. Let B = {β ij } be a set of real-valued bisector-functions over R d satisfying L.C.C. and A.C. Let V be any vector space of real functions over R d that contains B and constant functions.

If N is the dimension of V , the diagram defined by B is the pullback by some continuous function of an affine diagram in dimension N -1.

More explicitly, there exist a set C = {ψ ij • X + c ij } of oriented affine hyperplanes of R N -1 satisfying I.C. and A.C. and a continuous function φ :

R d → R N -1 such that for all i = j, bi ij = {x ∈ R d , β ij (x) ≤ 0} = φ -1 {y ∈ R N -1 , ψ ij (x) ≤ c ij }.
Proof. Let (γ 0 , . . . , γ N -1 ) be a basis of V such that γ 0 is the constant function equal to 1. Consider the evaluation application,

φ : x ∈ R d → (γ 1 (x), . . . , γ N -1 (x)) ∈ R N -1 .
If point x belongs to some b i ij , we have β ij (x) < 0. Furthermore, there exists real coefficients λ 0 ij , . . . , λ

N -1 ij such that β ij = N -1 k=0 λ k ij γ k . The image φ(x) of x thus belongs to the affine half-space B i ij of R N -1 of equation N -1 k=1 λ k ij X k < -λ 0 ij .
In this way, we can define all the affine half-spaces B i ij of R N -1 for i = j: B ij is an oriented affine hyperplane with normal vector (λ 1 ij , . . . , λ

ij = {x ∈ R d , β ij (x) ≤ 0} = φ -1 {y ∈ R N -1 , B ij (x) ≤ -λ n ij } (1.6)
Finally, let us prove that A.C. is also satisfied. Lemma 6 states that the B ij have exactly two inverse assignments satisfying A.C. Furthermore, Equation 1.6 implies that any of these two assignments defines an assignment for the b ij that also satisfies A.C. It follows that if the current assignment did not satisfy A.C., there would be more than two assignments for the b ij that satisfy A.C. This proves that A.C. is also satisfied by the B ij and concludes the proof.

We can now use Theorem 13 and specialize Theorem 14 to the specific case of diagrams whose bisectors are hyperspheres or hyperquadrics, or, more generally, to the case of diagrams whose class of bisectors spans a finite dimensional vector space.

Theorem 15. Any abstract diagram of R d with spherical bisectors such that the corresponding degree 2 polynomials satisfy L.C.C. is a Möbius diagram.

Proof. Since the spherical bisectors satisfy L.C.C., we can apply Theorem 14 and Theorem 13. Function φ of Theorem 14 is simply the lifting mapping x → (x, x 2 ), and we know from Theorem 13 that our diagram can be obtained as a power diagram pulled-back by φ. That is to say δ i (x) = Σ i (φ(x)), where Σ i is a hypersphere in R d+1 .

Another way to state this transformation is to consider the diagram with spherical bisectors in R d as the projection by φ -1 of the restriction of the power diagram of the hyperspheres

Σ i to the paraboloid φ(R d ) ⊂ R d+1 of equation x d+1 = x 2 .
Assume that the center of Σ j is (u j 1 , . . . , u j d+1 ), and that the squared radius of Σ j is w j . We denote by Σ j the power to Σ j . Distance δ j can be expressed in terms of these parameters:

δ j (x) = Σ j (φ(x)) = 1≤i≤d (x i -u j i ) 2 + (( 1≤i≤d x 2 i ) -u j d+1 ) 2 -w j .
Subtracting from each δ j the same term ( 1≤i≤D x 2 i ) 2 leads to a new set of distance functions that define the same minimization diagram as the δ j . In this way, we obtain new distance functions which are exactly the ones defining Möbius diagrams.

This proves that any diagram whose bisectors are hyperspheres can be constructed as a Möbius diagram.

The proof of the following theorem is similar to the previous one:

Incremental Voronoi Algorithms

Incremental constructions consist in adding the objects one by one in the Voronoi diagram, updating the diagram at each insertion. Incremental algorithms are well known and highly popular for constructing Euclidean Voronoi diagrams of points and power diagrams of spheres in any dimension. Because the whole diagram can have to be modified at each insertion, incremental algorithms have a poor worst-case complexity. However most of the insertions result only in local modifications and the worst-case complexity does not reflect the actual complexity of the algorithm in most practical situations. To provide more realistic results, incremental constructions are analyzed in the randomized framework. This framework makes no assumption on the input object set but analyzes the expected complexity of the algorithm assuming that the objects are inserted in random order, each ordering sequence being equally likely. The following theorem, whose proof can be found in many textbooks (see e.g. [START_REF] Boissonnat | Algorithmic Geometry[END_REF]) recalls that state-of-the-art incremental constructions of Voronoi diagrams of points and power diagrams have an optimal randomized complexity. .

Owing to the linearization techniques of Sec. 1.5, this theorem yields complexity bounds for the construction of linearizable diagrams such as Möbius, anisotropic or Apollonius diagrams. Incremental constructions also apply to the construction of Voronoi diagrams for which no linearization scheme exists. This is for instance the case for the 2-dimensional Euclidean Voronoi diagrams of line segments. The efficiency of the incremental approach merely relies on the fact that the cells of the diagram are simply connected and that the 1-skeleton of the diagram, (i. e. the union of its edges and vertices) is a connected set. Unfortunately, these two conditions are seldom met except for planar Euclidean diagrams. Let us take Apollonius diagrams as an illustration. Each cell of an Apollonius diagram is star shaped with respect to the center of the associated sphere and is thus simply connected. In the planar case, Apollonius bisectors are unbounded hyperbolic arcs and the 1-skeleton can easily be made connected by adding a curve at infinity. The added curve can be seen as the bisector separating any input object from an added fictitious object. In 3-dimensional space, the skeleton of Apollonius diagrams is not connected: indeed, we know from Sect. 1.4.3 that the faces of a single cell are in 1-1 correspondence with the faces of a 2-dimensional Möbius diagram and therefore may include isolated loops.

As a consequence, the rest of this section focuses on planar Euclidean diagrams. After some definitions, the section recalls the incremental construction of Voronoi diagrams, outlines the topological conditions under which this approach is efficient and gives some examples. The efficiency of incremental algorithms also greatly relies on the availability of some point location data structure to answer nearest neighbor queries. A general data structure, the Voronoi hierarchy, is described at the end of the section. The last subsection lists the main predicates involved in the incremental construction of Voronoi diagrams.

Planar Euclidean diagrams

To be able to handle planar objects that possibly intersect, the distance functions that we consider in this section are signed Euclidean distance functions, i.e. the distance δ i (x) from a point x to an object o i is:

δ i (x) = min y∈ōi y -x , if x ∈ o -min y∈ō c i y -x , if x ∈ o
where ōi is the closure of o i and ōc i the closure of the complement of o i . Note that the distance used to define Apollonius diagrams matches this definition. Then, given a finite set O of planar objects and o i ∈ O, we define the Voronoi region of o i as the locus of points closer to o i than to any other object in O

V (o i ) = {x ∈ R 2 : δ i (x) ≤ δ j (x), ∀o j ∈ O}.
Voronoi edges are defined as the locus of points equidistant to two objects O and closer to these two objects than to any other object in O, and Voronoi vertices are the locus of points equidistant to three or more objects and closer to these objects than to any other object in O. The Voronoi diagram Vor(O) is the planar subdivision induced by the Voronoi regions, edges and vertices.

The incremental construction described below relies on the three following topological properties of the diagram that are assumed to be met for any set of input objects:

1. The diagram is assumed to be a nice diagram, i. e. a diagram in which edges and vertices are respectively 0 and 1-dimensional sets. 2. The cells are assumed to be simply connected. 3. The 1-skeleton of the diagram is connected.

Owing to Euler formula, Properties 1 and 2 imply that the Voronoi diagram of n objects is a planar map of complexity O(n). Property 3 is generally not granted for any input set. Think for example of a set of points on a line. However, in the planar case, this condition can be easily enforced as soon as Properties 1 and 2 are met. Indeed, if the cells are simply connected, there is no bounded bisector and the 1-skeleton can be connected by adding a curve at infinity. The added curve can be seen as the bisector separating any input object from an added fictitious object. The resulting diagram is called the compactified version of the diagram.

Incremental Construction

We assume that the Voronoi diagram of any input set we consider is a nice diagram with simply connected cells and a connected 1-skeleton. Each step of the incremental construction takes as input the Voronoi diagram Vor(O i-1 ) of a current set of objects O i-1 and an object o i ∈ O i-1 , and aims to construct

the Voronoi diagram Vor(O i ) of the set O i = O i-1 ∪ {o i }.
In the following, we note V (o, O i-1 ) the region of an object o in the diagram Vor(O i-1 ) and

V (o, O i ) the region of o in Vor(O i ).
We note Skel(O i-1 ) the 1-skeleton of Vor(O i-1 ). Let x be a point of Skel(O i-1 ). We note N (x, O i-1 ) the nearest neighbors of x in O i-1 , i.e. the subset of objects of O i-1 that are closest to x. The point x is said to be in conflict with o i if x is closer to o i than to N (x, O i-1 ). Hence, the part of the skeleton that conflicts with o i , called the conflict skeleton for short, is exactly the intersection of the skeleton Skel(O i-1 ) with the region of o i in Vor(O i ). See Fig. 1.11.

The conflict skeleton is a subgraph of Skel(O i-1 ) and the endpoints of this subgraph are the vertices of the region V (o i , O i ). If V (o i , O i ) is not empty, the conflict skeleton is not empty either. Indeed, an empty conflict skeleton would imply that

V (o i , O i ) is included in a single region V (o, O i-1 ) of the diagram Vor(O i-1
) and the region V (o, O i ) would not be simply connected. Furthermore, the following lemma, due to Klein et al. [START_REF] Klein | Randomized incremental construction of abstract Voronoi diagrams[END_REF], proves that the conflict skeleton is a connected subgraph of Skel(O i-1 ). Lemma 7. If, for any input set, the Voronoi diagram is a nice diagram with simply connected regions and a connected 1-skeleton, the conflict skeleton of an additional object is connected.

Proof. We use the above notation and assume, for a contradiction, that the conflict skeleton of o i , which is Skel(O i-1 ) ∩ V (o i , O i ), consists of several disjoint connected components Sk 1 , Sk 2 , . . . , Sk . Each connected component Sk j has to intersect the boundary of the new region V (o i , O i ), otherwise Skel(O i-1 ) would not be connected, a contradiction. Then, if ≥ 2, there exists a path C in V (o i , O i ) connecting two points x and y on the boundary of V (o i , O i ) and separating Sk 1 and Sk 2 , see Fig. 1.10. The path C does not intersect Skel(O i-1 ) and is therefore included in the region V (o, O i-1 ) of some object o of O i-1 . Since, points arbitrarily close to x and y but outside

V (o i , O i ) belong to V (o, O i-1
), x and y can be joined by a simple path

D in V (o, O i ) ⊂ V (o, O i-1 ). The simple closed curve C ∪ D is contained in V (o, O i-1
) and encloses Sk 1 or Sk 2 , which contradicts the fact that Skel(O i-1 ) is connected.

Once the conflict skeleton is known, the Voronoi diagram Vor(O i-1 ) can be updated, leading to Vor(O i ). This is done by Procedure 2.

Procedure 3 describes a step of the incremental construction.

In the sequel, the incremental construction is analyzed in the randomized setting. It is assumed that each object has constant complexity, which implies that each operation involving a constant number of objects is performed in constant time. Because the conflict skeleton is connected, Substep 2 can be performed by traversing the graph Skel(O i-1 ) in time proportional to the number of edges involved in the conflict skeleton. These edges will be either deleted or shortened in the new diagram. Substep 3 takes time proportional

Procedure 2 Updating the Voronoi diagram

Input: Vor(Oi-1), Skel(Oi-1)

1. Create a new vertex at each endpoint of the conflict skeleton; 2. Remove vertices, edges and portions of edges that belong to Skel(Oi-1) ∩ V (oi, Oi); 3. Connect the new vertices s as to form the boundary of the new region.

Output: Vor(Oi)

Procedure 3 A step of the incremental algorithm Input: Vor(Oi-1) and a new object oi 1. Find a first point x of Skel(Oi-1) in conflict with oi ; 2. Compute the whole conflict skeleton ; 3. Update Vor(Oi-1) into Vor(Oi) using Procedure 2; 4. Update the location data structure.

Output: Vor(Oi)
to the number of edges involved in the conflict skeleton plus the number of edges of V (o i , O i ). The latter are the new edges. Hence Substeps 2 and 3 take time proportional to the number of changes in the 1-skeleton. Because each edge in the skeleton is defined by four objects and because the complexity of the Voronoi diagram of n objects is O(n), a standard probabilistic analysis (see e.g. [START_REF] Boissonnat | Algorithmic Geometry[END_REF]) shows that the expected number of changes at each step of the incremental algorithm is O(1). The overall randomized complexity of the algorithm is O(n).

The costs of Substeps 1 and 4 depend of course on the type of the input objects and of the location data structure. In Sect. 1.6.3, we described a location data structure, called the Voronoi hierarchy, that can be used in the case of disjoint convex objects. The Voronoi hierarchy allows to detect a first conflict in randomized time O(log 2 n). At each step, the data structure is updated in time O(m log 2 n) where m is the number of changes in the diagram. Because the expected number of changes at each step is O(1), the expected cost for updating the hierarchy is O(log 2 n). This yields the following theorem.

Theorem 18. The incremental construction of the planar Euclidean Voronoi diagram of n disjoint convex objects with constant complexity takes O(n log 2 n) expected time.

Note that the incremental construction described here is on-line, meaning that the algorithm does not need to know the whole set of objects right from the beginning. If the whole set of objects is known in advance, localization can be made easy and the maintenance of a location data structure is no longer required [START_REF] Alt | The Voronoi diagram of curved objects[END_REF]. The idea consists in picking one witness point inside each object and in building first the Voronoi diagram of witness points. In a second phase, each witness point is replaced in turn by the corresponding object. When replacing a witness point by the corresponding object, any point on the boundary of the cell of the witness point belongs to the conflict skeleton. The algorithm is no longer on-line but its randomized complexity is reduced to O(n log n).

Voronoi Diagrams of Line Segments

The above incremental construction applies to the Voronoi diagram of disjoint line segments. Indeed, in the case of disjoint line segments, the bisector curves are unbounded simple curves, each composed of at most seven line segments and parabolic arcs (see e.g. [START_REF] Boissonnat | Algorithmic Geometry[END_REF]Chap. 19]). Hence, Voronoi vertices and edges are respectively 0 and 1-dimensional sets. Furthermore, each region in the diagram is weakly star shaped with respect to its generating segment, meaning that the segment joining any point in the region to its closest point on the associated segment is included in the region. It follows that Voronoi regions are simply connected.

If the segments are allowed to share endpoints, the Voronoi diagram exhibits 2-dimensional Voronoi edges, hence violating the definition of nice Voronoi diagrams. A way to circumvent this problem consists in considering that each segment is composed of three distinct objects: the two endpoints and the open segment. If the two endpoints of a segment are inserted in the diagram prior to the open segment, the incremental construction encounters no 2-dimensional bisecting region and the algorithm presented above can be used.

Voronoi Diagrams of Curved Segments

Voronoi diagrams of disjoint curved segments have been studied by Alt, Cheong and Vigneron [START_REF] Alt | The Voronoi diagram of curved objects[END_REF]. Alt, Cheong and Vigneron introduce the notion of harmless curved segments defined as follows. A curved segment is said to be convex when the region bounded by the curved segment and the line segment joining its endpoints is convex. A spiral arc is a convex curved segment with monotonously increasing curvature. A harmless curved segment is either a line segment or a circular arc or a spiral arc. It can be shown that, if the input curved segments are split into harmless sub-segments and if each open curved sub-segment and its two endpoints are considered as three distinct sites, the Voronoi diagram is a nice diagram with simply connected regions. The incremental construction paradigm described above therefore applies.

Voronoi Diagrams of Convex Objects

The case of disjoint smooth convex objects is quite similar to the case of disjoint segments. The bisecting curves between two such objects is a 1dimensional curve. Furthermore, each Voronoi region is weakly star shaped with respect to the medial axis of its object, hence simply connected. Therefore, the Voronoi diagram can be built using the incremental algorithm. Note that the Voronoi diagram of disjoint smooth convex objects could also be obtained by applying the incremental algorithm to the curved segments forming the boundaries of the objects. However, this approach requires the subdivision of the boundary of each object into harmless parts and yields a Voronoi diagram which is a refinement of the diagram of the input objects.

If we still assume the objects to be smooth and convex but allow them to intersect, things become more difficult. Karavelas and Yvinec [29] have shown that the Voronoi regions remain simply connected if and only if the objects of O are pseudo-disks, meaning that the boundaries of any two objects of O intersect in at most two points. The above incremental algorithm can be adapted to work in this case. However, since the distance is a signed distance, some sites may have an empty region, which makes the algorithm slightly more complicated. Note that this may only happen when some of the objects are included in the union of others. The algorithm has to check that the new object has a non-empty region and must handle the case where the insertion of a new object causes the region of some other object to vanish. Karavelas and Yvinec [29] showed that there is no use to maintain a location data structure in this case because each insertion takes linear time anyway.

The algorithm can be generalized to the case of convex objects with piecewise smooth pseudo-circular boundaries. As in the case of segments, the main problem comes from the fact that sharp corners on the boundaries of objects yield 2-dimensional bisectors. This problem can be handled as in the case of line segments and planar curves by considering each corner as an object on its own.

The Voronoi Hierarchy

The first step when inserting a new object o i consists in finding one point of the current skeleton Skel(O i-1 ) in conflict with o i . If the objects do not intersect, this is done by searching the object o of O i-1 nearest to a point x of o i . Indeed, if the objects do not intersect, x belongs to the region of o in Vor(O i-1 ) and to the region of o i in Vor(O i ). Therefore o i has to be a neighbor of o in Vor(O i ) and some point on the boundary of V (o, O i-1 ) is in conflict with o i . If the objects intersect, things are slightly more complicated but nearest object queries can still be used to find out whether the new object is hidden and, if not, to find a first conflict (see [29]).

Let us describe now how to find the object o of a set O closest to a query point x. A simple strategy is to perform a walk in the Voronoi diagram Vor(O). The walk starts at any region of the diagram. When the walk visits the region V (o) of an object o it considers in turn each of the neighboring regions. If one of the neighbors of o, say o , is closer to x than o, the walk steps to the region V (o ). If none of the neighbors of o in Vor(O) is closer to x than o, then o is the object closest to x and the walk ends. Because the distance between x and the objects of the visited regions is decreasing, the walk cannot loop and is bound to end. However, the walk may visit all the regions before ending. The Voronoi hierarchy [29] is a randomized data structure that makes this strategy more efficient. The Voronoi hierarchy can be considered as a 2-dimensional version of the skip lists introduced by Pugh [44] and generalizes the Delaunay hierarchy described in [START_REF] Devillers | The Delaunay hierarchy[END_REF].

For a set of objects O, the Voronoi hierarchy H V (O) is a sequence of Voronoi diagrams Vor(Θ ), = 0, . . . , L, built for subsets of O forming a hierarchy, i.e,

O = Θ 0 ⊇ Θ 1 ⊇ • • • ⊇ Θ L .
The hierarchy H V (O) is built together with the Voronoi diagram Vor(O) according to the following rules:

1. Every object of O is inserted in Vor(Θ 0 ) = Vor(O); 2. An object o that has been inserted in Vor(Θ ), is inserted in Vor(Θ +1 ) with probability β.

To answer nearest object queries, the Voronoi hierarchy works as follows. Let us call θ the object of Θ closest to the query point x. First, a simple walk is performed in the top-most diagram to find θ L . Then, at each level = L -1, . . . , 0, a simple walk is performed in Vor(Θ ) from θ +1 to θ (see Fig. 1.12).

x Fig. 1.12. Locating x using the Voronoi hierarchy.

It is easy to show that the expected size of H V (O) is O( n 1-β ), and that the expected number of levels in H V (O) is O(log 1/β n). Moreover, the following lemma proves that the expected number of steps performed by the walk at each level is constant. Lemma 8. Let x be a point in the plane. Let θ (resp. θ +1 ) be the object closest to x in Θ (resp. Θ +1 ). Then the expected number of Voronoi regions visited during the walk in Vor(Θ ) from θ +1 to θ is O(1/β).

Proof. The objects whose regions are visited at level are closer to x than θ +1 . Consequently, if, among the objects of Θ , θ +1 is the k-th closest object to x, the walk in Vor(Θ ) performs at most k steps.

Let us show that θ +1 is the k-th closest object to x in Θ with probability β(1β) k-1 . Such a case occurs if and only if the two following conditions are satisfied:

1. Object θ +1 has been inserted in Θ +1 2. None of the k -1 objects of Θ that are closer to x than θ +1 has been inserted in Θ +1 .

The first condition occurs with probability β and the second with probability (1β) k-1 . Let n l be the number of objects in Θ . The expected number N of objects that are visited at level is bounded as follows:

N ≤ n k=1 k(1 -β) k-1 β < β ∞ k=1 k(1 -β) k-1 = 1 β ,
We still have to bound the time spent in each of the visited regions. Let o be the site of a visited region in Vor(Θ ). It is not efficient to consider in turn each neighbor o of o in Vor(Θ ) and compare the distances from x to o and o . Indeed, since the complexity of each region in the Voronoi diagram Vor(Θ ) may be Ω(n ), this would imply that the time spent at each level of the hierarchy is O(n), yielding a total of O(n) time per insertion. To avoid this cost, a balanced binary tree is attached to each Voronoi region in the Voronoi Exercise 24. Show that the planar Euclidean Voronoi diagram of n points can be computed on line in O(n log n) time (Hint: in the case of points, using a Delaunay hierarchy instead of the Voronoi hierarchy, nearest neighbor queries can be answered in O(log n). Upon insertion the structure is updated in randomized time O(log n). See [START_REF] Devillers | The Delaunay hierarchy[END_REF]) Exercise 25. Show that the planar Euclidean Voronoi diagram of n line segments can be computed on line in O(n log n) time. See e.g. [START_REF] Boissonnat | Algorithmic Geometry[END_REF] for a solution.

Exercise 26. Show that the planar euclidean diagram of n disjoint convex objects can be computed using predicates that involve only four objects.

Exercise 27. Provide a detailed description of the predicates of the incremental Voronoi diagram construction and a way to implement them efficiently for various types of simple objects (e.g., line segments, circles). See [31,30].

Exercise 28 (2D Abstract Voronoi diagrams). Klein et al. [START_REF] Klein | Randomized incremental construction of abstract Voronoi diagrams[END_REF] have defined abstract Voronoi diagrams in dimension 2 using bisecting curves. Each bisecting curve b ij is assumed to be an infinite curve separating the plane in two regions affected respectively to o i and o j and the Voronoi regions are defined as in Sect. 1.5.1. Klein et al. assume that the affectation fulfills the assignment condition but they do not assume the incidence condition. Instead they assume that each pair of bisecting curves intersect in only a finite number of connected components and that the interior of Voronoi regions are pathconnected. Show that, under Klein

Medial Axis

In this section, we introduce the concept of Medial Axis of a bounded set Ω, which can be seen as an extension of the notion of Voronoi diagram to infinite sets. Interestingly, it is possible to construct certified approximations of the medial axis of quite general sets efficiently. One approach to be described consists in sampling the boundary of Ω and then computing an appropriate subset of the Voronoi diagram of the sample which approximates the medial axis. Hence the problem of approximating the medial axis of Ω boils down to sampling the boundary of Ω, a problem that is closely related to mesh generation (see Chap. ??). Other informations on the medial axis can be found in Chap. ??.

Medial Axis and Lower Envelope

The medial axis of an open set Ω, denoted by M(Ω), is defined as the the set of points of Ω that have more than one nearest neighbor on the boundary of Ω. Nearest refers in this section to the Euclidean distance although the results may be extended to other distance functions. A medial sphere σ is a sphere centered at a point c of the medial axis and passing through the nearest neighbors of c on ∂Ω. Those points where σ is tangent (in the sense that σ does not enclose any point of ∂Ω) to ∂Ω are called the contact points of σ.

The concept of medial axis can be considered as an extension of the notion of Voronoi diagram to infinite sets. Let o be a point of the boundary of Ω and δ o be the distance function to o defined over

Ω ∀x ∈ Ω : δ o (x) = x -o .
The lower envelope of the infinite set of functions δ o is defined as

∆ -= inf o∈∂Ω δ o .
Following what we did for Voronoi diagrams (Sect. 1.2), we define, for any point x ∈ Ω, its index set I(x) as the set of all o such that ∆ -(x) = δ o (x). The set of points x such that |I(x)| > 1 constitutes the medial axis of Ω.

Computing the medial axis is difficult in general. If Ω is defined as a semialgebraic set, i.e. a finite collection of algebraic equations and inequalities, M(Ω) is also a semi-algebraic set that can therefore be computed using techniques from real algebraic geometry [START_REF] Attali | Stability and computation of medial axes: a state of the art report[END_REF]8]. This general approach, however, leads to algorithms of very high complexity. Theorem 1 can also be used but, still, working out the algebraic issues is a formidable task. Effective implementations are currently limited to simple objects. If Ω is a planar domain bounded by line segments and circular arcs, one can apply the results of Sect. 1.6. Further results can be found in [24].

An alternative and more practical approach consists in departing from the requirement to compute the medial axis exactly. In Sect. 1.7.2, we describe a method that approximates the medial axis of an object by first sampling its boundary, and then computing and pruning the Voronoi diagram of the sample.

Approximation of the Medial Axis

Approximating the medial axis of a set is a non trivial issue since sets that are close for the Hausdorff distance may have very different medial axes. This is illustrated in Fig. 1.13. Let S be a closed curve, Ω = R 2 \ S and P be a finite set of points approximating S. As can be seen on the figure, the skeleton of the Voronoi diagram of P is far from the medial axis of Ω since there are many long branches with no counterpart in M(Ω). These branches are Voronoi edges whose dual Delaunay edges are small (their lengths tend to 0 when the sampling density increases). In other words, the medial axis is not continuous under the Hausdorff distance. Notice however that if we remove the long branches in Fig. 1.13, we obtain a good approximation of the medial axis of S. This observation will be made precise in Lemma 10 below. It leads to an approximate algorithm that first sample S and extract from the Euclidean Voronoi diagram of the sample a sub-complex that approximates the medial axis of S. Given an open set Ω and a point x on the medial axis of Ω, we define D(x) as the diameter of the smallest closed ball containing the contact points of the medial sphere centered at x. We define the λ-medial axis of Ω, denoted M λ (Ω) as the subset of the medial axis of Ω consisting of points x such that D(x) ≥ λ.

Let Ω and Ω be two open sets and let S and S denote their boundaries. We assume that S and S are compact and that their Hausdorff distance d H (S, S ) is at most ε: any point of S is at distance at most ε from a point of S and vice versa. Notice that we do not specify S nor S to be finite or infinite point sets. For convenience, we rename medial axis of S (resp. S ) and write M(S) (resp. M(S )) the medial axis of R 3 \ S. Similarly, we rename medial axis of S and write M(S )) the medial axis of R 3 \ S.

The following lemma says that the λ-medial axis of S is close to the medial axis of S , provided that λ is sufficiently large. It should be emphasized that close here refers to the one-sided Hausdorff distance: the medial axis of S is not necessarily close to the λ-medial axis of S although, by exchanging the roles of S and S , the lemma states that the λ -medial axis of S is close to the medial axis of S for a sufficiently large λ . We will go back to this point later.

We say that a ball is S-empty if its interior does not intersect S. The sphere bounding a S-empty ball is called a S-empty sphere. Let σ be the maximal S -empty sphere tangent to σ at y . σ is tangent to S at at least two points. Let c be its center and r its radius. and, since, by assumption, r > 2ε, we get cc (r -2εh) ≤ 2ε(rε). Moreover, by assumption, l ≥ εr(1 -ε r ) which implies r -2εh ≥ 0. Indeed, h = r 2 -4l 2 ≤ r 2 -4εr + 4ε 2 = r -2ε.

We then deduce

c -c ≤ 2ε(r -ε) r -h -2ε ≤ 2εr r - √ r 2 -4l 2 -2ε .
We then get

c -c r ≤ 2ε(r + √ 2 -4l 2 -2ε) (r -2ε) 2 -(r 2 -4l 2 ) ≤ εr l 2 -εr + ε 2 .
The same bound plainly holds for |1 -r r |.

We consider the case where S is a surface of R 3 and S is a finite set of points on S at Hausdorff distance at most ε from S. To avoid confusion, we rename S as P. As already noticed, M(P) is the 1-skeleton of the Voronoi diagram Vor(P), called simply the Voronoi diagram of P in this section. The set of contact points of a medial sphere centered at a point c ∈ M(P) is the set of closest points of c in P. Any point in the relative interior of a face of Vor(P) has the same closest points in P. It follows that the λ-medial axis of P is the subset of the faces of Vor(P) whose contact points cannot be enclosed in a ball of diameter λ.

Lemma 10 then says that any Delaunay sphere of Del(P) passing through two sample points that are sufficiently far apart, is close to a medial sphere of S (for the Hausdorff distance). We have therefore bounded the one-sided Hausdorff distance from the λ-medial axis (Voronoi diagram) of an ε-sample P of S to the medial axis of S, when λ is sufficiently large with respect to ε. If we apply the lemma the other way around, we see that, for sufficiently large λ , M λ (S) is close to M(P). However, as observed above (Fig. 1.13), we cannot hope to bound the two-sided Hausdorff distance between M(S) and M(P).

The above lemma can be strengthened as recently shown by Chazal and Lieutier [14,[START_REF] Attali | Stability and computation of medial axes: a state of the art report[END_REF]. They proved that the λ-medial axis of S is close to the λmedial axis of P for a sufficiently large λ and some positive λ that depends on λ and ε. More precisely, let D be the diameter of S and k a positive constant. They showed that there exist three functions of ε, k(ε) = 15 √ 2 Consider now a family of point sets P ε parametrized by ε such that d H (S, P ε ) ≤ ε and let ε tends to 0. Because M η (S) tends to M(S) when η tends to 0, we deduce from the above inequalities, that lim The above discussion provides an algorithm to approximate the medial axis of S within any specified error (see Algorithm 4). The main issue is therefore to compute a sample of points on S (step 1). If S is a surface of R 3 , one can use a surface mesh generator to mesh S and take for P the vertices of the mesh. Various algorithms can be found in Chap. ?? and we refer to that chapter for a thorough description and analysis of these algorithms. Especially attractive in the context of medial axis approximation, are the algorithms that are based on the 3-dimensional Delaunay triangulation since we get the Voronoi diagram of the sample points (step 2) at no additional cost. An example obtained with the surface mesh generator of Boissonnat and Oudot [12]

Algorithm 4 Approximation of the Medial Axis

Voronoi Diagrams in Cgal

The Computational Geometry Algorithms Library Cgal [1] offers severals packages to compute Voronoi diagrams. Euclidean Voronoi diagrams of points and power diagrams are represented through their dual Delaunay and regular triangulations. Cgal provides Delaunay and regular triangulations in R 2 and R 3 . The implementation is based on a randomized incremental algorithm using a variant of the Voronoi hierarchy described in Sect. 1.6. Delaunay triangulations are also provided in higher dimensions.

The library also contains packages to compute Voronoi diagrams of line segments [27] and Apollonius diagrams in R 2 [START_REF] Karavelas | Dynamic additively weighted voronoi diagrams in 2d[END_REF]. Those packages implement the incremental algorithm described in Sect. 1.6. A prototype implementation of Möbius diagrams in R 2 also exists. This prototype computes the Möbius diagram as the projection of the intersection of a 3-dimensional power diagram with a paraboloid, as described in Sect. 1.4.1. This prototype also serves as the basis for the developement of a Cgal package for 3-dimensional Apollonius diagrams, where the boundary of each cell is computed as a 2-dimensional Möbius diagram, following the results of Sect. 1.4.3 [START_REF] Boissonnat | Convex hulls and Voronoi diagrams of additively weighted points[END_REF]. See Fig. 1.8.

Applications

Euclidean and affine Voronoi diagrams have numerous applications we do not discuss here. The interested reader can consult other chapters of the book, most notably Chap. ?? on surface meshing and Chap. ?? on reconstruction. Other applications can be found in the surveys and the textbooks mentionned in the introduction.

Additively and multiplicatively weighted distances arise when modeling growing processes and have important applications in biology, ecology and other fields. Consider a number of crystals, all growing at the same rate, and all starting at the same time : one gets a number of growing circles. As these circles meet, they draw a Euclidean Voronoi diagram. In reality, crystals start growing at different times. If they still grow at the same rate, they will meet Euclidean Voronoi diagrams of non punctual objects find applications in robot motion planning [37,26]. Medial axes are used for shape analysis [25], for computing offsets in Computer-Aided Design [19], and for mesh generation [START_REF] Price | Hexahedral mesh generation by medial surface subdivision: Part I, solids with convex edges[END_REF]42,47]. Medial axes are also used in character recognition, road network detection in geographic information systems, and other applications.

  The equivalence classes R d / ≡ are relatively open sets that cover R d . Their closures are called the faces of the minimization diagram of F (see Fig. 1.1).

Fig. 1. 1 .

 1 Fig. 1.1. The lower envelope of a set of univariate functions. The minimization diagram is drawn on the horizontal line with the corresponding indices. The face of index {1} consists of two components.

Fig. 1. 2 .

 2 Fig. 1.2. The Voronoi diagram of a set of 9 points.

Theorem 2 .

 2 The faces of the Euclidean Voronoi diagram Vor(P) of a set of points P are the vertical projections of the faces of the convex polyhedron V(P).Exercise 3. Consider the maximization diagram obtained by projecting the faces of h+ p1 ∩ • • • ∩ h + pn vertically.Characterize the points that belong to a face of this diagram in terms of the distance to the points of P.

Fig. 1 . 3 .

 13 Fig. 1.3. The polyhedron V(P), with one of its faces projected onto R d .

Fig. 1. 4 .

 4 Fig. 1.4. The Delaunay triangulation of a point set (in bold) and its dual Voronoi diagram (thin lines).

Fig. 1. 5 .

 5 Fig. 1.5. The polar hyperplane of a sphere.

Theorem 3 .

 3 The combinatorial complexity of the Voronoi diagram of n points of R d and of their Delaunay triangulation is Θ n d+1 2

  For a set P of n points uniformly distributed in a ball of R d , the combinatorial complexity of the Voronoi diagram of P is O(n) where the constant depends on the dimension d [20]. Other results are known for other point distributions [3, 5, 23].

Fig. 1

 1 Fig. 1.6. A power diagram.

Theorem 5 .

 5 The combinatorial complexity of the power diagram of n hyperspheres of R d and of its dual regular triangulation are Θ n d+1 2

  The minimization diagram of the δ i is called the Möbius diagram of ω and noted Möb (ω) (see Fig. 1.4.1).

Fig. 1. 7 .

 7 Fig. 1.7. A Möbius diagram.

  Given a set ω of n Möbius sites of R d+1 , the restriction of their Möbius diagram to a hypersphere S d is called a Möbius diagram on S d . It can be shown that such a diagram is also the restriction of a power diagram of hyperspheres of R d+1 to S d (Exercise 14) . We define spherical diagrams on S d as the diagrams on S d whose bisectors are hyperspheres of S d and that satisfy two properties detailed in Sect. 1.5.1 and 1.5.2. See Exercise 16 for more details on these conditions. This exercise proves that the restriction of a Möbius diagram, i.e. a Möbius diagram on S d , is a spherical diagram. Let us now prove the converse: any spherical diagram on S d is a Möbius diagram on S d . Let h be a hyperplane of R d+1 . The stereographic projection that maps S d to h maps any spherical diagram D on S d to some spherical diagram D on h. Theorem 15 implies that this D is in fact a Möbius diagram. Exercise 12 shows that D, which is the image of D by the inverse of the stereographic projection, is the restriction of some power diagram of R d+1 to S d . Exercise 14 then proves that it is indeed a Möbius diagram. Exercise 12. Show that the linearization theorem and its corollary still hold if one replaces the paraboloid Q by any hypersphere of R d+1 and the vertical projection by the corresponding stereographic projection. Exercise 13. Show that the intersection of a Möbius diagram in R d with a k-flat or a k-sphere σ is a Möbius diagram in σ. Exercise 14. Show that the restriction of a Möbius diagram of n Möbius sites to a hypersphere Σ ⊂ R d (i.e. a Möbius diagram on Σ) is identical to the restriction of a power diagram of n hyperspheres of R d with Σ, and vice versa. Exercise 15. The predicates needed for constructing a Möbius diagram are those needed to construct Pow(Σ) and those that decide whether a face of Pow(Σ) intersects Q or not. Write the corresponding algebraic expressions. Exercise 16. Explain how the two conditions A.C. and L.C.C. presented in Sect. 1.5.1 and 1.5.2 are to be adapted to the case of spherical diagrams on a sphere (Hint: consider the L.C.C. condition as a pencil condition, and define a pencil of circles on a sphere as the intersection of a pencil of hyperplanes with this sphere). Note that the restriction of a Voronoi diagram (affine or not) to S d always satisfies this adapted version of A.C. and prove that the restriction of an affine Voronoi diagram to S d satisfies L.C.C. so that Exercise 14 allows to conclude that the restriction of a Möbius diagram to S d , i.e. Möbius diagram on S d , is a spherical diagram.

Theorem 9 .

 9 The Voronoi diagram of n anisotropic sites of R d can be computed in time O(n D+1 2

Fig. 1

 1 Fig. 1.8. The Apollonius diagram of a set of circles. Compare with the power diagram of the same set of circles in Fig. 1.6.

d 2 +1

 2 log n). The intersection involved in Step 3 can be computed in time proportional to the number of faces of the power diagram of the Σ i 's, which is O(n d 2 +1 ). We have thus proved the following theorem due to Aurenhammer [6]: Theorem 10. The Apollonius diagram of a set of n hyperspheres in R d has complexity O(n d 2 +1 ) and can be computed in time O(n d 2 +1 log n).

Fig. 1

 1 Fig. 1.9. A cell in an Apollonius diagram of hyperspheres projects vertically onto a Möbius diagram in σ0.

  It follows that the radial projection of the boundary of A(σ 0 ) onto σ 0 is a Möbius diagram on σ 0 . Such a Möbius diagram on σ 0 can be computed by constructing the restriction of the power diagram of n hyperspheres of R d with the hypersphere σ 0 (see Exercise 14). Theorem 11. Let S be a set of n hyperspheres in R d . The worst-case complexity of a single Apollonius region in the diagram of n hyperspheres of R d is Θ(n d+1 2

  The bisector b ij subdivides R d into two open regions: one, b i ij , consisting of the points of R d that are closer to o i than to o j , and the other one, b j ij , consisting of the points of R d that are closer to o j than to o i . We can then define the Voronoi region of o i as the intersection of the regions b i ij for all j = i. The union of the closures of these Voronoi regions covers R d . Furthermore, if we assume that the bisectors are (d -1)-manifolds, the Voronoi regions then have disjoint interiors and we can define the closed region associated to b

  for a contradiction, that x ∈ b k ik . It follows that x ∈ bi ik , but x cannot lie on b ik , because this would imply that x ∈ b ik ∩b ij , which does not intersect b k jk . Hence, x ∈ b i ik and therefore, x ∈ b ij ∩ b k jk ∩ b i ik . We can then find x in the neighborhood of x such that x ∈ b j ij ∩ b k jk ∩ b i ki , contradicting Lemma 2. Let us now prove that b ij ∩ bk jk ⊂ bk ik . We have proved the inclusion for b ij ∩ b k jk . It remains to prove that b ij ∩ b jk ⊂ bk ik which is trivially true, by I.C. The two other inclusions are proved in a similar way.

Lemma 5 .

 5 For a given set B satisfying I.C. and assuming that we never have b ij ⊂ b ik for j = k, there are at most two ways of labeling the connected components of each R d \ b ij as b i ij and b j ij such that A.C. is verified. Proof. First assume that the sides b 1 12 and b 2 12 have been assigned. Consider now the labeling of the sides of b 1i , for some i > 2: let x be a point in the non empty set b 2i \ b 12 . First assume that x ∈ b 1 12 . Lemma 3 then implies that x ∈ b 1 1i . Conversely, if x ∈ b 2 12 , x ∈ b i 1i . In both cases, the assignment of the sides of b 1i is determined.

  e. ∃(λ, µ) ∈ R 2 β ki = λβ ij + µβ jk (L.C.C.) Note that L.C.C. implies I.C. and that in the case of affine bisectors L.C.C. is equivalent to I.C. Furthermore, it should be noted that, in the case of Voronoi diagrams, the bisector-functions defined as β ij = δ iδ j obviously satisfy L.C.C. We now prove that we can view diagrams satisfying L.C.C. as diagrams that can be linearized. Definition 1. A diagram D of n objects in some space E is said to be a pullback of a diagram D of m objects in space F by a function φ : E → F if m = n and if, for any distinct i, j, we have b

Theorem 16 .

 16 Any abstract diagram of R d with quadratic bisectors such that the corresponding degree 2 polynomials satisfy L.C.C. is an anisotropic Voronoi diagram. Exercise 23. Explain why, in Theorem 15, it is important to specify which bisector-functions satisfy L.C.C. instead of mentioning only the bisectors (Hint: Theorem 12 implies that there always exist some bisector-functions with the same zero-sets that satisfy L.C.C.)

Theorem 17 .

 17 The Euclidean Voronoi diagram of n points in R d and the power diagram of n spheres in R d can be constructed by an incremental algorithm in randomized time O n log n + n d+1 2

Fig. 1. 10 .

 10 Fig. 1.10. For the proof of Lemma 7.

Fig. 1. 11 .

 11 Fig. 1.11. Incremental construction of the Voronoi diagram of disjoint line segments.

  hierarchy. The tree attached to the region V (o) of o in Vor(Θ ) includes, for each Voronoi vertex v of V (o i ), the ray [v o , v) issued from the projection of v onto the boundary ∂o of o that passes through v. The rays are sorted by directions. Similarly, we associate to the query point x the ray [x o , x) where x o is the projection of x onto ∂o. When V (o) is visited, ray [x o , x) is located in the associated tree and we get the two rays [v o , v) and [w o , w) immediately before and after [x o , x). Let o be the neighbor of o whose cell is incident to v and w. We compare the distances from x to o and o . If o is closer to x than o, the walk steps to o . Otherwise, we know that o is the object of Θ closest to x and the walk halts (see [29] for details). Hence, visiting the Voronoi region of o i in Vor(Θ ) reduces to querying the tree and comparing the distances from x to o and o which takes O(log n ) time. Lemma 9. Using a hierarchy of Voronoi diagrams, nearest neighbor queries can be answered in expected time O(log 2 n). It has been shown [29] that the expected cost of updating the Voronoi hierarchy when inserting an object is O(log 2 n).

  et al. assumptions, the transitivity relation of Lemma 4 is satisfied and that Voronoi regions are simply connected. Exercise 29. Let O be a set of planar convex objects that may intersect and may not form a pseudo-circle set. Show that the Voronoi diagram Vor(O) may exhibit disconnected Voronoi regions. Propose an extension of the incremental algorithm to build the restriction Vor(O) ∩ U c of the Voronoi diagram Vor(O) to the region U c which is the complement of the union of the objects. The solution can be found in [29]. Exercise 30. Describes the geometric predicates required to implement the incremental algorithm. Provide algebraic expressions for the case of circles and line segments [31, 22].

Fig. 1. 13 .

 13 Fig. 1.13. On the left side, a closed curve S and the medial axis of Ω = R 2 \ S. On the right side, a dense sample E of points and its Voronoi diagram. The medial axis of R 2 \ E (which is the 1-skeleton of Vor(E)) is very different from the medial axis of Ω.

Lemma 10 .= x-y 4 ≥

 104 Let σ be a S-empty sphere centered at c, of radius r, intersecting S in two points x and y. If ε < r 2 and l def εr(1 -ε r ), there exists an S -empty sphere tangent to S in two points whose center c and radius r are such that |1 -r r | and c -c r are at most δ = εr l 2 -εr+ε 2 . Proof. Let σ be the maximal S -empty sphere centered at c and let r be its radius (seeFig. 1.14). The Hausdorff distance between S and S being at most ε, we have |rr | ≤ ε. Let y be a point of σ ∩ S .

Fig. 1 . 4 x -y 2 = r 2 +

 1422 Fig. 1.14. S is the continuous curve. x and y belong to S

4 √ D 3 ε and k (ε) = k 4 √ 2 √ 4 √

 4424 D 3 ε, such that M k(ε) (S) ⊂ M k (ε) (P) ⊕ B Dε ⊂ M k (ε) (S) ⊕ B Dε .Here, B r denotes the ball centered at the origin of radius r, and ⊕ the Minkowski sum.

ε→0d 4 √

 4 H (M(S), M 10 9D 3 ε (P ε )) = 0.

4 √

 4 Input: A surface S and a positive real ε 1. Sample S so as to obtain a sample P such that dH (S, P) ≤ ε; 2. Construct the Voronoi diagram of P; 3. Remove from the diagram the faces for which the diameter of the set of contact points is smaller than 10 9D 3 ε.Output: A PL approximation of M(S)

  is shown in Fig. 1.15. Exercise 31. Let O be a bounded open set. Show that M(O) is a retract of O (and therefore has the same homotopy type as O) [39].

Fig. 1 .

 1 Fig. 1.15. Two λ-medial axes of the same shape, with λ increasing from left to right, computed as a subset of the Voronoi diagram of a sample of the boundary (courtesy of Steve Oudot).

Fig. 1 .

 1 Fig. 1.16. A cell in an Apollonius diagram of spheres.

We may translate the hyperplanes vertically in order to have a non-empty intersection, or we may consider imaginary spheres with negative squared radii.

Acknowledgments

We thank D. Attali, C. Delage and M. Karavelas with whom part of the research reported in this chapter has been conducted. We also thank F. Chazal and A. Lieutier for fruitful discussions on the approximation of the medial axis.