Georg Osang
email: georg.osang@ist.ac.at

Mael Rouxel-Labbé
email: mael.rouxel.labbe@geometryfactory.com

Monique Teillaud
email: monique.teillaud@inria.fr

Mael Rouxel

Generalizing CGAL Periodic Delaunay Triangulations

Keywords: 2012 ACM Subject Classification Theory of computation → Computational geometry Delaunay triangulation, lattice, algorithm, software, experiments Digital Object Identifier 10.4230/LIPIcs.ESA.2020.37

Even though Delaunay originally introduced his famous triangulations in the case of infinite point sets with translational periodicity, a software that computes such triangulations in the general case is not yet available, to the best of our knowledge. Combining and generalizing previous work, we present a practical algorithm for computing such triangulations. The algorithm has been implemented and experiments show that its performance is as good as the one of the CGAL package, which is restricted to cubic periodicity.

Introduction

Delaunay triangulations are one of the most prominent structures in computational geometry. While nowadays many applications use Delaunay triangulations of finite point sets in Euclidean space, Delaunay originally introduced the notion in the context of infinite point sets with translational periodicity [START_REF] Delaunay | Sur la sphère vide. À la mémoire de Georges Voronoï[END_REF]. Such periodic point sets are abundant in fields such as crystallography and material sciences; thus their communities would benefit from software that computes Delaunay triangulations of infinite periodic point sets in Euclidean space. Specifically, given a d-dimensional lattice and its associated translation group, the orbits of a given finite point set in R d with respect to this translation group define a periodic point set.

Our aim is to compute a finite representation of the periodic Delaunay triangulation of such a periodic point set, specifically a projection of the triangulation onto the flat d-torus that is the quotient space of R d under the action of the translation group.

The first algorithm for this problem was already outlined in 1997 [START_REF] Dolbilin | Periodic Delone tilings[END_REF], yet a robust and efficient implementation for this problem does not exist to date as far as we know. The Voro++ library [START_REF] Rycroft | Voro++: A three-dimensional Voronoi cell library in c++[END_REF] is focused on crystallographic applications in 3 dimensions; however it is limited to orthogonal lattices. Zeo++ [START_REF] Thomas F Willems | Algorithms and tools for high-throughput geometry-based analysis of crystalline porous materials[END_REF] extends its functionality to arbitrary 3-dimensional lattices. Both libraries compute the Voronoi cell of a given input point as an intersection of half-spaces, by searching for other points around it that have an influence on its Voronoi cell. The combinatorics of the Delaunay triangulation cannot be easily accessed. The CGAL library [START_REF]The CGAL Project[END_REF] provides packages for periodic Euclidean Delaunay triangulations in 2D and 3D, which currently are limited to the integer lattice, referred to as the square and cubic setting and 2 and 3 dimensions, respectively [START_REF] Kruithof | 2D periodic triangulations[END_REF]6,5]. We propose an addition to CGAL that extends this functionality to arbitrary lattices.

The algorithm by Dolbilin and Huson [START_REF] Dolbilin | Periodic Delone tilings[END_REF] creates 3 d copies of each input point and computes their finite Delaunay triangulation, from which a representation of the periodic Delaunay triangulation is extracted. The CGAL algorithm [4,7] computes the triangulation in a finitely-sheeted covering space of the d-torus. It is based on the classical incremental algorithm by Bowyer and Watson [3,[START_REF] David F Watson | Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes[END_REF], and requires that the triangulation be a simplicial complex at any given time. Let us quickly recall that a triangulation is a simplicial complex, or is simplicial for short, if each of its simplices consists of a set of distinct vertices, and the intersection of any two simplices is either empty or a simplex. Operating directly on the d-torus does not guarantee this, see Figure 1a. Thus, in the cubic setting, a 3 d -sheeted cover is used (Figure 1b) until sufficiently many points have been inserted to guarantee that the triangulation in the 1-sheeted cover is a simplicial complex. Unfortunately, a 3 d -sheeted cover is not sufficient for more general periodic point sets: As the 3 d copies of each point have to be inserted iteratively into the 3 d -sheeted cover, simpliciality can be violated in the intermediate stages of point insertion, see Figure 1c. While there always exists a finitely-sheeted covering space [7] that ensures simpliciality, the number of sheets might be prohibitively large. Thus, we propose a different approach.

(a) The intersection of the two red edges is not a simplex but a set of two vertices. Thus the triangulation is not simplicial.

(b) The 9-sheeted cover guarantees a simplicial triangulation in the square setting.

(c) In the non-square setting, incremental point insertion of the 9 copies into the 9-sheeted cover can violate simpliciality. Here, the first of 9 copies is being inserted. Overview After formally defining the problem in Section 2, we propose an algorithm (Section 3) for periodic Delaunay triangulations that combines two different approaches and consists of two phases, both of which use Bowyer-Watson's algorithm. While for 2dimensional spaces algorithms based on flips circumvent the simpliciality requirement [START_REF] Despré | Flipping geometric triangulations on hyperbolic surfaces[END_REF], we stick to Bowyer-Watson's algorithm as it easily generalizes to 3 (and higher) dimensions. Furthermore it enables an efficient, clean, and easily maintainable implementation. The first phase of our algorithm (Section 3.2) refines the algorithm by Dolbilin and Huson [START_REF] Dolbilin | Periodic Delone tilings[END_REF], and its implementation details are based on some new results. It uses 3 d copies of each input point, regardless of the lattice, and computes a finite Euclidean Delaunay triangulation on this point set, from which a representation of the Delaunay triangulation on the d-torus is obtained. Once a simpliciality criterion is met, our algorithm switches to the second phase (Section 3.4), which conceptually follows the CGAL implementation of the cubic case [4]. It operates directly on the d-torus, maintaining only one copy of each input point, and thus provides better insertion running times than phase 1. A first version of our open-source implementation in 2D and 3D is publicly available.1 This implementation is expected to be an integral part of CGAL in a near future release. Experiments (Section 4) show similar performances as the CGAL implementation restricted to cubic lattices [5]. We close with a discussion of future extensions in Section 5.

Preliminaries

Let us recall various notions [START_REF] Engel | Geometric crystallography: an axiomatic introduction to crystallography[END_REF][START_REF] Conway | Low-dimensional lattices VI: Voronoi reduction of three-dimensional lattices[END_REF] that are employed throughout the algorithm. Let

B = {b 1 , b 2 , . . . , b d } be a basis of R d . The point set Λ := { d i=1 z i b i : z ∈ Z d } is called a lattice,
and B is its lattice basis. The lattice Λ is associated with the translation group Γ consisting of the translations φ λ : R d → R d mapping the origin to λ, for each λ ∈ Λ. The group Γ acts on R d and each translation of Γ maps Λ onto itself. We denote the length of the shortest non-zero lattice vector as sv(Λ). For a given lattice basis B, we call

B sup := B ∪ {b 0 } with b 0 = - d i=1 b i its superbase. A superbase is obtuse if for any pair b i and b j , b i , b j ≤ 0.
A basis is reduced if its superbase is obtuse [START_REF] Engel | Geometric crystallography: an axiomatic introduction to crystallography[END_REF]Definition 4.4]. This notion is defined in such a way that we can easily compute sv(Λ) and Dirichlet domains. The Dirichlet domain of a lattice point λ ∈ Λ is the region of λ in the Voronoi tesselation of Λ, or more formally dom(λ, Λ) := {p ∈ R d : p -λ ≤ p -ν ∀ν ∈ Λ}. It is a convex polytope, and we call the lattice Λ generic if each vertex of dom(0, Λ) is incident to the Dirichlet domains of exactly d other lattice points. For 2-dimensional generic lattices the Dirichlet domains are hexagons, for 3-dimensional generic lattices they are combinatorially equivalent to truncated cubes.

For a lattice vector λ, str(λ) = {p ∈ R d : -0.5 ≤ p,λ λ,λ < 0.5} is an infinite half-open strip that contains the subspace orthogonal to λ through the origin. Then dom(0, Λ) is the closure of the intersection of these strips for all non-zero lattice vectors. However, as dom(0, Λ) only has a finite number of facets, there must be a finite subset of strips whose closed intersection yields dom(0, Λ). Let V be the minimal set of lattice vectors (together with their negates) such that the closure of v∈V str(v) is dom(0, Λ). The vectors in V are commonly called Voronoi-relevant vectors. Each of them is a normal vector of a facet of the Dirichlet domain of 0, and thus there are at most 2(2 d -1) Voronoi relevant vectors [START_REF] Engel | Geometric crystallography: an axiomatic introduction to crystallography[END_REF]Theorem 3.6]. Let V + V -be a partition of V such that if v ∈ V + , then -v ∈ V -, and vice versa. For a fixed choice of V + (and implicitly V -), we define the canonical domain dom (0, Λ) := v∈V + str(v) (Figure 2). Its closure is dom(0, Λ), and its images under Γ, denoted dom (λ, Λ) := φ λ (dom (0, Λ)) for λ ∈ Λ, form a partition of R d . With kΛ for k ∈ Z referring to the lattice with basis {kb 1 , . . . , kb d }, we note that dom (0, kΛ) is dom (0, Λ) scaled by a factor of k.

For d ≤ 3, if we have a reduced lattice basis B with its superbase B sup , then V is a subset of { v∈S v : S ⊂ B sup , S = ∅, S = B sup } [8, Theorems 3 and 8], with equality if the lattice is generic. Note that sv(Λ) can be obtained as the length of the shortest vector in V.

Delaunay triangulations

The Delaunay triangulation Del(X) of an input point set X ⊂ R d is a collection of simplices up to dimension d whose vertex set is X and each d-simplex (which we refer to as a cell) corresponds to a set of d + 1 points whose open circumscribing ball does not contain any other points of X. We call the (d -1)-simplices of Del(X) its facets.

Given a lattice Λ and a finite set of points X, we get the periodic point set ΓX := {φ λ (x) : x ∈ X and λ ∈ Λ} consisting of the elements of the orbits of X under Γ. ΓX is globally invariant under Γ. Then Del(ΓX) is the periodic Delaunay triangulation of the infinite point set ΓX. Note that we can ignore degeneracies in ΓX by using the symbolic perturbation provided by CGAL [START_REF] Devillers | Perturbations for Delaunay and weighted Delaunay 3D triangulations[END_REF]: it is translation-invariant, so, degeneracies are triangulated in a consistent way, which ensures that the computed Del(ΓX) is actually invariant under Γ. The orbit space R d /Γ is a flat torus, and we denote its projection map as π : R d → R d /Γ. The torus triangulation Del(ΓX)/Γ is the projection of Del(ΓX) into R d /Γ. Using dom (0, Λ) as a geometric representation of the torus, we can use X 0 := ΓX ∩ dom (0, Λ) as canonical representatives of the vertex set of Del(ΓX)/Γ. While Del(ΓX)/Γ gives us a finite representation of Del(ΓX), unlike Del(ΓX) it is not necessarily simplicial (see Figure 1a).

