Voronoi diagram of orthogonal polyhedra in two and three dimensions

TL;DR

This paper presents a novel subdivision-based algorithm for constructing exact Voronoi diagrams of orthogonal polyhedra in 2D and 3D under the max-norm, avoiding full-dimensional elements and providing a skeletal representation.

Contribution

It introduces an original subdivision paradigm with a bounding-volume hierarchy for efficient, exact Voronoi diagram computation of orthogonal polyhedra, including an open-source Julia implementation.

Findings

Algorithm is adaptive and comparable in complexity to previous methods.

Provides a numerically stable, open-source implementation in Julia.

Complexity depends on site number, cell size, and vertex bounds.

Abstract

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. Under a mild assumption it is $O(n / \Delta + 1 / \Delta^2)$ in 2D or $O(n\alpha^2 / \Delta^2 +1 / \Delta^3)$ in 3D, where $n$ is the number of sites, namely edges or facets resp., $\Delta$ is the maximum cell size for the subdivision to stop, and $\alpha$ bounds vertex cardinality per facet. We also provide a numerically stable, open-source implementation in Julia, illustrating the practical nature of our algorithm.

Figures

Equations

$$D_s(\cdot ):{\mathbb{R}}^d\to \mathbb{R}:p\mapsto \{\begin{array}{l}{\mu }_{\infty }(p,s),\text{ if }p\in {\mathcal{Z}}^{+}(s),\\ \infty ,\text{ otherwise.}\end{array}$$

$$V_{\mathcal{D}}(s)=\{p\in \mathcal{P}\mid D_s(p)<\infty \text{ and }\forall s^{\prime }\in \mathcal{S}\setminus s{D}_s(p)<D_{s^{\prime }}(p)\}.$$

$$\widetilde{\phi}(\mathcal{C}):=\big{\{}s\in\mathcal{S}\mid\mathcal{Z}^{+}(s)\cap\mathcal{C}\neq\emptyset,\text{ and }\mu_{\infty}(p_{{}_{\mathcal{C}}},s)\leq 2r_{{}_{\mathcal{C}}}+\delta_{{}_{\mathcal{C}}}\big{\}},$$

$${\mu }_{\infty }(p_{{}_{\mathcal{C}}},s)$$

$$\le 2{\mu }_{\infty }(p_{{}_{\mathcal{C}}},p)+{\mu }_{\infty }(p_{{}_{\mathcal{C}}},s^{\ast })\le 2r_{{}_{\mathcal{C}}}+{\mu }_{\infty }(p_{{}_{\mathcal{C}}},s^{\ast }).$$

$$r=\frac{n}{2}+2(h-1)$$

Full text

1 Introduction

Orthogonal shapes are ubiquitous in numerous applications including raster graphics and VLSI design. We address Voronoi diagrams of 2- and 3-dimensional orthogonal shapes. We focus on the $L_{\infty}$ metric which is used in the relevant applications and has been studied much less than $L_{2}$.

A Voronoi diagram partitions space into regions based on distances to a given set $\mathcal{S}$ of geometric objects in $\mathbb{R}^{d}$. Every $s\in\mathcal{S}$ is a Voronoi site (or simply a site) and its Voronoi region under metric $\mu$, is $V_{\mu}(s)=\{x\in\mathbb{R}^{d}\mid\mu(s,x)<\mu(x,s^{\prime}),\text{ }s^{\prime}\in\mathcal{S}\setminus s\}.$ The Voronoi diagram is the set $\mathcal{V}_{\mu}(\mathcal{S})=\mathbb{R}^{d}\setminus\bigcup_{s\in\mathcal{S}}V_{\mu}(s)$, consisting of all points that attain their minimum distance to $\mathcal{S}$ by at least two Voronoi sites. For general input, the Voronoi diagram is a collection of faces of dimension $0,1,\dots,d-1$. A face of dimension $k$ comprises points equidistant to at least $d+1-k$ sites. Faces of dimension 0 and 1 are called Voronoi vertices and Voronoi edges respectively. 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.

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

The $L_{\infty}$ or Chebyshev distance is a metric defined on a vector space where the distance between two vectors is the greatest of their differences along any coordinate dimension. For two points $x,y\in\mathbb{R}^{d}$, with standard coordinates $x_{i}$ and $y_{i}$ respectively, their $L_{\infty}$ distance is $\mu_{\infty}(x,y)=\max_{i}\{|x_{i}-y_{i}|\}$. The $L_{\infty}$ distance of $x$ to a set $S\subset\mathbb{R}^{d}$ is $\mu_{\infty}(x,S)=\inf\{\mu_{\infty}(x,y)\mid y\in S\}$.

The $\mathbf{L_{\infty}}$ Voronoi diagram uses the $L_{\infty}$ metric to measure distances. In Figure 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 $V$ is the number of polyhedral vertices, the combinatorial complexity of our Voronoi diagrams equals $O(V)$ in 2D [18] and $O(V^{2})$ in 3D [4]. In 3D, it is estimated experimentally to be, in general, $O(V)$ [16].

Related work in 2D concerns $L_{\infty}$ Voronoi diagrams of segments. In [18], they introduce an $O(n\log n)$ sweep-line algorithm, where $n$ is the number of segments; they offer a robust implementation for segments with $O(1)$ number of orientations. Another algorithm implemented in library CGAL [6] is incremental. The $L_{\infty}$ Voronoi diagram of orthogonal polyhedra (with holes) is addressed in [16] in view of generalizing the sweep-line paradigm to 3D: in 2D it runs in $O(n\log n)$ as in [18], and in 3D the sweep-plane version runs in $O(kV)$, where $k=O(V^{2})$ 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 $L_{\infty}$ Voronoi diagram of rectilinear polygons [2]. An algorithm for the straight skeleton of a simple polygon (not necessarily rectilinear) has complexity $O(V^{\frac{17}{11}+\varepsilon})$ for fixed $\varepsilon>0$ [9]. For $x$-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 $L_{\infty}$ 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 $O(\min\{V^{2}\log V,k\log^{O(1)}V\})$, where $k=O(V^{2})$ 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 [13].

For polyhedral objects it is common to consider as Voronoi sites the vertices, edges and higher-dimensional faces (when $d\geq 3$) that form their boundary [21, 10, 14]. However, the combination of the underlying metric and the structure of the input polyhedron may cause the appearance of full-dimensional faces in the Voronoi diagram. Under this event, the Voronoi diagram cannot serve as a skeletal representation of the input shape and the design of an effective algorithm for its computation is significantly more difficult. For example, the $L_{2}$ Voronoi diagram of a polygon in the plane, under the standard definition, contains two dimensional regions of points equidistant to a reflex vertex and its adjacent edges. A usual approach towards eliminating the appearance of these two-dimensional regions is to consider as Voronoi sites the open edges and the polygon vertices and define by convention 1-dimensional bisectors between them. A similar concept is this of defining the ‘cones of influence” [12] or “zones” [22] of each site.

Full-dimensional faces on the Voronoi diagram in our setting are very frequent [16]; under $L_{\infty}$, when two points have same coordinate value, their bisector is full dimensional. Conventions have been adopted, to ensure bisectors between sites are not full-dimensional [18, 16, 6]. We address this issue in the next section.

The literature on subdivision algorithms for Voronoi diagrams is vast, e.g. [5, 8, 22] and references therein. Our work is closely related to [22, 5]. These algorithms are quite efficient, since they adapt to the input, and rather simple to implement. However, 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 $L_{\infty}$ distance of points to the boundary. The resulting Voronoi diagram, for general input, is $(d-1)-$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 $L_{\infty}$ 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 $O(n/\Delta+1/\Delta^{2})$ in 2D, where $n$ the number of sites (edges) and $\Delta$ the separation bound (maximum edge length of cells that guarantees termination). This bound is to be juxtaposed to the worst-case bound of $O(n\log n)$ of previous methods. In 3D, it is $O(n\,\alpha^{2}/\Delta^{2}+1/\Delta^{3})$ where $\alpha$ bounds vertex cardinality per facet (typically constant). Under a further assumption (Remark 2) this bound becomes $O(V/\Delta^{2}+1/\Delta^{3})$ whereas existing worst-case bounds are quasi-quadratic or cubic in $V$. $\Delta$ 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 conclude with some remarks, examples and implementation details.

2 Basic definitions and properties

We introduce useful concepts in general dimension. Let $\mathcal{P}$ be an orthogonal polyhedron of full dimension in $d$ dimensions, whose boundary consists of $n$ simply connected (without holes) facets; these are edges or flats in 2D and 3D, resp. Note that $\mathcal{P}$ includes the shape’s interior and boundary. Now $\mathcal{S}$ consists of the closed facets that form the boundary of $\mathcal{P}$ including all facets of the interior polyhedra. There are as many such polyhedra as the genus.

Under the $L_{\infty}$, the Voronoi diagram $\mathcal{V}_{\infty}(\mathcal{S})$ often contains full-dimensional regions (e.g. Figure 2(a)). Aiming at a $(d-1)$-dimensional Voronoi diagram, we will define an appropriate set of distance functions and restrict their minimization diagram to $\mathcal{P}$.

Let $V_{\infty}(s)$ denote the Voronoi region of site $s$ under the $L_{\infty}$ metric. Lemma 1 gives a property of standard $L_{\infty}$ Voronoi diagram preserved by Def. 1.

Lemma 1 ().

Let $s\in\mathcal{S}$. For every point $p\in V_{\infty}(s)$ it holds that $\mu_{\infty}(p,s)=\mu_{\infty}(p,\text{aff}(s))$, where aff$(s)$ is the affine hull of $s$.

Proof.

Assume w.l.o.g. that $s\subset\{x\in\mathbb{R}^{d}\mid x_{k}=c\}$, $k\in[d]$ and $c\in\mathbb{R}$. If $\mu_{\infty}(p,s)\neq\mu_{\infty}(p,\text{aff}(s))$, then $\mu_{\infty}(p,s)=\inf\{\max_{i\in[d]\setminus k}\{|p_{i}-q_{i}|\}\mid\forall q\in s\}$ and there is $q\in\partial s$ such that $\mu_{\infty}(p,s)=\mu_{\infty}(p,q)=|p_{j}-q_{j}|$, $j\in[d]\setminus k$. To see this, suppose on the contrary that $q$ is in the interior of $s$. Then, we can find $q^{\prime}\in s$ $\varepsilon$-close to $q$ such that $|p_{j}-q_{j}^{\prime}|=|p_{j}-q_{j}|-\varepsilon\Rightarrow\mu_{\infty}(p,q^{\prime})=\mu_{\infty}(p,q)-\varepsilon$, for any $\varepsilon>0$. This leads to a contradiction. Therefore, there is a site $s^{\prime}\neq s$ with $q\in s^{\prime}$. Since $p\in V_{\infty}(s)$, then $\mu_{\infty}(p,s)<\mu_{\infty}(p,s^{\prime})\leq\mu_{\infty}(p,q)$; contradiction. ∎

For $s\in\mathcal{S}$ let $\mathcal{H}(s)$ be the closed halfspace of $\mathbb{R}^{d}$ induced by aff$(s)$ such that for every $p\in s$ there exists a point $q\in\mathcal{H}(s)$ : $q\in int(\mathcal{P})$ and $\mu_{\infty}(p,q)<\epsilon$, $\forall\epsilon>0$. We define the (unoriented) zone of $s$ as $\mathcal{Z}(s):=\{p\in\mathbb{R}^{d}\mid\mu_{\infty}(p,s)=\mu_{\infty}(p,\text{aff}(s))\}$. The oriented zone of $s$ is $\mathcal{Z}^{+}(s):=\mathcal{H}(s)\cap\mathcal{Z}(s)$ (Figure 3). We associate to $s$ the distance function

[TABLE]

The minimization diagram of $\mathcal{D}=\{D_{s}\mid s\in\mathcal{S}\}$ restricted to $\mathcal{P}$ yields a Voronoi partitioning. The Voronoi region of $s$ with respect to $D_{s}(\cdot)$ is

[TABLE]

Definition 1.

The Voronoi diagram of $\mathcal{P}$ w.r.t. $\mathcal{D}$ is $\mathcal{V}_{\mathcal{D}}(\mathcal{P})=\mathcal{P}\setminus\bigcup_{s\in\mathcal{S}}V_{\mathcal{D}}(s)$.

This means one gets the Voronoi diagram of Figure 2(b). Clearly $\mathcal{P}\subset\bigcup_{s\in\mathcal{S}}\mathcal{Z}^{+}(s)$ (Figure 2(a)). Denoting by $\overline{X}$ the closure of a set $X$, then $V_{\infty}(s)\subseteq V_{\mathcal{D}}(s)\subseteq\overline{V_{\infty}(s)}\subseteq\mathcal{Z}^{+}(s)$. The bisector of $s,s^{\prime}\in\mathcal{S}$ w.r.t. $\mathcal{D}$ is $\text{bis}_{\mathcal{D}}(s,s^{\prime})=\{x\in\mathbb{R}^{d}\mid D_{s}(x)=D_{s^{\prime}}(x)<\infty\}$. Then $\text{bis}_{\mathcal{D}}(s,s^{\prime})\subset\text{affbis}(s,s^{\prime})$, where $\text{affbis}(s,s^{\prime})$ denotes the $L_{\infty}$ (affine) bisector of $\text{aff}(s),\text{aff}(s^{\prime})$. In 2D (resp. 3D) if sites have not the same affine hull, bisectors under $\mathcal{D}$ lie on lines (resp. planes) parallel to one coordinate axis (resp. plane) or to the bisector of two perpendicular coordinate axes (resp. planes). Although the latter consists of two lines (resp. planes), $\text{bis}_{\mathcal{D}}$ 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 dimensional region of the former diagram is ‘to the limit’ assigned to the site with the highest priority (Figure 4(b)). Such a perturbation always exists, both for 2D [16, Lem. 13] and 3D [16, Lem. 31].

Set $X$ is weakly star shaped with respect to $Y\subseteq X$ if $\forall x\in X$, $\exists y\in Y$ such that the segment $(x,y)$ belongs to $X$.

Lemma 2 ().

For every $s\in S$, $\overline{V_{\mathcal{D}}(s)}$ is weakly star shaped with respect to $s$.

Proof.

Let $p\in\overline{V_{\mathcal{D}}(s)}$ and $\rho=D_{s}(p)=\mu_{\infty}(p,q)$, $q\in s$. The open ball $B_{\infty}(p,\rho)$ centered at $p$ with radius $\rho$ is empty of sites. For $t\in(0,1)$, let $w=tp+(1-t)q$ on the line segment $(p,q)$. Then, since for every $i\in[d]$ it is $|w_{i}-q_{i}|=t|p_{i}-q_{i}|$, it holds that $w\in\mathcal{Z}^{+}(s)$ and $D_{s}(w)=t\rho$. If $w\not\in\overline{V_{\mathcal{D}}(s)}$ there is a site $s^{\prime}$ such that $D_{s^{\prime}}(w)<t\rho$. But $B_{\infty}(w,t\rho)\subseteq B_{\infty}(p,\rho)$ and $s^{\prime}$ intersects $B_{\infty}(p,\rho)$, leading to a contradiction. ∎

Therefore, since every $s$ is simply connected, from Lemma 2 $\overline{V_{\mathcal{D}}(s)}$ is simply connected and $V_{\mathcal{D}}(s)$ 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 $>d+1$, the vertex is degenerate. Lemma 3 is nontrivial: in metrics like $L_{2}$ degree is arbitrarily large. For $d=2,3$ this bound is tight [16].

Lemma 3 ().

(a) The maximum degree of a Voronoi vertex is less than or equal to $2^{d}d$. (b) When $d=2$, a Voronoi vertex cannot have degree 7.

Proof.

(a) Consider the vertex placed at the origin; $2^{d}$ orthants are formed around the vertex. To obtain the maximum number of Voronoi regions in each orthant, we count the maximum number of Voronoi edges in the interior of an orthant that have this Voronoi vertex as endpoint; at most one such edge can exist in each orthant. Since these Voronoi edges are equidistant to $d$ sites, result follows.

(b) Let $v^{*}=(x^{*},y^{*})$ be a Voronoi vertex of degree 7. Since 7 Voronoi edges meet at $v^{*}$, due to symmetry, we examine the two cases of Figure 5. When the configuration of Voronoi regions around the vertex is like in Figure 5a, then $s_{1}$ is a horizontal segment and $s_{2},s_{7}$ are vertical. Then, $\text{aff}(s_{2}),\text{aff}(s_{7})\subset\{(x,y)\in\mathbb{R}^{2}\mid x>x^{*}\}$. Since $v^{*}\in\mathcal{Z}^{+}(s_{2})\cap\mathcal{Z}^{+}(s_{7})$ and is equidistant to both $s_{2}$ and $s_{7}$, the affine hulls of $s_{2},s_{7}$ coincide. Then, whichever is the orientation of $s_{1}$, the affine bisectors of $s_{1},s_{2}$ and $s_{1},s_{7}$ cannot meet like in Figure 5a. When like in Figure 5b, since $b_{3}$ is vertical, $s_{1}$ is vertical. But since $b_{1}$ is horizontal, $s_{1}$ must be horizontal; a contradiction. ∎

3 Subdivision algorithm in two dimensions

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

3.1 Subdivision Phase

We consider subdivision cells as closed. Given cell $\mathcal{C}$, let $\phi(\mathcal{C})$ be the set of sites whose closed Voronoi region intersects $\mathcal{C}$: \phi(\mathcal{C})=\big{\{}s\in\mathcal{S}\mid\overline{V_{\mathcal{D}}(s)}\cap\mathcal{C}\neq\emptyset\big{\}}. For point $p\in\mathcal{P}$ we define its label set $\lambda(p)=\{s\in\mathcal{S}\mid p\in\overline{V_{\mathcal{D}}(s)}\}$. When $p\in\mathcal{P}^{c}$, where $\mathcal{P}^{c}$ is the complement of $\mathcal{P}$, then $\lambda(p)=\emptyset.$ Using the definition of label sets, one can alternatively write $\phi(\mathcal{C})$ as $\phi(\mathcal{C})=\bigcup_{p\in\mathcal{C}}\lambda(p)$.

Intuitively, storing $\phi(\mathcal{C})$ for every cell of the subdivision would ensure that upon termination and when the cardinality of this set for every leaf cell is a small constant, we could determine the topology of the Voronoi diagram inside each of them in constant time. The computation of $\phi(\mathcal{C})$ is hereditary, since $\phi(\mathcal{C})\subseteq\phi(\mathcal{C}^{\prime})$, if $\mathcal{C}^{\prime}$ is the parent of $\mathcal{C}$. But it is rather costly; given $\phi(\mathcal{C}^{\prime})$ with $|\phi(\mathcal{C}^{\prime})|=\kappa$, it takes $O(\kappa^{2})$ to compute $\phi(\mathcal{C})$, since the relative position of $\mathcal{C}$ to the bisector of every pair of sites in $\phi(\mathcal{C}^{\prime})$ must be specified. Alternatively, following the work of [22, 5], instead of implementing the exact predicate $\phi(\cdot)$, we compute an approximate one. We denote by $p_{{}_{\mathcal{C}}}$, $r_{{}_{\mathcal{C}}}$ the center and the $L_{\infty}$-radius of $\mathcal{C}$ and define the active set of $\mathcal{C}$ as:

[TABLE]

where $\delta_{{}_{\mathcal{C}}}=\min_{s}D_{s}(p_{{}_{\mathcal{C}}})\text{, if }p_{{}_{\mathcal{C}}}\in\mathcal{P}$, and 0 otherwise. We now explain how $\widetilde{\phi}$ approximates $\phi$ by adapting [5, Lem.2], where $\widetilde{\phi}$ appears as a soft version of $\phi$.

Lemma 4 ().

(a) For every cell $\mathcal{C}$, $\phi(\mathcal{C})\subseteq\widetilde{\phi}(\mathcal{C})$. (b) For a sequence of cells $(\mathcal{C})_{i}$ monotonically convergent to point $p\in\mathcal{P}$, $\widetilde{\phi}(\mathcal{C}_{i})=\phi(p)$ for $i$ large enough.

Proof.

(a) If $\phi(\mathcal{C})=\emptyset$, assertion follows trivially. Let $s\in\phi(\mathcal{C})$ and $p\in\mathcal{C}\cap\overline{V_{D}(s)}$. It holds that $D_{s}(p)\leq D_{s^{\prime}}(p)\Rightarrow\mu_{\infty}(p,s)\leq\mu_{\infty}(p,s^{\prime})$ for every $s^{\prime}\in\mathcal{S}$. We distinguish two cases according to the position of $p_{{}_{\mathcal{C}}}$ relatively to $\mathcal{P}$. If $p_{{}_{\mathcal{C}}}\in\mathcal{P}$ and $p_{{}_{\mathcal{C}}}\in\overline{V_{\mathcal{D}}(s^{*})}$, then:

[TABLE]

Otherwise, if $p_{{}_{\mathcal{C}}}\not\in\mathcal{P}$, since $\mathcal{C}\cap\mathcal{P}\neq\emptyset$, there is a site $s^{\prime}$ intersecting $\mathcal{C}$ such that $\mu_{\infty}(p,s^{\prime})\leq r_{{}_{\mathcal{C}}}$. Therefore, $\mu_{\infty}(p_{{}_{\mathcal{C}}},s)\leq\mu_{\infty}(p_{{}_{\mathcal{C}}},p)+\mu_{\infty}(p,s)\leq\mu_{\infty}(p_{{}_{\mathcal{C}}},p)+\mu_{\infty}(p,s^{\prime})\leq 2r_{{}_{\mathcal{C}}}.$

(b) There exists $i_{0}\in\mathbb{N}$ such that for $i\geq i_{0}$ $\mathcal{C}_{i}\cap\mathcal{P}\neq\emptyset$. Therefore, for $i\gg i_{0}$, since $p_{{}_{\mathcal{C}_{i}}}\rightarrow p$ and $r_{{}_{\mathcal{C}_{i}}}\rightarrow 0$, for every $s\in\mathcal{C}_{i}$, (a) implies that $s\in\lambda(p)=\phi(p)$. Since $\phi(p)\subseteq\widetilde{\phi}(\mathcal{C}_{i})$, result follows.

∎

The previous lemma is crucial in proving correctness of the algorithm. It is also worth mentioning that a preliminary version of $\widetilde{\phi}$ first appeared in [17].

One can easily verify $\widetilde{\phi}(\mathcal{C})\subseteq\widetilde{\phi}(\mathcal{C}^{\prime})$, therefore the complexity of computing $\widetilde{\phi}(\mathcal{C})$ is linear in the size of $\widetilde{\phi}(\mathcal{C}^{\prime})$. The algorithm proceeds as follows: For each subdivision cell we maintain the label sets of its corner points and of its central point, and $\widetilde{\phi}$. The subdivision of a cell stops whenever at least one of the termination criteria below holds (checked in turn). Upon subdivision, we propagate $\widetilde{\phi}$ 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 $M$ be the maximum degree of a Voronoi vertex $(M\leq 8)$.

3.1.1 Termination criteria:

(T1) $\mathcal{C}\subseteq V_{\mathcal{D}}(s)$ for some $s\in\mathcal{S}$; (T2) $int(\mathcal{C})\cap\mathcal{P}=\emptyset$; (T3) $|\widetilde{\phi}(\mathcal{C})|\leq 3$; (T4) $|\widetilde{\phi}(\mathcal{C})|\leq M$ and the sites in $\widetilde{\phi}(\mathcal{C})$ define a unique Voronoi vertex $v\in\mathcal{C}$.

When (T1) holds, $\mathcal{C}$ 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 $O(1)$ time the diagram’s topology in $\mathcal{C}$ since there are $\leq 3$ Voronoi regions intersected. (T4) stops cell subdivision if it contains a single degenerate Voronoi vertex. The process is summarized in Alg. 1.

Theorem 3.1 ().

Algorithm 1 halts.

Proof.

Consider an infinite sequence of boxes $\mathcal{C}_{1}\supseteq\mathcal{C}_{2}\supseteq\dots$ such that none of the termination criteria holds. Since (T1) and (T2) do not hold for any $\mathcal{C}_{i}$ with $i\geq 1$, the sequence converges to a point $p\in\mathcal{V}_{\mathcal{D}}(\mathcal{P})$. From Lemma 4(b), there exists $i_{0}\in\mathbb{N}$ such that $\widetilde{\phi}(\mathcal{C}_{i_{0}})=\phi(p)=\lambda(p)$. Since $|\lambda(p)|\leq 8$, (T4) will hold. ∎

Lemma 5 ().

For a subdivision cell $\mathcal{C}$, let $v_{1},\dots,v_{4}$ its corner vertices. For $s\in\mathcal{S}$, $\mathcal{C}\subseteq V_{\mathcal{D}}(s)$ if and only if $v_{1},\dots,v_{4}\in V_{\mathcal{D}}(s)$.

Proof.

Let $v_{1},\dots,v_{4}\in V_{\mathcal{D}}(s)$ and $p\in\mathcal{C}$. Then, $p\in\mathcal{Z}^{+}(s)$, since $\mathcal{Z}^{+}(s)$ is convex in 2D and $v_{1},\dots,v_{4}\in\mathcal{Z}^{+}(s)$. For $i=1,\dots,4$ the open ball $B_{i}:=B_{\infty}(v_{i},\mu_{\infty}(v_{i},s))$ is empty of sites. Since $B_{\infty}(p,\mu_{\infty}(p,\text{aff}(s)))\subset\cup_{i\in[4]}B_{i}$ it holds that $\mu_{\infty}(p,\mathcal{P})\geq\mu_{\infty}(p,\text{aff}(s))=\mu_{\infty}(p,s)$. On the other hand, $\mu_{\infty}(p,\mathcal{P})\leq\mu_{\infty}(p,s)$. So, if $p\not\in V_{\mathcal{D}}(s)$ there is a site $s^{\prime}$ s.t. $D_{s^{\prime}}(p)=D_{s}(p)$ and $p\in\mathcal{V}_{\mathcal{D}}(\mathcal{P})$. Therefore, since Voronoi regions are simply connected and Voronoi edges are straight lines, $p$ must be on the boundary of $\mathcal{C}$. The two possible configurations are shown in Figure 6(a) and are contradictory; for the first, use an argument similar to that of Lemma 3(b). For the second, notice that this cannot hold since the cell is square. We conclude that $\mathcal{C}\subseteq V_{\mathcal{D}}(s)$. The other direction is trivial. ∎

Lemma 6.

For a subdivision cell $\mathcal{C}$ it holds that $int(\mathcal{C})\cap\mathcal{P}=\emptyset$ if and only if $\lambda(p_{{}_{\mathcal{C}}})=\emptyset$ and for every $s\in\widetilde{\phi}(C)$ it holds that $s\cap int(C)=\emptyset$.

Proof.

Suppose that $\lambda(p_{{}_{\mathcal{C}}})=\emptyset$, i.e. $p_{{}_{\mathcal{C}}}\not\in\mathcal{P}$, and that $s\cap int(C)=\emptyset$, for every $s\in\widetilde{\phi}(C)$. If is $int(\mathcal{C})\cap\mathcal{P}\neq\emptyset$, since the center of the cell is not in $\mathcal{P}$, a segment (site) on the boundary of $\mathcal{P}$ must intersect $int(\mathcal{C})$. This site belongs to $\widetilde{\phi}(C)$, leading to a contradiction. Proof of the other direction is trivial. ∎

Hence one decides (T1) by checking the vertices’ labels. (T2) is valid for $\mathcal{C}$ iff $\lambda(p_{{}_{\mathcal{C}}})=\emptyset$ and $\forall s\in\widetilde{\phi}(\mathcal{C})$, $s\cap int(\mathcal{C})=\emptyset$. Fot (T4), the presence of a Voronoi vertex in $\mathcal{C}$ is verified through constructor VoronoiVertexTest: given $\mathcal{C}$ with $|\widetilde{\phi}(\mathcal{C})|\geq 3$, the affine bisectors of sites in $\widetilde{\phi}(\mathcal{C})$ are intersected. If the intersection point is in $\mathcal{C}$ and in $\mathcal{Z}^{+}(s)$ for every $s\in\widetilde{\phi}(\mathcal{C})$ then it is a Voronoi vertex. We do not need to check whether it is in $\mathcal{P}$ or not; since (T1) fails for $\mathcal{C}$, if $v\not\in\mathcal{P}$, there must be $s$ intersecting $\mathcal{C}$ such that $v\not\in\mathcal{Z}^{+}(s)$: contradiction.

3.2 Reconstruction Phase

We take the quadtree of the subdivision phase and output a planar straight-line graph (PSLG) $G=(V,E)$ representing the Voronoi diagram of $\mathcal{P}$. $G$ 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.

The reconstruction can be briefly described as follows: 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 $|\widetilde{\phi}(\cdot)|\geq 2$ that do not satisfy (T1) nor (T2).

3.2.1 Cell with two active sites.

When $\widetilde{\phi}(\mathcal{C})=\{s_{1},s_{2}\}$, $\mathcal{C}$ intersects $V_{\mathcal{D}}(s_{1})$ or $V_{\mathcal{D}}(s_{2})$ or both; at this point $\mathcal{C}$ cannot be wholly contained in $V_{D}(s_{1})$ or $V_{D}(s_{2})$, since (T1) is not satisfied. The intersection of $\text{bis}_{\mathcal{D}}(s_{1},s_{2})$ with the cell, when non empty, is part of the Voronoi diagram: for each $p\in\text{bis}_{\mathcal{D}}(s_{1},s_{2})\cap\mathcal{C}$ it holds that $D_{s_{1}}(p)=D_{s_{2}}(p)$ and $\lambda(p)\subseteq\widetilde{\phi}(\mathcal{C})=\{s_{1},s_{2}\}$. Therefore $p\in\overline{V_{\mathcal{D}}(s_{1})}\cap\overline{V_{\mathcal{D}}(s_{2})}$.

Remark 1.

If there is no Voronoi vertex in $\mathcal{C}$ and $p_{1},p_{2}\in\text{bis}_{\mathcal{D}}(s_{1},s_{2})\cap\mathcal{C}$ for $s_{1},s_{2}\in\widetilde{\phi}(\mathcal{C})$, then $p_{1}p_{2}\subset\text{bis}_{\mathcal{D}}(s_{1},s_{2})$.

Since $\text{bis}_{\mathcal{D}}(s_{1},s_{2})\subset\text{affbis}(s_{1},s_{2})$ we intersect the affine bisector with the boundary of the cell. An intersection point $p\in bis_{\infty}(\text{aff}(s_{1}),\text{aff}(s_{2}))$ is in $bis_{\mathcal{D}}(s_{1},s_{2})$ iff $p\in\mathcal{Z}^{+}(s_{1})\cap\mathcal{Z}^{+}(s_{2})$ (Lemma 1). If intersection points are both in $\mathcal{Z}^{+}(s_{1})\cap\mathcal{Z}^{+}(s_{2})$ (Figure 7(a)), we introduce a bisector node in the middle of the line segment joining them, labeled by $\{s_{1},s_{2}\}$. When only one intersection point is in $\mathcal{Z}^{+}(s_{1})\cap\mathcal{Z}^{+}(s_{2})$, then $s_{1},s_{2}$ must intersect in $\mathcal{C}$ (Figure 7(b)). We insert a bisector node at their intersection point labeled by $\{s_{1},s_{2}\}$.

3.2.2 Cell with 3 active sites or more.

When $|\widetilde{\phi}(\mathcal{C})|=3$ and the VoronoiVertexTest finds a vertex in $\mathcal{C}$ or when $|\widetilde{\phi}(\mathcal{C})|\geq 4$ (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 $\mathcal{P}$ in $\mathcal{C}$, bisector nodes are introduced and connected to the vertex node (Figure 8).

If no Voronoi vertex is in $\mathcal{C}$, 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.

Once the graph nodes are fixed they have to be connected appropriately. The only graph edges added so far are those that are completely contained in a subdivision cell. The remaining graph edges must cross two subdivision cells. We apply “dual marching cubes” [19] to enumerate pairs of neighboring cells in time linear in the size of the quadtree: cells are neighboring if they share a facet.

The traversal method involves two recursive functions: the faceProc and the edgeProc (see Figure 10). At first we call the faceProc for the bounding box of the input polygon. This function, for every internal node of the quadtree, recursively calls itself for each of its children and the edgeProc for every pair of neighboring children. When the faceProc reaches a leaf it terminates. When the edgeProc reaches two adjacent leaves, then the corresponding cells share a facet. Let $v_{1},v_{2}$ be graph nodes in neighboring cells. We connect them if and only if:

• •

$v_{1},v_{2}$ are bisector nodes and $\lambda(v_{1})=\lambda(v_{2}).$

• •

$v_{1}$ is a bisector node, $v_{2}$ is a Voronoi vertex node and $\lambda(v_{1})\subset\lambda(v_{2})$.

• •

$v_{1}$, $v_{2}$ are Voronoi vertex nodes, $\lambda(v_{1})\cap\lambda(v_{2})=\{s,s^{\prime}\}$ and $v_{1}v_{2}\subset\mathcal{P}$.

See Figure 9 for an example where $v_{1},v_{2}$ are Voronoi vertex nodes with $\lambda(v_{1})\cap\lambda(v_{2})=\{s,s^{\prime}\}$ and $v_{1}v_{2}\not\subset\mathcal{P}$.

Theorem 3.2 (Correctness).

The output graph is isomorphic to $\mathcal{V}_{\mathcal{D}}(\mathcal{P}).$

Proof.

We need to prove that the nodes in the graph are connected correctly. Let neighboring cells $\mathcal{C}_{1},\mathcal{C}_{2}$ and $v_{1},v_{2}$ graph nodes in each of them respectively. If $v_{1},v_{2}$ are bisector nodes and $\lambda(v_{1})=\lambda(v_{2})$, then the line segment $v_{1}v_{2}$ is in $bis_{{}_{\mathcal{D}}}(s_{1},s_{2})$, for $s_{1},s_{2}\in\lambda(v_{1})$, and on the Voronoi diagram (Rem. 1). If $v_{1}$ is a bisector node and $v_{2}$ is a Voronoi vertex node s.t. $\lambda(v_{1})\subseteq\lambda(v_{2})$, then $v_{1}v_{2}\subset\text{bis}_{\infty}(s_{1},s_{2})$. If the segment $v_{1}v_{2}$ is not on the Voronoi diagram then, there is a Voronoi vertex node different than $v_{2}$ in $\mathcal{C}_{1}$ or $\mathcal{C}_{2}$; contradiction. At last, let $v_{1}$ and $v_{2}$ be Voronoi vertex nodes such that their labels have two sites in common, say $s,s^{\prime}$, and the edge $v_{1}v_{2}\subset\mathcal{P}$. Vertices $v_{1},v_{2}$ are both on the boundary of $\overline{V_{\mathcal{D}}(s)}\cap\overline{V_{\mathcal{D}}(s^{\prime})}$. Since $v_{1}v_{2}\subset\mathcal{P}$, if it does not coincide with the Voronoi edge equidistant to $s,s^{\prime}$, then both $v_{1},v_{2}$ must also be on the boundary of a Voronoi region other than $\overline{V_{\mathcal{D}}(s)}$ and $\overline{V_{\mathcal{D}}(s^{\prime})}$. This leads to a contradiction. ∎

3.3 Primitives, Data-structures, Complexity

Assuming the input vertices are rational, Voronoi vertices are rational [18]. 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 $\widetilde{\phi}$ 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.

3.3.1 Primitives.

Membership in $\mathcal{H}(s)$ is trivial to decide, thus we focus on predicates that decide membership in $\mathcal{Z}(s)$. Given $p\in\mathbb{R}^{2}$ and $s\in\mathcal{S}$, let $pr_{\text{aff}(s)}(p)$ the projection of $p$ to aff$(s)$ and $I_{p,s}$ the $1d-$interval on $\text{aff}(s)$ centered at $pr_{\text{aff}(s)}(p)$ with radius $\mu_{\infty}(p,\text{aff}(s))$. inZone($p,s$) decides if $p\in\mathcal{Z}(s)$; this holds if and only if $I_{p,s}\cap s\neq\emptyset$ (Figure 11).

Given $\mathcal{C}$,$s\in\mathcal{S}$, ZoneInCell($s,\mathcal{C}$) decides if $\mathcal{Z}(s)\cap\mathcal{C}\neq\emptyset$. For this evaluation see Lemmma 7 and Figure 12.

Lemma 7.

Let $s\in\mathcal{S}$, $f_{1},f_{2}$ the two facets of $\mathcal{C}$ parallel to aff$(s)$, $\rho_{i}=\mu_{\infty}(f_{i},\text{aff}(s))$ for $i=1,2$ and $p_{{}_{\mathcal{C}}}^{\prime}=pr_{\text{aff}(s)}(p_{{}_{\mathcal{C}}}).$ Then, $\mathcal{Z}(s)\cap\mathcal{C}\neq\emptyset$ iff $\exists i\in\{1,2\}$ s.t. $B_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},r_{{}_{\mathcal{C}}}+\rho_{i})\cap s\neq\emptyset$.

