Region-expansion for the Voronoi diagram of 3D spheres

Computer-Aided Design 38 (2006) 417–430 www.elsevier.com/locate/cad

Region-expansion for the Voronoi diagram of 3D spheres Donguk Kim a, Deok-Soo Kim b,* a b

Voronoi Diagram Research Center, Hanyang University, 17 Haengdang-dong, Seongdong-gu, Seoul 133–791, South Korea Department of Industrial Engineering, Hanyang University, 17 Haengdang-dong, Seongdong-gu, Seoul 133–791, South Korea Received 18 August 2005; accepted 22 November 2005

Abstract Given a set of spheres in 3D, constructing its Voronoi diagram in Euclidean distance metric is not easy at all even though many mathematical properties of its structure are known. This Voronoi diagram has been known for many important applications from science and engineering. In this paper, we characterize the Voronoi diagram of spheres in three-dimensional Euclidean space, which is also known as an additively weighted Voronoi diagram, and propose an algorithm to construct the diagram. Starting with the ordinary Voronoi diagram of the centers of the spheres, the proposed region-expansion algorithm constructs the desired diagram by expanding the Voronoi region of each sphere, one after another. We also show that the whole Voronoi diagram of n spheres can be constructed in O(n3) time in the worst case. q 2006 Elsevier Ltd. All rights reserved. Keywords: Voronoi diagrams; Region-expansion; Proximity

1. Introduction Since its initial introduction in the early 20th century, the Voronoi diagram has been used in many disciplines of science and engineering due to its natural descriptive and manipulative capability. The ordinary Voronoi diagram for a point set and its construction have been studied extensively and the properties are well-known in two and higher dimensions [25]. However, the construction of the Voronoi diagram for spheres in the Euclidean distance metric, often referred to as an additively weighted Voronoi diagram [25], has not been explored sufficiently even though it has significant potential impacts on diverse applications in both science and engineering [1,9,23,27,29]. For example, the structural analysis of protein or RNA requires an efficient computational tool to analyze spatial structures among atoms [9,14,15,27]. In the design of new material and the analysis of its properties, a similar analysis is fundamental as well [21,24,28]. However, due to the lack of appropriate algorithms and stable running codes for the Voronoi diagrams of spheres, most applications have instead

* Corresponding author. Tel.: C82 2 2220 0472; fax: C82 2 2292 0472. E-mail addresses: [email protected] (D. Kim), dskim@ hanyang.ac.kr (D.-S. Kim).

0010-4485//$ - see front matter q 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.cad.2005.11.007

adapted an ordinary Voronoi diagram of points, a power diagram, or an a-hull. Unlike other diagrams, few reports are available on the construction of the Voronoi diagram of spheres. Aurenhammer discussed the transformation of the computation of the Voronoi diagram of d-dimensional spheres to that of (dC1)-dimensional power diagram obtained from the convex hull in (dC2)dimension [2]. Will [30] wrote a comprehensive PhD thesis dedicated to the computation of Voronoi regions in the Voronoi diagram of three-dimensional spheres. Will showed that the Voronoi region of a sphere has a Q(n2) combinatorial complexity and proposed the lower envelope algorithm, which takes an O(n2 log n) expected time for a single Voronoi region, where n is the number of spheres. By implementing the proposals by both himself and Aurenhammer, Will also provided experimental results on various data sets. Gavrilova, in her PhD thesis, reported several properties of the Voronoi diagram for spheres in arbitrary dimensions, including the shape of Voronoi regions, the nearest neighbors and the emptysphere property [7,8]. Luchnikov et al. [22] proposed a practical idea of tracing edges, which is a simple yet powerful way to obtain the desired diagram. Recently, Kim et al. [12–15] reported the detailed characterization of the edge-tracing algorithm and its full implementation for constructing the whole Voronoi diagram with the discussions on various applications including the analysis of protein structures. They showed that the whole Voronoi diagram can be constructed in O(n3) time in the

418

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

worst-case. Boissonnat and Karavelas [3] reported an algorithm to compute a Voronoi region using the convex hull of spheres transformed by inversion. The convex hull topology, then, provides the topology of region boundary. Since the computation of the convex hull takes O(n2) time in the worstcase, their algorithm also takes O(n3) time for the whole Voronoi diagram for n spheres in dimension three. In this paper, we will discuss several properties of the Voronoi diagram for three-dimensional spheres and present a new, simpler yet more powerful, algorithm for the construction of the diagram. The proposed region-expansion algorithm first constructs the Voronoi diagram of the centers of spheres. Considering the centers as shrunken spheres, the algorithm selects each sphere, one after another, and grows to its full size. While a sphere grows in size, the Voronoi region corresponding to the sphere expands simultaneously. Hence, the regionexpansion algorithm for three-dimensional spheres extends its precursor in the plane [16,17]. Repeating this process for all spheres, the correct topology of the desired Voronoi diagram can be obtained. In this research, the Voronoi diagram is represented as a variation of the radial-edge data structure since the diagram is a cell structure, one of the common nonmanifold models, bounded by faces where a face may have an arbitrary number of vertices and it can even have some topological holes inside. Details of the representation issues for the topology of Voronoi diagram are discussed in [4]. Extensive studies have been conducted on the topology of non-manifold models [5,6,18–20]. It is worth mentioning that the difficult edge-disconnectedness in the edge-graph of the Voronoi diagram is not a special case in the proposed algorithm. We also show in this paper that the whole Voronoi diagram of n spheres can be constructed in KO(n3) in the worst-case. This paper is organized as follows. In Section 2, some definitions and terminology are provided. In Section 3, overall scheme of the region-expansion is discussed. In Sections 4–6, we present the classification, detection, and time of events, respectively. In Section 7, we describe how to handle the topological structures associated with infinity. After discussing the disconnected case in Section 8, we provide two algorithms and their time complexity in Section 9. In Section 10, (b)

(a)

experimental results are provided and then we give concluding remarks in Section 11. 2. Terminology and definitions Let SZ{s1,s2,.,sn} be a set of spheres siZ(ci,ri), where ciZ(xi,yi,zi) and ri denote the center and the radius of si, respectively. We assume that no sphere is completely contained inside another although intersections between spheres are allowed. Associated with each sphere si, there is a corresponding Voronoi region VR(si) for si, where VR(si)Z {pjd(p,ci)Kri%d(p,cj)Krj,isj} and d($) denotes L2 Euclidean distance between two points. Then, VD(S)Z{VR(s 1), VR(s2),.,VR(sn)} is called the Voronoi diagram for S. We call si a generator of VR(si). Note that the topology of the Voronoi diagram in this paper can be represented as GZ (V,E,F) where V, E and F denote sets of Voronoi vertices, Voronoi edges and Voronoi faces, respectively. Suppose that we choose a generator which will expand to its full size starting from a point at its center. Then, the generator and the corresponding Voronoi region being expanded are called an expanding generator and an expanding region, respectively. Voronoi vertices on the boundary of an expanding region are called on-vertices Von, and the others are called offvertices Voff. The Voronoi edges on the expanding region are called on-edges Eon, and the edges which have no on-vertex are called off-edges Eoff. The other edges are called radiatingedges Erad, which are further categorized into two groups: (i) edges with one on-vertex and one off-vertex, and (ii) edges with on-vertices at both ends. Similarly, faces are also classified as follows: on-faces Fon, off-faces Foff, and radiating-faces Frad. Let a vertex sphere be an empty sphere simultaneously tangent to the generators defining a Voronoi vertex, and therefore its center is identical to the Voronoi vertex. We say that a sphere is empty if its interior does not intersect any generator. In this paper, we assume that the spheres are in general positions so that the degree of a Voronoi vertex is four and a Voronoi edge is the common intersection among three Voronoi faces. For the convenience of notation, we will omit the term Voronoi unless it is necessary. Similarly, the use of the term Euclidean will also be minimized.