Algorithm

The input to our algorithm is a lattice basis B for Λ and a set of points X defining the periodic point set. The output is an object representing Del(ΓX)/Γ. This object provides a uniform interface that, regardless of the internal state of our algorithm, allows the user to access the properties of Λ as well as the torus triangulation Del(ΓX)/Γ. For Λ this includes the reduced basis B = {b 1 , . . . , b d }. For Del(ΓX)/Γ this includes the set X 0 of canonical representatives for its vertex set, and its cells. Cells are not solely defined by their vertex set, but have additional geometric information attached. Specifically, each vertex of a cell is represented as a point x from X 0 with an associated offset o = (o 1 , . . . , o d ), which is an integer vector. The geometric location of the vertex then is x

+ d i=1 o i b i .
In alignment with other CGAL triangulations, we also provide access to simplices of lower dimensions, represented with associated vertex offsets akin to cells, as well as neighborhood relations. These include querying a cell for its adjacent cells, or querying a vertex for its incident cells or lower dimensional simplices. While many steps generalize, we restrict our focus to 2-and 3-dimensional triangulations, which are the most widely used.

Internally, our algorithm operates in two phases, which use two different data structures, respectively: The first phase maintains a finite Euclidean Delaunay triangulation while the second phase maintains a triangulation of Del(ΓX)/Γ. Let us now recall the result that is crucial to the first phase of our algorithm. For a point set X, let X 3 := dom (0, 3Λ) ∩ ΓX, i.e. all periodic copies of X that lie within the Dirichlet domain of 0 scaled by a factor of 3. We call a cell of Del(X 3 ) a periodic cell if it is also a cell of the periodic triangulation Del(ΓX) (see Figure 3a).

Proposition 1 [START_REF] Dolbilin | Periodic Delone tilings[END_REF]Lemma 3.4]). Given a point set X, each cell of Del(X 3 ) that has at least one vertex in dom (0, Λ) is a periodic cell. Furthermore the set of these cells contains at least one periodic copy of each cell of Del(ΓX)/Γ.

After some preprocessing that essentially consists in computing the canonical domain, the first phase internally maintains Del(X 3 ) using the CGAL packages for Euclidean Delaunay triangulations [START_REF] Jamin | 3D triangulations[END_REF][START_REF] Yvinec | 2D triangulation[END_REF]. For this we need to develop a systematic way of computing X 3 and obtaining the interface for Del(ΓX)/Γ from the internal data structure Del(X 3 ).

(a) Cells in blue are guaranteed to be periodic cells by Proposition 1. Green: non-periodic cells. Black vertices: Λ.

0 0 0 1 1 1 1 0 -1 -1 -1 1 -1 0 0 -1 1 -1 2 1 1 2 -1 -2 -2 -1 (b)
The set of canonical cells in blue. Each copy of dom(0, Λ) is labeled with its offset.

Figure 3

In blue, the points X3 with their Delaunay triangulation Del(X3). The large hexagon is dom(0, 3Λ), while the smaller ones are dom(0, Λ) and (middle, right) its periodic copies. The reduced lattice basis is drawn in red (left, middle).

Once we can guarantee that Del(ΓX)/Γ is simplicial and will remain so for any future point insertions, we switch to the second phase. Simpliciality guarantees that the Bowyer-Watson algorithm can be used directly on Del(ΓX)/Γ. For this purpose we leverage the CGAL machinery for periodic triangulations from the cubic setting [START_REF] Kruithof | 2D periodic triangulations[END_REF]5] and enhance its underlying data structures to work for the generic setting. As we keep only one representative for each vertex, inserting points in this phase is more efficient than in the first phase.

The remainder of this section describes the two phases and the transition between them, as well as the preprocessing of the lattice.

Lattice preprocessing

We first compute the reduced lattice basis B, which allows us to compute V, the face normals of the canonical domain. We use V to represent dom (0, Λ) and check for containment of a point within dom (0, Λ), which helps obtaining the canonical copy of each input point. To find all periodic copies of a canonical point x that lie within dom (0, 3Λ), it is sufficient to find all points λ such that the domain dom (λ, Λ) intersects dom (0, 3Λ), and then check E S A 2 0 2 0 37:6 Generalizing CGAL Periodic Delaunay Triangulations φ λ (x) for containment within dom (0, 3Λ) for these points. As the set of these points λ is independent of x, we compute it before phase 1.

Lattice reduction

Many lattice related problems, such as the shortest non-zero vector problem (SVP), are believed to be hard in general [START_REF] Van | Another NP-complete problem and the complexity of computing short vectors in a lattice[END_REF]2]. However in low dimensions, we can use two classical iterative algorithms to solve lattice reduction. Let B be a lattice basis and B sup = {b 0 , . . . , b d } its superbase, as defined in Section 2. Define c ij := b i , b j .

2-dimensional reduction

A 2-dimensional lattice basis is Lagrange-reduced [13, Section 4.2] if 0 ≤ 2|c 12 | ≤ c 11 ≤ c 22 . We can negate b 2 if necessary so that c 12 ≤ 0. Then {b 1 , b 2 , b 0 := -b 1 -b 2 } form an obtuse superbase, because c 01 = b 1 , -b 1 -b 2 = -c 11 -c 12 ≤ 0 due to 2|c 12 | ≤ c 11 ,
b 2 := b 2 -sb 1 , with s = 1 if c 12 > 0 and s = -1 otherwise. Vector b 2 is shorter than b 2 since c 22 = c 22 + c 11 -2sc 12 .
Since there are only finitely many pairs of lattice vectors within a ball of any given radius, a finite number of applications of this procedure will yield a Lagrange-reduced basis.

3-dimensional reduction

We outline Selling's algorithm [START_REF] Engel | Geometric crystallography: an axiomatic introduction to crystallography[END_REF]Section 4.4], an iterative algorithm that obtains an obtuse superbase in 3 dimensions. In each step of the algorithm, if there is a c ij > 0, the algorithm returns a new superbase B sup defined via b

i := -b i , b j := b j , b h := b h + b i and b k := b k + b i
where h and k are the remaining two indices different from i and j. For any basis B, let σ(B) := b∈Bsup ||b|| 2 . Notice that in each step, σ(B ) = σ(B ) -2c ij . In particular, σ(B ) < σ(B ). This fact, together with the fact that there are only finitely many lattice points in a ball of radius σ(B ) (and thus only finitely many quadruplets of vectors whose square magnitudes sum up to at most σ(B )) guarantees termination of the algorithm.

Intersecting domains

For a given canonical point x ∈ X 0 , each of its periodic copies y ∈ Γx can be obtained as φ λ (x) for some translation φ λ ∈ Γ. The corresponding lattice point λ can be uniquely written as d i=1 o i b i for some o ∈ Z d and we write and λ(o) := λ. We then call o =: o(y) the offset of y and φ λ(o) the translation associated with o. For each point x ∈ X 0 we need to determine the 3 d offsets for which φ λ(o) (x) is within dom (0, 3Λ). Fortunately, these offsets have to come from a fixed set of offsets that only depends on the lattice basis. Notice that if φ λ(o) (dom (0, Λ)) does not intersect dom (0, 3Λ), then φ λ(o) (x) cannot be in dom (0, 3Λ) for x ∈ X 0 . Thus we only need to check those offsets for which φ λ(o) (dom (0, Λ)) ∩ dom (0, 3Λ) = ∅. Lemma 2. There are at most 4 d -2 d + 1 translates of dom (0, Λ) that have non-empty intersection with dom (0, 3Λ).

Proof. If dom (λ, Λ) intersects dom (0, 3Λ) for some λ ∈ Λ, then λ must be inside dom (0, 4Λ). There are 4 d lattice points within dom (0, 4Λ). For each v ∈ V + , the lattice vectors -2v and 2v are on the boundary of dom (0, 4Λ), however only -2v is one of those 4 d points within dom (0, 4Λ). For these lattice points -2v on the boundary however the intersection between dom (-2v, Λ) and dom (0, 3Λ) is empty. Thus only

4 d -|V + | = 4 d -(2 d -1) translates of dom (0, Λ) intersect dom (0, 3Λ).
In general, we can find this set of offsets via a breadth-first search: We define a graph on Λ with each lattice point λ connected to λ + v for v ∈ V. We start the search at lattice point 0 and terminate once we have found 4 d -2 d + 1 offsets, or in the case of non-generic lattices once all 4 d lattice points inside dom (0, 4Λ) have been reached (with the ones on the boundary being discarded).

In 2 dimensions, there is in fact a fixed set S of 13 offsets such that for any lattice, dom (λ(o), Λ) only intersects dom (0, 3Λ) when o ∈ S (Figure 3b). Lemma 3. Given a 2-dimensional lattice with basis B and an offset o, if o ∈ {(0, 0), (-1, -1), (0, 1), (1, 0), (-1, 0), (0, -1), (1, 1), (-1, -2), (1, 2), (-2, -1), (2, 1), (-1, 1), (1, -1)}, then dom (λ(o), Λ) ∩ dom (0, 3Λ) = ∅.

Proof. Each dom (λ, Λ) adjacent to dom (0, Λ) has all but one of its facets in the interior of dom (0, 3Λ). Except for those adjacent via one of those facets, all dom (λ, Λ) at graph distance 2 from 0 intersect dom (0, 3Λ). The ones that don't are exactly the 6 domains dom (2v, Λ) for v ∈ V. As there are 19 domains at distance at most 2, we have found 13 domains that intersect dom (0, 3Λ). As this is the maximum possible number by Lemma 2, we have found all of them.

In 3 dimensions, the same argument yields a set of 51 fixed offsets. However, by Lemma 2, the total number of domains that intersect dom (0, 3Λ) is 57. Fortunately, the remaining 6 offsets come from a fixed set of 24 offsets that is independent of the lattice. Claim. Domains dom (ν, Λ) of graph-distance at most (k -2) from some λ = kv for an integer k ≥ 2 and v ∈ V cannot intersect dom (0, 3Λ).

Proof. We have