Proof.

$\mathcal{Z}(s)\cap\mathcal{C}\neq\emptyset$ iff $\mathcal{Z}(s)\cap f_{i}\neq\emptyset$ for at least one $i\in\{1,2\}$: Let $p\in\mathcal{Z}(s)\cap C$ s.t. $p\not\in f_{1}\cup f_{2}$ and $\text{pr}_{f_{i}}(p)$ be the projection of $p$ on $f_{i}$. There exists $i\in\{1,2\}$ s.t. $\mu_{\infty}(\text{pr}_{f_{i}}(p),\text{aff}(s))>\mu_{\infty}(p,\text{aff}(s))$. Then, $\text{pr}_{f_{i}}(p)\in\mathcal{Z}(s)$. It holds that $\mathcal{Z}(s)\cap f_{i}\neq\emptyset$ iff $B_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},r_{{}_{\mathcal{C}}}+\rho_{i})\cap s\neq\emptyset$: Let $q\in\mathcal{Z}(s)\cap f_{i}$ and $q^{\prime}$ its projection on $\text{aff}(s)$. Then $\mu_{\infty}(q^{\prime},s)\leq\mu_{\infty}(q,\text{aff}(s))=\rho_{i}$ and $\mu_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},q^{\prime})\leq r_{{}_{\mathcal{C}}}$. We deduce that $B_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},r_{{}_{\mathcal{C}}}+\rho_{i})\cap s\neq\emptyset$, since $\mu_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},s)\leq\mu_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},q^{\prime})+\mu_{\infty}(q^{\prime},s)\leq r_{{}_{\mathcal{C}}}+\rho_{i}.$ For the inverse direction, let $B_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},r_{{}_{\mathcal{C}}}+\rho_{i})\cap s\neq\emptyset$ and $q^{\prime}$ in $s$ s.t. $\mu_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},q^{\prime})\leq r_{{}_{\mathcal{C}}}+\rho_{i}$. Let $q$ be its projection on $\text{aff}(f_{i})$. If $q\in f_{i}$ we are done. Otherwise, $q$ is at $L_{\infty}$ distance from $f_{i}$ equal to $\mu_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},q^{\prime})-r_{{}_{\mathcal{C}}}$, attained at a boundary point $q^{\prime\prime}\in f_{i}$. Then, $\rho_{i}\leq\mu_{\infty}(q^{\prime\prime},s)\leq\mu_{\infty}(q^{\prime\prime},q^{\prime})=\max\{\rho_{i},\mu_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},q^{\prime})-r_{{}_{\mathcal{C}}}\}=\rho_{i}.$ It follows that $q^{\prime\prime}\in\mathcal{Z}(s)$.∎

