Jean-Daniel Boissonnat

Mathijs Wintraecken

Triangulating stratified manifolds I: a reach comparison theorem

Keywords: 1998 ACM Subject Classification Computational Geometry and Object Modeling Keywords and phrases Reach, Manifolds, surfaces, comparison theory Digital Object Identifier 10.4230/LIPIcs.xxx.yyy.p

In this paper, we define the reach for submanifolds of Riemannian manifolds, in a way that is similar to the Euclidean case. Given a d-dimensional submanifold S of a smooth Riemannian manifold M and a point p ∈ M that is not too far from S we want to give bounds on local feature size of exp -1 p (S). Here exp -1 p is the inverse exponential map, a canonical map from the manifold to the tangent space. Bounds on the local feature size of exp -1 p (S) can be reduced to giving bounds on the reach of exp -1 p (B), where B is a geodesic ball, centred at c with radius equal to the reach of S. Equivalently we can give bounds on the reach of exp -1 p • exp c (B c ), where now B c is a ball in the tangent space T c M, with the same radius. To establish bounds on the reach of exp -1 p • exp c (B c ) we use bounds on the metric and on its derivative in Riemann normal coordinates.

This result is a first step towards answering the important question of how to triangulate stratified manifolds.

Introduction and motivation

Motivation: Triangulating stratified manifolds

Triangulating of manifolds with boundary, submanifolds of Riemannian manifolds and stratified manifolds, is an important problem because of applications in dynamical systems, physics, and chemistry. Examples of stratified manifolds in applications include conformation spaces of molecules, such as discovered for cyclo-octane [START_REF] Martin | Topology of cyclo-octane energy landscape[END_REF].

There already exists a vast literature on triangulating surfaces in R 3 , see for example [START_REF] Amenta | A simple algorithm for homeomorphic surface reconstruction[END_REF][START_REF] Mourrain | Meshing of surfaces[END_REF][START_REF] Cazals | Delaunay triangulation based surface reconstruction[END_REF][START_REF] Cheng | Delaunay Mesh Generation[END_REF][START_REF] Dey | Curve and Surface Reconstruction: Algorithms with Mathematical Analysis[END_REF], as well as results on the triangulation of general manifolds in arbitrary dimensional Euclidean space [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF][START_REF] Cheng | Manifold reconstruction from point samples[END_REF]. The triangulation of solid objects with boundary and stratified solid objects in R 3 has been the object of study in [START_REF] Tamal | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF][START_REF] Tamal | Localized Delaunay refinement for volumes[END_REF][START_REF] Dey | Localized Delaunay refinement for piecewise-smooth complexes[END_REF][START_REF] Oudot | Meshing Volumes Bounded by Smooth Surfaces[END_REF][START_REF] Rineau | Meshing Volumes bounded by Piecewise Smooth Surfaces[END_REF], and Chapter 15 of [START_REF] Cheng | Delaunay Mesh Generation[END_REF].

Algorithms have been proposed for separating the strata of stratified manifolds [START_REF] Bendich | Inferring local homology from sampled stratified spaces[END_REF]; the resulting strata being manifolds with boundaries, see also [START_REF] Bendich | Stratification learning through homology inference[END_REF]. The triangulation of Riemannian manifolds also has received significant interest in the last few years, see [START_REF] Boissonnat | Delaunay triangulation of manifolds[END_REF][START_REF] Dyer | Riemannian simplices and triangulations[END_REF].

The reach

The reach is a key concept used in triangulations and related problems, see for example the overviews [START_REF] Mourrain | Meshing of surfaces[END_REF][START_REF] Cheng | Delaunay Mesh Generation[END_REF][START_REF] Dey | Curve and Surface Reconstruction: Algorithms with Mathematical Analysis[END_REF], because it not only gives bounds on the local curvature, but gives more global information, such as how close different parts of the set S lie to one another.

In Euclidean space, the reach of a subset S is defined as the minimal distance of a point on S to the medial axis (i.e. the set of points that have more than one closest point on S).

The reach of subsets of Euclidean space was introduced by Federer [START_REF] Federer | Curvature measures[END_REF]. The reach is denoted by rch. The same concept was called the condition number by Niyogi, Smale and Weinberger [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF]. The local feature size, a local version of the reach and called the reach at a specific point by Federer, was introduced into the Computational geometry community by Amenta and Bern [START_REF] Amenta | Surface reconstruction by Voronoi filtering[END_REF]. Remarkably, the authors of both [START_REF] Amenta | Surface reconstruction by Voronoi filtering[END_REF] and [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF] seem to have been unaware of the results of Federer. In this paper, we'll generalize the reach to submanifolds of Euclidean space to submanifolds of arbitrary Riemannian manifolds. Defining and understanding such a generalization of the reach to this setting is essential to be able to leverage the results on triangulations in Euclidean space we mentioned above. We'll explain this in more detail now.

Let M be a Riemannian manifold. We denote by T p M the tangent space at p and remind ourselves that the exponential map exp p maps vectors emanating from p in T p M to geodesics on M emanating from p with equal length, while also preserving the angles between these.

The reach of a submanifold of a Riemannian manifold can be defined in much the same way as in Euclidean space, roughly as the distance to the medial axis (i.e. the set of points that have more than one closest point on S). In the manifold setting some technical extra condition on the reach, involving the geometry of M is also necessary, see Section 2.

Suppose that we have bounds on the reach of a submanifold1 S of a Riemannian manifold M. We now want to consider a point p ∈ M that is not too far from S. In fact, for the sake of the argument, we'll assume that S lies well inside a convex geodesic ball centred at p. This assumption guarantees that exp -1 p (S) is well defined and we can talk about the (global) reach of exp -1 p (S). Working with more local assumptions (such as the local feature size) should be possible, but makes the argument even more technical. Now we would like to give bounds on the reach of exp -1 p (S) in T p M. With bounds on the reach, we can use results on triangulations in Euclidean space to triangulate in T p M and previous results on distortion Riemannian simplices under the exponential map to map the triangulation back to the manifold.

In preparation for the larger problem of triangulating Riemannian manifolds with boundary, submanifolds of a Riemannian manifold, and stratified Riemannian manifolds, this paper focuses on the essential prerequisite of understanding the reach or local feature size of exp -1 p (S) in terms of the reach of the submanifold S and the geometry of the manifold M in which S is embedded. Along the way, we will give an introduction to some fundamental results in Riemannian geometry, such as (a simplified version of) the Toponogov comparison theorem.

Future work

As we mentioned above, this work is part of a larger project on the triangulation of several kinds of Riemannian manifolds with substructure, for which a significant amount of work needs to be done. This is also reflected in the title of this paper.

To be able to go from this result to the triangulation of stratified manifolds we also need to extend and/or generalize results on The triangulation of piecewise smooth manifolds, and stratified manifolds in three dimensional Euclidean space. The results on the triangulation of Riemannian manifolds to accommodate constraints. This will be reported on in forthcoming papers.

2

The reach Similar to the definition for manifolds embedded in Euclidean space we define the reach for C ∞ compact Riemannian (sub)manifolds S ⊂ M. We shall denote the normal space of S at a point x by N x S, and the bundle by N M.

Definition 1 (The reach). We let the medial axis ax M (S) be the set of points in M that do not have a unique closest point on S, with respect to the Riemannian metric. We denote the projection of a point x in M on the closest point on S by π S (x). The pre-reach prch M (S) is then the shortest distance between ax M (S) and S. We now define the reach rch M (S) to be

min{prch M (S), ι M }, (1)
where ι M is the injectivity radius of M.

The injectivity radius is the largest radius r such that for any x, exp x restricted to the open ball centred at x with radius r in T x M is a diffeomorphism onto its image. Adding this bound to our definition of the reach is essential, because we would like to have that a tangent ball to S is indeed a topological ball.

Remark. In the mathematics literature, the closure of the medial axis is often called the cut locus in this context, see [START_REF] Chazal | Stability and finiteness properties of medial axis and skeleton[END_REF]Section 2] or [START_REF] Lieutier | Any open bounded subset of R n has the same homotopy type as it medial axis[END_REF]. We have chosen to use the terminology and concept most common in the computational geometry community in view of future applications.

We now give the result, which is an extension of Theorem 4.8.8 of Federer [START_REF] Federer | Curvature measures[END_REF] to the setting of Riemannian manifolds: Theorem 2 (Tubular neighbourhood). Let B NpS (r), be the ball of radius r centred at p in the normal space N p S ⊂ T p M of a C ∞ manifold with reach rch M (S), where r < rch M (S). For every point x ∈ exp p (B NpS (r)), we have π S (x) = p.

The proof of Theorem 2 is not so difficult and mainly follows Hirsch [START_REF] Hirsch | Differential Topology[END_REF]Section 4.5] and Spivak [30, Appendix I of Chapter 9] with some variations.

