André Lieutier

Boris Thibert

Geodesic as limits of geodesics on PL-surfaces

Keywords: subdivision surfaces, triangulations, PL-surfaces, geodesics, shortest paths, convergence

published or not. The documents may come

Geodesics as limits of geodesics on PL-surfaces

Introduction

A geodesic is usually defined as a curve on a surface that is locally a shortest path. Geodesics appear naturally in several applications, among which we can mention: i) The modelling of the human heart: the heart left ventricle can be modelled by a family of embedded surfaces; a muscular fiber of the central region of the left ventricle has 1 particular properties and can be considered as a geodesic of one of those surfaces [START_REF] Mourad | Description topologique de l'architecture fibreuse et modélisation mécanique du myocarde[END_REF][START_REF] Streeter | Jr Gross morphology and fiber geometry of the heart[END_REF]. ii) In the fabrication of composite parts by filament winding, the filament must idealy wind along geodesics [START_REF] Blais | Modélisation de déviateurs CRT : empilement du fil et tassage[END_REF]. iii) Finally the computation of radar cross sections involves the simulation of creeping ray which follow geodesics of the object [START_REF] Bouche | Méthodes asymptotiques en électromagnétisme (French) Mathématiques & Applications[END_REF][START_REF] Dessarce | Calculs par lancer de rayons[END_REF]. In this context, and since piecewise linear 2-manifolds (denoted by PL-surfaces in the following) are widely used for surface modelling, it is natural to consider the modelling of geodesics on surfaces and in particular on PL-surfaces.

We distinguish the geodesics from the more restricted class of shortest paths. A shortest path is a curve on a surface that is connecting two points and whose length does not decrease if it is perturbed (without moving the two extremities). A geodesic is a curve on a surface whose length does not decrease if it is pertubed in a small neighborhood of any point. A shortest path is clearly a geodesic, but the converse is not true (for example, a great circle is a geodesic but not a shortest path of the sphere).

There exist several algorithms that build shortest paths on PLsurfaces [START_REF] Kimmel | Finding shortest paths on surfaces[END_REF][START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF][START_REF] Maekawa | Computation of shortest paths on free-form parametric surfaces[END_REF][START_REF] Peyré | Geodesic Computations for Fast and Accurate Suface Remeshing and Parametrization[END_REF]. Concerning the geodesics, Pham-Trong and her coauthors have also proposed an algorithm that builds geodesics on PL-surfaces [START_REF] Pham-Trong | Pseudo-geodesics on three-dimensional surfaces and pseudo-geodesic meshes[END_REF]. In particular, they have also considered a sequence (T n ) n∈N of PL-surfaces defined by the De Casteljau subdivision for Bezier surfaces that is converging to a Bezier surface S. On each T n , they build a geodesic C n whose sequence converges to a curve C. The natural question is then to wonder whether C is a geodesic of S.

The convergence of geodesics has already been studied in the case of shortest paths by Hildebrandt et al. [START_REF] Hildebrandt | On the Convergence of Metric and Geometric Properties of Polyedral Surfaces[END_REF] and Memoli et al. [START_REF] Memoli | Distance functions and geodesics on submanifolds of R d and point clouds[END_REF]. They show that if a sequence (T n ) n∈N of PL-surfaces converges in Hausdorff distance to S, if the normals of T n also converge to the normals of S, then the limit curve of a sequence of shortest paths is a shortest path of S. However, this result does not hold anymore for geodesics: we provide in this paper (Section 4) a sequence (T n ) n∈N of PL-surfaces whose both distance and angular limit is a plane. However a sequence of geodesics C n ⊂ T n converges toward a limit curve C which is not a straight line of the plane, and thus not a geodesic of the plane.

It is worth noting that the result of convergence of Hildebrandt et al. [START_REF] Hildebrandt | On the Convergence of Metric and Geometric Properties of Polyedral Surfaces[END_REF] cannot be used in some applications: for example, in the modelling of the human heart, the curves modelling the fibers are closed and are not shortest paths [START_REF] Mourad | Description topologique de l'architecture fibreuse et modélisation mécanique du myocarde[END_REF]. Furthermore, this result cannot be used to validate the algorithm given in [START_REF] Pham-Trong | Pseudo-geodesics on three-dimensional surfaces and pseudo-geodesic meshes[END_REF]: indeed, Pham Trong and her coauthors build a sequence of geodesics that are not shortest paths in general.

The main result of this paper deals with convergence for geodesics. More precisely, we suppose, as for the result with shortest paths given in [START_REF] Hildebrandt | On the Convergence of Metric and Geometric Properties of Polyedral Surfaces[END_REF], that the sequence (T n ) n∈N of PL-surfaces converges in Hausdorff distance and in normals to a smooth surface S. In addition, we also suppose that there exist two constants K 1 and K 2 independant on n such that the length of the edges of T n is greater than K1 2 n and the maximal angle between the normals of T n and the normals of S is less than K2 2 n . In other words, the rate of convergence of the length of the edges cannot be faster than the rate of convergence of the normals.

This result is then applied to PL-surfaces T n that follow subdivision schemes, such as for example schemes for bicubic B-splines, or Catmull-Clark schemes, or Bezier surfaces. In particular, our results validate the algorithm of Pham-Trong and her coauthors [START_REF] Pham-Trong | Pseudo-geodesics on three-dimensional surfaces and pseudo-geodesic meshes[END_REF] that builds geodesics on subdivision surfaces.

It is interesting to note that shortest paths and geodesics do not deal with notions of the same order: the notion of shortest path relies on the notion of length, which is a quantity related to the first derivative. However, since a geodesic is defined locally, it depends on the infinitesimal variation of the length, which is a notion of second order.

In this paper, we focus on the problem of convergence of geodesics. In Section 2, we recall the main definitions. In Section 3, we recall the result of convergence of [START_REF] Hildebrandt | On the Convergence of Metric and Geometric Properties of Polyedral Surfaces[END_REF][START_REF] Memoli | Distance functions and geodesics on submanifolds of R d and point clouds[END_REF] for shortest paths. In Section 4 we give a counter example showing that the situation is more complicated for geodesics, and we also give the main result of convergence. We show in Section 5 that this result can be applied to several subdivision schemes (and in particular to the algorithm proposed in [START_REF] Pham-Trong | Pseudo-geodesics on three-dimensional surfaces and pseudo-geodesic meshes[END_REF]). The last section proves all these results. In the following, we will refer to triangulations instead of PL-surfaces.

Definitions

Smooth surfaces

In the following, a smooth surface means a C 2 surface which is regular, oriented, compact with or without boundary. We have the following proposition [START_REF] Federer | Curvature measures[END_REF] Proposition 1 Let S be a smooth compact surface of R 3 . Then there exists an open set U S of R 3 containing S and a continuous map ξ from U S onto S satisfying the following: if p belongs to U S , then there exists a unique point ξ(p) realizing the distance from p to S (ξ is nothing but the orthogonal projection onto S).

This proposition allows to introduce the following notion introduced by H. Federer [START_REF] Federer | Curvature measures[END_REF]: The reach of a surface S is the largest r > 0 for which ξ is defined on the open tubular neighborhood U r (S) of radius r of S.

Triangulations

A triangulation T is a connected topological 2-manifold made of a finite union of triangles of R 3 , such that the intersection of two triangles is either empty, or equal to a vertex, or equal to an edge.

Curves

In the following a curve C means a lipschitz parametrized curve C : [0, 1] → R 3 . Its length is denoted by l(C). Similarly, for 0 ≤ t a < t b ≤ 1 we denote by l(C, t a , t b ) the length of the curve C restricted to [t a , t b ]. As a particular case of Rademacher theorem, we know that C is differentiable almost everywhere. Moreover, it is the integral of its derivative. Whenever it exists, we denote by C ′ (t) the derivative of C at t. • We say that C has a uniform parametrization if it satisfies for almost every 3 with uniform parametrization is said to be a geodesic of a lipschitz surface M (M can be a triangulation or a smooth surface) if it locally minimizes the length, i.e. if for every t ∈ [0, 1], there exists 0 ≤ t a ≤ t ≤ t b ≤ 1 (where t a < t if t > 0 and

t ∈ [0, 1] C ′ (t) = l(C). • A curve C : [0, 1] → M ⊂ R
t b > t if t < 1), such that any lipschitz curve C : [0, 1] → M such that C(0) = C(t a ) and C(1) = C(t b ) satisfies l(C, t a , t b ) ≤ l( C).
The geodesic is said to be interior if for every t ∈ [0, 1], C(t) is interior to the surface M . The geodesic C is a shortest path, if the length of any curve on M connecting C(0) and C(1) has a length greater than l(C).

Properties of geodesics on triangulations

Let C be a polygonal curve of a triangulation T . Then C is a geodesic of T if and only if:

• C it is a straight line on each triangle. (The vertices p of C then belongs to the edges of T or are vertices of T .)

• If p belongs to the interior of an edge of T , then the incident and refracted angles of C at p are equal (see Figure 1).

• If p is a vertex of T , then C separates the set of the triangles of T containing p into two connected regions r 1 and r 2 (see Figure 1). If one denotes by α r1 p the sum of the angles α r1 i of the triangles of region r 1 at p (resp. by α r2 p the angles α r2 i of the triangles of region r 2 at p), one has

α r1 p ≥ π and α r2 p ≥ π. (1)
Remark that if a geodesic traverses a vertex p of T , then the sum of the angles of the triangles of T at the vertex p is greater than 2π: α r1 p + α r2 p ≥ 2π. We can also notice that if α r1 p + α r2 p > 2π, then the geodesic containing the vertices q and p is not unique (see Figure 1): two distinct polygonal curves C 1 and C 2 containing p and q that satisfy Equation (1) are geodesics.

Convergence of shortest paths

In this section, we recall a positive result of convergence for shortest paths, that can be found in [START_REF] Hildebrandt | On the Convergence of Metric and Geometric Properties of Polyedral Surfaces[END_REF] (a similar result can also be found in

p q r α r 1 1 α r 1 2 α r 1 m C region r 1
Figure 1: A geodesic C passing through a vertex of a 3D triangulation [START_REF] Memoli | Distance functions and geodesics on submanifolds of R d and point clouds[END_REF]). Proposition 2 states that if a sequence (T n ) n∈N of triangulations tends to a surface S in distance and in normals, then a sequence of shortest paths of T n is converging to a shortest path of S. Although this proposition has already been proved, for the sake of completeness, we give here a proof.

Proposition 2 Let (T n ) n∈N be sequence of triangulations that converges in the Hausdorff sense to a smooth C 2 surface S. Let (d n ) n∈N be sequence of real numbers converging to 0. We suppose that a) for every n, the restriction ξ n : T n → S of the map ξ to T n is bijective;

b) for every m ∈ T n , the angle between any triangle ∆ containing m and the tangent plane Π S ξ(m) of S at ξ(m) is smaller than d n ; Let a and b be two points of S, and C n : [0, 1] → R 3 be a shortest path on T n between ξ -1 n (a) and ξ -1 n (b), with uniform parametrization. Then • there exists a subsequence (C n k ) of (C n ) that is uniformly converging to a curve C;

• if a subsequence of (C n ) is converging to a limit curve C, then C is a shortest path of S.

Proof of Proposition 2

First remark that the convergence in distance and in normals imply that there exists a sequence (ǫ n ) n∈N converging to 0 such that for every curve γ n of T n we have (see [START_REF] Morvan | Unfolding Of Surfaces[END_REF] or [START_REF] Hildebrandt | On the Convergence of Metric and Geometric Properties of Polyedral Surfaces[END_REF]):

(1 -ǫ n )l(γ n ) ≤ l(ξ n • γ n ) ≤ (1 + ǫ n )l(γ n ).
Let us denote by C a shortest path between a and b on S. We then have :

l(ξ n • C n ) ≤ (1 + ǫ n )l(C n ) ≤ (1 + ǫ n )l(ξ -1 n • C) ≤ 1 + ǫ n 1 -ǫ n l( C),
which finally gives

l( C) ≤ l(ξ n • C n ) ≤ 1 + ǫ n 1 -ǫ n l( C). (2)
In particular l(ξ n • C n ) converges to l( C) and is then bounded. Then there exists a constant k independant of n such that for every t, t ′ ∈ (0, 1)

ξ n • C n (t) -ξ n • C n (t ′ ) ≤ l(ξ n • C n , t, t ′ ) ≤ (1 + ǫ n ) l(C n , t, t ′ ) = (1 + ǫ n ) l(C n )|t -t ′ | ≤ k |t -t ′ |.
The set {ξ n • C n , n ∈ N} is then equicontinuous (see [START_REF] Dieudonné | Éléments d'analyse. Tome I: Fondements de l'analyse moderne[END_REF] for example for more details). Arzela-Ascoli Theorem [START_REF] Dieudonné | Éléments d'analyse. Tome I: Fondements de l'analyse moderne[END_REF] then implies that a subsequence (ξ

n k • C n k ) k∈N of (ξ n • C n ) n∈N is
(C) ≤ lim inf l(C n ). However, Equation (2) implies that lim C n = lim ξ n • C n = l( C). We then have l(C) ≤ l( C) which means that C is a shortest path.

