A Simple Algorithm for Higher-Order Delaunay Mosaics and Alpha Shapes

We present a simple algorithm for computing higher-order Delaunay mosaics that works in Euclidean spaces of any finite dimensions. The algorithm selects the vertices of the order-k mosaic from incrementally constructed lower-order mosaics and uses an algorithm for weighted first-order Delaunay mosaics as a black-box to construct the order-k mosaic from its vertices. Beyond this black-box, the algorithm uses only combinatorial operations, thus facilitating easy implementation. We extend this algorithm to compute higher-order \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}α-shapes and provide open-source implementations. We present experimental results for properties of higher-order Delaunay mosaics of random point sets.


Introduction
Order-k Voronoi tessellations are a generalization of ordinary Voronoi tessellations.Instead of each domain corresponding to a single point in the input, A ⊆ R d , each order-k domain corresponds to a subset, Q ⊆ A, of size k, and consists of the set of points x ∈ R d for whom the points in Q are the closest k points within A. Its dual is the order-k Delaunay mosaic.We will formally define both in Section 2. Order-k Voronoi tessellations were introduced by [21] as a data structure for fast k closest point queries, namely in time O(k + log n) with n = #A.A less direct application is the computation of the distance-to-measure introduced in [5] and related to k closest point search in [13].Furthermore, certain subcomplexes of the order-k Delaunay mosaic realize the order-k α-shapes introduced in [14].Order-k α-shapes are a generalization of α-shapes [9] used to approximate the shape of a point set.Unlike ordinary α-shapes and depending on the parameter k, they exhibit robustness to noisy data points.
In the plane, the number of domains in the order-k Voronoi tessellation or, equivalently, the number of vertices in the order-k Delaunay mosaic is Θ(k(n − k)); see [15,21].For dimensions d ≥ 3, this number can vary significantly depending on the way the input points are distributed.The upper bound of O(k ) on the total size of the first k higher-order Delaunay mosaics [6] is tight, while the lower bound of Ω(k d n) [17] is only conjectured.For individual order-k Delaunay mosaics, the complexity is poorly understood.The problem is closely related to the (d + 1)-dimensional k-set problem.Specifically, the points in A ⊆ R d can be mapped to equally many points in R d+1 such that the order-k domains in R d correspond to k-sets in R d+1 , see e.g.[6].
The first algorithm to compute order-k Voronoi tessellations and Delaunay mosaics in the plane was described by Lee in [15].The algorithm computes the Voronoi tessellations one by one, in increasing order and in time O(k 2 n log n).Mulmuley [17] extended this algorithm beyond two dimensions, computing the first k levels in a special (d + 1)-dimensional hyperplane arrangement, which implicitly yield the order-k Voronoi tessellations and Delaunay mosaics in time O(s log n + k d n 2 ), in which s denotes the output size.Mulmuley [18] later described another algorithm, which instead adds hyperplanes one by one, and runs in time ) for d ≥ 3, which equals the worst-case output size.For d = 2, the expected runtime is O(k 2 n log n k ).Another incremental algorithm with similar complexity for d ≥ 3 has been described by Agarwal et al. [1].
In this paper, we describe a new algorithm for computing order-k Delaunay mosaics in Euclidean space of any finite dimension that stands out in its simplicity.It employs an algorithm for weighted first-order Delaunay mosaics, and otherwise uses only combinatorial operations.It thus benefits from highly optimized existing implementations and, if desired, can build upon their use of exact arithmetic.Its complexity depends on the complexity of the algorithm used for weighted Delaunay mosaics.Assuming it is linear in its output size, then out algorithm is also linear in its output size.We implement this algorithm and run it on various point sets, shedding light on the size and other properties of order-k Delaunay mosaics.In particular, we compare the total size of the first k Delaunay mosaics of random point sets with the (tight) worst-case upper bound, and we study the size of individual order-k Delaunay mosaics, for which no tight bounds are known in general.As far as we are aware, no such experimental investigations have been performed in the past, possibly due to the lack of a practical algorithm.We extend our algorithm to compute the radius function on an order-k Delaunay mosaic, which gives us the subcomplexes realizing order-k α-shapes.Open-source implementations of our algorithm are available [19,20].
Our algorithm makes use of the rhomboid tiling [10], which we will introduce in Section 2 alongside other necessary definitions.We will explore the combinatorial properties of this tiling and, by proxy, the properties of order-k Delaunay mosaics in Section 3. Using these results, we explain our algorithm in Section 4. We present experimental results obtained with two implementations of this algorithm in Section 5. Section 6 introduces a radius function on the order-k Delaunay mosaics and a way to compute it to yield order-k α-shapes.We close with a discussion of possible extensions and optimizations in Section 7.

