A probabilistic approach to reducing the algebraiccomplexity of computing Delaunay triangulations

A probabilistic approach to reducing the algebraic
complexity of computing Delaunay triangulations

Jean-Daniel Boissonnat 1    Ramsay Dyer 1    Arijit Ghosh 1
11 Geometrica, INRIA, Sophia Antipolis, France Jean-Daniel.Boissonnat@inria.fr
Abstract

Computing Delaunay triangulations in involves evaluating the so-called in_sphere predicate that determines if a point lies inside, on or outside the sphere circumscribing points . This predicate reduces to evaluating the sign of a multivariate polynomial of degree in the coordinates of the points . Despite much progress on exact geometric computing, the fact that the degree of the polynomial increases with makes the evaluation of the sign of such a polynomial problematic except in very low dimensions. In this paper, we propose a new approach that is based on the witness complex, a weak form of the Delaunay complex introduced by Carlsson and de Silva. The witness complex is defined from two sets and in some metric space : a finite set of points on which the complex is built, and a set of witnesses that serves as an approximation of . A fundamental result of de Silva states that if . In this paper, we give conditions on that ensure that the witness complex and the Delaunay triangulation coincide when is a finite set, and we introduce a new perturbation scheme to compute a perturbed set close to such that . Our perturbation algorithm is a geometric application of the Moser-Tardos constructive proof of the Lovász local lemma. The only numerical operations we use are (squared) distance comparisons (i.e., predicates of degree 2). The time-complexity of the algorithm is sublinear in . Interestingly, although the algorithm does not compute any measure of simplex quality, a lower bound on the thickness of the output simplices can be guaranteed.

Keywords.

Delaunay complex, witness complex, relaxed Delaunay complex, distance and incircle predicates, simplex quality, Lovász local lemma

1 Introduction

The witness complex was introduced by Carlsson and de Silva [cds-teuwc-04] as a weak form of the Delaunay complex that is suitable for finite metric spaces and is computed using only distance comparisons. The witness complex is defined from two sets and in some metric space : a finite set of points on which the complex is built, and a set of witnesses that serves as an approximation of . A fundamental result of de Silva [vds-wdt-08] states that if is the entire Euclidean space , and the result extends to spherical, hyperbolic and tree-like geometries. The result has also been extended to the case where is a smoothly embedded curve or surface of  [aem-ww-2007]. However, when the set of witnesses is finite, the Delaunay triangulation and the witness complexes are different and it has been an open question to understand when the two structures are identical. In this paper, we answer this question and present an algorithm to compute a Delaunay triangulation using the witness complex.

We first give conditions on that ensure that the witness complex and the Delaunay triangulation coincide when is a finite set (Section 3). Some of these conditions are purely combinatorial and easy to check. In a second part (Section 4), we show that those conditions can be satisfied by slightly perturbing the input set . Our perturbation algorithm is a geometric application of the Moser-Tardos constructive proof of the general Lovász local lemma. Its analysis uses the notion of protection of a Delaunay triangulation that we have previously introduced to study the stability of Delaunay triangulations [bdg-stabdt-2014].

Our algorithm has several interesting properties and we believe that it is a good candidate for implementation in higher dimensions.

1. Low algebraic degree.

The only numerical operations used by the algorithm are (squared) distance comparisons (i.e., predicates of degree 2). In particular, we do not use orientation or in-sphere predicates, whose degree depends on the dimension and are difficult to implement robustly in higher dimensions.

2. Efficiency.

Our algorithm constructs the witness complex of the perturbed set in time sublinear in . See Section 5.

3. Simplex quality and Delaunay stability.

Differently from all papers on this and related topics, we do not compute the volume or any measure of simplex quality. Nevertheless, through protection, a lower bound on the thickness of the output simplices can be guaranteed (see Theorem 26), and the resulting Delaunay triangulation is stable with respect to small metric or point perturbations [bdg-stabdt-2014].

4. No need for coordinates.

We can construct Delaunay triangulations of points that come from some Euclidean space but whose actual positions are unknown. We simply need to know the interpoint distances.

5. A thorough analysis.

