Voronoi diagram of orthogonal polyhedrain two and three dimensions

Voronoi diagram of orthogonal polyhedra
in two and three dimensions

Ioannis Z. Emiris Department of Informatics and Telecommunications
National and Kapodistrian University of Athens, Greece ATHENA Research and Innovation Center, Maroussi, Greece
   Christina Katsamaki Department of Informatics and Telecommunications
National and Kapodistrian University of Athens, Greece

Voronoi diagrams are a fundamental geometric data structure for obtaining proximity relations. We consider collections of axis-aligned orthogonal polyhedra in two and three-dimensional space under the max-norm, which is a particularly useful scenario in certain application domains. We construct the exact Voronoi diagram inside an orthogonal polyhedron with holes defined by such polyhedra. Our approach avoids creating full-dimensional elements on the Voronoi diagram and yields a skeletal representation of the input object. We introduce a complete algorithm in 2D and 3D that follows the subdivision paradigm relying on a bounding-volume hierarchy; this is an original approach to the problem. The complexity is adaptive and comparable to that of previous methods, namely linear in the number of sites, namely edges or facets resp. We also provide a numerically stable, open-source implementation in Julia, illustrating the practical nature of our algorithm.

max norm axis-aligned rectilinear straight skeleton subdivision method numeric implementation

1 Introduction

Orthogonal shapes are ubiquitous in numerous applications including raster graphics and VLSI design. Here we focus on the metric which is used in such applications and has been studied much less than . We address 2- and 3-dimensional Voronoi diagrams.

A Voronoi diagram partitions space into regions based on distances to a given set of geometric objects in . Every is a Voronoi site and its Voronoi region under metric , is

The Voronoi diagram is the set , consisting of all points that attain their minimum distance to by at least two Voronoi sites. For general input, the Voronoi diagram is a collection of faces of dimension . A face of dimension comprises points equidistant to at least sites. The union of Voronoi edges and vertices is the 1-skeleton. Equivalently, a Voronoi diagram is defined as the minimization diagram of the distance functions to the sites. The diagram is a partitioning of space into regions, each region consisting of points where some function has lower value than any other function.

Figure 1: Voronoi diagram of a rectilinear polygon with 2 holes.

In this paper, we study Voronoi diagrams in the interior of an orthogonal polyhedron. It may have arbitrarily high genus with holes defined by axis-aligned orthogonal polyhedra, not necessarily convex. The faces of all polyhedra meet at right angles, and their edges are aligned with the axes of a coordinate system. Facets are simply connected (without holes) for simplicity. The sites are the facets on the boundary of all polyhedra.

The distance of two points under is and the distance of to a set is . In Fig. 1, the Voronoi diagram111computed by our software and visualized with Axl viewer. of a rectilinear polygon with 2 holes is shown in blue. Our algorithm follows the Subdivision Paradigm and handles 2D and 3D sites. It reads in a region bounding all input sites and performs a recursive subdivision into cells (using quadtrees or octrees). Then, a reconstruction technique is applied to produce an isomorphic representation of the Voronoi diagram.

1.0.1 Previous Work.

If is the number of polyhedral vertices, the combinatorial complexity of our Voronoi diagrams equals in 2D [13] and in 3D [4]. In 3D, it is estimated experimentally to be, in general, [12].

Related work in 2D concerns Voronoi diagrams of segments. In [13], they introduce an sweep-line algorithm, where is the number of segments; they offer a robust implementation for segments with number of orientations. Another algorithm implemented in library CGAL [6] is incremental. The Voronoi diagram of orthogonal polyhedra (with holes) is addressed in [12] in view of generalizing the sweep-line paradigm to 3D: in 2D it runs in as in [13], and in 3D the sweep-plane version runs in , where is the number of events.