Definitions
Given a locally finite set, The order-k Delaunay mosaic is the cell complex dual to Vor k (A), denoted Del k (A).To realize the mosaic geometrically, we usually use the average of the points in Q as the location of the corresponding vertex in R d .In a few instances we use the sum rather than the average, for convenience.Figure 1 shows an example for k = 2.As we will see shortly, in d ≥ 3 dimensions, the order-k Delaunay mosaic is not necessarily simplicial even if the points are in general position.Assuming A is in general position, [10] established the existence of a tiling in R d+1 whose horizontal integer slices are the Delaunay mosaics.We recall the definition of the tiling and its most important properties.Let A ⊆ R d be locally finite and in general position.We construct a rhomboid, , for each partition A = A in A on A out for which there exists a sphere S in R d such that all points in A in lie inside S, all points in A on lie on S, and all points in A out lie outside S. Whenever convenient, we write A in ( ) = A in , A on ( ) = A on , and A out ( ) = A out to indicate the correspondence.Due to general position of A, we have #A on ≤ d + 1.A combinatorial vertex of is a collection of points that contains A in and is contained in A in ∪ A on , and we write for the collection of combinatorial vertices of .Setting y a = (a, −1) ∈ R d+1 , for every a ∈ A, we draw the rhomboids in R d+1 by mapping every combinatorial vertex to y Q = q∈Q y q , in which y ∅ = 0, by convention.The (d + 1)-st coordinate of y Q is therefore −#Q, and we call #Q the depth of the vertex.The geometric realization of a rhomboid is the convex hull of the locations of its combinatorial vertices, which is a rhomboid.We refer to A in ( ) as the anchor vertex of .The rhomboid tiling of A, denoted Rho(A), is the collection of thus defined rhomboids.By assumption of general position, every face of a rhomboid is again defined by a sphere as described above and thus belongs to the rhomboid tiling.As proved in [10], any two rhomboids are either disjoint or intersect in a common face, which implies that the rhomboid tiling is a complex embedded in R d+1 ; see Figure 2 for an example.The following properties have been observed in [10]: Proposition 1 (Rhomboid Tiling).Let A ⊆ R d be locally finite and in general position.1. Rho(A) is dual to an arrangement of hyperplanes in R d+1 ; 2. the slice of Rho(A) at depth k is the order-k Delaunay mosaic of A, scaled by a factor k.
The hyperplane arrangement will be introduced in Section 3.2.We elaborate on the second property: that each cell of the order-k Delaunay mosaic is a slice of some rhomboid.Combinatorially, each rhomboid is a cube and, again combinatorially, each cell of Del k (A) is a slice orthogonal to the cube diagonal that passes through a non-empty set of the vertices.For the (d + 1)-cube, there are d + 2 such slices, which we index from top to bottom by the generation 0 ≤ g ≤ d + 1.The g-th slice passes through d+1 g vertices, so we have a vertex at generations g = 0, d + 1, a d-simplex at generations g = 1, d, and some other d-dimensional polytope at generations 2 ≤ g ≤ d − 1.In d + 1 = 3 dimensions, we have a vertex, a triangle, another triangle, and another vertex, see Figure 4; but already in d + 1 = 4 dimensions, the middle slice is not a simplex; see Figure 3.We remark that in addition to the order-k Delaunay mosaic, also the degree-k Delaunay mosaic, which is the dual of the degree-k Voronoi tessellation [11], can be obtained as a slice of Rho(A) at depth k − 1 2 .Figure 2 Left: the rhomboid tiling of five points in R 1 .The highlighted rhomboid defined by Ain = {c} and Aon = {b, d} is the convex hull of the points yc, y {b,c} , y {c,d} , and y {b,c,d} .The horizontal line at depth k intersects the tiling in the order-k Delaunay mosaic.Right: the dual hyperplane arrangement.Following the dotted lines connecting the points of A on the horizontal axis to the paraboloid, we find the corresponding tangent hyperplanes.The highlighted rhomboids of dimension j = 0, 1, 2 are dual to the highlighted cells of dimension 2 − j in the arrangement.