(c)

3. Region-expansion process vs. events e

(e)

(d)

e

(f)

e′

Fig. 1. Region-expansion process in 2D.

Suppose that we are given with an ordinary Voronoi diagram for the centers of spheres. Then, each sphere is associated with a polyhedral Voronoi region, which is an incorrect Voronoi region for the sphere itself. Given a set of spheres, we view as each sphere grows, or expands, to its full size starting from a point at its center. While the expanding generator grows, the corresponding region expands simultaneously. If we can keep the topology among vertices, edges, faces, and regions for the intermediate Voronoi diagram correct and consistent, the complete Voronoi diagram can be computed by repeating the expansion process, generator by generator. If it is necessary, the related geometry can also be appropriately

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

updated in the expansion process. Hence, we call this process region-expansion. In Fig. 1, an example of the region-expansion process in a plane is illustrated. Fig. 1(a) shows seven circle generators where the radii of four of the seven circles are all zeroes, and the Voronoi diagram of all center points is shown in Fig. 1(b). Now, we choose a generator, shown as the shaded circle in Fig. 1(c), and continuously grow it from a point at its center. Then, its Voronoi region, whose boundary is depicted by the solid curves in Fig. 1(c), expands as the radius of the generator increases. Fig. 1(d) shows the resulting Voronoi diagram after completing region-expansion of this generator, and therefore its Voronoi region corresponds to the full-grown generator with respect to the center points in the neighborhood. Given this augmented Voronoi diagram, we perform the region-expansion of another generator, for example, the left most circle in the generator set. Fig. 1(e) shows the correct intermediate Voronoi diagram after completing the regionexpansion of the second generator. Repeating the regionexpansion process for all generators will eventually produce the desired Voronoi diagram as shown in Fig. 1(f). It is obvious that growing a generator, regardless of its dimension, always expands its Voronoi region and therefore increases the area or volume of the region. Based on this observation, it can be easily shown that each on-vertex, during the region-expansion, moves away from the initial expanding region by tracing the radiating-edge attached to the vertex. Similarly, each on-edge moves away from the initial expanding region by tracing a corresponding radiating-face. Note that a sufficiently small growth of an expanding generator leaves the combinatorial structure of the diagram unchanged and causes changes only in the geometries of vertices, edges and faces associated with the boundary of the expanding region. Fig. 1(c) shows such a case: two boundaries of the expanding generator have identical topological structures with slight differences in the geometries of the vertices and edges in the boundary. However, certain changes may occur in the topological structure at some point of time in the expansion process. For example, suppose that an on-vertex moves along a radiatingedge to meet the corresponding off-vertex on the edge. Then, the radiating-edge degenerates to a point and disappears afterwards. We call an event for such a situation which causes changes in the combinatorial structure of the Voronoi diagram during the region-expansion. For example, refer again to the two-dimensional case in Fig. 1. As shown in Fig. 1(c), a sufficiently small increase in the radius of the generator does not affect the topological structure of its region. However, the radius increase makes radiatingedges shorter than before as the region-expansion progress continues. The edge e in Fig. 1(c) is shorter than its counterpart in Fig. 1(b). At some point of time, the edge e will degenerate to a point and then will be flipped to e 0 as shown in Fig. 1(d). Hence, this topology change is called an event. Now, we present observations on the expansion process. Since Voronoi regions are always star-shaped, the on-faces on the expanding region never intersect other faces in the Voronoi diagram except at the boundaries. Suppose that an on-face does

419

intersect an off-face. Then, the Voronoi region, which has both on-face and off-face on its boundary, should have genus one; however, this contradicts the star-shaped property of a Voronoi region. Hence, the next lemma follows. Lemma 1. An on-face intersects another face only at an onedge during the region-expansion process, if it does intersect at all. This lemma states that the intersection between faces never occur at the interior of the faces, but at the boundaries of the faces. Hence, it is sufficient to consider only vertices and edges in the neighborhood to detect new topological changes in the expansion process. Furthermore, it is only necessary to consider edges on radiating-faces since the on-vertices are always constrained to move along the radiating-edges and the on-edges are similarly constrained to stay on the radiatingfaces. Since only the moving on-vertices or moving on-edges cause events, the next simple yet important theorem immediately follows. Theorem 1. When a region is expanding, the upcoming event occurs at vertices or edges on the radiating-faces. Note, however, that any edge except on-edges can be associated with events if the size of the expanding generator is sufficiently large. Since a generator has a prescribed size, only a subset of the events can be realized until the sphere expands to its full size.

4. Classification of events and their actions To devise the region-expansion algorithm based on events, we classify the events into four types and discuss the actions to be taken for each event. 4.1. States of edges As explained earlier, an event occurs only when moving onvertices or on-edges meet other vertices or edges. Before the discussion on the details of events, we first investigate the configurations among vertices and edges and assign a state to an appropriate edge in the configuration. The configurations that an event occurs are classified into the following three cases: † an on-vertex meets an off-vertex, † an on-vertex meets another on-vertex, and † an on-edge meets an off-edge. Note that each configuration accompanies a topological change for an edge. The first configuration corresponds to the case in which a radiating-edge starts to shrink from its initial on-vertex and disappears at its off-vertex. In this case, we assign an end-out state to the radiating-edge. Note that there is only one radiating-edge incident to an on-vertex, and this first configuration applies when the radiating-edge has an off-vertex at its other end. However, a radiating-edge may have two

420

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

on-vertices at its both ends and this case leads to the next configuration. The second configuration is when a radiating-edge has two on-vertices and both moving on-vertices meet at a point in the middle of the radiating-edge. Hence, this radiating-edge disappears afterwards. In this case, we assign a mid-out state to the radiating-edge. The third configuration is the case where a moving on-edge and an off-edge tangentially meet each other at a point in the middle of both edges. In this case, the off-edge is split into two pieces and we assign a split state to the off-edge. Note that the on-edge is also split into two pieces. In the region-expansion process, therefore, different configurations result in different edge states and we claim the next observation. Lemma 2. In the region-expansion process, the changes in the combinatorial structure can be identified according to the states of related Voronoi edges.