To decide if $s\cap\mathcal{C}\neq\emptyset$ and if $s\cap int(\mathcal{C})\neq\emptyset$, we use isIntersecting($s,\mathcal{C}$) and isStrictlyIntersecting($s,\mathcal{C}$) respectively. Design is trivial. All these predicates are computed in $O(1)$ time.

Deciding whether a site belongs to $\widetilde{\phi}(\mathcal{C})$, is done using the predicates:

• •

isIntersecting: to decide whether the site intersects the cell centered at $p_{{}_{\mathcal{C}}}$ with radius $2\cdot r_{{}_{\mathcal{C}}}+\delta_{{}_{\mathcal{C}}}$.

• •

ZoneInCell: to decide whether the zone of the site intersects $\mathcal{C}$.

Computing label sets.

If $p\in\mathcal{P}\cap\mathcal{C}$ then its closest sites are in $\widetilde{\phi}(\mathcal{C})$. Deciding if $p\in\mathcal{P}$ is done by LocationTest, which identifies position based on the sites that intersect $\mathcal{C}$: among these we select those with minimum $L_{\infty}$ distance to $p$ and for whom inZone($p,s$) is true. If a convex (resp. concave) corner with respect to the interior of $\mathcal{P}$ is formed by these sites then $p\in\mathcal{P}$ iff it belongs to the intersection (resp. union) of the oriented zones. If no corner is formed or even if $\mathcal{C}$ is not intersected by any site, decision is trivial. This takes $O(|\widetilde{\phi}(\mathcal{C})|)$ time. In more details, we define:

• •

$T:=\{s\in\mathcal{S}\mid s\cap\mathcal{C}\neq\emptyset\}$

• •

$d(p,T):=\min\{\mu_{\infty}(p,s)\mid s\in T\text{ and }p\in\mathcal{Z}(s)\}$

• •

$R(p,T):=\{s\in T\mid\mu_{\infty}(p,s)=d(p)\text{ and }p\in\mathcal{Z}(s)\}$

• •

$R^{\prime}(p,s,T)=\{s^{\prime}\in R(p,T)\mid s\cap s^{\prime}\neq\emptyset\}$, $s\in R(p,T)$.

If $T=\emptyset$ then the cell is contained in $\mathcal{P}$. Otherwise, even if $p\not\in\mathcal{P}$ there always exists a site $s$ intersecting $\mathcal{C}$ with $p\in\mathcal{Z}(s)$. Therefore $d(p,T)$ is well defined and $|R(p,T)|\geq 1$ for every $p\in\mathcal{C}$. We pick a site $s\in R(p,T)$ and then perform the LocationTest for $R^{\prime}(p,s,T)$ that returns true if and only if $p\in\mathcal{P}$; this test identifies the position of $p$ relative to $\mathcal{P}$ using one site if $|R^{\prime}(p,s,T)|=1$ or a corner formed by two sites in $R^{\prime}(p,s,T)$ otherwise. Note that $R^{\prime}(p,s,T)$ may contain 1 or 2 sites.