ν = kv + v 1 + • • • + v k-2 for some v i ∈ V, i = 1, . . . , k -2. Then w, v = k v, v + v, v 1 + • • • + v, v k-2 ≥ k v, v -(k -2) v, v = 2 v, v
due to the observation above. This means that ν is not in the interior of str(4v) and thus the interior of the 4-scaled domain. Therefore dom (ν, Λ) ∩ dom (0, 3Λ) = ∅.

For a generic 3-dimensional lattice, the Voronoi relevant vectors correspond to the non-empty subsets of B sup of size at most 3, independent of the actual vectors of B sup . Thus the graph we defined on Λ is isomorphic for all generic lattices. For non-generic lattices it is isomorphic to a subgraph of the graph for generic lattices; thus it is sufficient to restrict our attention to generic lattices.

We enumerate all lattice points at graph distance 3 from the origin. For a generic lattice, this yields 110 points, however all but 24 of them can be written as 3v or 3v + w for Voronoi relevant vectors v and w. The set of the 24 corresponding offsets is O 3 , while O ≤2 are those 51 offsets that the argument from Lemma 3 yields. All lattice points at graph distance 4 can be written as 4v + kw for k ≤ 2 and Voronoi relevant vectors v and w, and thus their domains and consequently any domains at higher distances cannot intersect dom (0, 3Λ).

Phase 1

Phase 1 maintains the Euclidean Delaunay triangulation of X 3 . For each new point to be inserted, we first find its canonical copy. Then we can compute its periodic copies that are contained within dom (0, 3Λ) using Lemmas 3 and 4. These are then inserted into Del(X 3 ). To provide user access to the cells of Del(ΓX)/Γ, we define a notion of canonical cell in Del(X 3 ) to get a representative for each cell of Del(ΓX)/Γ.

Canonical points