4.2. Classification of events Events are defined based on states. Consider the case where two radiating-edges have an identical off-vertex at one of their ends. In such a case, it is possible for both radiating-edges to have end-out states simultaneously and therefore both may disappear at the same time and location. For example, consider the case that two radiating-edges ei and ej have edge-out states simultaneously, and both ei and ej are incident to an identical off-vertex voff. In this case, both ei and ej should simultaneously disappear at the moment when two moving on-vertices corresponding to ei and ej meet at voff. This case maps to a two-end-event. Usually, however, only one radiating-edge having end-out state may disappear at a vertex, which we call a one-end-event. In this case, the changes in the combinatorial structure are different from the changes for the two-end-event case. Therefore, we treat these two cases as different events and they are handled in different ways. We note here that there can be more than two radiatingedges that are incident to an identical off-vertex. However, it can be shown that not more than two radiating-edges that are incident to an identical off-vertex may simultaneously disappear at a certain time in the region-expansion process. Hence, the one-end-event and two-end-event are all possibilities that correspond to the end-out state. On the other hand, a mid-out state and a split state map directly to a mid-event and a split-event, respectively. In the case of a mid-out state, for example, there is no such peculiarity like the end-out state since the topological changes in the related edge only occur at a point in the middle of the edge, and no other edge can share the point. Hence, we call this a midevent since the only topological change in this case is the disappearance of the edge with a mid-out state. Similarly, the topological change due to the split of an edge with a split-state is called a split-event. Hence, the next lemma follows. Lemma 3. In the region-expansion process, there are four event types: an one-end-event, a two-end-event, a mid-event, and a split-event. 4.3. Actions of events Depending on the events, different actions should be taken to maintain the correctness and consistency of the topology. Fig. 2 shows four types of events. The left and right columns show the situations before and after each event, respectively. In the figure, the white and black circles denote off-vertices and on-vertices, respectively. In addition, the symbols in the parentheses denote the edges or faces on the boundary of an expanding region, i.e. on-edges or on-faces.

Fig. 2. Events in the region-expansion: (a) one-end-event, (b) two-end-event, (c) mid-event, and (d) split-event.

4.3.1. One-end-event In the case of one-end-event in Fig. 2(a), the event edge e (a radiating-edge), shown as the thicker edge, disappears after the event is properly processed and three new on-edges are born to form a new on-face f1. In addition, both end vertices of

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

the event edge disappear, and three new on-vertices for the new on-face are born. 4.3.2. Two-end-event In a two-end-event, as shown in Fig. 2(b), the two radiatingedges e 0 and e 00 are removed, which have the end-out states simultaneously and incident to an identical off-vertex. Then, the radiating-face f1 which is bounded by both event edges also disappears. Therefore, three edges constituting the face disappear as well, but a new on-edge is born. Checking whether the end event for an edge is indeed one-end-event or two-end-event can be done simply looking at the event type of next edge in an event queue Q. 4.3.3. Mid-event Fig. 2(c) shows the case of a mid-event. In the figure, e denotes the edge with a mid-out state and e3 denotes an on-edge defining and bounding a radiating-face f1 with e. When this event is processed, the radiating-face f1 disappears, and two onfaces f2 and f3 are merged into one on-face f. Therefore, two edges e10 and e100 are merged into one edge e1. The edges e20 and e200 are similarly handled. Therefore, this event decreases one on-face and two on-edges from the topology under construction. 4.3.4. Split-event A split-event differs from the others in the sense that the event edge does not disappear. As shown in Fig. 2(d), the event edge e does not disappear but is divided into two edges e 0 and e 00 . After processing the current split-event, e 0 and e 00 become new radiating edges so that they will be tested for their event in the future. Similarly, the corresponding on-edge e5 is divided into two edges e50 and e500 . In this case, two new on-vertices, and two new on-edges are born as shown in the figure. In addition, a new on-face f bounded by two new on-edges is also born. Therefore, processing a split-event increases two faces (one on-face and one radiating-face) and two on-vertices. The cardinality change for edges is, however, a little different from the other events. After this event is appropriately handled, three on-edges and two radiating-edges are created, but decreasing one off-edge, which is the current event edge. s3

s1

s2

Fig. 3. Subdivision of the space.

421

5. Detection of events The type of an event can be determined by observing the state of an edge. Note that every edge, except on-edges, corresponds to one event via an edge state. Therefore, every event can be detected by the state of the edges triggering the event. In this section, we will discuss how to identify the states of edges and how to detect the corresponding events from states. 5.1. Edge states via vertex states As explained earlier, we have shown that an edge, except on-edges, has either one of the three types of states: end-out, mid-out, and split state. We found that these edge states can be determined by observing the conditions at the end vertices, which we call vertex states in the following. It can be easily shown that there can be only two situations at a vertex of an edge: the first case is when the edge starts to shrink from the vertex, and the second case is when the edge disappears at the vertex. We assign a ‘C’ state to the vertex of the former case, and a ‘K’ state to the vertex of the latter case. Suppose that an edge e is assigned with an end-out state. Then, e should have one ‘C’ vertex and one ‘K’ vertex since e starts to shrink at one vertex and e will eventually disappear at the other vertex in the region-expansion process. In the case that e is assigned with a mid-event, e starts to shrink from both ends. Therefore, in this case, e consists of two vertices both with ‘C’ states. If e has a split state, both vertices of e should have ‘K’ states. Recall that the vertices in the Voronoi diagram are all of degree four. Since a vertex state of a vertex, as defined in the above, is defined with respect to an edge incident to the vertex, each vertex is associated with four states corresponding to the incident edges. The four states may or may not be identical. The state of a vertex v can be determined as follows. By definition, a vertex sphere SðvÞ corresponds to v and has to be empty for it to be valid. In other words, if the interior of SðvÞ intersects one or more generators, then SðvÞ is not valid and therefore the local topology around v corresponding to SðvÞ is not correct. This may happen when an expanding generator touches the current SðvÞ and then intersects the interior of SðvÞ in the next moment during the region-expansion process. Therefore, if we properly observe how the expanding generator touches SðvÞ and if we look into the edges incident to v, we can correctly determine the states of v with respect to all incident edges. It is important to note that the states of vertices in the neighborhood of an expanding generator completely and correctly determine all events in the region-expansion. Since the region-expansion process always consistently and correctly maintains the topology, we only need to check the relation between SðvÞ and the expanding generator. Let fixed tangent points be the tangent points between SðvÞ and three generators si, sj, and sk which define an edge e incident to a vertex v. Suppose that sl is the fourth generator defining v. Let a reference tangent point be the fourth tangent

422

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