Combinatorial Properties
As proved in [4], the order-k Delaunay mosaic is the projection of the boundary complex of a convex polyhedron in R d+1 .To explain this construction, we define the lift of a ∈ R d as the point lift(a) = (a, a 2 ) ∈ R d+1 .For each k-tuple Q ⊆ A, we take the barycenter of their lifts, 1 k q∈Q lift(q), and obtain the order-k Delaunay mosaic as the vertical projection of the lower faces of the convex hull of these barycenters.Equivalently, we can interpret each barycenter of lifts as a weighted point in R d and get the order-k Delaunay mosaic as the weighted order-1 Delaunay mosaic of the weighted points.By itself, this approach does not scale well with k since there are #A k such barycenters.Most barycenters, however, are irrelevant as they do not contribute to the lower faces of the convex hull.If we could, somehow, identify the relevant barycenters without wasting time on the irrelevant ones, this procedure would efficiently construct the cells of the order-k Delaunay mosaic by computing the weighted first-order Delaunay mosaic.We will see how this can be done in Section 3.2.
In d ≥ 3 dimension, not all cells of Del k (A) are simplicial, even if the points in A are in general position.The cells carry important information, which for some applications is essential and cannot be easily recovered from a triangulation.This poses an additional challenge because most algorithms for computing convex hulls or weighted first-order Delaunay mosaics return a triangulated version of the correct mosaic.As explained in the following section, we address this issue by predicting the cells from their corresponding rhomboids.

Predicting Cells
Given a cell σ in the order-k Delaunay mosaic, the following lemma identifies the rhomboid, , that σ is a slice of; see Figure 4 for an illustration.Write V (σ) for the set of combinatorial vertices whose locations are the vertices of σ.Clearly, V (σ) ⊆ V ( ).

Lemma 2. Let ∈ Rho(A) be a rhomboid and σ ∈ Del
Since the depth of a vertex is determined by its cardinality, and the vertices of a slice are by definition all at the same depth, the vertices of the generation-g slice all satisfy #Q − #A in ( ) = g.The intersection of all g-subsets of A on ( ) is of course empty, which implies that the intersection of the combinatorial vertices of the slice is A in ( ).Furthermore, V (σ) = V ( ) for every slice σ of with generation g ≥ 1.The union of all g-subsets of A on ( ) is A on ( ) itself, and thus A on ( ) = V (σ) \ A in ( ).Finally, the generation of σ is the difference in depth of the anchor vertex, A in ( ), and the slice defining σ.The depth of σ is k and the depth of A in ( ) is its cardinality, which completes the proof.
If all of our order-k Delaunay cells are triangulated-e.g.due to being the output of a weighted first-order Delaunay algorithm-we cannot directly apply Lemma 2. Indeed, if τ is a simplex that is part of a triangulation of a slice σ of a rhomboid , then V (τ ) and V (τ ) do not necessarily equal A in ( ) and A in ( ) ∪ A on ( ).We can, however, still identify whether τ is a first-generation slice of and thus in fact is equal to σ.Using Lemma 2, we can then obtain .Proof.Let σ be the d-cell in Del k (A) that contains τ in its triangulation, and assume σ is a generation-g slice of .From Lemma 2, we know that A in ( ) ⊆ v for all v ∈ V (σ), and #A in (σ) = k − g.The remaining g points in every v are from A on ( ).We have V (τ ) ⊆ V (σ) with #V (τ ) = d + 1.So for τ to consist of vertices whose common intersection is of size k − 1, there need to be d + 1 distinct g-subsets of A on ( ) that all have g − 1 points in common.However, as #A on ( ) = d + 1, this is not possible unless g = 1.