Almost all papers in Computational Geometry rely on oracles to evaluate predicates exactly and assume that the complexity of those oracles is . Our (probabilistic) analysis is more precise. We only use predicates of degree 2 (i.e. double precision) and the analysis fully covers the case of non generic data.

1.1 Previous work

Millman and Snoeyink [ms-cpvd-2010] developed a degree-2 Voronoi diagram on a grid in the plane. The diagram of points can be computed using only double precision by a randomized incremental construction in expected time and expected space. The diagram also answers nearest neighbor queries, but it doesn’t use sufficient precision to determine a Delaunay triangulation.

Our work has some similarity with the -geometry introduced by Salesin et al. [DBLP:conf/compgeom/SalesinSG89]. For a given geometric problem, the basic idea is to compute the exact solution for a perturbed version of the input and to bound the size of this implicit perturbation. The authors applied this paradigm to answering some geometric queries such as deciding whether a point lies inside a polygon. They do not consider constructing geometric structures like convex hulls, Voronoi diagrams or Delaunay triangulations. Also, the required precision is computed on-line and is not strictly limited in advance as we do here.

Further developments have been undertaken under the name of controlled perturbation [dh-cp-2010]. The purpose is again to actually perturb the input, thereby reducing the required precision of the underlying arithmetic and avoiding explicit treatment of degenerate cases. A specific scheme for Delaunay triangulations in arbitrary dimensions has been proposed by Funke et al. [funke2005]. Their algorithm relies on a careful analysis of the usual predicates of degree and is much more demanding than ours.

1.2 Notation

In the paper, and denote sets of points in where denotes the standard flat torus . is finite but we will only assume that is closed in . The points of are called the witnesses and the points of are called the landmarks. To keep the exposition simpler, boundary issues will be handled in the full version of the paper (see Section 6).

We say that is an -sample if for any there is a with . If is a finite set, then it is a -sample for for all greater or equal to some finite . The parameter is called the sampling radius of . We further say that is -net if for all . We call the sparsity ratio of . Note that, for any two points of a -net, we have , which implies .

A simplex is a finite set. We always assume that contains a non-degenerate -simplex, and we demand that . This bound on avoids the topological complications associated with the periodic boundary conditions.

The volume of Euclidean -ball with unit radius will be denoted by .

2 Delaunay and witness complexes

Definition 1 (Delaunay center and Delaunay complex)

A Delaunay center for a simplex is a point that satisfies

The Delaunay complex of is the complex consisting of all simplexes that have a Delaunay center.

Note that is at equal distance from all the vertices of . A Delaunay simplex is top dimensional if is not the proper face of any Delaunay simplex. The affine hull of a top dimensional simplex has dimension . If is top dimensional, the Delaunay center is the circumcenter of which we denote . We write for the circumradius of .

Delaunay [delaunay1934] showed that if the point set is generic, i.e., if no empty sphere contains points on its boundary, then is a triangulation of (see the discussion in Section 3), and any perturbation of a finite set is generic with probability 1. We refer to this as Delaunay’s theorem.

We introduce now the witness complex that can be considered as a weak variant of the Delaunay complex.

Definition 2 (Witness and witness complex)

Let be a simplex with vertices in , and let be a point of . We say that is a witness of if

The witness complex is the complex consisting of all simplexes such that for any simplex , has a witness in .

Figure 1: Blue points are the witnesses and the black points are the landmarks . The simplex is witnessed by .

Observe that the only predicates involved in the construction of are (squared) distance comparisons, i.e. polynomials of degree in the coordinates of the points. This is to be compared with the predicate that decides whether a point lies inside, on or outside the sphere circumscribing a -simplex which is a polynomial of degree .

3 Identity of witness and Delaunay complexes

In this section, we make the connection between Delaunay and witness complexes more precise. We start with de Silva’s result [vds-wdt-08]:

Theorem 3

, and if then .

If is generic, we know that is embedded in by Delaunay’s theorem. It therefore follows from Theorem 3 that the same is true for . In particular, the dimension of is at most .

3.1 Identity from protection

When is not the entire space but a finite set of points, the equality between and no longer holds. However, by requiring that the -simplices of be -protected, a property introduced in [bdg-stabdt-2014], we are able to recover the inclusion , and establish the equality between the Delaunay complex and the witness complex with a discrete set of witnesses.