When the diagram is restricted in the interior of a polygon or polyhedron, it serves as a skeletal representation. A skeleton reduces the dimension of the input capturing its boundary’s geometric and topological properties. In particular, straight skeletons are very related to the Voronoi diagram of rectilinear polygons [2]. An algorithm for the straight skeleton of a simple polygon (not necessarily rectilinear) has complexity for fixed [9]. For -monotone rectilinear polygons, a linear time algorithm was recently introduced [7]. In 3D, an analogous equivalence of the straight skeleton of orthogonal polyhedra and the Voronoi diagram exists [4] and a complete analysis of 3D straight skeletons is provided in [3]. Specifically for 3D orthogonal polyhedra, in [4] they offer two algorithms that construct the skeleton in , where is the number of skeleton features. Both algorithms are rather theoretical and follow a wavefront propagation process. Recently, the straight skeleton of a 3D polyhedral terrain was addressed [11].

A Voronoi diagram can contain full-dimensional faces, as part of a bisector. Under , when two points have same coordinate value, their bisector is full dimensional (Fig. 1(a)). Conventions have been adopted, to ensure bisectors between sites are not full-dimensional [13, 12, 6]. We address this issue in the next section. Subdivision algorithms for Voronoi diagrams are numerous, e.g. [5, 8, 15]; our work is related to [15, 5]. These algorithms are quite efficient, since they adapt to the input, and rather simple to implement. None exists for our problem.

1.0.2 Our contribution.

We express the problem by means of the minimization diagram of a set of algebraic functions with restricted domain, that express the distance of points to the boundary. The resulting Voronoi diagram, for general input, is dimensional. We focus on 2D and 3D orthogonal polyhedra with holes, where the resulting Voronoi diagram is equivalent to the straight skeleton. We introduce an efficient and complete algorithm for both dimensions, following the subdivision paradigm which is, to the best of our knowledge, the first subdivision algorithm for this problem. We compute the exact Voronoi diagram (since bisectors are linear). The output data structure can also be used for nearest-site searching.

The overall complexity is output-sensitive, which is a major advantage. Under Hypothesis 1, which captures the expected geometry of the input as opposed to worst-case behaviour, the complexity is in 2D, where the number of sites (edges) and the separation bound (maximum edge length of cells that guarantees termination). This bound is to be juxtaposed to the worst-case bound of of previous methods. In 3D, it is where bounds vertex cardinality per facet (typically constant). Under a further assumption (Remark 1) this bound becomes whereas existing worst-case bounds are quasi-quadratic or cubic in . is measured under appropriate scaling for the bounding box to have edge length 1. Scaling does not affect arithmetic complexity, but may be adapted to reduce the denominators’ size in rational calculations. The algorithm’s relative simplicity has allowed us to develop a numerically stable software in Julia222https://gitlab.inria.fr/ckatsama/L_infinity_Voronoi/, a user-friendly platform for efficient numeric computation; it consists of about 5000 lines of code and is the first open-source code in 3D.

The rest of this paper is organized as follows. The next section provides structural properties of Voronoi diagrams. In Sect. 3 we introduce our 2D algorithm: the 2D and 3D versions share some basic ideas which are discussed in detail in this section. In particular, we describe a hierarchical data structure of bounding volumes, used to accelerate the 2D algorithm for certain inputs and is necessary for the efficiency of the 3D algorithm. Then we provide the complexity analysis of the 2D algorithm. In Sect. 4 we extend our algorithm and analysis to 3D. In Sect. 5 we give some examples and implementation details.

2 Basic definitions and properties

We introduce useful concepts in general dimension. Let be an orthogonal polyhedron of full dimension in dimensions, whose boundary consists of simply connected facets (without holes); these are edges or flats in 2D and 3D, resp. Note that includes the shape’s interior and boundary. Now consists of the closed facets that form the boundary of including all facets of the interior polyhedra. There are as many such polyhedra as the genus. Let denote the Voronoi region of site under the metric. Lem. 2 gives a property preserved by Def. 1.

Figure 2: Voronoi diagrams (in red): (a) standard, under , (b) under Def. 1.

lemmaprop Let . For every point it holds that , where aff is the affine hull of .


Let w.l.o.g. , and . Supposing that , then and there is s.t. . Then, there is at least one site s.t. . Since , then ; contradiction. ∎