point between SðvÞ and sl. Let a touch point denote the tangent point between SðvÞ and an expanding generator. Let H(e,v) be the hyperplane passing through the three fixed tangent points. Let HC(e,v) denote the open half-space containing the reference tangent point and limited by H(e,v). Then, HK(e,v) denotes the opposite open half-space. Hence, the whole space is divided into three subspaces: HK(e,v), H(e,v), and HC(e,v). Fig. 3 shows a two-dimensional analogy of the previous descriptions for three-dimension. In this figure, the circle centered at v is a vertex sphere. The hyperplane H(e,v) is defined as the line passing through two fixed tangent points on s1 and s2. The reference tangent point, the tangent between the vertex sphere and s3, is located in HC(e,v). If an expanding generator touches SðvÞ at a point in HC(e,v), then five spheres (si, sj, sk, sl, and the expanding generator) are simultaneously tangent to SðvÞ. After the moment of tangency, the expanding generator starts to intrude into SðvÞ to break its emptiness. To keep the emptiness of SðvÞ valid, v starts to move along the edge e. Note that SðvÞ moves and changes its radius according to the moving v, and moving SðvÞ should keep the tangency with si, sj, and sk which define e. In the planar case, shown in Fig. 3, e is defined by the generators s1 and s2. On the other hand, if an expanding generator touches SðvÞ at a point in HK(e,v), v should start to move further away from HK(e,v). However, there is no point on e that v can move since v itself is one of the end vertices of the edge e. Hence, the edge e cannot exist in its current topology structure and has to disappear after the moment of tangency. The above observations lead to the following lemma which facilitates the assignment of a ‘C’ or ‘K’ state to the vertex v with respect to the edge e. Lemma 4. Let an edge e be defined by two vertices v and v 0 . If a touch point exists in HC(e,v), the e starts to shrink at the vertex v. If a touch point exists in HK(e,v), the edge e disappears at the vertex v. If an edge e has both ‘C’ and ‘K’ vertices, e starts to shrink at the ‘C’ vertex and disappears at the ‘K’ vertex during the region-expansion. Therefore, e should have an end-out state. However, if the event time at the ‘K’ vertex is earlier than that of the ‘C’ vertex, then the edge has a split state. If both vertices of e are in ‘C’ states, e starts to shrink from both end vertices to disappear at a point in the middle of e. Hence, the edge e has a mid-out state. Suppose that both vertices of an edge e have ‘K’ states. Then, e should disappear at both end vertices. Note that this observation implies that the edge e should first go through a state of shrinking before it disappears. Since neither vertex contributes to the shrinking, there has to be a point in the middle of e from which the shrinkage starts. This constraint can be satisfied only by splitting e at a point in the middle of e. Hence, the edge e has a split state. In the general cases of the Voronoi diagram for spheres, there can be a case in which an edge does not have any end vertex. Suppose that a small sphere is located in the convex hull of two larger spheres. Then, the edge defined by these three spheres can be an ellipse without any vertex in its own

mathematical definition. In [12,13,15], we call this the big-brothers problem. The details on the big-brothers problem will be discussed in Section 8, and we focus here on the state of the elliptic edge. An elliptic edge e without any vertex can have nothing but a split state since the only first event that can be associated with e is the tangential touch between e and a moving on-edge of the expanding region. Hence, what should happen to e in the next moment is to split e at this tangent point. However, caution should be exercised here since e does not split into two pieces, but two end vertices are created on e so that e is now a part of the ellipse. The next theorem summarizes the above discussions how to detect edge states from vertex states, and provides the basis for the detection of events. Theorem 2. The state of an edge e is determined by the states of its two end vertices as follows: † If e has both ‘C’ and ‘K’ vertex states, e has neither an end-out or split state. † If e has two ‘C’ vertex state, e has a mid-out state. † If e has two ‘K’ vertex state, e has a split state. † If e has no vertex, e has a split state. 5.2. Event detection via edge states The type of event associated with an edge can be determined by the state of the edge. As was stated before, if an edge has an end-out state, the event corresponding to the edge is either oneend-event or two-end-event. Hence, we need to discriminate these two cases for an edge with an end-out state. When an edge has a mid-out or a split state, the edge directly corresponds to a mid-event or a split-event, respectively. 5.2.1. One-end-event and two-end-event If an edge e has an end-out state, the event corresponding to e occurs at an end vertex v of e. Note that v is an off-vertex. At the moment when an expanding generator touches the vertex sphere SðvÞ of v, the topology structure around v is not valid any more and has to be changed. Recall that the state of a vertex is defined with respect to an edge, and v is of degree four. Hence, v is associated with four vertex states: one for each edge incident to v. Suppose that some edges incident to v contribute ‘K’ vertex states to v. As explained earlier, these edges will disappear after the moment when an expanding generator touches SðvÞ. The other incident edges, i.e. the edges contributing to ‘C’ vertex state to v, start to shrink at this moment. It is important to note that there can be either one or two edges contributing the ‘K’ vertex state to v. Depending on the number of such edges, different changes may occur in the topology at v. In Fig. 2(a), only the edge e contributes the ‘K’ vertex state to the vertex v and therefore e disappears when the moving onvertex (shown as the black dot) meets v. The other edges e1, e2, and e3 incident on v start to shrink at that moment. We call such a change in the topology a one-end-event.

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

In Fig. 2(b), on the other hand, two edges e 0 and e 00 are incident on v, and they contribute the ‘K’ vertex state to v. In this case, these two edges will disappear simultaneously when two on-vertices (shown as the black dots) simultaneously meet v. As shown in the figures, this topological change is different from the one-end-event case and we call it a two-end-event. An important question is whether three or more incident edges can contribute ‘K’ vertex states to v. The answer is that such a situation is impossible. If three edges do contribute ‘K’ values to v, then the Voronoi region corresponding to these three edges disappears when the edges disappear. However, this cannot be realized in the region-expansion process since a generator sphere has always its Voronoi region unless the generator is fully contained in another. 5.2.2. Mid-event Suppose that there is an edge e having a mid-out state. When e disappears in the course of the region-expansion process, a mid-event occurs. Since there is no any other edge that can disappear at the same time and at the same place with e the edge with a mid-out state triggers a mid-event. Hence, a midevent can be detected by an edge having a mid-out state. Fig. 2(c) illustrates an example of a mid-event where the edge e has a mid-out state and both two on-vertices (in the black dots) have ‘C’ vertex states. In the case of a mid-event, one face always disappears with the edge e with a mid-out state, e.g. f1 in Fig. 2(c). To handle a mid–event correctly, such a face has to be appropriately identified. The disappearing face has the following properties: (i) the face is incident to e, and (ii) the face is bounded by two edges. Therefore, one of the three faces incident on e is the desired face. If there is only one face satisfying the above two properties, this face is the desired one. Unfortunately, however, it is possible that more than one face satisfy both properties, and in such a case we have to find the correct one among possibly three faces. Imagine the moment that an edge e, which has a mid-out state, disappears. At this moment, e degenerates to a point q and the moving vertex sphere SðqÞ becomes fixed at q. Note

423