Definition 4 (-protection)

We say that a simplex is -protected at if

We say that is -protected when each Delaunay -simplex of has a -protected Delaunay center. In this sense, -protection is in fact a property of the point set and we also say that is -protected. If is -protected for some unspecified , we say that is protected (equivalently is generic). We always assume since it is impossible to have a larger if is a -sample. The following lemma is proved in [2014arXiv1410.7012B]. For a subcomplex , we define as the subcomplex consisting of all the simplices that have at least one vertex in (note that this definition departs from common usage). We will use as a shorthand for .

Lemma 5 (Inheritance of protection)

Let be a -net and suppose . If every -simplex in is -protected, then all simplices in are at least -protected where . Specifically, any Delaunay -simplex is -protected at a point which is the barycenter of the circumcenters of a subset of -cofaces of .

The following lemma is an easy consequence of the previous one.

Lemma 6 (Identity from protection)

Let be a -net with . If all the -simplices in are -protected and is an -sample for with , then .

Proof.

By Theorem 3, we have

We now prove the other inclusion. Let be a simplex in and set

Let be a face of . Then, by Lemma 5, is -protected at a point such that

For any , any and any , we have

and

Hence, is a witness of . Since , there must be a point in which witnesses . Since this is true for all faces , we have . ∎

We end this subsection with a result proved in [bdg-stabdt-2014, Lemma 3.13] that will be useful in Section 5. For any vertex of a simplex , the face oppposite is the face determined by the other vertices of , and is denoted by . The altitude of in is the distance from to the affine hull of . The altitude of is the minimum over all vertices of of . A poorly-shaped simplex can be characterized by the existence of a relatively small altitude. The thickness of a -simplex is the dimensionless quantity that evaluates to if and to otherwise, where denotes the diameter of , i.e. the length of its longest edge.

Lemma 7 (Thickness from protection)

Suppose is a -simplex with circumradius less than and shortest edge length greater than or equal to . If every -face of is also a face of a -protected -simplex different from , then the thickness of satisfies

In particular, suppose , where is a -net, and every -simplex in is -protected, then every -simplex in is -thick.

3.2 A combinatorial criterion for identity

The previous result will be useful in our analysis but does not help to compute from since the -protection assumption requires knowledge of . A more useful result in this context will be given in Lemma 10. Before stating the lemma, we need to introduce some terminology and, in particular, the notion of good links.

Figure 2: Pinched simplicial torus is an example of -pseudo-manifold.

A complex is a -pseudo-manifold if it is a pure -complex and every -simplex is the face of exactly two -simplices.

Definition 8 (Good links)

Let be a complex with vertex set . We say has a good link if is a -pseudo-manifold. If every has a good link, we say has good links.

For our purposes, a simplicial complex is a triangulation of if it is a -manifold embedded in . We observe that a triangulation has good links.

Lemma 9 (Pseudomanifold criterion)

If is a triangulation of and has the same vertex set, then if and only if has good links.

Proof.

Let be the common vertex set of and , and let be a point in . Since is a triangulation, is a simplicial -sphere, which is a manifold. A -pseudo-manifold cannot be properly contained in a connected -manifold simplicial complex , because a -simplex that does not belong to cannot share a -face with a simplex of : Every -simplex of is already the face of two -simplices of . Since any two -simplices in can be connected by a path confined to the relative interiors of and -simplices, it follows that a -simplex must also belong to . For an arbitrary simplex , we have that is a face of a -simplex , because is a pure -complex (it is a manifold). It follows that and therefore .

Therefore, since , if has good links we must have . It follows that for every , and therefore . If does not have good links it is clearly not equal to , which does. ∎

We can now state the lemma that is at the heart of our algorithm. It follows from Theorem 3, Lemma 9, and Delaunay’s theorem:

Lemma 10 (Identity from good links)

If is generic and the vertices of have good links, then .

4 Turning witness complexes into Delaunay complexes