Figure 3: , , for segment .

For let be the closed halfspace of induced by aff such that for every there exists a point : and , . We define the (unoriented) definition zone of as . The oriented definition zone of is (Fig. 3). We associate to the distance function

The minimization diagram of restricted to yields a Voronoi partitioning. The Voronoi region of with respect to is

Definition 1.

The Voronoi diagram of with respect to is the set .

This means one gets the Voronoi diagram of Fig. 1(b). Clearly and . The bisector of w.r.t.  is . Then , where affbis denotes the (affine) bisector of . In 2D (resp. 3D) if sites have not the same affine hull, bisectors under lie on lines (resp. planes) parallel to one coordinate axis (resp. planes) or to the bisector of two perpendicular coordinate axes (resp. planes). Although the latter consists of two lines (resp. planes), lies only on one, and it can be uniquely determined by the orientation of the zones. Degeneracy of full-dimensional bisectors, between sites with the same affine hull, is avoided by infinitesimal perturbation of corresponding sites. This is equivalent to assigning priorities to the sites; the full dimesional region of the former diagram is ‘to the limit’ assigned to the site with the highest priority (Fig. 3(b)). Such a perturbation always exists, both for 2D [12, Lem. 13] and 3D [12, Lem. 31].

Figure 4: (a) 2D Voronoi diagram for polygon with colinear edges. (b) 1D Voronoi diagram after infinitesimal perturbation of edges, where .

Set is weakly star shaped with respect to if , such that the segment belongs to .


lemmastar For every , is weakly star shaped with respect to .


Let and , . The open ball centered at with radius is empty of sites. For , let on the line segment . For every it holds that and therefore . Then, and . If there is a site such that and must intersect ; a contradiction. ∎

Therefore, since every is simply connected, from Lem. 2 is simply connected and is also simply connected. Let the degree of a Voronoi vertex be the number of sites to which it is equidistant. If the degree is , the vertex is degenerate. Lem. 2 is nontrivial: in metrics like degree is arbitrarily large. For this bound is tight [12].


lemmadegree (a) The maximum degree of a Voronoi vertex is less than or equal to . (b) When , a Voronoi vertex cannot have degree 7.


(a) Consider the vertex placed at the origin; orthants are formed around the vertex. We count the maximum number of Voronoi regions that can be contained in each orthant and share this Voronoi vertex. Equivalently, we count the number of Voronoi edges that are in the interior of each orthant and have this Voronoi vertex as endpoint; at most one such edge can exist in each orthant. Since these Voronoi edges are equidistant to sites, result follows.

(b) Let be a Voronoi vertex of degree 7. Since 7 Voronoi edges meet at , due to symmetry, we examine the two cases of Fig. 5. When the configuration of Voronoi regions around the vertex is like in Fig. 5a, then is a horizontal segment and are vertical. Then, . Since and is equidistant to both of them, the affine hulls of coincide. Then, whichever is the orientation of , the affine bisectors of and cannot meet like in Fig. 5a. When like in Fig. 5b, since is vertical, is vertical. But since is horizontal, must be horizontal; a contradiction. ∎

Figure 5: The two cases in proof of Lemma 2.

3 Subdivision algorithm in two dimensions

Given manifold rectilinear polygon , i.e. every vertex being shared by exactly two edges, the input consists of and a box bounding . Non-manifold vertices can be trivially converted to manifold with an infinitesimal perturbation. Subdivision algorithms include two phases. First, recursively subdivide to 4 identical cells until certain criteria are satisfied, and the diagram’s topology can be determined in inside each cell. The diagram is reconstructed in the second phase.

3.1 Subdivision Phase

We consider subdivision cells as closed. Given cell , let be the set of sites whose closed Voronoi region intersects : . For point we define its label set . When then The computation of is hereditary, since , if is the parent of . But it is rather costly; given with , it takes to compute , since the relative position of to the bisector of every pair of sites in must be specified. Alternatively, we denote by , the center and the -radius of and define the active set of as:

where , and 0 otherwise. We now explain how approximates by adapting [5, Lem.2], where appears as a soft version of .


lemmayap (a) If then . (b) For a sequence of cells monotonically convergent to point , for large enough.


(a) Let and . It holds that for every . We distinguish two cases according to the position of relatively to . If and , then:

Otherwise, if , since , there is a site intersecting such that . Therefore,

(b) There exists such that for . Therefore, for , since and , for every , (a) implies that . Since , result follows.

One can easily verify , therefore the complexity of computing is linear in the size of . The algorithm proceeds as follows: For each subdivision cell we maintain the label sets of its corners and of its central point, and . The subdivision of a cell stops whenever at least one of the termination criteria below holds (checked in turn). Upon subdivision, we propagate and the label sets of the parent cell to its children. For every child we compute the remaining label sets and refine its active set. Let be the maximum degree of a Voronoi vertex .

3.1.1 Termination criteria:

(T1) for some ; (T2) ; (T3) ; (T4) and the sites in define a unique Voronoi vertex .

When (T1) holds, is contained in a Voronoi region so no part of the diagram is in it. (T2) stops the subdivision when the open cell is completely outside the polygon. If (T3) holds, we determine in the diagram’s topology in since there are Voronoi regions intersected. (T4) stops cell subdivision if it contains a single degenerate Voronoi vertex. The process is summarized in Alg. 1.

1:root bounding box of
2: root
3:while  do
4:      pop
5:     Compute and the label sets of the vertices and the central point.
6:     if (T1) (T2) (T3) (T4) then
7:         return
8:     else
9:         Subdivide into
11:     end if
12:end while
Algorithm 1 Subdivision2D()

theoremthmsub Algorithm 1 halts.


Consider an infinite sequence of boxes such that none of the termination criteria holds. Since (T2) does not hold , the sequence converges to a point . From Lem. 3.1(b), there exists such that . Since , (T4) will hold. ∎


lemmainpred If vertices for some then .


Let . Then, , since is convex in 2D, and . If there is a site s.t. . However, Voronoi regions are simply connected so is on the boundary of . There is a choice of such so that is not parallel to the cell facet on which is. Since lies on the affine bisector and is not a vertex of the cell, there exists in the interior of such that ; a contradiction. ∎

Hence one decides (T1) by checking the vertices’ labels. (T2) is valid for iff and , . Fot (T4), the presence of a Voronoi vertex in is verified through constructor VoronoiVertexTest: given with , the affine bisectors of sites in are intersected. If the intersection point is in and in for every then it is a Voronoi vertex. We do not need to check whether it is in or not; since (T1) fails for , if , there must be intersecting such that : contradiction.

3.2 Reconstruction Phase

We take the quadtree of the subdivision phase and output a planar straight-line graph (PSLG) representing the Voronoi diagram of . is a (vertex) labeled graph and its nodes are of two types: bisector nodes and Voronoi vertex nodes. Bisector nodes span Voronoi edges and are labeled by the two sites to which they are equidistant. Voronoi vertex nodes correspond to Voronoi vertices and so are labeled by at least 3 sites. We visit the leaves of the quadtree and, whenever the Voronoi diagram intersects the cell, bisector or vertex nodes are introduced. By connecting them accordingly with straight-line edges, we obtain the exact Voronoi diagram and not an approximation. We process leaves with that do not satisfy (T1) nor (T2).

3.2.1 Cell with two active sites.

When , intersects or or both. The intersection of with the cell, when non empty, is part of the Voronoi diagram: for each it holds that and . Therefore .


lemmaLtwosites If there is no Voronoi vertex in and for , then .

Since we intersect the affine bisector with the boundary of the cell. An intersection point is in iff . If intersection points are both in , we introduce a bisector node in the middle of the line segment joining them, labeled by . When only one intersection point is in , then must intersect in . Introduce a bisector node at the corner labeled by .

3.2.2 Cell with 3 active sites or more.

Figure 6: A Voronoi vertex node connected with two bisector nodes.

