Randomized Incremental Construction of Delaunay Triangulations of Nice Point Sets

Randomized incremental construction (RIC) is one of the most important paradigms for building geometric data structures. Clarkson and Shor developed a general theory that led to numerous algorithms which are both simple and efficient in theory and in practice. Randomized incremental constructions are usually space-optimal and time-optimal in the worst case, as exemplified by the construction of convex hulls, Delaunay triangulations, and arrangements of line segments. However, the worst-case scenario occurs rarely in practice and we would like to understand how RIC behaves when the input is nice in the sense that the associated output is significantly smaller than in the worst case. For example, it is known that the Delaunay triangulation of nicely distributed points in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {E}}^d$$\end{document}Ed or on polyhedral surfaces in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {E}}^3$$\end{document}E3 has linear complexity, as opposed to a worst-case complexity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta (n^{\lfloor d/2\rfloor })$$\end{document}Θ(n⌊d/2⌋) in the first case and quadratic in the second. The standard analysis does not provide accurate bounds on the complexity of such cases and we aim at establishing such bounds in this paper. More precisely, we will show that, in the two cases above and variants of them, the complexity of the usual RIC is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(n\log n)$$\end{document}O(nlogn), which is optimal. In other words, without any modification, RIC nicely adapts to good cases of practical value. At the heart of our proof is a bound on the complexity of the Delaunay triangulation of random subsets of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varepsilon }$$\end{document}ε-nets. Along the way, we prove a probabilistic lemma for sampling without replacement, which may be of independent interest.