Convergence of geodesics

The main result of this paper is given in this section. We first show in Section 4.1, by building a counter-example, that the previous result of convergence for shortest paths, does not hold anymore for geodesics in general. Then, in Section 4.2, we give a general result of convergence for geodesics.

Counter-example

The counter example shows a sequence (T n ) n∈N of triangulations whose both distance and angular limit is a plane. However a sequence of geodesics C n ⊂ T n converges toward a limit curve C which is not a straight line of the plane, and thus not a geodesic of the plane. We can notice that the triangulations T n and the curves C n satisfy all the assumptions of the result of Proposition 2, except that C n is a geodesic (and not a shortest path). However the limit curve C is not a geodesic. This counter example implies that the convergence in distance and in normals of T n to S is not sufficient to expect a result of convergence for geodesic. Detail of the construction: Figure 2 shows the triangulation T n for some n. The triangulation overlaps with the horizontal plane Π H of equation z = 0 outside the large circle and inside the small one. In the ring between the two circles, it is made of 4 n identical small "roof shaped" bumps detailed on the right of Figure 2

∈ N, β 1 + β 2 + β 3 + β 4 > π.
Observe that the sequence (T n ) n∈N of triangulations converges toward the plane in the Hausdorff sense. Furthermore the normals of T n tend to the normals of the plane. The shortest path on T n between the point a n and b n is the straight line (dotted line on Figure 2). However, the line a n m n n n b n , wraped around the circle between m n and n n is a local minimum, that is a geodesic, between a n and b n (see Section 2). These geodesics converge in the Hausdorff sense toward a curve C composed of two lines and an arc of circle.

Convergence toward a geodesic

The main result of this paper is the following theorem. It states that if a sequence (T n ) n∈N of triangulations is converging to a smooth surface S, then, under reasonable assumptions, a sequence of geodesics of T n is converging to a geodesic of S.

Theorem 1 Let S be a smooth surface of R 3 , r denote the reach of S, and (T n ) be a sequence of triangulations. Let K, K, θ min be positive constants and (d n ) a sequence of real numbers converging to 0. Suppose that for every n: The proof of this theorem is given in Section 6.

Remark 1 If a sequence (T n ) n∈N of triangulations satisfies the conditions of Theorem 1, then the valency of the vertices of T n is uniformly bounded by a constant independant on n. Indeed, the convergence of the normals implies that the sum of the angles at a vertex p n tends to 2π when n tends to infinity. Combined with the fact that all the angles are greater than θ min , that implies that the valency of p n is uniformly bounded.

Remark 2 Remark that the result of Theorem 1 still holds if the sequence (C n ) n∈N does not converge to a curve C, but if we suppose that (l(C n )) n∈N is bounded. Indeed, in that case, we can show (as in the proof of Proposition 2) that the family {C n , n ∈ N} is equicontinuous. The Arzela-Ascoli theorem (see [START_REF] Dieudonné | Éléments d'analyse. Tome I: Fondements de l'analyse moderne[END_REF] for example) then implies that a subsequence (C n k ) k∈N of (C n ) n∈N uniformly converges to a curve C. Theorem 1 can then be applied to (C n k ) k∈N .

Application to subdivision surfaces

The previous theorem can be easily applied to subdivision surfaces. In this section, we first give a general corollary, Corollary 1, that can be easily applied to several subdivision schemes. We then show that the result of convergence for geodesics holds when the triangulations follow a subdivision scheme for either bicubic B-splines, or Catmull-Clark schemes, or Bezier surfaces. We first need to give a few definitions.

Let (P n ) n∈N be a sequence of parametrized triangulations

P n : [0, 1] 2 → R 3 that is converging to a paramatrized smooth surface f : [0, 1] 2 → R 3 . The parameter domain [0, 1] 2 of each P n can be triangulated so that P n is linear on each triangle of [0, 1] 2 .
-We say that the parameter domain of P n is a triangular grid if its vertices are p i,j n = i 2 n , j 2 n (where i, j ∈ {0, ..2 n }) and the edges are p i,j n p i+1,j n , p i,j n p i,j+1 n and p i,j n p i+1,j+1 n . -We say that (P n ) n∈N uniformly converges to a function f with rate of convergence 1 2 n if:

∃N ∈ N, ∃K ∈ R, n > N ⇒ sup (u,v)∈[0,1] 2 P n (u, v) -f (u, v) ≤ K 2 n .
-We say that (P n ) n∈N uniformly converges in derivative to f with rate of convergence 1 2 n if there exists K > 0 and N ∈ N, such that for any n > N : sup i∈{0,..,2 n -1} j∈{0,..,2 n }

2 n P n i + 1 2 n , j 2 n -P n i 2 n , j 2 n - ∂f ∂u i 2 n , j 2 n ≤ K 2 n ,
and sup i∈{0,..,2 n } j∈{0,..,2 n -1}

2 n P n i 2 n , j + 1 2 n -P n i 2 n , j 2 n - ∂f ∂v i 2 n , j 2 n ≤ K 2 n .
Corollary 1 Let (P n ) n∈N be a sequence of parametrized triangulations

P n : [0, 1] 2 → R 3 and f : [0, 1] 2 → R 3 be a parametrized surface of class C 2 , such that: a) the parameter domain of each P n is a triangular grid, b) (P n ) n∈N uniformly converges to f with rate of convergence 1 2 n , c) (P n ) n∈N uniformly converges in derivative to f with rate of con- vergence 1 2 n , d) f is regular, i.e. ∀(u, v) ∈ [0, 1] 2 ∂f ∂u (u, v) ∧ ∂f ∂v (u, v) = 0. Let (C n ) n∈N be a sequence of polygonal curves C n : [0, 1] → R 3 with uniform parametrization such that C n is an interior geodesic of P n . If (C n ) n∈N converges in the sup norm sense toward a curve C which is interior to S, then C is of class C 2 and is a geodesic of S.
The proof of this corollary is given in Section 6.7.

The previous corollary can be easily applied to some subdivision schemes. As an example, we give the following corollary concerning subdivision scheme for bicubic B-spline (see for example [START_REF] Farin | Curves and surfaces for computer aided geometric design, a practical guide[END_REF] for details on B-splines). First remark that the subdivision scheme for bicubic B-splines generates a sequence of quadrangulations. Each quadrangulation can be considered as a triangulation by dividing each quadrangle into two triangles.

Corollary 2 Let (P n ) n∈N be a sequence of triangulations (or quadrangulations) defined by the subdivision scheme for bicubic B-spline that is converging to a regular B-spline f (i.e. satisfies assumption d) of Corollary 1).

Let (C n ) n∈N be a sequence of polygonal curves C n : [0, 1] → R 3 with uniform parametrization such that C n is an interior geodesic of P n . If (C n ) n∈N converges in the sup norm sense toward a curve C which is interior to S, then C is of class C 2 and is a geodesic of S.
Proof The assumptions a) and d) of Corollary 1 are clearly satisfied. Let u ∈ [0, 1]. The polygonal curve v → P n (u, v) follows the subdivision scheme for cubic B-spline and is uniformly converging to the function f u : v → f (u, v). More precisely, there exists K u ∈ R such that for every n ∈ N one has (see Theorem 4.12 of [START_REF] Dyn | Subdivision schemes in geometric modelling[END_REF] or Corollary 3.3 of [START_REF] Dyn | Subdivision schemes in CAGD[END_REF]):

sup v∈[0,1] P n (u, v) -f (u, v) ≤ K u 2 n .
In fact, one can show that K u does not depend on u: let us denote by P i,j = P 0 (u i , v j ) the poles of the initial control mesh (where i, j ∈ {0, ..M }). The k th pole P u,k of f u is the evaluation at u of the B-spline whose control net is P 0,k ,...P M,k . That implies that P u,k belongs to the convex hull of P 0,k ,...P M,k and then to the convex hull of the P i,j .

Since K u only depends on the maximum distance between two poles of the control net of f u (see the proof of Theorem 4.11 of [START_REF] Dyn | Subdivision schemes in geometric modelling[END_REF]), it only depends on the diameter of the convex hull of the P i,j and does not depend on u.

This implies that (P n ) n∈N uniformly converges toward f . Assumption b) is then proved. Similarly, assumption c) is also proved on the divided difference scheme (see for example [START_REF] Warren | Subdivision Methods for Geometric Design: A constructive Approach[END_REF] for details on divided difference schemes).

Remark 3 Corollary 2 directly implies that the result of convergence still holds for Catmull-Clark schemes if the limit curve C does not contain extraordinary points of S (see [START_REF] Catmull | Recursively generated B-spline surfaces on arbitrary topological surfaces[END_REF] for details on the Catmull-Clark scheme). By extraordinary point on the limit surface, we precisely mean the limit of the sequence of vertices corresponding to an extraordinary point of the triangulation through successive subdivisions. Indeed, after a sufficient number of iterations, the curve traverses only a finite number of bicubic B-splines patches where Corollary 2 can be applied.

More precisely, the curve C is at a distance greater than µ > 0 from all the extraordinary points. Let us denote by (P n ) n∈N the sequence of triangulations (or quadrangulations) defined by Catmull-Clark scheme. Then, by compacity, there exists a finite number of triangulations V 1 n ,...V p n and 0 = t 0 < t 1 < t 2 < ... < t p = 1 such that: i) each V i n follows a subdivision scheme for bicubic B-spline (as in Corollary 2) and is converging to

S i ⊂ S; ii) for every i, C n ([t i-1 , t i ]) ⊂ V i n and C([t i-1 , t i ]) ⊂ S i ; iii) V i n ∩ V i+1 n ∩ C n ([t i-1 , t i ]
) is homeomorphic to a connected curve. We then apply Corollary 2 on each V i n .

Remark 4 A proof similar to the one of Corollary 2 shows that this result also holds for Bezier surfaces and their successive control nets defined by the De Casteljau algorithm.

Proof of Theorem 1

The aim of this section is to prove Theorem 1. The proof being quite long, we first give an overview of each subsection.

6.1 We introduce the definitions needed for the proof, but not usefull for the statement of Theorem 1.

6.2 We give basic lemmas about angles and orthogonal projections onto planes. These lemmas are used in the rest of the proof.

6.3

We give results concerning geodesics on triangulations. We first bound from above the angle of deviation β 3D dev (p n ) of the geodesic at a vertex p n . Intuitively, this means that if the triangulation is "almost planar", then the geodesic is "not turning too much". We then also bound from above the angle of deviation β T S dev (p n ) of the projection of the geodesic onto a plane. Intuitively, if the plane is "almost parallel" to the triangles of the triangulation, then the projection of the geodesic is "turning much less" than the geodesic itself.

6. [START_REF] Clarke | Optimization and NonSmooth Analysis[END_REF] We bound the number of intersections of a polygonal curve C n of T n with the edges of a triangulation T n . Intuitively, we show that if the curve is "not too long"and if it does "not turn too much", then its number of intersection is "small".

6.5 There is no underlying triangulation in this section. The result gives a sufficient condition on a sequence of curves so that its limit curve is of class C 1,1 .

6.6

We give here the core of the proof of Theorem 1, that is using the previous sections.

6.7

We prove in this section Corollary 1.

Preliminary definitions

Let ǫ be smaller than the reach of S. Let T be a triangulation such that ξ induces an injection from T to S. Let denote by R(T ) the set of polygonal curves C of T that are linear on each triangle of T , to be more precise, if τ is a triangle (a triangle is defined here as a closed simplex, i.e. containing its boundary edges and vertices) of T , the image of each connected component of {t ∈ [0, 1] , C(t) ∈ τ } is a line segment: geodesics on T trivially satisfy this condition. Notice that this condition allows the curve to visit more than once a given triangle τ but, in this case, has to visit the interior of other triangle between two successive visits of τ . Let C ∈ R(T ) be a polygonal curve that belongs to the tubular neighborhood V ǫ (S) of radius ǫ of S. In the following, if m ∈ S, we denote by P S m the orthogonal projection onto the tangent plane of S at the point m.

• The total curvature of C is given by:

T C 3D (C) = p vertex of C β 3D dev (p),
where β 3D dev (p) is the deviation angle of C at the vertex p (see Figure 3). Similarly, for 0 ≤ t a < t b ≤ 1 we denote by T C 3D (C, t a , t b ) the total curvature of the curve C restricted to [t a , t b ].

• The tangent total curvature of C with respect to S is defined by

T C S T an (C) = p vertex of C β T S dev (p),
where β T S dev (p) is the deviation angle of P S ξ(p) (C) at the vertex ξ(p) (see Figure 3). • Let ♯C be the number of intersections between C and the edges of T . More precisely, if one denotes by E the set of edges of T and N CC (X) the number of connected components of a set X:

♯C = e∈E N CC (C([0, 1]) ∩ e)
Notice that each time C transversally crosses an edge away from a vertex, ♯C is increased of 1 and each time C crosses a vertex "generically" (that is without following an edge), ♯C is increased of the vertex valency. If it follows an edge from one vertex to the other, it crosses two vertices but the edge is counted once. Similarly, for 0 ≤ t a < t b ≤ 1 we denote by ♯(C, t a , t b ) the number of intersections between the curve C restricted to [t a , t b ] and the edges of T .

Basic lemmas about planes in R 3

In this section, we prove several very usefull basic lemmas.

Lemma 1 There exists K > 0 such that for every planes Π and Π 1 and for every vectors u and v of Π, one has:

|∠(u, v) -∠(P 1 (u), P 1 (v))| ≤ K ∠(u, v) α 2 ,
and also

|∠(u, v) -∠(P 1 (u), P 1 (v))| ≤ K ∠(P 1 (u), P 1 (v)) α 2 ,
where P 1 is the orthogonal projection onto Π 1 , and α is the angle between Π and Π 1 .

Proof

• We put θ = ∠(u, v) and θ 1 = ∠(P 1 (u), P 1 (v)). We clearly have

P 1 (u) ≤ u ≤ 1 cos α P 1 (u) .
The same inequality holds with v. Furthermore a simple calculus gives

Area(∆) = 1 cos α Area(P 1 (∆)).
Now by using the fact that

sin θ = Area(∆) 2 u v and sin θ 1 = Area(P 1 (∆)) 2 P 1 (u) P 1 (v) , we have cos α ≤ sin θ 1 sin θ ≤ 1 cos α ,
and then

| sin θ -sin θ 1 | = O(α 2 ) sin θ = O(α 2 ) θ. (3)
• We put X = u u and Y = v v . Then X -P 1 (X) ≤ sin α ≤ α. Furthermore we have π 2 -α ≤ ∠ (P 1 (X) -X, Y ) ≤ π 2 + α, which implies that | < X -P 1 (X), Y > | ≤ α 2 .
Similarly, we also have

| < Y -P 1 (Y ), X > | ≤ α 2 which implies that P 1 (X) P 1 (Y ) cos θ 1 =< P 1 (X), P 1 (Y ) > =< X, Y > +O(α 2 ) = cos θ + O(α 2 ).
We then have:

| cos θ -cos θ 1 | = O(α 2 ). (4)
• Suppose now that θ ∈ 0, π 4 . If α is small enough, then Equation (3) implies that θ p ∈ 0, 3π 8 and that

|θ 1 -θ| ≤   sup x∈[0,sin 3π 8 ] | arcsin ′ (x)|   | sin θ 1 -sin θ| = O(α 2 ) θ.
Remark that the same results holds if θ ∈ 3π 4 , π 2 .

• Suppose now that θ ∈ π 4 , 3π 4 . If α is small enough, then Equation (4) implies that

|θ 1 -θ| ≤   sup x∈[0,cos π 8 ] | arccos ′ (x)|   | cos θ-cos θ 1 | = O(α 2 ) = O(α 2 ) θ.
The second inequation of Lemma 1 is a direct consequence of the first one (with a larger constant).

2 sin ∠(u, v) 2 ≤ u -v min( u , v ) ≤ 2 sin ∠(u, v) 2 + 2ǫ 1 -ǫ .
Proof Suppose that v ≤ u . We then have

u -v v ≥ a v = 2 sin ∠(u, v) 2 .
We also have

u -v v ≤ a + b v = 2 sin ∠(u, v) 2 + u -v v ≤ 2 sin ∠(u, v) 2 + 2ǫ 1 -ǫ .
Lemma 3 For some constant K, if P 1 and P 2 are the respective projections on two planes Π 1 and Π 2 with ∠ (Π 1 , Π 2 ) = θ ≤ 1 10 and if u and v are two unit vectors and γ > 0 a number such that:

• ∠ (Π 2 , u) ≤ 1 10 • ∠ (Π 2 , v) ≤ 1 10 • ∠ (Π 1 , u) ≤ γ ≤ 1 10 • ∠ (Π 1 , v) ≤ γ ≤ 1
10 , one has:

∠ (P 2 (u), P 2 (v)) ≤ K ∠ (P 1 (u)), P 1 (v)) + sin θ ∠ (u, v) + γ 2 .
Proof Let δ = vu The affine projection P 1 induces a corresponding projection between vectors which is also denoted P 1 . δ can be splited in δ = P 1 (δ) + δ with P 1 (δ) and δ respectively parallel and orthogonal to Π 1 , which entails P 1 (δ) = P 1 (δ 1 ) and δ ≤ δ . One has P 2 (δ) = P 2 (P 1 (δ)) + P 2 ( δ) and:

P 2 (δ) ≤ P 2 (P 1 (δ)) + P 2 ( δ) ≤ P 1 (δ) + P 2 ( δ) ≤ P 1 (δ) + sin θ δ ≤ P 1 (δ) + sin θ δ .
From ∠ (Π 2 , u) ≤ 1 10 and ∠ (Π 2 , v) ≤ 1 10 , one has:

min ( P 2 (u) , P 2 (v) ) ≥ cos 1 10
.

Using twice Lemma 2 we have:

2 sin ∠ (P 2 (u), P 2 (v)) 2 ≤ 1 cos 1 10 P 2 (δ)
and

P 1 (δ) ≤ 2 sin ∠ (P 1 (u), P 1 (v)) 2 + O(γ 2 ).
This gives:

2 sin ∠ (P 2 (u), P 2 (v)) 2 ≤ 1 cos 1 10 2 sin ∠ (P 1 (u), P 1 (v)) 2 + O(γ 2 ) + sin θ δ .
Using that δ = 2 sin ∠(u,v)

2
and that, for any angle β ∈ [0, π] one has 2 π β ≤ 2 sin β 2 ≤ β we get the above bound on ∠ (P 2 (u), P 2 (v)).

Lemma 4 There exists K > 0 such that for every planes Π 1 and Π 2 and for every unit vectors u and v of R 3 such that ∠(u, Π i ) ≤ π 4 and ∠(v, Π i ) ≤ π 4 , we have:

|∠(P 1 (u), P 1 (v)) -∠(P 2 (u), P 2 (v))| ≤ K γ,
where P 1 and P 2 denote the respective projections on Π 1 and Π 2 , and γ denotes the angle between Π 1 and Π 2 .

Proof We put P = P 1 -P 2 , θ = ∠(u, v), θ 1 = ∠(P 1 (u), P 1 (v)) and

θ 2 = ∠(P 2 (u), P 2 (v)). Case 1: θ ∈ π 10 , π -π 10
Let X be a vector of R 3 . We put X 1 = P 1 (X) and X ⊥ 1 = X -X 1 . We then have P (X) = P (X 1 ) + P (X ⊥ 1 ) ≤ P (X 1 ) + P (X ⊥ 1 ) .

On the other hand

P (X 1 ) = X 1 -P 2 (X 1 ) ≤ sin γ X 1 ≤ γ X , and
P (X ⊥ 1 ) = P 2 (X ⊥ 1 ) ≤ sin γ X ⊥ 1 ≤ γ X ,
which implies that P (X) ≤ 2 γ X . We then have:

| P 1 (u) -P 2 (u) | ≤ P (u) ≤ 2 γ and | P 1 (v) -P 2 (v) ≤ 2 γ.
Furthermore, the fact that

P i (u) ≥ 1 √ 2 , P i (v) ≥ 1 √ 2 and cos θ i = P i (u) 2 + P i (v) 2 -P i (u -v) 2 2 P i (u) P i (v) ,
implies that there exists K > 0 such that

| cos θ 1 -cos θ 2 | ≤ K γ.
Then there exists K > 0 such that

|θ 1 -θ 2 | ≤ K γ. Case 2: θ ∈ 0, π 10
We denote by ∆ a triangle with edges u, v and uv and by θ ′ and θ ′′ the two other angles of ∆. Remark that θ ′ = θ ′′ ∈ π 2 -π 20 , π 2 . We also denote by θ ′ i and θ ′′ i the two angles of the triangle P i (∆). We then

have |θ ′ 1 -θ ′ 2 | ≤ K γ and |θ ′′ 1 -θ ′′ 2 | ≤ K γ. Then |θ 1 -θ 2 | = |(π -θ ′ 1 -θ ′′ 1 ) -(π -θ ′ 2 -θ ′′ 2 )| ≤ 2 K γ. Case 3: θ ∈ π -π 10 , π Since π -θ ∈ 0, π 10 , we have |θ 1 -θ 2 | = |(π -θ 1 ) -(π -θ 2 )| ≤ K γ.

Majoration of the deviation angles of a geodesic

Proposition 3 There exists K 1 , such that for every n: if C n is a geodesic of T n and p n is a vertex of C n , then we have:

β 3D dev (p n ) ≤ K 1 α n , and β T S dev (p n ) ≤ K 1 α 2 n ,
where α n is the maximal angle between all the triangles of T n containing p n and Π S ξ(pn) . Proof Case 1: p n is not a vertex of T n We denote by ∆ 1 and ∆ 2 the two triangles containing p n and by e their common edge. We denote by Π 1 and Π 2 the planes containing respectively ∆ 1 and ∆ 2 . We consider the following unit vectors:

pn angle α Π 1 ∆ 2 ∆ 1 - → V -→ V 1 -→ V 2 -→ U 1 -→ U 2 C n β β
-→ V is colinear to e; -→ U 1 and -→ U 2 are colinear to C n respectively in the planes Π 1 and Π 2 and oriented with the orientation of curve C n ;

-→ V 1 is the vector in the plane Π 1 orthogonal to -→ V ; -→ V 2 is the vector in the plane Π 2 orthogonal to -→ V and -→ n 1 is the vector normal to P 1 . We denote by α the angle between ∆ 1 and ∆ 2 , and β is the incident and refracted angle of C n at p n (see Figure 5). One has:

-→ V 2 = cos α -→ V 1 + sin α -→ n 1 and -→ U 1 = cos β -→ V + sin β -→ V 1 -→ U 2 = cos β -→ V + sin β -→ V 2 .
Which gives:

-→ U 2 - -→ U 1 = sin β (cos α -1) -→ V 1 + sin α -→ n 1 then 2 sin β 3D dev (p) 2 = -→ U 2 - -→ U 1 ≤ 1 2 α 2 + α ≤ 2 α ≤ 2α n . (5)
Then, there exists K 1 such that

β 3D dev (p) ≤ K 1 α n .
Similarly, one has, P S ξ(pn) being linear:

P S ξ(pn) -→ U 2 -P S ξ(pn) -→ U 1 = sin β (1 -cos α)P S ξ(pn) -→ V 1 + sin αP S ξ(pn) ( -→ n 1 ) .
That gives :

P S ξ(pn) -→ U 2 -P S ξ(pn) -→ U 1 ≤ sin β (1 -cos α) P S ξ(pn) -→ V 1 + sin α P S ξ(pn) ( -→ n 1 )
We know that:

P S ξ(pn) ( -→ n 1 ) ≤ sin α n ≤ α n .
Which gives:

P S ξ(pn) -→ U 2 -P S ξ(pn) -→ U 1 ≤ 1 2 α 2 + α α n ≤ 2α n 2 . ( 6
)
Lemma 2 then implies:

β T S dev (p n ) ≤ K 1 α 2 n .
Case 2: p n is a vertex of T n • In this case, the sum of all the angles of T n at p n is necessarily greater

p n α r 1 1 α r 1 2 α r 1 m C n projection region r 1 Π S ξ(pn) α r 1 pn Figure 6: Proof of Proposition 3 -Case 2
than 2π (see Section 2.4). To be more precise, the curve C n separates the set of the triangles of T n containing p n into two connected regions r 1 and r 2 (see Figure 6). If one denotes by α r1 p the sum of the angles α r1 i of the triangles of region r 1 at p n (resp. by α r2 p the angles α r2 i of the triangles of region r 2 at p n ), since C n is a geodesic, one has:

α r1
pn ≥ π and α r2 pn ≥ π.

By Lemma 1, the angular defect |2π -(α r1 p + α r2 p )| is less than 2Kπα 2 n . We then have π ≤ α r1 p ≤ π + 2Kπα 2 n . Let denote by α r1 i the angle of the projection onto Π S ξ(pn) of α r1 i . We denote by α r1 p the sum of the α r1 i . By using again Lemma 1, we have α r1 pα r1 p ≤ 2Kπα 2 n and then

β T S dev (p n ) = α r1 p -π ≤ α r1 p -α r1 p + α r1 p -π ≤ 4Kπα 2 n .
• The curve C n is included in the neighborhood of p n in two triangles ∆ 1 and ∆ 2 . We denote by Π 1 the plane containing ∆ 1 and by P 1 the orthogonal projection onto Π 1 . Let -→ U 1 and -→ U 2 denote respectively the two unitary vectors colinear to C n in ∆ 1 and in ∆ 2 , oriented with the orientation of curve C n . We put

-→ U 2 ⊥ = P 1 ( -→ U 2 ) - -→ U 2 .
Similarly as before, by using Lemma 1 and the projection P 1 onto Π 1 , we have