that SðqÞ is now simultaneously tangent to four generators (three generators defining e, and an expanding generator). In this situation, four tangent points on SðqÞ are coplanar and therefore cocircular. Note that this circle is not necessarily a great circle on SðqÞ. Let P be a plane where these four tangent points and the circle are placed as shown in Fig. 4. Suppose that p1, p2, and p3 denote three tangent points associated with the generators s1, s2, and s3 defining e, respectively, and p4 denotes the tangent point associated with the expanding generator. Without loss of generality, we assume that p1 through p4 are placed in the counter-clockwise orientation on the circle. Let aij denote the arc starting at pi and ending at pj in the counter-clockwise orientation as well. Then, aij corresponds to the Voronoi face defined by two generators si and sj. Therefore, if p4 is located on the arc aij, the Voronoi face defined by si and sj cannot exist and therefore disappears. 5.2.3. Split-event The case of a split-event is quite similar to the case of the mid-event in the sense that the event occurs at a point in the middle of an edge, not at a vertex. Therefore, it is not necessary to find other edges incident to the edge of a split state. However, even in the split-event, we should find a radiatingface and the on-edge on the radiating-face, which will cause the split-event. 6. Event time Once an event type is determined, it is necessary to know when the event takes place during the region-expansion process so that we can handle the event at the right time. In the region-expansion algorithm, events are appropriately ordered and handled according to their event time. There can be two approaches for such an ordering. The first approach is to compute the absolute event time that an event occurs. Note that the absolute event time is defined as the radius of an expanding generator when an event occurs. Once event times are computed, they are sorted in ascending order. Note that the region-expansion algorithm handles events in a manner similar to the discrete event simulation. The second approach is to find local precedence relationships for each event. These local precedences eventually lead to the finding of a global ordering among all events. 6.1. Absolute event time An event time is defined as the radius of the expanding generator when the event occurs. Hence, the current time in the process is given by the current radius of the expanding generator.

Fig. 4. Out-face detection in the mid-event.

6.1.1. One-end-event and two-end-event In the case of one-end-event and two-end-event, the onvertex of the event edge is moving along the edge toward the off-vertex on the edge. At the time when the on-vertex meets the off-vertex, the edge degenerates to a point resulting in an

424

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

end-event. Therefore, the event time of the event edge can be computed as the distance between the center of the expanding generator and the surface of the vertex sphere of the off-vertex. 6.1.2. Mid-event and split-event While the event time of an end-event can be easily computed, the event time for the mid- or split-event cannot be computed easily. Without loss of generality, we translate four generators (three generators defining the event edge, and an expanding generator) so that the expanding generator is centered at the origin. Let (xi,yi,zi) and riR0, iZ1,2,and,3, be the transformed centers and the radii of three generators defining the edge, respectively. Then, the equi-distant point (x,y,z) on the edge at a fixed time t in the expansion process can be formulated as the following system of equations: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxKx1 Þ2 C ðyKy1 Þ2 C ðzKz1 Þ2 Z r C r1 (1) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxKx2 Þ2 C ðyKy2 Þ2 C ðzKz2 Þ2 Z r C r2

(2)

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxKx3 Þ2 C ðyKy3 Þ2 C ðzKz3 Þ2 Z r C r3

(3)

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x 2 C y 2 C z2 Z r C t

(4)

where r is the distance between the equi-distant point and the generators. Gavrilova used a similar formulation to detect the topology change in the Voronoi diagram of spheres in a dynamic setting [8]. Note that r is negative if the equi-distant point is located inside three generators defining the target edge. Rearrangement of the above system yields the followings. x1 x C y1 y C z1 z C r1 r Z ðx21 C y21 C z21 Kr12 C 2rt C t2 Þ=2

(5)

x2 x C y2 y C z2 z C r2 r Z ðx22 C y22 C z22 Kr22 C 2rt C t2 Þ=2

(6)

x3 x C y3 y C z3 z C r3 r Z ðx23 C y23 C z23 Kr32 C 2rt C t2 Þ=2

(7)

x2 C y2 C z2 Z ðr C tÞ2

(8)

r C r1 O 0

(9)

r C r2 O 0

(10)

r C r3 O 0

(11)

Fig. 5. Graph representation of the four types of event.

r C tO 0

(12)

From Eqs. (5)–(7), (x,y,z) can be found in the function of r and t and substituting them in Eq. (8) produces a standard form of a quadratic equation in r, f(r):Zar2CbrCcZ0. In this equation, the coefficient a is a quadratic in t, b is cubic in t, and c is quartic in t. For t to be a valid event time for the mid- and split-event, f(r) should have a double root for r, where t is fixed as a constant. This double root corresponds to the case when the two vertices of the edge e in Fig. 2(c) meet at a point in the middle of e. The case where f(r) has two distinct roots corresponding to that of the two vertices is distinct. Therefore, the discriminant of f(r), denoted by dr(t):Zb2K 4ac, has to be zero. If dr(t)O0, two different r values can be obtained. If dr(t)!0, on the other hand, there is no real r value and therefore there is no point equi-distant from four generators. It turns out that the discriminant dr(t):Zb2K4ac is a quartic polynomial in t. Therefore, if we can find the roots of dr(t), some of them in the range of (0,r0] become candidates for the desired event time, where r0 is the predefined maximum radius of the expanding generator. From the equation f(r)Z0, we can obtain the corresponding r value as Kb2a if t is given since a and b are polynomials in t. Therefore, among the candidate t values, the one which satisfies the constraints (9)–(12) may be the event time of mid- or split-event. While there are two valid t values, the one with a smaller value is the desired one. 6.2. Relative event time For each event in the region-expansion, we can define precedence relationships between this event and some of the other events. There is a unique precedence relationship between the event of an event edge and the event of another edge incident to the event edge. 6.2.1. One-end-event Refer to Fig. 2(a). After a one-end-event occurs at the vertex v, the events that are associated with edges e1, e2, and e3 will occur. Note these edges are incident on v. Therefore, there are definite precedence relationships, three in total in the one-endevent, between the event of e and the events of the edges incident on v. However, the precedences among the events of e1, e2, and e3 are not known. Note that each subfigure in Fig. 5 matches each case in Fig. 2. Note again that the arrow in Fig. 5 denotes the precedence relationship between two consecutive events. The box in the figure denotes an event where the upper- and lowercells denote the event type and the location where an event occurs, respectively. Fig. 5(a), therefore, indicates that a one-end-event corresponding to an edge e occurs at vertex v, and edges e1, e2, and e3 trigger the events associated with themselves after this event occurs. It can be shown that all types of events, except a splitevent, may follow the one-end-event since all of e1, e2, and e3 have ‘C’ vertex states at v.

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