Let, as before, be a finite set of landmarks and a finite set of witnesses. In this section, we intend to use Lemma 10 to construct , where is close to , using only comparisons of (squared) distances. The idea is to first construct the witness complex which is a subcomplex of (Theorem 3) that can be computed using only distance comparisons. We then check whether using the pseudomanifold criterion (Lemma 9). While there is a vertex of that has a bad link (i.e. a link that is not a pseudomanifold), we perturb and the set of vertices , to be exactly defined in Section 4.3 (See Eq. 1) that are responsible for the bad link , and recompute the witness complex for the perturbed points. We write for the set of perturbed points at some stage of the algorithm. Each point is randomly and independently taken from the so-called picking ball . Upon termination, we have . The parameter , the radius of the picking balls, must satisfy Eq. (4) to be presented later. The steps are described in more detail in Algorithm 1. The analysis of the algorithm relies on the Moser-Tardos constructive proof of Lovász local lemma.

  Input:, , , ,
  while a vertex of has a bad link do
     perturb and the points in (defined in Eq. 1)
     update
  end while
  Output:
Algorithm 1 Delaunay triangulation from witness complex

4.1 Lovász local lemma

The celebrated Lovász local lemma is a powerful tool to prove the existence of combinatorial objects [the-probabilistic-method-alon-spencer]. Let be a finite collection of “bad” events in some probability space. The lemma shows that the probability that none of these events occur is positive provided that the individual events occur with a bounded probability and there is limited dependence among them.

Lemma 11 (Lovász local lemma)

Let be a finite set of events in some probability space. Suppose that each event is independent of all but at most of the other events , and that for all . If

( is the base of the natural logarithm), then

Assume that the events depend on a finite set of mutually independent variables in a probability space. Moser and Tardos [mt-cplll-2010] gave a constructive proof of Lovász lemma leading to a simple and natural algorithm that checks whether some event is violated and randomly picks new values for the random variables on which depends. We call this a resampling of the event . Moser and Tardos proved that this simple algorithm quickly terminates, providing an assignment of the random variables that avoids all of the events in . The expected total number of resampling steps is at most .

4.2 Three geometric lemmas

We will be using the following geometric lemmas.

We will need the following technical result to prove Lemma 15.

Lemma 12

If is a -sample of , the circumradius of any simplex in is less than . The witnesses of are at a distance less than from the closest vertex of and at a distance less than from any vertex of .

Proof.

By Theorem 3, is a simplex of and, since is a -sample of , the circumradius is less than . Let be a witness of and the vertex of closest to . Since is a -sample, . The diameter of is thus at most , and for any vertex . ∎

Lemma 13

If is a -net of and , then is an -net, where

Figure 3: For Lemma 14.
Lemma 14

(See [boissonnat:hal-00806107, Lemma 5.2]) Let be a hypersphere of of radius centered at and the spherical shell . Let, in addition, denote any -ball of radius . We have

where denotes the volume of a Euclidean -ball with unit radius.

4.3 Correctness of the algorithm

We first bound the density radius and the sparsity ratio of any perturbed point set . We will write and , and assume .

We refer to the terminology of the Lovász local lemma. Our variables are the points of which are randomly and independently taken from the picking balls , .

The events are associated to points of , the vertices of . We say that an event happens at when the link of in is not good, i.e., is not a pseudomanifold. We know from Lemma 6 that if is a vertex of and is not good, then there must exist a -simplex in that is not -protected for . We will denote by

  • the set of points of that can be in

  • the set of points of that can violate the -protected zone for some -simplex in

  • Since computing and exactly is difficult, we will rather compute the set

    (1)

    Using triangle inequality, Lemma 13, and the facts that and , we can show that

    where is the set of points in that correspond to points in .

  • the set of -simplices with vertices in that can belong to

The probability that is not good is at most the probability that one of the simplices of , say , has its -protecting zone violated by some point of . Write for the probability that belongs to the -protection zone of the -simplex . We have

(2)

The following lemmas provide upper bounds on , , and .

Lemma 15

and .

Proof.

Since any two points in are at least apart, a volume argument then gives the announced bound. Referring to (1) we have

We thus have

Lemma 16

An event is independent of all but at most other events.

Proof.

Two events and are independent if . Referring to defined in (1), this implies and are independent if . Hence an upper bound on the number of events that may not be independent from the event at is given by a volume argument as in the proof of Lemma 15. Here we count points since the events are associated to points:

Lemma 17

.

Proof.

The volume of Euclidean -ball with unit radius will be denoted by . The probability is the ratio between the volume of the portion of the picking region that is in the spherical shell and the volume of the picking region. Hence, using [boissonnat:hal-00806107, Lem. 5.2] and the crude bound , we obtain

Using Eq. (2), and Lemmas 15 and 17, we conclude that

An event depends on at most other events. Hence, to apply the Lovász Local Lemma 11, it remains to ensure that . In addition, we also need that to be able to apply Lemma 6. We thefore need to satisfy the following inequality:

(3)

Observe that , , and depend only on and . We conclude that the conditions of the Lovász local lemma hold if the parameter satisfies

(4)

Hence, if is sufficiently small, we can fix so that Eq. (4) holds. The algorithm is then guaranteed to terminate. By Lemma 10, the output is .

It follows from Moser-Tardos theorem that the expected number of times the "while-loop" in Algorithm 1 runs is and since , we get that the number of point perturbations performed by Algorithm 1 is on expectation, where the constant in the depends only on and not on the sampling parameters or . We sum up the results of this section in

Theorem 18

Under Eq. (4), Algorithm 1 terminates and outputs the Delaunay triangulation of some set whose distance to is at most . The number of point perturbations performed by the algorithm is .

5 Sublinear algorithm

When the set is generic, is embedded in and is therefore -dimensional. It is well known that the -skeleton of can be computed in time using only distance comparisons [DBLP:journals/algorithmica/BoissonnatM14]. Although easy and general, this construction is not efficient when is large.

In this section, we show how to implement an algorithm with execution time sublinear in . We will assume that the points of are located at the centers of a grid, which is no real loss of generality. The idea is to restrict our attention to a subset of , namely the set of full-leaf-points introduced in Section 5.2. These are points that may be close to the circumcentre of some -simplex. A crucial observation is that if a -simplex has a bounded thickness, then we can efficiently compute a bound on the number of its full-leaf-points. This observation will also allow us to guarantee some protection (and therefore thickness) on the output simplices, as stated in Theorem 26 below.

The approach is based on the relaxed Delaunay complex, which is related to the witness complex, and was also introduced by de Silva [vds-wdt-08]. We first introduce this, and the structural observations on which the algorithm is based.

5.1 The relaxed Delaunay complex

The basic idea used to get an algorithm sublinear in is to choose witnesses for -simplices that are close to being circumcentres for these simplices. With this approach, we can in fact avoid looking for witnesses of the lower dimensional simplices. The complex that we will be computing is a subcomplex of a relaxed Delaunay complex:

Definition 19 (Relaxed Delaunay complex)

Let be a simplex. An -center of is any point such that

We say that is an -Delaunay center of if

The set of simplices that have an -Delaunay centre in is a simplicial complex, called the -relaxed Delaunay complex, and is denoted .

We say is an -witness for if

We observe that is an -Delaunay center if and only if it is an -center and also an -witness.

Lemma 20

The distance between an -Delaunay center for and the farthest vertex in is less than . In particular, .

If and is a Delaunay center of , then any point in is a -Delaunay centre for . Thus .

If, for some , all the -simplices in have a -protected circumcentre, then we have that , and with , it follows (Lemma 20) that , and is itself -protected and has good links.

Reviewing the analysis of the Moser-Tardos algorithm of Section 4.3, we observe that the exact same estimate of that serves as an upper bound on the probability that one of the simplices in is not -protected at its circumcenter, also serves as an upper bound on the probability that one of the simplices in is not -protected at its circumcenter, provided that (using the diameter bound of from Lemma 20), which we will assume from now on. We can therefore modify Algorithm 1 by replacing by .

We now describe how to improve this algorithm to make it efficient. For our purposes it will be sufficient to set . In order to obtain an algorithm sublinear in , we will not compute the full but only a subcomplex we call . The exact definition of will be given in Section 5.3, but the idea is to only consider -simplices that show the properties of being -thick for some parameter to be defined later (Eq. 5 and 6). This will allow us to restrict our attention to points of that lie near the circumcentre. As explained in Section 5.3, this is done without explicitly computing thickness or circumcentres.