∠(P 1 ( -→ U 1 ), P 1 ( -→ U 2 )) ≤ 4Kπα 2 n .
We also have:

-→ U 1 - -→ U 2 = -→ U 1 -P 1 ( -→ U 2 ) + -→ U 2 ⊥ = P 1 ( -→ U 1 ) -P 1 ( -→ U 2 ) + -→ U 2 ⊥ .
Since P 1 ( -→ U 1 ) = 1 and P 1 ( -→ U 2 ) ≥ cos α n , by Lemma 2 we have:

2 sin ∠( -→ U1, -→ U2) 2 = -→ U 1 - -→ U 2 ≤ P 1 ( -→ U 1 ) -P 1 ( -→ U 2 ) + -→ U 2 ⊥ ≤ 2 sin ∠(P1( -→ U1),P1( -→ U2)) 2 + O(α n 2 ) + sin(α n ) ≤ ∠(P 1 ( -→ U 1 ), P 1 ( -→ U 2 )) + O(α n ) = O(α n ).
This implies that there exists K 1 such that Proposition 4 There exists a constant K 2 , such that for any curve C n ∈ R(T n ), one has:

β 3D dev (p n ) = ∠( -→ U 1 , -→ U 2 ) ≤ K 1 α n . 6.4 Majoration of ♯C n a n 1 a n mp n ξ(p n ) γ 1 γ 2
♯(C n ) ≤ K 2 1 + T C T S T an (C n ) + 2 n l(C n ) .

Proof of Proposition 4

If p n is a vertex of T n , let us denote by Cell(p n ) the set of points of the triangulation T n that are closer to p n than to the other vertices of T n . Let η min denote the length of the smallest edge of T n , and by θ min the smallest angle of the triangulation T n . First remark that the smallest altitude is larger than η min sin(θ min ). We put l min = ηmin sin(θmin) 4 . • Let us first consider a curve C n ∈ R(T n ) that satisfies:

l(C n ) ≤ l min and T C T S T an (C n ) ≤ π.
We are going to show that the number ♯C n of intersections between C n and the edges of T n is bounded by a constant independant of n. Case 1: There exists a vertex p n of T n such that the distance from p n to C n is less than l min . Since the length of the curve C n is less than l min , every point m of C n is a distance less than 2l min from the vertex p n . This implies that

C n ⊂ Cell(p n ).
By definition, the curve C n of R(T n ) "crosses" every intersected edge. This implies that C n either contains the vertex p n or is turning around the vertex p n without changing the sense (in the clockwise sense or in the counter-clockwise sense).

-If p n ∈ C n , then C n follows 0, 1 or 2 edges and only contains the vertex p n . This implies that ♯C n is less than the valence of p n which is uniformly bounded from above by a constant V (see Remark 1).

-If C n is turning around p n , we are going to show that C n cannot intersect three times the same edge e. If C n is intersecting twice the same edge e, then the discrete Gauss-Bonnet formulae of the projecton P S ξ(pn) (C n ) of the curve C n onto the plane Π S ξ(p) implies that (see Figure 7)

T C 3D (P S ξ(pn) (C n )) + γ 1 + γ 2 ≥ 2π, where γ 1 ∈ [0, π] and γ 2 ∈ [-π, 0]. We the have |γ 1 + γ 2 | ≤ π and: T C 3D (P S ξ(pn) (C n )) ≥ 2π -(γ 1 + γ 2 ) ≥ π.
Let us denote by a n 1 ,...a n mp n consecutive vertices of P S ξ(pn) (C n ) such that a n 1 and a n mp n belong to the same edge e. Remark that m pn -1 is equal to the valency of p n and thus is less than V . Since the angle between Π S ξ(a n i ) and Π S ξ(pn) is less than 1 2 n , Lemma 4 implies that:

T C T S T an (C n ) ≥ T C 3D (P S ξ(pn) (C n )) - K m pn 2 n ≥ π - K V 2 n .
Let us now suppose that C n is intersecting three times an edge e. Then we have that:

T C T S T an (C n ) ≥ 2π - 2K V 2 n > π,
for n large enough, which contradicts the assumption made on the curve C n . This implies that the number ♯C n of intersections between C n and T n is bounded from above by twice the valency of p n , which is less than 2V . Then ♯C n is uniformly bounded from above. Case 2: The distance between C n and all the vertices of T n is larger than l min . Let ∆ n = p n q n r n denote a triangle of T n that is intersected by C n . The intersection is a segment [a n , b n ] and we denote by θ pn the angle at p n (see Figure 8). We can suppose that p n a n ≤ p n b n and we have: Then the number ♯C n of intersection between C n and the edges of T n is less than:

a n b n ≥ 2 sin θ pn 2 p n a n ≥ 2 sin θ min 2 l min . p n q n r n a n b n θ pn
1 + l(C n ) 2 sin θmin 2 l min ≤ 1 + 1 2 sin θmin 2 .
In the two cases, ♯C n is bounded by a constant independant on n. That implies that there exists a constant K such that:

♯(C n ) ≤ K.
• Let us now consider any curve C n ∈ R(T n ). The curve C n can be subdivided in N curves C 1 n ,...,C N n that satisfy:

l(C i n ) ≤ l min and T C T S T an (C i n ) ≤ π, where N ≤ l(C n ) l min + T C T S T an (C n ) π + 1.
Then there exists a constant K 2 such that:

♯(C n ) ≤ N K ≤ K 2 1 + T C T S T an (C n ) + 2 n l(C n ) .

A sufficient condition for the regularity of the limit curve

Proposition 5 Let (C n ) n∈N be a sequence of polygonal curves C n : [0, 1] → R 3 , with uniform parametrization that converges toward a non constant curve C in the sup norm sense. If

∃k 1 , k 2 , ∀t a , t b ∈ [0, 1] T C 3D (C n , t a , t b ) ≤ k 1 2 n + k 2 l(C n , t a , t b ),
then the curve C has uniform parametrization, is of class C 1,1 and has curvature bounded by k 2 . Moreover l(C) = lim n→∞ l(C n ) and for any t 0 ∈ (0, 1):

lim n→∞ dC n dt + (t 0 ) = dC dt (t 0 ).
Proof One first proves a few lemmas.

Lemma 5 In the conditions of Proposition 5, for any θ > 0, there is an integer number N such that for any n, m ≥ N and any t ∈ [0, 1):

∠ dC n (t) dt + , dC m (t) dt + < θ. ( 7
)
Proof Recall that, because the curves C n have uniform parametrization, dCn dt + (t) is constant on [0, 1]. And, as C is non constant and

C n → C in the sup norm sense, let t < t ′ be such that C(t ′ ) = C(t). Then, for some N 0 , n ≥ N 0 ⇒ C -C n ∞ < 1 4 C(t ′ ) -C(t) , which entails: C n (t ′ ) -C n (t) > 1 2 C(t ′ ) -C(t) .
Therefore, for l = 1 2 C(t ′ ) -C(t) , for any n ≥ N 0 , dCn dt + (t) > l > 0. Let N ≥ N 0 be such that:

k 1 2 N < θ 16 , (8)
and:

∀n ≥ N, C n -C ∞ < θ 2 128k 2 . (9)
Notice that if Lemma 5 holds for θ < min π 2 , k 2 l , it holds in general. We proceed by contradiction. Let us assume that the assertion of the lemma does not hold for some t ∈ [0, 1 2 ], and θ < min π 2 , k 2 l . More precisely, let us assume that, for some t ∈ [0, 1 2 ] and n, m ≥ N , one has:

∠ dC n (t) dt + , dC m (t) dt + ≥ θ. (10)
Without loss of generality, one can assume that:

dC n dt + (t) ≥ dC m dt + (t) . (11)
We consider:

t ′ = t + θ 16k 2 dCn dt + (t) . Notice that θ < k 2 l entails t ′ < 1. One has l(C n , t, t ′ ) ≤ θ
16k2 which gives, with the assumption in Proposition 5 together with inequation (8):

T C 3D (C n , t, t ′ ) ≤ k 1 2 n + k 2 l(C n , t, t ′ ) ≤ θ 8 .
Therefore, one has, for any τ ∈ [t, t ′ ]:

dC n dt + (τ ) - dC n dt + (t) ≤ 2 sin( θ 16
) dC n dt + , and:

t ′ t dC n dt + (τ )dτ -(t ′ -t) dC n dt + (t) = t ′ t dC n dt + (τ ) - dC n dt + (t) ≤ (t ′ -t) dC n dt + (τ ) - dC n dt + (t) ≤ (t ′ -t)2 sin( θ 16 ) dC n dt + = 2 sin( θ 16 ) θ 8k 2 .
Which gives:

t ′ t dC n dt + (τ )dτ -(t ′ -t) dC n dt + (t) ≤ θ 2 64k 2 . ( 12
)
And similarly, using inequation [START_REF] Federer | Curvature measures[END_REF], one has:

t ′ t dC m dt + (τ )dτ -(t ′ -t) dC m dt + (t) ≤ θ 2 64k 2 . (13)
On another hand, inequations [START_REF] Federer | Curvature measures[END_REF] and [START_REF] Farin | Curves and surfaces for computer aided geometric design, a practical guide[END_REF] entail:

(t ′ -t) dC n dt + (t) -(t ′ -t) dC m dt + (t) ≥ sin θ (t ′ -t) dC n dt + (t) ,
which gives:

(t ′ -t) dC n dt + (t) -(t ′ -t) dC m dt + (t) ≥ sin θ θ 8k 2 > θ 2 16k 2 . (14)
But:

t ′ t dCn dt + (τ )dτ - t ′ t dCm dt + (τ )dτ = C n (t ′ ) -C n (t) -C m (t ′ ) + C m (t) ≤ 2 C n -C ∞ ,
which gives, using (9):

t ′ t dC n dt + (τ )dτ - t ′ t dC m dt + (τ )dτ ≤ θ 2 32k 2 .
This last equation can not hold together with inequations ( 14), ( 12) and [START_REF] Kimmel | Finding shortest paths on surfaces[END_REF]. Therefore, inequation [START_REF] Farin | Curves and surfaces for computer aided geometric design, a practical guide[END_REF] does not hold and the lemma is proved for t ∈ [0, 1 2 ]. By reparametrizattion of the curves by t → 1t, one gets the same property for t ∈ [ 1 2 , 1] but expressed with the left derivatives. However, in the condition of the lemma, the left and right derivatives satisfy:

∀n, ∀t, dC n dt -(t) - dC n dt + (t) ≤ k 1 2 n .
This entails, using that ∀n ≥ N 0 , dCn dt + (t) > l > 0 that both have same uniform limit angle.

Lemma 5 gives a Cauchy sequence condition on the angle allows to derive the following lemma on Cauchy sequence conditions on right derivatives.

Lemma 6

In the conditions of Proposition 5, the sequence of dCn dt + (t) is bounded by some number l. Moreover, it is a Cauchy sequence, that is, for any ǫ > 0, there is an integer number N such that for any n, m ≥ N and any t ∈ [0, 1):

dC n (t) dt + - dC m (t) dt + < ǫ. (15)
Proof We first claim that, for any β > 0, there is N such that, for any n, m ≥ N and t ∈ [0, 1 2 ], one has:

dC n (t) dt + - dC m (t) dt + ≤ β max dC n (t) dt + , dC m (t) dt + . (16)
As in the proof of Lemma 5, we consider N 0 and l > 0 such that ∀n ≥ N 0 , dCn dt

+ (t) > l > 0. We consider θ < min π 2 , k 2 l such that (1 -cos θ) < β 2 and N 1 ≥ N 0 such that k 1 2 N1 < θ 2 . ( 17
)
We consider again N ≥ N 1 , using Lemma 5, such that, ∀n, m ≥ N :

∠ dC n (t) dt + , dC m (t) dt + < θ (18)
and

C n -C m ∞ ≤ θβ 8k 2 . (19)
We consider some t ∈ [0, 1 2 ] and n, m ≥ N , and we will prove that (16) holds. Without loss of generality, one can assume that dCn(t)

dt + ≥ dCm(t) dt +
. As in the proof of Lemma 5, we consider the interval [t, t ′ ]:

t ′ = t + θ 2k 2 dCn dt + (t) , (20)
which gives (t ′t) dCn dt + (t) = θ 2k2 and we get from (17) that, ∀τ ∈ [t, t ′ ]:

∠ dC n dt + (τ ), dC n dt + (t) ≤ θ. (21)
We consider the unitary vector e = dCn dt + (t) dCn dt + (t) . One has, using Inequality [START_REF] Streeter | Jr Gross morphology and fiber geometry of the heart[END_REF]:

< e, C n (t ′ ) -C n (t) > = < e, t ′ t dC n dt + (τ )dτ > = t ′ t < e, dC n dt + (τ ) > dτ ≥ cos θ (t ′ -t) dC n dt + .
Similarly, using ( 18) and ( 21), one gets:

< e, C m (t ′ ) -C m (t) > = < e, t ′ t dC m dt + (τ )dτ > = t ′ t < e, dC m dt + (τ ) > dτ ≤ (t ′ -t) dC m dt + .
Using [START_REF] Peyré | Geodesic Computations for Fast and Accurate Suface Remeshing and Parametrization[END_REF] we get:

|< e, C n (t ′ ) -C n (t) > -< e, C m (t ′ ) -C m (t) >| ≤ θβ 4k 2 ,
and the three last inequalities sum up in:

cos θ (t ′ -t) dC n dt + ≤ (t ′ -t) dC m dt + + θβ 4k 2 .
Dividing both terms by (t

′ -t) = θ 2k2 dCn dt + (t) gives: cos θ dC n dt + (t) ≤ dC m dt + + β 2 dC n dt + (t) ,
which gives:

dC n dt + (t) - dC m dt + ≤ ( β 2 + 1 -cos θ) dC n dt + (t) .
The fact that (1 -cos θ) < β 2 proves (16). Property [START_REF] Memoli | Distance functions and geodesics on submanifolds of R d and point clouds[END_REF] easily proves the lemma. Indeed, taking β = 1 2 gives that, for some N :

n ≥ N ⇒ dC n dt + (t) - dC N dt + (t) ≤ 1 2 max dC n dt + (t) , dC N dt + (t) ,
which entails:

n ≥ N ⇒ dC n dt + (t) ≤ 2 dC N dt + (t) .
Therefore, the sequence of dCn dt + (t) is bounded by some number l and ( 16) entails that, for any β > 0, one has:

n ≥ N ⇒ dC n dt + (t) - dC N dt + (t) ≤ β l.
This implies that dCn dt + (t) is a Cauchy sequence. This fact, combined with Lemma 5, implies that dCn dt + (t) is a Cauchy sequence. The fact that for every t ∈ (0, 1), one has l(C n ) = dCn dt + (t) , implies by Lemma 6, that the sequence of lengths l(C n ) is converging to L. For two points a and b, d(a, b) denotes the euclidean distance between a and b. First one proves the following: Lemma 7 In the conditions of Proposition 5, for any ǫ with 0 < ǫ < 1 10 , there is η and an integer number N such that if 0 < t bt a ≤ η, then for any t, t ′ such that t a ≤ t < t ′ ≤ t b and for any n ≥ N :

C n (t ′ ) -C n (t) t ′ -t - C n (t b ) -C n (t a ) t b -t a ≤ ǫ C n (t b ) -C n (t a ) t b -t a (22)
and:

l(C n )(t b -t a )(1 -ǫ 2 ) ≤ d (C n (t a ), C n (t b )) ≤ l(C n )(t b -t a ). ( 23
)
Proof of Lemma 7

We put η = ǫ 4k2 l and N such that k1 2 N ≤ ǫ 4 . If 0 < t bt a ≤ η, we then have:

l(C n , t a , t b ) = l(C n ) (t b -t a ) ≤ l ǫ 4 k 2 l = ǫ 4 k 2 .
We then have

T C 3D (C n , t a , t b ) ≤ k 1 2 n + k 2 l(C n , t a , t b ) ≤ ǫ 4 + ǫ 4 = ǫ 2 .
First remark that the right derivative of C n exists everywhere. In the following, we denote it by dCn dt + . Furthermore, since C n is differentiable almost everywhere, the Lebesgues integral of its derivative is equal to the Lebesgues integral of dCn dt + . For any t 1 , t 2 ∈ [t a , t b ], one has:

∠ dC n dt + (t 1 ), dC n dt + (t 2 ) ≤ ǫ 2 .
Then, since C n has uniform parametrization, one has:

∀t 1 ∈ [0, 1] dC n dt + (t 1 ) = l(C n ),
and therefore, for any t 1 , t 2 ∈ [t a , t b ], one has:

dC n dt + (t 2 ) - dC n dt + (t 1 ) ≤ 2 l(C n ) sin ǫ 4 .
That implies that for any t 1 ∈ [t a , t b ] and t, t ′ with t a ≤ t < t ′ ≤ t b :

C n (t ′ ) -C n (t) t ′ -t - dC n dt + (t 1 ) = 1 t ′ -t t ′ t dC n dt + (τ ) - dC n dt + (t 1 ) dτ ≤ 2l(C n ) sin ǫ 4 .
Again:

C n (t ′ ) -C n (t) t ′ -t - C n (t b ) -C n (t a ) t b -t a = 1 t b -t a t b ta C n (t ′ ) -C n (t) t ′ -t - dC n dt + (τ ) dτ ≤ 1 t b -t a t b ta C n (t ′ ) -C n (t) t ′ -t - dC n dt + (t 1 ) dτ ≤ 2l(C n ) sin ǫ 4 . (24)
On the other hand:

C n (t b ) -C n (t a ) t b -t a 2 = 1 t b -t a t b ta dC n dt + (τ )dτ . 1 t b -t a t b ta dC n dt + (τ )dτ = 1 (t b -t a ) 2 t b ta t b ta dC n dt + (τ 1 ) • dC n dt + (τ 2 )dτ 2 dτ 1 .

This gives

l(C n ) 2 cos ǫ 2 ≤ C n (t b ) -C n (t a ) t b -t a 2 ≤ l(C n ) 2 and l(C n ) cos ǫ 2 ≤ C n (t b ) -C n (t a ) t b -t a ≤ l(C n ). (25)
Using ǫ < 1 10 , Equations ( 24) and (25) imply Equation [START_REF] Warren | Subdivision Methods for Geometric Design: A constructive Approach[END_REF]; Equation (25) proves Equation (23).

We recall that L = lim sup n→∞ l(C n ). Now, by using the fact that C n -C ∞ → 0, Lemma 7 gives that, for any ǫ with 0 < ǫ < 1 10 , there is η such that, if 0 < t bt a ≤ η:

C(t ′ ) -C(t) t ′ -t - C(t b ) -C(t a ) t b -t a ≤ ǫ C(t b ) -C(t a ) t b -t a .
Remark that:

C n (t b ) -C n (t a ) t b -t a ≤ l(C n ).
Which entails, by taking the limit of both sides as n → ∞:

C(t b ) -C(t a ) t b -t a ≤ L.
We have then:

C(t ′ ) -C(t) t ′ -t - C(t b ) -C(t a ) t b -t a ≤ ǫL. ( 26
)
For a given real number t 0 ∈ [0, 1], and an integer number j ≥ 1 we introduce the closed subset K j (t 0 ) as:

K j (t 0 ) = Closure C(t ′ ) -C(t) t ′ -t | max 0, t 0 - 1 j ≤ t < t ′ ≤ min 1, t 0 + 1 j .
From Equation (26) K j (t 0 ) is bounded for j large enough, and therefore compact. It is obviously not empty and decreasing for the inclusion: K j (t 0 ) ⊃ K j+1 (t 0 ). From compactness, the set K(t 0 ) = ∩ j≥1 K j (t 0 ) is not empty and Equation ( 26) entails that it must be contained in a ball of radius ǫL for arbitrarily small ǫ which entail that K(t 0 ) is a single point:

K(t 0 ) = {ψ(t 0 )}.
Again, Equation (26) entails:

∀ǫ > 0, ∃h > 0, max (0, t 0 -h) ≤ t < t ′ ≤ min (1, t 0 + h) ⇒ C(t ′ ) -C(t) t ′ -t -ψ(t 0 ) < ǫ.
Therefore, for 0 < t 0 < 1, ψ(t 0 ) is the derivative of C at t 0 . In fact, the expression above is stronger: it states that C is strictly differentiable (see [START_REF] Clarke | Optimization and NonSmooth Analysis[END_REF] page 30), which entails that the derivative function is continuous. For t 0 = 0 (resp. t 0 = 1 ) this shows that C has a right (resp. left) derivative at 0 (resp. 1). We have proven so far that C is of class C 1 .

Lemma 8

The sequence of right derivatives dCn dt + uniformly converges to dC dt . In other words, for every ǫ >, there exists N ∈ N such that:

n > N ⇒ ∀t ∈ (0, 1) dC n dt + (t) - dC dt (t) ≤ ǫ.
Proof Proof of Lemma 8

We know by Lemma 6 that dCn dt + uniformly converges. We only have to show that dCn dt + (t 0 ) converges to dC dt (t 0 ) for any t 0 ∈ (0, 1). From Lemma 7, for any ǫ such that 0 < ǫ < 1 10 , there is h 0 > 0 and an integer number N such that, for n ≥ N and h ≤ h 0 , one has:

C n (t 0 + h) -C n (t 0 -h) 2h - dC n dt + (t 0 ) ≤ ǫ.
Since C is of class C 1 , we can take h small enough to have:

C(t 0 + h) -C(t 0 -h) 2h - dC dt (t 0 ) ≤ ǫ and let N ≥ N be such that n ≥ N ⇒ C n -C ∞ ≤ ǫh. One gets: C n (t 0 + h) -C n (t 0 -h) 2h - C(t 0 + h) -C(t 0 -h) 2h ≤ 2ǫh 2h = ǫ,
which gives, for any n ≥ N :

dC n dt + (t 0 ) - dC dt (t 0 ) ≤ 3ǫ.
This is true for arbitrary small ǫ which proves the claim.

Lemma 8 allows to apply Lebesgues' dominated convergence theorem which entails that, for 0 ≤ t a < t b ≤ 1:

lim n→∞ l(C n , t a , t b ) = l(C, t a , t b ), (27)
and in particular:

lim n→∞ l(C n ) = L = l(C).
The assumptions of Proposition 5 entail that:

∀t a , t b ∈ [0, 1] ∠ dC n dt + (t a ), dC n dt + (t b ) ≤ k 1 2 n + k 2 l(C n , t a , t b ).
This, together with Lemma 8 and Equation ( 27) entail:

∀t a , t b ∈ [0, 1] ∠ dC dt (t a ), dC dt (t b ) ≤ k 2 l(C, t a , t b ),
which proves that k 2 bounds the curvature of C. Moreover, Lemma 8 entails that, for any t 0 with 0 < t 0 < 1:

dC dt (t 0 ) = L = l(C).
It follows that C has uniform parametrization and dC dt is K-Lipschitz with K = L 2 k 2 .

Proof of Theorem 1

We take the notations of Theorem 1. Let α n be the smallest real number such that for every m ∈ T n , the angle between any triangle ∆ containing m and the tangent plane Π S ξ(m) of S at ξ(m) is smaller than α n . By assumption, we have

α n ≤ K 2 n .
Step 1: Let us consider a given interval [t a , t b ] ⊂ [0, 1]. By definition, we have:

T C T S T an (C n , t a , t b ) = p vertex of Cn([ta,t b ]) β T S dev (p).
Propositions 3 and 4 imply that:

T C T S T an (C n , t a , t b ) ≤ K 1 (2α n ) 2 ♯ (C n , t a , t b ) (by Prop. 3) ≤ K 1 K 2 (2α n ) 2 1 + T C T S T an (C n , t a , t b ) + 2 n l(C n , t a , t b ) (by Prop. 4) ≤ 4K 2 4 n K 1 K 2 1 + T C T S T an (C n , t a , t b ) + 2 n l(C n , t a , t b ) .
Therefore, for some constants K 3 and K 4 independant of t a and t b , one has:

T C T S T an (C n , t a , t b ) ≤ K 3 2 n l(C n , t a , t b ) + K 4 4 n . (28)
Step 2: We have:

T C 3D (C n , t a , t b ) = p vertex of Cn(ta,t b ) β 3D dev (p) ≤ K 1 (2α n ) ♯ (C n , t a , t b ) (By Prop. 3) ≤ 2K 1 K 2 α n 1 + T C T S T an (C n , t a , t b ) + 2 n l(C n , t a , t b ) (By Prop. 4) ≤ 2KK1K2 2 n 1 + T C T S T an (C n , t a , t b ) + 2 n l(C n , t a , t b )
. Equation (28) then implies that there exist constants K and K ′ :

T C 3D (C n , t a , t b ) ≤ K l(C n , t a , t b ) + K ′ 1 2 n , (29)
and Proposition 5 then implies that the curve C is of class C 1,1 .

Step 3

We consider a point p 0 = C(t 0 ) and we recall that P S p0 denotes the projection on the plane Π S p0 tangent to S at p 0 . In this step, we are are going to prove the following lemma: Lemma 9 For any t 0 ∈ (0, 1), the projection P S C(t0) • C of C on the plane tangent to S at C(t 0 ) is twice derivable at t 0 and:

d 2 P S C(t0) • C dt 2 ˛t = t 0 = 0

Proof of Lemma 9

Let α be such that 0 < α < 1 and let r > 0 be the reach of the surface S. We recall the following proposition (see [START_REF] Federer | Curvature measures[END_REF] page 435):

Proposition 6 In the ball B(p 0 , αr), the map ξ is 1 1-α -Lipschitz. Let s be such that s l = αr 2 . Using the fact that C n has uniform parametrization and length upper bounded by l, the respective lengths of the arcs C n ([t 0s, t 0 ]) and C n ([t 0 , t 0 + s]) are smaller than αr 2 and therefore one has:

C n ([t 0 -s, t 0 + s]) ⊂ B C n (t 0 ), αr 2 .
Let us now take n large enough such that C n -C ∞ ≤ αr 2 . We then have:

C n ([t 0s, t 0 + s]) ⊂ B(p 0 , αr).

Using Proposition 6, it follows that the curves ξ • C n ([t 0s, t 0 ]) and ξ • C n ([t 0 , t 0 + s]) have length bounded by αr 1-α . Therefore, since the curvature of S is bounded by 1 r , if one denotes by Π S ξ(Cn(t)) the tangent planes to S at the point ξ (C n (t)), one has, for any t ∈ [t 0s, t 0 + s]:

2 sin ∠ Π S ξ(Cn(t)) , Π S ξ(Cn(t0)) 2 ≤ α 1 -α . ( 30
)
We consider now the sequence of curves C n which are the projections P S p0 (C n ) and which converge toward the projection P S p0 •C of the curve C on the plane Π S p0 . We consider the arc curve P S p0 (C ([t 0s, t 0 + s])) in the plane Π S p0 , for s = αr e l . Let t n i be the parameter of the i th vertex p n i = C n (t n i ) of C n . We apply Lemma 3, taking the unit vectors along dCn dt just before and just after the vertex p n i for the vector u and v of the proposition, and the projections P S ξ(p n i ) and P S p0 respectively for the projections P

β Π0 i ≤ K β T S dev (p n i ) + sin θ β 3D dev (p n i ) + 1 4 n ,
where θ is the angle between the planes Π S ξ(p n i ) and Π S p0 and satisfies sin θ ≤ 2 sin θ 2 ≤ α 1-α from Equation 30. One has then:

β Π0 i ≤ K β T S dev (p n i ) + α 1 -α β 3D dev (p n i ) + 1 4 n .
Therefore, by summing over all the vertices p n i , one has:

T C 3D P S p0 • C n , t 0s, t 0 + s ≤ K ′ T C T S T an (C n , t 0s, t 0 + s) + ♯ (C n , t 0s, t 0 + s)

4 n + α 1 -α
T C 3D (C n , t 0s, t 0 + s) .

Equation (28) gives:

T C T S T an (C n , t 0s, t 0 + s) ≤

K 3 2s 2 n + K 4 4 n .
Equation (29) implies:

T C 3D (C n , t 0 -s, t 0 + s) ≤ K 2s l + K ′ 1 2 n .
Proposition 4 implies that:

♯ (C n , t 0s, t 0 + s) ≤ K 2 1 + T C T S T an (C n , t 0s, t 0 + s) + 2 n 2s l .

Then, by combining all these results and using the fact that α = e ls r , we have, for some constant K 5 , K 6 and K 7 :

T C 3D P S p0 • C n , t 0 -s, t 0 + s ≤ K 5 4 n + K 6 s 2 n + K 7 s 2 .
The vector C ′ (t 0 ) belongs to the tangent plane of S at C(t 0 ). Furthermore, the fact that Df (0, 0) = 0 implies that D 2 f (γ(t 0 )).(γ ′ (t 0 ), γ ′ (t 0 )) = II C(t0) (C ′ (t 0 )),

where II C(t0) is the second fundamental form of S at the point C(t 0 ) [START_REF] Do | Differential geometry of curves and surfaces[END_REF]. We then have

C ′′ (t 0 ) = II C(t0) (C ′ (t 0 )) N S C(t0) .
In this expression, C ′′ (t 0 ) depends continuously on t 0 . That implies that C is of class C 2 in t 0 . The function C is then of class C 2 . Lemma 9 then implies that C has zero geodesic curvature, and then is a geodesic [START_REF] Do | Differential geometry of curves and surfaces[END_REF].

Proof of Corollary 1

We first need to check that the assumptions a) to e) of Theorem 1 are satisfied. The uniform convergence of P n to f clearly implies assumption a). Now, since the map ξ realises the distance to S, for every m = P n (u, v) ∈ P n , one has ξ(m)m ≤ f (u, v) -P n (u, v) which implies assumption b). By using the regularity of f and by compacity, we have that: Let ∆ n be a triangle of P n . The vertices of ∆ n are for example of the form p n = P n i 2 n , j 2 n q n = P n i + 1 2 n , j 2 n and r n = P n i + 1 2 n , j + 1 2 n . Now the fact that (P n ) n∈N uniformly converges in derivative to f implies that there exists N and K such that for every n > N , the lengths of p n q n , p n r n and q n r n are greater than 1 2 n m -K

2 n
and less than

1 2 n M + K 2 n .
If N is large enough we then have:

1 2 n m 2 ≤ p n q n ≤ 1 2 n 2 M, 1 2 n m 2 ≤ p n r n ≤ 1 2 n 2 M, 1 2 n m 2 ≤ r n q n ≤ 1 2 n 2 M, (32)
In particular, assumption d) is proved. The fact that P n uniformly converges in derivative to f also implies that there exists N and K such that for every n > N :

4 n--→ q n r n ∧ --→ q n p n - ∂f ∂u ∧ ∂f ∂v i 2 n , j 2 n ≤ K 2 n . ( 33
)
8

Conclusion and future works

The main result of this work gives sufficient conditions for a sequence of geodesics on PL-surfaces to converge toward a geodesic on a smooth limit surface. We believe this is a significant step toward an effective notion of geodesic: indeed, the usual definition of geodesic is not effective because it relies on the notions of smooth curves and surfaces and on the pointwise curvature which can not be exactly represented on computers. Our main theorem states that the usual notion of geodesic coincides with the limit of a sequence of PL-curves that can be represented (at least if one restricts ourselves to PL-surfaces with rational vertices coordinates). Therefore, by using our result, a realistic algorithm can output a sequence of curves whose limit is a geodesic of a smooth surface. Notice that, given a smooth surface with bounded curvature, there exists a sequence of PL-surfaces converging to it (and that matches the conditions of our theorem). However, in order to completely get the effective notion of geodesic, one still has to quantify the rate of convergence of this sequence of curves. We also believe that our result could be improved by relaxing the condition on the edge lengths: indeed, in the counter-example the lengths decrease with the order 1 4 n with respect to a decrease rate of 1 2 n of the angular convergence. We believe that it is possible to improve the theorem between the K 2 n condition of the theorem and the 1 4 n of the counter-example.

Another possible improvement of the result is to suppose that the limit surface is of class C 1,1 (instead of C 2 ). Notice that such a generalisation would be very usefull for some subdivision surfaces with extraordinary points. Indeed, at extraordinary points, the limit surface of some subdivision surfaces is only of class C 1,1 . We proved (for example for the Catmull-Clark scheme) that if the limit curve of a sequence of geodesics does not contain extraordinary points, then it is a geodesic. We believe that the result still holds if the limit curve contains extraordinary points.

  a) T n belongs to the tubular neighborhood U r (S) of radius r of S; b) for every m ∈ T n , mξ(m) ≤ d n ; c) for every m ∈ T n , the angle between any triangle ∆ containing m and the tangent plane Π S ξ(m) of S at ξ(m) is smaller than K 2 n ; d) the lengths of the edges of T n are greater than e K 2 n ; e) all the angles of T n are greater than θ min ; Let (C n ) n∈N be a sequence of polygonal curves C n : [0, 1] → R 3 with uniform parametrization such that C n is an interior geodesic of T n and ξ(C n ) does not intersect the boundary of S. If (C n ) n∈N converges toward a curve C in the sup norm sense, then C is of class C 2 and is a geodesic of S.

Figure 3 :

 3 Figure 3: Deviation angle of the curve and of its projection

Figure 4 :Lemma 2

 42 Figure 4: Proof of Lemma 2

Figure 5 :

 5 Figure 5: Proof of Proposition 3 -Case 1

Figure 7 :

 7 Figure 7: Projection of C n onto the plane Π S ξ(pn) (Proposition 4 -case 1)

Figure 8 :

 8 Figure 8: Proof of Proposition 4 -case 2

  uniformly converging to a curve C. Remark now that for every t ∈ (0, 1), the pointξ n (C n (t)) is the closest point of S to C n (t). That implies that C n (t)ξ n • C n (t)is smaller than the Hausdorff distance between T n and S, which tends to 0. Therefore, the subsequence (C n k ) k∈N is also uniformly converging to C.The fact that C n uniformly converges to C implies that l

  . The points d 1 n , d 2 n , p n and m n are on the plane Π H while the points t 1 n and t 2 n stand at some height above the

	an
	geodesic Cn
	mn	d 1 n	mn
		t 1 n	t 2 n
			β1
		d 2 n	β2
			pn
			β3
			β4
	nn
	bn		α = 2π 4 n
	Figure 2: Triangulation T n and geodesic C n seen from above: on the left we
	see the whole surface; the region in the dashed rectangle is depicted on the
	right
	plane. The faces d 1 n t 1 n d 2 n , p n t 2 n m n , d 1 n t 1 n t 2 n m n and d 2 n t 1 n t 2 n p n are planar and all make a slope s n with Π H . If we take s n = 3 2 n , one has, for each n

Let t ∈ [t 0s, t 0 + s]. Lemma 2 implies that there exists K 8 such that:

We finally have:

(31) Let ǫ > 0. By Lemma 8, there exists N such that for every n > N and for every u ∈ (0, 1):

We then have for every u ∈ (0, 1):

By using Equation (31), we have:

We finally get:

This allows to conclude the proof.

Step 4 Let t 0 ∈ (0, 1). In a neighborhood of C(t 0 ), the surface S can be parametrized by its tangent plane at C(t 0 ): in an appropriated frame with origin C(t 0 ), the surface is parametrized by (x, y)

where U 0 is a neighborhood of (0, 0), f is a function of class C 2 that satisfies f (0, 0) = 0 and Df (0, 0) = 0. For every t close enough to t 0 , we put γ(t) = P S C(t0) (C(t)). The function γ is of class C 1,1 , and by Lemma 9, we know that it is twice differentiable in t 0 and that γ ′′ (t 0 ) = 0. We have in a neighborhood of t 0 ,

.

The function C is then twice differentiable in t 0 and we have

.

If we take N such that K 2 n ≤ m 2 , we then have

Together with Equation (32), that implies that

The angle ∠( --→ q n p n , --→ q n r n ) is then lower bounded by a constant independant on n. The same result holds with the other angles of T n , which proves assumption e).

We now denote by N Tn qn a unitary normal of the triangle p n q n r n and N S i,j the normal of S at f ( i 2 n , j 2 n ). By using Equation (33) and Lemma 2, we have:

2 sin θmin,m " .

Let now m n ∈ ∆ n . Then we need to bound the angle ∠(N Tn qn , N S ξ(mn) ). There exists K such that

Then, by using Proposition 6, we have:

which implies that there exists k such that:

This result holds for all the triangles of T n . Assumption c) is then proved. We now only need to prove that ξ(C n ) is interior to S for n large enough: the curve C is interior to S. The compacity of S and C implies that the distance from C to the boundary of S is more than η > 0. The curve ξ(C n ) clearly tends to C. That implies that for n large enough ξ(C n ) is an interior curve.

Acknowledgements

We acknowledge Cédric Gérot for his advises concerning subdivision surfaces.


Huseyin Ugur Yildiz

Vehbi Cagri Gungor

Bulent Tavli

A Hybrid Energy Harvesting Framework for Energy Efficiency in Wireless Sensor Networks Based Smart Grid Applications

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A Hybrid Energy Harvesting Framework for Energy

Efficiency in Wireless Sensor Networks Based Smart Grid Applications Huseyin Ugur Yildiz * , Vehbi Cagri Gungor † , and Bulent Tavli ‡ * TED University, 06420, Ankara, Turkey, e-mail: hugur.yildiz@tedu.edu.tr † Abdullah Gul University, 38080, Kayseri, Turkey, e-mail: cagri.gungor@agu.edu.tr ‡ TOBB University of Economics and Technology, 06560, Ankara, Turkey, e-mail: btavli@etu.edu.tr Abstract-In smart grid applications, Wireless Sensor Networks (WSNs) which consist of battery limited sensor nodes are used on critical equipments of power distribution grids for monitoring purposes. WSN nodes have tight energy constraints hence it is important to reduce the energy consumption of sensor nodes due to harsh propagation characteristics of smart grid environment. One possible way to reduce the energy consumption is to utilize transmission power control where transmission powers are adjusted according to channel conditions. Another technique is to employ energy harvesting schemes to provide additional power for nodes by using environmental energy sources. Solar and electromagnetic energies are two possible environmental energy sources in outdoor substation environments. Solar energy can be efficiently exploited in a sunny day. On the other hand, electromagnetic energy can be used at any time. In this work, we propose a hybrid energy harvesting model that exploits both solar and electromagnetic energies and develop a Mixed Integer Programming (MIP) method to minimize the energy dissipation of sensor nodes. By using the MIP framework, we quantify the impact of the proposed hybrid energy harvesting model as well as transmission power control on the energy saving of nodes.