We start with Lemma 19 of Spivak [START_REF] Spivak | A comprehensive introduction to differential geometry: Volume I[END_REF]Chapter 9]:

Lemma 3. Let X be a compact metric space and X 0 ⊂ X a closed subset. Let f : X → Y be a local homeomorphism such that f | X0 is injective. Then there exists a neighbourhood U of X 0 such that f | U is injective.
With this lemma we can prove an embedding result, for which we have to make the following definition: Definition 4. Let N B denote the -neighbourhood in N S of S, that is the neighbourhood of all points closer than to S, where we identify S ⊂ N S via the zero section. Moreover write f : N B → M, for the map defined by sending S ⊂ N S to S and fiberwise sending N p S to exp p (N p S). Theorem 5. There exists an > 0 such that f : N B → M is a (global) homeomorphism onto its image, that is an embedding.

The proof in the appendix combines arguments from Section 4.5 of [START_REF] Hirsch | Differential Topology[END_REF] with Appendix I of Chapter 9 of [START_REF] Spivak | A comprehensive introduction to differential geometry: Volume I[END_REF], where small variations of this statement can be found.

We now define N S(r) fiberwise as those points in N p S that are at a distance r from p. We refer to r as the radius of the tubular neighbourhood. For any smooth manifold S embedded in M and x ∈ M, we know that the geodesic from x to π S (x) is normal to S at π S (x), as a direct consequence of the Gauss lemma (the Gauss lemma that refers to the exponential map). It follows that, for all such that f :

N B → M is a homeomorphism, each point in f (N S(r)) is a distance r from M, for all 0 < r < .
Because S is compact, so is N S(r) as is the closed -neighbourhood N B . Clearly f is continuous. In general, a continuous bijection from a compact to a Hausdorff space is a homeomorphism. This means that the map f from both N S(r) and the closed rneighbourhood N B r to their images are homeomorphisms, that is embeddings, if and only if f from both N S(r) and N B r are injective.

Because N B r is closed and N S(r) is the boundary of N B r , we have that if f from N S(r) and N B r to their respective images are a homeomorphisms, that is embeddings, there is an η > 0 such that the same holds for N S(r + η) and N B r+η . This means that we can make the neighbourhood larger, until a critical radius r where N S(r ) no longer gets embedded by f . Here we assume that r < ι M . Moreover, as we have just seen this is equivalent to the map not being injective, so at least two points x and y are mapped to the same image. Because we assume that r < ι M , we cannot have that x, y ∈ N q S for some q ∈ Su, if it were, then ι M = r . This means that the image of x, which is also the image of y, under f does not have a unique closest point on the manifold.

We now have the following: Lemma 6. The reach r is the smallest radius, with r ≤ ι M , such that f restricted to N S(r) and N B r are no longer homeomorphisms to their images. For all r ≤ rch M (S) and all points in f (N S(r)) are a distance r from S.

In particular we have proven Theorem 2.

Note that Theorem 2 immediately yields, Corollary 7. Let S be a submanifold M and p ∈ S. Any open ball B that is tangent to S at p and whose radius r satisfies r ≤ rch M (S) does not intersect M.

Proof. Let r < rch M (S). Suppose that the intersection of M and the open ball is not empty, then the π S (c) = p contradicting Theorem 2. The result for r = rch M (S) now follows by taking the limit.

Bounds on the reach

As mentioned, the main goal of this paper is to understand the reach or local feature size of exp -1 p (S) in terms of the reach of the submanifold S and the geometry of the manifold M in which S is embedded.

Approach

The lower bounds on the reach of exp -1 p (S), in terms of the reach of S and geometric properties of the manifold, follow by considering tangent balls to S that do not contain points of S, from various viewpoints, namely from the manifold and from various tangent spaces. The argument will be explained in detail using Figure 1.

M S c q T c M T p M exp -1 p exp -1 c exp c • exp -1 p c exp -1 c (q) exp -1 c (S) p exp -1 c (p) exp -1 p (S) exp -1 p • exp c (B c ) exp -1 p (q) B B c = exp -1 c (B)
Figure 1 An overview of the approach.

Let B be a geodesic ball in M centered at c, tangent to S at a point q and with radius the reach of S. As in Euclidean space and by Corollary 7, B has an empty intersection with S. B is indicated in red in the top of Figure 1.

The exponential map exp c gives coordinates on a neighbourhood of the manifold, as does exp p . We will use the inverse of the exponential maps exp c and exp p , and write B c = exp -1 c B and B p = exp -1 p B for the images (in T c M and in T p M respectively) of the geodesic ball B. Note that B c is an Euclidean ball which implies that its reach is equal to its radius, which is also the (geodesic) radius of B. See bottom left part of Figure 1. Note also that B c is tangent to S c = exp -1 c (S) and that B p is tangent to S p = exp -1 p (S). The composition exp -1 p • exp c gives a transformation between the two coordinate neighbourhoods, as indicated by the arrow from the bottom left to right in Figure 1.

Thanks to the Toponogov comparison theorem and a higher order variant, we have bounds on the metric as well as on the derivatives of the metric, both expressed in the coordinates induced by the exponential maps exp p and exp c . The bounds on the metric and its derivatives can then be used to give bounds on the first and second order derivatives of the transformation exp -1 p • exp c , the bottom arrow in Figure 1. Thanks to a result by Federer [START_REF] Federer | Curvature measures[END_REF] one can find a bound on the reach after the transformation, based on the reach of the original. This gives bounds on the reach (in T p M) of B p , that is on the reach of exp where the minimum over B is a minimum over geodesic balls in M of radius rch M (S) that are tangent to S at q. For the local feature size one takes q fixed instead of minimizing over S.

Overview

The outline of this section is as follows:

In section 3.2 we focus on Riemann normal coordinate systems. Thanks to standard comparison theorems such as the Toponogov comparison theorem [START_REF] Heintze | A general comparison theorem with applications to volume estimates for submanifolds[END_REF][START_REF] Karcher | Riemannian comparison constructions[END_REF], we are able to give bounds on the metric in Riemann normal coordinates. The Riemann normal coordinates of a neighbourhood of p are those coordinates that are found by lifting the metric to the tangent space at p via the exponential map. The work by Kaul [START_REF] Kaul | Schranken für die Christoffelsymbole[END_REF] provides us bounds on the Christoffel symbols in the Riemann normal coordinate neighbourhood, and thus indirectly bounds on the derivative of the metric. In Section 3.3 we'll study the coordinates transformation exp -1 c • exp -1 p , the bottom arrow in Figure 1. In Section 3.3.1, we first see how we can go from bounds on the metric in the Riemann normal coordinates to bounds on the coordinate transformation exp -1 c • exp -1 p . In Section 3.3.2, we'll be applying a result by Federer [START_REF] Federer | Curvature measures[END_REF], that will yields the reach.

Bounds on the metric and Kaul's bound on the Christoffel symbols

In this section we review the bounds on the metric and Kaul's bounds on the Christoffel symbols, see [START_REF] Kaul | Schranken für die Christoffelsymbole[END_REF]. The expressions for these bounds have been simplified by Von Deylen [31, Section 6], at the cost of weakening the bounds. We shall make use of his simplification.

Here and throughout d denotes the dimension of M and we adopt the Einstein summation convention, that is, if an index occurs twice (once as an upper and once as a lower index) we sum over this index without writing a summation sign. We denote by g ij (x) the metric in this coordinate system at x and δ ij is the standard Euclidean metric. The inverses of the metrics are denotes by g ij (x) and δ ij , respectively. The norm with respect to the Riemannian metric is denoted by | • | g while the norm with respect to the Euclidean metric is denoted by | • | E . Distances on M will be denoted by d M . As before, ι M is the injectivity radius.

The Christoffel symbols are

Γ κ µν = 1 2 g κλ (∂ µ g λν + ∂ ν g λµ -∂ λ g µν ),
where ∂ µ denotes the partial derivative with respect to the coordinate x µ The Christoffel symbols are used to express the covariant derivative of a vector field v κ in local coordinates,

∇ ν v κ = ∂ ν v κ + Γ κ µν v µ .
R σ µνλ will denote the Riemann curvature tensor in local coordinates. We refer to Do Carmo [START_REF] Carmo | Riemannian Geometry[END_REF] and Spivak [START_REF] Spivak | A comprehensive introduction to differential geometry: Volume I[END_REF], as some of the standard texts introducing these concepts.

We shall now assume that the curvature and its derivative are bounded, that is in any orthonormal coordinate system,

|R σ µνλ | ≤ R max (2)
|∇ κ R σ µνλ | ≤ R ∇ max . ( 3
)
We now consider the Riemann normal coordinates around p:

