6.3. Generalization of the Voronoi Diagram

6.3.1. Voronoi Diagram and Delaunay Triangulation in 3D

The Voronoi diagram in 3D for a set of points S subdivides the 3D space into cells of the same neighborship and can be represented by a graph in 3D. For convenience, we assume that the points are in general position; for example, no five points lie on a common sphere. The corresponding data structure is given by a DCEL in 3D. We will see that the incremental construction approach is a simple and effective construction scheme in the 3D case.

Formal description. Formally, we extend the notations of the 2D case. The set S denotes the set { p ₁ , p ₂ , , p _n } of n points in 3D. The 3D bisector Bis( p _i , p _j ) of two points p _i , p _j ˆˆ S is defined to be a plane through the midpoint of the line segment p _i p _j that is perpendicular to the segment p _i p _j . A Voronoi region VoR( p _i , S ) is given by the intersection of the half-spaces bounded by bisectors Bis( p _i , p _j ) for all j ‰ i . As an intersection of half-spaces, the boundary of VoR( p _i , S ) consists of faces, edges, and vertices, and represents a 3D convex polyhedron.

Complexity. The complexity is as follows .

Obviously, the number of faces of VoR( p _i , S ) is at most n ˆ 1. Each half-space defined by Bis( p _i , p _j ) will contribute to at most one face.
Applying the Eulerian formula, the number of edges and vertices of VoR( p _i , S ) is in O ( n ).
Counting the complexity of all regions , the total number of components of the diagram VD( S ) in 3D is in O ( n ² ).
Unfortunately, there are configurations of S with ( n ² ) components; see [Dewdney and Vranch 77].
Fortunately, in practical situations, the complexity of VD( S ) in higher dimensions is linear. For example, if the points of S are drawn uniformly at random in a unit ball, the expected size of the diagram is in O ( n ); see [Dwyer 91].

Altogether, we can assume that, in practical situations, VD( S ) is a geometric graph in 3D of linear size.

Delaunay triangulation DT( S ) in 3D. In analogy to the 2D case, the Delaunay triangulation DT( S ) in 3D is defined to be the geometric dual of VD( S ) with vertices from S .

For two points p _i and p _j that share a common face in the VD( S ), there is an edge p _i p _j in DT( S ).
Every edge of the VD( S ) subdivides three faces. Thus, for every edge in the VD( S ), there is a triangle in DT( S ).
At a vertex v of VD( S ), four points of S have the same distance to v and the four points share four edges with each other. Therefore, there is a tetrahedron for the vertex v .

Note that we have assumed general position. Altogether, DT( S ) is a graph in 3D consisting of tetrahedra.

Equivalently, DT( S ) may be defined extending the empty circle property. If the circumsphere of four points of S does not contain another point of S inside or on the boundary, the corresponding tetrahedron of the four points is a Delaunay tetrahedron. The circumcenters of these empty spheres are just the vertices of VD( S ). In analogy to the triangulation in 2D, the Delaunay triangulation DT( S ) is a partition of the convex hull of S into tetrahedra, provided that S is in general position.

There are robust implementations of the empty sphere property (see Chapter 9), and the general position condition can be easily fulfilled for the test in 3D (see Section 9.5). Therefore, it is highly recommended to extend the incremental construction scheme of the previous section to 3D.

Simple incremental construction in 3D. We consider an incremental construction of the dual DT( S ). Some algorithms for the incremental construction of VD( S ) in 3D also exist; see [Watson 81], [Field 86], [Tanemura et al. 83], or [Inagaki et al. 92].

In the dual environment, analogously to the 2D case, we have to perform edge flips for removing all tetrahedra of DT _i _ˆ ₁ that are in conflict with p _i . That is, the circumsphere of the tetrahedra ( p _j , p _k , p _l , p _m ) of four points p _j , p _k , p _l , and p _m contains p _i . In analogy to the 2D case, this can happen only if p _j , p _k , p _l , p _m , and p _i are in convex position.