When and the VoronoiVertexTest finds a vertex in or when (a vertex has already been found), we introduce a Voronoi vertex node at the vertex, labeled by corresponding sites. In the presence of corners of in , bisector nodes are introduced and connected to the vertex node (Fig. 6).

If no Voronoi vertex is in , we repeat the procedure described in previous paragraph for each pair of sites. Even if a bisector node is found, it is not inserted at the graph if it is closer to the third site.

3.2.3 Connecting the graph nodes.

The remaining graph edges must cross two subdivision cells. We apply “dual marching cubes” [14] to enumerate pairs of neighboring cells in time linear in the size of the quadtree: cells are neighboring if they share a facet. Let be graph nodes in neighboring cells. We connect them iff:
   are bisector nodes and
   is a bisector node, is a Voronoi vertex node and .
   , are Voronoi vertex nodes, and


[Correctness]theoremreconstruction The output graph is isomorphic to


Let neighboring cells and graph nodes in each of them respectively. If are bisector nodes and , then the line segment is in and on the Voronoi diagram (Lem. 3.2.1). If is a bisector node and is a Voronoi vertex node s.t. , then . If not on the Voronoi diagram then, there is a Voronoi vertex node different than in or ; contradiction. At last, if and are Voronoi vertex nodes, then are adjacent iff their labels differ in two Voronoi sites, say . They are connected as long as coincides with the Voronoi edge equidistant to . ∎

3.3 Primitives, Data-structures, Complexity

Assuming the input vertices are rational, Voronoi vertices are rational [13]. Computing Voronoi vertices, and intersections between affine bisectors and cell facets require linear operations, distance evaluations and comparisons. Therefore, they are exact. The above operations, computing and deciding site-cell intersection are formulated to allow for a direct extension to 3D. In the sequel we discuss design of predicates, computation of label sets and construction of a Bounding Volume Hierarchy. We also provide a complexity analysis of the algorithm.

Figure 7: Test performed by inZone().

Membership in is trivial to decide, thus we focus on predicates that decide membership in . Given and , let the projection of to aff and the interval on centered at with radius . inZone() decides if ; this holds iff (Fig. 7). Given ,, ZoneInCell() decides if . For this evaluation see Lem. 3.3 and Fig. 8.


lemmazoneinc Let , the two facets of parallel to aff, for and Then, iff s.t. .


iff for at least one : Let s.t. and be the projection of on . There exists s.t. . Then, . It holds that iff : Let and its projection on . Then and . We deduce that , since For the opposite, let and in s.t. . Let be its projection on . If we are done. Otherwise, is at distance from equal to , attained at a boundary point . Then, It follows that . ∎

Figure 8: Illustration of test performed by ZoneInCell

To decide if and if , we use isIntersecting() and isStrictlyIntersecting() respectively. Design is trivial. All these predicates are computed in .

3.3.1 Computing label sets.

If then its closest sites are in . Deciding if is done by LocationTest, which identifies position based on the sites that intersect : among these we select those with minimum distance to and for whom inZone() is true. If a convex (resp. concave) corner w.r.t. the interior of is formed by these sites then iff it belongs to the intersection (resp. union) of the oriented definition zones. If no corner is formed or even if is not intersected by any site, decision is trivial. This takes .

3.3.2 Bounding Volume Hierarchy.

We decompose into a collection of rectangles such that any two of them have disjoint interior. We construct a kd-tree on the reflex vertices of the polygon, splitting always at a vertex. An orthogonal polygon with holes, has reflex vertices. The kd-tree subdivides the plane into at most regions. Every terminal region contains a nonempty collection of disjoint rectangles. Let be the maximum number of such rectangles. Using this decomposition, we construct a Bounding Volume Hierarchy (BVH) [10]. It is a tree structure on a set of objects stored at the leaves along with their bounding volume while internal nodes store information of their descendants’ bounding volume. Two important properties are minimal volume and small overlap. They are achieved by using the Axis Aligned Bounding Box (AABB) as bounding shape and building the BVH by bottom-up traversing the constructed kd-tree: at every leaf of the kd-tree we compute the AABB of its rectangles (namely a terminal bounding box) and for every internal node we compute the AABB of its children. The bounding volumes of a node’s children intersect only at their boundary. Space complexity is linear in tree size.