x : (x 1 , . . . , x d ) → exp p (x i E i )
for some orthonormal basis E i of T M . We now have the following simplification of the Toponogov comparison theorem:

Lemma 8 (Lemma 6.8 of [START_REF] Von Deylen | Numerische Approximation in Riemannschen Mannigfaltigkeiten mithilfe des Karcher'schen Schwerpunktes[END_REF]).

If d M (p, x) ≤ r, with R max r 2 ≤ π 2 4 and r ≤ ι M 2 , then |g ij (x) -δ ij | ≤ R max r 2 .
Moreover the result of Kaul [START_REF] Kaul | Schranken für die Christoffelsymbole[END_REF] simplifies to: Lemma 9 (Lemma 6.9 of [START_REF] Von Deylen | Numerische Approximation in Riemannschen Mannigfaltigkeiten mithilfe des Karcher'schen Schwerpunktes[END_REF]).

If d M (p, x) ≤ r, with R max r 2 ≤ π 2 4 and r ≤ ι M 2 , then |g κλ (x)Γ κ µν (x)v µ w ν u λ | ≤ 10R max r + 5R ∇ max r 2 , for all v, u, w ∈ T x M such that |u| g = |v| g = |w| g = 1.
We now note that thanks to Lemma 8, we have that

||v| 2 g -|v| 2 E | = g ij v i v j -δ ij v i v j ≤ ij |g ij -δ ij ||v i v j | ≤ ij R max r 2 |v i ||v j | = R max r 2 i 1|v i | 2 ≤ R max r 2 i 1 i |v i | 2 , = dR max r 2 |v| 2 E ,
where we used the Cauchy-Schwartz inequality, and made the summation explicit at several points for reasons of clarity. Combining this with Lemma 9, we see that

|g κλ (x)Γ κ µν (x)| ≤ 10R max r + 5R ∇ max r 2 (1 -R max d r 2 ) 3/2 .

Using that

∂ ν g κµ (x) = g κλ (x)Γ λ µν + g µλ (x)Γ λ κν
and taking absolute values, we find the following corollary

Corollary 10. If d M (p, x) ≤ r, with R max r 2 ≤ π 2 4 and r ≤ ι M 2 , then |∂ ν g κµ (x)| ≤ 20R max r + 10R ∇ max r 2 (1 -R max d r 2 ) 3/2 .
We now recall two results from linear algebra:

Let E be a d × d-matrix, then 1 √ d E ∞ ≤ E 2 ≤ √ d E ∞ ,
see for example (2.3.11) of [START_REF] Golub | Matrix Computations[END_REF].

If G = I + E, where G and E are d × d-matrices, I denotes the identity matrix and E 2 ≤ 1, then

I -G -1 2 ≤ E 2 1 -E 2 ,
see for example [START_REF] Horn | Matrix Analysis[END_REF]Section 5.8].

If now, E ∞ ≤ c, with √ dc < 1, then I -G -1 ∞ ≤ dc 1 - √ dc .
With this result and Lemma 8, we immediately have the following

Corollary 11. If d M (p, x) ≤ r, with √ dR max r 2 ≤ 1 and r ≤ ι M 2 , then |g ij -δ ij | ≤ dR max r 2 1 - √ dR max r 2 .
We'll also make use of the following result [START_REF] Horn | Matrix Analysis[END_REF]Corollary 6.3.4], which we'll formulate as a lemma, Lemma 12. Let E be an d × d-matrix, and G = I + E, with I the identity matrix. If λ is an eigenvalue of G, then |λ -1| ≤ E 2 .

We now immediately have

Corollary 13. If d M (p, x) ≤ r, with R max r 2 ≤ π 2 4 and r ≤ ι M 2 , any eigenvalue λ of g ij (x) now satisfies |λ -1| ≤ √ dR max r 2 .

From bounds on the metrics to bounds on the coordinate transformations

The starting point of this section is the following: We are given a metric in Riemann normal coordinates at two different points. We want to study the coordinate transformation between these coordinates systems, based on our knowledge of the metric in these two coordinate systems.

In fact, we assume we have bounds on the first and second order derivatives of the metric in both coordinate systems. These bounds yield bounds on the first and second derivatives of the coordinates transformation. It is easy to see why this is so, by considering the limit case: Suppose both metrics are the Euclidean metric, then the transformation from one coordinate system to the other is a combination of a rotation and translation.

From bounds on the coordinate transformation, a result of Federer [START_REF] Federer | Curvature measures[END_REF] gives bounds on the reach of exp -1 p (B) where B is a geodesic ball centred at c with radius r, assuming that p and c are not too far from each other. Here we emphasize that the radius of the ball B equals the radius of exp -1 c (B). We will consider a coordinate transformation from a coordinate system x to a coordinate system y. Because the emphasis is on coordinate transformations we'll follow a different convention in this section, and only this section, and use Latin indices. We'll use the indices a, b, c, e, f for y-coordinate system and the indices i, j, k, l, m for x and write

y a = T a i x i + Q a ij x i x j + O(x 3
). Here we assumed that the coordinate systems are chosen such that the origins are mapped to one another, which can be done without loss of generality. Q a ij is symmetric in i and j. We have

g ij (x) = g ab (y(x)) ∂y a ∂x i ∂y b ∂x j g ij (x) = g ij (0) + (∂ k g ij )(0)x k + O(x 2 ) g ab (y) = g ab (0) + (∂ c g ab )(0)y c + O(y) 2 y c = T c m x m + O(x 2 ) ∂y a ∂x i = T a i + Q a ik x k + O(x 2 )
Combining these gives

g ij (0) = g ab (0)T a i T b j , (4)
(∂ k g ij )(0) = ∂ c g ab T c k T a i T b j + g ab (Q a ik T b j + T a i Q b jk )

Bounds on the transformations

With the concepts and notations developed in the previous section, we can give bounds on the transformation in terms of bounds on the metric and its derivatives. We find the following.

Lemma 14. Let x and y be two coordinates systems for the same point on the manifold and assume that the metric and its derivatives in these coordinates systems are as follows

g ij (x) = g ij (0) + (∂ k g ij )(0)x k + O(x 2 ) g ab (y) = g ab (0) + (∂ c g ab )(0)y c + O(y) 2 .
We assume further that any eigenvalue λ of g ij is bounded by |λ -1| ≤ A, the any eigenvalue λ of g ab is bounded by | λ -1| ≤ B, and the entries of g ef are bounded from above by 1 + C.

Moreover we assume that for all i, j, k we also have that |∂ k g ij | ≤ ∂g max,x , and for all a, b, c, that |∂ a g bc | ≤ ∂g max,y . Now we have that the coordinate transformation between x and y,

y a = T a i x i + Q a ij x i x j + O(x 3
), satisfies the following constraints: The Lipschitz constants, or the metric distortion of the linear approximation T , are bounded by

√ 1 + A/ √ 1 -B and its inverse by √ 1 + B/ √ 1 -A, and |Q a ij | ≤ 3d 2 ∂g max,x (1 + C) √ 1 + B √ 1 -A + d 3 ∂g max,y (1 + A)(1 + C) 1 -B .
The proof of this statement can be found in the appendix.

Using Federer's estimates: from bounds on the coordinate transformation to bounds on the reach

In this section, we are finally able to give the bounds on the reach by applying Theorem 4.19 of Federer [START_REF] Federer | Curvature measures[END_REF]:

Theorem 15 (Federer). Let S be a subset of R d with rch(S) > t > 0, and s > 0. If

f : {x | d(x, S) < s} → R d is a C 2 diffeomorphism such that |D f | ≤ M |D f -1 | ≤ N |D 2 f | ≤ P,
where D denotes the derivative, then

rch( f (S)) ≥ min{sN -1 , (M t -1 + P ) -1 N -2 }.
We can now combine this result with the estimates of the previous sections. We want to investigate how an empty tangent ball to S transforms under the exponential map. Because a geodesic ball is also an Euclidean ball in the tangent space of its centre (lifted via the exponential map), this is equivalent to giving bounds on the reach of this ball under the map exp -1 p • exp c . The bounds on the reach under the map exp -1 p • exp c use almost all previous results in this paper: In particular the bounds on the metric and its derivatives are given in Lemma 8, and Corollaries 10, 11, and 13, while Lemma 14 tells us how to go from bounds on the metric to bounds on the coordinate transformation. Federer's result now gives us the reach after the transformation.

Our main result now reads:

Theorem 16. Let M be a smooth d-dimensional Riemannian manifold whose curvatures are bounded as follows:

|R σ µνλ | ≤ R max ( 2
)
|∇ κ R σ µνλ | ≤ R ∇ max . ( 3
)
Suppose that r p and r c are the radii of geodesic balls centred at p and c respectively such that B(c, 2r c ) ⊂ B(p, r p ), where we have made the centres and radii explicit, √ dR max r 2 p ≤ 1 and r p ≤ ι M 2 . Then the reach of exp -1 p (B(c, r c )) ⊂ T p M is bounded. Specifically,

rch(exp -1 p (B(c, r c ))) ≥ min{r c N -1 , (M r -1 c + P ) -1 N -2 },
where

M = √ 1 + A √ 1 -B N = √ 1 + B √ 1 -A P = 12∂g max,x (1 + C) √ 1 + B √ 1 -A + 8∂g max,y (1 + A)(1 + C) 1 -B , with A = √ dR max r 2 c B = √ dR max r 2 p . C = dR max r 2 p 1 - √ dR max r 2 p ,
and ∂g max,x = ∂g max (r c ), ∂g max,y = ∂g max (r p ) with

∂g max (r) = 20R max r + 10R ∇ max r 2 (1 -R max d r 2 ) 3/2
Remark. Our theorem also gives a bound on the reach (or local feature size, if we want to concentrate on the local case) of a submanifold under exp -1 p by simply applying the theorem to every tangent ball (or every tangent ball tangent to a given point, for the local feature size). Here we naturally assume that all the tangent balls satisfy the conditions of the theorem.

A Appendix: Proofs

Proof of Theorem 5. We first note that the fiberwise restriction of f to B NpS (r) is a diffeomorphism for each r < ι M . Because M is C ∞ , f is smooth and we can consider the derivative T (p,0) f at a point (p, 0) ∈ N S. The tangent space splits as follows T (p,0) N S = T p S ⊕N p S. T (p,0) f is the identity if restricted to the tangent space of M as well as to a fiber. This gives that f is an immersion and in particular a local homeomorphism on its image because the codimension is zero. Due to Lemma 3, f : N B → M is injective for some sufficiently small > 0.

Proof of Lemma 14. We write G for the matrix g ij (0), and we assume that the eigenvalues λ i , are not far from 1, that is |λ i -1| ≤ A for all i and some A ≥ 0. Similarly, we write G for g ab (0), and assume that its eigenvalues λi are bounded by | λi -1| ≤ B for some B ≥ 0. We'll also write G = o t A D A o A with o A the orthogonal matrix that diagonalizes G and D A the diagonal matrix with the eigenvalues of G on the diagonal, that is diag(λ i ). We let S A denote the matrix with square roots √ λ i on the diagonal, that is We shall now consider the quadratic term. We start with

S A = diag( √ λ i ). And similarly, G = o t B D B o B and S B = diag( λi ). Now (4) gives G = T t GT o t A S A IS A o A = T t o t B S B IS B o B T o t A S t A IS A o A = T t o t B S t B IS B o B T I = S -t A o A T t o t B S t B IS B o B T o t A S -1 A , which means that S B o B T o t A S -1 A = o,
(∂ k g ij )(0) = ∂ c g ab T c k T a i T b j + g ab (Q a ik T b j + T a i Q b jk ).
Reshuffling and permuting the indices and changing the sign for the last equation gives

(∂ k g ij )(0) -∂ c g ab T c k T a i T b j = g ab Q a ik T b j + g ab T a i Q b jk (∂ i g jk )(0) -∂ c g ab T c i T a j T b k = g ab Q a ji T b k + g ab T a j Q b ki -(∂ j g ki )(0) + ∂ c g ab T c j T a k T b i = -g ab Q a kj T b i -g ab T a k Q b ij .
Adding the terms yields

(∂ k g ij )(0) + (∂ i g jk )(0) -(∂ j g ki )(0) -∂ c g ab T c k T a i T b j -∂ c g ab T c i T a j T b k + ∂ c g ab T c j T a k T b i = g ab Q a ik T b j + g ab T a i Q b jk + g ab Q a ji T b k + g ab T a j Q b ki -g ab Q a kj T b i -g ab T a k Q b ij = 2g ab Q a ik T b j ,
and thus

(∂ k g ij )(0) + (∂ i g jk )(0) -(∂ j g ki )(0) (T -1 ) j e g ef -∂ c g ae T c k T a i g ef -∂ c g eb T c i T b k g ef + ∂ e g ab a k T b i g ef = (∂ k g ij )(0) + (∂ i g jk )(0) -(∂ j g ki )(0) -∂ c g ab T c k T a i T b j -∂ c g ab T c i T a j T b k + ∂ c g ab T c j T a k T b i (T -1 ) j e g ef =2g ab Q a ik T b j (T -1 ) j e g ef =Q f ik ,
The idea now is the following: If we assume that the left hand side of the previous equation is close to zero (this is in line with Corollary 10 because we assume that the derivatives of the metric are not too large if the neighbourhood is not too small), g ab is close to δ ab , and T b j is close to a rotation, all entries of Q a ij have to be close to zero too. Let us now assume that for all k, i, j we have that

|∂ k g ij | ≤ ∂g max,x ,
and for all a, b, c, that |∂ a g bc | ≤ ∂g max,y .

We'll also assume that entries of g ef are bounded in absolute value by 1 + C. We will use that for a tensor U µν , with |U µν | ≤ U max , and the coordinates of vectors v µ , w µ are bounded by |v µ | ≤ v max and |w µ | ≤ w max , we have that

|T µν v µ w ν | ≤ µ,ν T max v max v max = d 2 T max v max v max (5)
where we made the summation explicit, as well as its obvious generalization. Thanks to the triangle inequality we now have

|Q f ik | = ((∂ k g ij )(0) + (∂ i g jk )(0) -(∂ j g ki )(0)) (T -1 ) j e g ef -∂ c g ae T c k T a i g ef -∂ c g eb T c i T b k g ef + ∂ e g ab T a k T b i g ef ≤| ((∂ k g ij )(0) + (∂ i g jk )(0) -(∂ j g ki )(0)) (T -1 ) j e g ef | + |∂ c g ae T c k T a i g ef | + |∂ c g eb T c i T b k g ef | + |∂ e g ab T a k T b i g ef | ≤3d 2 ∂g max,x (1 + C) √ 1 + B √ 1 -A + d 3 ∂g max,y (1 + A)(1 + C) 1 -B

- 1 p

 1 • exp c (B c ). This is indicated in the bottom right of Figure 1, by the green Euclidean ball inside the red deformed ball. Consider the Euclidean ball b p ⊆ B p ⊂ T p M with radius rch(B p ) that is tangent to ∂B p at exp -1 p (q), where ∂ indicates the boundary. The Euclidean ball b p is also tangent to S p at exp -1 p (q) and, since B is empty of points of S, B p does not intersect S p nor does b p . The reach of S p is lower bounded by the minimum radius of any such Euclidean ball b p , that is rch(S p ) ≥ min q∈S min B rch(B p ) = min q∈S min B radius(b p ),

  with o an orthogonal transformation. This in turn implies that T = o t B S -1 B oS A o A is close to an orthogonal transformation, if A and B are close to zero. The Lipschitz constant of a composition of function is the product of the Lipschitz constants and thus the Lipschitz constant of T is bounded by (1-A)(1-B) and (1 + A)(1 + B) respectively. Because we have that for any vector |v µ | ≤ |v|, where the first | • | should be read as an absolute value and the second as the norm, and A ij = e t i Ae j , where the e i denote basis vectors, the entries of T are bounded by from above by √ 1 + A/ √ 1 -B. By the same argument we have that the entries of T -1 are bounded by √ 1 + B/ √ 1 -A.

The result for a stratum would be similar.

Acknowledgements

We thank Eddie Aamari, David Cohen-Steiner, Isa Costantini , Fred Chazal, Ramsay Dyer, and Alef Sterk for discussion and Pierre Pansu for encouragement. The research leading to these results has received funding from the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement No. 339025 GUDHI (Algorithmic Foundations of Geometry Understanding in Higher Dimensions).


Moussa Déthié Sarr

François Delobel
email: delobel@sancy.clermont-universite.fr

Michel Misson
email: misson@sancy.clermont-universite.fr

Ibrahima Niang
email: iniang@ucad.sn

Automatic and dynamic network establishment for linear WSNs

Keywords:

The deployment of Wireless Sensor Networks (WSNs) in the mostly linear large structures, such as rivers, pipelines, etc, suffers from being the worst case for classical addressing schemes, e.g. Distributed Address Assignment Mechanism (DAAM) or stochastic addressing for Zigbee. Using DAAM for physical topologies composed of long lines of nodes connected together wastes addresses and generates orphan nodes.

We show in this paper the inherent limitations of classical (DAAM, stochastic) and specialized (usually cluster-orientated) addressing schemes for Linear WSNs. DiscoProto, is an addressing and routing scheme which builds a logical network corresponding to a corresponding physical linear network without any knowledge of physical topology. In this paper, we show thanks to a realistic simulation using Castalia (Omnet based simulator) that DiscoProto avoids waste of addresses and allows a high association ratio. We also propose a dynamic version of our protocol called Dynamic DiscoProto in the second part of the paper. Dynamic DiscoProto allows to add new nodes or new branches in an existing linear network.

Introduction

Wireless sensor networks (WSNs) are networks formed by a collection of sensor nodes with limited resources (power, bandwidth, memory, CPU and communication range) [START_REF] Akyildiz | A survey on sensor networks[END_REF][START_REF] Yick | Wireless sensor network survey[END_REF]. Because sensor nodes used to perform these applications have limited resources, the main challenges in WSNs are scalability, self-organization and low energy consumption.

WSNs are used to monitor many applications fields (rivers, health monitoring [START_REF] Baker | Wireless sensor networks for home health care[END_REF], volcanoes monitoring [START_REF] Lorincz | Deploying a wireless sensor network on an active volcano[END_REF], railways [START_REF] Zimmerling | Localized power-aware routing in linear wireless sensor networks[END_REF], bridge structure [START_REF] Kim | Health monitoring of civil infrastructures using wireless sensor networks[END_REF]).

Most applications are based on the logical topologies defined in WSNs standards (star, mesh, cluster-tree) [START_REF]Computer Society. Part 15.4: Wireless Medium Access Control (MAC) and Physical Layer (PHY) Specifications for Low-Rate Wirelees Personal Area Networks (WPANs)[END_REF][START_REF]ZigBee Specification. ZigBee Standard Organization[END_REF].

However, existing network formation protocols and nodes addressing assignment schemes [START_REF] Ameer | A survey on clustering algorithms for wireless sensor networks[END_REF][START_REF] Yick | Wireless sensor network survey[END_REF][START_REF]ZigBee Specification. ZigBee Standard Organization[END_REF] proposed in WSNs1 are not suitable for applications monitoring pipelines, rivers, bridges, tunnels or railway structures because these protocols assume that the nodes can be easily organized in a mesh or cluster-tree logical topology. In fact, these application fields imply a mostly linear deployment [START_REF] Jawhar | Linear wireless sensor networks: Classification and applications[END_REF][START_REF] Jawhar | A framework for pipeline infrastructure monitoring using wireless sensor networks[END_REF][START_REF] Meng-Shiuan Pan | Address assignment and routing schemes for zigbee-based long-thin wireless sensor networks[END_REF] to monitor a mostly linear infrastructure.

Many classical WSNs addressing techniques are inefficient in this case because linear structures are often the worst case of their addressing techniques (by opposition to a well-balanced tree or a stochastic deployment). Moreover, in large cluster-trees or linear networks, it is difficult to provide an addressing mechanism covering all nodes because of the limitation due to the network parameters (Cm, Rm, Lm), the small (16-bits) size of logical addresses and the waste of address. Stochastic addressing has high probability to generate address conflicts on large networks and requires large routing tables.

To avoid those limitations, specialized protocols and mechanisms are proposed in [START_REF] Jawhar | A framework for pipeline infrastructure monitoring using wireless sensor networks[END_REF][START_REF] Zimmerling | Localized power-aware routing in linear wireless sensor networks[END_REF][START_REF] Meng-Shiuan Pan | Address assignment and routing schemes for zigbee-based long-thin wireless sensor networks[END_REF], we showed in [START_REF] Déthié Sarr | Automatic discovery of topologies and addressing for linear wireless sensors networks[END_REF] the drawbacks of these protocols and mechanisms for large linear WSNs.

The goal of our paper is to propose a robust network formation mechanism not only for linear but also for partially linear wireless network topologies 2 . This mechanism provides:

• Linear network formation functionality without any prior knowledge of physical topology, without limitations of the logical topologies and without manual configuration of any kind,

• Efficient hierarchical routing and addressing based on exact size of subtrees,

• Dynamic node management after network formation phase (adding or removing nodes, or whole branches), preserving most of the benefits of hierarchical routing,

• Address block reallocation to ensure the reliability of LWSNs in harsh environment to keep the mechanism robust.

The paper is organized in the following manner. In section II, we present the related works in WSNs. In section III we discuss of the state of the art of addressing mechanisms suitable for linear networks in ZigBee cluster-tree. We introduce our automatic discovery of topology and addressing scheme for linear WSNs in Section IV. In section V, we present our proposition, a Dynamic DiscoProto mechanism. In section VI, we describe our implementation in Castalia simulator. Section VII presents the results of our simulation and section VII conclude this paper with a summary of our contribution and perspectives of future works.

2 State of the art 2.1 Cluster-tree network formation in ZigBee IEEE 802.15.4 standard [START_REF]Computer Society. Part 15.4: Wireless Medium Access Control (MAC) and Physical Layer (PHY) Specifications for Low-Rate Wirelees Personal Area Networks (WPANs)[END_REF] introduces two different categories of nodes. A Full Function Device (FFD) which is a sensor node with router capability and a Reduce Function Device can only be used as end-device node ZigBee [START_REF]ZigBee Specification. ZigBee Standard Organization[END_REF] defines three roles of nodes: a ZigBee Coordinator, ZigBee Routers (ZR) that participates in routing messages for mesh and cluster-tree topologies and ZigBee End Devices (ZED) to provide sensing data capabilities. FFD devices can be used as ZC, ZR or ZED nodes. However, RFD devices can only be used as ZED nodes.

In ZigBee cluster-tree network formation [START_REF]ZigBee Specification. ZigBee Standard Organization[END_REF][START_REF] Cunha | On the use of ieee 802.15.4/zigbee as federating communication protocols for wireless sensor networks[END_REF], the parent node, for every new child to associate, chooses an address and assigns it to its new child. A PAN coordinator uses parameters such as the maximum number of children per node (Cm), the maximum number of children routers (Rm), and the maximum tree depth (Lm) to perform a Distributed Address Assignment Mechanism (DAAM) and to guarantee a unique and hierarchical addressing. For this, each parent node uses the Cskip function to calculate the size of address sub-block depending on its depth (d) in the network.

In this cluster-tree network, the network address block is calculated with the Equation ( 1)):