Similar to the 2D case, we make use of a 3D tetrahedral DAG and find the tetrahedron, say ( p _j , p _k , p _l , p _m ), of DT _i _ˆ1 that contains p _i . With similar arguments, we first insert the edges p _i p _j , p _i p _k , p _i p _l , and p _i p _m , thus obtaining three new tetrahedra; see Figure 6.9.

Figure 6.9: The tetrahedron (p _j , p _k , p _l , p _m ) contains p _i and four edges p _i p _j , p _i p _k , p _i p _l , and p _i p _m are inserted. For every face of P( ˜… (p _i )), we check the outer tertrahedra; in this case, (p _j , p _l , p _m , q) with the opposite p _i .

Now we obtain P ( ˜… ( p _i )), the set of outer triangles of ( p _j , p _k , p _l , pm ) opposite p _i . In the 2D case, P ( ˜… ( p _i )) was a set of edges opposite p _i . Similar to the 2D case, we go on with p _i and with P ( ˜… ( p _i )). In the 2D case, every edge of P ( ˜… ( p _i )) may subdivide two triangles. Analogously, every triangle of P ( ˜… ( p _i )) may subdivide two tetrahedra as follows. Every triangle of P ( ˜… ( p _i )) can be a ground face G of a tetrahedron with a point q outside P ( ˜… ( p _i )) and opposite p _i , and G builds a tetrahedron with p _i . Together with p _i , we have to check whether the four points of the outer tetrahedron and p _i are in conflict. That is, does the circumcenter of G and q contain p _i ? This cannot happen if the five points are not in convex position, since p _i lies outside the tetrahedron of G and q in this case. This fact can be proven similar to the 2D case; see Figure 6.5. If the five points are in convex position, the triangle G is called flippable .

If G is flippable, we apply the circumsphere test. If G is no longer Delaunay, i.e., the circumsphere test returns true, then we can either remove an edge or we can insert an edge. Note that in the 2D case, we replace one edge by another, so here we have one of the main differences between 2D and 3D. The nature of the 3D edge flip depends on the situation; see Figure 6.10. In the first case, we replace two tetrahedra by three tetrahedra in the tetrahedral decomposition. In the second case, we replace three tetrahedra by two tetrahedra. We refer to the first case as a 2 3 edge flip; the second case is denoted by a 3 2 edge flip. We change the number of tetrahedra of the five points by +1 or ˆ 1. Obviously, other edge flippings are impossible !

Figure 6.10: The first figure shows the 2 3 edge flip, whereas the second figure shows a 3 2 edge flip. The ground face G is flippable because the endpoints of G, p _i , and q are in convex position. New triangles (shaded) of P( ˜… (p _i )) opposite p _i are found after the flipping.

As shown in Figure 6.10, it might happen that T is in conflict with p _i and the corresponding vertex q opposite p _i belongs to the ground face G of another triangle of P ( ˜… ( p _i )). This is another significant difference between 2D and 3D. In this case, there is an outer tetrahedra with two triangles stemming from P ( ˜… ( p _i )) that has to vanish . Note that P ( ˜… ( p _i )) can be non-convex; see Figure 6.10.

If an edge flip is performed, new ground faces (triangles) in P ( ˜… ( p _i )) might become available. A triangle T in P ( ˜… ( p _i )) is uniquely determined by the fact that p _i is the apex of a tetrahedron in the current triangulation and T lies opposite p _i ; see Figure 6.9.

In principle, the presented algorithm can be implemented in analogy to the 2D case; we store the history of a tetrahedral subdivisions in a DAG. Now, a DAG node has four children for the triangles of the corresponding tetrahedron, a triangle t uniquely separates two tetrahedra T ₁ and T ₂ , and so on. If there are no more flippable tetrahedra beyond the triangles of the current P ( ˜… ( p _i )) opposite p _i , the process stops and DT _i is constructed .

Unfortunately, not every order of edge flips in a single insertion step might finally succeed. Fortunately, [Rajan 91] showed that the flipping process will always result in a current Delaunay triangulation if the flips are done in an appropriate order. ^[1] This result was used for an O ( n ² ) implementation by [Joe 91b].