6.2.2. Two-end-event As shown in Fig. 2, a two-end-event occurs at vertex v which is triggered by two edges e 0 and e 00 . Then, the events associated with edges e1 and e2, which are incident on v, will occur. Similarly to the one-end-event, we can make precedence relationships associated with a two-end-event as shown in Fig. 5(b). As in a one-end-event, the possible type of the following events is either a one-end-event, a two-end-event, or a mid-event. 6.2.3. Mid-event Fig. 5(c) shows the precedence relationship due to a midevent. Recall that a mid-event occurs at a point in the middle of an edge while the one-end-event or the two-end-event occurs at a vertex. Hence, the event, which precedes a mid-event, must have occurred at the two end vertices of the edge. Therefore, there should be two events, corresponding to the two vertices and following the edge, preceding a mid-event and this is shown in Fig. 2(c). In the case of a mid-event, no other event directly follows the mid-event since all edges associated with this event become on-edges after processing the mid-event, and an on-edge does not have an event. 6.2.4. Split-event In the case of a split-event, there is no preceding event, only following events. The event edge splits into two new edges, and they are again associated with their own events. Note that these new edges are always radiating-edges and each new edge can have either an end-out state or a mid-out state. Therefore, the event immediately following a split-event can be one-endevent, two-end-event, or mid-event. 6.2.5. Representation of event precedences in DAG From the events and their precedence relationships, we can define a directed acyclic graph (DAG) where nodes are events and arcs are the precedence relationships between events. A DAG is uniquely defined for each expanding generator when it is to be expanded, and the DAG does not change during the whole process of region-expansion for a generator. A DAG can be constructed by verifying the events associated with edges and appropriately stitching nearby events with directed arcs while following the appropriate edges starting from those on the initial radiating-faces of the generator. Since generators have predefined radii, the search range to verify events is determined by the size of an expanding generator. In general, this search space is rather limited in the near neighborhood of an expanding generator. Hence, the number of nodes and arcs of DAG can be determined for an expanding generator prior to the actual processing of the region-expansion. Once a DAG is constructed, the topological ordering among all nodes can be made in a linear time of the number of nodes in the worst-case [10]. Note that the ordering in DAG may not be unique, and thus the sequence of events obtained from the relative event time approach may be different from the one obtained from the absolute event time approach. However, it is guaranteed that the topology of the resulting Voronoi diagram

425

after expanding the region for a generator to its full size can be consistently maintained even though numerical errors exist. Provided a sufficient level of numerical accuracy, the results of both absolute and relative approaches are correct. Given a DAG representation of the precedence among events, there can be multiple valid orderings of events. If an ordering satisfies the precedence requirement, any ordering will produce an identical result for the expansion of a generator. For example, refer to Fig. 5(b). In this case, both events before the arcs e 0 and e 00 should be processed before the current event. The order of processing e1 and e2 is not important unless it conflicts with some other node later. Since the topological ordering among all nodes can be found in the linear time for the number of nodes, we can appropriately order all of the events for the current region-expansion process according to the precedence relationships represented as the DAG in the linear time of the number of nodes in the worstcase. Note that the number of nodes in DAG can be, in the worst-case, O(n2), which is also the number of edges in the Voronoi diagram. Hence, the next lemma follows. Lemma 5. The region-expansion of a generator can be done in O(n2) time in the worst-case, where n is the number of generators.

7. Handling of infinity Infinity has to be appropriately handled in the computation of a Voronoi diagram. We treat infinity by introducing one, and only one, fictitious generator, which represents the infinity. Then, each of the spheres contributing to the convex hull boundary shares a fictitious face with the fictitious generator. We call this face an infinite face. Hence, there are the same number of infinite faces as the number of spheres on the convex hull boundary. In addition, the infinite face is bounded by a number of infinite edges, and every infinite edge consists of two infinite vertices. Note that the infinite vertices, edges, and faces altogether represent the convex hull of spheres. The infinite vertices, edges, and faces correspond to the planar, conic, and spherical parts of the convex hull boundary, respectively. Topology changes related to infinity can also be handled in a similar way as ordinary edges. For an infinite vertex vN, one ordinary edge and three infinite edges are incident to vN. Therefore, only three generator spheres are associated with vN. 7.1. Ordinary edge with an infinite vertex To determine the state of an edge incident to vN, regardless whether it is ordinary or infinity, the vertex sphere SðvNÞ at vN can be regarded as a common tangent plane to the three generators defining the edge. Note that there are eight possible tangent planes on three generators. Among them, we select the one that coincides with the convex hull boundary. To determine the state of vN with respect to an ordinary edge incident to vN, we need to know where the touch point (the tangent point between an expanding generator and SðvNÞ)

426

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

is located in the plane SðvNÞ. Depending on the position of the touch point in the plane, a different state is assigned to vN. If the touch point is in the circle passing through three fixed tangent points (the tangent points between three generators and SðvNÞ), the ordinary edge has to disappear and the expanding generator will contribute to the convex hull after the contact. Therefore, the region inside the circle can also be denoted by HK(e,vN) as the usual cases, where e is an ordinary edge. As a result, the infinite vertex vN has a ‘K’ state. Suppose that the touch point is located outside the circle. Then, the tangent sphere SðvNÞ), which coincides with a planar face on the convex hull boundary, should obviously disappear. In this case, a new tangent sphere can be defined among the expanding generator and the three generators corresponding to the edge. Therefore, outside of the triangle can be regarded as HC(e,vN) and vN obtains a ‘C’ state. 7.2. Infinite edge As explained earlier, an infinite edge eN corresponds to a conic part of the convex hull boundary and is related with two generators s1 and s2. Let H(eN,vN) be a line passing through two tangent points between SðvNÞ and both s1 and s2. Let HC(e,NvN) denote the halfplane limited by H(eN,vN) on SðvNÞ and contains the

reference tangent point. Then, HK(eN,vN) denotes the other side of HC(eN,vN) on the plane. Suppose that a touch point lies in HK(eN,vN). In this case, the conic face corresponding to s1 and s2 on the convex hull boundary cannot exist since no common tangent plane can be defined between s1 and s2 after the contact. Hence, eN cannot exist and disappears, and therefore vN is assigned with a ‘K’ state as is usually done. Similarly, when a touch point is in HC(eN,vN), the corresponding planar face on the convex hull boundary starts to change while remaining tangent with s1 and s2. Therefore, the edge eN has a ‘C’ state at vN. Depending on the states of both end vertices, the event type of an infinite edge is determined by Theorem 2 in exactly the same way it is usually done. It is interesting to observe that an identical set of state types apply to both infinite edges as well as ordinary edges. 7.3. Event time An event time of an infinite edge can be computed slight differently depending on the event type. If an infinite edge eN has a one-end-event or two-end-event, the event time can be calculated as the distance between the expanding generator and the vertex sphere SðvNÞ if eN contributes a ‘K’ state to vN. Recall that SðvNÞ is a plane. In the case that eN has a mid-event, the event time can be obtained as the distance between the center of the expanding generator and the farthest supporting tangent plane of two generators s1 and s2 corresponding to the edge. In the case of a split-event, on the other hand, the distance between the center of expanding generator and the closest supporting tangent plane between s1 and s2 becomes the event time of the edge. 8. Disconnectedness in the edge graph

Fig. 6. Example of the disconnected case.

Suppose that a Voronoi diagram is represented as a graph consisting of Voronoi vertices and Voronoi edges. It is known that the Voronoi diagrams for points can be represented by a graph of a single connected component. It is possible, however, that the graph for the Voronoi diagram of spheres may consist of multiple connected components. Suppose that some small spheres are contained in the convex hull of two larger spheres. Then, the graph formed among these spheres may be disconnected from the others since there may be no edge going to the infinity. This situation is called a big-brothers problem in the Voronoi diagram of spheres [13]. Fig. 6(a) shows an example of the disconnected case. In this example, two mid-sized spheres are located in-between two larger spheres. In between those two mid-sized spheres, there is another small sphere. Four other relatively large spheres are not related to the big-brothers problem in this example. Note that the edges associated with the mid-sized spheres are disconnected from the edges defined among larger spheres. Thus, the Voronoi face between the two larger spheres contains a topological hole (or sometimes called as an island) inside the face. Fig. 6(b) is a magnified image of Fig. 6(a). As shown in

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