Amax = Cskip(0) × Rm + Cm -Rm (1)
Where:

Cskip(d) = 1 + Cm × (Lm -d -1), if Rm = 1 1+Cm-Rm-Cm×Rm Lm-d-1 1-Rm , otherwise
When a ZED node joins the network and becomes the n th child end device of the parent node at the depth (d), its address is calculated with equation [START_REF] Yick | Wireless sensor network survey[END_REF].

A child = A parent + Rm × Cskip(d) + n (2)
However, if it is a ZR node which joins the network and becomes the n th child router of the parent node at depth (d), then its address is calculated with equation 3:

A child = A parent + 1, if n = 1 (3)
A child = A parent + (n -1) × Cskip(d), if n > 1
Cm, Rm, Lm parameters have to be manually defined before the construction phase of the network according to its shape and strongly impact the network formation process. Those values can be optimally chosen only for perfectly balanced trees, with the same number of sons for each node. In fact, in real situation, it is difficult to predict good values for parameters Cm, or Lm. Poorly chosen values causes problems of non-associated nodes (orphan nodes) in large networks. In [START_REF] Park | Distributed borrowing addressing scheme for zigbee/IEEE 802.15.4 wireless sensor networks[END_REF][START_REF] Meng-Shiuan Pan | The orphan problem in zigbee wireless networks[END_REF], the authors propose a mechanism to reduce the number of non-associated nodes.