It was shown in [Guibas et al. 92b] that the expected number of tetrahedra occurring during an incremental construction of the Delaunay triangulation in 3D is bounded by O ( n ² ). A drawback is that the number of final tetrahedra might be in O ( n ), whereas the expected number of edge flips remains to be in O ( n ² ). In [Edelsbrunner and Shah 96], the above results were collected and the randomized incremental construction was extended to regular triangulations.

It remains to show how to compute the flipping order after P ( ˜… ( p _i )) was initialized . Among all flippable tetrahedra opposite p _i having a ground face in P ( ˜… ( p _i )), we compute an order of the candidates by the following elegant idea. The order can be easily computed and incorporated into the DAG algorithm.

The convex hull in ^d ⁺ ¹ and the Delaunay triangulation in defined by empty circle predicates are identical with respect to a lifting transformation; see [Brown 79] or [Edelsbrunner 87]. As will be shown in Section 9.4.2, we can project points p = ( p ₁ , p ₂ , , p _d ) in onto the paraboloid P in ^d ⁺ ¹ by

see Figure 6.11 and Figure 9.19 for d = 2. The lower convex hull of the transformed points on the paraboloid in ^d ⁺ ¹ can be projected back into , and we get the Delaunay triangulation by the empty circle property; see Theorem 9.30 in Section 9.4.2.

Figure 6.11: The projection of points from 2D onto the paraboloid in 3D. The lower convex hull on the paraboloid can be projected back to 2D and represents the Delaunay triangulation in 2D. In general, the Delaunay triangulation in

is equivalent to the convex hull in

^d ⁺¹ .

Let us assume that a new point p _i in is inserted, and let us consider the situation on the convex hull. We use a kinetic interpretation. On the hull, we can assume that the new point comes from the interior of a surface and moves down to its final position on the paraboloid. During this movement, a sequence of triangles of the lower convex hull on the paraboloid becomes visible from the moving point; that is, the moving point lies below the hyperplane passing through the corresponding triangles . These triangles can no longer belong to the lower convex hull; they must vanish successively. They disappear in a fixed order, and this is exactly the order we would like to have for a flipping insertion process. Let us compute this order now.

Every triangle on the lower convex hull lies in a unique hyperplane; see also Section 9.4.2. A triangle on the convex hull becomes visible if the moving point meets the intersection of the triangles hyperplane with the projection line of p _i in dimension d + 1.

In detail in 3D, for the new point p _i = ( , , ) in tetrahedron t the point in 4D moves along a projection line l = {( , , , w ) w ˆˆ } parallel to the fourth axis. The point moves from a point

in a 3D hyperplane, stemming from the projection of the vertices t , to a point

on the paraboloid. Here, C _t = ( , , ) is the circumcenter of the tetrahedron t , and R _t denotes the radius of the circumsphere of t in 3D. Let W denote the fourth coordinate. For the neighboring tetrahedra of P ( ˜… ( p _i )) in 3D, we consider the corresponding projected tetrahedra on the lower convex hull on the paraboloid in 4D. The hyperplane of a projected tetrahedron t has an intersection ( , , , w _t ) with l .

We compute all intersections w _t with the line l for every neighboring tetrahedron in P ( ˜… ( p _i )) and sort them in decreasing order. This gives the flipping order. If , we can neglect the tetrahedron t . The details of the computations can be found in [Rajan 91] and for a tetrahedron t on the boundary of P ( ˜… ( p _i )) with circumcenter C _t = ( , , ) and circumradius R _t , we can show

We can compute the values w _t for every tetrahedron around P ( ˜… ( p _i )); sorting the values gives an appropriate flipping order. Altogether, we get the following result.

Theorem 6.7. The Delaunay triangulation of a set of n points in 3D can be constructed incrementally in expected time O ( n ² ) , using expected quadratic space. The average is taken over the different orders of inserting the n sites.

6.3.2. Constrained Voronoi Diagrams

For many applications, it is not sufficient to consider neighborship relations without any constraints. For example, a small distance between two cities may be unimportant if there is an obstacle (natural or artificial) that blocks their sight to each other. To overcome such problems, the concept of constrained Voronoi diagrams was introduced by [Lee and Lin 86].