Introduction
The randomized incremental construction (RIC) is an algorithmic paradigm introduced by Clarkson and Shor [12], which has since found immense applicability in computational geometry, e.g. [27,28]. The general idea is to process the input points sequentially in a random order, and to analyze the expected complexity of the resulting procedure. The theory developed by Clarkson and Shor is quite general and led to numerous algorithms which are simple and efficient, both in theory and in practice. On the theory side, randomized incremental constructions are usually optimal in space and time in the worst case, as exemplified by the construction of convex hulls, Delaunay triangulations, and arrangements of line segments. Randomized incremental constructions are also known to perform very efficiently in practice, which, together with their simplicity, makes them the most popular candidates for implementations. Not surprisingly, the cgal library includes several randomized incremental algorithms, e.g. for computing Delaunay triangulations.
Experimental evidence has shown that randomized incremental constructions often work much better than the worst-case analysis suggests, which is fortunate since worst-case situations are rare in applications. This paper aims at understanding how randomized incremental constructions behave, when the input is nice in the sense that the associated construction is significantly smaller than in the worst case.
We need a model of good point sets to describe the input data and analyze the algorithms. This will be done through the notion of ε-nets, which have a long and rich history since their introduction in the 1950s in the works of Kolmogorov and others on functional analysis and topological vector spaces (see e.g. [31]). ε-Nets have become ubiquitous in many theoretical as well as applied areas, from geometry and functional analysis to probability theory and statistics, where they are often used as countable or finite approximations of continuous spaces.
When we work with such a hypothesis of "nice" distribution of the points in space, a volume counting argument ensures that the local complexity of the Delaunay triangulation around a vertex is bounded by a constant (dependent only on the dimension). Unfortunately, to be able to control the complexity of the usual randomized incremental algorithms [3,10,12,15], it is not enough to control the final complexity of the Delaunay triangulation. We need to control also the complexity of the triangulation of random subsets. One might expect that a random subsample of size k of an ε-net is also an ε -net for ε = ε d √ n/k. Actually this is not quite true, it may happen with reasonable probability that a ball of radius O(ε ) contains (log k/ log log k) points or that a ball of radius (ε d √ log k) does not contain any point. For the convenience of the reader, we briefly sketch the proofs in the appendix (Lemma A.1). However, it can only be shown that such a subsample is an (ε /(log (1/ε))-covering and an (ε log(1/ε))-packing, with high probability. Thus this approach can transfer the complexity of an ε-net to the one of a random subsample of an ε-net but with an extra multiplicative factor of (log 1/ε) = (log n). It follows that, in the two cases we consider, the standard analysis does not provide accurate bounds on the complexity of the (standard) randomized incremental construction. Our results are based on proving that, in expectation, the above bad scenarios occur rarely, and the algorithm achieves optimal time complexity.

Related work The Delaunay triangulations of nicely-distributed points have been
Extensions We also give some extensions of our results for periodic spaces. Our extensions are in four directions: (i) a more general notion of well-distributed point sets, the (ε, κ)-samples; (ii) a different notion of subsampling-the Bernoulli or i.i.d. sample where each point is selected to be in Y independently of the others, with probability q = s/n, (iii) a more general class of spaces-Euclidean d-orbifolds, and (iv) a more general class of metrics-those having bounded-distortion with respect to the Euclidean metric. Precisely, for all the above cases, we show that the Delaunay triangulation of a random subsample has a linear size in expectation. We believe that our methods should work for an even larger class of spaces, though this might require more delicate handling of boundary effects and other features specific to the metric space under consideration.
Outline The rest of the paper is as follows. In Sect. 2, we define the basic concepts of Delaunay triangulation, ε-net, flat torus, and random samples. We state our results in Sect. 3. In Sect. 4, we bound the size of the Delaunay triangulation of a uniform random sample of a given size extracted from an ε-net on the flat torus T d . In Sect. 5, we analyze the case when the uniform random subsample is drawn from an ε-net on a polyhedral surface in E 3 . In Sect. 6, we use the size bounds established in Sects. 4 and 5, to compute the space and time complexity of the randomized incremental construction for constructing Delaunay triangulations of ε-nets. Finally, in Sect. 7, we state and prove some extensions. Proofs missing in the main sections are given in the appendix.

Notations
We denote by ( p, r ), B( p, r ) and B[ p, r ], the sphere, the open ball, and the closed ball of center p and radius r respectively. For x ∈ E 2 , y ≥ 0, D(x, r ) denotes the disk with center x and radius r , i.e., the set of points {y ∈ E 2 : y − x < r }, and similarly D[x, r ] denotes the corresponding closed disk. The volume of the unit Euclidean ball of dimension d is denoted V d and the area of the boundary of such a ball is denoted S d−1 . It is known that For an event E in some probability space , we use 1 [E] to denote the indicator variable 1 [E] = 1 [E] (ω) which is 1 whenever ω ∈ E, and 0 otherwise. We use [n] to mean the set {1, 2, . . . , n}. Given a discrete set A, (A) denotes its cardinality and, for k ∈ Z + , A k denotes the set of k-sized subsets of A. Given an event A in some probability space, P[A] denotes the probability of A occurring. For a random variable Z in a probability space, E[Z ] denotes the expected value of Z . Lastly, e = 2.7182 . . . denotes the base of the natural logarithm.

"-Nets
A set X of n points in a metric space M, is an ε-packing if any pair of points in X are at least distance ε apart, and an ε-cover if each point in M is at distance at most ε from some point of X . The set X is an ε-net if it is an ε-cover and an ε-packing simultaneously. The definition of an ε-net applies to any metric space, but in the case of Euclidean metric, some additional properties can be proven. We shall use · to denote the Euclidean 2 norm. The following lemmas are folklore. Lemma 2.1 (Maximum packing size) Any packing of the ball of radius r ≥ ρ in dimension d by disjoint balls of radius ρ/2 has a number of balls smaller than (3r /ρ) d .
Proof Consider a maximal set of disjoint balls of radius ρ/2 with centers inside the ball B(r ) of radius r . Then the balls with the same centers and radius ρ cover the ball B(r ) (otherwise it contradicts the maximality). By a volume argument we get that the number of balls is bounded from above by Lemma 2.2 (Minimum cover size) Any covering of a ball of radius r in dimension d by balls of radius ρ has a number of balls greater than (r /ρ) d .

Proof
The volume argument gives a lower bound of For d ∈ Z + , the flat d-torus T d is the compact quotient group E d /Z d , with addition as the group action. More generally, for k ∈ Z + , the flat torus of length k is T d

Lemma 2.3 (ε-Net size bounds)
Given ε ∈ (0, 1/2], let X be an ε-net over the flat torus T d . Then, ( Proof Observe that, by the minimum distance property of the points in X , the balls of radius ε/2 centered at each point in X are disjoint, and by a volume argument there can be at most The balls of radius ε centered at each point in X cover the space thus their number is at least This completes the proof of the lemma.

Delaunay Triangulation
For simplicity of exposition and no real loss of generality, all finite point sets considered in this paper will be assumed to be in general position, i.e., no set of d + 2 points lie on a sphere. Given a set X in some ambient topological space, the Delaunay complex of X is the (abstract) simplicial complex with vertex set X which is the nerve of the Voronoi diagram of X , that is, a simplex σ (of arbitrary dimension) belongs to Del(X ) iff the Voronoi cells of its vertices have a nonempty common intersection. Equivalently, σ can be circumscribed by an empty ball, i.e., a ball whose bounding sphere contains the vertices of σ and whose interior contains no points of X .
The Delaunay complex is a triangulation if it triangulates the ambient space, or more precisely, the Delaunay complex Del(X ) of a point set X over an ambient space M, is said to be a Delaunay triangulation of M if there exists a homeomorphism between Del(X ) and M. The combinatorial complexity of a Delaunay triangulation is the total number of simplices of all dimensions, contained in the triangulation. The combinatorial complexity of a Delaunay triangulation is at most 2 d times the number of maximum-dimensional simplices contained in it.
Given a set X in some ambient space M, with its Delaunay complex Del(X ), the star of a subset S ∈ X , or star(S), is the set of all simplices in Del(X ) that are incident to at least one point in S. For a point p ∈ X , we shall use the shorthand expression star( p) to mean star({ p}).
Given topological spaces S and C, and a continuous map π : C → S, C is a covering space of S if π is such that for every point x ∈ S, there exists an open neighborhood U of x, such that the pre-image π −1 (U ) is a disjoint union of open neighborhoods in C, each of which is homeomorphically mapped onto U by π . A covering C of S is m-fold or m-sheeted if the cardinality of the pre-image of each point x ∈ S under the covering map is m. For example, T d k forms a k d -sheeted covering space of T d , with the covering map x → x mod 1, the modulus operation being defined coordinate-wise. Caroli and Teillaud [11] showed the following. Note that the above theorem implies that the Delaunay triangulation of any finite point set in T d always exists in the 3 d -sheeted covering of T d .
A key property of ε-nets is that their Delaunay triangulations have linear size.
Proof. Let us first bound the number of d-dimensional simplices in Del(X ). Observe that the circumradius of any d-simplex in Del(X ) cannot be greater than ε, since this would imply the existence of a ball in T d of radius at least ε, containing no points from X . Therefore given a point p ∈ X , any point that lies in a Delaunay simplex incident to p, must be at most of distance 2ε from p. Again by a volume argument, the number of such points is at most

Randomized Incremental Construction and Random Subsamples
For the algorithmic complexity aspects, we state a version of a standard theorem for the RIC procedure (see e.g. [14]). We first need a necessary condition for the theorem. When a new point p is added to an existing triangulation, a conflict is defined to be a previously existing simplex whose circumball contains p.  A subset Y of set X is a uniform random sample of X of size s if Y is any possible subset of X of size s with equal probabilities. In case the multiplicity of a point in X is greater than 1, the sample counts only one copy of the point; all other copies are present in Y if and only if the original point is present.
In order to work with uniform random samples, we shall prove a lemma on the uniformly random sampling distribution or sampling without replacement, which is stated below, and will be a key probabilistic component of our proofs. The lemma provides a bound on the probability of a non-monotone compound event, that is, if the event holds true for a fixed set of k points, there could exist supersets as well as subsets of the chosen set for which the event does not hold. This may well be of general interest, as most natural contiguity results with Bernoulli (i.e., independent) sampling are for monotone events. Lemma 2.8 Let a, b, c ∈ Z + with 2b ≤ a ≤ c and t ≤ c. Let C be a set, and B and T two disjoint subsets of C. If A is a random subset of C, chosen uniformly from all subsets of C having size a, the probability that A contains B and is disjoint from T , is at most where a, b, c are the cardinalities of A, B, and C respectively, and the cardinality of T is at least t.

Proof
The total number of ways of choosing the random sample A is c a . The number of ways of choosing A so that B ⊂ A and T ∩ A = ∅, is c−b−t a−b . Therefore the required probability is where in the last step observe that for the product for each i the term (1 − i/a) in the numerator is smaller than the corresponding term (1 − i/c) in the denominator, since a ≤ c. A similar observation holds for the product

Results
Random samples of ε-nets in T d The following theorem gives a constant bound on the expected complexity of star( p) for the Euclidean metric on the flat torus T d .
Polyhedral surfaces in E 3 A polyhedral surface S in E 3 is a collection of a finite number of polygons F ⊂ S, called facets, which are pairwise disjoint or meet at a vertex or along an edge. We show that the expected complexity of the Delaunay triangulation of a uniformly random subsample of an ε-net on a polyhedral surface is linear in the size of the subsample: , X be an ε-net on a fixed polyhedral surface S, with C facets having total area A and total length L along its boundaries, with n points and let Y ⊂ X be a random sub-sample of X having size s. Then the Delaunay triangulation Algorithmic bounds We next use the above combinatorial bounds to get the space and time complexity of the randomized incremental construction of the Delaunay triangulation of an ε-net on the flat d-torus or on a polyhedral surface in E 3 .

Euclidean Metric on T d
In this section, we prove that a subsample Y of a given size s, drawn randomly from an ε-net X ⊂ T d , has a Delaunay triangulation in which the star of any given vertex has a constant expected complexity. Hence, the expected complexity of the triangulation is linear in the size of the subsample. The constant of proportionality is bounded by 2 cd 2 , where c is a constant independent of ε and d.

Existence of Delaunay triangulation Del(Y):
In order to ensure we always have the Delaunay complex embedded in T d , we shall use Theorem 2.4. Accordingly, we get two different regimes of candidate simplices in the triangulation. When the circumradius of a candidate simplex σ ∈ star( p) is at most 1/4, then the simplex lies in a ball of radius at most 1/2 with center p. By Theorem 2.4, in this regime the Delaunay complex star( p) embeds in the one-sheeted covering of T d . Therefore, for a fixed set of vertices, there is a unique circumball. When the circumradius is greater than 1/4, the simplex is contained in a ball of radius greater than 1/2 around p, and therefore star( p) embeds in the 3 d -sheeted covering of T d , i.e., T d 3 . In this case, each vertex has 3 d copies, and so for a given choice of d vertices together with p, one can have (3 d ) d circumballs.
Proof framework Now we set up the formal proof. We first focus on bounding the expected number of d-simplices in star( p). Recall that n := (X ). Define q := (s − 1)/(n − 1) and δ : To bound the expected number of d-simplices, we shall consider the probability of presence in Del(Y), of d-simplices in X incident to p and having radius in the interval I k , as k ranges over Z + .
Throughout this section, we shall use σ to mean a d-simplex incident to p, with circumcenter c σ and circumradius r σ . To count the number of d-simplices in star( p) with circumradius in I k , let S p (k) denote the set of d-simplices having p as a vertex, with circumradius r σ ∈ I k . Set s p (k) := (S p (k)). For σ ∈ S p (k), let P p (k) denote the upper bound on the probability that σ appears in Del(Y), that is, Finally, let Z p (k) denote the number of d-simplices σ ∈ S p (k) such that σ ∈ Del(Y). The main lemma in the proof is a bound on the expected complexity of the star of p, in terms of s p (k) and P p (k).
Proof. For a simplex σ ∈ S p (k), let 1 [σ ] be the indicator random variable which is 1 if σ ∈ star( p), and zero otherwise. Then Z p (k) := σ ∈S p (k) 1 [σ ] , and (star( p)) = k≥0 Z p (k). Taking expectations over the random choice of Y, we get It only remains therefore, to establish bounds on s p (k) and P p (k) as functions of k, and finally to bound the sum k≥0 s p (k) P p (k). Following the earlier discussion on the existence of the Delaunay triangulation Del(Y), we shall split the sum k≥0 E[Z p (k)] into the two regimes, 0 ≤ k ≤ k max and k > k max , where k max denotes log 2 (1/4δ).
In this regime, the circumradii r σ of the candidate d-simplices are at most 1/4 since, by the definition of k max , we have r σ ≤ 2 k max = 1/4. Therefore every set of d + 1 vertices in Y, has a unique circumball. We begin by establishing a bound on P p (k). Let n k denote the minimum number of points of X in the interior of a circumball of a d-simplex σ ∈ S p (k), over all σ ∈ S p (k): n k := min σ ∈S p (k) (B(c σ , r σ ) ∩ X ). First, we bound n k from below using Lemma 2.2.

Lemma 4.2
Let σ be a simplex incident to p, having circumradius r σ ∈ I k , k ≥ 0. Then Proof When k ≤ k max , the radius of the circumball of a simplex σ is at most 2 k max δ ≤ 1/4 < 1/2. Applying Theorem 2.4, we work in the one-sheeted covering of T d . Using the fact that X is an ε-covering, we apply Lemma 2.2 to get that n k ≥ (2 k−1 δ/ε) d . Now applying Lemma 2.8, we can bound P p (k).
Proof The simplex σ can be a Delaunay simplex only if (i) the set of its vertices is included in the subsample Y, and (ii) all points in B(c σ , r σ ) ∩ X are excluded from Y. The idea is therefore, to use Lemma 4.2 to bound the number of points in B(c σ , r σ ) ∩ X from below by n k , and then upper bound the probability that all these points are excluded from Y.
This suggests applying Lemma 2.8, with the universe having c = n − 1 elements, sample size a = s − 1, included set having b = d elements, and excluded set having t = n k elements. We verify first that the conditions of the lemma are satisfied, i.e., where the equality was by the substitution q = (s − 1)/(n − 1), and the last inequality followed from the fact that Next, we shall upper bound s p (k) from above. We first state a simple observation. Proof This follows simply from the triangle inequality. For the first statement, we have that for any p ∈ σ , p, p ≤ p, c σ + c σ , p = 2r σ . The second statement follows by replacing the above inequalities with strict inequalities for the points in the open ball B(c σ , r σ ).
Now we can bound s p (k) using the above observation together with Lemma 2.1.
Next, using the above bounds on values s p (k) and n k , we shall bound on the sum in the following three lemmas.
Proof Substituting the bounds on s p (0) and P p (0) proven in Lemmas 4.2, 4.3, and 4.5, we have where in the second inequality we used the definition of δ to get q(δ/ε) d = 2d, and in the last inequality we used Stirling's approximation d d /d! ≤ e d , and the fact that 2e < 6.
. Now from Lemmas 4.2, 4.3, and 4.5, we have that for all k ≥ 0, Observe that k≥0 s p (k) P p (k) ≤ k≥0s p (k)P p (k). In order to bound k≥0 s p (k) P p (k), it therefore suffices to simply bound k≥0s p (k)P p (k). For the rest of the proof, therefore, we shall focus on bounding this sum.
From the analysis, it will be easy to see that the close neighbor vertices will form the bulk of the d-simplices in E [star( p)], and the contribution of vertices farther from p will decrease exponentially with the distance.
When k = 0 In this case, the ratiõ which is at most (2/e) 6 for d ≥ 2. Therefore for all 0 ≤ k ≤ k max − 1, we have that Thus, the sum k max k=0 E[Z p (k)] is upper bounded by the sum of a geometric progression with leading terms p (0)P p (0) ≤ 6 d 2 +d and common ratio (2/e) 6 , which is at most (1 − (2/e) 6 ) −1 · 6 d 2 +d . In other words, most of the Delaunay neighbors of p are the points nearest to p.
Lastly, we shall show that the expected number of d-simplices with exceptionally large circumradii, i.e., k>k max E[Z p (k)], is at most a small constant.
Case II: d-simplices with large circumradii k > k max . In this regime, the circumradii of the candidate d-simplices are greater than 1/4. Therefore, by Theorem 2.4, we shall work in the 3 d -sheeted covering of T d .
Proof From Lemma 2.3, we have that n = (X ) ≤ 2 −d d d/2 ε −d . Therefore, by Lemma 2.2, any ball B of radius at least 2 k max δ ≥ 1/4, has at least Here, since 2 k+1 δ > 1/2, we shall use Theorem 2.4 and work in the 3 d -sheeted covering of T d . The maximum number of d-tuples which can possibly form a Delaunay d-simplex with p, is at most n−1 d ≤ (n−1) d /d! Since we are working in the 3 d -sheeted covering space, each vertex of a simplex σ ∈ S p (k) can be chosen from one of at most 3 d copies in the covering space. Thus, each simplex in S p (k) yields less than 3 d 2 possible Delaunay spheres in T d . Therefore, the expected number of d-simplices having radius at least 2 k max δ, is at most (1) is a decreasing function of s, and it is easy to check that the value at s 0 is smaller than 5. The lemma follows.

Polyhedral Surfaces in E 3
In this section, we introduce a partition of the sub-sample Y into boundary and interior points, to do a case analysis of the expected number of edges in the Delaunay triangulation Del(Y), depending on whether the end-points of a potential Delaunay edge are boundary or interior points, and whether they lie on the same facet or on different facets.

Main ideas
Our overall strategy will be to mesh the proofs of Attali-Boissonnat [4] and Theorem 3.1. Briefly, Attali and Boissonnat reduce the problem to counting the Delaunay edges of the point sample, which, up to a constant factor, is equal to the number of simplices (see the proof of Theorem 3.2). They then count the edges by distinguishing between boundary and interior points of a facet. For boundary points, they allow all possible edges. For interior points, the case of edges with end-points on the same facet is easy to handle, while geometric constructions are required to handle the case of end-points on different facets, or that of edges with one end-point in the interior and another on the boundary.
However, we shall need to introduce some new ideas to adapt our previous methods to this setting. Firstly, for any pair of vertices, there are infinitely many balls containing the pair on the boundary, and as soon as one of these balls is empty (with respect to the points in Y), the pair appears as an edge in the Delaunay triangulation. This is handled using a geometric construction (see Lemma 5.9). Basically, the idea is to build a constant-sized packing of a sphere centered on a given point, using large balls, such that any sphere of a sufficiently large radius which passes through the point, must contain a ball from the packing.
Secondly, since we no longer have a regular point distribution on the boundaries of S, boundary effects could penetrate deep into the interior. To handle this, we introduce the notion of levels of a surface, instead of the fixed strip around the boundary used in [4], and use a probabilistic, rather than deterministic, classification of boundary and interior points. The new classification is based on the level of a point and the radius of the largest empty disk passing through it.
Recall the definitions of X , Y, and S from Theorem 3.2. We shall use κ to denote the maximum number of points of a given point set in a disk of radius 2ε. When X is an ε-net, κ is at most 6 d = 36 (using Lemma 2.1 with r = 2ε, d = 2, and ρ = ε). We define the sampling density q := s/n and the sampling radius δ := ε/ √ q. For a curve , l( ) denotes its length. For a subset of a surface R ⊂ S, a(R) denotes the area of R. For sets A, B ⊂ E 3 , A ⊕ B denotes the Minkowski sum of A and B, i.e., the set {x + y : x ∈ A, y ∈ B}. For convenience, the special case A ⊕ B(0, r ) shall be denoted by A ⊕ r .
We next present some general lemmas, which will be needed in the proofs of the main lemmas.

Level sets, boundary points and interior points
We now introduce some definitions which will play a central role in the analysis. First we define the notion of levels. Given facet F ∈ S and k ≥ 0, define the level set L ≤k := F ∩ (∂ F ⊕ 2 k δ), L =k := L ≤k \ L ≤k−1 . For x ∈ X , the level of x, denoted Lev(x), is k such that x ∈ L =k . Let L ≤k (X ), L =k (X ) denote L ≤k ∩X , L =k ∩X respectively. Note that for x ∈ L =k , k ≥ 1, Next, we define a bi-partition of the point set into boundary and interior points. Given a point x ∈ F, x is a boundary point or x ∈ Bd F (Y), if (i) d(x, ∂ F) ≤ δ, or (ii) there exists an empty disk (w.r.t. Y) of radius greater than 2 k−1 δ, whose boundary passes through x, where k is the minimum positive integer such that d( The above bi-partition induces a classification of potential Delaunay edges, depending on whether the end-points are boundary or interior points. Let E 1 denote the set of edges {x, y} in Del(Y) such that x, y ∈ Bd S (Y). Let E 2 denote the set of edges {x, y} in Del(Y) such that x, y ∈ Int F (Y) for some F ∈ S. Let E 3 denote the set of edges {x, y} in Del(Y) such that x, y ∈ Int S (Y) and x ∈ F, y ∈ F , F, F different facets in S. Let E 4 denote the set of edges Recall that the polyhedral surface S has C facets, with total area of its faces A, and total length of its boundaries L. We have the following lemmas, to be proven in Sect. 5.2.

Lemma 5.4 E [ (E 4 )] ≤ O(1) · κ 2 L 2 s/A.
Given the above lemmas, the proof of Theorem 3.2 follows easily. [4,Sect. 4], by Euler's formula, the number of tetrahedra t(Del(Y)) in the Delaunay triangulation of S, is at most e(Del(Y)) − (Y) = e(Del(Y)) − s, where e(Del(Y)) is the number of edges in the Delaunay triangulation. Observe that the number of triangles in Del(Y) is also at most 4 times the number of tetrahedra. Therefore to bound the complexity of Del(Y) up to constant factors, it suffices to count the edges of Del(Y). Next, observe that any point x ∈ Y is either a boundary or an interior point, that is, Bd S (Y) Int S (Y) = Y. An edge in Del(Y), therefore, can be either between two points in Bd S (Y), or two points in Int S (Y), or between a point in Bd S (Y) and another in Int S (Y). The case of a pair of points in Int S (Y) is further split based on whether the points belong to the same facet of S or different facets. Thus using the above exhaustive case analysis, the proof follows simply by summing the bounds.

Proof of Theorem 3.2 As in
Before proving Lemmas 5.1-5.4, we first present

Some Technical Results
The following geometric and probabilistic lemmas prove certain properties of ε-nets on polyhedral surfaces and random subsets, as well as exploit the notion of boundary and interior points to get an exponential decay for boundary effects penetrating into the interior.

Proposition 5.5 [4] Let F be a facet of S. For any Borel set R ⊂ F, we have a(R)
and therefore, Proposition 5.6 [4] Let F be a facet of S, a curve contained in F, and j ∈ N. Then The following fact can be easily verified using a greedy strategy. The next lemma bounds the number of points of the original net X , in a given level of S (defined with respect to the sampling radius δ), in terms of the total length L of the facet boundaries, together with the sampling radius δ and some parameters related to the net X .
Proof. The first inequality is obvious, as L =k ⊆ L ≤k . The proof of the second inequality follows by applying Proposition 5.6 over the boundaries of the facets. For a fixed facet F ∈ S, we get Summing over all F ∈ S, we get As mentioned previously, our goal is to bound the expected number of edges appearing in Del(Y) by bounding the probability that a given pair of vertices in Y appears as an edge. However, for a given pair of vertices there are infinitely many balls containing this pair on the boundary, and if any of these balls contains a point of Y in their interior, then the pair must form an edge in Del(Y). Bounding the above probability therefore seems to require us to forbid infinitely many events from occurring. The next couple of lemmas show how to restrict this to forbidding finitely many events, and obtain an exponentially decreasing bound on the probability of having long Delaunay edges.
Given a point x on some facet F ∈ S with supporting plane P, we call a disk D ⊂ P a supporting disk for x, if the boundary ∂ D x and Y ∩ int D = ∅, that is, the interior of the disk D does not contain any point of the random sample Y. The following lemma shows that given any point x on a facet F ∈ S with supporting plane P, there exists a finite collection of sufficiently large discs in F, such that any large supporting disc for x must contain a disc from K x .

Lemma 5.9
Let F be a facet of S with supporting plane P, and x ∈ F with Lev(x) > 0. Then there exists a collection K x of at most c B = 7 disks in F such that (i) each D ∈ K x is contained in F; (ii) each D ∈ K x has radius r 0 /4, where r 0 = 2 k δ and k ∈ N is such that 0 ≤ k < Lev(x); and (iii) any supporting disk for x of radius at least r 0 , contains at least one disk in K x .
Proof Let D 0 ⊂ P be a supporting disk for x, with center y ∈ P and radius r ≥ r 0 . Let D 1 = D(x , r 0 ) be the unique disk with center x on the line x y, radius r 0 , and having x ∈ ∂ D 1 . Note that Consider 2 = (x, r 0 /2), and let p = x x ∩ 2 , that is, the point p lies on the line x x , at distance r 0 /2 from x (and therefore from x as well). We shall build a minimal covering K of the circle 2 , by disks centered in 2 , having radius r 0 /4 (see Fig. 2 for the construction). From Lemma 5.7, we get (K) = 7. Let D ∈ K be a disk in the covering. Then by the triangle inequality, D ⊂ D (x, r 0 /2 + r 0 /4) ⊂ D(x, r 0 ). As before, by the definitions of Lev(x) and r 0 , this implies D ⊂ F. Thus K satisfies conditions (i) and (ii) of the lemma. Further, since K is a covering of 2 , there exists D p ∈ K such that p ∈ D p . Therefore, the disk D p ⊂ D 1 ⊂ D 0 , and D p ⊂ F. Thus D p ∈ K satisfies condition (iii). Now taking K x = K completes the proof of the lemma.
The following lemma gives a bound on the probability of having a large supporting disk, which falls exponentially as the radius of the disk increases. This will be crucial to our main proof, as it will be used in the proofs of all Lemmas 5. 1-5.4. The idea is to use the construction of the previous lemma, together with the probabilistic bounds from Lemma 2.8.

Lemma 5.10 (Decay lemma)
Note that we bound directly the joint probability of a collection of t vertices, simultaneously being chosen in Y and each of them having a supporting disk. This slight technical complication is required due to the fact that all these simultaneously occurring events are mutually dependent, as the random sample Y is chosen with a fixed size rather than choosing points independently. This is where the probabilistic Lemma 2.8 in handling such compound events proves to be so versatile.

Proof of Lemma 5.10
Firstly, consider the case where k max = 0, i.e., all the k i 's are zero. In this case we simply upper bound the probability of the event E, by the probability of including all the points x 1 , . . . , x t in Y. By Lemma 2.8, this is at most q t .
We now come to the case when k max > 0. Since for all i ∈ [t], k i < Lev(x i ), we can apply Lemma 5.9 for each i,with k = k i , to conclude that for each i, there exists a collection K i of at most c B disks of radius r * i /4, such that any supporting disk for x i having radius greater than r * i , must contain some disk D = (B 1 , . . . , B t ) ∈ T denote a generic element of T . Taking the union bound over the set T gives the following: Let j := arg max i∈[t] k i , so that k max = k j . Now, the event E requires the set x 1 , . . . x t to be in the sample Y, and the interiors of the disks D * i to be free from points in Y. In particular, the disk D * j should not contain any points in Y. Therefore applying Lemma 2.8 for the universe C = X , the random sample A = Y, the included subset B = {x 1 , . . . , x t }, and the excluded subset Z = D * j ∩ X of size at least z = π(r * j /4) 2 /(4πε 2 ), we get where in the last step we used the fact that q = s/n, δ = ε √ n/s, r * j = 2 k max δ, and set c 1 = c t B . This completes the proof of the lemma. The next lemma gives upper and lower bounds on the number of points in the original set X , as well as the point sample Y, contained in a disc in some facet of S, as a function of the disc radius.
Lemma 5.11 (Growth lemma) Given any point x ∈ S in a facet F, and 0 ≤ k < Lev(x), we have Proof By the definition of Lev(x), we have that D(x, 2 Lev(x)−1 δ) ⊂ F. Now the statement follows by the application of Proposition 5.5: where δ := ε √ n/s. This gives the first statement of the lemma, with q = s/n. The second statement follows simply by taking expectation.

Proofs of Lemmas 5.1-5.4
The proofs of Lemmas 5.1 and 5.2 now follow from the decay and growth lemmas, together with similar ideas as for the flat torus case.
Proof of Lemma 5.1 Let x 1 , x 2 ∈ Bd S (Y). To bound the expected number of edges in E 1 , we simply bound the number of pairs (x 1 , x 2 ) ∈ Bd S (Y) · Bd S (Y). Let l 1 := Lev(x) and l 2 := Lev(y), and let l : . By definition, if l = 0, then x 1 , x 2 ∈ Bd S (Y). For l ≥ 1, we get that x 1 ∈ Bd S (Y) and x 2 ∈ Bd S (Y) only if there exists a disk of radius at least 2 l−1 δ passing through x 1 or x 2 , and containing no points of Y. Therefore to bound the probability that (x 1 , x 2 ) ∈ (Bd S (Y)) 2 , we can apply the decay Lemma 5.10, with t = 2, for i ∈ {1, 2}. We get where c 2 = c 2 /4 = 2 −9 . Summing over all choices of levels of x 1 and x 2 , we have By symmetry, it is enough to assume without loss of generality that l 1 ≥ l 2 , i.e., l = l 1 . Thus, Applying (6) and the level size Lemma 5.8, we get Using the definitions of q and δ, together with Proposition 5.5, and writing the terms outside the summation as N 1 , we obtain The summation can be bounded using Lemma A.2 to get Now substituting c 2 = 2 −9 gives E[ (E 1 )] ≤ 2 · 10 4 c 1 s · 4π(9κ L) 2 /A.

Proof of Lemma 5.2
Let l stand for min {Lev(x), Lev(y)}. Observe that if l = 0, then either x or y is a boundary point, and hence we can assume l ≥ 1. Let x = arg min z∈∂ F d(z, x), and y = arg min z∈∂ F d(z, y), i.e., x is the closest point to x in ∂ F, and similarly for y . By the definition of Bd S (Y), observe that By the growth Lemma 5.11, the expected number of Delaunay neighbors y of a point x such that d(x, y) ≤ δ is at most E[D(x, δ) ∩ Y] ≤ q · (4/q) = 4. Thus the expected number of edges in E 2 from pairs (x, y) with x, y ∈ Int F (Y) for some F ∈ S, and d(x, y) ≤ δ, is at most 4 · (Y). For longer-distance edges, let k ≥ 1 be such that 2 k−1 δ ≤ d(x, y) ≤ 2 k δ. Taking t = 2, x 1 = x, x 2 = y, k 1 = k − 1, and k 2 ≤ k 1 , and applying the decay Lemma 5.10, we get that where c 2 = c 2 /4. Summing over all possible choices of l ≥ 1 and k ≤ l, we get where in the second step we applied the growth Lemma 5.11, and in the third step we bounded the sum k≥0 2 2k exp (−c 2 2 2k ) using Lemma A.2 and the fact that q = s/n = ε 2 /δ 2 . Note that c 3 ≤ 4 · 10 3 , and c 4 := c 1 c 3 ≤ 2 · 10 5 .
For the proofs of Lemmas 5.3 and 5.4, we need some more geometric ideas of [4].

Proof of Lemma 5.3 Let x, y ∈ Int S (Y)
, where x ∈ F and y ∈ F , for some F, F ∈ S. Let F be fixed. To analyze this case, we shall first give a geometric construction of [4], and state an observation from its proof. ) Let P and P denote the supporting planes of the facets F and F , respectively. Let P B be the bisector plane of P and P . We denote by x ∈ P the reflection of x ∈ P with respect to P B , and similarly, by y ∈ P the reflection of y ∈ P . Let B = B xyx y be the smallest ball in E 3 passing through x, y, x , y , having intersections D 1 = B ∩ P and D 2 = B ∩ P with P and P , respectively.

Construction 5.12
Attali and Boissonnat [4] observed that Therefore, if there exists a ball B ∈ E 3 such that x, y ∈ ∂ B, and Y ∩ int B = ∅, Observe that, as in Case II, we have x ∈ D(y, 2 Lev(y) δ), since otherwise y ∈ Bd S (Y). (Note that our definition of boundary points allows us to ignore the fact that x is not necessarily a point in X .) Further, the set {x ∈ D(y, 2 k δ)}, 0 < k < Lev(y), is bounded in size by (D(y , 2 k δ) ∩ F ∩ Y). The rest of the analysis for the fixed facet F , therefore, follows as in Case II. Summing over all F ∈ S \ {F}, we get E[ (E 3 )] ≤ c 4 (C − 1)κs.
Before proving Lemma 5.4, we briefly describe a construction which will be central to our analysis. ) Let P be a plane and Z be a finite set of points. To each point x ∈ Z , assign the region V (x) = V x (Z ) ⊂ P of points y ∈ P such that the sphere tangent to P at y and passing through x encloses no point of Z .

Construction 5.14
We summarize some conclusions of Attali-Boissonnat regarding the construction. Proofs can be found in [4].

Proposition 5.15 (i) V is a partition of P.
(ii) For each x ∈ Z, V (x) is an intersection of regions that are either disks or complements of disks. (iii) The total length of the boundary curves in V is equal to the total length of the convex boundaries.
Proof The proofs of (i) and (ii) are easy. We prove (iii). Consider a point x ∈ Z , and let V (x) be the region corresponding to x in V. By (ii), V (x) = D∈D x D ∩ D ∈C xD , where D x is a set of disks and C x is a set of complements of disks in the plane P. Let y ∈ ∂ V (x). Then if y ∈ D∈D x D, there exists D 1 ∈ D x such that y ∈ ∂ D 1 , and so y is part of a convex segment in ∂ V (x). Otherwise, there existsD 2 ∈ C x such that y ∈ ∂D 2 . In this case, let V (z), z ∈ Z , denote the region such that y ∈ ∂ V (z). Then D 2 ⊃ V (z), and therefore y belongs to a convex segment in ∂ V (z). Thus, every point y ∈ ∂ V (x) is convex either for V (x) or for a neighboring region of V (x), and so the total length of the convex boundary curves in V gives the total length of all the boundary curves.
For the rest of this subsection, we shall apply Construction 5.14 for the plane P and the points in Bd S (Y) as Z . Let T := Int F (Y) for some facet F ∈ S. Given x ∈ Z , y ∈ P \ V (x), let k y = k y (x) denote the least k ≥ 0 such that y ∈ ∂ V (x) ⊕ 2 k δ. Also, x ∈ V (x), by the definition of V (x). Therefore we get x ∈ D y ∩ V (x).
Proof Suppose {x, y} ∈ E 4 . Then there exists a ball B ∈ E 3 with x, y ∈ ∂ B, and Y ∩ int B = ∅. Therefore D y := B ∩ P also satisfies Y ∩ int D y = ∅. By Proposition 5. 16 we have that D y ∩ V (x) = ∅. Therefore, y ∈ V (x) ⊕ 2r y , where r y is the radius of D y . But since y ∈ Int(F), we have that any disk having y on its boundary and containing no point of Y in its interior can have radius at most 2 Lev(y)−1 δ. Therefore r y ≤ 2 Lev(y)−1 δ. Now taking k y such that 2 k y δ = 2r y , we get k y ≤ Lev(y). Now we partition the pairs of vertices {x, y} ∈ E 4 with x ∈ Bd S (Y), depending on whether y ∈ V F (x) or y ∈ ∂ V F (x) ⊕ 2 k y δ. That is, given a facet F ∈ S, let E 4 (Int(F)) denote the set of edges {x, y} ∈ E 4 with y ∈ int V F (x), and E 4 (Bd(F)) denote the set of edges in E 4 with y ∈ ∂(V F (x)) ⊕ 2 k δ, for k ∈ [0, k y ]. Define E 4 (Int) := F∈S E 4 (Int(F)) and E 4 (Bd) := F∈S E 4 (Bd(F)), respectively.

Proof of Lemma 5.4
The proof follows from Lemmas 5.18 and 5.19, which bound the expected number of edges in E 4 (Int) and E 4 (Bd) respectively.
is a fixed set of points. The number of pairs contributing to E 4 Since the last bound holds for any choice of Y , taking expectation over all choices we get Proof To compute the expected value of E 4 (Bd(S)), fix a face F ∈ S. Consider a pair of points x, y ∈ X such that y ∈ F. Let E x,y denote the event {x, y} ∈ E 4 (Bd(F)). The value of E 4 (Bd) is the number of x, y ∈ X such that E x,y occurs. Taking expectations, Observe that E x,y occurs only if (i) x ∈ Bd S (Y) and (ii) k y (x) ≤ Lev(y), by applying Construction 5.14, for the plane P, Z = Bd S (Y), T = Y ∩ P, and using Proposition 5.16. By Lemma 5.17, k y (x) ∈ [0, Lev(y)]. Let P k 1 ,k 2 denote the probability that {x, y} ∈ E 4 (Bd(F)), with Lev(x) = k 1 and k y (x) = k 2 . Equation (7) can be rewritten in terms of k 1 and k 2 as Applying the decay Lemma 5.10 with t = 2, x 1 = x, x 2 = y, k 1 = max {0, k 1 − 1} (since x ∈ Bd S (Y)), and k 2 = max {0, k 2 − 1}, we get where k * := max {0, k 1 − 1, k 2 − 1}, and f (k * ) = 0 if k * = 0 and c 2 2 2k * otherwise, with c 2 = c 2 /4.
As in the proof of Lemma 5.1, we shall use symmetry to handle the cases k 1 ≥ k 2 and k 2 > k 1 together. We get By the level size Lemma 5.8, we get that (L =k 1 ∩ X ) ≤ 9κ L2 k 1 δ/ε 2 . Using Proposition 5.6, we get ({∂ V (x) ⊕ 2 k 2 δ} ∩ Y) ≤ 9κ 2 k 2 δ · l(∂ V (x))/ε 2 . By Proposition 5.15 (iii), each boundary in the partition V is convex for some x ∈ Bd S (Y). Therefore we need to sum l(∂ V (x)) only over the convex curves in ∂ V (x), x ∈ Bd S (Y). The length of these curves is at most l(∂ F). Thus we get Using Lemma A.2, the above summation is bounded by a constant. This results in

Randomized Incremental Construction (Proof of Theorem 3.3)
In this section, we show how Theorems 3.1 and 3.2 imply bounds on the computational complexity of constructing Delaunay triangulations of ε-nets. Our main tool shall be Theorem 2.7. However, we need to show first that Condition 2.6 holds. The standard proof of this (see e.g. [10,12], and also the discussion in [

Euclidean Orbifolds and Bounded-Distortion Metrics
In this section, we shall give some extensions of Theorems 3.1 and 3.3. The proofs of our theorems follow by finding covering spaces of bounded multiplicity where the Delaunay complex can be embedded, and generalizing Lemmas 4.5-4.8 to such spaces. Given a space S, ε ∈ [0, 1], and κ ∈ Z + , an (ε, κ)-sample is a set of points for which any ball of radius ε in S, contains at least one point and at most κ points. The proof of Theorem 3.2 goes through for (ε, κ)-samples as well, using a generic κ in place of the fixed value used in the proof. Similarly the proofs of Theorems 3.1 and 3.3 can also be translated to the (ε, κ) setting. We have Theorem 7.1 Theorems 3.1, 3.2, and 3.3 hold when the point set X is an (ε, κ)-sample.
Our combinatorial results for uniformly random subsets of nicely-distributed point sets also hold when the subsets are chosen by independent sampling. The only change in the proof of Theorem 3.1 for independent sampling is when computing P p (k). We make use of the fact that points are selected independently, to get directly that P p (k) ≤ q d (1 − q) n k ≤ q d exp (−qn k ). A similar adjustment is needed in the proof of Theorem 3.2. The rest of the proofs follow as before. Thus we get Theorem 7.2 Theorems 3.1 and 3.2 also hold in the case when the random sample is an i.i.d. sample with probability parameter q = s/n.
Coming to our results for Delaunay triangulations of Euclidean d-manifolds and embedded metrics in T d , we need a few definitions first.
Euclidean d-orbifolds A d-dimensional Bieberbach group G is a discrete group of isometries acting on E d . A d-orbifold E d /G is the compact quotient space (i.e., collection of orbits) of E d acted on by a d-dimensional Bieberbach group G. When the group action is free (i.e., has no fixed points), the d-orbifold is a closed Euclidean d-manifold. Every Euclidean d-manifold is the quotient space of some d-Bieberbach group acting on E d [7,32]. For Euclidean d-orbifolds we have: Proof In this case, the existence of the covering space follows from the algorithmic version of Bieberbach's theorem [6] by Sect. 4]. Now, observe that in the Euclidean metric, the diameter · (Z d P) of the largest ball not containing any point from the set Z d P is at most √ d, with equality holding when (P) = 1, and that (S ) ≤ (S) for any S ⊇ S, since adding points can only decrease the diameter of the largest empty ball. Therefore, in the metric d, we have that By the previously observed bound on · (Z d P), we have

Conclusion and Remarks
In this paper, we analyzed the behavior of the usual RIC algorithm for the Delaunay triangulation of nice point sets, focusing on the cases where the ambient space is the flat d-torus or a polyhedral surface in E 3 . Similar questions can be asked for other spaces where the Delaunay triangulation is known to have low complexity for "nice" point sets.
We leave for further research a more general analysis of RIC of Delaunay triangulations of cases such as polyhedral surfaces in higher dimensions, as well as extending the techniques developed in this paper to the RIC of other geometric problems.
"bins". Now choosing k points from the n points in the ε-net is akin to throwing k packets into k bins-with high probability, the maximum load of a bin is (log k/ log log k), i.e., there exists a ball of radius ε with (log k/ log log k) points.
(ii) Similarly, there are (ε / log k) = (k/ log k) balls of radius ε d √ log k that cover the ε-net, and choosing k points from this covering can be thought of as an instance of the coupon collector problem: we are allowed k draws (i.e., we choose k points), each draw is one of the "coupons" (i.e., the balls in the covering) with equal probability, and we want the probability that at least one coupon has not been collected, i.e., at least one ball from the covering has no point chosen from it. Now the well-known lower tail of the number of draws taken to collect all coupons in the coupon collector problem, gives that with high probability, for k = (k/ log k) coupons, and at most O(k log k ) = k draws, all coupons will not be collected, that is, some ball of radius (ε d √ log k) will be empty.
Proof The proof follows from elementary calculus. The maximum term is when 2 an = 1/b, i.e., n = log 2 (1/b)/a, and evaluates to 2 an exp (−b2 an ) = 1/eb, and terms decrease exponentially afterwards. The sum is therefore upper bounded by twice the maximum term, times the number of terms before the maximum, which is the claimed bound.