Figure 1: DAAM problem

On the other hand, with DAAM, it is difficult to expand network area particularly in Linear WSNs. Indeed, if we consider trees with high depth (close to linear), the addressing space of ZigBee (2 16 = 65535), limits the possible configurations (Cm, Rm, Lm) to (2, 2, 15), [START_REF] Baker | Wireless sensor networks for home health care[END_REF][START_REF] Baker | Wireless sensor networks for home health care[END_REF][START_REF] Ameer | A survey on clustering algorithms for wireless sensor networks[END_REF], [START_REF] Lorincz | Deploying a wireless sensor network on an active volcano[END_REF][START_REF] Lorincz | Deploying a wireless sensor network on an active volcano[END_REF][START_REF]Computer Society. Part 15.4: Wireless Medium Access Control (MAC) and Physical Layer (PHY) Specifications for Low-Rate Wirelees Personal Area Networks (WPANs)[END_REF]. In linear WSNs, we consider topologies whith depths that can exceed 15 hops.

DIstributed Borrowing Addressing Scheme for ZigBee/IEEE 802.15.4 Wireless Sensor Networks

The authors in [START_REF] Park | Distributed borrowing addressing scheme for zigbee/IEEE 802.15.4 wireless sensor networks[END_REF] propose the Distributed Borrowing Addressing (DIBA) mechanism to solve association failure in distributed address assignment mechanism used in ZigBee. All router nodes in this network broadcast beacon frames including the Available Address Count (AAC), which is the number of children that the routers can add.