the figure, it is possible that the edge graph defined between the mid-sized spheres is again disconnected from the edge graph due to the smaller sphere. Proper handling of the big-brothers problem in a single frame work of an algorithm for the Voronoi diagram is usually difficult as stated in [13,22]. In the proposed region-expansion algorithm, however, this is not a special case at all, and therefore can be handled as just a usual case of all disconnected edges without any extra computation or code. Fig. 7 illustrates the situation before a related event occurs (Fig. 7(a)) and after the event has occurred (Fig. 7(b)). Note that the edge e is the event edge to be investigated in this example and the face f1 is not an on-face but a radiating-face. When a transition is made from Fig. 7(a) to (b), the event associated with the event edge e is always a mid-event. When a mid-event is associated with an edge e, there are two on-faces f2 and f3 for the expanding generator and these faces are connected to the two vertices of e. In general, f2 and f3 are topologically distinct and are merged into a single face as shown in Fig. 2(c). If, however, f2 and f3 are topologically identical faces, the merging of f2 and f3 via the mid-event results in the disconnected edge graph as illustrated in Fig. 7(b). Note that it is also possible to reverse the above process in the middle of a region-expansion, i.e. from Fig. 7(b) to (a). As before, this is not an extra-ordinary case as well, and can be handled by as usual since a split-event corresponds to the event edge, either e1 or e2 in Fig. 7(b). Even though there is no such a disconnected case in the final Voronoi diagram for spheres, some disconnected cases may happen during the region-expansion process. Note that the disconnected edges are connected via split-events. 9. Algorithm The region-expansion algorithm for a single Voronoi region can be described by the following steps. The whole Voronoi diagram for all spheres can be obtained by repeating the following algorithm for all generators, starting with the ordinary Voronoi diagram for the centers of all spheres. In this section, we provide two algorithms: one is based on the absolute event time, and the other one is based on the relative event time (using DAG). Note that an appropriate topology data structure can support the iteration in a coherent scheme. (a)

427

Algorithm 1. (Expand region of siZ(ci,ri)) 1. Find all the edges bounding radiating-faces and insert them into a set E. 2. For each edge er2E\Eon, determine its event time t and state. If t!ri, insert er into the event queue Q which implements a priority queue. 3. Pop an edge eq from Q and determine the corresponding event. If eq has end-out state, check if eq causes one-endevent or two-end-event by watching the next edges in Q with the same event time. 4. Perform an appropriate action for the detected event. After treatment, find new edges bounding new radiating-faces and insert them into a set Enew. Perform Step 2 for all edges in the edge set Enew. 5. Repeat Steps 3 and 4 until Q is empty. Algorithm 2. (Expand region of si using DAG) 1. Initialization Find all the initial radiating-edges and insert them into a queue Q. Find all the split stated edges bounding radiating-faces, and insert them into Q. 2. DAG construction While Q is not empty, do the followings: Pop an edge eq from Q. If the event time t of eq is greater than ri or eq is already processed, do nothing. Else if the state of eq is end-out, † determine if the corresponding event is one-endevent or two-end-event, † insert these incident edges to the event vertex except the event edge(s) into Q, † find split stated edges bounding the faces which will be radiating-faces after the event, and insert them into Q, and † construct a DAG node and arcs for the detected event. Else if the state of eq is mid-out, construct a DAG node and connect arcs for a mid-event. Else if the state of eq is split, † find split stated edges bounding the faces which will be radiating-faces after the event, and insert them into Q, and (b)

Fig. 7. Disconnected case.

428

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

Fig. 8. Examples of the Voronoi diagram for spheres: (a) Voronoi diagram for 15 spheres. The Voronoi region of the sphere in the center is shaded, (b) Voronoi diagram for a protein data (PDB [26] code: 3rec) consisting of 63 atoms (25 Hs, 20 Cs, 10 Os, 7 Ns, and 1 P).

† construct a DAG node and two arcs for a split-event, and insert two split edges into Q. 3. Make a topological ordering from DAG and perform appropriate actions for each event by the sequence in the ordering. Algorithm 1 takes O(n2 log n) time in the worst-case to expand a Voronoi region for an expanding generator, starting from a center point of the generator to a complete sphere. Note that there can be O(n2) number of edges in the Voronoi diagram of spheres as well as points in three dimension, and sorting all events is necessary in the algorithm. (Note that the combinatorial complexity of a single Voronoi region is known as Q(n2) [29].) Therefore, the whole Voronoi diagram can be constructed by Algorithm 1 in O(n3 log n) time in the worst-case. Algorithm 2, however, takes O(n2) time in the worst-case since all events in the region-expansion for a generator are stored in DAG and an appropriate ordering of all events is identified in the linear time of the nodes in the DAG. Note that the number of nodes in the DAG is O(n2) time in the worstcase. Therefore, the whole Voronoi diagram can be constructed by Algorithm 2 in O(n3) time in the worst-case. We want to mention here that reducing the sizes of all generators by the smallest one at the outset of the algorithm results in much faster computation, even though the reduction is a constant factor.

Hence, the computations for the other geometries for Voronoi edges and faces are not included. Note that most applications require only the topologies of Voronoi diagrams. However, we can compute the geometries for Voronoi edges and faces as suggested in [13] if necessary. Recall that the worst-case time complexity of Algorithm 1 is O(n3 log n) while it is O(n3) for the edge-tracing algorithm. However, the empirical experiment shows a rather surprising result as shown in the figure. Fig. 9(b) shows the computation time required by Algorithm 1 only shows a roughly linear increase. This result is mainly due to the fact that proteins consist of atoms whose size difference is not too big and the atoms are distributed somewhat regularly. Therefore, the

10. Experiments and discussions We have fully implemented the proposed algorithm and tested it with various protein data models available in the Protein Data Bank (PDB) [26]. Fig. 8 shows the computed Voronoi diagrams for two models: one for a hypothetical system consisting of 15 spheres with three different radii, and the other for a protein model coded 3rec in PDB which consists of 63 atoms. For the visual convenience, only the Voronoi edges are displayed. We have compared the performance of the proposed regionexpansion algorithm (Algorithm 1) against the previously reported edge-tracing algorithm [13] and the result is shown in Fig. 9(a). The statistics in this figure reflects only the time required by the computation of topologies of the Voronoi diagrams and includes the computation of Voronoi vertices.

Fig. 9. Computation time: (a) comparison between the Edge-Tracing Algorithm and Region-Expansion Algorithm and (b) running time of the RegionExpansion Algorithm.

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