We discuss the case of a Voronoi diagram of points in the plane with additional line segment obstacles. Formally, let S be a set of n points in the plane, and let L be a set of non- crossing line segments spanned by S . By induction, it is easy to show that L contains at most 3 n ˆ 6 line segments. The segments in L may be considered as obstacles.

The constrained distance between two points x and y takes visibility information between the points with respect to the segments of L into account. The constrained distance is defined by

where d ( x, y ) denotes the Euclidean distance between x and y . The distance function gives rise to the constrained Voronoi diagram of S and L , ConstrVD( S, L ) for short. Regions of sites that are close but not visible from each other are separated by segments in L ; an example is shown in Figure 6.12.

Figure 6.12: The constrained Voronoi diagram ConstrVD(S, L) and its dual.

The exact dual of VD( S, L ) may no longer be a full triangulation of S (even if S is included). For example, in Figure 6.12, in the dual of VD( S, L ) including L , the face given by p ₁ , p ₁₀ , p ₉ , and p ₅ is not a triangle. The edges p ₅ p ₃ and p ₉ p ₃ of DT( V ) are not part of the dual of ConstrVD( S, L ).

Fortunately, we can modify ConstrVD( S, L ) in order to dualize it into a near to Delaunay triangulation DT ( S, L ), which includes L .

For every line segment l , we proceed as follows. All sites belonging to clipped regions lying to the left of a segment l build a Voronoi diagram to the right of l and vice versa; see Figure 6.12 for an example. The neighborship within these new diagrams leads to additional edges, which are inserted into the dual of ConstrVD( S, L ). For example, in Figure 6.12, in the diagram of p ₁ , p ₅ , p ₉ , and p ₁₀ extended to the right of l , the sites p ₅ and p ₁₀ finally become neighbors, and the corresponding edge in the dual of ConstrVD( S, L ) triangulates the corresponding graph. By construction, the endpoints of a given line segment l are always neighbors in the new diagrams, and therefore, the line segments themselves are inserted.

Algorithms for computing the constrained Voronoi diagram VD( S, L ) and the constrained Delaunay triangulation DT ( S, L ) have been proposed by several authors: Lee and Lin [Lee and Lin 86], Wang and Schubert [Wang and Schubert 87], Chew [Chew 89], and Wang [Wang 93]. Good implementation schemes can be found in [Seidel 88] and [Kao and Mount 92]. An application of DT ( S, L ) to quality mesh generation can be found in [Chew 93].

6.3.3. Types of Generalizations

In this section, we simply list some of the most famous generalization schemes and show examples of some intuitive ones.

First, Voronoi diagrams can be considered in other spaces that use different distance functions . For the L ₁ and the L ₂ metric, similar construction schemes are available. The bisectors of two sites are no longer line segments (see Figure 6.13 for an example), but they consist of at most three segments and divide the plane into distinct half-spaces also. This means that in a brute force manner, we can simply compute the intersections of half-spaces for computing the Voronoi region of a single site. This results in a simple O ( n ³ ) construction algorithm. Optimal O ( n log n ) construction schemes for the L _p norms were considered by several authors; see [Hwang 79, Lee 80, Lee and Wong 80]. Examples of Voronoi diagrams for L ₁ and L ₂ are shown in Figure 6.14.

Figure 6.13: Bisectors under L ₁ and L ₂ distance functions. One can build the bisector by the intersections of scaled copies of the unit circle.

Figure 6.14: Voronoi diagrams under L ₁ and L ₂ distance functions.

More generally , we can consider convex distance functions . A convex distance function is defined by a convex set C with fixed base point O . The distance between two points p and q is given by the following transformation (see also Figure 6.15). Translate C onto p and consider the ray through q starting from p . The ray hits C in a unique point q , and the distance with respect to the convex set C is defined by

Figure 6.15: Bisectors under convex distance functions.