Identifying Vertices
Given a triangulation of the order-k Delaunay mosaic, we just saw how to identify its first-generation cells.From these, we can obtain the corresponding rhomboids and their higher-generation slices.We shall now prove that if we have triangulations of the order-j Delaunay mosaics, for all j < k, we can assemble the complete vertex set of the order-k Delaunay mosaic by taking slices at depth k obtained from first-generation cells at lower depths.We note that this only holds in the unweighted setting.
To prepare the proof of this result, we recall the definition of the hyperplane arrangement postulated by Proposition 1.For each point a ∈ A, write f a : R d → R for the affine map defined by The collection of such hyperplanes decomposes R d+1 into convex cells, which we call the hyperplane arrangement of A, denoted Arr(A); see Figure 2. The cells in the arrangement are intersections of hyperplanes and closed half-spaces.More formally, for each cell there is an ordered three-partition A = A in A on A out such that the cell consists of all points This three-partition is the key to establishing the bijection between the cells of Arr(A) and the rhomboids of Rho(A) that proves the duality claimed in Proposition 1.We call top-dimensional cells of Arr(A) chambers; they satisfy A on = ∅.The depth of a chamber is #A in or, equivalently, the number of hyperplanes that are above this chamber; it equals the depth of the dual vertex in Rho(A).To see the aforementioned relationship between the arrangement and the higher-order Voronoi tessellations, we observe that the chamber in Arr(A) with A in = Q vertically projects to dom(Q).We can therefore construct Vor k (A) by computing and projecting all chambers whose ordered three-partitions satisfy #A in = k; see [8, Chapter 13] or [11].
We call a chamber γ a bowl if only one of its facets bounds it from above or, equivalently, if there is only one chamber at the next lower depth that shares a facet with γ.We call the hyperplane that contains this facet the lid of the bowl.Lemma 4. A hyperplane that is a lid of a bowl at depth 1 is not a lid of any other chambers.
Proof.Let γ be a bowl at arbitrary depth, and let P be its lid.Every other hyperplane that contains a facet of γ bounds γ from below.The top facet of γ is the only part of P that is above all of these hyperplanes; that is: all other parts of P are below at least one of the other hyperplanes.This implies that every other bowl with lid P has at least one other hyperplane above it, and is thus of depth at least 2. Now assume γ is at depth 1.If there were another bowl γ with lid P , then the above argument would yield that all other bowls are at depth at least 2, contradicting our assumption on γ.Thus γ has to be the unique bowl with lid P .
With this lemma, we are ready to state and prove the main combinatorial insight that motivates our algorithm.In a nutshell, it says that the first-generation cells form clusters without interior vertices.In R 2 , this is equivalent to saying that these clusters have outerplanar 1-skeletons.
Theorem 5. Let A ⊆ R d be locally finite and k ≥ 2. Then every vertex in Del k (A) is vertex of some d-cell of generation g ≥ 2.
Proof.In the unweighted setting, each hyperplane is tangent to the paraboloid P and contains a facet of the unique depth-0 chamber.Thus, each hyperplane is the lid to a chamber at depth 1.As this is true for every hyperplane, all chambers of depth 2 or higher have no lids by Lemma 4. This means that any chamber of depth at least 2 has at least two upper facets.Because the upper boundary is connected, there are two upper facets that meet in a (d − 1)-face, the dual rhomboid of this face has dimension 2, and its bottom vertex is dual to the chamber.Thus we can obtain this vertex, v, knowing the other three vertices of the 2-rhomboid.
Any 2-rhomboid is a face of some (d + 1)-dimensional rhomboid, , which thus contains v at generation at least 2, i.e. v has depth at least #A in ( ) + 2. Knowing A in ( ) and A on ( ), we obtain this vertex via Equation (1).