Lemma 8.

Let $p\in\mathcal{C}$. The following equivalence holds: $p\in\mathcal{P}\text{ if and only if }T=\emptyset\text{ or {LocationTest}}(p,R^{\prime}(p,s,T))\text{ returns {true} for some }s\in R(p,T)$

Proof.

When $T=\emptyset$ the equivalence is obvious. When $T\neq\emptyset$, we distinguish two cases according to the cardinality of $R^{\prime}(p,s,T)$: When $|R^{\prime}(p,s,T)|=1$ it is straightforward that $p\in\mathcal{P}\Leftrightarrow p\in\mathcal{Z}^{+}(s)$. When $|R^{\prime}(p,s,T)|=2$, the configuration of sites in $R^{\prime}(p,s)$ is like in Figures 13(a) and 13(b). The test returns true if and only if $p\in\mathcal{P}$ in any case. ∎

3.3.2 Rectangular Decomposition and Bounding Volume Hierarchy.

We decompose $\mathcal{P}$ into a collection of rectangles such that any two of them have disjoint interior. The resulting decomposition may not be optimal with respect to the number of rectangles but allows for an immediate construction of a data structure (a Bounding Volume Hierarchy) that can answer point and rectangle-intersection queries on the decomposition’s rectangles efficiently. This data structure can accelerate the 2D algorithm for certain inputs and is essential for efficiency of the 3D version of the algorithm.

