Convex Hull and Voronoi Diagram of Additively ... - Semantic Scholar
Convex Hull and Voronoi Diagram of Additively Weighted Points
Jean-Daniel Boissonnat
Christophe Delage
April 12, 2005
Abstract We provide a complete description of dynamic algorithms for constructing convex hulls and Voronoi diagrams of additively weighted points of Rd . We present simple algorithms and provide a complete description of all the predicates. The algorithms have been implemented in R3 and experimental results are reported. Our implementation follows the CGAL design and, in particular, is made both robust and efficient through the use of filtered exact arithmetic.
1
Introduction
In this paper, we provide a complete description of dynamic algorithms for constructing convex hulls and Voronoi diagrams of additively weighted points of Rd . The algorithms have been implemented in R3 and experimental results are reported. Our motivation comes from the fact that weighted points can be considered as hyperspheres when the weights are positive and is twofold. On one hand, spheres are non linear objects and, besides the combinatorial and algorithmic questions, numerical and robustness issues deserve a careful investigation, which has not been fully done yet. On the other hand, spheres are objects of major concern in various fields, most notably structural biology, and effective implementations of basic geometric algorithms for spheres are needed. We first revisit the problem of computing the convex hull of n weighted points of Rd . This problem has already been solved optimally [5, 6]. We present a simpler fully dynamic algorithm and provide a complete description of all the predicates for any d. We then consider the construction of additively weighted Voronoi diagrams. It is known that the construction of such diagrams reduces to intersecting a power diagram in Rd 1 with half-cones [3]. We are not aware of robust implementations of this algorithm. Other algorithms have been recently designed and implemented in the planar case [8, 13]. In R3 , we are only aware of two prototype implementations, one by Will [14] for computing a single cell, and one by Kim et al. [12] that computes the entire diagram. None of these implementations is provably robust. Moreover, the algorithm by Kim et al. assumes that the graph of the edges of each cell is connected, which is not true in general. We apply our result on the construction of the convex hull of additively weighted points to the construction of a Voronoi cell in the Voronoi diagram of n additively weighted points. The construction, which makes use of inversion, is close to the algorithm of [6]. The main contribution of this work is to provide a full analysis of the predicates involved, a thorough treatment of the degenerate cases, and a CGAL implementation. Our predicates, when specialized to the planar case (d 2), are simpler and of lower degree than the best predicates known so far
[email protected] [email protected]
1
[11, 2]. Our implementation follows the CGAL design and, in particular, is made both robust and efficient through the use of filtered exact arithmetic. The paper is organized as follows. In section 2, we establish a new correspondence between convex hulls of additively weighted points in Rd and power diagrams of spheres of Rd , from which we deduce an algorithm to construct such hulls. In section 3, we recall a similar correspondence for a cell in the Voronoi diagram of additively weighted points and present an algorithm for constructing such a cell. In section 4, we describe the predicates. In section 5, we show how to handle the degenerate cases. In section 6, we report on experimental results. Finally, we conclude in section 7. In the sequel, a weighted point, or site for short, of Rd is a pair p, w where p is a point of d R and w is a real number. When w is positive, we also call a weighted point a hypersphere. Given a set n spheres Σ σ1 , . . . , σn of Rd , σi ci , ri , the power of a point x Rd to σi is:
d σ , x x c r and the power diagram of Σ, noted P Σ , is the subdivision of R consisting of the n cells P σ , . . . , P σ where P σ x R , d σ , x d σ , x , j 1, . . . , n . P σ We write P σ , . . . , σ . . . P σ . When non empty, P σ , . . . , σ is a face of P Σ , of dimension d k if the spheres are in general position.
2 i
2
P
i
i
d
1
n
i
0
2
k
d
i
P
0
P
j
0
k
k
Convex Hull of Additively Weighted Points
s0 , . . . , sn be a set of weighted points of Rd . We write si pi , wi , i 0, . . . , n. We Let S consider first the case where all the weights wi are non negative, i.e. the sites are spheres. The convex hull of S, CH S , is the smallest closed convex subset of Rd containing all the spheres of S. A supporting hyperplane H of S is a hyperplane tangent to at least one of the spheres of S, and such that all the spheres of S lie in the same half-space limited by H. A face of CH S of circularity k, 0 k d, is the portion of CH S that consists of the points whose supporting hyperplanes are tangent to a same subset of d k spheres. For d 3, faces of circularity 0 are planar faces tangent to three spheres, faces of circularity 1 are conical patches tangent to two spheres, and faces of circularity 2 are spherical patches contained in some si .
2.1
From power diagrams to convex hulls of spheres
Let Π be a supporting hyperplane tangent to k spheres, and m be the unit normal vector of Π that is oriented away from S. As there is exactly one supporting hyperplane that has a given oriented normal, m defines Π exactly. Π is a supporting hyperplane tangent to s1 , . . . , sk if and only if (cf Fig. 1)
p p m p p
m
i
1
i
1
w1 w1
wi , 1 wi , k
i i
k
(1)
n
(2)
We rewrite 2 as follows:
m
p p w w m p p m p w 2m p p p 2p m 2m p p p p 2w m p p 2w , m p 1
i
1
m2
2 1
1
1
2
1
2 1
2 1
1
1
1
i
i
i
i
2
i
2
2
2 i
2 i
2 i
i
2wi
Lemma 1 The k-faces of P σ , . . . , σ S are in 1-1 correspondence with the faces of circularity k 1 of CH S .
1 rewrites the same way. Thus, denoting ri2 p2i 2wi , σi pi , ri and Σ σ0 , . . . , σ k , this is equivalent to m being in the open d k 1 -face of σ1 , . . . , σk in the power diagram P Σ . As m belongs to the unit sphere S x R: x 1 , we have proven
1
n
Π r1 m
ri
c1
ci
p p
Figure 1: Here, m
i
1
w1
wi .
The above construction works also when some or all wi are negative. Although a geometric interpretation in terms of convex hull is then missing, the result of our construction is called the convex hull of the weighted points, or AWCH for short, by analogy to the case of positive weights.
2.2
Algorithm
According to lemma 1, constructing the convex hull of the set S of additively weighted points (or σ0 , . . . , σ n sites) s0 , . . . , sn reduces to computing the intersection of the power diagram of Σ with the unit sphere S. The affine hull of a face f is denoted aff f . We say that a k-face f is a sub-face of a k 1 -face f (or equivalently that f is a super-face of f ) when f f . For a face f and a sub-face f of f , H f, f denotes the halfspace of aff f bounded by aff f that contains f . For instance, when f is a 1-face (a line segment), f is one of its endpoints, and H f, f is a ray. We present a static algorithm and an incremental algorithm. We assume that the si are in general position so that we do not have any degeneracy. Degeneracies will be considered in section 5.
Static algorithm. The algorithm first constructs the power diagram P Σ and then determines, for each face f of P Σ , whether f intersects S or not. The result is stored in tag f :
• tag f
if and only if aff f is outside S,
• tag f
if and only if f does not intersect S but aff f intersects S,
• tag f
if and only if f intersects S,
• tag f
if and only if f is inside S.
Assuming we know tag f for each sub-face f of f , we compute tag f as follows:
1. If aff f does not intersect S, then tag f ,
2. else, if for each sub-face f of f , tag f , then, by convexity of f , tag f ,
3
3. else, if there is a sub-face f of f such that tag f of f , tag f .
or tag f , then, by connexity
4. else, if for each sub-face f of f , tag f , and aff f interior of f intersects S and tag f ,
S
H f, f , then the relative
5. else, tag f .
Assume wlog that f P σ0 , . . . , σk 1 and f the two following geometric predicates.
P σ0 , . . . , σk . This algorithm needs to evaluate
k-RadicalIntersection f determines whether aff f is outside S or not.
k-RadicalSide f, f determines whether aff f intersects S and aff f is outside S.
S
H f, f , assuming that aff f
Computing all the tags clearly takes time proportional to the size of the diagram. It follows that d our algorithm computes the additively weighted convex hull of n sites in time O n log n n 2 , which is worst-case optimal.
Incremental algorithm. We use an incremental algorithm for constructing the power diagram, and we attach to each face f the number num f of its sub-faces that are tagged . The insertion of a new sphere in the power diagram creates some new faces, deletes some faces, and replaces some faces by smaller ones. In this last case, a face f is replaced by a smaller face f¯ f that is incident to the cell of the new sphere. f¯ is called a cut face. Notice that a cut face of dimension less than d has exactly one new subface. We denote by m the number of deleted faces plus the number of new faces. 1. For a new face f , num f is easily computed by looking at all its sub-faces. This can be done in time proportional to m. 2. We now update num f¯ for a cut face f¯. We set num f¯ num f . Then num f¯ is decremented by the number of the sub-faces of f that are deleted, and updated according to the tag of the new and cut sub-faces of f¯. This can be simply done in the following way. When the tag of a face f changes, or f is deleted, we update num f for each super -face f of f . Updating num f¯ therefore takes O m time.
3. We update tag f¯ for a cut k-face f¯, k 1. We compute tag f¯ from num f¯ and tag f , where f is the new sub-face f of f¯. As f¯ is a cut face, tag f¯ differs from tag f only when tag f
. If num f¯ is positive, then tag f¯ . Now, if num f¯ is 0, only the relative interior of f¯ can intersect S. Hence
• if tag f , then we update tag f¯ according to the outcome of k-RadicalSide f¯, f • otherwise, as the set of the sub-faces of f¯ is connected (f¯ is of dimension at least 2), f¯ has the same tag as f .
For a cut 1-face f¯, we just compute tag f¯ directly. Updating the tag of a cut face takes a constant time. The incremental additively weighted convex hull algorithm has therefore the same complexity as the incremental power diagram. When the sites are inserted in random order, the expected d time complexity of our algorithm is therefore O n log n n 2 .
4
Under some assumptions, this bound can be improved. First, according to our experiments (see section 6), the number of spheres with a non empty cell in the power diagram is usually proportional to the number of points h on the additively weighted convex hull. In that case, the running time for n insertions is:
O n log h
d
h 2 ,
Moreover, h can be much smaller than n. It is known that the convex hull of a set of n points uniformly distributed inside a sphere of R3 has O n points on its convex hull. The same result holds trivially for spheres with the same radius. In R3 , assuming that the number of cells in the power diagram is proportional to h and that h is O n , the complexity of our algorithm is O n log n .
3
Additively Weighted Voronoi Diagram
The additively weighted distance, denoted d , from a point m of Rd to a site si d
si , m
pi
m
p ,w i
i
is:
wi .
s , . . . , s of sites, the additively weighted Voronoi cell of s , V s is: m R j, d s , m d s , m . V s It is possible that V s ∅, but this only happens when m R , d s , m d s , m , for Considering the set S
0
i
i
n
i
d
i
j
i
d
j
i
some sj . In that case, we say that si is hidden by sj . When dealing with spheres (i.e. ri , rj 0), si is hidden by sj when si sj . The additively weighted Voronoi diagram, or AWVD for short, of S, noted V S , is the cell complex whose d-cells are the V si . The construction of a single cell of the diagram reduces to computing an additively weighted convex hull, via an inversion. More precisely, the cell of s0 in V S is combinatorially equivalent to the additively weighted convex hull of S s0 , s1 , . . . , sn , where s0 is centered at the origin pi p0 2 wi w0 2 . This and has weight 0, si is centered at piαip0 and has weight piαip0 , αi scheme only works when s0 is not hidden by, nor does it hide any other si . See [6] for details. Using the construction of a single cell as a subroutine, we can compute the whole diagram in a fully dynamic manner. All the Voronoi cells are stored in a data structure with pointers between the corresponding elements. In this data structure, Γ s denotes the set of neighbors of s and Γs s the set of neighbors of s that are also neighbors of s . In order to keep track of the hidden sites, we keep, for each site s, a list hidden s of the sites s hides. We now describe the three main ingredients needed to update an additively weighted Voronoi diagram.
Localization. Given a point m of Rd and a starting site s, this procedure, called Locate m, s , returns the cell of the diagram in which m lies. This is done by means of a simple walk: • if a neighbor s is closer to m than s, jump to s and iterate, • otherwise, m belongs to the cell of s. This walk is linear in the size of the diagram, but can be improved using a hierarchical data structure as in [8], or using an auxiliary location data structure like a Kd-tree. This localization algorithm requires only one predicate: given two sites s0 and s1 , determine if a query point m is closer to s0 than to s1 . This predicate is called SideOfBisector. 5
Insertion. The insertion procedure needs to decide whether a site hides another site : we call this predicate IsTrivial. To avoid ambiguities when considering diagrams of different sets of Algorithm 1 Inserts a site in an additively weighted Voronoi diagram Insert s : Locate c, s0 , where c is the center of s 1: sc 2: if IsTrivial s, sc then 3: hidden sc hidden sc s 4: else 5: queue sc ; mark sc 6: while queue ∅ do 7: let si queue; queue queue si 8: if IsTrivial si , s then 9: hidden s hidden s si 10: R Γ si 11: else 12: N Γ si 13: insert si in V s and s in V si 14: R N Γ si Γsi s 15: end if 16: R r R r is not marked 17: queue queue R; mark R 18: end while 19: unmark hidden s Γ s 20: end if
sites, we introduce the following notation. Given a set of sites S and s, s S, we denote VS s the cell of s in the additively weighted Voronoi diagram of S, and VS s, s VS s VS s . Definition 1 A site s
S is hidden in S if and only if VS s
Definition 2 A site s S is in conflict with s
∅.
S if and only if VS s
V s s . S
Notice that only the non-hidden sites of S may be in conflict with s.
Definition 3 The conflict graph of a site s S is G X, E where X ∅. in conflict with s, and xy E if and only if VS x, y V S s s
S is the set of sites
In other words, the conflict graph of s is the dual of the portion of the cells and d the Voronoi diagram that will disappear after when inserting s. Lemma 2 The conflict graph is connected.
1 -faces of
Proof. Given two sites x and y in the conflict graph G of some s, we take px in VS x VS s s and py in VS y VS s s . As VS s s is arc connected, we can follow a path from px to py in VS s s , and each time we cross a d 1 -face of the diagram, we follow the corresponding edge of G. This gives a path from x to y in G. The first insertion (i.e. the insertion in an empty diagram) is easy, we just create a new cell, insert in it, and insert the new site in V . Once there is a least one site in the diagram, we use algorithm 3 which essentially performs a breadth-first walk on the conflict graph, handling the newly hidden sites along the way. We are now ready to prove the following.
6
∅ and a new site s S, the procedure
Theorem 1 Given the AWVD of a set of sites S Insert s computes the AWVD of S s . Proof. When inserting s, there are two cases.
1. s will be hidden. This is correctly detected by the test of line 2. Indeed, s is hidden if and only if it is hidden by sc (which is the closest site to s), see [9], Lemma 1. 2. s is not hidden. Assuming that si is in conflict with s, then R (just before line 16) is the set of neighbors of si in the conflict graph of s. Also, sc is in conflict with s. Thus, lines 5 to 18 perform a breadth-first walk on the conflict graph G of s. As G is connected (lemma 2), all of its vertices are visited. For each si , if it is hidden by s, it is added to hidden s (line 9), and si is not inserted in the Voronoi cells of s. Removing a site s from a diagram is straightforward:
Removal.
• remove s from the cells adjacent to V
s,
• let s1 , . . . , sk be the neighbors of s; for all 1 to rebuild the hole made by the removal of s,
i, j
k, i
j, insert s into the cell of s i j
• insert all the sites s was hiding.
Complexity. An insertion in the Voronoi diagram of n sites performs O n insertions in O n additively weighted convex hulls, plus a localization, which can be implemented in O n time. The overall time to construct the Voronoi diagram of n sites is therefore:
n log n
O n n Which gives, for d
d
n 2
O n2 log n
d
n 2
1
.
3, O n3 . This bound can be improved under the following assumptions:
• O n sites appear on the convex hull,
• the sites have O 1 neighbors,
• the underlying power diagram of every Voronoi cell has a O s non-hidden points, where s is the number of neighbors of the cell. Those assumptions are not too restrictive, they happen to be satisfied on a variety of input data (see section 6.) Under those assumptions, each insertion in the diagram involves: one localization, O 1 insertions in Voronoi cells, and O 1 insertions in the convex hull. Using a hierarchy as in [8] gives an average localization running time of O log n , under the current assumptions; O n insertions in the convex hull cost O n log n ; and O n insertions in Voronoi cells of O 1 size cost O n . Hence, under the above assumptions, the incremental construction of the AWVD of n sites can be done in O n log n time. This analysis is corroborated by our experiments (see section 6.)
4
Predicates
We consider a set S s0 , . . . , sn of sites, where si is centered at pi , and has weight wi . If each input data is a b-bit integer, the size of each monomial occuring in a predicate is upper bounded by 2 b 1 d . Moreover, let v be the number of variables that occur in a predicate; for
7
the predicates considered in this paper, v is a constant. It follows that a predicate of degree d requires precision p d b 1 log v . Here, the predicates are polynomials in the unknowns pi , wi , and the algebraic degree of each of them is given. In addition to the predicates mentioned in section 2.2 and 3, we need the two well-known predicates that are needed to construct the power diagram: Orientation and PowerTest. Predicates IsTrivial and SideOfBisector are detailed in [10] for d 2, and are straightforward to extend to arbitary dimension. IsTrivial is of degree 2, and SideOfBisector is of degree 4. Basic linear algebra provides explicit formulas for the other predicates (see the appendice.) A summary of the degree of all the predicates is given below. Notice that k d is a special case for the k-RadicalIntersection predicate, leading to a lower degree than the k d case.
Algorithm Dimension IsTrivial SideOfBisector Orientation PowerTest 1-RadicalIntersection 2-RadicalIntersection 3-RadicalIntersection k-RadicalIntersection, 1 d-RadicalIntersection 1-RadicalSide 2-RadicalSide 3-RadicalSide k-RadicalSide, 1 k d Maximum degree
AWCH 3 d 3
2
2 3 2 4 k
3 4 2 8 6
d d
2 8 12 4k 2d 1 3 5 2k 1 4d 4
d
1 3
1
1 3 5
AWVD 3 d 3 2 2 4 4 5 d 2 6 d 3 6 6 16 16 10 20 4k 8 2d 4 3 3 7 7 9 9 2k 3 16 4d 4
2 2 4 4 5 6 8
3 7
4
8
8
The predicates of the algorithm for the additively weighted Voronoi diagram of [8], detailed in [11], have a maximal degree of 16 for d 2. Here we have degree 8 for d 2, and 16 for d 3.
5
Degenerate Cases
5.1
Additively Weighted Convex Hull
In this section, we show how to handle degeneracies in the algorithm for the convex hull of additively weighted points. We call a case degenerate when some predicate returns 0, instead of “positive” or “negative”. A simple way of dealing with these cases is to carefully choose a non-zero sign to be returned by a predicate when it evaluates to 0. • Predicates Orientation or PowerTest return zero when Σ is not in general position. Any standard perturbation scheme will work for us.
• When predicate k-RadicalIntersection returns 0, some d k -flat is tangent to S (it intersects S but not the open ball bounded by S.) We can consider that this d k -flat lies outside S.
• Predicate k-RadicalSide f, f returns 0 when the projection of the origin on aff f is on aff f . In that case, both aff f and aff f intersects S, or both do not intersect S. As predicate k-RadicalSide is only called when aff f intersects S and aff f does not, this predicate is never called on degenerate inputs.
8
5.2
Additively Weighted Voronoi Diagram
In the case of the additively weighted Voronoi diagram, the previous perturbation scheme does not work. Indeed, as we compute the Voronoi cells separately, we not only need to resolve degeneracies in each cell, but also to ensure that consistent decisions are taken when we compute the neighboring cells.
Definition 4 A set of sites S is called degenerate when there exists k 1 sites of S s0 , . . . , sk , 1 k d, such that V s0 , . . . , sk is not empty and is of dimension strictly less than d k. When an input is non-degenerate, any small enough perturbation will let the combinatorial structure of the Voronoi diagram unchanged. For a face f , we define L f s S : f V s . If a degeneracy s0 , . . . , sk is minimal (i.e. if s0 , . . . , sk si is not degenerate for 0 i k) then perturbing the weight of any site in L V s0 , . . . , sk will remove the degeneracy. Given a degenerate input s0 , . . . , sk , finding a minimal degeneracy is easy : w.l.o.g. we ∅, V s , . . . , s check if s0 , . . . , sk 1 is still degenerate. As V s0 , . . . , sk 0 k 1 has dimension at least 1. A degeneracy of dimension m 1 in the intersection between the unit sphere and a power diagram can only appear if two m 1 -faces of the power diagram are equal, which can be tested with the Orientation and PowerTest predicates. Now, we can handle the degeneracies in our algorithm by means of symbolic perturbations. When faced with a predicate that returns zero on s0 , . . . , sk , we find a minimal degeneracy s0 , . . . , sk and perturb one site, in L V s0 , . . . , sk , say si , i.e. we replace wi by wi . The predicates are now polynomials in and we need to evaluate the sign of the non-zero coefficient of smallest degree. Notice that this scheme does not increase the degree of the predicates since the perturbation is linear. It can occur that some predicate returns zero, and the Voronoi diagram is not degenerate, and thus some site gets perturbed unnecessarily. To ensure consistency from one cell to another, we just need to choose the site to perturb in a way that is independent of the site whose cell is under construction when we detect the degeneracy. One way to do that is to choose the smallest site according to some global ordering (for instance, the lexicographical order on the centers.)
6
Experimental results
A prototype implementation of the convex hull of weighted points, using the algorithm presented here, has been done in R3 . It uses CGAL 3.1, mainly its 3D regular triangulations (see [7]), which are the duals of the power diagrams. While not yet fully optimized, this implementation already follows the CGAL standard of genericity and robustness. It is parametrizable through the use of C++ templates and can dynamically filter the predicates to avoid problems in degenerate (or near-degenerate) cases. At the time of writing, the 3D regular triangulation is not yet compatible with the hierarchical localization structure of CGAL, thus, localization in the regular triangulation is done by a simple O k walk (where k is the size of the triangulation.) Thus, further speed-up of our implementation is expected when a fast localization in regular triangulations will be available. The running time of Figure 2 are obtained on a Athlon at 1333MHz, with 133MHz DDRSDRAM memory and 256Ko of L2 cache. The degenerate input is a set of sites randomly chosen in a cube, with their weight equal to their height, so that all the sites are tangent to the lower face of the cube. We have also implemented our algorithm to construct additively weighted Voronoi diagram in R3 . Running times (on the same machine) are shown in Figure 3. Benchmarks were done
9
18
18 increasing weight decreasing weight random
14
increasing weight decreasing weight random
16 Running time in seconds
Running time in seconds
16
12 10 8 6 4
14 12 10 8 6 4
2
2
0
0 1
2
3 4 5 6 7 8 9 Number of input spheres (in thousands)
10
0
5 10 15 20 25 30 Number of input spheres (in thousands)
degenerate input
35
non-degenerate input
Figure 2: Additively Weighted convex hull benchmarks, all using filtered predicates, for various input, and insertion order. using dynamically filtered predicates. On this non-degenerate input, there is, roughly, a factor 5 slow-down compared to double-precision arithmetic. The input here comes from a direct application of our algorithm. The sites have their centers on a surface (here, the hand model) and the weights are of the form lfskx where lfs x is an approximation of the local feature size of the surface at x, and k is a parameter. This kind of diagram has been used to efficiently compute a sizing field, which is useful for 3D meshing (see [1] for details). On Figure 3, k 1 on the left and k 0.3 on the right.
1000
1200
Running time in seconds
800
Running time in seconds
increasing weight decreasing weight random
900
700 600 500 400 300 200
increasing weight decreasing weight random
1000 800 600 400 200
100 0
0 0
5 10 15 20 25 30 Number of input spheres (in thousands)
35
0
no hidden site
5 10 15 20 25 30 Number of input spheres (in thousands)
35
8,369 non-hidden sites
Figure 3: Additively weighted Voronoi diagram benchmarks using filtered predicates for various input sizes, numbers of hidden sites, and insertion orders. Like many other incremental algorithms, this one is sensitive to the insertion order. Figure 3 shows clearly that it is best to insert sites in the order of decreasing weights. The reason is that a site with a larger weight tends to have more neighbors, and thus, tends to take longer to insert. The difference is even greater when there are many hidden spheres. A screenshot is shown in Figure 4, where one cell is represented by meshing its boundary. Our implementation computes all the edges of the cell (they are faces of circularity 1 in the underlying convex hull), and sample them. Then, for each face of the cell (ie. each face of circularity 2 in the convex hull) is approximated using the meshing algorithm of [4]. We plan to have the code included in the cgal library soon. 10
Figure 4: A screenshot of one cell of a 3D additively weighted Voronoi diagram. The size of power diagrams. When running the AWCH algorithm on random spheres, with centers uniformly distributed in a ball, and radii uniformly chosen in some interval, we observe a remarkable phenomenon: almost all the spheres that do not contribute to the convex hull are hidden (i.e. have empty cells) in the underlying power diagram. This also occur when the AWCH are cells of an additively weighted Voronoi diagram. This phenomenon has been observed on many input sets: in no Voronoi cells of this example, the number of cells in the power diagram is more than seven times the number of neighbors of the Voronoi cell. Moreover, for only 1% of the Voronoi cells, the number of cells in the underlying power diagram is more than twice the number of neighbors in the Voronoi diagram. Although this observation does not hold in general —it is possible to construct n spheres such that only O 1 of them contribute to the convex hull, while all of them appear in the power diagram— this makes our algorithms efficient in practice.
7
Conclusion
We have presented fully robust implementations of two algorithms for constructing the convex hull and the Voronoi diagram of additively weighted points (and spheres). To the best of our knowledge, no certified algorithms existed previously. This work does not settle the main open question in this area : what is the combinatorial complexity of the Voronoi diagram of n additively weighted points in Rd ? Tight bounds are only known for d 2 and odd dimensions. However, experimenting with our code may provide new insights such as the one mentionned in section 6 that eventually will help improving the combinatorial bounds.
References [1] P. Alliez, D. Cohen-Steiner, M. Yvinec, and M. Desbrun. Variational tetrahedral meshing, 2005. To appear at SIGGRAPH. [2] Francois Anton. Columbia, 2004.
Voronoi diagrams of semi-algebraic sets.
Ph.d. thesis, University of British
[3] F. Aurenhammer and H. Imai. Geometric relations among Voronoi diagrams. Geom. Dedicata, 27:65–75, 1988.
11
[4] J. D. Boissonnat and S. Oudot. Provably good surface sampling and approximation. In Proc. 1st Symp. on Geometry Processing, pages 9–18, 2003. [5] Jean-Daniel Boissonnat, Andr´e C´er´ezo, Olivier Devillers, Jacqueline Duquesne, and Mariette Yvinec. An algorithm for constructing the convex hull of a set of spheres in dimension d. Comput. Geom. Theory Appl., 6:123–130, 1996. [6] Jean-Daniel Boissonnat and Menelaos Karavelas. On the combinatorial complexity of Euclidean Voronoi cells and convex hulls of d-dimensional spheres. In Proc. 14th ACM-SIAM Sympos. Discrete Algorithms (SODA), pages 305–312, 2003. [7] The CGAL Manual, 2004. Release 3.1. [8] Menelaos Karavelas and Mariette Yvinec. Dynamic additively weighted voronoi diagrams in 2d. In Proc. 10th European Symposium on Algorithms, pages 586–598, 2002. [9] Menelaos Karavelas and Mariette Yvinec. Dynamic additively weighted voronoi diagrams in 2d. Rapport de recherche 4466, INRIA, 2002. [10] Menelaos I. Karavelas and Ioannis Z. Emiris. Predicates for the planar additively weighted Voronoi diagram. Technical Report ECG-TR-122201-01, INRIA Sophia-Antipolis, 2002. [11] Menelaos I. Karavelas and Ioannis Z. Emiris. Root comparison techniques applied to computing the additively weighted Voronoi diagram. In Proc. 14th ACM-SIAM Sympos. Discrete Algorithms (SODA), pages 320–329, 2003. [12] Deok-Soo Kim, Youngsong Cho, Donguk Kim, Jonghwa Bhak, and Sung-Hoon Lee. Euclidean voronoi diagram of 3d spheres and applications to protein structure analysis. In Kokichi Sugihara, editor, 1st International Symposium on Voronoi Diagrams in Science and Engineering, 2004. [13] Deok-Soo Kim, Donguk Kim, and Kokichi Sugihara. Updating the topology of the dynamic voronoi diagram for spheres in euclidean d-dimensional space. Computer-Aided Design, 18:541–562, 2001. [14] Hans-Martin Will. Fast and efficient computation of additively weighted Voronoi cells for applications in molecular biology. In Proc. 6th Scand. Workshop Algorithm Theory, volume 1432 of Lecture Notes Comput. Sci., pages 310–321. Springer-Verlag, 1998.
12
A
Predicates
In order to avoid doing the calculi twice, we give formulas that are valid for the predicates of both the convex hull and the Voronoi diagram algorithms. To do so, we consider sites centered at αqii , with weights αωii . In the case of the convex hull, we take: qi
pi ,
ωi
wi ,
αi
1
and, if we want to compute the cell of sj in the AWVD, we take : qi
pi
pj ,
ωi
wi
wj ,
αi
qi2 1
j
ωi2 when i when i
j
qi 2 c2i 2 αωii , σi ci , ri . In addition, We denote dα the degree of We denote ci αi and ri an αi (dα is 0 in the case of a convex hull, and 2 for the additively weighted Voronoi diagram.) Also, note that if si and sj do not hide each other, then pi pj 2 wi wj 2 0 (see [6]), thus, we always have αi 0. In the following, we will be primarily interested in the sign of expressions rather than in sgn their values. For the sake of convenience, we will write A B, when A and B have the same sign.
A.1
The Orientation Predicate
1 -tuple of points c , . . . , c 1 1 1 1 c
c
This predicate is therefore of degree d d .
The orientation of the d
0
d
is:
q0 α0
0
d
qd αd
α0
sgn
q0
αd qd
.
α
A.2
The PowerTest Predicate
The power test of the d
c2 0
1 c0
r02
1 cd
This is a degree d
A.3
c2d
2 -tuple of weighted points σ , . . . , σ , σ is: 1
1 1 α 1
q c ω c r 2
2 2 0
t
rd2 dα
2 t
qd αd ωd αd
q0 α0 ω0 α0
2 t
d
qt αt ωt αt
t
0
sgn
0
0
αd αt qd qt ωd ωt
1 predicate.
The k-RadicalIntersection Predicate
Given a face of a power diagram f , defined by s0 , . . . , sk , this predicate determines if aff f is entierly outside the unit sphere. We want to know where the origin projects onto aff f , comparing the norm of this projection π to 1 will give the k-RadicalIntersection predicate. π is found by resolving the linear system Ak X bk , keeping only the coordinates corresponding to m in X, where:
0 0 .. .
Ak
0 α0 .. .
0
0
α0
α
k
Id
Qtk
Qk
0kk
1 1
,
X
u m λ0 .. . λk
αk 13
0 .. .
,
bk
0 ω0 .. .
ωk
where 0kk is the zero k k matrix, Id is the identity d d matrix, and Qk q0 , . . . , qk is a d k 1 matrix. The details are in appendix A.5. The matrix Ak where the ith column is i replaced by bk is denoted Ak . Then, the ith coordinate of the projection πk of O onto aff f is:
i 1
Ak
i πk
Ak
The degree of Ak can be obtained by developing the determinant along the first line and first column, giving a sum of term of the form:
Id
M2 0
αi αj M1
,
where M1 and M2 are Qtk and Qk , minus one line and one column, respectively. This term rewrites to αi αj M1 M2 . M1 M2 is a k k matrix, whose coefficient are of degree 2, therefore,
Ak is of degree 2 k dα . Clearly, Ak has the same degree as Ak . For k d, the projection is easier to compute. As aff f degenerates to a single point, the projection of O is this same point. When k d, the projection π of O is the solution of the system Ad X bd , where:
i 1
Ad
X
α0 .. . Qtd , αd
u m
,
bd
ω0 .. . . ωd
Ad is then a d 1 d 1 matrix, and has degree (as well as Ad k 1 is also a special case :
i 1
π1
α 0 ω1 α 0 q1
α 1 ω0 α 0 q1 α 1 q0 2
) d
dα .
α 1 q0 .
The k-RadicalIntersection predicate is equivalent to finding the sign of:
1
2
d Aki 1 Ak 2 i 1
sgn
Ak
2
d i Ak i 1
d . When k
When k d, this polynomial has degree 4 k When k 1, we have:
1 2
α
.
d, it has degree 2 d
d . α
1
π12
This is a degree 2 dα
A.4
1
1
α 0 ω1 α 0 q1
α 1 ω0 2 α 1 q0 2
sgn
α 0 q1
α 1 q0
α ω 2
0 1
α 1 ω0 2 .
predicate.
The k-RadicalSide Predicate
We are given a face of a power diagram f and f a sub-face of f . W.l.o.g. f and f are defined by σ0 , . . . , σk 1 and σ0 , . . . , σk , respectively. Assuming that aff f intersects the unit sphere, and that aff f does not, this predicate determines if H f, f intersects the unit sphere. We
14
want to know on which side of P σ0 , σk the orthogonal projection onto aff f of O lies. Using the notations of the previous section, we test the sign of:
where xi
j
Developing the 2 Ak
1
Ak
1
α0
Ak
i Ak 1
Ak 1
α k 1
i
α0 x0 i i Ak α k xk
1
0
i
Ak 1
0 α0 .. .
α0
1
d i 1
d i xk i 1 0 0 .. .
k
α0 αk
α0
ω0 ωk
α k
α0 αk
ω0 A ωk k
1
2 determinants gives:
which can be rewritten to:
1
i x0 i xk
is the j th coordinate of qi . This expressions has the same sign as:
d i 1
0
ω k Ak
Id
0 0 .. .
Qtk 1
0 ω0 .. .
α
d i x0 i 1
αk
Qk
i
αk
0kk
0
0
0
Qtk
α k 1
ω 0 Ak
Id
0 α0 .. .
1
Ak 1
0 0 .. .
k 1
ωk 1 ωk 0
qkt
0
α0
1
1
q0t
0
1
0 0 .. .
α0
α Qk
0 ω0 .. .
.
k 1
1
0kk
ωk 1 ω0 0
0
Now, taking the term of α0 from the second determinant and putting it in the first one gives:
0 0 .. .
1
k
Ak
1
0 α0 .. .
α0
αk 1
αk
0
0
0 0 .. .
Id
Qtk
α0
Qk
0 0 .. .
k 1
0 0 α1 .. .
0
αk
0kk
1
ωk 1 ωk 0
qkt
0 ω0 .. .
1
α
0 Id q0t q1t .. .
α k 1
0
0
qkt 1 q0t
0 0 .. .
α0
0 ω0 ω1 .. . ωk 1 ω0 0
α
0 0 .. .
1
k
Ak
1
0 α0 .. .
α k
0
0 Id
0 0 .. .
Qtk
0 ω0 .. . ωk
15
α0
α
k 1
Qk
1
0kk
1
Qk
1
.
0kk
0
There are two identical lines in the second determinant: this determinant is therefore null. As α0 0, we finally get:
k 1
Using the same technique used for computing the degree of Ak 1 , we find that this predicate is degree 2 k dα 1. When k 1, this predicate degenerates to finding on which side of P σ0 , σ1 the origin lies. This is given by the sign of: α0 ω1 α 1 ω0
which is a degree dα
A.5
1 polynomial.
Projection onto a k-Flat of the Origin
In this section, we compute the projection of O onto a k-flat F defined by σ0 , . . . , σk . F is defined by the following equations:
c r . m c r .. . m c r m c r c 2 , we can rewrite the preceeding system of equation to: α q α q m αω αω .. . . α q α q m αω αω m
2
0
pi αi
and ri2
1
2 0
2
0
In our case, ci
2 0
2
2 1
2
2 k
k
ωi αi
2 i
0 1
1 0
0 k
k 0
0 1
1 0
0 k
k 0
Now, the orthogonal projection of O onto F lies on F and on the hyperplane generated by the vectors α0 q1 α1 q0 , . . . , α0 qk αk q0 , thus it is the solution of:
λ1 α 0 q 1
α 1 q0
α 0 qk Putting u
q0 m ω 0 α0
and λ0
1 α0
m m αq m
k i 1 α i qi ,
k 0
u, m α ,q k
α 0 ω1
α 0 ωk
.. .
m αω m α 1 ω0
.
k 0
we can rewrite the previous system as follows: 0 0
α 0 λ0 α k λk λ0 q0 λk qk m α0 , q 0 u, m
0
λk α 0 qk α k q0 α 0 q1 α 1 q0
k
16
ω0
ωk
.. .
Jean-Daniel Boissonnat
Christophe Delage
April 12, 2005
Abstract We provide a complete description of dynamic algorithms for constructing convex hulls and Voronoi diagrams of additively weighted points of Rd . We present simple algorithms and provide a complete description of all the predicates. The algorithms have been implemented in R3 and experimental results are reported. Our implementation follows the CGAL design and, in particular, is made both robust and efficient through the use of filtered exact arithmetic.
1
Introduction
In this paper, we provide a complete description of dynamic algorithms for constructing convex hulls and Voronoi diagrams of additively weighted points of Rd . The algorithms have been implemented in R3 and experimental results are reported. Our motivation comes from the fact that weighted points can be considered as hyperspheres when the weights are positive and is twofold. On one hand, spheres are non linear objects and, besides the combinatorial and algorithmic questions, numerical and robustness issues deserve a careful investigation, which has not been fully done yet. On the other hand, spheres are objects of major concern in various fields, most notably structural biology, and effective implementations of basic geometric algorithms for spheres are needed. We first revisit the problem of computing the convex hull of n weighted points of Rd . This problem has already been solved optimally [5, 6]. We present a simpler fully dynamic algorithm and provide a complete description of all the predicates for any d. We then consider the construction of additively weighted Voronoi diagrams. It is known that the construction of such diagrams reduces to intersecting a power diagram in Rd 1 with half-cones [3]. We are not aware of robust implementations of this algorithm. Other algorithms have been recently designed and implemented in the planar case [8, 13]. In R3 , we are only aware of two prototype implementations, one by Will [14] for computing a single cell, and one by Kim et al. [12] that computes the entire diagram. None of these implementations is provably robust. Moreover, the algorithm by Kim et al. assumes that the graph of the edges of each cell is connected, which is not true in general. We apply our result on the construction of the convex hull of additively weighted points to the construction of a Voronoi cell in the Voronoi diagram of n additively weighted points. The construction, which makes use of inversion, is close to the algorithm of [6]. The main contribution of this work is to provide a full analysis of the predicates involved, a thorough treatment of the degenerate cases, and a CGAL implementation. Our predicates, when specialized to the planar case (d 2), are simpler and of lower degree than the best predicates known so far
[email protected] [email protected]
1
[11, 2]. Our implementation follows the CGAL design and, in particular, is made both robust and efficient through the use of filtered exact arithmetic. The paper is organized as follows. In section 2, we establish a new correspondence between convex hulls of additively weighted points in Rd and power diagrams of spheres of Rd , from which we deduce an algorithm to construct such hulls. In section 3, we recall a similar correspondence for a cell in the Voronoi diagram of additively weighted points and present an algorithm for constructing such a cell. In section 4, we describe the predicates. In section 5, we show how to handle the degenerate cases. In section 6, we report on experimental results. Finally, we conclude in section 7. In the sequel, a weighted point, or site for short, of Rd is a pair p, w where p is a point of d R and w is a real number. When w is positive, we also call a weighted point a hypersphere. Given a set n spheres Σ σ1 , . . . , σn of Rd , σi ci , ri , the power of a point x Rd to σi is:
d σ , x x c r and the power diagram of Σ, noted P Σ , is the subdivision of R consisting of the n cells P σ , . . . , P σ where P σ x R , d σ , x d σ , x , j 1, . . . , n . P σ We write P σ , . . . , σ . . . P σ . When non empty, P σ , . . . , σ is a face of P Σ , of dimension d k if the spheres are in general position.
2 i
2
P
i
i
d
1
n
i
0
2
k
d
i
P
0
P
j
0
k
k
Convex Hull of Additively Weighted Points
s0 , . . . , sn be a set of weighted points of Rd . We write si pi , wi , i 0, . . . , n. We Let S consider first the case where all the weights wi are non negative, i.e. the sites are spheres. The convex hull of S, CH S , is the smallest closed convex subset of Rd containing all the spheres of S. A supporting hyperplane H of S is a hyperplane tangent to at least one of the spheres of S, and such that all the spheres of S lie in the same half-space limited by H. A face of CH S of circularity k, 0 k d, is the portion of CH S that consists of the points whose supporting hyperplanes are tangent to a same subset of d k spheres. For d 3, faces of circularity 0 are planar faces tangent to three spheres, faces of circularity 1 are conical patches tangent to two spheres, and faces of circularity 2 are spherical patches contained in some si .
2.1
From power diagrams to convex hulls of spheres
Let Π be a supporting hyperplane tangent to k spheres, and m be the unit normal vector of Π that is oriented away from S. As there is exactly one supporting hyperplane that has a given oriented normal, m defines Π exactly. Π is a supporting hyperplane tangent to s1 , . . . , sk if and only if (cf Fig. 1)
p p m p p
m
i
1
i
1
w1 w1
wi , 1 wi , k
i i
k
(1)
n
(2)
We rewrite 2 as follows:
m
p p w w m p p m p w 2m p p p 2p m 2m p p p p 2w m p p 2w , m p 1
i
1
m2
2 1
1
1
2
1
2 1
2 1
1
1
1
i
i
i
i
2
i
2
2
2 i
2 i
2 i
i
2wi
Lemma 1 The k-faces of P σ , . . . , σ S are in 1-1 correspondence with the faces of circularity k 1 of CH S .
1 rewrites the same way. Thus, denoting ri2 p2i 2wi , σi pi , ri and Σ σ0 , . . . , σ k , this is equivalent to m being in the open d k 1 -face of σ1 , . . . , σk in the power diagram P Σ . As m belongs to the unit sphere S x R: x 1 , we have proven
1
n
Π r1 m
ri
c1
ci
p p
Figure 1: Here, m
i
1
w1
wi .
The above construction works also when some or all wi are negative. Although a geometric interpretation in terms of convex hull is then missing, the result of our construction is called the convex hull of the weighted points, or AWCH for short, by analogy to the case of positive weights.
2.2
Algorithm
According to lemma 1, constructing the convex hull of the set S of additively weighted points (or σ0 , . . . , σ n sites) s0 , . . . , sn reduces to computing the intersection of the power diagram of Σ with the unit sphere S. The affine hull of a face f is denoted aff f . We say that a k-face f is a sub-face of a k 1 -face f (or equivalently that f is a super-face of f ) when f f . For a face f and a sub-face f of f , H f, f denotes the halfspace of aff f bounded by aff f that contains f . For instance, when f is a 1-face (a line segment), f is one of its endpoints, and H f, f is a ray. We present a static algorithm and an incremental algorithm. We assume that the si are in general position so that we do not have any degeneracy. Degeneracies will be considered in section 5.
Static algorithm. The algorithm first constructs the power diagram P Σ and then determines, for each face f of P Σ , whether f intersects S or not. The result is stored in tag f :
• tag f
if and only if aff f is outside S,
• tag f
if and only if f does not intersect S but aff f intersects S,
• tag f
if and only if f intersects S,
• tag f
if and only if f is inside S.
Assuming we know tag f for each sub-face f of f , we compute tag f as follows:
1. If aff f does not intersect S, then tag f ,
2. else, if for each sub-face f of f , tag f , then, by convexity of f , tag f ,
3
3. else, if there is a sub-face f of f such that tag f of f , tag f .
or tag f , then, by connexity
4. else, if for each sub-face f of f , tag f , and aff f interior of f intersects S and tag f ,
S
H f, f , then the relative
5. else, tag f .
Assume wlog that f P σ0 , . . . , σk 1 and f the two following geometric predicates.
P σ0 , . . . , σk . This algorithm needs to evaluate
k-RadicalIntersection f determines whether aff f is outside S or not.
k-RadicalSide f, f determines whether aff f intersects S and aff f is outside S.
S
H f, f , assuming that aff f
Computing all the tags clearly takes time proportional to the size of the diagram. It follows that d our algorithm computes the additively weighted convex hull of n sites in time O n log n n 2 , which is worst-case optimal.
Incremental algorithm. We use an incremental algorithm for constructing the power diagram, and we attach to each face f the number num f of its sub-faces that are tagged . The insertion of a new sphere in the power diagram creates some new faces, deletes some faces, and replaces some faces by smaller ones. In this last case, a face f is replaced by a smaller face f¯ f that is incident to the cell of the new sphere. f¯ is called a cut face. Notice that a cut face of dimension less than d has exactly one new subface. We denote by m the number of deleted faces plus the number of new faces. 1. For a new face f , num f is easily computed by looking at all its sub-faces. This can be done in time proportional to m. 2. We now update num f¯ for a cut face f¯. We set num f¯ num f . Then num f¯ is decremented by the number of the sub-faces of f that are deleted, and updated according to the tag of the new and cut sub-faces of f¯. This can be simply done in the following way. When the tag of a face f changes, or f is deleted, we update num f for each super -face f of f . Updating num f¯ therefore takes O m time.
3. We update tag f¯ for a cut k-face f¯, k 1. We compute tag f¯ from num f¯ and tag f , where f is the new sub-face f of f¯. As f¯ is a cut face, tag f¯ differs from tag f only when tag f
. If num f¯ is positive, then tag f¯ . Now, if num f¯ is 0, only the relative interior of f¯ can intersect S. Hence
• if tag f , then we update tag f¯ according to the outcome of k-RadicalSide f¯, f • otherwise, as the set of the sub-faces of f¯ is connected (f¯ is of dimension at least 2), f¯ has the same tag as f .
For a cut 1-face f¯, we just compute tag f¯ directly. Updating the tag of a cut face takes a constant time. The incremental additively weighted convex hull algorithm has therefore the same complexity as the incremental power diagram. When the sites are inserted in random order, the expected d time complexity of our algorithm is therefore O n log n n 2 .
4
Under some assumptions, this bound can be improved. First, according to our experiments (see section 6), the number of spheres with a non empty cell in the power diagram is usually proportional to the number of points h on the additively weighted convex hull. In that case, the running time for n insertions is:
O n log h
d
h 2 ,
Moreover, h can be much smaller than n. It is known that the convex hull of a set of n points uniformly distributed inside a sphere of R3 has O n points on its convex hull. The same result holds trivially for spheres with the same radius. In R3 , assuming that the number of cells in the power diagram is proportional to h and that h is O n , the complexity of our algorithm is O n log n .
3
Additively Weighted Voronoi Diagram
The additively weighted distance, denoted d , from a point m of Rd to a site si d
si , m
pi
m
p ,w i
i
is:
wi .
s , . . . , s of sites, the additively weighted Voronoi cell of s , V s is: m R j, d s , m d s , m . V s It is possible that V s ∅, but this only happens when m R , d s , m d s , m , for Considering the set S
0
i
i
n
i
d
i
j
i
d
j
i
some sj . In that case, we say that si is hidden by sj . When dealing with spheres (i.e. ri , rj 0), si is hidden by sj when si sj . The additively weighted Voronoi diagram, or AWVD for short, of S, noted V S , is the cell complex whose d-cells are the V si . The construction of a single cell of the diagram reduces to computing an additively weighted convex hull, via an inversion. More precisely, the cell of s0 in V S is combinatorially equivalent to the additively weighted convex hull of S s0 , s1 , . . . , sn , where s0 is centered at the origin pi p0 2 wi w0 2 . This and has weight 0, si is centered at piαip0 and has weight piαip0 , αi scheme only works when s0 is not hidden by, nor does it hide any other si . See [6] for details. Using the construction of a single cell as a subroutine, we can compute the whole diagram in a fully dynamic manner. All the Voronoi cells are stored in a data structure with pointers between the corresponding elements. In this data structure, Γ s denotes the set of neighbors of s and Γs s the set of neighbors of s that are also neighbors of s . In order to keep track of the hidden sites, we keep, for each site s, a list hidden s of the sites s hides. We now describe the three main ingredients needed to update an additively weighted Voronoi diagram.
Localization. Given a point m of Rd and a starting site s, this procedure, called Locate m, s , returns the cell of the diagram in which m lies. This is done by means of a simple walk: • if a neighbor s is closer to m than s, jump to s and iterate, • otherwise, m belongs to the cell of s. This walk is linear in the size of the diagram, but can be improved using a hierarchical data structure as in [8], or using an auxiliary location data structure like a Kd-tree. This localization algorithm requires only one predicate: given two sites s0 and s1 , determine if a query point m is closer to s0 than to s1 . This predicate is called SideOfBisector. 5
Insertion. The insertion procedure needs to decide whether a site hides another site : we call this predicate IsTrivial. To avoid ambiguities when considering diagrams of different sets of Algorithm 1 Inserts a site in an additively weighted Voronoi diagram Insert s : Locate c, s0 , where c is the center of s 1: sc 2: if IsTrivial s, sc then 3: hidden sc hidden sc s 4: else 5: queue sc ; mark sc 6: while queue ∅ do 7: let si queue; queue queue si 8: if IsTrivial si , s then 9: hidden s hidden s si 10: R Γ si 11: else 12: N Γ si 13: insert si in V s and s in V si 14: R N Γ si Γsi s 15: end if 16: R r R r is not marked 17: queue queue R; mark R 18: end while 19: unmark hidden s Γ s 20: end if
sites, we introduce the following notation. Given a set of sites S and s, s S, we denote VS s the cell of s in the additively weighted Voronoi diagram of S, and VS s, s VS s VS s . Definition 1 A site s
S is hidden in S if and only if VS s
Definition 2 A site s S is in conflict with s
∅.
S if and only if VS s
V s s . S
Notice that only the non-hidden sites of S may be in conflict with s.
Definition 3 The conflict graph of a site s S is G X, E where X ∅. in conflict with s, and xy E if and only if VS x, y V S s s
S is the set of sites
In other words, the conflict graph of s is the dual of the portion of the cells and d the Voronoi diagram that will disappear after when inserting s. Lemma 2 The conflict graph is connected.
1 -faces of
Proof. Given two sites x and y in the conflict graph G of some s, we take px in VS x VS s s and py in VS y VS s s . As VS s s is arc connected, we can follow a path from px to py in VS s s , and each time we cross a d 1 -face of the diagram, we follow the corresponding edge of G. This gives a path from x to y in G. The first insertion (i.e. the insertion in an empty diagram) is easy, we just create a new cell, insert in it, and insert the new site in V . Once there is a least one site in the diagram, we use algorithm 3 which essentially performs a breadth-first walk on the conflict graph, handling the newly hidden sites along the way. We are now ready to prove the following.
6
∅ and a new site s S, the procedure
Theorem 1 Given the AWVD of a set of sites S Insert s computes the AWVD of S s . Proof. When inserting s, there are two cases.
1. s will be hidden. This is correctly detected by the test of line 2. Indeed, s is hidden if and only if it is hidden by sc (which is the closest site to s), see [9], Lemma 1. 2. s is not hidden. Assuming that si is in conflict with s, then R (just before line 16) is the set of neighbors of si in the conflict graph of s. Also, sc is in conflict with s. Thus, lines 5 to 18 perform a breadth-first walk on the conflict graph G of s. As G is connected (lemma 2), all of its vertices are visited. For each si , if it is hidden by s, it is added to hidden s (line 9), and si is not inserted in the Voronoi cells of s. Removing a site s from a diagram is straightforward:
Removal.
• remove s from the cells adjacent to V
s,
• let s1 , . . . , sk be the neighbors of s; for all 1 to rebuild the hole made by the removal of s,
i, j
k, i
j, insert s into the cell of s i j
• insert all the sites s was hiding.
Complexity. An insertion in the Voronoi diagram of n sites performs O n insertions in O n additively weighted convex hulls, plus a localization, which can be implemented in O n time. The overall time to construct the Voronoi diagram of n sites is therefore:
n log n
O n n Which gives, for d
d
n 2
O n2 log n
d
n 2
1
.
3, O n3 . This bound can be improved under the following assumptions:
• O n sites appear on the convex hull,
• the sites have O 1 neighbors,
• the underlying power diagram of every Voronoi cell has a O s non-hidden points, where s is the number of neighbors of the cell. Those assumptions are not too restrictive, they happen to be satisfied on a variety of input data (see section 6.) Under those assumptions, each insertion in the diagram involves: one localization, O 1 insertions in Voronoi cells, and O 1 insertions in the convex hull. Using a hierarchy as in [8] gives an average localization running time of O log n , under the current assumptions; O n insertions in the convex hull cost O n log n ; and O n insertions in Voronoi cells of O 1 size cost O n . Hence, under the above assumptions, the incremental construction of the AWVD of n sites can be done in O n log n time. This analysis is corroborated by our experiments (see section 6.)
4
Predicates
We consider a set S s0 , . . . , sn of sites, where si is centered at pi , and has weight wi . If each input data is a b-bit integer, the size of each monomial occuring in a predicate is upper bounded by 2 b 1 d . Moreover, let v be the number of variables that occur in a predicate; for
7
the predicates considered in this paper, v is a constant. It follows that a predicate of degree d requires precision p d b 1 log v . Here, the predicates are polynomials in the unknowns pi , wi , and the algebraic degree of each of them is given. In addition to the predicates mentioned in section 2.2 and 3, we need the two well-known predicates that are needed to construct the power diagram: Orientation and PowerTest. Predicates IsTrivial and SideOfBisector are detailed in [10] for d 2, and are straightforward to extend to arbitary dimension. IsTrivial is of degree 2, and SideOfBisector is of degree 4. Basic linear algebra provides explicit formulas for the other predicates (see the appendice.) A summary of the degree of all the predicates is given below. Notice that k d is a special case for the k-RadicalIntersection predicate, leading to a lower degree than the k d case.
Algorithm Dimension IsTrivial SideOfBisector Orientation PowerTest 1-RadicalIntersection 2-RadicalIntersection 3-RadicalIntersection k-RadicalIntersection, 1 d-RadicalIntersection 1-RadicalSide 2-RadicalSide 3-RadicalSide k-RadicalSide, 1 k d Maximum degree
AWCH 3 d 3
2
2 3 2 4 k
3 4 2 8 6
d d
2 8 12 4k 2d 1 3 5 2k 1 4d 4
d
1 3
1
1 3 5
AWVD 3 d 3 2 2 4 4 5 d 2 6 d 3 6 6 16 16 10 20 4k 8 2d 4 3 3 7 7 9 9 2k 3 16 4d 4
2 2 4 4 5 6 8
3 7
4
8
8
The predicates of the algorithm for the additively weighted Voronoi diagram of [8], detailed in [11], have a maximal degree of 16 for d 2. Here we have degree 8 for d 2, and 16 for d 3.
5
Degenerate Cases
5.1
Additively Weighted Convex Hull
In this section, we show how to handle degeneracies in the algorithm for the convex hull of additively weighted points. We call a case degenerate when some predicate returns 0, instead of “positive” or “negative”. A simple way of dealing with these cases is to carefully choose a non-zero sign to be returned by a predicate when it evaluates to 0. • Predicates Orientation or PowerTest return zero when Σ is not in general position. Any standard perturbation scheme will work for us.
• When predicate k-RadicalIntersection returns 0, some d k -flat is tangent to S (it intersects S but not the open ball bounded by S.) We can consider that this d k -flat lies outside S.
• Predicate k-RadicalSide f, f returns 0 when the projection of the origin on aff f is on aff f . In that case, both aff f and aff f intersects S, or both do not intersect S. As predicate k-RadicalSide is only called when aff f intersects S and aff f does not, this predicate is never called on degenerate inputs.
8
5.2
Additively Weighted Voronoi Diagram
In the case of the additively weighted Voronoi diagram, the previous perturbation scheme does not work. Indeed, as we compute the Voronoi cells separately, we not only need to resolve degeneracies in each cell, but also to ensure that consistent decisions are taken when we compute the neighboring cells.
Definition 4 A set of sites S is called degenerate when there exists k 1 sites of S s0 , . . . , sk , 1 k d, such that V s0 , . . . , sk is not empty and is of dimension strictly less than d k. When an input is non-degenerate, any small enough perturbation will let the combinatorial structure of the Voronoi diagram unchanged. For a face f , we define L f s S : f V s . If a degeneracy s0 , . . . , sk is minimal (i.e. if s0 , . . . , sk si is not degenerate for 0 i k) then perturbing the weight of any site in L V s0 , . . . , sk will remove the degeneracy. Given a degenerate input s0 , . . . , sk , finding a minimal degeneracy is easy : w.l.o.g. we ∅, V s , . . . , s check if s0 , . . . , sk 1 is still degenerate. As V s0 , . . . , sk 0 k 1 has dimension at least 1. A degeneracy of dimension m 1 in the intersection between the unit sphere and a power diagram can only appear if two m 1 -faces of the power diagram are equal, which can be tested with the Orientation and PowerTest predicates. Now, we can handle the degeneracies in our algorithm by means of symbolic perturbations. When faced with a predicate that returns zero on s0 , . . . , sk , we find a minimal degeneracy s0 , . . . , sk and perturb one site, in L V s0 , . . . , sk , say si , i.e. we replace wi by wi . The predicates are now polynomials in and we need to evaluate the sign of the non-zero coefficient of smallest degree. Notice that this scheme does not increase the degree of the predicates since the perturbation is linear. It can occur that some predicate returns zero, and the Voronoi diagram is not degenerate, and thus some site gets perturbed unnecessarily. To ensure consistency from one cell to another, we just need to choose the site to perturb in a way that is independent of the site whose cell is under construction when we detect the degeneracy. One way to do that is to choose the smallest site according to some global ordering (for instance, the lexicographical order on the centers.)
6
Experimental results
A prototype implementation of the convex hull of weighted points, using the algorithm presented here, has been done in R3 . It uses CGAL 3.1, mainly its 3D regular triangulations (see [7]), which are the duals of the power diagrams. While not yet fully optimized, this implementation already follows the CGAL standard of genericity and robustness. It is parametrizable through the use of C++ templates and can dynamically filter the predicates to avoid problems in degenerate (or near-degenerate) cases. At the time of writing, the 3D regular triangulation is not yet compatible with the hierarchical localization structure of CGAL, thus, localization in the regular triangulation is done by a simple O k walk (where k is the size of the triangulation.) Thus, further speed-up of our implementation is expected when a fast localization in regular triangulations will be available. The running time of Figure 2 are obtained on a Athlon at 1333MHz, with 133MHz DDRSDRAM memory and 256Ko of L2 cache. The degenerate input is a set of sites randomly chosen in a cube, with their weight equal to their height, so that all the sites are tangent to the lower face of the cube. We have also implemented our algorithm to construct additively weighted Voronoi diagram in R3 . Running times (on the same machine) are shown in Figure 3. Benchmarks were done
9
18
18 increasing weight decreasing weight random
14
increasing weight decreasing weight random
16 Running time in seconds
Running time in seconds
16
12 10 8 6 4
14 12 10 8 6 4
2
2
0
0 1
2
3 4 5 6 7 8 9 Number of input spheres (in thousands)
10
0
5 10 15 20 25 30 Number of input spheres (in thousands)
degenerate input
35
non-degenerate input
Figure 2: Additively Weighted convex hull benchmarks, all using filtered predicates, for various input, and insertion order. using dynamically filtered predicates. On this non-degenerate input, there is, roughly, a factor 5 slow-down compared to double-precision arithmetic. The input here comes from a direct application of our algorithm. The sites have their centers on a surface (here, the hand model) and the weights are of the form lfskx where lfs x is an approximation of the local feature size of the surface at x, and k is a parameter. This kind of diagram has been used to efficiently compute a sizing field, which is useful for 3D meshing (see [1] for details). On Figure 3, k 1 on the left and k 0.3 on the right.
1000
1200
Running time in seconds
800
Running time in seconds
increasing weight decreasing weight random
900
700 600 500 400 300 200
increasing weight decreasing weight random
1000 800 600 400 200
100 0
0 0
5 10 15 20 25 30 Number of input spheres (in thousands)
35
0
no hidden site
5 10 15 20 25 30 Number of input spheres (in thousands)
35
8,369 non-hidden sites
Figure 3: Additively weighted Voronoi diagram benchmarks using filtered predicates for various input sizes, numbers of hidden sites, and insertion orders. Like many other incremental algorithms, this one is sensitive to the insertion order. Figure 3 shows clearly that it is best to insert sites in the order of decreasing weights. The reason is that a site with a larger weight tends to have more neighbors, and thus, tends to take longer to insert. The difference is even greater when there are many hidden spheres. A screenshot is shown in Figure 4, where one cell is represented by meshing its boundary. Our implementation computes all the edges of the cell (they are faces of circularity 1 in the underlying convex hull), and sample them. Then, for each face of the cell (ie. each face of circularity 2 in the convex hull) is approximated using the meshing algorithm of [4]. We plan to have the code included in the cgal library soon. 10
Figure 4: A screenshot of one cell of a 3D additively weighted Voronoi diagram. The size of power diagrams. When running the AWCH algorithm on random spheres, with centers uniformly distributed in a ball, and radii uniformly chosen in some interval, we observe a remarkable phenomenon: almost all the spheres that do not contribute to the convex hull are hidden (i.e. have empty cells) in the underlying power diagram. This also occur when the AWCH are cells of an additively weighted Voronoi diagram. This phenomenon has been observed on many input sets: in no Voronoi cells of this example, the number of cells in the power diagram is more than seven times the number of neighbors of the Voronoi cell. Moreover, for only 1% of the Voronoi cells, the number of cells in the underlying power diagram is more than twice the number of neighbors in the Voronoi diagram. Although this observation does not hold in general —it is possible to construct n spheres such that only O 1 of them contribute to the convex hull, while all of them appear in the power diagram— this makes our algorithms efficient in practice.
7
Conclusion
We have presented fully robust implementations of two algorithms for constructing the convex hull and the Voronoi diagram of additively weighted points (and spheres). To the best of our knowledge, no certified algorithms existed previously. This work does not settle the main open question in this area : what is the combinatorial complexity of the Voronoi diagram of n additively weighted points in Rd ? Tight bounds are only known for d 2 and odd dimensions. However, experimenting with our code may provide new insights such as the one mentionned in section 6 that eventually will help improving the combinatorial bounds.
References [1] P. Alliez, D. Cohen-Steiner, M. Yvinec, and M. Desbrun. Variational tetrahedral meshing, 2005. To appear at SIGGRAPH. [2] Francois Anton. Columbia, 2004.
Voronoi diagrams of semi-algebraic sets.
Ph.d. thesis, University of British
[3] F. Aurenhammer and H. Imai. Geometric relations among Voronoi diagrams. Geom. Dedicata, 27:65–75, 1988.
11
[4] J. D. Boissonnat and S. Oudot. Provably good surface sampling and approximation. In Proc. 1st Symp. on Geometry Processing, pages 9–18, 2003. [5] Jean-Daniel Boissonnat, Andr´e C´er´ezo, Olivier Devillers, Jacqueline Duquesne, and Mariette Yvinec. An algorithm for constructing the convex hull of a set of spheres in dimension d. Comput. Geom. Theory Appl., 6:123–130, 1996. [6] Jean-Daniel Boissonnat and Menelaos Karavelas. On the combinatorial complexity of Euclidean Voronoi cells and convex hulls of d-dimensional spheres. In Proc. 14th ACM-SIAM Sympos. Discrete Algorithms (SODA), pages 305–312, 2003. [7] The CGAL Manual, 2004. Release 3.1. [8] Menelaos Karavelas and Mariette Yvinec. Dynamic additively weighted voronoi diagrams in 2d. In Proc. 10th European Symposium on Algorithms, pages 586–598, 2002. [9] Menelaos Karavelas and Mariette Yvinec. Dynamic additively weighted voronoi diagrams in 2d. Rapport de recherche 4466, INRIA, 2002. [10] Menelaos I. Karavelas and Ioannis Z. Emiris. Predicates for the planar additively weighted Voronoi diagram. Technical Report ECG-TR-122201-01, INRIA Sophia-Antipolis, 2002. [11] Menelaos I. Karavelas and Ioannis Z. Emiris. Root comparison techniques applied to computing the additively weighted Voronoi diagram. In Proc. 14th ACM-SIAM Sympos. Discrete Algorithms (SODA), pages 320–329, 2003. [12] Deok-Soo Kim, Youngsong Cho, Donguk Kim, Jonghwa Bhak, and Sung-Hoon Lee. Euclidean voronoi diagram of 3d spheres and applications to protein structure analysis. In Kokichi Sugihara, editor, 1st International Symposium on Voronoi Diagrams in Science and Engineering, 2004. [13] Deok-Soo Kim, Donguk Kim, and Kokichi Sugihara. Updating the topology of the dynamic voronoi diagram for spheres in euclidean d-dimensional space. Computer-Aided Design, 18:541–562, 2001. [14] Hans-Martin Will. Fast and efficient computation of additively weighted Voronoi cells for applications in molecular biology. In Proc. 6th Scand. Workshop Algorithm Theory, volume 1432 of Lecture Notes Comput. Sci., pages 310–321. Springer-Verlag, 1998.
12
A
Predicates
In order to avoid doing the calculi twice, we give formulas that are valid for the predicates of both the convex hull and the Voronoi diagram algorithms. To do so, we consider sites centered at αqii , with weights αωii . In the case of the convex hull, we take: qi
pi ,
ωi
wi ,
αi
1
and, if we want to compute the cell of sj in the AWVD, we take : qi
pi
pj ,
ωi
wi
wj ,
αi
qi2 1
j
ωi2 when i when i
j
qi 2 c2i 2 αωii , σi ci , ri . In addition, We denote dα the degree of We denote ci αi and ri an αi (dα is 0 in the case of a convex hull, and 2 for the additively weighted Voronoi diagram.) Also, note that if si and sj do not hide each other, then pi pj 2 wi wj 2 0 (see [6]), thus, we always have αi 0. In the following, we will be primarily interested in the sign of expressions rather than in sgn their values. For the sake of convenience, we will write A B, when A and B have the same sign.
A.1
The Orientation Predicate
1 -tuple of points c , . . . , c 1 1 1 1 c
c
This predicate is therefore of degree d d .
The orientation of the d
0
d
is:
q0 α0
0
d
qd αd
α0
sgn
q0
αd qd
.
α
A.2
The PowerTest Predicate
The power test of the d
c2 0
1 c0
r02
1 cd
This is a degree d
A.3
c2d
2 -tuple of weighted points σ , . . . , σ , σ is: 1
1 1 α 1
q c ω c r 2
2 2 0
t
rd2 dα
2 t
qd αd ωd αd
q0 α0 ω0 α0
2 t
d
qt αt ωt αt
t
0
sgn
0
0
αd αt qd qt ωd ωt
1 predicate.
The k-RadicalIntersection Predicate
Given a face of a power diagram f , defined by s0 , . . . , sk , this predicate determines if aff f is entierly outside the unit sphere. We want to know where the origin projects onto aff f , comparing the norm of this projection π to 1 will give the k-RadicalIntersection predicate. π is found by resolving the linear system Ak X bk , keeping only the coordinates corresponding to m in X, where:
0 0 .. .
Ak
0 α0 .. .
0
0
α0
α
k
Id
Qtk
Qk
0kk
1 1
,
X
u m λ0 .. . λk
αk 13
0 .. .
,
bk
0 ω0 .. .
ωk
where 0kk is the zero k k matrix, Id is the identity d d matrix, and Qk q0 , . . . , qk is a d k 1 matrix. The details are in appendix A.5. The matrix Ak where the ith column is i replaced by bk is denoted Ak . Then, the ith coordinate of the projection πk of O onto aff f is:
i 1
Ak
i πk
Ak
The degree of Ak can be obtained by developing the determinant along the first line and first column, giving a sum of term of the form:
Id
M2 0
αi αj M1
,
where M1 and M2 are Qtk and Qk , minus one line and one column, respectively. This term rewrites to αi αj M1 M2 . M1 M2 is a k k matrix, whose coefficient are of degree 2, therefore,
Ak is of degree 2 k dα . Clearly, Ak has the same degree as Ak . For k d, the projection is easier to compute. As aff f degenerates to a single point, the projection of O is this same point. When k d, the projection π of O is the solution of the system Ad X bd , where:
i 1
Ad
X
α0 .. . Qtd , αd
u m
,
bd
ω0 .. . . ωd
Ad is then a d 1 d 1 matrix, and has degree (as well as Ad k 1 is also a special case :
i 1
π1
α 0 ω1 α 0 q1
α 1 ω0 α 0 q1 α 1 q0 2
) d
dα .
α 1 q0 .
The k-RadicalIntersection predicate is equivalent to finding the sign of:
1
2
d Aki 1 Ak 2 i 1
sgn
Ak
2
d i Ak i 1
d . When k
When k d, this polynomial has degree 4 k When k 1, we have:
1 2
α
.
d, it has degree 2 d
d . α
1
π12
This is a degree 2 dα
A.4
1
1
α 0 ω1 α 0 q1
α 1 ω0 2 α 1 q0 2
sgn
α 0 q1
α 1 q0
α ω 2
0 1
α 1 ω0 2 .
predicate.
The k-RadicalSide Predicate
We are given a face of a power diagram f and f a sub-face of f . W.l.o.g. f and f are defined by σ0 , . . . , σk 1 and σ0 , . . . , σk , respectively. Assuming that aff f intersects the unit sphere, and that aff f does not, this predicate determines if H f, f intersects the unit sphere. We
14
want to know on which side of P σ0 , σk the orthogonal projection onto aff f of O lies. Using the notations of the previous section, we test the sign of:
where xi
j
Developing the 2 Ak
1
Ak
1
α0
Ak
i Ak 1
Ak 1
α k 1
i
α0 x0 i i Ak α k xk
1
0
i
Ak 1
0 α0 .. .
α0
1
d i 1
d i xk i 1 0 0 .. .
k
α0 αk
α0
ω0 ωk
α k
α0 αk
ω0 A ωk k
1
2 determinants gives:
which can be rewritten to:
1
i x0 i xk
is the j th coordinate of qi . This expressions has the same sign as:
d i 1
0
ω k Ak
Id
0 0 .. .
Qtk 1
0 ω0 .. .
α
d i x0 i 1
αk
Qk
i
αk
0kk
0
0
0
Qtk
α k 1
ω 0 Ak
Id
0 α0 .. .
1
Ak 1
0 0 .. .
k 1
ωk 1 ωk 0
qkt
0
α0
1
1
q0t
0
1
0 0 .. .
α0
α Qk
0 ω0 .. .
.
k 1
1
0kk
ωk 1 ω0 0
0
Now, taking the term of α0 from the second determinant and putting it in the first one gives:
0 0 .. .
1
k
Ak
1
0 α0 .. .
α0
αk 1
αk
0
0
0 0 .. .
Id
Qtk
α0
Qk
0 0 .. .
k 1
0 0 α1 .. .
0
αk
0kk
1
ωk 1 ωk 0
qkt
0 ω0 .. .
1
α
0 Id q0t q1t .. .
α k 1
0
0
qkt 1 q0t
0 0 .. .
α0
0 ω0 ω1 .. . ωk 1 ω0 0
α
0 0 .. .
1
k
Ak
1
0 α0 .. .
α k
0
0 Id
0 0 .. .
Qtk
0 ω0 .. . ωk
15
α0
α
k 1
Qk
1
0kk
1
Qk
1
.
0kk
0
There are two identical lines in the second determinant: this determinant is therefore null. As α0 0, we finally get:
k 1
Using the same technique used for computing the degree of Ak 1 , we find that this predicate is degree 2 k dα 1. When k 1, this predicate degenerates to finding on which side of P σ0 , σ1 the origin lies. This is given by the sign of: α0 ω1 α 1 ω0
which is a degree dα
A.5
1 polynomial.
Projection onto a k-Flat of the Origin
In this section, we compute the projection of O onto a k-flat F defined by σ0 , . . . , σk . F is defined by the following equations:
c r . m c r .. . m c r m c r c 2 , we can rewrite the preceeding system of equation to: α q α q m αω αω .. . . α q α q m αω αω m
2
0
pi αi
and ri2
1
2 0
2
0
In our case, ci
2 0
2
2 1
2
2 k
k
ωi αi
2 i
0 1
1 0
0 k
k 0
0 1
1 0
0 k
k 0
Now, the orthogonal projection of O onto F lies on F and on the hyperplane generated by the vectors α0 q1 α1 q0 , . . . , α0 qk αk q0 , thus it is the solution of:
λ1 α 0 q 1
α 1 q0
α 0 qk Putting u
q0 m ω 0 α0
and λ0
1 α0
m m αq m
k i 1 α i qi ,
k 0
u, m α ,q k
α 0 ω1
α 0 ωk
.. .
m αω m α 1 ω0
.
k 0
we can rewrite the previous system as follows: 0 0
α 0 λ0 α k λk λ0 q0 λk qk m α0 , q 0 u, m
0
λk α 0 qk α k q0 α 0 q1 α 1 q0
k
16
ω0
ωk
.. .