Algorithm
We outline our algorithm in this section; its correctness follows from the results of the previous sections.We compute the Delaunay mosaics one by one in sequence of increasing order.For Del 1 (A), the vertex set is the set A of input points.Whenever we have the vertex set of Del j (A), we compute its (triangulated) d-cells using an off-the-shelf algorithm for weighted Delaunay triangulations.We use Lemma 3 to identify the first-generation d-cells, while discarding all other cells.From each first-generation cell, we obtain and save the higher-generation d-cells and vertices defined by the same rhomboid.These will appear in Delaunay mosaics of higher orders.By Theorem 5, once we have processed all Del j (A) for j < k, we will have obtained the complete vertex set of Del k (A) in the process, thus allowing our algorithm to continue until we have all Delaunay mosaics up to the desired order.Algorithm 1 is a more formal write-up of the above outline, and Figure 4 visualizes the process.A dimension-agnostic python implementation and a 2-and 3-dimensional C++ implementation using CGAL [22] are available at [19,20].If we store all first-generation cells with their anchor vertices, we can use this algorithm to implicitly construct the rhomboid tiling, as done in [20].
To get a handle on the runtime of the algorithm, we consider the two steps used to compute the order-k Delaunay mosaic after finishing the construction of the first k −1 mosaics.The first step is geometric and invokes the black-box algorithm to construct the weighted Delaunay mosaic from which we get vertices and cells of (unweighted) higher-order Delaunay mosaics.The runtime of this step depends on the runtime of the black box algorithm, which in many cases is output-dependent.The second step is combinatorial and determines, for each output simplex from the first step, whether it is first generation, in which case it is a genuine cell of the mosaic.Assuming constant dimension, d, identifying whether an order-(k − 1) cell is of first generation and, in this case, obtaining the higher-generation cells takes time O(k).Thus, for a given k, the combinatorial step takes time O(kC k ), in which C k is the number of d-cells of Del k (A).With each vertex being represented as a k-tuple of points, this is linear in the output size, assuming we store each cell naively as a set of its vertices.If the runtime of each black-box invocation were linear in the output size, the total runtime for producing the first k higher-order Delaunay mosaics would thus be linear in the output size as well.
In practice, it is more efficient to store a cell as a set of pointers or indices to its vertices, only requiring space O(kV k + C k ), with V k denoting the number of vertices of Del k (A).Using this representation, the combinatorial step is not linear in the output size unless the number of cells of Del k (A) is linear in the number of vertices.