Rectangular decomposition.

It is known [20] that the minimum number of rectangles in a partition of a polygon with $n$ vertices and $h$ holes is $n/2+h-g-1$, where $g$ is the maximum number of nonintersecting chords that can be drawn either horizontally or vertically between reflex vertices. To the direction of minimizing the number of rectangles, the optimal algorithm has $O(n^{3/2}\log n)$ time complexity [15].

However, for our purpose, the decomposition does not need to be optimal. We observe that drawing an axis-parallel edge, from every reflex vertex results to a rectangular decomposition. Thus, we construct a kd-tree on the reflex vertices of the polygon, splitting always at a vertex. Assuming the polygon has $r$ reflex vertices, the kd-tree subdivides the plane into at most $r+1$ regions. Every terminal region contains a disjoint collection of rectangles (nonempty). We denote by $t$ the maximum number of rectangles in a terminal region.

Lemma 9.

For an orthogonal polygon on $n$ vertices with $h$ holes, let $r$ be the number of its reflex vertices. It holds that:

[TABLE]

Proof.

A simple orthogonal polygon without holes on $N$ vertices has $N/2-2$ reflex vertices. So, if we denote by $n_{0}$ the number of the outer contour’s vertices, and by $n_{i}$ the respective number for every hole ($i\in[h]$) we have that: r=n_{0}/2-2+\sum_{i=1}^{h}\big{(}n_{i}/2+2\big{)}=\frac{n}{2}+2(h-1). ∎

Bounding Volume Hierarchy.

A Bounding Volume Hierarchy (BVH) [11] 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.

In our setting, the geometric objects are the axis aligned rectangles obtained by the decomposition. Consequently, the most appropriate bounding shape is the Axis Aligned Bounding Box (AABB) offering a good tradeoff between minimal volume and simplicity of representation. As for minimizing the overlap of bounding boxes at the same level in the hierarchy, we group rectangles in a specific way in order to accomplish that:

The BVH is built in a bottom-up manner, by traversing the kd-tree previously constructed and adding some additional information to its nodes. 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 find the AABB of its two children. In that way, the bounding volumes of a node’s children intersect only at their boundary (see Figure 14 for a simple illustration). Space complexity is linear in tree size.

Rectangle-Intersection queries: Given query rectangle $Q$ the data structure reports all rectangles in the decomposition overlapping with $Q$. Starting from the root, for every internal node, we check whether $Q$ intersects its bounding rectangle or not. In the latter case the data structure reports no rectangles. In the former, we check the position of $Q$ 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 $Q$ relative to every rectangle in it. Let $k$ be the number of terminal bounding boxes intersected by $Q$. Following [1], we count the number of internal nodes visited on each level of the tree and show:

Theorem 3.3 ().

Rectangle intersection queries are answered in $O(k\lg r+kt)$.

Proof.

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

Point queries: Given $p\in\mathbb{R}^{2}$, we report on the rectangles of the decomposition in which $p$ 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 $O(\lg r+t)$.

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 $\Delta$ be the maximum value s.t. for every cell of edge length $\geq\Delta$ at least one termination criterion holds. Then, the maximum tree height is $L=O(\lg(1/\Delta))$. Let $\beta$ be the minimum distance of two Voronoi vertices, and $\gamma$ the relative thickness of $\mathcal{P}^{c}$, i.e. the minimum diameter of a maximally inscribed $L_{\infty}$-ball in $\mathcal{P}^{c}$, where $\mathcal{P}^{c}$ is the complement of $\mathcal{P}$.

Lemma 10.

Separation bound $\Delta$ is $\Omega(\min\{\gamma,\beta\})$, where the asymptotic notation is used to hide some constants.

Proof.

The algorithm mainly subdivides cells that intersect $\mathcal{V}_{\mathcal{D}}(\mathcal{P})$, since a cell inside a Voronoi region or outside $\mathcal{P}$ is not subdivided (Termination criteria (T1), (T2)). Most subdivisions occur as long as non neighboring Voronoi regions are “too close”. Consider $\mathcal{C}$ centered at $p_{{}_{\mathcal{C}}}\in V_{\mathcal{D}}(s)$ and $s^{\prime}\in\widetilde{\phi}(\mathcal{C})\setminus\phi(\mathcal{C})$, with $V_{\mathcal{D}}(s),V_{\mathcal{D}}(s^{\prime})$ non neighboring. For $r_{{}_{\mathcal{C}}}<\frac{\mu_{\infty}(p_{{}_{\mathcal{C}}},s^{\prime})-\mu_{\infty}(p_{{}_{\mathcal{C}}},s)}{2}$ site $s^{\prime}$ is not in $\widetilde{\phi}(\mathcal{C})$. It holds that $\mu_{\infty}(p_{{}_{\mathcal{C}}},s^{\prime})-\mu_{\infty}(p_{{}_{\mathcal{C}}},s)\geq\zeta(s,s^{\prime})$, where $\zeta(s,s^{\prime})=\min\{\mu_{\infty}(p,q)\mid p\in\overline{V_{\mathcal{D}}(s)},q\in\overline{V_{\mathcal{D}}(s^{\prime})}\}$, i.e. the minimum distance of the closure of the two Voronoi regions. When $\overline{V_{\mathcal{D}}(s)},\overline{V_{\mathcal{D}}(s^{\prime})}$ are connected with a Voronoi edge, $\zeta(s,s^{\prime})=\Omega(\beta)$. When a minimum cell size of $\Omega(\beta)$ is not sufficient for $s^{\prime}$ to not belong in $\widetilde{\phi}(\mathcal{C})$, then there is a hole between $\overline{V_{\mathcal{D}}(s)},\overline{V_{\mathcal{D}}(s^{\prime})}$ and $\Delta$ is $\Omega(\gamma)$ in this case. ∎

This lower bound is tight: in Figure 15(a) for $\Delta=0.8125\beta$, and in Figure 15(b) for $\Delta=\gamma$. Next we target a realistic complexity analysis rather than worst-case. For this, assume the site distribution in $\mathcal{C}_{0}$ is “sufficiently uniform”. Formally:

Uniform Distribution Hypothesis (UDH): For $L_{\infty}$ balls $A_{1}\subseteq A_{0}\subset\mathcal{C}_{0}$, let $N_{0}$ (resp. $N_{1}$) be the number of sites intersecting $A_{0}$ (resp. $A_{1}$). We suppose ${N_{1}}/{N_{0}}=O({vol(A_{1})}/{vol(A_{0})})$, where $vol(\cdot)$ denotes the volume of a set in $\mathbb{R}^{d}$, $d$ being the dimension of $\mathcal{C}_{0}$.

Theorem 3.4.