number of events actually happened during the region expansion of a single generator could be in a small constant number. This in fact verifies that the observations by Halperin and Overmars [11] for protein models are correct. It is worth noting that the order of expanding spheres in the region-expansion process may influence the amount of computation. This needs to be further investigated in the future for the performance optimization of the algorithm. To make the region-expansion algorithm robust against numerical errors, the analyses should be done for some basic numerical computations such as the computations of Voronoi vertices and event times. In geometric algorithms, such numerical errors may cause inconsistency in the computed topological structure to result in an undesired output as an infinite loop or a crash in the middle of computation. In the region-expansion algorithm, there are several numerical computations. First, vertex spheres and vertex states need to be exactly computed since the vertex spheres are used in the decision of vertex states and in the computation of the event times for one-end-event and two-end-event. Second, the computation of the event times for mid-event and split-event requires solving a quartic polynomial equation, which requires high computational cost. Moreover, the coefficients of the quartic polynomial consist of many terms. In our implementation, for example, each coefficient consists of roughly 1000 terms where each term requires a few multiplications. Hence, the exact computation for this equation would require a significant amount of computation as well as memory.

11. Conclusion In this paper, we have presented a region-expansion algorithm for constructing a Voronoi diagram for spheres by handling events from the ordinary Voronoi diagram for the centers of the spheres. The algorithm can compute the whole Voronoi diagram in O(n3) time in the worst-case where n is the number of spheres. In addition, most of the algorithm has been fully implemented and the runtime performance shows a linear pattern for protein data. The edge-disconnectedness problem in the Voronoi diagram can be incorporated in a single framework of the region-expansion algorithm. Acknowledgements This research was supported by the Creative Research Initiatives from the Ministry of Science and Technology in Korea. 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.

429

[3] Boissonnat J-D, Karavelas MI. On the combinatorial complexity of Euclidean Voronoi cells and convex hulls of d-dimensional spheres. In: Proceedings of the 14th Annual ACM-SIAM Symposium on Discrete Algorithms; 2003. p. 305–12. [4] Cho Y, Kim D, Kim D-S. Topology representation for the Voronoi diagram of 3D spheres. Int J CAD/CAM 2005; 5(1). [5] Choi GH, Han SH, Lee HC. Optional storage of non-manifold information. Trans Soc CAD/CAM Eng 1997;2(3):150–60. [6] Choi Y. Vertex-based boundary representation of non-manifold geometric models. PhD thesis, Carnegie-Mellon University, USA; 1989. [7] Gavrilova M. Proximity and applications in general metrics. PhD thesis, Department of Computer Science, The University of Calgary, Calgary, Canada; 1998. [8] Gavrilova M, Rokne J. Updating the topology of the dynamic Voronoi diagram for spheres in Euclidean d-dimensional space. Comput Aided Geom Des 2003;20(4):231–42. [9] 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. [10] Goodrich MT, Tamassia R. Data Structures and Algorithms in Java. 2nd ed. New York: Wiley; 2001. [11] Halperin D, Overmars MH. Spheres, molecules, and hidden surface removal. In: Proceedings of The 10th ACM Symposium on Computational Geometry; 1994, pp. 113–122. [12] Kim D-S, Cho Y, Kim D. Edge-tracing algorithm for Euclidean Voronoi diagram of 3D spheres. In: Proceedings of the 16th Canadian Conference on Computational Geometry; 2004. p. 176–9. [13] Kim D-S, Cho Y, Kim D. Euclidean Voronoi diagram of 3D balls and its computation via tracing edges. Comput Aided Des 2005;37(13):1412–24. [14] Kim D-S, Cho Y, Kim D, Cho C-H. Protein sructure analysis using Euclidean Voronoi diagram of atoms. In: Proceedings of the International Workshop on Biometric Technologies (BT2004); 2004. p. 125–9. [15] Kim D-S, Cho Y, Kim D, Kim S, Bhak J, Lee S-H. Euclidean Voronoi diagrams of 3D spheres and applications to protein structure analysis. In: Proceedings of the 1st International Symposium on Voronoi Diagrams in Science and Engineering (VD2004); 2004. p. 137–44. [16] Kim D-S, Kim D, Sugihara K. Voronoi diagram of a circle set from Voronoi diagram of a point set: I. Topology. Comput Aided Geom Des 2001;18:541–62. [17] Kim D-S, Kim D, Sugihara K. Voronoi diagram of a circle set from Voronoi diagram of a point set: II. Geometry. Comput Aided Geom Des 2001;18:563–85. [18] Lee K. Principles of CAD/CAM/CAE systems. Reading, MA: Addison Wesley; 1999. [19] Lee SH, Lee K. Compact boundary representation and generalized Euler operators for non-manifold geometric modeling. Trans Soc CAD/CAM Eng 1996;1(1):1–19. [20] Lee SH, Lee K. Partial entity structure: a compact boundary representation for non-manifold geometric modeling. ASME J Comput Inf Sci Eng 2001;1(4):356–65. [21] Luchnikov VA, Medvedev N, Naberukhin YI, Schober HR. VoronoiDelaunay analysis of normal modes in a simple model glass. Phys Rev B 2000;62:3181–9. [22] Luchnikov VA, Medvedev NN, Oger L, Troadec J-P. Voronoi-Delaunay analysis of voids in systems of nonpherical particles. Phys Rev E 1999; 59(6):7205–12. [23] Montoro JCG, Abascal JLF. The Voronoi polyhedra as tools for structure determination in simple disordered systems. J Phys Chem 1993;97(16): 4211–5. [24] Naberukhin YI, Voloshin VP, Medvedev NN. Geometrical analysis of the structure of simple liquids: percolation approach. Mol Phys 1991;73: 917–36. [25] Okabe A, Boots B, Sugihara K, Chiu SN. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. 2nd ed. Chichester: Wiley; 1999. [26] RCSB Protein Data Bank Homepage. http://www.rcsb.org/pdb/ [27] Richards FM. The interpretation of protein structures: total volume, group volume distributions and packing density. J Mol Biol 1974;82: 1–14.

430

D. Kim, D.-S. Kim / Computer-Aided Design 38 (2006) 417–430

[28] Sastry S, Corti DS, Debenedetti PG, Stillinger FH. Statistical geometry of particle packings. I. Algorithm for exact determination of connectivity, volume, and surface areas of void space in monodisperse and polydisperse sphere packings. Phys Rev E 1997;56:5524–32. [29] Voloshin VP, Beaufils S, Medvedev NN. Void space analysis of the structure of liquids. J Mol Liq 2002;96–97:101–12. [30] Will H-M. Computation of additively weighted Voronoi cells for applications in molecular biology. PhD thesis, Swiss Federal Institute of Technology, Zurich; 1999. Donguk Kim is a senior researcher in Voronoi Diagram Research Center at Hanyang University, Seoul, Korea. He received his B.S., M.S. and Ph.D. degrees from Hanyang University in 1999, 2001 and 2004, respectively. His research interests include computational geometry, geometric modeling and their applications in the molecular biology.

Deok-Soo Kim is a professor in Department of Industrial Engineering, Hanyang University, Korea. Before he joined the university in 1995, he worked at Applicon, USA, and Samsung Advanced Institute of Technology, Korea. He received a B.S. from Hanyang University, Korea, an M.S. from the New Jersey Institute of Technology, USA, and a Ph.D. from the University of Michigan, USA, in 1982, 1985 and 1990, respectively. His current research interests mainly lie in the theory and applications of Voronoi diagram while he has been interested in various geometric problems. He is current the director of Voronoi Diagram Research Center supported by the Ministry of Science and Technology, Korea.