Experimental Results
In 2 dimensions, the number of cells in the (order-1) Delaunay mosaic is always linear in the number of input points, while in d ≥ 3 dimensions, the size of the mosaic depends on the input set itself-and not just its cardinality-and ranges from Ω(n) to O(n Size in 3 dimensions.To shed light on the size range of order-k Delaunay mosaics, we compute them for a few 3-dimensional point sets relevant to these bounds.Note that for order-k Delaunay mosaics the number of vertices varies as well.Figure 5 shows the numbers of vertices and 3-dimensional cells for all higher-order Delaunay mosaics of four sets of size n = 200 each: points on the moment curve, points sampled on the torus (with major radius 1 and minor radius 0.5 obtained by uniformly sampling the angles of its parametrization), points uniformly sampled inside the unit ball, and a point set in convex position forming each other, with roughly three times as many cells as vertices.Other than in Figure 5, we therefore omit the information about the vertices and show only the plots for the cells.The moment curve and polytope sets are both in convex position.Nevertheless, the size of the mosaic for the moment curve grows large faster for small k, and reaches its peak at k ≈ n/3, while for the polytope the peak is at k ≈ n/2.Notice how a faster rise also goes along with an earlier decay.This is a consequence of the fact that the total size of all order-k Delaunay mosaics together-or, equivalently of the rhomboid tiling-only depends on the input size, n, and not on the relative position of the input points.
Size increase for small order.Looking more closely at the growth for small k relative to the input size, we observe that the polytope and unit ball exhibit linear growth while the size of the mosaic seems to grow quadratically for the moment curve, see Figure 6.This is consistent with the bounds on first-order Delaunay mosaics mentioned earlier.For the torus, the size seems to grow slightly superlinearly, which is again consistent with the O(n log n) bound for smooth surfaces mentioned above.
Variance.To probe whether the above figures are representative, we investigate the variance in number of cells for the polytope and the unit ball.As shown in Figure 7, the variance is particularly small for the polytope, and it is considerably larger of the unit ball.Curiously, the variance dips at k = n/2.Generations.We also investigate the distribution of cells of different generations.All point sets exhibit a pattern similar to that in Figure 8, with the fraction of first-generation cells decreasing and the fraction of d-th-generation cells increasing as the order grows.The change is most prominent for small and large k, while the fractions remain almost constant in the range k ≈ n/2, provided n is significantly larger than the dimension d.
Curse of dimensionality.Like many geometric structures, order-k Delaunay mosaics are subject to the dimensionality curse.Figure 9 shows how the size of order-k Delaunay mosaics behaves for point sets in different dimensions.
Vertex degrees.Order-k Delaunay mosaics in R 3 exhibit an interesting distribution of vertex degrees for random point sets; see Figure 10.The distribution looks like the sum of two distributions-with the second one only covering values 2 modulo 3-and is exhibited for  all k except very small and very large ones.We do not know the reason for vertices being frequently incident to 5, 8, 11, . . .d-cells, but suspect these numbers correspond to geometric configurations of cells of different generations, such as three octahedra sharing a common vertex with two tetrahedra.
Clusters.First-generation cells of any order-k Delaunay mosaic come in clusters connected by shared facets.We investigate the distribution of their sizes, leaving the discussion of their potential algorithmic significance for later.Figure 11 shows cluster size distributions in R 3 for different orders.For very small k, the distribution depends on how the points are sampled, while for all other k, the cluster sizes seem to follow an exponential distribution.The decay rate increases with k and seems to be linked to the fraction of first-generation cells.It culminates in all clusters being singletons for k = n − 3.For k > n − 3, there are no more first-generation cells.

Order-k Alpha Shapes
Beyond order-k Delaunay mosaics, our algorithm can be extended to compute order-k alpha shapes, as introduced in [14].To this end, the rhomboid tiling is endowed with a radius function on its rhomboids [10].It is inherited by the Delaunay mosaics, which are slices of the rhomboid tiling, and their sublevel sets with respect to this radius function are complexes that geometrically realize the order-k α-shapes.In this section, we recall the definition of the radius function from [10], and present an efficient way of computing it.
To get started, we note that the radius function needs a representation for every rhomboid in the tiling, but the algorithm in Section 4 computes only the top-dimensional rhomboids.This is easily remedied by noticing that the dimension of a rhomboid is k = #A on ( ) and its 3 k faces correspond to the different ways of partitioning A on ( ) into three sets.For the remainder of the discussion, assume that we have a representation for the rhomboids of all dimensions 0 ≤ j ≤ d + 1 in Rho(A).Each j-dimensional rhomboid, ∈ Rho(A), corresponds to a (d + 1 − j)-dimensional cell in the dual arrangement, * ∈ Arr(A).We introduce With slight abuse of notation, we write P t for the graph of this function.This graph is the paraboloid P 0 dropped down vertically by a distance t 2 .We define the squared radius function R 2 : Rho(A) → R, which maps a rhomboid to the minimum t such that P t has a non-empty intersection with * .We call a sphere constrained by if it encloses A in ( ), passes through all points of A on ( ), and has no other points of A inside.Letting S min ( ) be the smallest such sphere, we get an alternative interpretation of the radius function: Lemma 6. R 2 ( ) equals the squared radius of S min ( ).
Proof.The proof of Theorem 1 of [10] establishes a map from points in the Arr(A) to spheres: a point y = (x, z) ∈ R d × R below the paraboloid P is mapped to the sphere, S, with center x and squared radius x 2 − 2z.Importantly, if * is the cell in the dual arrangement whose interior contains y, then S is constrained by , which is the rhomboid dual to * .Now let t = R 2 ( ), and let r 2 be the squared radius of S min ( ).By definition, t is the smallest value for which P t contains a point y ∈ * .The aforementioned map maps y to a sphere constrained by , thus r 2 ≤ t.When reversing this map, S min ( ) is mapped to a point of * .As t was the smallest value for which P t touches * , we have t ≤ r 2 .Thus the squared radius of S min ( ) equals R 2 ( ).
To compute this radius function, we first get the smallest sphere constrained by a rhomboid.While Welzl's algorithm [23] for smallest enclosing sphere can be adapted to this task, it takes O(n) with n = #A for each such sphere computation.To improve on this bound, we recall that Lemma 3 of [10] establishes that rhomboids of the same radius value come in intervals [ min , max ] := { ∈ Rho(A) | min ⊆ ⊆ max } whose lower bound, min , is a vertex.To identify the vertex v that a rhomboid forms an interval with, we need to identify its vertex with the same radius value.By Lemma 6 this means the radii of S min ( ) and S min (v) have to be the same, and it is not difficult to see that the spheres S min ( ) and S min (v) are in fact the same.As A on (v) = ∅ for any vertex v, the sphere achieving the radius value of v is defined solely by inclusions and exclusion constraints.Therefore all constraints of that require points of A on ( ) to be on the sphere need to be converted to inclusion and exclusion constraints without affecting the resulting sphere.We know that such constraints exist because the lower bound of the interval is a vertex.This observation gives rise to the following lemma.
Lemma 7. Let be a rhomboid that is an upper bound of an interval.Let A I ⊆ A on ( ) such that the smallest enclosing sphere S of A I that excludes A on ( ) \ A I is the same as the circumsphere of A on ( ).Then forms an interval with the vertex v = A in ( ) ∪ A I .
Proof.As is an upper bound of an interval, its sphere, S min ( ), is only supported by A on ( ).Indeed, if there were another point a ∈ A in ( )-or a ∈ A out ( )-on the surface of this sphere, then the rhomboid with A on ( ) = A on ( ) ∪ {a} and A in ( ) = A in ( ) \ {a}-or A out ( ) = A out ( ) \ {a}-would be a higher-dimensional rhomboid with the same sphere S min ( ) = S min ( ), contradicting that be an upper bound of an interval.
As S min ( ) is only supported by A on ( ), this means that S min ( ) is the same as the circumsphere of A on ( ), which by our assumption is the same as S. Now the inclusion and exclusion constraints of S are part of the constraint set for S min (v), but because S = S min ( ) it does in fact fulfill all the constraints of S min (v).Thus S min (v) = S = S min ( ), proving that they are in the same interval.
Algorithm.Assume is a j-rhomboid that is an upper bound of an interval.Let S be the circumsphere of A on ( ).For each point a ∈ A on ( ), we need to decide whether to impose an inclusion or exclusion constraint on it.Let S a be the circumsphere of A on ( ) \ {a}.If a is outside of S a , then imposing an exclusion constraint for a would yield S a rather than S, thus we add a to A I in order to impose an inclusion constraint for it.Similarly, if a is inside of S a , we have to impose an exclusion constraint for a and thus do not add it to A I .
While this is difficult for an individual rhomboid, it becomes straightforward if we compute all intervals in the rhomboid tiling.We know that all (d + 1)-rhomboids are upper bounds of intervals.After marking all rhomboids that are contained in such intervals, we know that all remaining unmarked d-rhomboids are upper bounds of intervals.Thus by processing the rhomboids in decreasing dimension, all unmarked rhomboids we encounter are upper bounds.

Discussion
This paper presents a simple algorithm for computing order-k Delaunay mosaics in Euclidean space of constant dimension.Implementations of the algorithm-in C++ for points in R 2 and R 3 and in python for points in R d -are provided [19,20].This software includes the application to the persistence of k-fold covers described in [10].The remainder of this section discusses this application and possible extensions and optimizations of our algorithm.
k-fold covers.The sublevel sets of the order-k Delaunay mosaics with respect to the radius function introduced in Section 6 are homotopy equivalent to k-fold covers of Euclidean balls.It follows that our algorithms facilitate the computation of persistence of these k-fold covers.Furthermore, the circumcenters of the spheres that are used in the computation of the radius function provide the geometric locations of the order-k Voronoi vertices and allow reconstructing the order-k Voronoi tessellation via duality.
Weighted setting.Our algorithm generalizes to points with real weights, but not easily.The main challenge is the extraction of the vertices of the order-k mosaic from lower-order mosaics.This extraction relies on Theorem 5, which does not hold for weighted points.Indeed, a crucial assumption in this theorem is that every lifted hyperplane is incident to the depth-0 chamber of the arrangement, and this property is generally violated for weighted points.This is the same assumption used in the prior dimension-agnostic algorithms [1,17,18].For sets of weighted points that satisfy this assumption, our algorithm and these prior algorithms still work.To overcome this limitation, we would need a way to detect all bowls in the arrangement, because they correspond to the vertices in the Delaunay mosaics our algorithm is not able to find.Identifying the bowls is an independent problem, and any solution to it can be combined with our algorithm.Once we know the bowls and add the corresponding vertices to the appropriate mosaics, our algorithm works as before.

Clusters of cells.
As mentioned in Section 5, first-generation cells in the order-k Delaunay mosaic are organized in clusters.To formally define them, consider the graph whose nodes are the cells and whose arcs are the shared facets (i.e. the 1-skeleton of the order-k Voronoi tessellation).A cluster is a connected component in the subgraph induced by the firstgeneration cells.It is not difficult to see that two such cells belong to a common cluster if and only if the corresponding rhomboids have the same anchor vertex.Let be one of these rhomboids and recall that the anchor vertex is A in ( ), which in this case is a collection of k − 1 points of A. Each combinatorial vertex of any cell in the cluster contains these k − 1 points, plus one additional point, which differentiates between these vertices.In other words, the cluster as a subcomplex of the order-1 Delaunay mosaic of these additional points.With this insight, we could replace the weighted Delaunay mosaic of the entire vertex set by multiple instances of unweighted Delaunay mosaics, namely one per cluster.This alternative strategy avoids the need to compute averages of points at the cost of extra book-keeping to group the vertex set of Del k (A) into clusters.We mention that in R 2 , the structure of each cluster satisfies the requirements that allow for the construction in time linear in the number of points [2].
Exact arithmetic.The CGAL software library [22] supports exact arithmetic by distinguishing between exact constructions and exact predicates.The latter are geometric tests with a true or false answer, such as whether or not a given point lies on a given sphere.By itself, the CGAL algorithm for weighted Delaunay triangulations requires exact predicates but no exact constructions.Our algorithm, on the other hand, computes averages of collections of input points, which are the locations of the vertices of the mosaic.This is an exact construction and indeed the only one needed to run our algorithm with exact arithmetic.In practice, exact constructions are a significant overhead with noticeable impact on the runtime, which would be nice to avoid.

Figure 1
Figure 1Superposition of the order-2 Voronoi tessellation (in black) and the order-2 Delaunay mosaic (in blue) of a set of six points in the plane.Each domain of the tessellation corresponds to two of these six points, and the corresponding vertex of the mosaic is the average of these two points.

Figure 3
Figure 3 Slices of a 4-dimensional rhomboid defined by Ain( ) = ∅ and Aon( ) = {a, b, c, d}.The non-trivial slices are a tetrahedron at generation g = 1, an octahedron at generation g = 2, and another tetrahedron at generation g = 3.

Figure 4
Figure 4 First-, second-, and third-order Delaunay mosaics of the set A = {a, b, c, d, e} in R 2 as slices of the 3-dimensional rhomboid tiling.For clarity, only two of the rhomboids are shown, with their first-generation slices in red and second-generation slices in dark blue.The rhomboids on the left and right are defined by Ain = {b}, Aon = {a, c, d} and Ain = ∅, Aon = {c, d, e}, respectively.

Lemma 3 .
A d-simplex, τ , in a triangulation of Del k (A) is a first-generation d-cell of Del k (A) if and only if the intersection of its combinatorial vertices is of size k − 1.

aFigure 5
Figure 5 Number of vertices (left) and of 3-dimensional cells (right) in the order-k Delaunay mosaics of four sets with n = 200 points in R 3 each.

Figure 6
Figure 6 Number of cells in the order-k Delaunay mosaics for small k in relation to the input size, for various 3-dimensional point sets.

Figure 7 Figure 8
Figure 7Variance of the number of 3-dimensional cells in the order-k Delaunay mosaics of randomly sampled points in convex position (left) and in a unit ball (right).The statistics of each plot are obtained from 30 sets of 50 points each.In black: the mean; in red: the range of one standard deviation around the mean; in grey: the range between the minimum and maximum.

Figure 9 Figure 10
Figure 9 Number of d-cells in the order-k Delaunay mosaic for 20 points (left) and 50 points (right) randomly sampled in the unit ball for different dimensions d.

Figure 11
Figure 11 From left to right: distribution of cluster sizes in Delaunay mosaics of order 2, 50, and 90 for 100 random points in the unit ball.40 41