Index Terms-smart grids, wireless sensor networks, energy harvesting, energy efficiency, mixed integer programming.

I. INTRODUCTION

Traditional power distribution grids spread over a large area that are monitored and diagnosed through wired communications which have high capital and operational expenditures [START_REF] Gungor | Opportunities and challenges of wireless sensor networks in smart grid[END_REF]. In order to overcome such high expenditures, grid monitoring can be performed by using wireless communications techniques. Wireless Sensor Networks (WSNs) are proven to be a viable solution for grid monitoring applications such as power fraud, fault detection, outage detection, power line automation, etc. The capital and operational expenditures of WSNs are cheaper than conventional wired communication techniques [START_REF] Sahin | Quality-of-service differentiation in single-path and multi-path routing for wireless sensor network-based smart grid applications[END_REF]. The integration of modern communication methods (such as WSNs) as well as monitoring, automation, and control abilities into power grid infrastructure to improve the reliability, productivity, efficiency, etc. creates the paradigm of smart grids (SGs) [START_REF] Bilgin | On the performance of multi-channel wireless sensor networks in smart grid environments[END_REF]. Although WSNs are proven to be a promising solution for SG systems, tight battery limitations of WSN nodes create a bottleneck for long term monitoring.

WSN based SG systems experience extremely harsh environmental challenges such as fading and interference caused by the generated electromagnetic fields [START_REF] Fadel | A survey on wireless sensor networks for smart grid[END_REF]. There are several ways to combat the extreme propagation characteristics of SGs. One possible way is to employ transmission control approaches. Increasing the transmission power levels globally in the network helps to reduce the packet errors. However, the price paid for this case is the increment in total energy dissipation of the network. Another way is to adjust the transmission power level on each link according to the channel conditions (i.e., link-level) which can, potentially, prevent the excessive energy dissipation in the network.

Although transmission power control approaches may reduce the overall energy consumption in the network, additional energy can be provided to nodes by using energy harvesting approaches. In the literature, some power sources available in the SG environment are solar and electromagnetic energies [START_REF] Han | A nonintrusive power supply design for self-powered sensor networks in the smart grid by scavenging energy from AC power line[END_REF]. Considering an outdoor SG application, solar energy can be harvested during the day time by using solar panels (or photovoltaic cells) as long as the weather is not cloudy [START_REF] Erol-Kantarci | Suresense: sustainable wireless rechargeable sensor networks for the smart grid[END_REF]. However, if the weather is cloudy or during the night time, then solar energy cannot be harvested efficiently. At this point electromagnetic field generated by the alternating current (AC) flowing through distribution lines can be used as a backup energy source to be exploited by using harvester devices that are either clamped or not clamped around the conductor [START_REF] Moghe | A scoping study of electric and magnetic field energy harvesting for wireless sensor networks in power system applications[END_REF]. Fluctuation of the solar energy generation rate due to the stochastic nature of weather conditions motivates us to develop a hybrid energy harvesting method by using the electromagnetic energy as an alternative energy source.

Our contributions in this work are itemized as follows:

1) We propose a hybrid energy harvesting model that exploits either solar [START_REF] Hsu | Adaptive duty cycling for energy harvesting systems[END_REF], [START_REF] Yao | Quality-driven energy-neutralized power and relay selection for smart grid wireless multimedia sensor based IoTs[END_REF] or electromagnetic energies [START_REF] Bhuiyan | A miniature energy harvesting device for wireless sensors in electric power system[END_REF]- [START_REF] Yuan | Magnetic field energy harvesting under overhead power lines[END_REF] depending on the time and probabilistic weather conditions. 2) We develop a Mixed Integer Programming (MIP) framework which is built on top a handshaking based linklayer model employing transmission power control approaches at link-level by using the physical layer characteristics of a well known SG environment (i.e., outdoor ISBN 978-3-903176-05-8 c 2018 IFIP 500 kV substation [START_REF] Kurt | Packet size optimization in wireless sensor networks for smart grid applications[END_REF]). The proposed MIP framework aims to minimize the energy dissipation of the maximum energy consuming node. This is critical because the lifetime of WSNs is usually limited by the lifetime of the the highest energy consuming node. 3) We quantitatively investigate the performance of the proposed hybrid energy harvesting model in terms of energy dissipation of sensor nodes by using the MIP framework through numerical evaluations. To the best of our knowledge, this is the first study that numerically explores the effects of decrement in energy consumption of sensor nodes when employing a hybrid energy harvesting approach by using a MIP framework that combats the harsh environmental conditions of SGs with the help of transmission power control.

II. SYSTEM MODEL A. Wireless Channel Model

In this work, we choose Tmote Sky motes as our WSN node platform [15]. We use log-normal shadowing channel model to account the propagation loss which has proven to be an accurate channel model [START_REF] Bilgin | On the performance of multi-channel wireless sensor networks in smart grid environments[END_REF]. The signal-to-noise ratio (i.e., γ ij (l) in dB) at the receiving node-j due to the transmission from node-i with power-l can be expressed as

γ ij (l) = P ant tx (l) -P L 0 + 10nlog 10 d ij d 0 -X σ P ant rx,ji (l) -P n . (1)
In this equation, P ant tx (l) is antenna transmit power when using power level-l (in dBm), P ant rx,ji (l) is the receive power (in dBm), P L 0 is the reference path loss value (in dB), n is the path loss exponent, d ij is the link distance (in meters), d 0 is the reference distance (in meters), X σ is a zero-mean Gaussian random variable with variance σ 2 denoting the shadowing random variable (in dB), and P n is the noise floor value which is taken as -94 dBm for Tmote Sky platforms [15].

We adopt two different SG propagation environments namely outdoor 500 kV substation line-of-sight (LOS) and non-line-of-sight (NLOS). Experimentally determined lognormal channel parameters for these environments are based on the field tests conducted at Georgia Power, Atlanta, GA, USA [START_REF] Gungor | Opportunities and challenges of wireless sensor networks in smart grid[END_REF]. Tmote Sky motes use Chipcon CC2420 radios and reported P ant tx (l) values as well as circuitry transmission power consumption values (i.e., P crc tx (l)) for this radio platform is given in Table I [15]. Tmote Sky platforms use offset quadrature phase shift keying as the modulation. After including processing gain costs (i.e., η = 16), the successful reception probability of data packet at node-j due to the transmission of node-i can be calculated as [14]

P s,L ij (l) = 1 -Q η × γ ij (l) L , (2)
where γ ij (l) is the signal-to-noise ratio in ordinary form and L is the data packet size in bits. Similarly, P f,L ij (l) = 1-P s,L ij (l) shows the failure reception probability. B. Link-Layer Model 1) Handshaking Model: We assume a slotted communication scheme at the link-layer where a single slot encompasses durations for M P -byte of data packet transmission (i.e., T tx (M P )), M A -byte of acknowledgement (ACK) packet transmission (i.e., T tx (M A )), propagation delay (T pd ), and guard times (2 × T grd ) which are applied at both the beginning and end of the active slot to prevent synchronization errors [START_REF] Kurt | Packet size optimization in wireless sensor networks for smart grid applications[END_REF]. We choose M P = 128 bytes, M A = 20 bytes, and data rate of Tmote Sky platforms (R) as 250 kbps [15]. Hence, the active slot time is calculated as:

T slot = [2 × T grd + T tx (M P ) + T pd + T tx (M A )] = 4.78 ms.
During an active slot a two-way handshaking should be performed successfully for a reliable communication performance. The successful handshaking is obtained when both M P -byte and M A -byte of data and ACK packets are successfully received at the intended nodes. Thus, the probability of successful handshaking is calculated as

P s,HS ij (l, k) = P s,M P ij (l) × P s,M A ji (k)
where l and k are the transmission power levels for data and ACK packets, respectively. Considering a stop-and-wait automatic repeat request scheme, each packet has to be retransmitted λ lk ij = 1 P s,HS ij (l,k) times on the average.

2) Transmitter Energy Model: Transmitter node spends P crc tx (l) × T tx (M P ) Joules of energy for transmitting M Pbyte of data packet. In the rest of the active slot time (i.e., T slot -T tx (M P ) seconds), the transmitter node consumes P crc rx × (T slot -T tx (M P )) Joules of energy to stay in receive/idle mode for ACK packet reception where P crc rx = 69 mW is the power consumption for reception [START_REF] Kurt | Packet size optimization in wireless sensor networks for smart grid applications[END_REF]. E P P = 12.66 µJ is the packet processing energy which is dissipated only once during a single slot [START_REF] Kurt | Packet size optimization in wireless sensor networks for smart grid applications[END_REF]. Including the effects of retransmissions with the factor of λ lk ij , the total energy dissipation of the transmitter node-i in a single slot is obtained as

E lk tx,ij = E P P + λ lk ij × [P crc tx (l) × T tx (M P ) + P crc rx × (T slot -T tx (M P ))]. (3)
3) Receiver Energy Model: Receiving node-j spends

P crc rx × [T slot -T tx (M A )]
Joules of energy for being in receive/idle mode to receive a data packet from the transmitting node-i. Upon receiving the data packet, node-j consumes P crc tx (k) × T tx (M A ) Joules for ACK packet transmission by using power level-k. Let ε denote the amount of energy dissipated for reception while considering a successful handshake. If the handshaking has failed due to the ACK

(l)×P f,M A ji (k) P s,HS ij (l,k)
is dissipated due to the retransmissions. Finally, if the handshaking error is caused by the failure of the data packet, then receiving node will be in receive/idle mode for the whole slot duration in which

P f,M P ij (l) P s,HS ij (l,k) ×[P crc rx ×T slot ]
Joules of energy is dissipated when considering retransmissions. Including the effects of packet processing, the total energy dissipation of the receiving node-j in a single slot is calculated as

E lk rx,ji = ε P crc rx [T slot -T tx (M A )] + [P crc tx (k) × T tx (M A )] + P s,M P ij (l)P f,M A ji (k) P s,HS ij (l, k) ε + P f,M P ij (l) P s,HS ij (l, k) (P crc rx × T slot ) + E P P . (4)

C. Energy Harvesting Methods

In this part, we present mechanisms for solar and electromagnetic energy harvesting methods utilized throughout this work in Sections II-C1-II-C3. Next, we describe the details of our proposed hybrid energy harvesting model in Section II-C4.

1) Solar Energy Harvester (SEH): For the solar energy harvesting (SEH) method, we assume that photovoltaic solar panels are installed on top of each WSN node. SEH profile given in [START_REF] Yao | Quality-driven energy-neutralized power and relay selection for smart grid wireless multimedia sensor based IoTs[END_REF] is adopted which is measured and averaged between the beginning of June to middle of August for 72 days [START_REF] Hsu | Adaptive duty cycling for energy harvesting systems[END_REF]. In Figure 1, we present the average solar harvested power (i.e., P SEH ) with respect to time. As seen from this figure, solar energy can be exploited between 5:00 a.m. and 7:00 p.m. and the maximum harvested power is measured as 60 mW at 12:00 p.m. 2) Conductor Winding Electromagnetic Harvester (CWEH): A conductor winding electromagnetic harvester (CWEH) is composed of an inductive multi-turn coil where its cross sectional area is surrounding the magnetic core material which is constructed either by using a single or multi layers of high permeability metal. In [START_REF] Bhuiyan | A miniature energy harvesting device for wireless sensors in electric power system[END_REF], multi-layered core structures are proven to have better performance than single layered structures. By using the theoretical analysis given in [START_REF] Bhuiyan | A miniature energy harvesting device for wireless sensors in electric power system[END_REF], the root mean square (RMS) voltage induced at coil terminals for multi-layered core structures can be expressed as

V coil = (N CW EH )µ 0 µ r ωf Iln 1 + h r , (5)
where N CW EH and r denote the number of turns in the coil and the distance between coil and the conductor. µ 0 = 4π10 -7 H/m is the permeability of free space and µ r is the relative permeability. ω and h are width and the thickness of the core material. Finally, I and f represent the RMS magnitude and the frequency of the current. Since µ r ignores the effect of the core gap, we modify µ r to include core gaps. In the following equation, µ rg shows the effective permeability accounting for the core gap which is calculated as

µ rg = lm A 1 A lg + 0.241+ 1 π ln b lg p + lm A µr . (6)
In this equation A = h × ω, p = 2(h + ω), b, l g , and l m are cross-sectional area of the core, perimeter of the cross section of the core, half the inner toroidal core perimeter, core gap length, and average core perimeter without the gap, respectively. l m can be calculated by using

l m = l a -l g = π(D o -D i ) ln ( Do Di ) -l g , (7)
where l a , D o , and D i = 2r = 2b/π parameters are average core perimeter, outer, and inner core diameters, respectively. If the magnetic core is considered as n p -layered, then µ rg can be improved as