As will be shown in Section 5.4 (Lemma 24), the modification of Algorithm 1 that computes instead of will terminate and output a complex with good links. However, this is not sufficient to guarantee that the output is correct, i.e., that . In order to obtain this guarantee, we insert an extra procedure check(), which, without affecting the termination guarantee, will ensure that the simplices of have -protected circumcenters for a positive . It follows that and, by Lemma 9, that . The modification of Algorithm 1 that we will use is outlined in Algorithm 2.

  input:, , , , ,
  
  compute:
  while a vertex of has a bad link or check do
     perturb and the points in  ( is defined in Section 4.3.)
     update
  end while
  output: with protection (See Lemma 25)
Algorithm 2 Protected Delaunay triangulation from

We describe the details of computing and of the check() procedure in the following subsections.

5.2 Computing relaxed Delaunay centers

We observe that the -Delaunay centers of a -simplex are close to the circumcenter of , provided that has a bounded thickness:

Lemma 21 (Clustered -Delaunay centers)

Assume that is a -sample. Let be a non degenerate -simplex, and let be an -center for at distance at most from the vertices of , for some constant . Then is at distance at most from the circumcenter of . In particular, if is an -Delaunay center for , then

Proof.

Using [bdg-stabdt-2014, Lemma 4.3] and , we get

The second assertion follows from this together with Lemma 20, and the assumption that . ∎

It follows from Lemma 21 that -Delaunay centers are close to all the bisecting hyperplanes

The next simple lemma asserts a kind of qualitative converse:

Lemma 22

Let be a -simplex and be the bisecting hyperplane of and . A point that satisfies , for any is a -center of .

Let be a -simplex of and let be the smallest box with edges parallel to the coordinate axes that contains . Then the edges of have length at most (Lemma 20). Any -Delaunay center for is at a distance at most from . Therefore all the -Delaunay centers for lie in an axis-aligned hypercube with the same center as and with side length at most . Observe that the diameter (diagonal) of is at most .

Our strategy is to first compute the -centers of that belong to and then to determine which ones are -witnesses for . Deciding if an -center is an -witness for can be done in constant time since is a -net and (Lemma 20).

We take . To compute the -Delaunay centers of , we will use a pyramid data structure. The pyramid consists of at most levels. Each level is a grid of resolution . The grid at level consists of the single cell, . Each node of the pyramid is associated to a cell of a grid. The children of a node correspond to a subdivision into subcells (of the same size) of the cell associated to . The leaves are associated to the cells of the finest grid whose cells have diameter .

A node of the pyramid that is intersected by all the bisecting hyperplanes of will be called a full node or, equivalently, a full cell. By our definition of , a cell of the finest grid contains an element of at its centre. The full-leaf-points are the elements of associated to full cells at the finest level. By Lemma 22, the full-leaf-points are -centers for . In order to identify the full-leaf cells, we traverse the full nodes of the pyramid starting from the root. Note that to decide if a cell is full, we only have to decide if two corners of a cell are on opposite sides of a bisecting hyperplane, which reduces to evaluating a polynomial of degree 2 in the input variables.

We will compute a bound on the number of full cells in the pyramid. Consider a level in the pyramid and the associated grid whose cells have diameter . Any point in a full cell of is at distance from the bisecting hyperplanes of and so, by Lemma 22, is a -center for . By construction of , the distance between a point and a vertex is bounded by

Therefore, it follows from Lemma 21 that all the full cells of are contained in a ball of radius centered at . Letting be a hypercube with the dimensions of a cell in , a volume argument gives a bound on , the number of full cells in :

where is the volume of the Euclidean ball of unit radius.

Observe that does not depend on . Hence we can apply this bound to all levels of the pyramid, from the root level to the leaves level. Since there are at most levels, we conclude with the following lemma:

Lemma 23

The number of full cells is less than

is also a bound on the time to compute the full cells.

5.3 Construction of

By Lemma 20, all the simplices incident to a vertex of are contained in , and it follows from the fact that is a -net that

In the first step of the algorithm, we compute, for each , the set , and the set of -simplices