Under UDH the algorithm complexity is $O(n/{\Delta}+1/\Delta^{2})$, where $n$ 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 $|\widetilde{\phi}(\mathcal{C}_{0})|=n$. The cardinality of active sets decreases as we move to the lower levels of the quadtree: Let $A(p,d,R)=\{q\in\mathbb{R}^{2}\mid d\leq\mu_{\infty}(p,q)\leq 2R+d\}$. For cell $\mathcal{C}$ and $s\in\widetilde{\phi}(\mathcal{C})$, $s\cap A(p_{{}_{\mathcal{C}}},\delta_{{}_{\mathcal{C}}},r_{{}_{\mathcal{C}}})\neq\emptyset$. Let $E=vol(A(p_{{}_{\mathcal{C}}},\delta_{{}_{\mathcal{C}}},r_{{}_{\mathcal{C}}}))$. For $\mathcal{C}_{1}$ a child of $\mathcal{C}$ and $s_{1}\in\widetilde{\phi}(\mathcal{C}_{1})$, $s_{1}\cap A(p_{{}_{\mathcal{C}_{1}}},\delta_{{}_{\mathcal{C}_{1}}},r_{{}_{\mathcal{C}_{1}}})\neq\emptyset$. Since $B_{\infty}(p_{{}_{\mathcal{C}}},\delta_{{}_{\mathcal{C}}})$ is empty of sites and may intersect with $A(p_{{}_{\mathcal{C}_{1}}},\delta_{{}_{\mathcal{C}_{1}}},r_{{}_{\mathcal{C}_{1}}})$, we let $E_{1}=vol(A(p_{{}_{\mathcal{C}_{1}}},\delta_{{}_{\mathcal{C}_{1}}},r_{{}_{\mathcal{C}_{1}}})\setminus(A(p_{{}_{\mathcal{C}_{1}}},\delta_{{}_{\mathcal{C}_{1}}},r_{{}_{\mathcal{C}_{1}}})\cap B_{\infty}(p_{{}_{\mathcal{C}}},\delta_{{}_{\mathcal{C}}})))$. We prove that in any combination of $\delta_{{}_{\mathcal{C}}},\delta_{{}_{\mathcal{C}_{1}}},r_{{}_{\mathcal{C}}}$ it is $E_{1}\leq E/2$. Under Hypothesis 1, a cell at tree level $i$ has $|\widetilde{\phi}(\mathcal{C}_{i})|=O(n/2^{i})$. Computation per tree level, is linear in sum of active sets’ cardinality, therefore summing over all levels of the tree, we find the that complexity of the subdivision phase is $O(n/\Delta)$. The complexity of the reconstruction phase is $O(\tilde{n})$, where $\tilde{n}$ is the number of leaf nodes in the quadtree, which is in turn $O(1/\Delta^{2})$. This allows to 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 $O(1)$, which is the typical case. Then, we prove the following.

Lemma 11.

We denote by $\mathcal{C}^{\prime}$ the parent of $\mathcal{C}$ in the subdivision. Using BVH accelerates the refinement of $\mathcal{C}$ if ${|\widetilde{\phi}(\mathcal{C}^{\prime})|}/{|\widetilde{\phi}(\mathcal{C})|}=\Omega(\lg n+t)$.

Proof.

A label set $\lambda(p)$ is determined by performing a point and a rectagle-intersection query; once the point is detected to lie inside a rectangle $R_{0}$ of a leaf $T_{0}$ we find an initial estimation $d_{0}$ of $\mu_{\infty}(p,\mathcal{P})$. Since the closest site to $p$ may be on another leaf, we do a rectangle-intersection query centered at $p$ with radius $d_{0}$. The closest site(s) to $p$ are in the intersected leaves. Thus, finding $\lambda(p)$ takes $O(k(p,d_{0})\lg r+k(p,d_{0})t)$ (Thm. 3.3), where $k(p,d_{0})=O(1)$ is the number of BVH leaves intersected by $\overline{B_{\infty}(p,d_{0})}$. Computing the sites in $\widetilde{\phi}(\mathcal{C})$ is accelerated if combined with a rectangle intersection query to find segments at $L_{\infty}$-distance $\leq 2r_{{}_{\mathcal{C}}}+\delta_{{}_{\mathcal{C}}}$ from $p_{{}_{\mathcal{C}}}$. Let $k_{{}_{\mathcal{C}}}$ be the maximum number of BVH leaves intersected by this rectangle-intersection query. We obtain a total refinement time for the cell equal to $O(k_{{}_{\mathcal{C}}}\lg r+k_{{}_{\mathcal{C}}}t)$. Since $k_{{}_{\mathcal{C}}}=O(|\widetilde{\phi}(\mathcal{C})|)$ and $r=O(n)$ the lemma follows. ∎

4 Subdivision algorithm in three dimensions

Let $\mathcal{P}$ be a manifold orthogonal polyhedron: every edge of $\mathcal{P}$ is shared by exactly two and every vertex by exactly 3 facets. For non-manifold input we first employ trihedralization of vertices, discussed in [16, 4]. Input consists of $\mathcal{S}$ and bounding box $\mathcal{C}_{0}$ of $\mathcal{P}$. 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 $s$, $\mathcal{Z}^{+}(s)$ is not necessarily convex and therefore the distance function $D_{s}(\cdot)$ cannot be computed in $O(1)$ time: it is not trivial to check membership in $\mathcal{Z}^{+}(s)$. 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 Chapter 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 (Theorem 3.3).

4.0.1 Subdivision.

The active set $\widetilde{\phi}$, $\phi$ and the label set of a point are defined as in to 2D. Most importantly, Lemma 4 is valid in 3D as well. The algorithm proceeds as follows: We recursively subdivide $\mathcal{C}_{0}$ 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 $\widetilde{\phi}$. Upon subdivision, we propagate $\widetilde{\phi}$ from a parent cell to its children for further refinement. We denote by $M$ the maximum degree of a Voronoi vertex $(M\leq 24)$.

3D Termination criteria: (T1’) $int(\mathcal{C})\cap\mathcal{P}=\emptyset$, (T2’) $|\widetilde{\phi}(\mathcal{C})|\leq 4$, (T3’) $|\widetilde{\phi}(\mathcal{C})|\leq M$ and the sites in $\widetilde{\phi}(\mathcal{C})$ define a unique Voronoi vertex $v\in\mathcal{C}$.

Subdivision is summarized in Alg. 2. (T1’) is valid for $\mathcal{C}$ iff $\lambda(p_{{}_{\mathcal{C}}})=\emptyset$ and $\forall s\in\widetilde{\phi}(C)$ it holds that $s\cap int(\mathcal{C})=\emptyset$. Detecting a Voronoi vertex in $\mathcal{C}$ 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.

(T1) used in 2D is omitted, for it is not efficiently decided: labels of the cell vertices cannot guarantee that $\mathcal{C}\subseteq V_{\mathcal{D}}(s)$. However, as the following lemma indicates, termination of the subdivision is not affected: cells contained in a Voronoi region whose radius is $<r^{*}$, where $r^{*}$ is a positive constant, are not subdivided.

Lemma 12 ().

Let $\mathcal{C}\subseteq V_{\mathcal{D}}(s)$. There exists $r^{*}>0$ s.t. if $r_{{}_{\mathcal{C}}}<r^{*}$ it holds that $\widetilde{\phi}(\mathcal{C})=\{s\}.$

Proof.

Let $s^{\prime}\in\mathcal{S}\setminus s$. We will prove that if $\mathcal{Z}^{+}(s^{\prime})\cap\mathcal{C}\neq\emptyset$, it holds that $\delta_{{}_{\mathcal{C}}}<\mu_{\infty}(p_{{}_{\mathcal{C}}},s^{\prime})$. Therefore there is $r(s^{\prime})>0$ such that $2r(s^{\prime})+\delta_{{}_{\mathcal{C}}}<\mu_{\infty}(p_{{}_{\mathcal{C}}},s^{\prime})$. Let $r^{*}$ be the minimum of these radii for every site different than $s$. When $r_{{}_{\mathcal{C}}}<r^{*}$, it holds that $\widetilde{\phi}(\mathcal{C})=\{s\}.$