µ rgn = 1 + (µ rg -1) n p × d 1 T . (8)
The parameters, n p , d 1 = h/n p , and T are number of core layers, thickness of single core layer, and thickness of entire core, respectively. Maximum power delivered to R coil load resistance can be defined as [START_REF] Yuan | Magnetic field energy harvesting under overhead power lines[END_REF]

P CW EH = V coil /2 2 R coil . (9)
The load resistance, R coil , is calculated as [START_REF] Yuan | Magnetic field energy harvesting under overhead power lines[END_REF]

R coil = 2π(0.5D o )(N CW EH )ρ + πN 2 CW EH d 2 wire ρ L c , (10)
where ρ, d wire , and L c are resistivity, diameter of the winding wire, and length of the core, respectively. For CWEH, we adopt the parameters given in [START_REF] Bhuiyan | A miniature energy harvesting device for wireless sensors in electric power system[END_REF], [START_REF] Yuan | Magnetic field energy harvesting under overhead power lines[END_REF] which are listed as: 3) Free Standing Electromagnetic Harvester (FSEH): Free standing electromagnetic harvesters (FSEHs) are used in cases when CWEHs cannot be applied directly due to the distance constraints between the conductor and the harvesting coil. Since FSEHs cannot be clamped on top of the conductor, they can be placed near power cables, busbars and transformers. Hence, the gap assumption stated in Equation 6 is not valid [START_REF] Roscoe | Harvesting energy from magnetic fields to power condition monitoring sensors[END_REF]. According to the theoretical analysis provided in [START_REF] Roscoe | Harvesting energy from magnetic fields to power condition monitoring sensors[END_REF], the RMS voltage becomes

N CW EH =
V coil = 2πf µ ef f (N F SEH )π(0.5D) 2 B, (11)
where N F SEH is the number of turns in the coil, B is the magnitude of the magnetic flux density parallel to the coil axis, D is the diameter of the cylindrical core, and 1 < µ ef f < µ r is the effective relative permeability. The maximum power can be harvested by using FSEHs (i.e., P F SEH ) is calculated by using Equation 9. For FSEH we used the parameters in [START_REF] Tanaka | Measurement of power frequency electric and magnetic fields near power facilities in several countries[END_REF], [START_REF] Roscoe | Harvesting energy from magnetic fields to power condition monitoring sensors[END_REF] as N F SEH = 10K, D = D o = 0.05 m, µ ef f = 23 H/m, and B = 10 µT rms , respectively. In Figure 2 we present P CW EH and P F SEH harvested powers with respect to I (in Figure 2a) and B (in Figure 2b) as functions of n p (in Figure 2a) and N F SEH (in Figure 2b) considering the parameters given above. Incrementing the number of layers (n p ) in the magnetic core and number of turns in coils (N F SEH ) increases the power harvested. Moreover, increasing the RMS magnitude of the current (I) and flux density (B) also increases the amount of harvested power.

4) Hybrid Energy Harvesting Model: In our hybrid energy harvesting model, we assume that each critical power substation has a sensor node for monitoring. In addition, substations are covered with SEH panels as well as CWEH and FSEH devices. During an hour interval, the weather is cloudy with probability P cloud independent of other 1-hour time intervals. If the weather is not cloudy during the day time (i.e., 5:00 a.m. to 7:00 p.m.), all sensor nodes utilize harvested solar energy according to the profile given in Figure 1. Otherwise, if the weather is cloudy or during the night time, CWEH or FSEH are used to harvest energy with equal probability. Harvested powers are converted to electrical energy and transferred to the battery of each node.

D. MIP Framework for Energy Minimization

In this section, we present the MIP framework that aims to minimize the energy dissipation of the highest energy consuming node (i.e., a MinMax problem), since the lifetime of WSNs is usually limited by the lifetime of the the highest energy consuming node, by considering the hybrid energy harvesting strategy as well as transmission power control. The network topology is considered as a directed graph G = (V, A) with V represents all nodes including the base station (node-1). We also define set W which includes all nodes except node-1 (i.e., W = V \ {1}). A = {(i, j) : i ∈ W, j ∈ V -i} is the ordered set of arcs. Note that the set A implies that no node sends data to itself and the base station cannot generate data. The number of data packets flowing from node-i to node-j at time intervalt is represented with the integer variable f t ij . We define the set T = {1, 2, . . . , 24} to represent hours which is used to divide a single day to 24 hours (e.g., t = 1 represents 12:00 a.m. to 1:00 a.m., t = 2 represents 1:00 a.m. to 2:00 a.m. and so on). We represent each hour as 720 rounds (i.e., N = 720) where a round duration is taken as 5 seconds (i.e., τ rnd = 5 seconds). The MIP model is given in Figure 3. The objective function of the MIP model minimizes the highest energy dissipating node (i.e., E). Constraints of the MIP model are given in 12-20.

Constraint 12 ensures that incoming flows (i.e., j∈V f t ij ), generated flows at each round (i.e., N × s i ), and outgoing flows (i.e., j∈W f t ji ) are balanced at each sensor node (i.e., ∀i ∈ W ) in each time interval (i.e., ∀t ∈ T ). Constraint 13 enforces all sensor nodes to convey their collected data to the base station at each time interval-t. Constraint 14 is defined Minimize E Subject to:

j∈V f t ij - j∈W f t ji = N × s i ∀i ∈ W, ∀t ∈ T (12) j∈W f t j1 = N × j∈W s j ∀t ∈ T (13)
τ t b,i = T slot   (i,j)∈A λ lk ij f t ij + (j,i)∈A λ lk ji f t ji   + N × T DA ∀i ∈ W, ∀t ∈ T (14)
j∈V

E lk tx,ij f t ij tx + j∈W E lk rx,ji f t ji rx + P slp × (N τ rnd -τ t b,i ) sleep + N × E DA acquisiton ≤ e t i + P t h,i × (N τ rnd -τ t b,i ) harvest ∀i ∈ W, ∀t ∈ T (15)
t∈T

e t i ≤ E ∀i ∈ W (16) T slot   j∈V λ lk ij f t ij + j∈W λ lk ji f t ji + (j,n)∈A λ lk jn f t jn I i,lk jn   ≤ N × τ rnd ∀i ∈ V, ∀t ∈ T (17)
I i,lk jn = 1 if P ant rx,ji (l) ≥ P sns or P ant rx,ni (k) ≥ P sns 0 o.w.

(18) to calculate the busy time of node-i (i.e., τ t b,i ) including the effects of retransmissions for each time interval-t. Note that T DA = 5 ms is the time to acquire a data packet [START_REF] Kurt | Packet size optimization in wireless sensor networks for smart grid applications[END_REF]. Lefthand side of Constraint 15 is used to calculate the total amount energy consumed for transmission, reception, sleep, and data acquisition at each node-i for each time interval-t. E DA = 57 µJ is the energy consumed for data acquisition and P slp = 3 µW is the sleep power [START_REF] Kurt | Packet size optimization in wireless sensor networks for smart grid applications[END_REF]. Right-hand side of this constraint bounds the total amount of energy dissipation to e t i plus extra energy provided by the hybrid energy harvesting method. In this definition, e t i variable denotes the amount of energy consumed by node-i at time interval-t. We assume that energy harvesting is performed when node-i is in the sleep mode which lasts N τ rnd -τ t b,i seconds at each time interval-t. In this constraint P t h,i denotes the power harvested at node-i in time interval-t which is derived in the previous section. Constraint 16 is used to minimize the maximum energy consuming node in a single day (i.e., t∈T e t i ). Constraint 17 is used for each node-i to limit the aggregated duration of incoming, outgoing, and interfering flows to the total round duration of one hour (i.e., N × τ rnd ). Interference mechanism is stated in 18 and works as follows: if node-i interferes the two-way data/ACK communications between node-j and node-n by using power levels-l and k as long as the reception power is greater than the nominal receiver sensitivity level (i.e., P sns = -94 dBm), the interference function for node-i will be 1. Finally, Constraints 19 and 20 show the non-negativity conditions of the variables. We employ a link-level transmission power control scheme to determine the optimal power level pair for data packet and ACK packet (i.e., {l o ij , k o ji }) transmissions that minimizes the total energy dissipation in each link (i.e., E lk tx,ij + E lk rx,ji ).

f t ij ≥ 0 ∀(i, j) ∈ A, ∀t ∈ T (19)

III. ANALYSIS

In this section, we present the performance of our hybrid energy harvesting model on the energy consumption of sensor nodes. We construct the wireless channel model (Section II-A), link-layer model (Section II-B), and hybrid energy harvesting model (Section II-C) in MATLAB. The MIP model (Section II-D) is built and solved in GAMS by using CPLEX solver.

We consider an outdoor 500 kV substation power grid environment where |V | = 49 nodes are placed near to substations. Sensor nodes are placed to form a square (grid) topology where the base station is located at the center of the power grid. We note the distance between any node pairs in the same line as d int . As in [START_REF] Gungor | Opportunities and challenges of wireless sensor networks in smart grid[END_REF], we consider both LOS and NLOS cases while adopting n and σ values as 2.42 and 3.12 dB for LOS case; 3.51 and 2.95 dB for NLOS case. We choose d int values as 25, 27, 30 meters for LOS case; 8, 9, 10 meters for NLOS case to account for harsh SG environments. On each substation, SEH, CWEH, and FSEH devices are installed and harvested energy, according to the principle stated in Section II-C4, is transferred to the batteries of WSN nodes.

In Figure 4, we show our results by using the MIP model. We present the energy dissipation of the highest energy consuming node (E) with respect to the probability of cloudy weather (P cloud ) for LOS (in Figure 4a) and NLOS cases (in Figure 4b) considering the aforementioned d int values. We also present the results for non-harvesting cases (i.e., w/o harvesting) in this figure as dashed lines by substituting P t h,i with zero. Our results reveal that E does not change with the P cloud when considering w/o harvesting. We observe that highest energy consumption of a node can be in the interval 56.19-66.49 Joules and 54.99-66.06 Joules for LOS and NLOS cases, respectively. We also note that E increases as d int increases due to the raise of the transmission energy (i.e., E lk tx,ij ) for a sparse network. The proposed hybrid energy harvesting model (i.e., w/ harvesting) can reduce E by 47.29% to 91.46% and 31.27% to 87.54% for LOS and NLOS cases, respectively. It is quite revealing that amount of energy harvested is greatly increased (hence reducing E) as P cloud decreases. When solar energy is fully exploited (i.e., P cloud = 0), 75.88%-91.46% and 68.51%-87.54% of energy savings can be obtained for LOS and NLOS cases, respectively. On the other hand, when P cloud = 1, only CWEH and FSEH are utilized which would yield 47.29%-81.40% and 31.27%-72.81% of energy savings for LOS and NLOS cases, respectively.

IV. CONCLUSION

In this work, a hybrid energy harvesting model that exploits both solar and electromagnetic energies is proposed for WSN based SG applications. The proposed hybrid harvesting method considers the stochastic nature of weather conditions into the account. An MIP framework is designed on top of a log-normal shadowing channel as the physical layer and handshaking based link layer model that can employ transmission power control on each link. The MIP framework is used as a benchmark for the analysis and minimizes the energy dissipation of the highest energy consuming WSN node. Quantitative analysis is performed on a square WSN topology deployed for an outdoor 500 kV substation system for two different channel cases (LOS and NLOS). Our results reveal that developed hybrid energy harvesting framework can reduce the maximum energy consuming node up to 91.46% when solar energy is fully benefited. Nevertheless, when the solar energy cannot be exploited, the proposed hybrid harvesting method can still decrease the energy consumption up to 81.40%.

Fig. 2 :

 2 Fig. 2: Harvested powers for CWEH (P CW EH ) and FSEH (P F SEH ) methods.

  280, ω = 50 mm, f = 60 Hz, I = 13.5 A, l g = 1 mm, µ r = 10 5 H/m, d wire = 0.14 mm, d 1 = 0.1016 mm, T = 0.78 mm, D o = 16.3 mm, D i = 14.32 mm, l m = 47.3 mm, n p = 5, L c = 50 mm, and ρ = 1.11 Ω/m.

Fig. 3 :

 3 Fig. 3: MIP framework.

Fig. 4 :

 4 Fig. 4: Energy dissipation of the most energy consuming node (E in Joules) as a function of P cloud for the outdoor 500 kV substation LOS and NLOS environments w/ and w/o energy harvesting on a grid topology considering several d int values.

TABLE I :

 I Transmission power consumption (P crc tx (l) in mW) and output antenna power (P ant tx (l) in dBm) for each power level (l) for the CC2420 radio platform[15].

	l	P crc tx (l) P ant tx (l)	l	P crc tx (l)	P ant tx (l)
	3	25.5	-25	19	41.7	-5
	7	29.7	-15	23	45.6	-3
	11	33.6	-10	27	49.5	-1
	15	37.5	-7	31	52.2	0

ACKNOWLEDGEMENT

The work of V. C. Gungor was supported by the Turkish National Academy of Sciences Distinguished Young Scientist Award Program (TUBA-GEBIP) under Grand no. V.G./TUBA-GEBIP/2013-14.