This is also true for the unit circles of L ₁ and L ₂ . The bisector of two points p _i and p _j is the locus of all points with the same distance to p _i and p _j . Therefore, we can use scaled copies of C translated to p _i and p _j . The locus of all intersections of the scaled copies for all scale factors represents the bisector. Therefore, it can be easily shown that for a convex polygon C the bisector of two points has complexity O ( C ).

Furthermore, we can consider different environments . We have already seen that the Voronoi diagram can be generalized to 3D space. Many other environments may be used. For example, we can consider a set of sites on a graph G = ( V, D ). The distance between two arbitrary points on the graph is given by the shortest distance along the graph. An example of the Voronoi diagram of a set of sites on a graph is shown in Figure 6.16. Construction and complexity results for various Voronoi diagrams on trees, and graphs can be found in [Hurtado et al. 04].

Figure 6.16: The Voronoi diagram of a set of sites on a graph.

Beyond metrics and environments, one can generalize the concept of Voronoi diagram with respect to properties of the sites .

Straightforwardly, one might consider more general sites such as line segments or polygonal chains. An example of a Voronoi diagram of a set of line segments is given in Figure 6.17. Many O ( n log n ) construction schemes were presented; see [Yap 87a], [Fortune 87], and [Klein et al. 93]. A simple randomized incremental construction scheme with expected O ( n log n ) running time is presented in [Alt and Schwarzkopf 95]. First, the Voronoi diagram of the endpoints of the line segment is computed. Then the line segment bisectors are inserted successively. The given approach works also for curved objects with properties similar to line segments.

Figure 6.17: The Voronoi diagram of a set of non-intersecting polygonal objects.

Furthermore, every site might have a certain weight that influences the distance function. That is, for two points p _i and p _j with real-valued weights w ( p _i ) and w ( p _j ), the weighted bisector is defined to be the locus of points q so that

if we apply the weights multiplicatively, or

if the weights are meant to be additive. An example of a Voronoi diagram with additive weights is given in Figure 6.18. For additive weights, the bisector is given by the arc of a hyperbola, because the difference of the distances d ( p _i , q ) ˆ d ( p _j , q ) is a constant w ( p _j ) ˆ w ( p _i ). Obviously, a big weight results in a smaller area.

Figure 6.18: The Voronoi diagram with additive weights.

Generalizations of Voronoi diagrams might consider also the intention of the cell subdivision . The normal Voronoi diagram subdivides the plane into cells of the same nearest neighborship. But we can ask also for a subdivision of the farthest neighborship. That is, a cell represents all points that belong to a site farther away than any other site. This subdivision is called the farthest point Voronoi diagram .

A more general concept of the neighborship subdivision is called the kth order Voronoi diagram . A k th order Voronoi diagram in 2D subdivides the plane into regions R dedicated to a set of k sites , , , such that for every point p of a region R , the first k nearest neighbors are , , , . Figure 6.19 shows an example for k = 3. In this setting, the ( n ˆ 1) order Voronoi diagram of n sites represents the farthest point Voronoi diagram.

Figure 6.19: The 3-order Voronoi diagram.

Furthermore, we can make use of colors for separating the point set P into sets of points P ₁ , P ₂ , P _k . Every set P _i has its own color , and holds. The color Voronoi diagram subdivides the plane into regions of nearest neighborship with respect to the sets P _i . The region of P _i contains all points q that have a nearest neighbor stemming from P _i .

Additionally, we can combine the given concepts. For example, the farthest color Voronoi diagram represents a subdivision of the plane into regions with respect to sets of points P ₁ , P ₂ , P _k , so that the region of P _i contains all points q that have a farthest neighbor from P _i . An example of a farthest color Voronoi diagram is given in Figure 6.20 with three colors. Note that the regions are no longer connected and that a set might have an empty region. The diagram has complexity O ( nk ) and can be computed in O ( nk log n ) using methods stemming from the famous book [Sharir and Agarwal 95] on Davenport-Schinzel sequence s; see also [Abellanas et al. 01a] and [Abellanas et al. 01b].

Figure 6.20: The farthest color Voronoi diagram.

^[1] The new site p _i might be in conflict with many tetrahedra along the boundary of P ( ˜… ( p _i )).