Euclidean Voronoi diagram of 3D balls and its ... - CiteSeerX
Computer-Aided Design 37 (2005) 1412–1424 www.elsevier.com/locate/cad
Euclidean Voronoi diagram of 3D balls and its computation via tracing edges Deok-Soo Kima,b,*, Youngsong Chob, Donguk Kimb a
Department of Industrial Engineering, Hanyang University, 17 Haengdang-dong, Seongdong-gu, Seoul 133-791, South Korea b Voronoi Diagram Research Center, Hanyang University, 17 Haengdang-dong, Seongdong-gu, Seoul 133-791, South Korea Received 20 November 2004; accepted 5 February 2005
Abstract Despite its important applications in various disciplines in science and engineering, the Euclidean Voronoi diagram for spheres, also known as an additively weighted Voronoi diagram, in 3D space has not been studied as much as it deserves. In this paper, we present an algorithm to compute the Euclidean Voronoi diagram for 3D spheres with different radii. The presented algorithm follows Voronoi edges one by one until the construction is completed in O(mn) time in the worst-case, where m is the number of edges in the Voronoi diagram and n is the number of spherical balls. As building blocks, we show that Voronoi edges are conics that can be precisely represented as rational quadratic Be´zier curves. We also discuss how to conveniently represent and process Voronoi faces which are hyperboloids of two sheets. q 2005 Elsevier Ltd. All rights reserved. Keywords: Euclidean Voronoi diagrams; Edge-tracing; Rational quadratic Be´zier; Apollonius problem
1. Introduction Since its initial definition and characterization by the Russian mathematician G. Voronoı¨ in 1907 and the first computational algorithm presented by Shamos and Hoey [30], the Voronoi diagram has been known as one of the central topics in many disciplines in science and engineering including computational geometry. Due to its natural descriptive and manipulative capability, the Voronoi diagram and its variations have been known as various names such as a Thiessen polygon [34], a medial axis transformation (MAT) [19], a symmetric axis transformation (SAT) [4,5], a skeleton [18], a proximity map [12], a Dirichlit tessellation, a thinning, etc. Even long time before Voronoı¨, an idea quite close to what is now known as a Voronoi diagram was used by Descartes to * Corresponding author. Address: Department of Industrial Engineering, Hanyang University, 17 Haengdang-dong, Seongdong-gu, Seoul 133-791, South Korea. Tel.: C82 2 2220 0472; fax: C82 2 2292 0472. E-mail addresses: [email protected] (D.-S. Kim), [email protected] (Y. Cho), [email protected] (D. Kim).
0010-4485//$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.cad.2005.02.013
describe the distribution of matters in the solar system and its environs [24]. We believe that Aurenhammer [3] and Okabe et al. [24] might be the best source of information for the comprehensive introduction to various Voronoi diagrams, their properties, algorithms and related problems. The ordinary Voronoi diagram for point set has been studied extensively and its properties are well-known in 2 and higher dimensions. However, Voronoi diagrams for circles in 2D and spheres in 3 or higher dimensions in the Euclidean distance measure, which we call Euclidean Voronoi diagrams, have not been explored as much as they deserve even though they have significant potential impacts on diverse applications in both science and engineering [1,11,22,28,35]. We want to note here that a Euclidean Voronoi diagram of circles or spheres is usually referred to as an additively weighted Voronoi diagram in the computational geometry community [24]. However, in this paper we will use the term a Euclidean Voronoi diagram of circles or spheres to emphasize its geometric interpretations and immediate practical applications in various scientific and engineering disciplines. We believe that a Euclidean Voronoi diagram for spheres in 3D can have more significant applications in various fields than its planar counterpart since most important geometric problems are 3D in nature. For example,
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
Fig. 1. Examples of the Euclidean Voronoi diagram of circles: (a) random circles, and (b) circles on a spiral.
the structural analysis of proteins or RNA’s requires an efficient computational tool to analyze the spatial properties among atoms [11,16,28]. In the design of new materials, a similar analysis is fundamental as well [15,23]. However, due to the lack of appropriate algorithms and stable running codes for Euclidean Voronoi diagrams of circles and spheres, most applications instead adapt an ordinary Voronoi diagram of points or a power diagram for the circles or spheres. For Euclidean Voronoi diagrams of circles in a plane, a simple yet fast and robust algorithm is available [13,14]. See Fig. 1 for the examples of computed output for circle sets. In Kim et al. [13,14], a Euclidean Voronoi diagram for circles is computed via an edge-flipping procedure, which takes O(n2) time in the worst-case where n is the number of circles, using the ordinary Voronoi diagram of center points as an initial starting diagram [13,14]. Note that this diagram has been used in various applications such as a cable packing problem [33]. Since the edge-flipping in 3D, however, lacks one-to-one correspondence between before and after flipping [24], a direct application of the edge-flipping concept may not be appropriate for the Euclidean Voronoi diagram for spheres in 3 or higher dimensions unless special cares are taken. Therefore, developing a good practical algorithm for a Euclidean Voronoi diagram of spheres in 3D is very important for practical reasons. Unlike other types of Voronoi diagrams or power diagrams, there are only very few previous works on the Euclidean Voronoi diagram of spheres. Aurenhammer discussed the possibility of transforming the computation of the Euclidean Voronoi diagram of spheres in d-dimension to one of (dC1)-dimensional power diagram construction by the intersections of (dC2)-dimensional half-spaces [2]. However, neither an algorithm nor an implementation is given for dR3. Gavrilova, in her Ph.D. Thesis, reported important properties of Euclidean Voronoi diagrams in arbitrary dimemsions, including shapes of the Voronoi regions, nearest neighbors and empty-sphere properties [8,9]. Will wrote a comprehensive Ph.D. thesis dedicated to the computation of the Voronoi regions for the Euclidean Voronoi diagram of spheres [36]. In his thesis, Will showed that a Voronoi region can be computed in Q(n2)
1413
Fig. 2. Examples of EVD(S): (a) an a-helix consisting of 67 atoms, and (b) PDB ID 1BH8 consisting of 1074 atoms (680 C’s, 181 N’s, 203 O’s, and 10 S’s).
time and actually implemented the algorithm so that the complete Voronoi diagram can be computed by constructing one Voronoi region after another. Luchnikov et al. proposed a practical idea of tracing edges which is simple yet powerful method for obtaining a Voronoi diagram of spheres in 3D [21]. However, a correct and efficient implementation of their algorithm is not so easy and the details of their algorithm are not available. Recently, we reported on the full implementation of an edge-tracing algorithm for constructing the Voronoi diagram with discussions of various applications including the analysis of protein structures [16,17]. In this paper, we will discuss several properties of the Euclidean Voronoi diagram for spheres in 3D and present the details of an algorithm based on the edge-tracing algorithm which computes the Voronoi diagram by tracing Voronoi edges. Given 3D balls, we describe how to trace the edges correctly and efficiently and discuss the representations of edges and faces. Having fully implemented the algorithm, we provide some experimental results as well. Shown in Fig. 2(a) is a Euclidean Voronoi diagram computed from a set of 67 atoms which describes an a-helix of a real protein. Fig. 2(b) illustrates a Voronoi diagram for a complete protein data with the PDB ID 1BH8 [40] which consists of 1074 atoms (680 C’s, 181 N’s, 203 O’s, and 10 S’s). In both Voronoi diagrams, faces are not intentionally rendered for the purpose of visualization. This paper is organized as follows. In Section 2, we define some terminology and notations used in this paper. Then, we discuss geometric primitives such as vertex positions, edge equations, and face equations of a Voronoi diagram for 3D balls in Section 3. The edge-tracing construction for the topology of the Voronoi diagram is described in Section 4. In Section 5, we conclude this paper with the discussions on the experimental performance of the proposed algorithm and some future works.
2. Terminology and definitions Let BZ{b1,b2,.,bn} be a set of generators for a Voronoi diagram where bi is a 3-dimensional spherical ball. Hence, biZ(ci,ri) where ciZ(xi,yi,zi) and ri denote the center and
1414
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
the radius of a ball bi, respectively. We assume that no ball is completely contained inside another ball even though intersections are allowed between balls. Associated with each ball bi, there is a corresponding Voronoi region VRi for bi, where VRi Z fpjdistðp; ci ÞK ri % distðp; cj ÞK rj ; i sjg. Then, EVD(B)Z{VR1, VR2,., VRn} is called a Euclidean Voronoi diagram for B. The topology of the Voronoi diagram in this paper is represented as an edge-graph GZ(V,E) where VZ{v1,v2,.} and EZ{e1,e2,.} are sets of Voronoi vertices and edges, respectively. Once G is available, the face set FZ{f1,f2,.} of Voronoi diagram can be easily constructed. In this paper, the ordinary L2-metric is used to define a Euclidean Voronoi diagram. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi In other words, 2 2 2 distðci ; cj ÞZ ðxi K xj Þ C ðyi K yj Þ C ðzi K zj Þ . As in ordinary point set Voronoi diagrams, Voronoi regions corresponding to balls on the boundary of the convex hull of B are unbounded. Other regions are bounded by a set of faces, called Voronoi faces, where a Voronoi face f is defined by two immediately neighboring balls. Note that a Voronoi face is always a connected-subset of hyperboloid of two sheets. A Voronoi face intersects another Voronoi face to form a Voronoi edge e which is a conic section. When Voronoi edges intersect, a Voronoi vertex v is defined. Note that there is always one-to-one correspondence between a Voronoi vertex v and a sphere S tangent to neighboring balls defining the vertex. Hence a tangent sphere is always centered at the corresponding vertex v. The tangent sphere S is called empty since it never intersects or contains any other ball except at the tangent points with its generating balls. In this paper, we assume that the degree of a Voronoi vertex is always four which means that no five balls are cotangent to an empty tangent sphere. We also assume that a Voronoi edge is always the common intersection among three Voronoi faces and therefore only three balls define an edge. Even with these assumptions, topological conflicts may occur due to numerical errors in the course of the Voronoi diagram construction [31,32]. It is known that the exact computation techniques can resolve such topological inconsistencies [31,32,37]. In this paper, however, we rather focus on describing the Voronoi diagram construction assuming that spheres are in general positions and leave the exact computation issue for the future research. Balls, in this paper, refer to spherical generators defining a Voronoi diagram while tangent spheres denote spheres corresponding to Voronoi vertices. For convenience of notation, we will ignore the term Voronoi and call them simply a vertex, edge or face unless further definition or explanation is necessary. In a similar fashion, the term Euclidean will also be ignored. Shown in Fig. 3 is a 2-dimensional analogy of a small Voronoi diagram of 3D balls. An edge e is oriented from a start vertex vs to an end vertex ve. The empty tangent spheres corresponding to vs and ve are called a start tangent sphere Ss and an end tangent sphere Se, respectively. An oriented edge
Fig. 3. Euclidean Voronoi diagram of spheres and related terminology.
e is always defined by three neighboring balls which surround e, and these balls are called gate balls bg1,bg2 and bg3 (In the figure illustrating its 2D analogy, only two gate balls are shown.) They are called gate balls since the triangle defined by the centers of the balls can be viewed as a gate which opens to the edge direction. The other two balls bs and be are called a start ball and an end ball since they uniquely correspond to the start and end vertices, respectively. To compute an ending vertex ve for an edge e, it is necessary to do computations with some balls in the given ball set B. Given three gate balls bg1, bg2 and bg3 and a start ball bs, it turns out that the computation of the tangent sphere Se, corresponding to an end vertex ve, requires to test if Se can be defined among four balls consisting of three gate balls bg1, bg2 and bg3 and each ball in KZB\{bg1,bg2,bg3}. We call K a candidate ball set for an edge e. In other words, in the example, KZ{bs,be,bi,bj}. Note that K contains the start ball bs as well (We will discuss this later). Associated with each edge, an angular distance q is defined as another useful distance measure. Given three gate balls, let bg1 be the smallest. Then, an angular distance q of a point p from a vertex vs is the angle :vscg1p, where cg1 is the center of bg1. Note that 0!q!2p and q is a directed distance where a positive direction is defined by the edge direction from the start vertex vs. When there are a number of angular distances defined by different balls, we call the one giving the minimum angle the closest to the start vertex. 3. Geometric primitives In the computation of a Voronoi diagram of circles in a plane, the computations of geometric primitives such as the coordinates of the vertices and the equations of the edges have been separated from the construction of the topology between vertices and edges [13,14]. In the current problem for spheres in 3D, we take a similar approach. In 3D, the geometry computation consists of vertices, edges, and faces while its 2D counterpart consists of vertices and edges. We want to mention here that the computation of vertices is inevitable in the construction of the topology. However, the computation of equations for edges and faces
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
can be avoided in the construction of a valid topology of a Voronoi diagram. Once a topology is constructed with the coordinates of vertices, they can be separately computed without much effort when it is necessary. 3.1. Voronoi vertices From the definition of a Voronoi diagram, a Voronoi vertex is the center of an empty sphere tangent to four nearby balls. Gavrilova and Rokne [9] recently reported an elegant algorithm to compute the sphere(s) tangent to the four balls in a general setting and we briefly summarize it here as follows. Suppose that b1,.,b4, denote four balls from which a Voronoi vertex, and therefore a corresponding tangent sphere, is to be computed. Note that the center of the tangent sphere is the equi-distant point from the four given balls. Without loss of generality, we can assume that all balls are reduced by the radius of the smallest ball b4. Then, b4 is reduced to a point. Transforming the four-ball system so that b4 coincides with the origin, the tangent sphere S with a center (x,y,z) and a radius r can be computed by solving the following system. Let a transformed shrunken ball bi have a center (xi,yi,zi) and a radius ri, for iZ1,.,4. ðx K x1 Þ2 C ðy K y1 Þ2 C ðz K z1 Þ2 Z ðr C r1 Þ2
(1)
ðx K x2 Þ2 C ðy K y2 Þ2 C ðz K z2 Þ2 Z ðr C r2 Þ2
(2)
ðx K x3 Þ2 C ðy K y3 Þ2 C ðz K z3 Þ2 Z ðr C r3 Þ2
(3)
x2 C y 2 C z2 Z r 2 :
(4)
By reducing the system to a simpler form, the first three equations are reduced to three linear equations in terms of r. In other words, each coordinate, x, y, or z, is expressed by a linear equation of an unknown variable r. Then, substituting these equations into the fourth equation results in a quadratic equation in terms of r. By looking at the rank of the linear system consisting of three linear equations, they showed that the four ball configuration can be categorized into one of four different cases with none, one, two or infinite tangent spheres. For a detailed explanation, the readers is referred to Gavrilova and Rokne [9]. 3.2. Voronoi edges Once the edge-graph GZ(V,E) of the topology of a Voronoi diagram is computed, it is often necessary to compute the actual equations of Voronoi edges. Unless it is a degenerate case, an edge is always defined as a locus of points equi-distant from the surfaces of three balls b1,b2, and b3. Similar to a Voronoi vertex, we shrink the balls by the radius of the smallest ball and transform the ball-system so that the smallest ball shrunken to a point coincides with the origin. Let the transformed shrunken balls have centers
1415
(xi,yi,zi) and radii ri, iZ1,2 and 3. Note that the smallest ball b3 now has a radius of zero and coincides with the origin. Then, the equation of an edge is the solution of the following three equations. ðx K x1 Þ2 C ðy K y1 Þ2 C ðz K z1 Þ2 Z ðr C r1 Þ2
(5)
ðx K x2 Þ2 C ðy K y2 Þ2 C ðz K z2 Þ2 Z ðr C r2 Þ2
(6)
x2 C y2 C z2 Z r 2 ;
(7)
where (x,y,z) is a point on the edge. Rearranging the above system yields the following. x1 x C y1 y C z1 z K ðx21 C y21 C z21 K r12 Þ=2 C r1 r Z 0
(8)
x2 x C y2 y C z2 z K ðx22 C y22 C z22 K r22 Þ=2 C r2 r Z 0
(9)
x2 C y2 C z2 Z r 2
(10)
Two linear equations in Eqs (8) and (9) show that they are planes and the solution of both equations is a line in general. Since the third equation is a sphere, the solution of the whole system is two points on the sphere in general. Hence, the edge equation can be obtained by a variable r. Looking back at the linear equations once more, bothffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 C y2 C z2 planes move at constant velocities of Kr = x 1 1 1 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 and Kr2 = x2 C y2 C z2 in the normal direction of both planes, respectively. This means that the trajectory of the moving intersection line defines another plane. Hence, its intersection with the third equation with a variable r, which is the resulting Voronoi edge, has to be a planar curve. Similarly, two latter equations together define a hyperboloid of two sheets. Therefore, a Voronoi edge is an intersection between a hyperboloid and a plane, and indeed is a conic curve. It can be also differently shown by the fact that a Voronoi edge satisfies the definition of a spine curve of Dupin cyclide [25,26]. Will has also shown that a Voronoi edge is a conic curve and further classified it by looking at the relative configuration among these balls [36]. Even though the type of the curve is known, its mathematical representation is another issue since applications usually require processing of these curves. Hence, a convenient curve representation is necessary for the processing. Since Voronoi edges are conics, they can be exactly represented in a form of a rational quadratic Be´zier curve defined as bðtÞ Z
w0 ð1 K tÞ2 p0 C 2w1 ð1 K tÞtp1 C w2 t2 p2 ; t 2½0; 1; w0 ð1 K tÞ2 C 2w1 ð1 K tÞt C w2 t2 (11)
where p0, p1, and p2 are the control points, and w0, w1, and w2 are the corresponding weights, respectively [7]. It is known that the above equation can be completed once all five of the following parameters are known: two end points, tangent vectors at both end points, and a point through
1416
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
becomes a tangent plane of B12 at v. As a result, a vector starting from v and ending at p has the same angle with vectors starting from v and ending at c1 and c2. In a similar manner, the angle bisecting plane between b2 and b3 becomes a tangent plane at v, and therefore every vector starting from v on the tangent plane has the same angle with vectors starting from v and ending at c2 and c3. Therefore, the intersection between the two tangent planes becomes the tangent line of a Voronoi edge e at v and it is equi-angular with the three vectors starting from v and ending at c1, c2, and c3. , Fig. 4. Proof of Lemma 1.
which the curve passes [7]. In our problem, two end points are already given as Voronoi vertices of the edge. Let b1, b2, and b3 be three gate balls for an edge e where r1%r2%r3, and v be a vertex from which we want to define a tangent vector of the edge. Suppose that ci is the center of a ball bi and ti Z ðci K vÞ=jjci K vjj. Then, the following lemma holds. Theorem 1. (Angle Trisectioning) The tangent vector t at a point v on a Voronoi edge e is equi-angular with three vectors t1,t2 and t3. Proof. An edge e is a common intersection among three bisecting surfaces, i.e. between b1 and b2, b2 and b3, and b3 and b1. Let us consider a bisector B12 between b1 and b2. Let v be an arbitrary point on e, and let a plane P bisect the angle :c1vc2 and p be a point on P other than v as shown in Fig. 4. In addition, we locate a point a on a line segment vc2 such that jjvKajjZjjvKc1jj. Then, jjpKc1jjZjjpKajj holds since Dvpa and Dvpc1 are congruent. Since Dpc2a forms a triangle, the triangular inequality holds, i.e. jjpK c2 jj% jjpK ajjC jjaK c2 jj. Since jjpKajjZjjpKc1jj and jjaKc2jjZr2Kr1, the above inequality can be rearranged as jjpK c2 jjK r2 % jjpK c1 jjK r1 . Therefore, P cannot intersect B12 at any point except at v unless r1Zr2, and the angle bisector P
Due to this lemma, the following corollary immediately holds and is presented without a proof. Corollary 2. The tangent vector t at v satisfies the following system of equations. ti $t Z tj $t Z tk $t;
(12)
where ts0. The last parameter, a passing point p, can be found as follows. Suppose that we define a plane P passing through the centers of three gate balls bg1, bg2 and bg3. Then, the intersection of P with the three balls results in three circles, Cg1, Cg2, and Cg3, on P as shown in Fig. 5(a). Therefore, the passing point p, which lies on this plane, has to be the center of a circle C tangent to these three circles. Hence, computing C is equivalent to solving The Apollonius 10th Problem in a plane P and its solution process is well-described in Kim et al. [14] and Rokne [29]. Once these five parameters are known, the complete rational quadratic Be´zier representation can be found by solving the following equation. Suppose that a curve is represented by the standard form of a rational quadratic Be´zier, i.e. w0Zw2Z1. Then, w1 can be obtained by t1 ; w1 Z pffiffiffiffiffiffiffiffiffi 2 t0 t2
Fig. 5. Passing point on a Voronoi edge solved by The Apollonius Tenth Problem: (a) one tangent circle, and (b) two tangent circles.
(13)
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
where t0, t1, and t2 are the barycentric coordinates of a passing point with respect to the triangle vertices b0, b1, and b2 [7]. Note that the barycentric coordinates are positive if a passing point is located inside a triangle Db0b1b2. In addition, w1 can be negative if the passing point is not inside the triangle. If p lies on a valid edge segment, w1 computed by the above equation can be directly applied. If p lies outside a valid edge segment, however, the correct rational quadratic Be´zier representation of the valid edge segment can be obtained by reversing the sign of w1 calculated by the above equation using p. There could be, however, even worse cases. When an edge is elliptic, there could be two such intersection points, p1 and p2, between P and the edge e as shown in Fig. 5(b). In this case, the situation can be either (i) both are on a valid edge segment, (ii) both are not on a valid edge segment, or (iii) only one of them is on a valid edge segment. The key part of this problem is how to tell which points, p1, p2, or both, are on a valid segment. Suppose we pick an arbitrary point between two points, for example p1. Then, checking if p1 lies on a valid edge segment can be easily done by testing if the angular distance of p1 from the start vertex vs is less than the angular distance of the end vertex ve from vs. Note that ve is already computed. If p1 is on a valid edge segment, the weight w1 computed by Eq. (13) can be directly used. If p1 is not on the valid segment of the edge, the sign of computed weight w1 should be reversed regardless of the computed sign of w1. From the above discussions for both cases, the following lemma can be concluded. Lemma 3. Suppose that the weight w1 is computed by Eq. (13). If a passing point p is on the valid edge segment, w1 can be directly used as the weight. Otherwise, its sign should be reversed.
1417
b Z K4ri2 C 4y2j K 4rj2 K 8yi yj C 8ri rj C 4y2i ; c Z 4z2i K 4rj2 K 4ri2 C 8ri rj C 4z2j K 8zi zj ; d Z K8xi yj K 8yi xj C 8xj yj C 8xi yi ; e Z 8yj zj C 8yi zi K 8yi zj K 8zi yj ; f Z 8xj zj K 8xi zj K 8zi xj C 8xi zi ; g Z K8xi ri rj K 8xj ri rj C 4xi z2j K 4x3j K 4x3i C 4xi rj2 C 4xi ri2 K 4xj z2j K 4xi y2i K 4xi z2i C 4xj ri2 C 4xi x2j C 4xi y2j K 4xj y2j C 4y2i xj C 4x2i xj C 4z2i xj C 4xj rj2 ; h Z K8yi ri rj K 8yj ri rj K 4y3i K 4y3j C 4y2i yj C 4yj ri2 C 4yj rj2 C 4yi z2j C 4yi y2j K 4x2j yj K 4x2i yi C 4yi rj2 K 4yj z2j C 4z2i yj C 4yi ri2 K 4yi z2i C 4x2i yj C 4yi x2j ; k Z 4zj ri2 C 4y2i zj C 4zi rj2 K 4y2i zi K 4z3i C 4zj rj2 K 4z3j K 4y2j zj K 8zi ri rj K 8zj ri rj C 4z2i zj C 4zi ri2 K 4x2i zi C 4x2i zj K 4x2j zj C 4zi x2j C 4zi z2j C 4zi y2j ; and l Z 2x2i y2i C 2x2i z2i K 2x2i x2j K 2x2i y2j K 2x2i z2j K 2x2i ri2 K 2x2i rj2 C 2y2i z2i K 2y2i x2j K 2y2i y2j K 2y2i z2j K 2y2i ri2 K 2y2i rj2 K 2z2i x2j K 2z2i y2j K 2z2i z2j K 2z2i ri2 K 2z2i rj2 C 2x2j y2j C 2x2j z2j C 4x2i ri rj C x4i C y4i C z4i C x4j C y4j C z4j C ri4 C rj4 C 4y2i ri rj C 4z2i ri rj C 4x2j ri rj C 4y2j ri rj
3.3. Voronoi faces
C 4z2j ri rj K 2x2j ri2 K 2x2j rj2 C 2y2j z2j K 2y2j ri2 K 2y2j rj2 Once the edge-graph GZ(V,E) is known, it is easy to tell which pair of balls define a Voronoi face. Two topologically neighboring balls bi and bj define a Voronoi face as jjp K ci jj K ri Z jjp K cj jj K rj
(14)
where ci and ri are the center and radius of ball bi. Then, the above equation defines a hyperboloid of two sheets. Rearranging the above equation results in the following implicit equation. ax2 Cby2 Ccz2 Cdxy Ceyz Cfzx Cgx Chy Ckz Cl Z 0; (15) where a Z 4x2j K 4rj2 K 8xi xj C 8ri rj C 4x2i K 4ri2 ;
K 2z2j ri2 K 2z2j rj2 K 4ri3 rj C 6ri2 rj2 K 4ri rj3 :
ð16Þ
Given a face equation, it is usually better to represent the faces in parametric form for reasons that are well-known in the geometric modeling community. Since the Voronoi face in our problem is always a hyperboloid, the natural choice for the parametric representation is a rational Be´zier surface with an appropriate degree and topology. It turns out that the topological shape of a Voronoi face in our problem form a polygon such as a triangle, rectangle, pentagon, etc. In addition, it can even have holes inside. Therefore, the Be´zier form cannot be used directly but needs approprite modifications. Hence, there are three alternative ways to represent a Voronoi face as a rational Be´zier form. If a rectangular topology is used for rational Be´zier surface representation, it is in general necessary to trim
1418
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
the surface for a non-rectangular shaped Voronoi face. It is known that trimming of a parametric surface takes a significant amount of computing time and memory. Most seriously, there is usually no way to avoid the gaps between two neighboring surfaces if trimming has to be employed. If a Voronoi face is subdivided into a number of triangular segments on a hyperboloid, each triangle can be represented as a rational Be´zier surface having triangular topology. If this approach is used in our problem, trimming is not necessary. However, it may be necessary to implement some geometric operations redundantly. For example, consider the rolling-ball blending surface between balls. It is known that this surface can be conveniently represented using a rational surface having rectangular topology. Hence, having both triangular and rectangular topologies for rational surfaces simultaneously does not seem to be a good way to represent Voronoi faces. If each triangle in the above is again subdivided into three rectangles, all of them can be represented as rational Be´zier surfaces having rectangular topology. While avoiding the problem discussed above, this approach requires a significantly larger amount of memory. While the choice of an appropriate surface representation remains to be made, principal uses of Voronoi diagrams in most applications suggest a different alternative. The main uses of Voronoi faces are the visualization of the faces themselves, the computation of volumes and boundary areas of Voronoi regions, etc. Suppose that two balls b1 and b2 define a Voronoi face where r1%r2. In addition, let pm be the mid point between the centers of b1 and b2. Then, we transform both balls so that pm is located at the origin and the centers of b1 and b2 are placed on the positive and negative Z-axis, respectively, as shown in Fig. 6. We call this position of two balls a canonical position. Let b1 and b2 be centered at (0,0,d) and (0,0,Kd), respectively. Then, the equation for the bisector surface between the balls in the canonical position is
Fig. 6. Voronoi face in the canonical position.
represented as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 C y2 C ðz K dÞ2 K r1 Z x2 C y2 C ðz C dÞ2 K r2 : (17) While the above equation defines a hyperboloid of two sheets, the valid Voronoi face between two balls is always one of either sheet, not both. Rearranging the above equation in an explicit form gives the following equation. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2r ðKr 2 C 4d2 Þð4x2 C 4y2 C 4d 2 K r 2 Þ z ZG ; (18) K4r 2 C 16d 2 where rZr2Kr1. As shown in the above equation, the equation is divided into the positive and negative parts, and the desired sheet of a hyperboloid is the positive part. In Fig. 6, for example, fV is the valid sheet and fI is the invalid one. When a Voronoi face is represented and manipulated in applications, it is always necessary to test the validity of the face. For example, when we have to triangulate a Voronoi face for a visualization purpose, we have to choose one of either points on the hyperboloid computed from the above equation when X and Y coordinates are given. Since we have transformed the Voronoi face in the canonical position, we can easily locate the correct surface by simply looking at the signs of points on the hyperboloid. If Voronoi faces are handled in this way, both evaluating a point on the face and verifying if a point is on the face become much simpler. 3.4. Visualization of Voronoi faces To illustrate how to process Voronoi faces in applications, we show how to render those faces as an example. Shown in Fig. 7 are examples similar to those presented in Luchnikov et al. [21]. Fig. 7(a) shows an example of an atom complex which consists of 15 atoms with identical radii, while Fig. 7(d) shows a similar example with three different radii. Their Voronoi diagrams are computed in Fig. 7(b) and (e), respectively. For convenience of visualization, we have shown only the Voronoi region for the yellow atom in the center while ignoring edges and faces emanating toward infinity. Fig. 7(c) and (f) show the Voronoi regions only. Note that the faces in Fig. 7(c) are planar and the associated edges are line segments like an ordinary point set Voronoi diagram since the sizes of atoms are identical. On the other hand, the faces in Fig. 7(f) are parts of hyperboloids and the edges are conics represented as rational quadratic Be´zier curves. To render a Voronoi region, it is necessary to triangulate the faces bounding the region so that a graphics library such as OpenGL can be used. To triangulate a Voronoi face, we need to sample some points on the valid part of the hyperboloidal surface associated with the Voronoi face. Since a Voronoi face is a subset of a hyperboloid of two
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
Fig. 7. Fifteen balls with three different radii and the corresponding Euclidean Voronoi diagram. Voronoi edges are represented as rational quadratic Be´zier curves while Voronoi faces are represented as implicit surfaces: (a) 15 balls with identical radii, (b) the corresponding Voronoi diagram where faces and edges going to infinity are ignored, (c) the Voronoi region for the yellow ball, (d) 15 balls with three different radii, (e) the corresponding Voronoi diagram, and (f) the Voronoi region for the yellow ball. Note that the polyhedron is the Voronoi region for the yellow ball which has 14 faces with the other balls.
sheets, such a sampling operation is not an easy task. After correctly sampling such an appropriate point set, it is also necessary to do a triangulation and then eliminate extraneous triangles. It turns out that performing these tasks directly in 3D is relatively difficult as well. Instead, we propose a different approach to take advantage of the surface representation presented in the previous subsection. Given a Voronoi face in the world coordinate system, there are two related balls. Suppose that we transform both balls in the canonical position. Then, the Voronoi face between the transformed balls is always a single-valued function w.r.t. X and Y. Note that a Voronoi face is bounded by a number of Voronoi edges represented as rational quadratic Be´zier curves. Therefore, the edges can be also easily transformed to the same canonical position. Let f be a Voronoi face defined by two balls bi and bj in the world coordinate system. Then, transforming bi and bj to the canonical position produces the corresponding Voronoi face fV as well in the canonical position. In the canonical position, the valid points on the correct sheet of hyperboloid for the face can be easily found. Let EZ{ekjkZ1,2,.} be an edge set bounding f in the world coordinate system where ek is represented as a rational quadratic Be´zier curve. Suppose now that Ec be the edge set of fV in the canonical position. Then, the projection of Ec as well as fV onto the XY-plane produces a polygon Poly bounded by a number of conics in a plane. Note that the transformation to a canonical position and projection of boundary curves can be done easily by simply processing the control points of edges since Be´zier curves and surfaces are affine invariant.
1419
Suppose that we sample a number of points in Poly bounded by Ec and a number of points from Ec itself. Then, we apply a constrained Delaunay triangulation using the polylines on the boundary as constraints and remove extraneous triangles exterior to the boundary Ec of the polygon. If we lift-up the sample points and the topology among the sample points onto fV, we can get the correct triangulation of Voronoi face fV in the canonical position. Hence, the back-transformation of fV to f produces the desired triangular mesh in 3D so that a graphics library such as OpenGL can be used. For the constrained Delaunay triangulation, we have used CGAL library [38]. Note that the polygon Poly is not necessarily convex. It may even contain some holes inside. Hence, we need to be a little bit careful to sample appropriate points inside Ec. The optimal distribution of sample points itself from both inside and boundary curves can be another issue to be studied further. From our experience, sampling points from uniform grids is quite satisfactory. While the constrained Delaunay triangulation in CGAL is too general for the given problem and therefore takes a significant amount of computation, our triangulation problem is quite a special case. Therefore, devising a specialized algorithm would be beneficial.
4. Topology construction by tracing edges Given the building blocks to compute the geometric part of a Voronoi diagram, it is necessary to devise an algorithm to construct its topology counterpart. The basic idea of the edge-tracing algorithm is quite simple. The algorithm first locates a true Voronoi vertex v0 by computing an empty tangent sphere defined by four appropriate nearby balls. Provided that v0 has been found, four edges e0, e1, e2, and e3 emanating from v0 can be easily identified and pushed into a stack called an Edge-stack. Hence, these edges have v0 as their starting vertices. After popping an edge from the stack, the algorithm computes the end vertex of the edge. Note that a vertex can be found by computing an appropriate empty tangent sphere from each of candidate balls plus three gate balls defining the popped edge. If an appropriate empty tangent sphere is found, the center of the sphere becomes the end vertex of the popped edge. Once the end vertex of currently popped edge is found, it is also possible to compute three more edges emanating from this new vertex. Hence, when an edge is created, its end vertex is used as the start vertex of another three new-born edges. (Some of the new-born edges may turn out to be invalid later on.) Note that these edges are also pushed into the Edge-stack. By repeating this process until the Edge-stack is empty, the computation of a Voronoi diagram of a connected graph is completed. Hence, this process is iterated as many times as the number of edges, which is O(n2) in the worst-case and O(n)
1420
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
on the average, where n is the number of balls. Even though the idea is simple yet powerful, designing a correct and efficient algorithm is not so easy. Hence, we elaborate the details of the algorithm step by step as follows. 4.1. Initialization The initialization process computes a tangent sphere from each combination of four balls and tests if the tangent sphere is empty or not. Since there are O(n4) number of fourball combinations and a tangent sphere from each combination should be tested for emptiness against nK4 other balls, a brute-force approach to find an initial vertex v0 can take as long as O(n5) time in the worst-case. A simple trick can help to avoid such a high worst-case time complexity. Suppose that we add four more fictitious balls to the initial ball set. These fictitious balls are intentionally configured so that there is only one empty tangent sphere among them. Let the Voronoi vertex corresponding to this tangent sphere be a fictitious vertex. Then, these fictitious balls are located in the neighborhood of the input ball set so that the Voronoi vertex among the fictitious balls can still remain as a Voronoi vertex for all balls including the fictitious ones. Then, we can safely start the process from this fictitious initial vertex, which takes O(1) time for the initialization. Tracing the edges from the fictitious initial vertex will eventually lead to a Voronoi vertex which is defined only among the input balls. Then, we choose this vertex as the true initial vertex to start the process to construct the topology of the Voronoi diagram by ignoring the four fictitious balls from the ball set and the computations done until now. We want to mention that the proposed approach of computing an initial vertex by adding four fictitious balls not only reduces the worst-case time complexity but also avoids a potentially difficult problem of dealing with the cases where there is none, one, or two possible tangent spheres corresponding to end vertex(es) for a given ball configuration. 4.2. Tracing edges via the edge-stack When an initial vertex v0 is located, we define four edges e0, e1, e2, and e3 of the Voronoi diagram emanating from v0 and push them into the Edge-stack. Hence, these edges have v0 as their starting vertices. Once the edges are stored in the Edge-stack, we pop one edge after another to trace. This tracing is done by completing the definition of a popped edge by computing its end vertex. Note that the end vertex can be found by computing an appropriate empty tangent sphere from 3 gate balls surrounding the popped edge plus one appropriate ball of the nK3 candidate balls. Suppose that a valid end vertex vi is computed for a popped edge. Then, we again define three more edges ei1,
ei2, and ei3 emanating from vi, and again push them into the Edge-stack. Hence, these new-born edges have vi as their starting vertices as well. However, some of these new-born edges may be determined in a later step to be redundant since their end vertices are already created. This will be explained in Section 4.4. By repeating this process until the Edge-stack is empty, the coordinates of vertices and the topology among vertices and edges can be completely identified. Hence, pushing and popping operations on edges with the Edge-stack lead to the complete edge-graph of a Voronoi diagram for balls. While the computation of vertices, each of which corresponds to a tangent sphere among the four balls, is inevitable in the process of constructing the edge-graph, the computation of edge and face geometries is not necessary in the process of edge-graph computation. 4.3. End vertex computation Computing an end vertex of a popped edge is equivalent to finding an end ball associated with the end vertex. Hence, we compute a tangent sphere from four balls, three defining balls of the popped edge and one candidate ball in the ball set, to find the end vertex. After a tangent sphere is computed, its emptiness is tested against all the other balls except the four balls defining the tangent sphere. Here, the cardinality of the candidate ball set is nK3, not nK4, consisting of all balls except only the three surrounding balls for the edge being traced. In other words, the start ball should be also considered as a candidate ball for the computation of the end vertex. This is because a certain configuration of four balls can define two tangent spheres, not only one, as described earlier. Hence, the next lemma follows without a proof. Lemma 4. For an edge e, and therefore three corresponding gate balls bg1, bg2 and bg3, and an input ball set B, the candidate ball set K is defined as KZB\{bg1,bg2,bg3}. Due to the lemma, we can design an algorithm to compute an end vertex given an edge e with its gate balls and the start ball. The following is a simple-minded first version of such an algorithm. Algorithm END-VERTEX-Q Step 1. For bk2K, where K is the candidate ball set. Step 1.1. Compute a tangent sphere Sk from bk and three gate balls Step 1.2. For bl2K\{bk}, Step 1.2.1. If Sk intersects bl then GOTO Step 1. End-for Step 1.3. Enqueue Sk into Q End-for Step 2. If Q is empty, the end vertex is going infinity. Step 3. Else, find S in Q closest to the start vertex, and use the center of S as the end vertex of the current edge.
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
In the brute-force algorithm presented above, the first For-loop in Step 1 iterates nK3 times to compute tangent spheres, and the nested For-loop in Step 1.2 iterates nK4 times for the emptiness test of each tangent sphere against other balls. Since the output from Step 1 is O(n) of empty spheres, Step 3 takes O(n) time in the worst-case. Therefore, algorithm END-VERTEX-Q takes O(n2) time in the worstcase and therefore the next lemma follows. Lemma 5. Using Algorithm END-VERTEX-Q, the edge tracing algorithm constructs a complete edge-graph of a Voronoi diagram of spheres in O(mn2) time, where m is the number of Voronoi edges and n is the number of balls. The above time complexity, however, can be improved by the following lemma. Lemma 6. Given a start vertex and gate balls for an edge e, the end vertex of e can be located in O(n) time in the worstcase. Proof. (See Fig. 8) Suppose that an edge e is being traced starting from a vertex vs, and a ball bi, together with three
1421
gate balls bg1, bg2 and bg3, successfully defines a tangent sphere Si. However, Si may or may not define a correct end vertex for e depending on the other balls. As shown in Fig. 8(a), Si may not be empty but intersect another ball bj, i!j. In this case, the center vei of Si cannot be the end vertex of the edge at all since the angular distance qi of vei from the start vertex vs is greater than qj associated with Sj (See Fig. 1(a)). Hence, vej, the center of tangent sphere Sj defined by the ball bj can be now considered as a candidate for the end vertex of the edge. In Fig. 8(b), bj does not intersect Si and defines another tangent sphere Sj. In this case, which is similar to the above case, it can be shown that qi!qj. Hence, we can simply ignore bj from further consideration. In Fig. 8(c), bj does not intersect Si and also defines another tangent sphere Sj. In this case, it can be shown that qi>qj. Hence, we now consider vej for candidate end vertex instead of vei. Therefore, a scanning of nK3 candidate balls once will eventually locate the valid end vertex, and this process takes O(n) time in the worst-case as well as average-case. , Based on the above lemma, we can devise a more efficient algorithm as follows. Algorithm END-VERTEX-L Step 1. Compute a tangent sphere S1 from three gate balls and b1 in the candidate ball set K. Step 2. S)S1. Step 3. For bk2K\{b1}, Step 3.1. Compute a tangent sphere Sk from three gate balls and bk Step 3.2. If Sk is closer to the start vertex than S, then S)Sk End-for Step 4. Use the center of S as the end vertex of current edge. Given the three gate balls for a popped edge, Algorithm END-VERTEX-L first computes a tangent sphere Si with an arbitrary candidate ball bi. Then, the algorithm selects another candidate ball bj from the candidate ball set and computes a tangent sphere Sj with bj and three surrounding balls. If bj intersects Si, the current Si is replaced by Sj. If not, we choose the tangent sphere between Si and Sj which is closer to the start vertex of the popped edge in their angular distances. Since all balls in the candidate set is scanned only once, this process runs in O(n) time in the both the worst and average cases and finds the correct end vertex. Hence, based on Lemma 6, the following theorem holds.
Fig. 8. End vertex of an edge.
Theorem 7. Using Algorithm END-VERTEX-L, the edge tracing algorithm constructs a complete edge-graph of a Euclidean Voronoi diagram of spheres in O(mn) time, where m is the number of edges and n is the number of vertices.
1422
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
4.4. Checking if a vertex is already computed
5. Discussions and conclusions Voronoi diagrams for spheres in 3-dimensional space in Euclidean distance have many applications in various disciplines in science and engineering. Even though the Euclidean Voronoi diagram of spheres can be so important, its algorithms have not been studied as much as it deserves. In this paper, we proposed an edge-tracing algorithm to compute the 3-dimensional Eulcidean Voronoi diagram for spheres which runs in O(mn) time in the worst-case where m is the number of Voronoi edges and n is the number of given spheres. As building blocks, we have shown that Voronoi edges are conics and can be represented as rational quadratic Be´zier curves. We also discussed how to conveniently represent and process Voronoi faces which are parts of hyperboloids of two sheets. Fig. 9 illustrates the run-time behavior of the current implementation. The X- and Y-axis are the number of atoms in the example file and the computation time in seconds, respectively. The experiments were performed on a Pentium IV with the clock speed of 3.0 GHz and 2 GB memory. In the figure, the lower curves denote accelerations using
computation time (sec.)
Suppose that a new vertex for an edge is obtained by tracing an edge. If the vertex has not previously been computed, then we can safely use the new vertex to complete the edge as an end vertex. However, if it already exists, we have to be careful with the new vertex since it is redundant. Instead of creating the new vertex for the end vertex of the edge being traced, we rather search the existing one to complete the edge definition since the existing one already has associated topology information partially determined. Whenever a new vertex is computed, therefore, it is necessary to check if this vertex is already created or not. For the search of a vertex, we have devised a table of computed vertices and named it the Vertex Index Dictionary (VIDIC). Each entry, which corresponds to a vertex, in the dictionary consists of indices to four balls defining the vertex and a pointer to the definition of the vertex. Note that the number of entries in the dictionary is bounded by the number of vertices of a Voronoi diagram and therefore is O(n2) in the worst-case. However, we believe that a tighter worst-case bound for the size of VIDIC much lower than O(n2) exists in the worst-case as well as average-case since a vertex can be permanently removed from VIDIC once it is connected with four edges. While a simple linear search can be used on VIDIC, we prefer a binary search which takes O(log n) time in the worst-case, which is possible by maintaining the sorted list of the entries. In addition, we want to mention here that hashing, which takes O(1) time on the average, is also possible on VIDIC.
35 30 25 20 15 10
5 0 0
200
400
600
800
1000
1200
the number of spheres original
plane filter
bucket
bucket + plane filter
Fig. 9. Computation time statistics for various models with different sizes.
geometric filters and show significant performance enhancements. We expect that further increase in speed can be achieved by optimizing the codes. While the numbers of vertices, edges, and faces are O(n2) in the worst-case, their average numbers are known as linear to the number of balls. In addition, if we devise an appropriate bucket, we expect the quadratic growth of runtime behavior can be reduced to linear since computing an end vertex can take only O(1) time on the average if the balls are distributed accordingly. For this claim to be mathematically justified, however, the distribution of balls should be carefully investigated and it leaves a challenge. While an ordinary atomic structure yields a Voronoi diagram which is usually a singly connected edge-graph, it is possible that the valid edge-graph of a Voronoi diagram be disconnected for the given input balls. Let two balls bi and bj be non-intersecting and closely located to define a Voronoi face fij as shown in Fig. 10. Suppose that a third ball bk, much smaller than the other two, be located in-between bi and bj so that all three balls altogether define a single elliptic edge e, and the other balls are located far away from the edge e so that no vertex can be defined on e. Hence the edge e has no explicit vertex on it. Then, bi and bk define a face fik, and bj and bk define another face fjk. In this case, fik and fjk intersect at the elliptic edge e. Note that the edge e now defines an island in the face fij and is disconnected from the other part of edge-graph. Hence, in this case, there
Fig. 10. Big brothers problem for a disconnected case.
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
is no edge connecting the small ball to the other balls outside. If this case happens in a part of larger Voronoi diagram, the edge-graph of the whole Voronoi diagram is not a single edge-graph but disconnected. In this case, the isolated edge is either circular or elliptic. If the center point of the small ball is located on the line passing through the centers of two big-brothers, then the edge is circular. If the center of small ball is a little bit off the line, then the edge becomes now elliptic. In this case, special care should be taken to handle the situation to come up with a complete Voronoi diagram. Current work leaves many points of future research issues regarding on the Voronoi diagram construction. The acceleration of running performance of the proposed algorithm can be achieved by employing some geometric filters so that a search space can be significantly reduced. Another improvement would be to devise an algorithm with a better worst-case time complexity. Besides, developing an easy-to-implement algorithm for the Voronoi diagram would be plausible. The most important issue for the construction of a Voronoi diagram for a large model is to guarantee the robustness of the algorithm. In the current research, however, we have mainly focused on the computation of Voronoi diagram assuming that the generators are in general position and we leave the robustness issue for the immediate future research topic. To construct the Voronoi diagram robustly, the exact computation technique should be incorporated into the algorithm. This can be achieved either by using an exact computation package such as CORE library [39] or devising an exact computation scheme specific to our problem. The former approach is easy to come up with a robust code while the latter approach is more efficient in the run-time. Based on this Voronoi diagram, we expect that several important problems in biological systems such as molecular interaction can be solved rather easily and efficiently. In particular, the computation of interface between proteins or protein vs. small molecules can be available with Voronoi diagrams. Moreover, It is noteworthy that all these problems can be solved in a single framework, the Euclidean Voronoi diagram for spheres.
Acknowledgements This research was supported by Creative Research Initiative grant from the Ministry of Science and Technology, Korea. Authors thank Prof. Chee Yap and Prof. Kokichi Sugihara for suggesting the possibility of linear time reduction for the end ball search. Authors also would like to thank Dr Sang-Min Park for her valuable comments and careful review.
1423
References [1] Angelov B, Sadoc J-F, Jullien R, Soyer A, Mornon J-P, Chomilier J. Nonatomic solvent-driven Voronoi tessellation of proteins: an open tool to analyze protein folds. Proteins: Struct Funct Genet 2002;49(4): 446–56. [2] Aurenhammer F. Power diagrams: properties, algorithms and applications. SIAM J Comput 1987;16:78–96. [3] Aurenhammer F. Voronoi diagrams—a survey of a fundamental geometric data structure. ACM Comput Surv 1991;23(3):345–405. [4] Blum H. A transformation for extracting new descriptors of shape. Proceedings of symposium models for perception of speech & visual form 1967 p. 362–80. [5] Blum H, Nagel RN. Shape description using weighted symmetric axis features. Pattern Recogn 1978;10:167–80. [7] Farin G. Curves and surfaces for computer-aided geometric design: a practical guide. 4th ed. San Diego: Academic Press; 1996. [8] Gavrilova M. Proximity and Applications in General Metrics. Ph.D. Thesis: The University of Calgary, Dept. of Computer Science, Calgary, AB, Canada; 1998. [9] Gavrilova M, Rokne J. Updating the topology of the dynamic Voronoi diagram for spheres in Euclidean d-dimensional space. Comput Aided Geometric Des 2003;20(4):231–42. [11] Goede A, Preissner R, Fro¨mmel C. Voronoi cell: new method for allocation of space among atoms: elimination of avoidable errors in calculation of atomic volume and density. J Comput Chem 1997; 18(9):1113–23. [12] Held M. On the computational geometry of pocket machining. Berlin: Springer; 1991. [13] Kim D-S, Kim D, Sugihara K. Voronoi diagram of a circle set from Voronoi diagram of a point set: I. Topology. Comput Aided Geometric Des 2001;18(6):541–62. [14] Kim D-S, Kim D, Sugihara K. Voronoi diagram of a circle set from Voronoi diagram of a point set: II. Geometry. Comput Aided Geometric Des 2001;18(6):563–85. [15] Kim D-S, Chung Y-C, Kim JJ, Kim D, Yu K. Voronoi diagram as an analysis tool for spatial properties for ceramics. J Ceram Process Res 2002;3(3):150–2. [16] Kim D-S, Cho Y, Kim D, Cho C-H. Protein structure analysis using Euclidean Voronoi diagram of atoms. Proceedings of international workshop on biometric technologies (BT 2004). Special forum on modeling and simulation in biometric technology 2004 p. 125–9. [17] Kim D-S, Cho Y, Kim D. Edge-tracing algorithm for Euclidean Voronoi diagram of 3D spheres. Proceedings of 16th canadian conference on computational geometry 2004 p. 176–9. [18] Kirkpatrick DG. Efficient computation of continuous skeletons. Proceedings of 14th IEEE symposium foundations of computer science 1979 p. 18–27. [19] Lee DT. Medial axis transformation of a planar shape. IEEE Trans Pattern Anal Machine Intell 1982;4:363–9. [21] Luchnikov VA, Medvedev NN, Oger L, Troadec J-P. VoronoiDelaunay analyzis of voids in systems of nonspherical particles. Phys Rev E 1999;59(6):7205–12. [22] Montoro JCG, Abascal JLF. The Voronoi polyhedra as tools for structure determination in simple disordered systems. J Phys Chem 1993;97(16):4211–5. [23] Ohno K, Esfarjani K, Kawazoe Y. Computational materials science from ab initio to Monte Carlo methods. Berlin: Springer; 1999. [24] Okabe A, Boots B, Sugihara K, Chiu SN. Spatial tessellations: concepts and applications of Voronoi diagrams. 2nd ed. Chichester: Wiley; 1999. [25] Paluszny M, Boehm W. General cyclides. Comput Aided Geometric Des 1998;15(7):699–710. [26] Pottman H. personal communication. 2004. [28] Richards FM. The interpretation of protein structures: total volume, group volume distributions and packing density. J Mol Biol 1974;82: 1–14.
1424
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
[29] Rokne J. Appolonius’s 10th problem. In: Arvo J, editor. Graphics Gems II. New York: Academic Press; 1991. p. 19–24. [30] Shamos MI, Hoey D. Closest-point problems. Proceedings of the 16th annual IEEE symposium on foundations of computer science 1975 p. 151–62. [31] Sugihara K. A simple method for avoiding numerical errors and degeneracy in Voronoi diagram construction. IEICE Transactions. Fundamentals 1992;E75-A:468–77. [32] Sugihara K, Iri M. Construction of the Voronoi diagram for ‘one million’ generators in single precision arithmetic. Proc IEEE 1992;80: 1471–84. [33] Sugihara K, Sawai M, Sano H, Kim D-S, Kim D. Disk packing for the estimation of the size of a wire bundle. Jpn J Ind Appl Math 2004; 21(3):259–78.
[34] Thiessen AH. Precipitation averages for large areas. Monthly Weather Rev 1911;39:1082–4. [35] Voloshin VP, Beaufils S, Medvedev NN. Void space analysis of the structure of liquids. J Mol Liq 2002;96–97:101–12. [36] Will H.-M. Computation of Additively Weighted Voronoi Cells for Applications in Molecular Biology. Ph.D. Thesis, ETH, Zurich; 1999. [37] Yap C. Robust geometric computation. In: Goodman JE, O’Rourke J, editors. Handbook of discrete and computational geometry. New York: CRC Press LLC; 1997. [38] Computational Geometry Algorithms Library (CGAL) Homepage. http://www.cgal.org/. [39] The CORE Homepage. http://www.cs.nyu.edu/exact/core/. [40] RCSB Protein Data Bank (PDB) Homepage. http://www.rcsb.org/pdb/.
Euclidean Voronoi diagram of 3D balls and its computation via tracing edges Deok-Soo Kima,b,*, Youngsong Chob, Donguk Kimb a
Department of Industrial Engineering, Hanyang University, 17 Haengdang-dong, Seongdong-gu, Seoul 133-791, South Korea b Voronoi Diagram Research Center, Hanyang University, 17 Haengdang-dong, Seongdong-gu, Seoul 133-791, South Korea Received 20 November 2004; accepted 5 February 2005
Abstract Despite its important applications in various disciplines in science and engineering, the Euclidean Voronoi diagram for spheres, also known as an additively weighted Voronoi diagram, in 3D space has not been studied as much as it deserves. In this paper, we present an algorithm to compute the Euclidean Voronoi diagram for 3D spheres with different radii. The presented algorithm follows Voronoi edges one by one until the construction is completed in O(mn) time in the worst-case, where m is the number of edges in the Voronoi diagram and n is the number of spherical balls. As building blocks, we show that Voronoi edges are conics that can be precisely represented as rational quadratic Be´zier curves. We also discuss how to conveniently represent and process Voronoi faces which are hyperboloids of two sheets. q 2005 Elsevier Ltd. All rights reserved. Keywords: Euclidean Voronoi diagrams; Edge-tracing; Rational quadratic Be´zier; Apollonius problem
1. Introduction Since its initial definition and characterization by the Russian mathematician G. Voronoı¨ in 1907 and the first computational algorithm presented by Shamos and Hoey [30], the Voronoi diagram has been known as one of the central topics in many disciplines in science and engineering including computational geometry. Due to its natural descriptive and manipulative capability, the Voronoi diagram and its variations have been known as various names such as a Thiessen polygon [34], a medial axis transformation (MAT) [19], a symmetric axis transformation (SAT) [4,5], a skeleton [18], a proximity map [12], a Dirichlit tessellation, a thinning, etc. Even long time before Voronoı¨, an idea quite close to what is now known as a Voronoi diagram was used by Descartes to * Corresponding author. Address: Department of Industrial Engineering, Hanyang University, 17 Haengdang-dong, Seongdong-gu, Seoul 133-791, South Korea. Tel.: C82 2 2220 0472; fax: C82 2 2292 0472. E-mail addresses: [email protected] (D.-S. Kim), [email protected] (Y. Cho), [email protected] (D. Kim).
0010-4485//$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.cad.2005.02.013
describe the distribution of matters in the solar system and its environs [24]. We believe that Aurenhammer [3] and Okabe et al. [24] might be the best source of information for the comprehensive introduction to various Voronoi diagrams, their properties, algorithms and related problems. The ordinary Voronoi diagram for point set has been studied extensively and its properties are well-known in 2 and higher dimensions. However, Voronoi diagrams for circles in 2D and spheres in 3 or higher dimensions in the Euclidean distance measure, which we call Euclidean Voronoi diagrams, have not been explored as much as they deserve even though they have significant potential impacts on diverse applications in both science and engineering [1,11,22,28,35]. We want to note here that a Euclidean Voronoi diagram of circles or spheres is usually referred to as an additively weighted Voronoi diagram in the computational geometry community [24]. However, in this paper we will use the term a Euclidean Voronoi diagram of circles or spheres to emphasize its geometric interpretations and immediate practical applications in various scientific and engineering disciplines. We believe that a Euclidean Voronoi diagram for spheres in 3D can have more significant applications in various fields than its planar counterpart since most important geometric problems are 3D in nature. For example,
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
Fig. 1. Examples of the Euclidean Voronoi diagram of circles: (a) random circles, and (b) circles on a spiral.
the structural analysis of proteins or RNA’s requires an efficient computational tool to analyze the spatial properties among atoms [11,16,28]. In the design of new materials, a similar analysis is fundamental as well [15,23]. However, due to the lack of appropriate algorithms and stable running codes for Euclidean Voronoi diagrams of circles and spheres, most applications instead adapt an ordinary Voronoi diagram of points or a power diagram for the circles or spheres. For Euclidean Voronoi diagrams of circles in a plane, a simple yet fast and robust algorithm is available [13,14]. See Fig. 1 for the examples of computed output for circle sets. In Kim et al. [13,14], a Euclidean Voronoi diagram for circles is computed via an edge-flipping procedure, which takes O(n2) time in the worst-case where n is the number of circles, using the ordinary Voronoi diagram of center points as an initial starting diagram [13,14]. Note that this diagram has been used in various applications such as a cable packing problem [33]. Since the edge-flipping in 3D, however, lacks one-to-one correspondence between before and after flipping [24], a direct application of the edge-flipping concept may not be appropriate for the Euclidean Voronoi diagram for spheres in 3 or higher dimensions unless special cares are taken. Therefore, developing a good practical algorithm for a Euclidean Voronoi diagram of spheres in 3D is very important for practical reasons. Unlike other types of Voronoi diagrams or power diagrams, there are only very few previous works on the Euclidean Voronoi diagram of spheres. Aurenhammer discussed the possibility of transforming the computation of the Euclidean Voronoi diagram of spheres in d-dimension to one of (dC1)-dimensional power diagram construction by the intersections of (dC2)-dimensional half-spaces [2]. However, neither an algorithm nor an implementation is given for dR3. Gavrilova, in her Ph.D. Thesis, reported important properties of Euclidean Voronoi diagrams in arbitrary dimemsions, including shapes of the Voronoi regions, nearest neighbors and empty-sphere properties [8,9]. Will wrote a comprehensive Ph.D. thesis dedicated to the computation of the Voronoi regions for the Euclidean Voronoi diagram of spheres [36]. In his thesis, Will showed that a Voronoi region can be computed in Q(n2)
1413
Fig. 2. Examples of EVD(S): (a) an a-helix consisting of 67 atoms, and (b) PDB ID 1BH8 consisting of 1074 atoms (680 C’s, 181 N’s, 203 O’s, and 10 S’s).
time and actually implemented the algorithm so that the complete Voronoi diagram can be computed by constructing one Voronoi region after another. Luchnikov et al. proposed a practical idea of tracing edges which is simple yet powerful method for obtaining a Voronoi diagram of spheres in 3D [21]. However, a correct and efficient implementation of their algorithm is not so easy and the details of their algorithm are not available. Recently, we reported on the full implementation of an edge-tracing algorithm for constructing the Voronoi diagram with discussions of various applications including the analysis of protein structures [16,17]. In this paper, we will discuss several properties of the Euclidean Voronoi diagram for spheres in 3D and present the details of an algorithm based on the edge-tracing algorithm which computes the Voronoi diagram by tracing Voronoi edges. Given 3D balls, we describe how to trace the edges correctly and efficiently and discuss the representations of edges and faces. Having fully implemented the algorithm, we provide some experimental results as well. Shown in Fig. 2(a) is a Euclidean Voronoi diagram computed from a set of 67 atoms which describes an a-helix of a real protein. Fig. 2(b) illustrates a Voronoi diagram for a complete protein data with the PDB ID 1BH8 [40] which consists of 1074 atoms (680 C’s, 181 N’s, 203 O’s, and 10 S’s). In both Voronoi diagrams, faces are not intentionally rendered for the purpose of visualization. This paper is organized as follows. In Section 2, we define some terminology and notations used in this paper. Then, we discuss geometric primitives such as vertex positions, edge equations, and face equations of a Voronoi diagram for 3D balls in Section 3. The edge-tracing construction for the topology of the Voronoi diagram is described in Section 4. In Section 5, we conclude this paper with the discussions on the experimental performance of the proposed algorithm and some future works.
2. Terminology and definitions Let BZ{b1,b2,.,bn} be a set of generators for a Voronoi diagram where bi is a 3-dimensional spherical ball. Hence, biZ(ci,ri) where ciZ(xi,yi,zi) and ri denote the center and
1414
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
the radius of a ball bi, respectively. We assume that no ball is completely contained inside another ball even though intersections are allowed between balls. Associated with each ball bi, there is a corresponding Voronoi region VRi for bi, where VRi Z fpjdistðp; ci ÞK ri % distðp; cj ÞK rj ; i sjg. Then, EVD(B)Z{VR1, VR2,., VRn} is called a Euclidean Voronoi diagram for B. The topology of the Voronoi diagram in this paper is represented as an edge-graph GZ(V,E) where VZ{v1,v2,.} and EZ{e1,e2,.} are sets of Voronoi vertices and edges, respectively. Once G is available, the face set FZ{f1,f2,.} of Voronoi diagram can be easily constructed. In this paper, the ordinary L2-metric is used to define a Euclidean Voronoi diagram. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi In other words, 2 2 2 distðci ; cj ÞZ ðxi K xj Þ C ðyi K yj Þ C ðzi K zj Þ . As in ordinary point set Voronoi diagrams, Voronoi regions corresponding to balls on the boundary of the convex hull of B are unbounded. Other regions are bounded by a set of faces, called Voronoi faces, where a Voronoi face f is defined by two immediately neighboring balls. Note that a Voronoi face is always a connected-subset of hyperboloid of two sheets. A Voronoi face intersects another Voronoi face to form a Voronoi edge e which is a conic section. When Voronoi edges intersect, a Voronoi vertex v is defined. Note that there is always one-to-one correspondence between a Voronoi vertex v and a sphere S tangent to neighboring balls defining the vertex. Hence a tangent sphere is always centered at the corresponding vertex v. The tangent sphere S is called empty since it never intersects or contains any other ball except at the tangent points with its generating balls. In this paper, we assume that the degree of a Voronoi vertex is always four which means that no five balls are cotangent to an empty tangent sphere. We also assume that a Voronoi edge is always the common intersection among three Voronoi faces and therefore only three balls define an edge. Even with these assumptions, topological conflicts may occur due to numerical errors in the course of the Voronoi diagram construction [31,32]. It is known that the exact computation techniques can resolve such topological inconsistencies [31,32,37]. In this paper, however, we rather focus on describing the Voronoi diagram construction assuming that spheres are in general positions and leave the exact computation issue for the future research. Balls, in this paper, refer to spherical generators defining a Voronoi diagram while tangent spheres denote spheres corresponding to Voronoi vertices. For convenience of notation, we will ignore the term Voronoi and call them simply a vertex, edge or face unless further definition or explanation is necessary. In a similar fashion, the term Euclidean will also be ignored. Shown in Fig. 3 is a 2-dimensional analogy of a small Voronoi diagram of 3D balls. An edge e is oriented from a start vertex vs to an end vertex ve. The empty tangent spheres corresponding to vs and ve are called a start tangent sphere Ss and an end tangent sphere Se, respectively. An oriented edge
Fig. 3. Euclidean Voronoi diagram of spheres and related terminology.
e is always defined by three neighboring balls which surround e, and these balls are called gate balls bg1,bg2 and bg3 (In the figure illustrating its 2D analogy, only two gate balls are shown.) They are called gate balls since the triangle defined by the centers of the balls can be viewed as a gate which opens to the edge direction. The other two balls bs and be are called a start ball and an end ball since they uniquely correspond to the start and end vertices, respectively. To compute an ending vertex ve for an edge e, it is necessary to do computations with some balls in the given ball set B. Given three gate balls bg1, bg2 and bg3 and a start ball bs, it turns out that the computation of the tangent sphere Se, corresponding to an end vertex ve, requires to test if Se can be defined among four balls consisting of three gate balls bg1, bg2 and bg3 and each ball in KZB\{bg1,bg2,bg3}. We call K a candidate ball set for an edge e. In other words, in the example, KZ{bs,be,bi,bj}. Note that K contains the start ball bs as well (We will discuss this later). Associated with each edge, an angular distance q is defined as another useful distance measure. Given three gate balls, let bg1 be the smallest. Then, an angular distance q of a point p from a vertex vs is the angle :vscg1p, where cg1 is the center of bg1. Note that 0!q!2p and q is a directed distance where a positive direction is defined by the edge direction from the start vertex vs. When there are a number of angular distances defined by different balls, we call the one giving the minimum angle the closest to the start vertex. 3. Geometric primitives In the computation of a Voronoi diagram of circles in a plane, the computations of geometric primitives such as the coordinates of the vertices and the equations of the edges have been separated from the construction of the topology between vertices and edges [13,14]. In the current problem for spheres in 3D, we take a similar approach. In 3D, the geometry computation consists of vertices, edges, and faces while its 2D counterpart consists of vertices and edges. We want to mention here that the computation of vertices is inevitable in the construction of the topology. However, the computation of equations for edges and faces
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
can be avoided in the construction of a valid topology of a Voronoi diagram. Once a topology is constructed with the coordinates of vertices, they can be separately computed without much effort when it is necessary. 3.1. Voronoi vertices From the definition of a Voronoi diagram, a Voronoi vertex is the center of an empty sphere tangent to four nearby balls. Gavrilova and Rokne [9] recently reported an elegant algorithm to compute the sphere(s) tangent to the four balls in a general setting and we briefly summarize it here as follows. Suppose that b1,.,b4, denote four balls from which a Voronoi vertex, and therefore a corresponding tangent sphere, is to be computed. Note that the center of the tangent sphere is the equi-distant point from the four given balls. Without loss of generality, we can assume that all balls are reduced by the radius of the smallest ball b4. Then, b4 is reduced to a point. Transforming the four-ball system so that b4 coincides with the origin, the tangent sphere S with a center (x,y,z) and a radius r can be computed by solving the following system. Let a transformed shrunken ball bi have a center (xi,yi,zi) and a radius ri, for iZ1,.,4. ðx K x1 Þ2 C ðy K y1 Þ2 C ðz K z1 Þ2 Z ðr C r1 Þ2
(1)
ðx K x2 Þ2 C ðy K y2 Þ2 C ðz K z2 Þ2 Z ðr C r2 Þ2
(2)
ðx K x3 Þ2 C ðy K y3 Þ2 C ðz K z3 Þ2 Z ðr C r3 Þ2
(3)
x2 C y 2 C z2 Z r 2 :
(4)
By reducing the system to a simpler form, the first three equations are reduced to three linear equations in terms of r. In other words, each coordinate, x, y, or z, is expressed by a linear equation of an unknown variable r. Then, substituting these equations into the fourth equation results in a quadratic equation in terms of r. By looking at the rank of the linear system consisting of three linear equations, they showed that the four ball configuration can be categorized into one of four different cases with none, one, two or infinite tangent spheres. For a detailed explanation, the readers is referred to Gavrilova and Rokne [9]. 3.2. Voronoi edges Once the edge-graph GZ(V,E) of the topology of a Voronoi diagram is computed, it is often necessary to compute the actual equations of Voronoi edges. Unless it is a degenerate case, an edge is always defined as a locus of points equi-distant from the surfaces of three balls b1,b2, and b3. Similar to a Voronoi vertex, we shrink the balls by the radius of the smallest ball and transform the ball-system so that the smallest ball shrunken to a point coincides with the origin. Let the transformed shrunken balls have centers
1415
(xi,yi,zi) and radii ri, iZ1,2 and 3. Note that the smallest ball b3 now has a radius of zero and coincides with the origin. Then, the equation of an edge is the solution of the following three equations. ðx K x1 Þ2 C ðy K y1 Þ2 C ðz K z1 Þ2 Z ðr C r1 Þ2
(5)
ðx K x2 Þ2 C ðy K y2 Þ2 C ðz K z2 Þ2 Z ðr C r2 Þ2
(6)
x2 C y2 C z2 Z r 2 ;
(7)
where (x,y,z) is a point on the edge. Rearranging the above system yields the following. x1 x C y1 y C z1 z K ðx21 C y21 C z21 K r12 Þ=2 C r1 r Z 0
(8)
x2 x C y2 y C z2 z K ðx22 C y22 C z22 K r22 Þ=2 C r2 r Z 0
(9)
x2 C y2 C z2 Z r 2
(10)
Two linear equations in Eqs (8) and (9) show that they are planes and the solution of both equations is a line in general. Since the third equation is a sphere, the solution of the whole system is two points on the sphere in general. Hence, the edge equation can be obtained by a variable r. Looking back at the linear equations once more, bothffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 C y2 C z2 planes move at constant velocities of Kr = x 1 1 1 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 and Kr2 = x2 C y2 C z2 in the normal direction of both planes, respectively. This means that the trajectory of the moving intersection line defines another plane. Hence, its intersection with the third equation with a variable r, which is the resulting Voronoi edge, has to be a planar curve. Similarly, two latter equations together define a hyperboloid of two sheets. Therefore, a Voronoi edge is an intersection between a hyperboloid and a plane, and indeed is a conic curve. It can be also differently shown by the fact that a Voronoi edge satisfies the definition of a spine curve of Dupin cyclide [25,26]. Will has also shown that a Voronoi edge is a conic curve and further classified it by looking at the relative configuration among these balls [36]. Even though the type of the curve is known, its mathematical representation is another issue since applications usually require processing of these curves. Hence, a convenient curve representation is necessary for the processing. Since Voronoi edges are conics, they can be exactly represented in a form of a rational quadratic Be´zier curve defined as bðtÞ Z
w0 ð1 K tÞ2 p0 C 2w1 ð1 K tÞtp1 C w2 t2 p2 ; t 2½0; 1; w0 ð1 K tÞ2 C 2w1 ð1 K tÞt C w2 t2 (11)
where p0, p1, and p2 are the control points, and w0, w1, and w2 are the corresponding weights, respectively [7]. It is known that the above equation can be completed once all five of the following parameters are known: two end points, tangent vectors at both end points, and a point through
1416
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
becomes a tangent plane of B12 at v. As a result, a vector starting from v and ending at p has the same angle with vectors starting from v and ending at c1 and c2. In a similar manner, the angle bisecting plane between b2 and b3 becomes a tangent plane at v, and therefore every vector starting from v on the tangent plane has the same angle with vectors starting from v and ending at c2 and c3. Therefore, the intersection between the two tangent planes becomes the tangent line of a Voronoi edge e at v and it is equi-angular with the three vectors starting from v and ending at c1, c2, and c3. , Fig. 4. Proof of Lemma 1.
which the curve passes [7]. In our problem, two end points are already given as Voronoi vertices of the edge. Let b1, b2, and b3 be three gate balls for an edge e where r1%r2%r3, and v be a vertex from which we want to define a tangent vector of the edge. Suppose that ci is the center of a ball bi and ti Z ðci K vÞ=jjci K vjj. Then, the following lemma holds. Theorem 1. (Angle Trisectioning) The tangent vector t at a point v on a Voronoi edge e is equi-angular with three vectors t1,t2 and t3. Proof. An edge e is a common intersection among three bisecting surfaces, i.e. between b1 and b2, b2 and b3, and b3 and b1. Let us consider a bisector B12 between b1 and b2. Let v be an arbitrary point on e, and let a plane P bisect the angle :c1vc2 and p be a point on P other than v as shown in Fig. 4. In addition, we locate a point a on a line segment vc2 such that jjvKajjZjjvKc1jj. Then, jjpKc1jjZjjpKajj holds since Dvpa and Dvpc1 are congruent. Since Dpc2a forms a triangle, the triangular inequality holds, i.e. jjpK c2 jj% jjpK ajjC jjaK c2 jj. Since jjpKajjZjjpKc1jj and jjaKc2jjZr2Kr1, the above inequality can be rearranged as jjpK c2 jjK r2 % jjpK c1 jjK r1 . Therefore, P cannot intersect B12 at any point except at v unless r1Zr2, and the angle bisector P
Due to this lemma, the following corollary immediately holds and is presented without a proof. Corollary 2. The tangent vector t at v satisfies the following system of equations. ti $t Z tj $t Z tk $t;
(12)
where ts0. The last parameter, a passing point p, can be found as follows. Suppose that we define a plane P passing through the centers of three gate balls bg1, bg2 and bg3. Then, the intersection of P with the three balls results in three circles, Cg1, Cg2, and Cg3, on P as shown in Fig. 5(a). Therefore, the passing point p, which lies on this plane, has to be the center of a circle C tangent to these three circles. Hence, computing C is equivalent to solving The Apollonius 10th Problem in a plane P and its solution process is well-described in Kim et al. [14] and Rokne [29]. Once these five parameters are known, the complete rational quadratic Be´zier representation can be found by solving the following equation. Suppose that a curve is represented by the standard form of a rational quadratic Be´zier, i.e. w0Zw2Z1. Then, w1 can be obtained by t1 ; w1 Z pffiffiffiffiffiffiffiffiffi 2 t0 t2
Fig. 5. Passing point on a Voronoi edge solved by The Apollonius Tenth Problem: (a) one tangent circle, and (b) two tangent circles.
(13)
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
where t0, t1, and t2 are the barycentric coordinates of a passing point with respect to the triangle vertices b0, b1, and b2 [7]. Note that the barycentric coordinates are positive if a passing point is located inside a triangle Db0b1b2. In addition, w1 can be negative if the passing point is not inside the triangle. If p lies on a valid edge segment, w1 computed by the above equation can be directly applied. If p lies outside a valid edge segment, however, the correct rational quadratic Be´zier representation of the valid edge segment can be obtained by reversing the sign of w1 calculated by the above equation using p. There could be, however, even worse cases. When an edge is elliptic, there could be two such intersection points, p1 and p2, between P and the edge e as shown in Fig. 5(b). In this case, the situation can be either (i) both are on a valid edge segment, (ii) both are not on a valid edge segment, or (iii) only one of them is on a valid edge segment. The key part of this problem is how to tell which points, p1, p2, or both, are on a valid segment. Suppose we pick an arbitrary point between two points, for example p1. Then, checking if p1 lies on a valid edge segment can be easily done by testing if the angular distance of p1 from the start vertex vs is less than the angular distance of the end vertex ve from vs. Note that ve is already computed. If p1 is on a valid edge segment, the weight w1 computed by Eq. (13) can be directly used. If p1 is not on the valid segment of the edge, the sign of computed weight w1 should be reversed regardless of the computed sign of w1. From the above discussions for both cases, the following lemma can be concluded. Lemma 3. Suppose that the weight w1 is computed by Eq. (13). If a passing point p is on the valid edge segment, w1 can be directly used as the weight. Otherwise, its sign should be reversed.
1417
b Z K4ri2 C 4y2j K 4rj2 K 8yi yj C 8ri rj C 4y2i ; c Z 4z2i K 4rj2 K 4ri2 C 8ri rj C 4z2j K 8zi zj ; d Z K8xi yj K 8yi xj C 8xj yj C 8xi yi ; e Z 8yj zj C 8yi zi K 8yi zj K 8zi yj ; f Z 8xj zj K 8xi zj K 8zi xj C 8xi zi ; g Z K8xi ri rj K 8xj ri rj C 4xi z2j K 4x3j K 4x3i C 4xi rj2 C 4xi ri2 K 4xj z2j K 4xi y2i K 4xi z2i C 4xj ri2 C 4xi x2j C 4xi y2j K 4xj y2j C 4y2i xj C 4x2i xj C 4z2i xj C 4xj rj2 ; h Z K8yi ri rj K 8yj ri rj K 4y3i K 4y3j C 4y2i yj C 4yj ri2 C 4yj rj2 C 4yi z2j C 4yi y2j K 4x2j yj K 4x2i yi C 4yi rj2 K 4yj z2j C 4z2i yj C 4yi ri2 K 4yi z2i C 4x2i yj C 4yi x2j ; k Z 4zj ri2 C 4y2i zj C 4zi rj2 K 4y2i zi K 4z3i C 4zj rj2 K 4z3j K 4y2j zj K 8zi ri rj K 8zj ri rj C 4z2i zj C 4zi ri2 K 4x2i zi C 4x2i zj K 4x2j zj C 4zi x2j C 4zi z2j C 4zi y2j ; and l Z 2x2i y2i C 2x2i z2i K 2x2i x2j K 2x2i y2j K 2x2i z2j K 2x2i ri2 K 2x2i rj2 C 2y2i z2i K 2y2i x2j K 2y2i y2j K 2y2i z2j K 2y2i ri2 K 2y2i rj2 K 2z2i x2j K 2z2i y2j K 2z2i z2j K 2z2i ri2 K 2z2i rj2 C 2x2j y2j C 2x2j z2j C 4x2i ri rj C x4i C y4i C z4i C x4j C y4j C z4j C ri4 C rj4 C 4y2i ri rj C 4z2i ri rj C 4x2j ri rj C 4y2j ri rj
3.3. Voronoi faces
C 4z2j ri rj K 2x2j ri2 K 2x2j rj2 C 2y2j z2j K 2y2j ri2 K 2y2j rj2 Once the edge-graph GZ(V,E) is known, it is easy to tell which pair of balls define a Voronoi face. Two topologically neighboring balls bi and bj define a Voronoi face as jjp K ci jj K ri Z jjp K cj jj K rj
(14)
where ci and ri are the center and radius of ball bi. Then, the above equation defines a hyperboloid of two sheets. Rearranging the above equation results in the following implicit equation. ax2 Cby2 Ccz2 Cdxy Ceyz Cfzx Cgx Chy Ckz Cl Z 0; (15) where a Z 4x2j K 4rj2 K 8xi xj C 8ri rj C 4x2i K 4ri2 ;
K 2z2j ri2 K 2z2j rj2 K 4ri3 rj C 6ri2 rj2 K 4ri rj3 :
ð16Þ
Given a face equation, it is usually better to represent the faces in parametric form for reasons that are well-known in the geometric modeling community. Since the Voronoi face in our problem is always a hyperboloid, the natural choice for the parametric representation is a rational Be´zier surface with an appropriate degree and topology. It turns out that the topological shape of a Voronoi face in our problem form a polygon such as a triangle, rectangle, pentagon, etc. In addition, it can even have holes inside. Therefore, the Be´zier form cannot be used directly but needs approprite modifications. Hence, there are three alternative ways to represent a Voronoi face as a rational Be´zier form. If a rectangular topology is used for rational Be´zier surface representation, it is in general necessary to trim
1418
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
the surface for a non-rectangular shaped Voronoi face. It is known that trimming of a parametric surface takes a significant amount of computing time and memory. Most seriously, there is usually no way to avoid the gaps between two neighboring surfaces if trimming has to be employed. If a Voronoi face is subdivided into a number of triangular segments on a hyperboloid, each triangle can be represented as a rational Be´zier surface having triangular topology. If this approach is used in our problem, trimming is not necessary. However, it may be necessary to implement some geometric operations redundantly. For example, consider the rolling-ball blending surface between balls. It is known that this surface can be conveniently represented using a rational surface having rectangular topology. Hence, having both triangular and rectangular topologies for rational surfaces simultaneously does not seem to be a good way to represent Voronoi faces. If each triangle in the above is again subdivided into three rectangles, all of them can be represented as rational Be´zier surfaces having rectangular topology. While avoiding the problem discussed above, this approach requires a significantly larger amount of memory. While the choice of an appropriate surface representation remains to be made, principal uses of Voronoi diagrams in most applications suggest a different alternative. The main uses of Voronoi faces are the visualization of the faces themselves, the computation of volumes and boundary areas of Voronoi regions, etc. Suppose that two balls b1 and b2 define a Voronoi face where r1%r2. In addition, let pm be the mid point between the centers of b1 and b2. Then, we transform both balls so that pm is located at the origin and the centers of b1 and b2 are placed on the positive and negative Z-axis, respectively, as shown in Fig. 6. We call this position of two balls a canonical position. Let b1 and b2 be centered at (0,0,d) and (0,0,Kd), respectively. Then, the equation for the bisector surface between the balls in the canonical position is
Fig. 6. Voronoi face in the canonical position.
represented as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 C y2 C ðz K dÞ2 K r1 Z x2 C y2 C ðz C dÞ2 K r2 : (17) While the above equation defines a hyperboloid of two sheets, the valid Voronoi face between two balls is always one of either sheet, not both. Rearranging the above equation in an explicit form gives the following equation. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2r ðKr 2 C 4d2 Þð4x2 C 4y2 C 4d 2 K r 2 Þ z ZG ; (18) K4r 2 C 16d 2 where rZr2Kr1. As shown in the above equation, the equation is divided into the positive and negative parts, and the desired sheet of a hyperboloid is the positive part. In Fig. 6, for example, fV is the valid sheet and fI is the invalid one. When a Voronoi face is represented and manipulated in applications, it is always necessary to test the validity of the face. For example, when we have to triangulate a Voronoi face for a visualization purpose, we have to choose one of either points on the hyperboloid computed from the above equation when X and Y coordinates are given. Since we have transformed the Voronoi face in the canonical position, we can easily locate the correct surface by simply looking at the signs of points on the hyperboloid. If Voronoi faces are handled in this way, both evaluating a point on the face and verifying if a point is on the face become much simpler. 3.4. Visualization of Voronoi faces To illustrate how to process Voronoi faces in applications, we show how to render those faces as an example. Shown in Fig. 7 are examples similar to those presented in Luchnikov et al. [21]. Fig. 7(a) shows an example of an atom complex which consists of 15 atoms with identical radii, while Fig. 7(d) shows a similar example with three different radii. Their Voronoi diagrams are computed in Fig. 7(b) and (e), respectively. For convenience of visualization, we have shown only the Voronoi region for the yellow atom in the center while ignoring edges and faces emanating toward infinity. Fig. 7(c) and (f) show the Voronoi regions only. Note that the faces in Fig. 7(c) are planar and the associated edges are line segments like an ordinary point set Voronoi diagram since the sizes of atoms are identical. On the other hand, the faces in Fig. 7(f) are parts of hyperboloids and the edges are conics represented as rational quadratic Be´zier curves. To render a Voronoi region, it is necessary to triangulate the faces bounding the region so that a graphics library such as OpenGL can be used. To triangulate a Voronoi face, we need to sample some points on the valid part of the hyperboloidal surface associated with the Voronoi face. Since a Voronoi face is a subset of a hyperboloid of two
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
Fig. 7. Fifteen balls with three different radii and the corresponding Euclidean Voronoi diagram. Voronoi edges are represented as rational quadratic Be´zier curves while Voronoi faces are represented as implicit surfaces: (a) 15 balls with identical radii, (b) the corresponding Voronoi diagram where faces and edges going to infinity are ignored, (c) the Voronoi region for the yellow ball, (d) 15 balls with three different radii, (e) the corresponding Voronoi diagram, and (f) the Voronoi region for the yellow ball. Note that the polyhedron is the Voronoi region for the yellow ball which has 14 faces with the other balls.
sheets, such a sampling operation is not an easy task. After correctly sampling such an appropriate point set, it is also necessary to do a triangulation and then eliminate extraneous triangles. It turns out that performing these tasks directly in 3D is relatively difficult as well. Instead, we propose a different approach to take advantage of the surface representation presented in the previous subsection. Given a Voronoi face in the world coordinate system, there are two related balls. Suppose that we transform both balls in the canonical position. Then, the Voronoi face between the transformed balls is always a single-valued function w.r.t. X and Y. Note that a Voronoi face is bounded by a number of Voronoi edges represented as rational quadratic Be´zier curves. Therefore, the edges can be also easily transformed to the same canonical position. Let f be a Voronoi face defined by two balls bi and bj in the world coordinate system. Then, transforming bi and bj to the canonical position produces the corresponding Voronoi face fV as well in the canonical position. In the canonical position, the valid points on the correct sheet of hyperboloid for the face can be easily found. Let EZ{ekjkZ1,2,.} be an edge set bounding f in the world coordinate system where ek is represented as a rational quadratic Be´zier curve. Suppose now that Ec be the edge set of fV in the canonical position. Then, the projection of Ec as well as fV onto the XY-plane produces a polygon Poly bounded by a number of conics in a plane. Note that the transformation to a canonical position and projection of boundary curves can be done easily by simply processing the control points of edges since Be´zier curves and surfaces are affine invariant.
1419
Suppose that we sample a number of points in Poly bounded by Ec and a number of points from Ec itself. Then, we apply a constrained Delaunay triangulation using the polylines on the boundary as constraints and remove extraneous triangles exterior to the boundary Ec of the polygon. If we lift-up the sample points and the topology among the sample points onto fV, we can get the correct triangulation of Voronoi face fV in the canonical position. Hence, the back-transformation of fV to f produces the desired triangular mesh in 3D so that a graphics library such as OpenGL can be used. For the constrained Delaunay triangulation, we have used CGAL library [38]. Note that the polygon Poly is not necessarily convex. It may even contain some holes inside. Hence, we need to be a little bit careful to sample appropriate points inside Ec. The optimal distribution of sample points itself from both inside and boundary curves can be another issue to be studied further. From our experience, sampling points from uniform grids is quite satisfactory. While the constrained Delaunay triangulation in CGAL is too general for the given problem and therefore takes a significant amount of computation, our triangulation problem is quite a special case. Therefore, devising a specialized algorithm would be beneficial.
4. Topology construction by tracing edges Given the building blocks to compute the geometric part of a Voronoi diagram, it is necessary to devise an algorithm to construct its topology counterpart. The basic idea of the edge-tracing algorithm is quite simple. The algorithm first locates a true Voronoi vertex v0 by computing an empty tangent sphere defined by four appropriate nearby balls. Provided that v0 has been found, four edges e0, e1, e2, and e3 emanating from v0 can be easily identified and pushed into a stack called an Edge-stack. Hence, these edges have v0 as their starting vertices. After popping an edge from the stack, the algorithm computes the end vertex of the edge. Note that a vertex can be found by computing an appropriate empty tangent sphere from each of candidate balls plus three gate balls defining the popped edge. If an appropriate empty tangent sphere is found, the center of the sphere becomes the end vertex of the popped edge. Once the end vertex of currently popped edge is found, it is also possible to compute three more edges emanating from this new vertex. Hence, when an edge is created, its end vertex is used as the start vertex of another three new-born edges. (Some of the new-born edges may turn out to be invalid later on.) Note that these edges are also pushed into the Edge-stack. By repeating this process until the Edge-stack is empty, the computation of a Voronoi diagram of a connected graph is completed. Hence, this process is iterated as many times as the number of edges, which is O(n2) in the worst-case and O(n)
1420
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
on the average, where n is the number of balls. Even though the idea is simple yet powerful, designing a correct and efficient algorithm is not so easy. Hence, we elaborate the details of the algorithm step by step as follows. 4.1. Initialization The initialization process computes a tangent sphere from each combination of four balls and tests if the tangent sphere is empty or not. Since there are O(n4) number of fourball combinations and a tangent sphere from each combination should be tested for emptiness against nK4 other balls, a brute-force approach to find an initial vertex v0 can take as long as O(n5) time in the worst-case. A simple trick can help to avoid such a high worst-case time complexity. Suppose that we add four more fictitious balls to the initial ball set. These fictitious balls are intentionally configured so that there is only one empty tangent sphere among them. Let the Voronoi vertex corresponding to this tangent sphere be a fictitious vertex. Then, these fictitious balls are located in the neighborhood of the input ball set so that the Voronoi vertex among the fictitious balls can still remain as a Voronoi vertex for all balls including the fictitious ones. Then, we can safely start the process from this fictitious initial vertex, which takes O(1) time for the initialization. Tracing the edges from the fictitious initial vertex will eventually lead to a Voronoi vertex which is defined only among the input balls. Then, we choose this vertex as the true initial vertex to start the process to construct the topology of the Voronoi diagram by ignoring the four fictitious balls from the ball set and the computations done until now. We want to mention that the proposed approach of computing an initial vertex by adding four fictitious balls not only reduces the worst-case time complexity but also avoids a potentially difficult problem of dealing with the cases where there is none, one, or two possible tangent spheres corresponding to end vertex(es) for a given ball configuration. 4.2. Tracing edges via the edge-stack When an initial vertex v0 is located, we define four edges e0, e1, e2, and e3 of the Voronoi diagram emanating from v0 and push them into the Edge-stack. Hence, these edges have v0 as their starting vertices. Once the edges are stored in the Edge-stack, we pop one edge after another to trace. This tracing is done by completing the definition of a popped edge by computing its end vertex. Note that the end vertex can be found by computing an appropriate empty tangent sphere from 3 gate balls surrounding the popped edge plus one appropriate ball of the nK3 candidate balls. Suppose that a valid end vertex vi is computed for a popped edge. Then, we again define three more edges ei1,
ei2, and ei3 emanating from vi, and again push them into the Edge-stack. Hence, these new-born edges have vi as their starting vertices as well. However, some of these new-born edges may be determined in a later step to be redundant since their end vertices are already created. This will be explained in Section 4.4. By repeating this process until the Edge-stack is empty, the coordinates of vertices and the topology among vertices and edges can be completely identified. Hence, pushing and popping operations on edges with the Edge-stack lead to the complete edge-graph of a Voronoi diagram for balls. While the computation of vertices, each of which corresponds to a tangent sphere among the four balls, is inevitable in the process of constructing the edge-graph, the computation of edge and face geometries is not necessary in the process of edge-graph computation. 4.3. End vertex computation Computing an end vertex of a popped edge is equivalent to finding an end ball associated with the end vertex. Hence, we compute a tangent sphere from four balls, three defining balls of the popped edge and one candidate ball in the ball set, to find the end vertex. After a tangent sphere is computed, its emptiness is tested against all the other balls except the four balls defining the tangent sphere. Here, the cardinality of the candidate ball set is nK3, not nK4, consisting of all balls except only the three surrounding balls for the edge being traced. In other words, the start ball should be also considered as a candidate ball for the computation of the end vertex. This is because a certain configuration of four balls can define two tangent spheres, not only one, as described earlier. Hence, the next lemma follows without a proof. Lemma 4. For an edge e, and therefore three corresponding gate balls bg1, bg2 and bg3, and an input ball set B, the candidate ball set K is defined as KZB\{bg1,bg2,bg3}. Due to the lemma, we can design an algorithm to compute an end vertex given an edge e with its gate balls and the start ball. The following is a simple-minded first version of such an algorithm. Algorithm END-VERTEX-Q Step 1. For bk2K, where K is the candidate ball set. Step 1.1. Compute a tangent sphere Sk from bk and three gate balls Step 1.2. For bl2K\{bk}, Step 1.2.1. If Sk intersects bl then GOTO Step 1. End-for Step 1.3. Enqueue Sk into Q End-for Step 2. If Q is empty, the end vertex is going infinity. Step 3. Else, find S in Q closest to the start vertex, and use the center of S as the end vertex of the current edge.
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
In the brute-force algorithm presented above, the first For-loop in Step 1 iterates nK3 times to compute tangent spheres, and the nested For-loop in Step 1.2 iterates nK4 times for the emptiness test of each tangent sphere against other balls. Since the output from Step 1 is O(n) of empty spheres, Step 3 takes O(n) time in the worst-case. Therefore, algorithm END-VERTEX-Q takes O(n2) time in the worstcase and therefore the next lemma follows. Lemma 5. Using Algorithm END-VERTEX-Q, the edge tracing algorithm constructs a complete edge-graph of a Voronoi diagram of spheres in O(mn2) time, where m is the number of Voronoi edges and n is the number of balls. The above time complexity, however, can be improved by the following lemma. Lemma 6. Given a start vertex and gate balls for an edge e, the end vertex of e can be located in O(n) time in the worstcase. Proof. (See Fig. 8) Suppose that an edge e is being traced starting from a vertex vs, and a ball bi, together with three
1421
gate balls bg1, bg2 and bg3, successfully defines a tangent sphere Si. However, Si may or may not define a correct end vertex for e depending on the other balls. As shown in Fig. 8(a), Si may not be empty but intersect another ball bj, i!j. In this case, the center vei of Si cannot be the end vertex of the edge at all since the angular distance qi of vei from the start vertex vs is greater than qj associated with Sj (See Fig. 1(a)). Hence, vej, the center of tangent sphere Sj defined by the ball bj can be now considered as a candidate for the end vertex of the edge. In Fig. 8(b), bj does not intersect Si and defines another tangent sphere Sj. In this case, which is similar to the above case, it can be shown that qi!qj. Hence, we can simply ignore bj from further consideration. In Fig. 8(c), bj does not intersect Si and also defines another tangent sphere Sj. In this case, it can be shown that qi>qj. Hence, we now consider vej for candidate end vertex instead of vei. Therefore, a scanning of nK3 candidate balls once will eventually locate the valid end vertex, and this process takes O(n) time in the worst-case as well as average-case. , Based on the above lemma, we can devise a more efficient algorithm as follows. Algorithm END-VERTEX-L Step 1. Compute a tangent sphere S1 from three gate balls and b1 in the candidate ball set K. Step 2. S)S1. Step 3. For bk2K\{b1}, Step 3.1. Compute a tangent sphere Sk from three gate balls and bk Step 3.2. If Sk is closer to the start vertex than S, then S)Sk End-for Step 4. Use the center of S as the end vertex of current edge. Given the three gate balls for a popped edge, Algorithm END-VERTEX-L first computes a tangent sphere Si with an arbitrary candidate ball bi. Then, the algorithm selects another candidate ball bj from the candidate ball set and computes a tangent sphere Sj with bj and three surrounding balls. If bj intersects Si, the current Si is replaced by Sj. If not, we choose the tangent sphere between Si and Sj which is closer to the start vertex of the popped edge in their angular distances. Since all balls in the candidate set is scanned only once, this process runs in O(n) time in the both the worst and average cases and finds the correct end vertex. Hence, based on Lemma 6, the following theorem holds.
Fig. 8. End vertex of an edge.
Theorem 7. Using Algorithm END-VERTEX-L, the edge tracing algorithm constructs a complete edge-graph of a Euclidean Voronoi diagram of spheres in O(mn) time, where m is the number of edges and n is the number of vertices.
1422
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
4.4. Checking if a vertex is already computed
5. Discussions and conclusions Voronoi diagrams for spheres in 3-dimensional space in Euclidean distance have many applications in various disciplines in science and engineering. Even though the Euclidean Voronoi diagram of spheres can be so important, its algorithms have not been studied as much as it deserves. In this paper, we proposed an edge-tracing algorithm to compute the 3-dimensional Eulcidean Voronoi diagram for spheres which runs in O(mn) time in the worst-case where m is the number of Voronoi edges and n is the number of given spheres. As building blocks, we have shown that Voronoi edges are conics and can be represented as rational quadratic Be´zier curves. We also discussed how to conveniently represent and process Voronoi faces which are parts of hyperboloids of two sheets. Fig. 9 illustrates the run-time behavior of the current implementation. The X- and Y-axis are the number of atoms in the example file and the computation time in seconds, respectively. The experiments were performed on a Pentium IV with the clock speed of 3.0 GHz and 2 GB memory. In the figure, the lower curves denote accelerations using
computation time (sec.)
Suppose that a new vertex for an edge is obtained by tracing an edge. If the vertex has not previously been computed, then we can safely use the new vertex to complete the edge as an end vertex. However, if it already exists, we have to be careful with the new vertex since it is redundant. Instead of creating the new vertex for the end vertex of the edge being traced, we rather search the existing one to complete the edge definition since the existing one already has associated topology information partially determined. Whenever a new vertex is computed, therefore, it is necessary to check if this vertex is already created or not. For the search of a vertex, we have devised a table of computed vertices and named it the Vertex Index Dictionary (VIDIC). Each entry, which corresponds to a vertex, in the dictionary consists of indices to four balls defining the vertex and a pointer to the definition of the vertex. Note that the number of entries in the dictionary is bounded by the number of vertices of a Voronoi diagram and therefore is O(n2) in the worst-case. However, we believe that a tighter worst-case bound for the size of VIDIC much lower than O(n2) exists in the worst-case as well as average-case since a vertex can be permanently removed from VIDIC once it is connected with four edges. While a simple linear search can be used on VIDIC, we prefer a binary search which takes O(log n) time in the worst-case, which is possible by maintaining the sorted list of the entries. In addition, we want to mention here that hashing, which takes O(1) time on the average, is also possible on VIDIC.
35 30 25 20 15 10
5 0 0
200
400
600
800
1000
1200
the number of spheres original
plane filter
bucket
bucket + plane filter
Fig. 9. Computation time statistics for various models with different sizes.
geometric filters and show significant performance enhancements. We expect that further increase in speed can be achieved by optimizing the codes. While the numbers of vertices, edges, and faces are O(n2) in the worst-case, their average numbers are known as linear to the number of balls. In addition, if we devise an appropriate bucket, we expect the quadratic growth of runtime behavior can be reduced to linear since computing an end vertex can take only O(1) time on the average if the balls are distributed accordingly. For this claim to be mathematically justified, however, the distribution of balls should be carefully investigated and it leaves a challenge. While an ordinary atomic structure yields a Voronoi diagram which is usually a singly connected edge-graph, it is possible that the valid edge-graph of a Voronoi diagram be disconnected for the given input balls. Let two balls bi and bj be non-intersecting and closely located to define a Voronoi face fij as shown in Fig. 10. Suppose that a third ball bk, much smaller than the other two, be located in-between bi and bj so that all three balls altogether define a single elliptic edge e, and the other balls are located far away from the edge e so that no vertex can be defined on e. Hence the edge e has no explicit vertex on it. Then, bi and bk define a face fik, and bj and bk define another face fjk. In this case, fik and fjk intersect at the elliptic edge e. Note that the edge e now defines an island in the face fij and is disconnected from the other part of edge-graph. Hence, in this case, there
Fig. 10. Big brothers problem for a disconnected case.
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
is no edge connecting the small ball to the other balls outside. If this case happens in a part of larger Voronoi diagram, the edge-graph of the whole Voronoi diagram is not a single edge-graph but disconnected. In this case, the isolated edge is either circular or elliptic. If the center point of the small ball is located on the line passing through the centers of two big-brothers, then the edge is circular. If the center of small ball is a little bit off the line, then the edge becomes now elliptic. In this case, special care should be taken to handle the situation to come up with a complete Voronoi diagram. Current work leaves many points of future research issues regarding on the Voronoi diagram construction. The acceleration of running performance of the proposed algorithm can be achieved by employing some geometric filters so that a search space can be significantly reduced. Another improvement would be to devise an algorithm with a better worst-case time complexity. Besides, developing an easy-to-implement algorithm for the Voronoi diagram would be plausible. The most important issue for the construction of a Voronoi diagram for a large model is to guarantee the robustness of the algorithm. In the current research, however, we have mainly focused on the computation of Voronoi diagram assuming that the generators are in general position and we leave the robustness issue for the immediate future research topic. To construct the Voronoi diagram robustly, the exact computation technique should be incorporated into the algorithm. This can be achieved either by using an exact computation package such as CORE library [39] or devising an exact computation scheme specific to our problem. The former approach is easy to come up with a robust code while the latter approach is more efficient in the run-time. Based on this Voronoi diagram, we expect that several important problems in biological systems such as molecular interaction can be solved rather easily and efficiently. In particular, the computation of interface between proteins or protein vs. small molecules can be available with Voronoi diagrams. Moreover, It is noteworthy that all these problems can be solved in a single framework, the Euclidean Voronoi diagram for spheres.
Acknowledgements This research was supported by Creative Research Initiative grant from the Ministry of Science and Technology, Korea. Authors thank Prof. Chee Yap and Prof. Kokichi Sugihara for suggesting the possibility of linear time reduction for the end ball search. Authors also would like to thank Dr Sang-Min Park for her valuable comments and careful review.
1423
References [1] Angelov B, Sadoc J-F, Jullien R, Soyer A, Mornon J-P, Chomilier J. Nonatomic solvent-driven Voronoi tessellation of proteins: an open tool to analyze protein folds. Proteins: Struct Funct Genet 2002;49(4): 446–56. [2] Aurenhammer F. Power diagrams: properties, algorithms and applications. SIAM J Comput 1987;16:78–96. [3] Aurenhammer F. Voronoi diagrams—a survey of a fundamental geometric data structure. ACM Comput Surv 1991;23(3):345–405. [4] Blum H. A transformation for extracting new descriptors of shape. Proceedings of symposium models for perception of speech & visual form 1967 p. 362–80. [5] Blum H, Nagel RN. Shape description using weighted symmetric axis features. Pattern Recogn 1978;10:167–80. [7] Farin G. Curves and surfaces for computer-aided geometric design: a practical guide. 4th ed. San Diego: Academic Press; 1996. [8] Gavrilova M. Proximity and Applications in General Metrics. Ph.D. Thesis: The University of Calgary, Dept. of Computer Science, Calgary, AB, Canada; 1998. [9] Gavrilova M, Rokne J. Updating the topology of the dynamic Voronoi diagram for spheres in Euclidean d-dimensional space. Comput Aided Geometric Des 2003;20(4):231–42. [11] Goede A, Preissner R, Fro¨mmel C. Voronoi cell: new method for allocation of space among atoms: elimination of avoidable errors in calculation of atomic volume and density. J Comput Chem 1997; 18(9):1113–23. [12] Held M. On the computational geometry of pocket machining. Berlin: Springer; 1991. [13] Kim D-S, Kim D, Sugihara K. Voronoi diagram of a circle set from Voronoi diagram of a point set: I. Topology. Comput Aided Geometric Des 2001;18(6):541–62. [14] Kim D-S, Kim D, Sugihara K. Voronoi diagram of a circle set from Voronoi diagram of a point set: II. Geometry. Comput Aided Geometric Des 2001;18(6):563–85. [15] Kim D-S, Chung Y-C, Kim JJ, Kim D, Yu K. Voronoi diagram as an analysis tool for spatial properties for ceramics. J Ceram Process Res 2002;3(3):150–2. [16] Kim D-S, Cho Y, Kim D, Cho C-H. Protein structure analysis using Euclidean Voronoi diagram of atoms. Proceedings of international workshop on biometric technologies (BT 2004). Special forum on modeling and simulation in biometric technology 2004 p. 125–9. [17] Kim D-S, Cho Y, Kim D. Edge-tracing algorithm for Euclidean Voronoi diagram of 3D spheres. Proceedings of 16th canadian conference on computational geometry 2004 p. 176–9. [18] Kirkpatrick DG. Efficient computation of continuous skeletons. Proceedings of 14th IEEE symposium foundations of computer science 1979 p. 18–27. [19] Lee DT. Medial axis transformation of a planar shape. IEEE Trans Pattern Anal Machine Intell 1982;4:363–9. [21] Luchnikov VA, Medvedev NN, Oger L, Troadec J-P. VoronoiDelaunay analyzis of voids in systems of nonspherical particles. Phys Rev E 1999;59(6):7205–12. [22] Montoro JCG, Abascal JLF. The Voronoi polyhedra as tools for structure determination in simple disordered systems. J Phys Chem 1993;97(16):4211–5. [23] Ohno K, Esfarjani K, Kawazoe Y. Computational materials science from ab initio to Monte Carlo methods. Berlin: Springer; 1999. [24] Okabe A, Boots B, Sugihara K, Chiu SN. Spatial tessellations: concepts and applications of Voronoi diagrams. 2nd ed. Chichester: Wiley; 1999. [25] Paluszny M, Boehm W. General cyclides. Comput Aided Geometric Des 1998;15(7):699–710. [26] Pottman H. personal communication. 2004. [28] Richards FM. The interpretation of protein structures: total volume, group volume distributions and packing density. J Mol Biol 1974;82: 1–14.
1424
D.-S. Kim et al. / Computer-Aided Design 37 (2005) 1412–1424
[29] Rokne J. Appolonius’s 10th problem. In: Arvo J, editor. Graphics Gems II. New York: Academic Press; 1991. p. 19–24. [30] Shamos MI, Hoey D. Closest-point problems. Proceedings of the 16th annual IEEE symposium on foundations of computer science 1975 p. 151–62. [31] Sugihara K. A simple method for avoiding numerical errors and degeneracy in Voronoi diagram construction. IEICE Transactions. Fundamentals 1992;E75-A:468–77. [32] Sugihara K, Iri M. Construction of the Voronoi diagram for ‘one million’ generators in single precision arithmetic. Proc IEEE 1992;80: 1471–84. [33] Sugihara K, Sawai M, Sano H, Kim D-S, Kim D. Disk packing for the estimation of the size of a wire bundle. Jpn J Ind Appl Math 2004; 21(3):259–78.
[34] Thiessen AH. Precipitation averages for large areas. Monthly Weather Rev 1911;39:1082–4. [35] Voloshin VP, Beaufils S, Medvedev NN. Void space analysis of the structure of liquids. J Mol Liq 2002;96–97:101–12. [36] Will H.-M. Computation of Additively Weighted Voronoi Cells for Applications in Molecular Biology. Ph.D. Thesis, ETH, Zurich; 1999. [37] Yap C. Robust geometric computation. In: Goodman JE, O’Rourke J, editors. Handbook of discrete and computational geometry. New York: CRC Press LLC; 1997. [38] Computational Geometry Algorithms Library (CGAL) Homepage. http://www.cgal.org/. [39] The CORE Homepage. http://www.cs.nyu.edu/exact/core/. [40] RCSB Protein Data Bank (PDB) Homepage. http://www.rcsb.org/pdb/.