A new node joining the network chooses the parent node which has the biggest AAC. However, if the parent node has no available address or if its tree depth is Lm, it uses the address borrowing mechanism to obtain an address for the new node : When the parent node P receives an Association Request from a new node and has reached the limit of Cm, it broadcasts an Address Borrowing Request to borrow addresses from its neighbor nodes. Some of them responds with an Address Borrowing Response including Available Address (AA) and AAC. 1: Routing tables of some nodes shown in Figure 2 In dense networks, Address Borrowing Response messages cause many collisions and can result in association failures. In addition, when all available addresses have been borrowed in 1-hop neighbors, the candidate node can no longer be associated because multi-hop borrowing is not allowed in DIBA. Moreover, DIBA increases the complexity of routing in the case of numerous address borrowing because purely hierarchical routing becomes insufficient (see Table 1).

Address Assignment and Routing Schemes for ZigBee-Based Long-Thin WSNs

In [START_REF] Meng-Shiuan Pan | Address assignment and routing schemes for zigbee-based long-thin wireless sensor networks[END_REF], the authors consider a special case of ZigBee cluster-tree where each cluster is composed of nodes deployed on line segment. In addition, each node in the Long-Thin WSNs (LTWSNs) has router capabilities and two special nodes need to be defined (see Figure 3):

• The cluster head node is the ZigBee router node connected to the cluster bridge node of the parent cluster.

• A cluster bridge node is the ZigBee router node connected to the cluster head node of the child cluster. In order to deploy this mechanism, the head and bridge of each cluster have to be manually defined. Then, a logical network is formed (see Figure 4) where each cluster in a physical network is represented by a node in a logical topology. For the node address assignment, the authors divide the 16-bit network address in two parts: The first part of m bits is used for cluster addressing and the second part of 16 -m bits for node addressing in a cluster.

Then, as in ZigBee Cluster-Tree Addressing, two parameters are fixed: CLm (Cluster maximum Length) and CCm (Cluster maximum Children) to compute the CCskip function defining the starting Cluster ID (CID) of each cluster in a logical network.

• if CCm > 1 CCskip(d) = 1 -CCm CLm-d 1 -CCm , (5)
Each cluster at depth d in the logical network assigns Cluster ID to its n th cluster child using Equation 6

CID Child = CID P arent + (n -1) × CCskip(d) + 1, (6)
if CCm >= 2

• If CCm = 1, the logical network is a segment line with no branch and cluster addresses are assigned sequentially.

After the CID assignment, each node in the cluster starts a network setup sending HELLO messages with 16-bit network address where the CID part is its cluster ID and the nodeID is 0. The network coordinator starts a nodeID address assignment broadcasting beacon with nodeID set to 0. If a node u in the neighborhood is not already associated, it has to respond with an Assoc Request to the beacon sender having the strongest signal strength. Each beacon sender v with an address (CID v , nodeID v ) sorts Assoc Request senders by signal strength in a List L. Then it sets nodeID = nodeID v + 1 and assigns sequentially nodeID u = nodeID u + 1 to Assoc Request senders having a same CID. After assigning nodeID for all received requests, the beacon sender either delegates this role of beacon sender to the bridge node or to the last node in L.

This approach suffers the same limitations as DAAM: The maximum number of clusters that can be addressed is 2 m and the maximum addressable node is 2 (16-m) . In addition, the network depth CLm is limited by the chosen parameter m (see the example in Equation ( 7)).

for CCm = 2, Amax ≤ 2 m ⇒ 2 × (2 CLm -1) ≤ 2 m (7) CLm ≤ ln (2 m-1 +1) ln 2
Finally, the head and bridge of each cluster have to be manually selected during the network deployment. Therefore, human intervention is also necessary to elect another node as cluster-head or bridge when one of this node fails. This manual selection method of the cluster-head or bridge can cause a network partitioning. In fact, all traffic of clusters is routed by a cluster-head to bridge node of its parent cluster.

Automatic Discovery of Topology and Addressing for Linear WSNs

In [START_REF] Déthié Sarr | Automatic discovery of topologies and addressing for linear wireless sensors networks[END_REF], the authors propose a discovery and a network formation mechanism adapted for linear network environment. We confirm this work with a realistic simulation and extend it to become dynamic.

Initial formation of network

[13] performs the network setup without knowledge of the topology parameters (maximum depth, maximum number of children by node, position of the node) and also provides a hierarchical addressing scheme to ensure efficient (and memory-saving) data routing.

In the discovery network phase, [START_REF] Déthié Sarr | Automatic discovery of topologies and addressing for linear wireless sensors networks[END_REF] propose a distributed algorithm called DiscoProto to generate a spanning tree minimizing several criteria. For this, let us consider in a graph a set of connected nodes CN and a set of unconnected nodes UN. Let c ∈ CN and u ∈ U N . DiscoProto incrementally constructs a connected network with the best link (c, u) having the local best cost. Note that the cost being a local optimum3 , several nodes can be added at the same time, providing they are far enough not to interact. So the formation phase of the network is actually done in parallel.

The choice of the best association father/son depends on three criteria (cost) sorted by importance:

• First, the number of common neighbors with c and u, common(c, u). This parameter represents the proximity of two nodes. Closest noded are associated first.

• Second, the number of sons that the connected node c has already associated, called sons(c). If two connected nodes c 1 and c 2 have the same number of common neighbors, i.e common(c 1 , u) = common(c 2 , u), u will be associated with the node with fewer sons between c1 and c2. This criterion ensures a good linearity of the network and keeps it as close as possible to the physical topology.

• Third, the sum of neighbors of c and u. If two nodes u 1 an u 2 have same two first criteria, the node with low connectivity will be associated first (fail first heuristic, as we will use a distributed greedy algorithm).

Those values are aggregated together as a single objective value by the following equation:

objective(c, u) =α × common(c, u) -sons(c) -β × (neighbors(c) + neighbors(u)), (8)
with α = 10 and β = 1 1000 .

To discover a linear network and to associate all nodes, each node goes thru five different states.

Neighborhood Discovery State

All awaken nodes broadcast Hello messages. Nodes that receive these messages update their neighbor tables. This step is important because computing the objective value for an association requires the knowledge of the neighborhood.

Wait for Father State

After neighborhood discovery, all nodes stay in this state until receiving a FatherOffer message, and then go into next step.

Son Collection State

In this state, connected nodes, i.e. nodes that are already associated, start to broadcast FatherOffer messages in the neighborhood. This message contains information needed to compute the objective value of the association (number of neighbors and number of sons). Each unassociated neighbor that received FatherOffer messages computes the association objective value with the given function (see Equation ( 8)). Then, it answers to the FatherOffer sender with unicast SonOffer message containing the objective value.

If no SonOffer messages are received from potential parents until the expiration of the timer for collect offers, connected nodes can retry sending FatherOffer message up to three times.

When the timer of a connected node sending FatherOffer expires, it selects the best SonOffer in its neighborhood and broadcasts its value using a ChallengeOffer message to the other connected nodes. This challenge is locally propagated (up to three hops) so that only the best father/son association in a neighborhood will be chosen. If a ChallengeOffer sender node receives an ChallengeOffer with the higher objective value from its neighborhood, then it cancels its potential father/son association and restarts a new Son Collection. At the expiration of ChallengerOffer timer, nodes that have not cancelled their potential father/son association start the association step.

After the challenge phase in Sons Collection state, the winner node sends an AssociationAccept message to the node selected as new son. This node answers with an AssociationAck and memorizes it as new father. When an AssociationAccept sender node does not receive an AssociationAck until expiration of the timer, it can retry up to three time before giving up the selected node and restarting a new Son Collection state. If a selected node receives an AssociationAccept when it is already connected, it sends an AssociationFailed to the sender node to cancel the useless association. Therefore, the latter gives up the selected node and restarts a new Son Collection state.

A new son node sets its potential father node as new father node and starts a new Son Collection state. When a FatherOffer sender did not receive SonOffer message after three repetitions, it stops collecting sons and goes to the Wait for SonSize state.

Wait for SonSize State

This state lasts until all the node's sons have provided the size of their subtree. Then the node computes the value size = n i=1 SonSize(i) + 1 and sends this value to its father node using PropaSons messages. When all nodes have left Wait for SonSize State, all the nodes in the network know the size of all their subtrees. A node leaving Wait for SonSize State enters the Wait for Address state.