Suppose that $\delta_{{}_{\mathcal{C}}}=\mu_{\infty}(p_{{}_{\mathcal{C}}},s^{\prime})=\mu_{\infty}(p_{{}_{\mathcal{C}}},q)$ for $q\in s^{\prime}$. Then $q\in\text{aff}(s)$ and $s^{\prime}$ cannot be a subset of $\text{aff}(s)$. So, $s^{\prime}$ is adjacent to $s$. If $\mathcal{Z}^{+}(s^{\prime})\cap\mathcal{C}=\emptyset$ then $s^{\prime}\not\in\widetilde{\phi}(\mathcal{C})$. If $\mathcal{Z}^{+}(s^{\prime})\cap\mathcal{C}\neq\emptyset$, since for adjacent sites it holds that $\mathcal{Z}^{+}(s)\cap\mathcal{Z}^{+}(s^{\prime})=\text{bis}_{\mathcal{D}}(s,s^{\prime})$, bisector intersects $\mathcal{C}$ which is a contradiction. Thus, there is no site $s^{\prime}$ with $\mathcal{Z}^{+}(s^{\prime})\cap\mathcal{C}\neq\emptyset$ and $\mu_{\infty}(p_{{}_{\mathcal{C}}},s^{\prime})=\delta_{{}_{\mathcal{C}}}$. ∎

Theorem 4.1 ().

Algorithm 2 halts.

4.0.2 Reconstruction.

We construct a graph $G=(V,E)$, representing the 1-skeleton of the Voronoi diagram. The nodes of $G$ are of two types, skeleton nodes and Voronoi vertex nodes, and are labeled by their closest sites. Skeleton nodes span Voronoi edges and are labeled by 3 or 4 sites. We visit the leaves of the octree and process cells with $|\widetilde{\phi}(\mathcal{C})|\geq 3$ and that do not satisfy (T1’). We introduce the nodes to the graph as in 2D. Graph edges are added between corners and Voronoi vertex nodes inside a cell. We run dual marching cubes (linear in the octree size) and connect graph nodes $v_{1},v_{2}$ located in neighboring cells, iff:

   $-$ $v_{1},v_{2}$ are skeleton nodes and $\lambda(v_{1})=\lambda(v_{2})$, or

   $-$ $v_{1}$ is a skeleton node, $v_{2}$ is a Voronoi vertex node and $\lambda(v_{1})\subset\lambda(v_{2})$, or

   $-$ $v_{1},v_{2}$ are Voronoi vertex nodes, $\lambda(v_{1})\cap\lambda(v_{2})=\{s,s^{\prime},s^{\prime\prime}\}$ and $v_{1}v_{2}\subset\mathcal{P}$.

Theorem 4.2 (Correctness).

The output graph is isomorphic to the 1-skeleton of the Voronoi diagram.

Primitives. Deciding membership in $\mathcal{H}(\cdot)$ is trivial. The predicates of Sect. 3.3 extend to 3D and the runtime of each is that of a rectangle-intersection query on the BVH constructed for the corresponding site at preprocessing: Let $pr_{\text{aff}(s)}(p)$ be the projection of $p$ to aff$(s)$ and $B_{p,s}$ the $2d-$box on $\text{aff}(s)$ centered at $pr_{\text{aff}(s)}(p)$ with radius $\mu_{\infty}(p,\text{aff}(s))$. Then, $p\in\mathcal{Z}(s)$ iff $B_{p,s}\cap s\neq\emptyset$. A query with $B_{p,s}$ is done by inZone($p,s$). For ZoneInCell($p,\mathcal{C}$) we do a query with $\overline{B_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},r_{{}_{\mathcal{C}}}+\rho_{i})}$ where $p_{{}_{\mathcal{C}}}^{\prime}=pr_{\text{aff}(s)}(p_{{}_{\mathcal{C}}})$, $\rho_{i}=\mu_{\infty}(f_{i},\text{aff}(s))$ and $f_{1},f_{2}$ the two facets of $\mathcal{C}$ parallel to aff$(s)$. Queries with $\overline{B_{\infty}(p_{{}_{\mathcal{C}}}^{\prime},r_{{}_{\mathcal{C}}})}$ are also performed by isIntersecting($s,\mathcal{C}$) and isStrictlyIntersecting($s,\mathcal{C}$). When computing label sets, LocationTest is slightly modified, since the corners used to identify the position of a point can also be formed by 3 sites.

Complexity. Under appropriate scaling so that the edge length of $\mathcal{C}_{0}$ be 1, if $\Delta$ is the separation bound, then the maximum height of the octree is $L=O(\lg(1/\Delta))$. The algorithm mainly subdivides cells intersecting $\mathcal{P}$, unlike the 2D algorithm that mainly subdivides cells intersecting $\mathcal{V}_{\mathcal{D}}(\mathcal{P})$, because a criterion like (T1) is missing. This absence does not affect tree height, since by Lemma 12 the minimum cell size is same as when we separate sites whose regions are non-neighboring (handled by Lemma 10). If $\beta$ is the minimum distance of two Voronoi vertices and $\gamma$ the relative thickness of $\mathcal{P}^{c}$, taking $\Delta=\Omega(\min\{\beta,\gamma\})$ suffices, as in 2D (Lemma 10).

Theorem 4.3.

Under UDH and if $n$ is the number of polyhedral facets, $\alpha$ the maximum number of vertices per facet and $t_{\alpha}$ the maximum number of rectangles in a BVH leaf, the algorithm’s complexity is $O(n\,\alpha(\lg\alpha+t_{\alpha})/\Delta^{2}+1/\Delta^{3})$.

Proof.

We sum active sets’ cardinalities of octree nodes, since refining a cell requires a number of rectangle-intersection queries linear in the size of its parent’s active set. Let $\mathcal{C}$ and its child $\mathcal{C}_{1}$. Any $s\in\widetilde{\phi}(\mathcal{C})$ satisfies $\delta_{{}_{\mathcal{C}}}\leq\mu_{\infty}(p_{{}_{\mathcal{C}}},s)\leq 2r_{{}_{\mathcal{C}}}+\delta_{{}_{\mathcal{C}}}$. We denote by $E$ the volume of the annulus $\{q\in\mathbb{R}^{2}\mid\delta_{{}_{\mathcal{C}}}\leq\mu_{\infty}(p_{{}_{\mathcal{C}}},q)\leq 2r_{{}_{\mathcal{C}}}+\delta_{{}_{\mathcal{C}}}\}$ and by $E_{1}$ the volume of the respective annulus for $\mathcal{C}_{1}$, minus the volume of the annulus’ intersection with $B_{\infty}(p_{{}_{\mathcal{C}}},\delta_{{}_{\mathcal{C}}})$. It is easy to show $E_{1}\leq E/2$. Under Hypothesis 1, we sum all levels and bound by $O(4^{L}n)$ the number of rectangle-intesection queries. Using Thm. 3.3, we find that the complexity of the subdivision phase is $O(n\,\alpha(\lg\alpha+t_{\alpha})/\Delta^{2})$. The complexity of the reconstruction phase is $O(\tilde{n})$, where $\tilde{n}$ is the number of leaf nodes in the octree, which is in turn $O(1/\Delta^{3})$. This allows to conclude. ∎

Remark 2.

$t_{\alpha}=O(\alpha)$ so the bound of Thm. 4.3 is $O(n\alpha^{2}/\Delta^{2}+1/\Delta^{3})$. Let $V$ be the number of input vertices. It is expected that $n\alpha=O(V)$; also $\alpha$ is usually constant. In this case, the complexity simplifies to $O(V/\Delta^{2}+1/\Delta^{3})$.

5 Implementation and concluding remarks

Our algorithms were implemented in Julia and use the algebraic geometric modeler Axl for visualization. They are available in https://gitlab.inria.fr/ckatsama/L_infinity_Voronoi. They are efficient in practice and their performance scales well in thousands of input sites. Some examples with runtimes are given in Figure 15 and 16. All experiments were run on a 64-bit machine with an Intel(R) Core(TM) i7-8550U CPU @1.80GHz and 8.00 GB of RAM.

Our complexity bounds rely on a hypothesis of ”sufficiently uniform” input. Even though Uniform Distribution Hypothesis might seem a strong assumption, note that our analysis does not encounter the width of the subdivision tree and considers it as a complete tree. We anticipate that in instances where our hypothesis does not hold, the “non-uniform” distribution of sites is captured by sudden increases or decreases in the width of the tree passing from a higher level to a lower one. We expect that there is a trade-off between the number of cells at level $i$ of the subdivision tree variating from $4^{i}$ and the cardinality of the active set of a cell being independent of the current level, so that the total complexity remains linear in the number of sites (facets). We plan to extend our research towards proving non trivial complexity bounds that exploit to the maximum degree the adaptivity of subdivision algorithms and do not rely on assumptions.

Acknowledgements. We thank Evanthia Papadopoulou for commenting on a preliminary version of the paper and Bernard Mourrain for collaborating on software. Both authors are members of AROMATH, a joint team between INRIA Sophia-Antipolis (France) and NKUA.