Rectangle-Intersection queries: Given query rectangle the data structure reports all rectangles in the decomposition overlapping with . Starting from the root, for every internal node, we check whether intersects its bounding rectangle or not. In the latter case the data structure reports no rectangles. In the former, we check the position of relative to the bounding boxes of the node’s children so as to decide for each one if it should be traversed or not. We continue similarly: when we reach a terminal bounding box, we check the position of relative to every rectangle in it. Let be the number of terminal bounding boxes intersected by . Following [1] we show:


theorembvh Rectangle intersection queries are answered in .


Let a rectangle-intersection query and an internal node of the BVH tree visited during the query. We distinguish two cases; first, the subtree rooted at contains a terminal bounding box that intersects . There are such nodes at each level. Otherwise, intersects with the bounding rectangle stored at but does not intersect any terminal bounding box of the subtree rooted at . There are at least two such terminal bounding boxes, say and . Since does not intersect there is a line passing through a facet of separating from . Similarly, there exists a line passing through a facet of that separates it from . W.l.o.g. there is a choice of such that and are distinct -if all the terminal bounding boxes of the subtree can be separated by the same line, then cannot intersect -. If are perpendicular, then their intersection also intersects . Since the bounding boxes of each level are strictly non-overlapping, every vertex of intersects a constant number of them (up to 4). So, there is a constant number of such nodes at a given level. When are parallel and no vertex of intersects , then the terminal bounding rectangles of the subtree can be partitioned to those separated by from and to those separated by from . For these distinct sets of terminal bounding boxes to be formed, there must occur a split of by a line parallel and in between . So there is a reflex vertex of the polygon in , causing this split. But ; a contradiction. So there are internal nodes visited at each level of tree. The visited leaf nodes correspond to the terminal bounding boxes that intersect and since each of them encloses at most rectangles, the additional amount of performed operations equals . Summing over all levels of the tree yields a total query complexity of . ∎

Point queries: Given , we report on the rectangles of the decomposition in which lies inside (at most 4 rectangles). When zero, the point lies outside the polygon. Since it is a special case of a rectangle-intersection query, the query time complexity is .

Complexity. Analysis requires a bound on the height of the quadtree. The edge length of the initial bounding box is supposed to be 1 under appropriate scaling. Let separation bound be the maximum value s.t. for every cell of edge length at least one termination criterion holds. Then, the maximum tree height is . Let be the minimum distance of two Voronoi vertices, and the relative thickness of , i.e. the minimum diameter of a maximally inscribed -ball in , where is the complement of .

Lemma 1.

Separation bound is , where the asymptotic notation is used to hide some constants.

Proof. The algorithm mainly subdivides cells that intersect . Most subdivisions occur as long as cells contain adjacent Voronoi vertices, i.e. joined with a Voronoi edge, or non neighboring Voronoi regions are “too close”. In the fist case, let be adjacent Voronoi vertices; it holds that . Consider , w.l.o.g. centered at , and let . Then the maximum cell radius such that is . Therefore, when the minimum cell size (edge length) in the subdivision is , adjacent Voronoi vertices are separated. For the latter case, consider with . When a minimum cell size of is not sufficient for an to not belong in , there is a hole between . Then, minimum cell size is . ∎

This lower bound is tight: in Fig. 8(a) for , and in Fig. 8(b) for . Next we target a realistic complexity analysis rather than worst-case. For this, assume the site distribution in is “sufficiently uniform”. Formally:

(a) Input consists of 28 sites and 164 cells are generated. Total time is 12.0 ms. Minimum cell size is .
(b) Input consists of 12 sites and 64 cells are generated. Total time is 5.6 ms. Minimum cell size is .
Figure 9: The 1-skeleton of the Voronoi diagram is shown in blue.

Hypothesis 1. For subsets , let (resp. ) be the number of sites intersecting (resp. ). We suppose , where denotes the volume of a set in , being the dimension of .