Wait for Address State

After receiving all son sizes, the coordinator computes appropriate address block (AB) for the entire network using this following Equation 9.

AB = ( n i=1 SonSize(i) + 1) × F skip (9)
As the addresses can be smoothly distributed in the network, spare addresses can be uniformly distributed (F skip coefficient in Equation ( 9)).

Every node (including the sink) uses PropaAddr and PropaAddrAck messages to dispatch address blocks to its sons. Each node that receives the PropaAddr message containing the range [a, b], chooses a as its own address. Then, it keeps addresses from [a + 1, a + F skip] as available addresses for future association. F skip is the number of available addresses. The remaining interval [a + F skip + 1, b] is divided between its N sons proportionally to the size of their subtree. Figure 5 shows an example of address distribution after the network formation. Table 2: Son address table of node 6 in Figure 5 3

.2 Routing algorithm

A routing algorithm based on our hierarchical addressing ensures efficient data routing. Each node in the network maintains a son address table containing for each son: its logical and physical addresses and the first and last logical addresses allocated in its subtree (see

memory size = son(C) × (3 × network address size (10)
+mac address size)

The worst care complexity of routing is trying all the rules up to the fourth one, with a complexity of numbers of entry.

Suppose that the node n with an address addr receives a data packet, it performs the following procedures:

1. If node n is the destination, the packet is accepted and transmitted to the application layer.

If not, the node verifies if the destination address belongs to one of the entries of the son address

table. If it is found, then the packet is sent to the son node corresponding to the matching entry.

3. Otherwise, the packet is sent to the father node.

Dynamic DiscoProto

In this section, we introduce a new mechanism for DiscoProto to dynamically add sparse nodes or branches to an existing linear WSNs. We also present an address reallocation mechanism in case nodes have exhaust their spare addresses. We consider that the network is already formed and that all connected nodes have received their own address block. A new node that needs to join the network starts by broadcasting HELLO messages until expiration of hello timer cycle. When an associated node i.e., a connected node with an address, receives a HELLO message, it updates its neighbor table and responds by broadcasting HELLO messages. At expiration of the hello timer, associated nodes which received HELLO messages from new node enter Son Collection state (described in section 3.1.3). After associating the new node, the selected father node (with logical address myaddr ) waits until receiving the size N of the new son subtree node before entering Addresses Propagation state.

Network join procedure

When the selected father receives the size of all its subtrees, it performs the following tasks:

1. If the father node has enough spare addresses (i.e., N ≤ R -nbrN ewSon) then it attributes the first available address myaddr + nbrN ewSon to the new son node, where nbrN ewSon is the number of new sons connected after the network setup.

2. Otherwise, it sends an AddressRequest message to the sink including N + 1 (see Figure 6). The sink answers with an AddressResponse message including the new address block [a, b] where a = lastEndBlock + 1 and b = lastEndBlock + F skip × (N + 1) + 1.

3. When the father receives the new address block range from the sink, it keeps the address range [a + 1, a + F skip + 1] for future associations and the remaining address range is dispatched between its newly added sons proportionally to their subtree size using PropaAddr.

The join procedure for new branches operates in the same manner as for new spare nodes (see Figure 7). However, the addressing mechanism only starts after the collection son procedure for all nodes in the branch. The first connected node of the branch waits to collect all son size subtree before sending it to its father to prevent useless messages and address block fragmentation. The father node of the new node performs the above procedures to distribute the address block to the new branch.

Address block reallocation

In our Linear WSNs, connected nodes can exhaust their available addresses by different manner:

• Connected nodes have distributed all available addresses to new nodes or branches,

• Connected nodes have not enough available addresses for new branches that need to join the network,

• New associated nodes have new candidate nodes or branches that need to join the network For all these cases, connected nodes send an AddressRequest message to the sink and the sink answers with AddressResponse. In order to maintain addressing, nodes traversed by the AddressResponse add a new entry in their routing tables. This entry contains start Block, end Block and the destination node of the AddressResponse (see 7f Note that requesting large blocks instead of multiple single addresses avoids bloating routing tables.

Routing Algorithm

The routing mechanism based on our hierarchical addressing mechanism and the information in the routing table is the same as described in subsection 3.2. Son address table maintains 16 bits prefix, stat Block and end Block addresses and 64 bits MAC addresses. The total memory size of son addresses table can be calculated as in Equation [START_REF] Jawhar | A framework for pipeline infrastructure monitoring using wireless sensor networks[END_REF]. Let son(C) be the sons number of node C and newAddrBlock the number of times that new blocks are forwarded to the subtree of node C. memory size = (son(C) + newAddrBlock) ×(3 × network address size [START_REF] Jawhar | A framework for pipeline infrastructure monitoring using wireless sensor networks[END_REF] +mac address size)

The Table 5 shows an example of entry of node 3.

Implementation in Castalia Simulator

In Automatic Discovery of Topology and Addressing for Linear WSNs [START_REF] Déthié Sarr | Automatic discovery of topologies and addressing for linear wireless sensors networks[END_REF], we validated the topology formation mechanism on a home-made simulator developed in Ruby, using a green thread per node. This environment had no realistic low layer (MAC and physical layer) and perfect propagation: Indeed, every message was received by all nodes located in range.

To validate the protocol from [START_REF] Déthié Sarr | Automatic discovery of topologies and addressing for linear wireless sensors networks[END_REF], realistic propagation conditions and collision have to be taken into account. Therefore, we propose an implementation of Automatic Discovery of Topology in the WSNs Castalia environment, which has a proper physical and MAC layer (S-MAC protocol [START_REF] Ye | An energy-efficient mac protocol for wireless sensor networks[END_REF]).

Castalia environment

Castalia [START_REF] Boulis | Castalia : A simulator for wireless sensor networks[END_REF][START_REF] Tselishchev | Experiences and lessons from implementing a wireless sensor network mac protocol in the castalia simulator[END_REF] is a simulator specially designed for WSNs based on OMNet++ platform. It implements advanced channel model and radio model, realistic modeling of RSSI and carrier sensing and it is designed for adaptation and expansion. Castalia provides many path loss models, including the log-normal shadowing model which is well adapted to WSNs [START_REF] Zuniga | Analyzing the transitional region in low power wireless links[END_REF]. It also has a useful energy model. Physical node parameters are based on real specifications of CC2410 (see Tables 6 and7) motes.

Name normal

Link cost computation

Castalia computes RSSI and LQI values during transmission between nodes in a network. We use this advantage to slightly improve the objective function in Equation ( 12) to compute link cost in network graph by adding RSSI measurements. If two potentials parents of the same node have the same number of common neighbors, the same number of associated sons, and the same number of neighbors, we select the potential parent with the better RSSI value.

objective(c, u) =α × common(c, u) -sons(c) -β × (neighbors(c) + neighbors(u)) (12)
+ × RSSI(c)

with α = 100 , β = 1 1000 , = 1 10000
Neighborhood discovery state is stopped when the timer expires or when the node receives a FatherOffer message from potential parents. During this state, all nodes collect RSSI of their neighbor nodes and put the nodes with the best value of RSSI in their neighbor table.

Evaluation and Results

In this section, we evaluate our protocol in Linear Wireless Sensors Networks (LWSNs) using Castalia Simulator.

Evaluation of initial network formation

The simulation scenarios use different random topologies with 50 to 300 nodes placed in areas of side 1km to 6km. The characteristics of our topology generator are the following.

• Only connected topologies are built.

• The distance between two nodes is 20m ± 0.2m, so that each node can reach at least two other nodes in a direction (Physical 2-vertex-connectivity for better network resilience).

• During node placement, the probability to create a branch coming out a node is 10% and the deviation angle is between 40 and 140 degrees or -40 and -140 degrees to avoid parallel lines.

• Flickering effect is between -10 and 10 degrees so that the lines are not perfect. To evaluate our protocol in Castalia simulator, we run for each network enumerated above 10.000 experiences with simulation parameters in Table 8. After checking the association rates, we will discuss influence of various factors on total association time (time for all the notes to be associated, including timeouts), either in the case of an initial network formation, or for the case of a new branch addition. Then, we will study the overhead linked to control messages and the distribution of the control messages.

Efficiency

The Figure 9 shows the comparison of association rates of our protocol DiscoProto and DAAM ZigBee for different network sizes. For DAAM, we use parameters Cm = 2, Rm = 2 and Lm = 15. Unlike DAAM, our mechanism associates and provides addresses for almost all nodes of the network node. With DAAM, a lot of orphan nodes appear because of maximum depth limited. Collision problems because of MAC layer in dense linear network prevent message delivery in son collection step and nodes stop sending FatherOffer after three attempts. Therefore, the association rate for the initial network is close to 100% but the few orphans nodes have other chances to get associated. But, even if a node misses the association phase, it will try to join network again during the network join procedure describes in Section 4.1.

Scalability

Figure 10 represents the average association time by network size. We note that the average time increases sublineary: our protocol has the capacity to perform node association in parallel manner which explains why it scales well.

Influence of number of disjunctions

In Figure 11, we show that the average association time decreases with number of disjunctions until a limit value depending on network size. We call disjunction a node in the logical topology with more than one son (usually because a branch has been created by the topology generator). At first, the increasing number of disjunctions allows DiscoProto to run in parallel. But, when the density raises, Figure 12 shows the average number of disjunctions (in the network) with regards to the number of branches actually generated. These figures confirmed the results obtained in [START_REF] Déthié Sarr | Automatic discovery of topologies and addressing for linear wireless sensors networks[END_REF] for the procedure of discovery and topology generation. In a realistic WSN environment, our implemented protocol setups a linear network as well as in the home-made Ruby simulator with poor physical and MAC layer.

Evaluation of the cost of addition of new nodes

We study the cost of association of new nodes to an already formed network using Dynamic DiscoProto. The metrics we are studying are not only the association time but also the cost in terms of number of association messages. The efficiency of the association is measured in terms of association success (percentage of successfully associated nodes). We study initial networks of various sizes (100, 200, 300, 400, 500 nodes), submitted to associations of several new branches (5, 10, 15, 20 and 25 new branches). Each new branch is 20 nodes long. The nodes of the new branch start up at random time after the end of initial network formation. Figure 13 shows that the number of new branches does not affect exponentially the association delay for the same network size.

Influence of initial network size on association time of new branches

Figure 14 shows that the association delay for new branches is not dependent on the network size.

Table 9 shows that the association rate before and after addition of new branches, with an initial network of 200 nodes. Dynamically adding branches does not lower the association rate (it is even better because orphans nodes of the network had time to dynamically associate afterward).

Conclusion

In this paper, we evaluate and improve DiscoProto, an efficient discovery mechanism for Linear WSNs topologies. In addition, we propose Dynamic DiscoProto to support addition of new nodes or new branches to an existing network while keeping the routing nearly hierarchical. The evaluation of our protocol in the Castalia simulator shows that it can be adapted for large Linear WSNs but can cause latency in dense WSNs with many overlapping branches. The implemented routing algorithm supports new address bloc reallocation and ensures efficient data routing for entire network. The next step is to implement a mechanism to detect and repair failing links (both upward and downward). The idea is to define different procedures used by a node to adopt another father without many messages exchanged and to ensure a fast convergence of routing tables. Another aspect we plan to study is to find the optimal number of hops for broadcasting ChallengerOffer message in the neighborhood and, more generally, the impact of different constants (timers of different messages) used in the Sons Collection Step on the association rate and average association time. We also plan to evaluate our protocol with MAC protocol dedicated to Linear WSNs proposed by authors in [START_REF] Hadji | Using a token approach for the MAC layer of linear sensor networks: Impact of the node position on the packet delivery[END_REF].

For the research perspective in the future, we plan to study the deployment of our protocol in real environment (river, pipeline, bridge) using tinyOS nodes [START_REF] Amjad | TinyOSnew trends, comparative views, and supported sensing applications: A review[END_REF] or nodes with Lora [START_REF] Mouhamed | Free space range measurements with semtech lora technology[END_REF] modulation and Mac and Network layers adapted for mesh networks. We could also propose that the pan-coordinator maintains a mapping of DiscoProto logical addresses to IPv6 addresses to use have our nodes as be a part of IoT [START_REF] Hsieh | ScriptIoT: A script framework for and internet-of-things applications[END_REF].

Figure 3 :Figure 4 :

 34 Figure 3: Role assignment in a LT-WSN

Figure 5 :

 5 Figure 5: Example of address distribution

( a )

 a New address block request (b) New address block response (c) Address block propagation

Figure 6 :

 6 Figure 6: Address block reallocation procedure

  (a) Add new branch of nodes (b) Add new branch of nodes (c) Subtree size collection (d) New address block request (e) New address block reallocation (f) Addresses propagation to newly added nodes

Figure 7 :

 7 Figure 7: new branch join procedure

Figure 8 Figure 8 :

 88 Figure 8 shows an example of generate linear topology

Figure 9 :

 9 Figure 9: Average association rate by network size

Figure 10 :

 10 Figure 10: Average association time by network size

Figure 11 :

 11 Figure 11: Average association time by disjunction

Figure 12 :

 12 Figure 12: Average number of disjunction by number of branches

Figure 13 :

 13 Figure 13: Average association delay of addition of new branches

Figure 14 :Table 9 :

 149 Figure 14: Average association time of 400 new nodes by network size

Figure 15

 15 Figure15shows that the number of control messages sent during network building grows up linearly with the number of branches because each potential node responds to each fatherOffer message diffusing in the neighborhood. Association messages are HELLO, FatherOffer, SonOffer, ChallengerOffer, AssociationAccept.

Figure 15 :

 15 Figure 15: Average number of messages with varying number of branches for a network of 200

  The node P that sends the request chooses an AA with the biggest AAC. If multiple nodes with same AAC are found, P selects the node with the highest AA. The neighbor node called LN (Lender Node) that has received a request message calculates the appropriate AA to lend by the following Equation (4), where num BA is the number of borrowed address.

	AA = Addr LN + (Rm -num BA) × Cskip(d) + 1,	(4)
		num BA < Rm
	Figure 2: DIBA mechanism on cluster-tree WSNs
	If the incoming node receives an available address for lending (called AAL), it also reserves the
	address block [AAL, AAL + Cskip(d) -1], where d is the depth of the node lending the address block
	(see Figure 2).
	Routing table of node 32
		Address Lender Borrower
	Parent	0	-	-
	Child 1	33	-	-
	Child 2	48	-	-
	Child 3	56	48	-
	Routing table of node 48
		Address Lender Borrower
	Parent	32	-	-
	Child 1	56	-	32
	Table

Table 3

 3

	).

Table 3 :

 3 Son address table formatThe memory used by the routing table is:

Table 4

 4

		).
	start Block end Block addr Node id Son
	9	11	9	4
	12	26	12	5
	27	32	9	4