For each point of X, finding its periodic copy that lies in dom (0, Λ) is equivalent to solving the closest vector problem (CVP), i.e. given x ∈ X, determining the lattice point that is closest to x. For arbitrary dimensions this problem is known to be NP-hard [START_REF] Van | Another NP-complete problem and the complexity of computing short vectors in a lattice[END_REF]. For the exact version of CVP, various iterative algorithms have been described [1, [START_REF] Sommer | Finding the closest lattice point by iterative slicing[END_REF][START_REF] Micciancio | A deterministic single exponential time algorithm for most lattice problems based on Voronoi cell computations[END_REF]. As we are only operating in 2 and 3 dimensions, any of them would suffice for us in practice in terms of running time, and we will describe the algorithm by Sommer et al [START_REF] Sommer | Finding the closest lattice point by iterative slicing[END_REF] due to its simplicity.

For a real number r, define round(r) to be the closest integer to r. If this integer is not unique, then it is the one with the smallest absolute value. This definition ensures convergence in cases where x is on the boundary of a Dirichlet domain of the lattice [START_REF] Sommer | Finding the closest lattice point by iterative slicing[END_REF]. Recall that V + is a set of normals of the facets of the canonical domain, and can be obtained from a reduced basis B (see Section 2). We first sort the vectors of V + by their magnitude. As dom (0, Λ) is the intersection of the strips str(v) for v ∈ V + , we need to find the periodic copy of x that is in all these str(v). We loop through the vectors v of V + and for each of them perform the following operations: Compute c := x,v v,v . If -0.5 ≤ c < 0.5, then already x ∈ str(v). If not, then we subtract round(c) • v from x. Note that after such a step, it is guaranteed that x ∈ str(v). However after modifying x for a longer vector v, it might happen that x is moved outside of str(v ) for some shorter vector v again. Therefore we need to repeatedly loop through V + and perform this operation until x ∈ str(v) for all v ∈ V + .

Notice that in each step where x is modified, the magnitude of x strictly decreases. Therefore this algorithm terminates in finite time, as the number of lattice vectors within a ball of given radius ||x|| is finite.

Extracting representative cells

Recall that our interface specifies access to the cells (and lower-dimensional simplices) of Del(ΓX)/Γ. Thus for each of its cells we have to provide one representative from Del(X 3 ). For a cell σ of either Del(ΓX) or Del(X 3 ), with V (σ) = {x 1 , . . . , x d }, we define its offset as the vector o(σ) = min x∈V (σ) {o(x)} where the minimum is taken lexicographically. Note that this definition differs from [4, Convention 3.3.1] where the coordinate-wise minimum is taken. If the offset of a cell is the 0-vector, we call it a canonical cell. Note that due to our definition a canonical cell always has a vertex with offset 0. Therefore by Proposition 1 a canonical cell of Del(X 3 ) is also periodic , see Figure 3b. For σ ∈ Del(ΓX), the translated cell σ -d i=1 o i b i is called its canonical representative. This means that for each class of cells in Del(ΓX) (or equivalently each cell of Del(ΓX)/Γ), there is a unique canonical representative in Del(X 3 ). Therefore we can get a set of representative cells by iterating over the cells of Del(X 3 ) and selecting those that are canonical.

Neighborhood relations

In accordance with our interface, we need to provide neighborhood relations for the vertices and cells of Del(ΓX)/Γ. These vertices and cells are represented by the canonical vertices and cells of Del(X 3 ), whose neighbors in Del(X 3 ) (which we have access to) may differ from the neighbors in Del(ΓX)/Γ (which we want to find). We outline how to get the neighors of a cell of Del(ΓX)/Γ from Del(X 3 ), and note that other neighborhood relations work in a conceptually similar way. Note that we do not store these relations explicity, but we compute them upon request and may cache them for future access.

Consider a canonical cell σ of Del(X 3 ) and a neighboring cell σ n . If σ n is canonical, it is the neighbor of σ in Del(ΓX). If σ n is not canonical but has a vertex in dom (0, Λ), then it is periodic and we return the canonical representative of this cell. However, it is possible that all vertices of σ n are outside dom (0, Λ), and thus σ n might not be a periodic cell at all. In that case, we need to consider the facet τ separating σ and its neighbor σ n . As τ is a facet of a canonical cell, it is also a periodic facet (i.e. a facet of Del(ΓX)). We compute τ 's offset o := o(τ ). Then τ := φ λ(-o) (τ ) is the canonical representative of τ and σ := φ λ(-o) (σ) as well as its neighbor σ n across τ are periodic cells because they share τ , which has a vertex in dom (0, Λ). If σ n is canonical, then is it the canonical representative of the neighbor of σ in Del(ΓX); if not, then its canonical representative is. Figure 3c illustrates this process.

In practice, we store the canonical representative of each vertex of X 3 . To obtain the canonical representative of a cell (or facet) σ, we need to choose a vertex x whose offset is minimal (lexicographically) among its vertices. This ensures that the periodic copy x of x in the canonical version of the cell is inside dom (0, Λ). Then one of the cells (or facets) incident to x in Del(X 3 ) is the canonical representative of σ.

Transition

Phase 2 is more efficient than phase 1 as it directly operates on Del(ΓX)/Γ; however, we cannot use it from the start as the Bowyer-Watson algorithm [3,[START_REF] David F Watson | Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes[END_REF] comes with some constraints. The Bowyer-Watson algorithm is an incremental algorithm inserting points one by one. For each new point x, it determines the conflict zone, which is the set of cells whose circumsphere contains x. All these cells are removed, and all boundary facets of the resulting hole are connected to x to fill in the hole with new cells. The algorithm requires this hole to be a topological d-ball, which is not always guaranteed for Del(ΓX)/Γ. However the following criterion is a sufficient condition for the Bowyer-Watson algorithm to work [7].

Lemma 5 ([7,Criterion 3.11]). If for every cell in Del(ΓX)/Γ the circumradius is smaller than 1 4 sv(Λ), then Del(ΓX )/Γ is simplicial for every X ⊇ X.

If this criterion is fulfilled, we can safely switch to phase 2. To detect at which point in phase 1 the criterion is fulfilled, we maintain a set S big of big cells, which are canonical cells whose circumradius is larger than or equal to 1 4 sv(Λ). We update S big during each point insertion.

Assume we wish to insert a new point x into the periodic triangulation of some point set X. On a high level, we need to remove the big cells in the conflict zone of x from S big , and then add the newly created big cells to S big . In practice, we have already computed the triangulation Del(X 3 ) and have to insert the periodic copies of x that are within dom (0, 3Λ), i.e. Γx ∩ dom (0, 3Λ), into Del(X 3 ). We first detect the conflict zone of x 0 , the canonical copy of x. For each big cell in the conflict zone we check if it has a canonical copy in S big , and remove that copy from S big if it does. While the conflict zone may contain non-periodic E S A 2 0 2 0 37:10 Generalizing CGAL Periodic Delaunay Triangulations cells, Lemma 6 guarantees that we capture a representative of each cell from Del(ΓX) that is in conflict with x. Next we insert all the copies Γx ∩ dom (0, 3Λ) into Del(X 3 ). Finally, for each big cell incident to x 0 in the resulting triangulation, we add its canonical copy to S big . Note that all these cells have a canonical copy by Proposition 1 as they have a vertex inside dom (0, Λ). Lemma 6. Let C be the conflict zone of some point x ∈ dom (0, Λ) with respect to Del(ΓX). Then the conflict zone of x with respect to Del(X 3 ) contains at least one periodic copy of each cell from C.

Proof. First observe that if a cell of Del(ΓX) has its circumcenter at the origin, then all its vertices must be within dom(0, Λ), because if the circumsphere contains a point outside dom(0, Λ), then it also contains its canonical copy. Via translation it follows that if a cell has its circumcenter in dom (0, 2Λ), then its vertices are in dom (0, 3Λ) and thus it is a cell of Del(X 3 ). So assume we have a cell σ in C whose circumcenter c is not within dom (0, 2Λ). Then there is a facet of dom (0, 2Λ) with respect to which c is outside. Let f be its face normal. Then x + f is also contained in the circumsphere of σ. Reversely, σ -f is a cell whose circumsphere contains x, and furthermore its circumcenter is closer to 0 than c, see Figure 4. As there are only finitely many periodic copies of c within a given distance from

c c -f 0 f Figure 4
The red point is in the conflict zone of the blue cell σ, which is not a cell of Del(X3) and whose circumcenter c is outside dom (0, 2Λ). However the light blue cell σ -f is a periodic copy of σ that is in Del(X3) and x is in its conflict zone. 0, after applying this process a finite number of times we eventually obtain a cell whose circumsphere contains x and whose circumcenter is in dom (0, 2Λ). This cell is a periodic copy of σ, is contained in Del(X 3 ) and thus also part of the conflict zone of x with respect to Del(X 3 ).

Once S big is empty, the criterion of Lemma 5 is fulfilled, and we can internally convert our triangulation from Del(X 3 ) to Del(ΓX)/Γ, which is maintained in phase 2. We initialize the periodic triangulation data structure with the set of canonical vertices X 0 and canonical cells obtained from Del(X 3 ), as well as the adjacency and incidence relations outlined earlier.

Phase 2

Phase 2 operates directly on the torus triangulation Del(ΓX)/Γ. Thus it maintains only one copy of each cell and vertex and point insertion is faster than in phase 1. The data structure it uses to represent Del(ΓX)/Γ is akin to the one used in CGAL for the cubic case, and closely resembles the interface we defined for our algorithm: Only X 0 is stored as vertex set, and each cell is represented by its vertices, with a vertex encoded as a pair (x, o) of a point x ∈ X 0 and an offset so that its geometric location is φ λ(o) (x). While in the cubic case each offset coordinate either takes the value 0 or 1 and offsets can be encoded in d bits, in our more general setting offsets can take any of the lattice-specific values from Lemma 2. Unlike in phase 1, most neighborhood relations for Del(ΓX)/Γ are already stored explicitly in the data structure: For each cell its adjacent cells are stored, and for each vertex one incident cell is stored. The remaining neighborhood relations required by our interface are obtained implicitly from the explicitly stored ones. With each point insertion, all stored neighborhood relations have to be updated accordingly.

Point insertion

To insert a new point x into Del(ΓX)/Γ using the Bowyer-Watson algorithm, first we compute the canonical copy x 0 . Then we locate the cell containing x 0 , via a traversal starting from an arbitrary cell. The conflict zone is computed via a search starting in the cell containing the point 0 . Whenever we are traversing cells, care has to be taken to maintain the correct offset of the affected cells relative to x 0 , and similarly when creating new cells to fill in the hole left by the deleted conflict zone. For the cubic case the details are described in [4, Section 3.3], and we omit the technical adjustments needed to make these steps work in the more general case.

4

Experimental results

Points until transition

We experimentally evaluated the number of points required until the criterion of Lemma 5 is fulfilled and the transition to phase 2 occurs. Lemma 5 requires all of R d /Γ to be covered by the balls of radius 1 4 sv(Λ) around X 0 . As the volume of the torus equals the volume of the Dirichlet domain of 0, denoted as vol(dom(0, Λ)), we expect the number of points until switching to phase 2 to be roughly proportional to vol(dom(0, Λ))/sv(Λ) d , assuming the points are sampled uniformly at random. To investigate this in 2 dimensions, we parametrize a 2-dimensional space of lattices. We call the parameters the elongation and the skew s. The basis of the lattice with elongation and skew s is b 1 = ( , 0) and b 2 = (s • /2, 1). With ≥ 1 and s between 0 and 1 we can parametrize all 2-dimensional lattices up to symmetry and scaling. Note that the skew affects sv(Λ) but not vol(dom(0, Λ)), while the elongation is proportional to vol(dom(0, Λ)) but does not affect the sv(Λ). Fixing the skew at 0 and varying the elongation (Figure 5a), we see that the number of random points needed until the phase switch appears to be proportional to the elongation. The same applies to the number of points until the resulting triangulation is simplicial for the first time. Figure 5b shows the same statistics for lattices of fixed area but varying skew. In addition we plot the inverse of sv(Λ)

2 for comparison, and observe that it behaves similarly albeit not entirely proportionally. In 3-dimensions, a parametrization of lattices up to symmetry and scaling needs 5 parameters, and thus an analysis like in 2D is impractical. As mentioned before, we expect that the number of points until the switch to phase 2 is proportional to vol(dom(0, Λ))/sv(Λ) 3 .

By comparing the two, Figure 6 confirms this relationship for 182 lattices whose basis vectors have randomly sampled direction and a random magnitude between 1 and 100.

Running times

We evaluate the running times of our algorithm for different 3-dimensional lattices, and compare them to the existing CGAL implementations. The data points collected are from Delaunay triangulations of 10 k uniformly sampled random points for k up to 7. Each data point is an average of 300 trials (10 trials for 10 7 points). The experiments were conducted on a laptop running Fedora 30 64-bits, with two 6-core Intel(R) i9-8950HK CPU clocked at 2.90GHz, and with 32GB of RAM. The CGAL kernel used was CGAL::Exact_predicates_inexact_constructions_kernel and CPU time was measured using CGAL's timer tools. The code was compiled using clang 8.0.0 with compilation flags -O3 and -DNDEBUG.

Table 1 shows a comparison between the CGAL implementation of Euclidean Delaunay triangulations [START_REF] Jamin | 3D triangulations[END_REF] (with random points uniformly sampled in the unit cube), periodic Delaunay triangulations in the cubic setting [5] and our algorithm for various lattices, including the cubic lattice. For each lattice, we also measured the average number of points until the switch to phase 2. As our algorithm and the one for the cubic setting are based on the same code base when operating directly on the torus triangulation, their runtimes are comparable for large point sets. It should be noted that when similar experiments were conducted in 2010 to compare the Euclidean and cubic periodic algorithms [4, Section 3.6.2], both were performing comparably. Since then, Euclidean Delaunay triangulations in CGAL have seen significant optimizations, which were not applied to periodic triangulations. This also explains why our algorithm is faster than the cubic periodic algorithm for point sets where phase 1 takes up a significant portion of the running time, as internally we use a Euclidean rather than a periodic triangulation.

Discussion

Weighted points Phase 1 of our algorithm readily generalizes to weighted point sets. In particular, Proposition 1 still holds for weighted point sets (see Appendix A for a proof). Phase 2 only works for weighted points under additional restrictions because Del(ΓX)/Γ can not be guaranteed to remain simplicial, in particular after inserting points with large weights. Such point insertions that break simpliciality can be prevented by requiring points to be inserted in decreasing order of weight, or like in the cubic implementation in CGAL by severely restricting the range of weights a point can have. Thus an implementation is subject to a tradeoff between flexibility (phase 1) and performance (phase 2).

Higher dimensions Conceptually, our algorithm generalizes to higher dimensions. The only step that is dimension-dependent and therefore requires significant adjustments is obtaining a representation of dom (0, Λ). While up to 3 dimensions every lattice has an obtuse superbase, this does not hold in higher dimensions. Therefore we would need to find a different way of obtaining the set of Voronoi relevant vectors. Complexity-wise, computing the canonical representative of a point is believed to be hard. Furthermore, the number of copies that we compute of each input point in phase 1 is 3 d . Therefore we expect the runtime of our algorithm to be exponential in the dimension, but still feasible in practice as long as the dimension is not too large.

Dummy points

For best performance, transition from phase 1 to phase 2 should occur as soon as possible. Intuitively, this happens when the input point distribution does not have large gaps. We can achieve this by first inserting additional dummy points, which are removed at the end (if possible) [7]. In practice, we can choose these points from a sufficiently fine hexagonal lattice (or body-centered cubic in 3 dimensions), such that the open spheres of radius 1 4 sv(Λ) around the points cover the entire torus.

Software distribution

We aim to provide a CGAL package for periodic Delaunay triangulations for arbitrary lattices in 2 and 3 dimensions. The current state of our implementation is available online (see Footnote 1). Both the 2-and 3-dimensional implementations have been integrated into the CGAL codebase and only require some additional refactoring and optimizing to adhere to CGAL's quality standards. Automated tests and documentation still have to be produced. An extension to weighted Delaunay triangulations, referred to as regular triangulations in CGAL, is planned for the future. 1 Running time (in seconds) of various Delaunay triangulation algorithms on random point sets of different sizes. Our algorithm ("Lattice") is evaluated for different lattices including the cubic lattice, face-centered cubic (FCC) lattice and two other lattices Λ1 and Λ2 with bases B1 = {(0.5, -0.5, 0.1), (-0.5, 0.5, 0.1), (0.5, 0.5, -0.1)} and B2 = {(1, 0, 0), (-0.5, √ 3/2, 0), (0, 0, 0.05)}. For each lattice we also record the average number of points until the switch to phase 2 (n switch ).

A Generalization to weighted points

For a weighted point x, we denote its weight with w x . For weighted points x, y ∈ X, their power distance is dist(x, y) := x -y -w x -w y . For an arbitrary point p ∈ R that is not part of X we define the power distance from x to p as dist(x, p) := x -p -w x (as if the weight of p was 0). We call two points x and y orthogonal if dist(x, y) = 0. For d + 1 points Q ⊆ X (if in general position) there is a unique point z with weight w z that is orthogonal to all points of Q. This point is called the orthocenter or power sphere of Q. If Q is a set of points forming a cell (simplex) in the Delaunay triangulation, then for all points x ∈ Q and y ∈ X \ Q it holds that dist(y, z) > dist(x, z) = 0. We call this the empty-sphere property. The converse holds as well. The perpendicular bisector of two weighted points x and y is the set of points p with dist(x, p) = dist(y, p).

Proposition 7 (Generalization of Lemma 3.2-3.4 from [START_REF] Dolbilin | Periodic Delone tilings[END_REF]). Given a set X of representatives for our point set, let X 3 := dom (0, 3Λ) ∩ ΓX, i.e. all periodic copies of these points that lie within the Dirichlet domain of 0 scaled by a factor of 3. Let T 0 be those cells of Del(X 3 ) that have at least one vertex in dom (0, Λ). Then the cells of T 0 are all part of the triangulation of ΓX. Furthermore, these cells contain at least one representative of each class of cells from Del(ΓX).

Proof. We will prove the proposition in 3 steps.

Claim 1. Assume one of our points is x 0 = 0 (with arbitrary weight). Then the orthocenter c T of any Delaunay cell that has x 0 as a vertex is within dom(0, Λ).

Proof. Assume not. Then there is a face F i of dom(0, Λ) such that F i separates x 0 from c T . Let f i be the corresponding translation lattice vector orthogonal to F i . Then x 0 + f i is strictly closer to c T than x 0 because it has the same weight as x 0 . This is a contradiction to the empty-sphere property.

Claim 2. For a cell containing x 0 as a vertex, all vertices are in dom(0, 2Λ).

Proof. Assume some vertex x is not. Consider the perpendicular bisector between x and x -f i where f i is the face normal vector orthogonal to the face F i of dom(0, Λ) with respect to which x is outside of dom(0, 2Λ). This bisector is separated from x 0 and the orthocenter c T by F i , the face of dom(0, Λ) that is parallel to this bisector. In particular, it follows from Claim 1 that c T is (strictly) on the x -f i side of the perpendicular bisector. So x -f i is closer to c T than x and thus the empty-sphere property is violated. See Figure 7 for reference. Claim 3. All cells σ having a vertex in dom (0, Λ) are entirely contained in dom (0, 3Λ).

Proof. Let x be a vertex of σ that is within dom (0, Λ). Shift the entire point set ΓX and its triangulation by the vector -x so that x now coincides with 0. Now from Claim 2 it follows that the other vertices of the shifted cell are within dom(0, 2Λ). Adding the vector +x to these shifted vertices we get back to the original setting, but because x ∈ dom (0, Λ) we also know now that these vertices are within dom (0, 3Λ), as dom (0, 3Λ) is the Minkowski sum of dom (0, Λ) and dom(0, 2Λ).

f i F i x 0 x x -f i c T Figure 7
A cell (dashed) with a vertex x that is outside of dom(0, 2Λ). The dotted line is the perpendicular bisector between x and x -fi.

Every cell of Del(ΓX) that has a vertex in dom (0, Λ) has all its vertices in X 3 . Furthermore, because it fulfils the empty-sphere property in ΓX, then it also fulfils this property in X 3 ⊂ ΓX, and thus is present in Del(X 3 ). As every point from X has a representative in dom (0, Λ), also each cell from Del(ΓX) has a representative in T 0 , proving the statement.

B Detailed proof of Lemma 4

Lemma 4 (extended version). Let O 3 := {(3, 2, 1), (2, 1, -1), (3, 1, 2), (2, -1, 1), (1, -1, -2), (1, -2, -1), (2, 3, 1), (1, 2, -1), (1, 3, 2), (-1, 2, 1), (-1, 1, -2), (-2, 1, -1), (2, 1, 3), (1, -1, 2), (1, 2, 3), (-1, 1, 2), (-1, -2, 1), (-2, -1, 1), (-1, -2, -3), (-1, -3, -2), (-2, -1, -3), (-3, -1, -2), (-2, -3, -1), (-3, -2, -1)} and O ≤2 := {(0, 0, 0), (-1, -1, -1), (0, 0, 1), (0, 1, 0), (1, 0, 0), (1, 1, 0), (1, 0, 1), (0, -1, -1), (0, 1, 1), (-1, 0, -1), (-1, -1, 0), (1, 1, 1), (0, 0, -1), (0, -1, 0), (-1, 0, 0), (1, 2, 0), (1, 0, 2), (-1, -2, -2), (2, 1, 0), (2, 0, 1), (1, -1, -1), (0, 1, 2), (-2, -1, -2), (0, 2, 1), (-1, 1, -1), (-2, -2, -1), (-1, -1, 1), (1, 1, 2), (-1, -1, -2), (1, 2, 1), (0, 1, -1), (-1, -2, -1), (0, -1, 1), (2, 1, 1), (1, 0, -1), (1, -1, 0), (-2, -1, -1), (-1, 0, 1), (-1, 1, 0), (1, 2, 2), (-1, 0, -2), (-1, -2, 0), (2, 1, 2), (0, -1, -2), (2, 2, 1), (1, 1, -1), (0, -2, -1), (1, -1, 1), (-2, -1, 0), (-2, 0, -1), (-1, 1, 1)}.

Given a 3-dimensional lattice and an offset o, there is a subset S of size 6 of O 3 such that if o ∈ O ≤2 ∪ S, then dom (λ(o), Λ) ∩ dom (0, 3Λ) = ∅.

Proof. As mentioned in the proof of Theorem 3.5 from [START_REF] Engel | Geometric crystallography: an axiomatic introduction to crystallography[END_REF], v 2 is on the boundary of dom (0, Λ) for every Voronoi-relevant vector v.

Claim 1. The orthogonal projection of v into the 1-dimensional subspace spanned by another Voronoi relevant vector w ∈ V + has magnitude less than w.

Proof. All facets of dom (0, Λ) are contained in str(w) for all w ∈ V + , so because v 2 is on the boundary of dom (0, Λ), it is contained in str(w). By definition of str(w) this implies | 1 2 v, w / w, w | ≤ 1 2 , from which it follows that | v, w | ≤ w, w .

Claim 2. Domains dom (ν, Λ) of graph-distance (k -2) from some λ = kv for an integer k and v ∈ V cannot intersect dom (0, 3Λ).

Proof. We have

ν = kv + v 1 + • • • + v k-2 for some v i ∈ V. Then w, v = k v, v + v, v 1 + • • • + v, v k-2 ≥ k v, v -(k -2) v, v = 2 v
, v due to Claim 1. This means that ν

Figure 1

 1 Figure 1 Representation of the projections of the periodic Delaunay triangulation into the 1-sheeted (left) and 9-sheeted cover (middle, right) of the 2-torus.

  dom (0, Λ) (red) is the intersection of the halfopen strips (blue) of V + = {v1, v2, v3}. (b) dom (0, Λ) of a 3-dimensional lattice (blue). Lattice basis B in black, V + ⊃ B in orange.

Figure 2

 2 Figure 2 Canonical domains and Voronoi relevant vectors for lattices in 2D and 3D.

  Finding the periodic neighbor of σ across the edge τ , which is σ n , the canonical representative of σ n .

  and similarly c 02 ≤ 0. If a basis is not Lagrange-reduced with 2|c 12 | > c 11 , we exchange b 2 for a shorter vector that forms a basis with b 1 : Let

Lemma 4 .

 4 There is a set O ≤2 of 51 offsets and a set O 3 of 24 offsets such that for any 3-dimensional lattice with basis B, there is a subset S of O 3 of size 6 such that if o ∈ O ≤2 ∪ S, then dom (λ(o), Λ) ∩ dom (0, 3Λ) = ∅. Proof. We use the observation that | v, w | ≤ v, v for any two Voronoi relevant vectors w, v ∈ V [13, Section 3.5]. Then the following claim holds:

  Results for fixed skew s = 0 and varying elongation . Results for = 4 and varying s. The inverse of the square length of the shortest non-zero lattice vector in black, with its y-axis on the right.

Figure 5

 5 Figure 5 For different lattices, the number of points inserted into a periodic triangulation until (blue) Del(ΓX)/Γ is simplicial for the first time, and (orange) Del(ΓX)/Γ fulfills Lemma 5. Each data point is the mean of 200 trials, and the bars represent the standard deviation.

Figure 6

 6 Figure 6Phase switching in 3D: Each point in the plot represents a lattice, with vol(dom(0, Λ))/sv(Λ) 3 on the x-axis and the number of points inserted until the switch to phase 2 occurred on the y-axis. Each y-value is obtained as the average over 200 point sets, each distributed uniformly at random.

2 0 2 0 37:14 Generalizing CGAL Periodic Delaunay Triangulations

  Adrian Bowyer. Computing Dirichlet tessellations. The computer journal, 24(2):162-166, 1981.

	References 1 Erik Agrell, Thomas Eriksson, Alexander Vardy, and Kenneth Zeger. Closest point search in lattices. IEEE Trans. Inform. Theory, 48(8):2201-2214, 2002. doi:10.1109/TIT.2002.800499. 2 Miklós Ajtai. The shortest vector problem in L2 is NP-hard for randomized reductions. In Proceedings of the Thirtieth Annual ACM Symposium on Theory of Computing, STOC '98, page 10-19, New York, NY, USA, 1998. Association for Computing Machinery. doi: 10.1145/276698.276705. Cubic [5] Lattice Lattice -Cubic Cubic FCC Λ1 Λ2 vol(dom(0,Λ)) sv(Λ) 3 -1.00 1.00 0.71 12.50 346.41 n switch -145 [4] 141 94 2519 89950 10 0 0.0000 0.0001 0.0002 0.0001 0.0003 0.0004 10 1 0.0000 0.0160 0.0033 0.0026 0.0044 0.0035 10 2 0.0004 0.1848 0.0461 0.0287 0.0460 0.0380 10 3 0.0049 0.5957 0.0858 0.0446 0.9812 0.6372 10 4 0.0487 0.9591 0.3832 0.1642 4.6602 16.5759 10 5 0.5679 4.8119 4.5153 2.7868 10.8139 362.7956 10 6 6.5974 93.9327 95.1447 51.5945 58.1568 517.9715 10 7 59.5152 2314.3618 2317.7867 1215.2648 1799.4515 2983.0943 3 Algorithm Euclidean [14] Table

E S A

https://members.loria.fr/Monique.Teillaud/CGAL_periodicDT_ESA20/ E S A

0 2 0

Manuel Caroli. Triangulating Point Sets in Orbit Spaces. PhD thesis, Université Nice Sophia Antipolis, 2010. URL: https://tel.archives-ouvertes.fr/tel-00552215.

Manuel Caroli, Aymeric Pellé, Mael Rouxel-Labbé, and Monique Teillaud. 3D periodic triangulations. In CGAL User and Reference Manual. CGAL Editorial Board, 5.0.2 edition, 2020. URL: http://doc.cgal.org/latest/Manual/packages.html#PkgPeriodic3Triangulation3.

6 Manuel Caroli and Monique Teillaud. 3D periodic triangulations. In CGAL User and Reference Manual. CGAL Editorial Board, 3.5 edition,

2009. 7 Manuel Caroli and Monique Teillaud. Delaunay triangulations of closed Euclidean d-orbifolds.Discrete & Computational Geometry, 55(4):827-853, 2016. URL: https://hal.inria.fr/

Georg Osang: This author is partially supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, grant no. 788183. Monique Teillaud: This author is partially supported by the grants ANR-17-CE40-0033 of the French National Research Agency ANR (project SoS) and INTER/ANR/16/11554412/SoS of the Luxembourg National Research fund FNR (https://members.loria.fr/Monique.Teillaud/collab/ SoS/).

is not in the interior of str(4v) and thus the interior of the 4-scaled domain. Therefore dom (ν, Λ) ∩ dom (0, 3Λ) = ∅.

, every lattice point λ can be written as a non-negative integer combination of three of these extended basis vectors, i.e.

This representation is unique up to permutation of basis vectors with the same coefficient.

Enumerating all dom (λ, Λ) at graph distance 3 from dom (0, Λ), we get lattice points λ with non-negative integer combinations of the following types:

All of these, except for type 3a + 2b + c, can be written as 3v or 3v + w for some Voronoi relevant vectors v, w ∈ V. This means that they are, or are adjacent to, a domain centered at 3v, and from Claim 2 it follows that they cannot intersect dom (0, 3Λ). A similar argument shows that none of the domains at graph distance 4 from dom (0, Λ) can intersect dom (0, 3Λ), and therefore also none at higher graph distance. Now the 51 offsets from O ≤2 are those that the argument from Lemma 3 yields, while O 3 contains the offsets corresponding to type 3a+2b+c. As by Lemma 2 there are 57 intersecting domains, only 6 of the offsets from O 3 can correspond to intersecting domains.


Xi Tian
email: tianxi@u.nus.edu

Qihang Zeng

Denys Nikolayev

John S Ho
email: johnho@nus.edu.sg

Conformal Propagation and Near-Omnidirectional Radiation with Surface Plasmonic Clothing

Keywords: Wearable antennas, spoof surface plasmons, smart textiles

Wireless technologies are essential components of wearable devices for applications ranging from connected healthcare to human-machine interfaces. Their performance, however, is hindered by the human body, which obstructs the propagation of wireless signals over a wide range of directions. Here, we demonstrate conformal propagation and near-omnidirectional radiation of wireless signals near the body using clothing integrated with spoof surface plasmonic structures. These structures, fabricated entirely from conductive textiles, induce wireless signals emitted by nearby antennas to propagate around the body as surface waves and radiate as propagating waves in directions otherwise obstructed by the body. We describe the procedure for designing textile-based spoof surface plasmon waveguides, radiating elements, and impedance matching sections for operation in the 2.4-2.45 GHz ISM band. Using a tissue phantom model of the human torso, we experimentally demonstrate 2.62 greater angular coverage by a dipole placed near the body compared to without the clothing.

I. INTRODUCTION

W IRELESSLY-connected wearable devices have many applications in emerging healthcare, communication, and human-computer interface technologies [START_REF] Han | Battery-free, wireless sensors for full-body pressure and temperature mapping[END_REF]. To create more functional devices, techniques to increase the performance of the radio-frequency wireless system (for example, Wi-Fi, Bluetooth, or 4G/5G) are essential because it represents a major constraint on the size, power consumption, and lifetime of the device. A major obstacle in the performance of wearable wireless technologies, however, is the presence of the human body. Because the human body obstructs the propagation of radio-frequency waves, the ability of radiofrequency signals to radiate into a wide range of directions around the body is greatly limited. This obstruction results in poor angular coverage that adversely impacts the reliability and energy efficiency of the wearable device.

Antenna design approaches have been extensively investigated to improve the radiation pattern emanating from devices worn on the body [START_REF] Lin | Characteristics of electric field and radiation pattern on different locations of the human body for in-body wireless communication[END_REF]- [START_REF] Lan | A free-access segmented coplanar waveguide for on-body communication[END_REF]. Classes of antennas developed for use in wearable devices include cavity slot antennas [START_REF] Haga | Characteristics of cavity slot antenna for body-area networks[END_REF], [START_REF] Moro | Wearable textile antenna in substrate integrated waveguide technology[END_REF], planar inverted-F antennas [START_REF] Gao | A wearable PIFA with an all-textile metasurface for 5 GHz WBAN applications[END_REF], [START_REF] Lin | A compact planar inverted-F antenna for 2.45 GHz on-body communications[END_REF], crosspatch antennas [START_REF] Bhattacharjee | A compact dual-band dual-polarized omnidirectional antenna for on-body applications[END_REF], [START_REF] Tak | A tapered double cross patch antenna for dual-band monopole-like radiation[END_REF], higher-mode antennas [START_REF] Paraskevopoulos | Higher-mode textile patch antenna with embroidered vias for on-body communication[END_REF], metallic button antennas [START_REF] Xiaomu | Wearable button antenna for dual-band WLAN applications with combined on and offbody radiation patterns[END_REF], [START_REF] Sanz-Izquierdo | Dualband wearable metallic button antennas and transmission in body area networks[END_REF], and steerable-beam antennas [START_REF] Islam | A 900 MHz beam steering parasitic antenna array for wearable wireless applications[END_REF]- [START_REF] Ha | Reconfigurable beam steering using a microstrip patch antenna with a U-slot for wearable fabric applications[END_REF]. Although these designs exhibit useful properties both on and off the body, they do not address the obstruction of wave propagation and remain limited in angular coverage when worn on the body. Another approach is to integrate the antennas into clothing, such as by fabricating them using conductive textile materials [START_REF] Da Conceicao Andrade | Reconfigurable textile-based ultra-wideband antenna for wearable applications[END_REF]- [START_REF] Hall | Antennas and propagation for body centric communications[END_REF]. By allowing multiple antennas to be placed on the body [START_REF] Lee | Omnidirectional vest-mounted body-worn antenna system for UHF operation[END_REF], this strategy can potentially circumvent obstruction and provide near-omnidirectional coverage around the body. However, current implementations rely on wired connections between the antennas and the signal source. These connectors present a major inconvenience to the user, and reduce the robustness of the clothing to washing and daily wear.

Recent advances in metamaterials have inspired the development of a class of thin, flexible waveguides that support spoof surface plasmons (SSP) at radio-frequencies [START_REF] Shen | Conformal surface plasmons propagating on ultrathin and flexible films[END_REF]. These structures have been used to create a wide range of microwave components and devices [START_REF] Chen | Wideband leaky-wave antennas loaded with gradient metasurface for fixed-beam radiations with customized tilting angles[END_REF]- [START_REF] Tian | Compact high-efficiency resonator antenna based on dispersion engineering of even-mode spoof surface plasmon polaritons[END_REF], including ultrawideband filters [START_REF] Wang | Supercompact and ultrawideband surface plasmonic bandpass filter[END_REF], [START_REF] Gao | An ultra-wideband surface plasmonic filter in microwave frequency[END_REF], ultrathin ring resonators [START_REF] Zhou | Spoof localized surface plasmons on ultrathin textured MIM ring resonator with enhanced resonances[END_REF], and high-isolation power dividers [START_REF] Zhou | Spoof surface plasmon polaritons power divider with large isolation[END_REF]. We have recently integrated such waveguides into clothing using conductive textile materials [START_REF] Tian | Wireless body sensor networks based on metamaterial textiles[END_REF]. Such clothing can significantly increase the efficiency of wireless transmission between wearable devices because of the excitation of surface waves by nearby antennas: wireless signals emitted by devices near the clothing can be induced to propagate conformally around the body rather than to radiate into the surrounding space. Importantly, these designs exploit the evanescent field associated with SSP modes to enable excitation by proximity to existing wireless devices without requiring connectors to the clothing. To interface with external devices, however, the integration of antennas and impedance matching sections is required to launch surface waves into radiating waves in the desired directions. Although their potential for integration with antennas was previously shown, the design of textile structures capable of redirecting radiation from a nearby devices to overcome obstruction by the body has not been previously explored.

In this paper, we demonstrate a surface plasmonic clothing that enables conformal propagation and near-omnidirectional radiation of wireless signals around the body. To overcome the obstruction of radiation by the body [Fig. 1(a)], our design incorporates several advances over prior work: (i) an optimized design procedure for textile-based SSP waveguides and patch antennas for launching surface waves, (ii) the integration of impedance matching sections, and (iii) detailed measurement and analysis of the radiation pattern generated around the body. We analyze the excitation of the surface modes and describe the procedure for designing such structures, including the SSP waveguide, radiating elements, and impedance matching section, for operation in the 2.4-2.5 GHz industrial, scientific, and medical (ISM) band. We fabricate the clothing using thin conductive textiles attached to a fabric substrate and experimentally demonstrate near-omnidirectional radiation around a tissue phantom with angular coverage 2.62 times greater than without the clothing.

II. CONCEPT AND DESIGN

A. Surface Plasmonic Clothing Concept Fig. 2 illustrates the design of the textile structure comprising of a SSP waveguide, matching sections, and radiating elements. When a radiating antenna is placed in close proximity to the SSP waveguide, the interaction between the source current distribution and the evanescent field of the surface mode induces a portion of the outgoing energy to propagate as surface waves. These surface waves propagate conformally along the waveguide, following the curvature of the body without reflection. Impedance matching sections then enable the surface waves to excite radiating elements placed around the body. These radiating elements launch the surface waves as radiation into the surrounding space, thereby enabling the radiation pattern emanating from the body to cover directions otherwise obstructed by the body. Crucially, the surface plasmonic structure affords control over the propagation and radiation of wireless signals through passive proximity to wirelessly-enabled wearable devices. The clothing enables the field radiated by a wearable device to be redirected around the body without physical connectors or modification of its built-in antenna.

Several design considerations are essential to implement the surface plasmonic clothing. First, the SSP waveguide must support surface plasmon-like modes over the operating frequency range. Second, the modes should be tailored to enable contactless excitation using antennas placed in close proximity to the body. Third, the transmission line needs to be interfaced to radiating elements that efficiently convert the surface wave into a propagating wave. Finally, the SSP structure must be fabricated using conductive textiles and integrated into clothing.

B. Surface Plasmonic Waveguides

We describe the design of the textile-based SSP waveguide in order to support surface modes with the desired properties over the operating frequency band. The design process uses a simplified model of the body comprising of a half-space in which the upper region (z > 0) is filled with freespace and a lower region (z < 0) with biological tissue of dielectric permittivity ε body . The textile structure is placed on the interface z = 0 and supports surface-plasmon-like modes propagating in the ±x directions.

We first introduce the design parameters by providing an approximate description of the surface modes. While the modes exhibit complex field distributions that depend on the detailed textile structure, the electric field can in general be approximated by the form

E ± (r) = p(y)   iα n /β 0 1   e ±iβx-αn|z| , (1)
where p(y) describes the structure-dependent field profile in the y direction, α n are the decay constant in the upper (n = 1) and lower (n = 2) regions, β is the propagation constant, and a time-dependency of e -iωt has been assumed [START_REF] Maier | Plasmonics: Fundamentals and Applications[END_REF]. In addition to the field distribution described in Eq. ( 1), the structure should also support surface-plasmon-like dispersion characteristics. In the simplest approximation, the propagation constant β follows the surface plasmon dispersion relation

β = k 0 ε 1 ε 2 /(ε 1 + ε 2 ),
where k 0 = ω/c is the wavenumber and ε 2 = 1-ω 2 p /ω 2 , with both ε 1 and ω p as fitting parameters. This relation results in a dispersion curve that lies to the right of the light line β = k 0 and approaches a horizontal asymptote at the frequency

ω sp = ω p / √ 1 + ε 1 .
The field distribution and dispersion relation together imply a bound surface mode that is analogous to an optical surface plasmon polariton at the metal-dielectric interface, although the material parameters ε 1 and ε 2 should be strictly regarded as fitting parameters without their usual physical interpretation.

The structure of the SSP textile consists of a planer patterned top layer, an intermediate fabric layer, and an unpatterned bottom layer as shown in Fig. 3(a). In the upper region z > 0, the mode has a decay constant α 1 given by α 1 = β 2 -k 2 0 . Since the dispersion relation requires β > k 0 , α 1 is purely real and the mode does not radiate into the upper region z > 0. The addition of the unpatterned bottom layer prevents the evanescent field from extending into the lower region α 2 → ∞, which introduces additional radiative losses into the body region due to coupling to modes with propagation constant β < ω √ ε body /c. While analytical expressions for dispersion relation supported by the SSP waveguides in unbounded free space are available [START_REF] Wan | Guiding spoof surface plasmon polaritons by infinitely thin grooved metal strip[END_REF], the increased complexity of our structure requires numerical simulation to obtain the surface mode characteristics.

To design the SSP structure, we use an eigenmode solver of CST Microwave Studio (Dassault Systems). Unit-cell analysis of the structure yielded the dispersion curves on Fig. 3(a) corresponding to a set of geometrical parameters. These dispersion curves exhibit the surface-plasmon-like characteristics with properties that can be tuned by varying the geometry. Using this approach, the geometrical parameters can be optimized to achieve a desired β and mode field distribution at a particular operating frequency. Fig. 3(b)-(c) show the electric field of the SSP mode for different geometrical parameters. The field is strongly confined within the width of the structure, but extends evanescently above the surface. The decay length α -1 1 is primarily controlled by the comb height h [START_REF] Shen | Conformal surface plasmons propagating on ultrathin and flexible films[END_REF], and is largely insensitive to the other comb parameters over a wide range of values [Fig. 3(d)]. The region between the top comb structure and bottom unpatterned layer is characterized by the thickness t and permittivity ε fabric of the underlying fabric, which are generally selected based on wearability considerations. The SSP mode exhibits a weak dependence on these parameters, which can be compensated by tuning h to obtain the required β. Using cotton fabric (ε fabric = 1.5) with thickness t = 2 mm, we obtained the structure parameters a = 4 mm, b = 6 mm, h = 21 mm, and p = 8 mm to achieve β = 0.46π cm -1 at 2.4 GHz. The dimensions of the structure represent feature sizes that are well within the tolerances of textile cutting processes.

C. Surface Mode Excitation

Owing to the evanescent field extending from the surface of the SSP textile, the surface modes can be wirelessly excited by nearby antennas without physical connectors. In contrast, most conventional microwave waveguides and transmission lines strongly confine the modal field within its structure and require connectors to efficiently interface with external devices. The surface excitation effect can be described by separating the electric field of the surface mode into forward (+) and backward (-) propagating components

E(r) = a + E + (r) + a -E -(r), (2)
where the a ± is the mode amplitude and E ± are normalized such that

1 2 |p(y)| 2 dy = α 1 α 2 /(α 1 + α 2 )
. When a source with current density J(r) is placed above the SSP clothing, the surface modes are excited with amplitude [34, Ch. 8.12.C]

a ± = - β 2ωε 0 J(r) • E ∓ (r) dr. (3)
Equation [START_REF] Mendes | On-body transmission performance of a novel dual-mode wearable microstrip antenna[END_REF] shows that excitation depends only on the overlap between the current density and the evanescent field above the surface. The current distribution can be provided by a radiating antenna placed in close proximity to the structure, in which case the excitation of surface modes is accompanied by suppression of radiation into the surrounding space. As the excitation occurs in the near-field of the antenna, far-field properties such as gain and directivity do not play a direct role in determining the amplitude of the excited mode. We characterize the efficiency of surface mode excitation at 2.4 GHz using a dipole antenna of resonant length. The excitation efficiency is defined as η = P ssp /P in , where P in is the input power into the dipole antenna and P ssp is the power leaving ports in the ±x directions of the SSP waveguide. Since the surface plasmon modes are transverse-magnetic, the transverse y components of the electric field is negligible [Fig. 3(b)] and the surface modes are predominantly excited when the dipole is oriented in the x or z directions. Fig. 4(a) shows the theoretically calculated η as a function of the distance d above the surface for a x-oriented dipole, using (3) and analytical expressions for the electric field and dipole current distribution. The excitation efficiency is high when d is comparable to the decay length α in the full-wave simulations in Fig. 4(b). This model yields an excitation efficiency in close agreement with analytical calculations, and shows that the power input into the antenna is transferred from radiation to surface modes [Fig. 4(c)].

To generate near-omnidirectional radiation, the excitation distance d should be fixed to provide the surface excitation necessary to radiate power equally in each direction. For the case of two secondary radiating elements-each excited by the surface wave propagating in the ±x direction-the excitation efficiency required is η = 2/3 such that the remaining 1/3 energy is directly radiated by the dipole. The designed SSP structure with h = 21 mm, as shown in Fig. 3(a), achieves this target efficiency when d is approximately 5 mm. In practice, d is largely determined by the design of the wearable device, such as thickness of the packaging around the antenna, and the manner by which it is worn on the body. The optimal d can be tuned over a wide range by varying the parameter h to control the decay constant α using the previously described procedure.

D. Radiating Elements and Matching Section

Conversion of surface waves into radiation is required to provide coverage in desired directions around the body. Previous studies of surface-to-propagation wave conversion for the SSP structure provide the basis for the design of such elements [START_REF] Zhang | High gain and high-efficiency millimeter-wave antenna based on spoof surface plasmon polaritons[END_REF], [START_REF] Du | Wideband fish-bone antenna utilizing oddmode spoof surface plasmon polaritons for endfire radiation[END_REF], although the structure must satisfy additional requirements in order to be integrated into clothing and achieve near-omnidirectional radiation. First, radiation from each element should be single-sided to prevent transmission into the human body, which results in losses and tissue heating. Second, the structure should be simple and planar to be compatible with textile fabrication processes and allow easy integration into clothing. Third, the directivity should be sufficient to provide the required angular coverage in each direction without interfering with adjacent elements.

To meet these criteria, we use microstrip patch antennas as the radiating elements due to their ease of fabrication, single-sided radiation characteristics, and low directivity. The interface between the SSP waveguide and the patch antenna uses an impedance matching section, while the ground plane is simply extended over the length of the patch antenna [Fig. 5(a)]. We design the impedance matching section using a standard gradient corrugated strip with a linear taper by varying the length of the matching section and identifying the minimum reflection to a port 50 Ω [START_REF] Wang | Supercompact and ultrawideband surface plasmonic bandpass filter[END_REF], [START_REF] Kianinejad | Low-loss spoof surface plasmon slow-wave transmission lines with compact transition and high isolation[END_REF]- [START_REF] Ma | Broadband and high-efficiency conversion from guided waves to spoof surface plasmon polaritons[END_REF]. The two conductive fabric layers are modeled using with a thickness of 0.1 mm and separated by a 2 mm dielectric substrate with ε fabric = 1.5 and loss tangent tan δ = 0.04 [START_REF] Salvado | Textile materials for the design of wearable antennas: A survey[END_REF] to model the textile substrate.

E. Integration into Clothing

To demonstrate control of radiation around the body, we design the textile structure shown in Fig. 6(a) for integration into clothing. The structure consists of a SSP waveguide with 35 unit cells that terminates at both ends with a matching section and a secondary radiating element. The entire structure wraps partially around the waist with a separation angle of θ, defined from the center of each radiating element, as illustrated in Fig. 6(b)-(c). When a radiating antenna is placed near the front of the body, energy is directed to propagate as surface waves around the body and radiate into backwards direction. As a result, the resulting radiation pattern can be nearly omnidirectional, even when the antenna is obstructed by the body. Although we use two secondary radiating elements to demonstrate this effect, the structure can be extended to incorporate additional radiating elements, each with higher directivity, in order to achieve greater control over the radiation pattern.

We characterize the radiation pattern generated when the structure wrapped around a cylindrical model of the human torso (40-cm diameter, ε body = 40) and a dipole antenna placed at distance d above surface. The radiation patterns in Fig. 6(e) show the angular coverage for d = 5-7 mm. The angular coverage is maximized when d = 6 mm, which is in close agreement with the distance corresponding to the required excitation efficiency in Fig. 4(c). The separation angle θ between the radiating elements affects the angular coverage as well. Compared with the -3 dB threshold widely used to measure the directivity of antennas in free-space, here we use a lower reference threshold of -10 dB in order to quantify the angular coverage around the body relative to the radiation pattern of a dipole placed in front of the body, which has a front-back ratio that exceeds 20 dB. As shown in Fig. 6(g)-(j), the radiation pattern is nearly omnidirectional for θ = 60 • -150 • . The angular coverage is maximized at θ = 120 • when the radiating elements and dipole antenna are equally distributed around the body. Fig. 7 compares the absolute gain of each antenna for θ = 120 • . These results demonstrate two key reasons why SSP clothing could improve the angular coverage. First, the angular redistribution of the radiated energy mitigates the body obstruction. Second, similar to the dipole with a ground plane in Fig. 7(b), the minimization of power absorption in tissues increases the total energy radiated into the surrounding space. For this reason, the SSP clothing enhances the angular coverage with minimal penalty on the broadside gain (< 5 dB, compared with dipole on the body without a ground plane).

We perform a full-wave simulation of the structure using a dipole antenna is placed on the center of the abdomen 6 mm above the transmission line. The stimulations use a voxel model of the human body (Donna, CST Microwave Studio), which captures the tissue heterogeneity of the body. Fig. 8(a)-(b) show the simulated electric field profile for the dipole antenna without and with the SSP device. A significant portion of the energy leaving the antenna can be seen to propagate as surface waves around the body, which radiate in the backward directions. As demonstrated in previous studies [START_REF] Shen | Conformal surface plasmons propagating on ultrathin and flexible films[END_REF], propagation of the SSP waveguide is highly robust to bending and deformation. Along the curvature of the body, simulations show that reflection losses are about 3.9% compared to the flat structure and propagation losses about 1.5%. clothing is 114 • at 2.4 GHz, compared to 278 • in free-space. When placed near the SSP clothing, however, the angular coverage is increased about 2.7 times to 314 • , providing nearly omnidirectional radiation. The simulated angular coverage results are summarized in Table I.

III. EXPERIMENTAL DEMONSTRATION

A. Experimental Setup

We experimentally demonstrate the operation of the clothing by measuring the radiation pattern emanating from a 20 cmdiameter cylindrical tissue phantom to mimic the human body. The surface plasmonic clothing is fabricated by using a laser cutter to pattern the structure from thin conductive textiles Cu/Ni nylon ripstop fabric sheets (Conductive Fabric 4713, Holland Shielding Systems). The textiles are then attached to a cotton fabric sheet using a fabric adhesive. The resulting device is free of connectors and is entirely made of textile, which enables it to conform to the body when worn and be robust to washing. To test the clothing, the device is wrapped conformally around the cylinder and excited using a wideband dipole antenna resonant at 2.4 GHz placed at a distance d above the surface. The entire setup is placed on a rotating stage in an anechoic chamber as shown in Fig. 9(a).

The tissue phantom comprises of a cylindrical polypropylene container (20-cm diameter) filled with an agarose gel to mimic the shape and electromagnetic properties human torso. The phantom solution uses agarose, water and sugar following the recipe described in Ref. [START_REF] Dabbagh | Tissue-mimicking gel phantoms for thermal therapy studies[END_REF]. We verified the phantom solution by measuring the dielectric permittivity spectrum with an open-ended dielectric probe. As shown in Fig. 9(b), the dielectric permittivity is approximately ε r = 54 within the 2.4-2.5 GHz ISM band, which is close to permittivity of muscle tissue.

B. S-parameter Measurements

We measured the scattering parameter S 11 to characterize the effect of proximity to the SSP structure on the power leaving the dipole.

C. Radiation Pattern Measurements

We measured the radiation patterns at 2.4 GHz and 2.45 GHz with and without the SSP clothing on the tissue phantom. Fig. 10(a)-(b) shows the gain normalized to the maximum in the xy-plane. In absence of the SSP clothing, the angular coverage at 2.4 GHz is 123 • and follows the typical broadside pattern associated with obstruction of propagation by the body. When the SSP clothing is placed on the phantom, however, energy is radiated in directions obstructed by the tissue phantom, resulting in a nearly omnidirectional angular coverage of up to 322 • . The measured radiation patterns exhibit a slight asymmetry due to imperfect alignment between the dipole antenna and body surface, as well as imprecision in the placement of the flexible fabric material against the tissue phantom. However, the measured angular coverage (Table II) is in close agreement with simulations based on the computational human body model (Table I). Notably, the near-omnidirectional pattern is achieved solely by proximity to the SSP clothing, without altering the transmit power or the antenna.

IV. ROBUSTNESS ANALYSIS

To study the effect of the surrounding environment on the textile structure, we evaluated SSP waveguide transmission when a human arm, modeled by a cylinder with a diameter of 5 cm with a uniform dielectric permittivity of ε body = 40, is placed over the structure. The arm represents a "worst case" obstruction since most other common environmental objects have dielectric permittivities that are less than or comparable to fabric. Fig. 11 These results suggest that careful placement of the matching sections at positions where the human arm usually locate, as is done in our design, can mitigate performance variations due to posture while allowing interactions with nearby antennas over a wide range of positions.

Fig. 11(c)-(d) show another effective approach to reduce the sensitivity of the SSP waveguide by adding an additional fabric layer over sections of the SSP waveguide where contact is expected. The additional fabric layer significantly reduces the effect of the arm, and can enable effective transmission ∆|S 21 | < 4 dB, which maintains significant enhancement in radiation in the backward direction, even at contact d 2 = 0 mm for a superstrate thickness of t = 4 mm.

In addition to its use for connector-free interfacing with nearby antennas, the SSP waveguide can also provide distinct advantages in robustness to mechanical damage. Fig. 11(e) illustrates a SSP waveguide and a microstrip line cut with a 0.5-mm gap to mimic mechanical tear during wearing. Simulations of the time-averaged electric field show that the surface wave propagating along the SSP waveguide can significantly traverse the gap, whereas the microstrip line incurs large reflection [Fig. 11(f)]. As shown in the Fig. 11(g), the transmission provided by the SSP waveguide over the gap is 9 dB greater than the microstrip line. Use of the SSP waveguide in regions vulnerable to tear may improve the robustness of the design in wearable applications.

V. CONCLUSION

We have demonstrated conformal propagation and omnidirection radiation of wireless signals emitted from the human body using clothing patterned with surface plasmonic structures. These structures, fabricated using conductive textiles, induce radio-frequency waves emitted by nearby antennas to propagate around the body as surface waves, where they are launched as radiation into the surrounding space. In contrast with prior work on textile-integrated antenna, the surface plasmonic design provides control over propagation and radiation through proximity to the antenna -connectors are not required. Analysis of surface wave excitation using standard waveguide excitation theory allows the amplitude of the excited surface mode to be predicted from the current distribution of the source. Using full-wave simulations and measurements, we verified that such textiles can overcome the obstruction of radiation by the human body and provide nearly omnidirectional angular coverage about 2.62 times larger than in absence of the structures. These results highlight the potential of clothing to engineer the electromagnetic environment around the body in order to enhance the integration of technology into everyday life.

Fig. 1 .

 1 Fig. 1. Function of the surface plasmonic clothing. (a) Illustration of spoof surface plasmonic (SSP) structure on clothing. (b) Without the structures, the body obstructs the propagation of radio-frequency waves. (c) The structure induces the emitted wireless signals to propagate as surface waves around the body (red arrow) and then radiate around the body (blue arrows).

Fig. 3 .

 3 Fig. 3. Textile SSP waveguide design. (a) Geometrical parameters of the structure and dispersion curves for varying h. The remaining parameters are a = 4 mm, b = 6 mm, p = 8 mm, t = 2 mm, and ε fabric = 1.5). (b) Electric field of the mode in the xy plane at 2.45 GHz for h = 21 mm. (c) Electric field of the mode in the yz plane for varying h and b. White dashed lines show the decay length 3α -1 . (d) Propagation constant β as a function of h, p, a, and t.

-1 1 ≈ 7 Fig. 4 .

 174 Fig. 4. Excitation of surface modes. (a) Theoretically calculated excitation efficiency of the surface modes for h = 19 mm. (b) Electric field profile when an x-oriented electric dipole (2.4 GHz resonant length) is placed at distance d above the surface. (c) Simulated excitation efficiency of surface modes as a function of d for varying waveguide width h.

Fig. 5 (

 5 b) shows the simulated electric field in the xzplane cutting through the center of the patch antenna when the structure is placed on a planar interface between air and tissue (ε r = 50). The device efficiently converts surface waves excited by energy input into from a waveguide port into radiation into the upper half-space with minimal reflection |S 11 | < 10 dB over the 2.4-2.45 GHz band [Fig. 5(c)]. Due to the ground plane, the presence of the body does not significantly affect the performance of the structure.

Fig. 5 .

 5 Fig. 5. Radiative elements and matching section. (a) Schematic of structure including the SSP waveguide, matching section, and antenna. (b) Electric field magnitude in the xz longitudinal plane for excitation of the SSP waveguide. (c) Simulated |S 11 | of the radiating element.

6 .

 6 Surface plasmonic clothing design. (a) Illustration of clothing integrated with the SSP waveguide and two secondary radiating elements for redirecting radiation from a device worn on the front of the body. (b) Top view of the flattened structure. (c) Side view of the curved structure. (d) Cross-section view of the structure on a simplified cylindrical model of the body. The angular separation between the radiating elements is θ and d is the distance of the dipole antenna (representing the wearable device) d from the surface. (e) Simulated radiation patterns for d = 5-7 mm. (f)-(j) Simulated radiation patterns for θ = 60-180 • . Shaded area indicates the angular coverage where the gain is higher than -10 dB.

Fig. 7 .

 7 Fig. 7. Antenna radiation and gain on a cylindrical model of the human torso. (a)-(b) Dipole without ground, (c)-(d) dipole with ground, (e)-(f) textilebased patch antenna, and (g)-(h) dipole with the SSP textile structure. The gain is normalized to the overall highest gain, which is 4.19 dB for the dipole with ground. The cylinder has a diameter of 40 cm and a uniform dielectric permittivity of ε body = 40.

Fig. 8 .

 8 Fig. 8. Antenna radiation in a computational human body model. (a)-(b) Simulated electric field in the transverse plane of the body (a) without and (b) with the SSP clothing when a dipole antenna (solid white arrow) is placed near the front of the body. Dashed white lines indicate the direction of wave propagation. (c)-(d) Radiation patterns with and without the surface plasmonic clothing at (c) 2.4 GHz and (d) 2.45 GHz. The shaded areas indicate the angular coverage where the gain is greater than -10 dB.

Fig. 9 .

 9 Fig. 9. Experimental setup. (a) Measurement setup using a 20-cm diameter cylindrical container filled with tissue-mimicking hydrogel. (b) Measured dielectric permittivity of hydrogel as a function of frequency. (c) Measured |S 11 | at different excitation distance d.

  Fig. 9(c) shows the |S 11 | spectrum for d varying from 3 mm to 20 mm. At close proximity (up to 6 mm), the |S 11 | curve exhibits two dips located at about 2.36 GHz and 2.45 GHz, and is lower than -10 dB across the ISM band. The dips correspond to the resonances of the patch antenna as well as the modes associated with the finite length of the SSP waveguide. At about d = 10 mm, the spectral dips vanish, indicating that the SSP mode is weakly excited. Beyond d = 15 mm, the S 11 spectrum is insensitive to d, which is consistent with the dipole directly radiating into the surrounding space without exciting surface modes.

  (a)-(b) illustrate the placement of the arm in close proximity to different parts of the SSP clothing. The change in transmission ∆|S 21 | (relative to |S 21 | in absence of the arm) becomes significant when the distance of the arm over the SSP waveguide section of the structure d 1 is less than 1 cm, and exceeds -3 dB when d 1 < 4 mm. However, the

Fig. 10 .

 10 Fig. 10. Measured radiation patterns with and without the SSP clothing. Radiation pattern measured at (a) 2.4 GHz and (b) 2.45 GHz. Shaded areas indicate the angular coverage where the gain is higher than -10 dB. Black circles indicate angular position of the secondary radiating elements.

Fig. 11 .

 11 Fig. 11. Robustness analysis. (a) Illustrations of an environmental object (arm) placed at a distance d 1 above different sections of the surface plasmonic clothing. (b) Change in transmission ∆|S 21 | (relative to no object) as a function of d 1 for different positions of the arm. (c) Illustration of a fabric superstrate with thickness of tc added above the SSP waveguide when the arm is placed above the surface of superstrate at a distance d 2 . (d) ∆|S 21 | as a function of d 2 for different superstrate thickness. (e) Schematics of microstrip line and SSP waveguide cut with a 0.5 mm gap. (f)-(g) Time-averaged electric field amplitude and ∆|S 21 | after cutting.

This work was supported by the National Research Foundation Singapore (NRFF2017-07) and Ministry of Education Singapore (MOE2016-T2-2-016 and MOE2016-T3-1-004). (Xi Tian and Qihang Zeng contributed equally to this work.Body Antenna SSP structure