Theorem 3.1.

Under Hypothesis 1 the algorithm complexity is , where is the total number of boundary edges (including any holes).

Proof. At each node, refinement and checking the termination criteria run in time linear in the size of its parent’s active set. At the root . The cardinality of active sets decreases as we move to the lower levels of the quadtree: Let . For cell and , . Let . For a child of and , . Since is empty of sites and may intersect with , we let . We prove that in any combination of it is . Under Hypothesis 1, a cell at tree level has . Computation per tree level, is linear in sum of active sets’ cardinality, therefore summing over all levels of the tree, we conclude. ∎

Queries in the BVH can be used to compute label sets and the active set of a cell. Assume the number of segments touching a rectangle’s boundary is , which is the typical case. Then, we prove the following.

Lemma 2.

We denote by the parent of in the subdivision. Using BVH accelerates the refinement of if .

Proof. A label set is determined by performing a point and a rectagle-intersection query; once the point is detected to lie inside a rectangle of a leaf we find an initial estimation of . Since the closest site to may be on another leaf, we do a rectangle-intersection query centered at with radius . The closest site(s) to are in the intersected leaves. Thus, finding takes (Thm. 3.3.2), where is the number of BVH leaves intersected by . Computing the sites in is accelerated if combined with a rectangle intersection query to find segments at -distance from . Let be the maximum number of BVH leaves intersected by this rectangle-intersection query. We obtain a total refinement time for the cell equal to . Since and the lemma follows. ∎

4 Subdivision algorithm in three dimensions

Let be a manifold orthogonal polyhedron: every edge of is shared by exactly two and every vertex by exactly 3 facets. For non-manifold input we first employ trihedralization of vertices, discussed in [12, 4]. Input consists of and bounding box of . An octree is used to represent the subdivision.

The main difference with the 2D case is that Voronoi sites can be nonconvex. As a consequence, for site , is not necessarily convex and therefore the distance function cannot be computed in : it is not trivial to check membership in . It is direct to extend the 2D algorithm in three dimensions. However, we examine efficiency issues.

For an efficient computation of the basic predicates (of Sect. 3.3), we preprocess every facet of the polyhedron and decompose it to a collection of rectangles. Then a BVH on the rectangles is constructed. The basic operation of all these predicates in 2D is an overlap test between an interval and a segment in 1D. In 3D, the analog is an overlap test between a 2D rectangle and a site (rectilinear polygon). Once the BVH is constructed for each facet, the rectangle-intersection query takes time logarithmic in the number of facet vertices (Thm 3.3.2).

4.0.1 Subdivision.

The active set , and the label set of a point are defined as in to 2D. Most importantly, Lem. 3.1 is valid in 3D as well. The algorithm proceeds as follows: We recursively subdivide into 8 identical cells. The subdivision of a cell stops whenever at least one of the termination criteria below holds. For each cell of the subdivision we maintain the label set of its central point and . Upon subdivision, we propagate from a parent cell to its children for further refinement. We denote by the maximum degree of a Voronoi vertex .

3D Termination criteria: (T1’) , (T2’) , (T3’) and the sites in define a unique Voronoi vertex .

Subdivision is summarized in Alg. 2. (T1’) is valid for iff and it holds that . Detecting a Voronoi vertex in proceeds like in 2D. A Voronoi vertex is equidistant to at least 4 sites and there is a site parallel to each coordinate hyperplane among them.

1:root bounding box of
2: root
3:while  do
4:      pop
5:     Compute the label set of central point and .
6:     if (T1’) (T2’) (T3’) then
7:         return
8:     else
9:         Subdivide into
11:     end if
12:end while
Algorithm 2 Subdivision3D()

theoremthreedhalts Algorithm 2 halts.

(T1) used in 2D is omitted, for it is not efficiently decided: labels of the cell vertices cannot guarantee that . However, when and with , we can prove that . Therefore we can choose a sufficiently small such that . The following lemma gives a lower bound on the maximum such radius:


lemmaLrad Let , and