Table 4 :

 4 Son address table of node 6 in Figure6c

	start Block end Block addr Node id Son
	6	26	6	3
	27	30	9	3
	33	44	33	12

Table 5 :

 5 Son address table of node 3 in Figure

Table 7 :

 7 CC2410 TX level specifications implemented in Castalia nodes

	Data rate		250 kbps
	Modulation type	PSK
	Bandwidth		20 MHz
	Noise bandwidth	194 MHz
	Noise floor		-100
	Sensitivity		-95
	Power consumption	62 mW
	Table 6: CC2410 RX specifications implemented in Castalia nodes
	Tx (dBm)	0	-1	-3	-5
	Tx (mW) 52.42 55.18 50.69 46.2

Table 8 :

 8 Parameters of simulation

	Simulation time (sec)	1800, 3600, 5200, 7200
	Numbers of runs	10000
	Propagation model	Shadowing (PE=2, 4, sigma=4, 0)
	Tx output power	0 dBm

Low-Power WAN such as LoRaWAN, Sigfox or Weightless could be a good solution to monitor such infrastructures if the nodes are in line-of-sight with the gateways, which is not the case for tunnels, rivers in canyons..

In this article, we consider LWSNs as mostly Linear Wireless Sensors Networks, WSNs constituted of long branches of nodes connected together.

The distributed algorithm is a relaxation of the non-distributed Prim-Jarnik algorithm and best costs are only compared in a portion of the graph.


est destinée au dépôt

Second-Order Statistics of Complex Signals

The pdf copy of the final published text can be obtained from to the author at the folowing address : bernard.picinbono@lss.supelec.fr Abstract-The second-order statistical properties of complex signals are usually characterized by the covariance function. However, this is not sufficient for a complete second-order description, and it is necessary to introduce another moment called the relation function. Its properties, and especially the conditions that it must satisfy, are analyzed both for stationary and nonstationary signals. This leads to a new perspective concerning the concept of complex white noise as well as the modeling of any signal as the output of a linear system driven by a white noise. Finally, this is applied to complex autoregressive signals, and it is shown that the classical prediction problem must be reformulated when the relation function is taken into consideration.

  Bernard Picinbono, Pascal Bondon To cite this version: Bernard Picinbono, Pascal Bondon. Second-Order Statistics of Complex Signals. IEEE Transactions on Signal Processing, 1997, IEEE Transactions on Signal Processing, 45 (2), pp.411 -420. hal-01661242 Second-Order Statistics of Complex Signals Bernard Picinbono, Fellow, IEEE and Pascal Bondon, Member, IEEE