Daniel Antunes
email: daniel.martins-antunes@univ-smb.fr

Jacques-Olivier Lachaud
email: jacques-olivier.lachaud@univ-smb.fr

Hugues Talbot
email: hugues.talbot@centralesupelec.fr

An Elastica-driven Digital Curve Evolution Model for Image Segmentation

Keywords: Multigrid convergence, Digital estimator, Curvature, Shape Optimization, Image Segmentation

Geometric priors have been shown to be useful in image segmentation to regularize results. For example, the classical Mumford-Shah functional uses region perimeter as prior. This has inspired much research in the last few decades, with classical approaches like the Rudin-Osher-Fatemi and most graph-cut formulations, which all use a weighted or binary perimeter prior. It has been observed that this prior is not suitable in many applications, for example for segmenting thin objects or some textures, which may have high perimeter/surface ratio. Mumford observed that an interesting prior for natural objects is the Euler Elastical model, which involves the squared curvature. In other areas of science, researchers have noticed that some physical binarization processes, like emulsion unmixing can be well-approximated by curvature-related flow like the Willmore flow. However, curvature-related flows are not easy to compute because curvature is difficult to estimate accurately, and the underlying optimisation processes are not convex. In this article, we propose to formulate a digital flow that approximates an Elastica-related flow using a multigrid-convergent curvature estimator, within a discrete variational framework. We also present an application of this model as a post-processing step to a segmentation framework.

Introduction

Geometric quantities are particularly useful as regularizers in low-level image analysis, especially when little prior information is known about the shape of interest. Length penalization is a well-behaved, general purpose regularizer and many models in the literature make use of it, from active contours [START_REF] Caselles | Geodesic active contours[END_REF] to level-set formulations [START_REF] Malladi | Image processing via level set curvature flow[END_REF][START_REF] Malladi | Shape modeling with front propagation: A level set approach[END_REF]. Discrete graph-based variational models have been particularly succesful to incorporate length penalization as a penalizer, while keeping the ability to extract a global optimum [START_REF] Boykov | Interactive graph cuts for optimal boundary amp; region segmentation of objects in n-d images[END_REF][START_REF] Appleton | Globally optimal geodesic active contours[END_REF].

However length regularization shows limitations when segmenting small or thin and elongated objects, as it tends to shrink solutions or yields disconnected solutions. Such drawbacks can be problematic in image segmentation, image restoration or image inpainting. It is thus rather natural to consider curvature (especially squared curvature) as a potential regularizer. Energies involving both length and squared curvature are often called the Elastica model (that dates back to Euler). The Willmore energy is its n-dimensional extension. These were brought to attention in computer vision by Mumford [31]. A similar concept is also present in the second-order regularizer of the original snake model [START_REF] Kass | Snakes: Active contour models[END_REF].Indeed, an explicit curvature term was often employed in early level-set methods, like [START_REF] Malladi | Image processing via level set curvature flow[END_REF][START_REF] Malladi | Shape modeling with front propagation: A level set approach[END_REF][START_REF] Ballester | Filling-in by joint interpolation of vector fields and gray levels[END_REF] for segmentation or inpainting. However, Casselles et al. in [START_REF] Caselles | Geodesic active contours[END_REF] observed that these terms as employed were not geometric, in the sense that they depend on the discretization parameters.

Existing works

The use of curvature for surface regularization in the case of the Willmore energy was studied in [START_REF] Bobenko | Discrete willmore flow[END_REF]. Although results look promising, authors used an inexact linearized version of the square curvature at every step.

One of the first successful uses of geometric, squared curvature in image processing is the inpainting algorithm described in [START_REF] Masnou | Level lines based disocclusion[END_REF], where authors evaluate the absolute curvature along the level lines of a simply-connected image to reconstruct its occluded parts. The non-intersection property of level lines induces the construction of an efficient dynamic programming algorithm. The curvature energy is nevertheless only coarsely approximated. We may quote [START_REF] Chan | Euler's elastica and curvature based inpaintings[END_REF] as another geometric inpainting method involving Elastica, which is optimized with level-sets.

In image segmentation, Zehiry and Grady have shown that injecting curvature as a regularizer can help recover thin and elongated objects [START_REF] El-Zehiry | Fast global optimization of curvature[END_REF]. Similarly Schoenemann et al. [START_REF] Schoenemann | Curvature regularity for region-based image segmentation and inpainting: A linear programming relaxation[END_REF] have considered the Elastica energy in their ratiobased image segmentation approach. In [START_REF] Schoenemann | The elastic ratio: Introducing curvature into ratio-based image segmentation[END_REF], Schoenemann et al. extend ratiocut approaches to segmentation with a curvature term. Their global optimization framework shows the considerable advantages of using such regularizer in common binary image segmentation tasks. However the time complexity of their algorithm is prohibitive, even if curvature is again coarsely approximated. In a more recent work, Lim et al. propose an edge-weighted Willmore-like flow to segment vertebra [START_REF] Lim | Introducing willmore flow into level set segmentation of spinal vertebrae[END_REF] using a level-set formulation. They observe improved results compared to [START_REF] Caselles | Geodesic active contours[END_REF] and other perimeter-based priors, but observe that a shape prior and a close-enough initialization must also be provided to achieve good results.

In fact, it is still a very challenging task to efficiently handle curvature in the context of image segmentation. State-of-art methods are in practice difficult to optimize and do not scale easily [START_REF] El-Zehiry | Fast global optimization of curvature[END_REF][START_REF] Schoenemann | Curvature regularity for region-based image segmentation and inpainting: A linear programming relaxation[END_REF][START_REF] Strandmark | Curvature regularization for curves and surfaces in a global optimization framework[END_REF][START_REF] Nieuwenhuis | Efficient squared curvature[END_REF]. In order to achieve reasonable running times, such approaches make use of coarse curvature estimations for which the approximation error is often unknown. Improving the quality of the curvature estimator has an important impact on the accuracy of the results, but is computationally too costly in these methods.

Segmentation methods using Elastica or Willmore energies can be formulated via the computation of the corresponding flow, i.e. following its gradient. Well-founded level-set methods involving Willmore energy [START_REF] Droske | A level set formulation for willmore flow[END_REF] have been used in medical image segmentation [START_REF] Lim | Introducing willmore flow into level set segmentation of spinal vertebrae[END_REF]. Another level-set formulation similar to the Chan-Vese model was proposed in [START_REF] Zhu | Image segmentation using euler's elastica as the regularization[END_REF]. Willmore flow can also be carried out by phase-field models (see [START_REF] Bretin | Phase-field approximations of the willmore functional and flow[END_REF]) but they are less suited to image segmentation since interfaces are blurred in such models. Threshold dynamics can also be considered for such energies [START_REF] Esedoglu | Threshold dynamics for high order geometric motions[END_REF] and has been proposed for image desocclusion [START_REF] Esedoglu | Threshold dynamics for shape reconstruction and disocclusion[END_REF]. However, flow-based methods will always deliver a local minimum, which heavily depend on initialization for good results. A similar approach is to formulate the Elastica as a regularization in a variational framework and use operator splitting to optimize it. Tai et al. [START_REF] Tai | A fast algorithm for euler's elastica model using augmented lagrangian method[END_REF] and more recently Deng et al. [START_REF] Deng | A new operator splitting method for the euler elastica model for image smoothing[END_REF] have proposed this approach, which is well-suited to inverse problem solving. The problem remains non-convex and the nature of the converged iterates are not studied, i.e. how close they are to a global optimum.

In contrast, researchers have sought convex relaxation of curvature-related formulation. In [START_REF] Goldluecke | Introducing total curvature for image processing[END_REF], Goldluecke and Cremers introduce the Total Curvature, based on the Menger-Melnikov curvature of a measure. The formulation is non-convex but can be relaxed into a well-correlated but non-tight convex one. The relaxed formulation can be used to solve inverse problems including segmentation. However computing a solution is still expensive, requiring signifiant parallelisation efforts with a GPU to achieve acceptable computing times. In [START_REF] Bredies | A convex, lower semicontinuous approximation of euler's elastica energy[END_REF], Bredies et al. study a convex approximation of the Elastica energy, using a lifting scheme. Thanks to the lifting, it is not restricted to binary sets. It shows promising results and can be used for image restoration in addition to inpainting and segmentation. The scheme is however neither tight nor exact, and required 4 dimensions for a 2D approximation. In addition, it does not carry over to nD. It also remains computationally expensive.

Discrete approaches have been used to tackle the Elastica difficulties. As noted before, [START_REF] Schoenemann | Curvature regularity for region-based image segmentation and inpainting: A linear programming relaxation[END_REF] have linearized a mesh-based approach, which was expanded upon by [START_REF] Strandmark | Curvature regularization for curves and surfaces in a global optimization framework[END_REF], while [START_REF] El-Zehiry | Fast global optimization of curvature[END_REF] have used a binary solver for a non-submodular approach on a regular square grid of pixels, which was expanded in [START_REF] El-Zehiry | Contrast driven Elastica for image segmentation[END_REF]. In [START_REF] Nieuwenhuis | Efficient squared curvature[END_REF] , a specific solver is developed to yield a better solution to an extension of the formulation in [START_REF] El-Zehiry | Fast global optimization of curvature[END_REF].

Motivation

In spite of the multiplicity of existing approaches, it seems none is currently completely satisfactory. While the sought solution is expected to be very smooth based on geometrical expectations, the mathematical properties of the related operators make this smoothness elusive. Non-convexity and the presence of high-order terms make Elastica difficult to optimize and to approximate, let alone compute efficiently. The computed solutions in the literature are rarely evaluated, and the complexity of the proposed solution is often high.

Consequently, we are interested in studying purely discrete variational segmentation models involving an Elastica energy. Recently, new estimators of curvature along digital contours have been proposed [START_REF] Roussillon | Accurate curvature estimation along digital contours with maximal digital circular arcs[END_REF][START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF][START_REF] Schindele | Multigrid convergence for the MDCA curvature estimator[END_REF] with the desirable multigrid convergence property. This motivates us to propose models in which they can be used successfully.

In this paper, we propose to investigate the use of the digital integral invariant curvature estimator [START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF] in a discrete variational segmentation framework. More precisely, we show how to incorporate it in a digital flow minimizing its squared curvature. This model is expressible as a discrete combinatorial optimization model, which we solve using a classical variant of QPBO algorithm [START_REF] Rother | Optimizing binary mrfs via extended roof duality[END_REF]. Our solution does not compute a global solution to the whole domain but a global solution to a band surrounding a given shape. Our approach thus leads to a digital flow similar to continuous Elastica, as shown by our experiments. We present an application of this model as a post-processing step in a segmentation framework, and demonstrate how it can improve standard results. This paper is an extended version of a previous work that appeared in [START_REF] Antunes | Digital curvature evolution model for image segmentation[END_REF]. We have included a deeper discussion on our optimization method and we present a more thorough comparison of our segmentation method with two other related works. Moreover, we have included a new local combinatorial optimization approach for elastica minimization, which show that our discrete model is well-founded. Finally, the code produced in this paper is freely available on github. 1

Outline

Section 2 reviews the concept of multigrid convergence and highlights its importance for the definition of digital estimators. We describe two convergent estimators used in this paper, one that approaches the length of elementary discrete contour elements, and the other that approaches the curvature of a discrete contour. They are used in the optimization model and in the definition of the digital elastica. Section 3 describes a local combinatorial optimization model suitable for several curvature estimators, that presents interesting results but is too costly in practice. Section 4 describes the proposed Elasticadriven curvature evolution model along with several illustrations of digital flows. Section 5 explains how to use this evolution model as a post-processing step in an image segmentation framework and compares it to two related approaches. We conclude and draw some perspectives to this work in Section 6.

Multigrid convergent estimators

Our objective is to delineate shapes within digital images with some priors related to continuous geometry. The goal of this section is to introduce the concept of multigrid convergence and its potential when analyzing digital images with variational models involving geometric quantities, with the constraint that only digitized shapes are observed.

A digital shape is the result of some quantization process over an object X lying in some continuous space of dimension 2 (here). For example, the Gauss digitization of a continuous subset X of the Euclidean plane R 2 with grid step h > 0 is defined as

D h (X) = X ∩ (hZ) 2 .
Given a shape X and its digitization D h (X), a digital estimator û for some geometric quantity u is intended to compute u(X) by using only the digitization. This problem is not well-posed, as the same digital object could be the digitization of infinitely many objects very different from X. Therefore, a characterization of what constitutes a good estimator is necessary.

Let u be some geometric quantity of X (e.g. tangent, curvature). We wish to devise a digital estimator û for u. It is reasonable to state that û is a good estimator if û(D h (X)) converges to u(X) as we refine our grid. For example, counting pixels is a convergent estimator for area (with a rescale of h 2 ); but counting boundary pixels (with a rescale of h) is not a convergent estimator for perimeter. Multigrid convergence is the mathematical tool that makes this definition precise. Given any subset Z of (hZ) 2 , we can represent it as a union of axis-aligned squares with edge length h centered on the point of Z. The topological boundary of this union of cubes is called h-frontier of Z. When Z = D h (X), we call it h-boundary of X and denote it by ∂ h X.

Definition 1 (Multigrid convergence for local geometric quantites)

A local discrete geometric estimator û of some geometric quantity u is (uniformly) multigrid convergent for the family X if and only if, for any X ∈ X, there exists a grid step h X > 0 such that the estimate û(D h (X), x, h) is defined for all x ∈ ∂ h X with 0 < h < h X , and for any x ∈ ∂X,

∀x ∈ ∂ h X with x -x ∞ ≤ h, û(D h (X), x, h) -u(X, x) ≤ τ X (h),
where τ X : R + \ {0} → R + has null limit at 0. This function defines the speed of convergence of û towards u for X.

For a global geometric quantity (e.g. perimeter, area, volume), the definition remains the same, except that the mapping between ∂X and ∂ h X is no longer necessary.

Multigrid convergent estimators provide a quality guarantee and should be preferred over non-multigrid convergent ones. In the next section, we describe two estimators that are important for our purpose.

Tangent and Perimeter Estimators

The literature presents several perimeter estimators that are multigrid convergent (see [START_REF] Coeurjolly | A comparative evaluation of length estimators of digital curves[END_REF][START_REF] Coeurjolly | Multigrid Convergence of Discrete Geometric Estimators[END_REF] for a review), but in order to define the digital Elastica we need a local estimation of length and we wish that integration over these local length elements gives a multigrid convergent estimator for the perimeter.

Definition 2 (Elementary Length) Let a digital curve C be represented as a sequence of grid vertices in a grid cell representation of digital objects (in a grid with step h). Further, let v be a multigrid convergent estimator for the unit tangent vector. The elementary length ŝ(e) at some oriented grid edge e ∈ C is defined as

ŝ(e) = h v(e) • e.
The integration of the elementary length along the digital curve is a multigrid convergent estimator for perimeter if one uses the λ-MST [START_REF] Lachaud | Fast, accurate and convergent tangent estimation on digital contours[END_REF] tangent estimator (see [START_REF] Lachaud | Non-Euclidean spaces and image analysis : Riemannian and discrete deformable models, discrete topology and geometry[END_REF]).

Integral Invariant Curvature Estimator

Generally speaking, an invariant is a function whose value is unaffected by the action of some group on the elements of the domain. Perimeter and curvature are examples of invariants for shapes on R 2 with respect to the euclidean group of rigid transformations. Definition of integral area invariant and its link with curvature is proven in [START_REF] Manay | Integral invariant signatures[END_REF].

Definition 3 (Integral area invariant) Let X ⊂ R 2 and B r (p) the ball of radius r centered at point p. Further, let 1 X (•) be the characteristic function of X. The integral area invariant σ X,r (•) is defined as

∀p ∈ ∂X, σ X,r (p) = Br(p) 1 X (x)dx.
The value σ X,r (p) is the area of the intersection of the ball B r (p) with shape X. By approaching the shape at point p ∈ X, one can rewrite the intersection area σ X,r (p) in the form of the Taylor expansion [START_REF] Pottmann | Integral invariants for robust geometry processing[END_REF]:

σ X,r (p) = π 2 r 2 - κ(X, p) 3 r 3 + O(r 4 ),
where κ(X, p) is the curvature of X at point p. By isolating κ we can define a curvature estimator

κ(p) := 3 r 3 πr 2 2 -σ X,r (p) , (1)
Such an approximation is convenient as one can simply devise a multigrid convergent estimator for the area.

Definition 4 Given a digital shape D ⊂ (hZ) 2 , the area estimatir Area(D, h) is defined as

Area(D, h) := h 2 Card (D) . (2)
It is well known that this area estimator is multigrid convergent (e.g. see [START_REF] Klette | Multigrid convergence of calculated features in image analysis[END_REF]).

In [START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF], the authors combine the approximation(1) and digital estimator (2) to define a multigrid convergent estimator for the curvature.

Definition 5 (Integral Invariant Curvature Estimator) Let D ⊂ (hZ) 2 a digital shape. The integral invariant curvature estimator is defined as

κr (D, x, h) := 3 r 3 πr 2 2 -Area (B r (x) ∩ D, h) .
This estimator is multigrid convergent with speed O(h

3 ) for radii chosen as r = Θ(h 1 3 ). This estimator is also robust to noise and can be extended to estimate the mean curvature of three dimensional shapes.

Elastica energy estimator

In the remaining of this article, we are considering in the minimization of the Elastica energy. Given a Euclidean shape X with smooth enough boundary, the Elastica is defined as

E(X) = ∂X (α + βκ 2 )ds, for α ≥ 0, β ≥ 0. ( 3
)
The digital version of Elastica energy approaching the continuous Elastica energy is defined as follows (it uses multigrid convergent estimators):

Ê(D h (X)) = e∈∂D h (X) ŝ(e) α + βκ 2 r (D h (X), ė, h) , (4)
where ė denotes the center of the edge e. In the expression above, we will substitute an arbitrary subset Z of Z 2 to D h (X) since the continuous shape X is unknown. In the following we omit the grid step h to simplify expressions (or, putting it differently, we assume that the shape of interest is rescaled by 1/h and we set h = 1).

Local combinatorial optimization

The objectives of this section is to show both the potential of minimizing a purely digital Elastica energy but also to underline the difficulties of using it in a global combinatorial optimization framework. Given a digital shape S (0) we describe a process that generates a sequence S (i) of shapes with non-increasing Elastica energy. The idea is to define a neighborhood of shapes N (i) to the shape S (i) and choose the element of N (i) with lowest energy. The process is suited for the integral invariant estimator but also for other curvature estimators, for example, MDCA [START_REF] Roussillon | Accurate curvature estimation along digital contours with maximal digital circular arcs[END_REF]. As a matter of fact, our experiments have shown that either estimators induce similar results.

Let S be a 2-dimensional digital shape. We adopt the cellular-grid model to represent S, i.e., pixels and its lower dimensional counterparts, linels and pointels, are part of S. In particular, we denote by ∂S the topological boundary of S, i.e., the connected sequence of linels such that for each linel we have one of its incident pixels in S and the other not in S.

Let d S : Ω → R be the Euclidean distance transformation with respect to shape S. The value d S (x) gives the Euclidean distance between x and the closest pixel in S.

Definition 6 (m-Ring Set) Given a digital shape S ∈ Ω, its distance transformation d S and natural numbers m > 0, the m-ring set of S is defined as

R m (S) := R - m (S) ∪ R + m (S) := {x ∈ Ω | m -1 < d S (x) ≤ m} ∪ {x ∈ Ω | m -1 < d S (x) ≤ m} .
Consider the following set of neighbor candidates to S:

U(S) = {D |D ⊂ R 1 (S) ∪ S and D is connected}.
Such set can be extremely large and its complete exhaustion is prohibitively expensive. Instead, we explore a subset of it with the help of n-glued curves.

An oriented closed curve C is a closed connected sequence of linels with a well-defined interior. A segment of C is a connected subsequence c ∈ C of its linels.

Definition 7 (Glued Curve) Given closed curves C 1 , C 2 agreeing with some orientation q, a glued curve is a closed curve (c 1 , 1 , c 2 , 2 ) with orientation q and c 1 ∈ C 1 , c 2 ∈ C 2 . The linels 1 , 2 are called junction linels.

Definition 8 (n-Glued Curve Set) Given closed curves C 1 , C 2 with same orientation, its set of n-glued curves is defined as Let S O = (S ∪ R + 1 (S)) and S I = (S \ R - 1 (S)), the neighborhood set to shape S is defined as

G n (C 1 , C 2 ) = {(c 1 , 1 , c 2 , 2 ) | |c 2 | = n},
N (S, N ) = 1≤n≤N int G n (∂S O , ∂S) ∪ int G n (∂S, ∂S I ) ,
where int(C) is the interior of the shape bounded by the close oriented curve C. Algorithm 11 describes the local combinatorial process and Figure 2 presents the digital curve evolution when executing this algorithm for two different shapes with N = 20. input : A digital set S; the maximum length of glued curves N ; the maximum number of iterations maxIt; and a stop condition tolerance 1 delta ←-tolerance +1; 2 while i < maxIt and delta > tolerance do

3 for X ∈ N (S (i) , N ) do 4 if Ê(X) < Ê(X ) then 5 X ←-X 6 end 7 end 8 i ←-i +1; 9 S (i) ←-X ; 10 delta ←-Ê(S (i-1) ) -Ê(S (i) ); 11 end
Algorithm 1: Local combinatorial optimization for elastica minimization.

The running time of algorithm 11 is summarized in table 1. All the experiments in this paper were executed on a five-core 3.4Ghz CPU and the number of pixels in the triangle, square and flower shapes are respectively 841, 1867, 521 for grid step h = 1.0. Although its use in practical applications is limited, we demonstrate that digital estimators are effective in their measurements and the flows evolve as expected. We observe that it is a complete digital approach, and we do not suffer from discretization and rounding problems, a common issue in continuous models. Furthermore we have checked that this approach works indifferently with Integral Invariant curvature estimator and Maximal Digital Circular Arc curvature estimator. So the convergence of the digital curvature estimator seems to be the cornerstone to get a digital curve behaving like a continuous Elastica. In the next section, we explore a more efficient approach.

Digital Curvature Flow

Given a digital shape S ∈ Ω, the idea is to evolve a flow S (i) where S (i+1) is derived from the minimization of the sum of the squared curvature at S (i) . In order to guarantee connectivity, we limit the optimization region to a subset of Ω, namely the inner pixel boundary of S (i) . We also use the notation |S| to denote the cardinal of a digital set S.

Discrete variational model for discrete curve evolution

We start by definition the optimization variables.

Definition 9 (Inner pixel boundary) Given a digital shape S embedded in a domain Ω, we define its inner pixel boundary set O(S) as

O(S) = { x | x ∈ S, |N 4 (x) ∩ S| < 4 } ,
where N 4 (x) denotes the 4-adjacent neighbor set of x (without x).

In other words, each flow iteration decides which pixels in the inner boundary are to be removed and which are to be kept. To simplify notation, we denote the optimization region of S (i) as O (i) , and the foreground set as i) . We recall the definition of our digital curvature estimator :

F (i) = S (i) \ O (
κ2 (y) = c 1 c 2 -|B r (y) ∩ S| 2 , ( 5
)
where c 1 = 3/r 6 and c 2 = πr 2 /2. We define O

x j + |F (i) r (y)| 2 = c 1 • C + 2 |F (i) r (y)| -c 2 xj ∈O (i)
r (y)

x j + xj ∈O (i) r (y) x 2 j + xj ,x k ∈O (i) r (y) j<k 2x j x k , where C = c 2 2 -2c 2 • |F (i) r (y)| + |F (i)
r (y)| 2 is a constant. By ignoring constants and multiplication factors and using the binary character of the variables, we define the following family of energies for given parameters α, β, γ ≥ 0.

E m (X, S) = x∈X αs(x) + y∈Rm(S) β (1/2 + |F (i) r (y)| -c 2 ) • xj ∈O (i) r (y) x j + xj ,x k ∈O (i) r (y) j<k x j x k + x∈X γg(S, x), (6)
where g(S, x) denotes a similarity term as distance to initial shape S or data fidelity term and s(x) denotes a length penalization term. Each choice of m generates a different flow, which is generaly described in the Digital Curve Evolution (DCE) algorithm 15.

input : A digital set S; The ring number m; Length(α), curvature(β) and similarity(γ) coefficients; the maximum number of iterations maxIt;

1 S (i) ←-S;

2 while i < maxIt do 3 //Expansion mode 4 if i is odd then 5 X (i) ←-arg min X E m (X, S (i) ); 6 S (i) ←-F (i) + X (i) ; 7 S (i) ←-S (i) ; 8 end 9 //Shrinking mode 10 else 11
X (i) ←-arg min X E m (X, S (i) ); 12
S (i) ←-F (i) + X (i) ; 13 end 14 i ←-i +1;
15 end Algorithm 2: Digital curvature evolution algorithm (DCE). Recall that the integral invariant estimator approaches curvature by computing the difference between half of the area of a chosen ball and the area of the intersection of this ball with the shape. In particular, regions of positive curvature have fewer pixels in their intersection set than on its complement w.r.t the estimation ball. This implies that variables in such regions are labeled with 1, as the unbalance grows otherwise. We attenuate curvature if we shift the center of the estimation ball towards the interior of the shape, which means removing the 1-labeled pixels. That is why we take the complement of the optimization solution.

The same reasoning applies for non-convex parts. Indeed, concave regions are convex in the shape complement. In the expansion mode we apply the same reasoning on the image complement, and by doing this we are able to handle concavities. It is called expansion mode because the optimization region, in this case, is the outer pixel boundary of the original shape. Table 2 sums up these arguments.

Length and similarity terms should be properly defined in order to comply with the complement step of the DCE algorithm. The length penalization is defined as

s(x) = xj ∈N4(x) t(x j ), where t(x j ) =      (x -x j ) 2 , if x j ∈ O (i) r (x -0), if x j ∈ F (i) r (x -1), otherwise (7)
shrink mode κ 0 κ ≥ 0 κ < 0 X x k = 1 x k ∈ {0, 1} x k = 0 S (i+1) ← S (i) \ X eroded prob. eroded unchanged expansion mode κ 0 κ ≥ 0 κ < 0 X xk = 1 xk ∈ {0, 1} xk = 0 S (i+1) ← S(i) \ X dilated prob. dilated unchanged
Table 2: Since the curvature is negated when reversing the curve (i.e. κ = -κ), this process can only shrink convex parts in shrink mode and expand concave parts in expansion mode.

We do not use similarity terms in this section, thus we postpone the definition of such terms until later in this article, with the description of how DCE can be used in an image segmentation framework.

In figure 4 we see some results for the DCE algorithm with m = 1. We observe a global evolution towards rounder shapes, but with several artifacts. We minimize the effects of a jaggy boundary by setting α > 0. Nonetheless, a higher estimation ball radius creates unstable shapes. In fact, the estimator is very sensitive in regions of low squared curvature, and it is precisely in those regions that spurious pixels are created. In the previous section we noticed that the algorithm produces shapes with many artifacts due to the small uncertainties of the estimator along regions of low squared curvature. We argue that, evaluating the estimation ball along outer ring sets we avoid those sensitive areas by focusing the optimization process only on regions with highest squared curvature value. In our experiments, the best results are obtained by executing DCE algorithm with m equal to r, where r is the estimation ball radius (see figure 5 ). We observe that Elastica may increase after some iterations if chosen radius is too large, as in the case of the triangle in figure 6 in which the flow converges to a single point. We conjecture that an appropriated value for the radius should be given by the shape reach. The produced flow has no difficulties in handling changes on topology, and it presents different speeds for regions with low and high curvature values, as illustrated in figure 7.

Optimization method

Let f be a function of n binary variables with unary and pairwise terms, i.e.

f (y 1 , • • • , y n ) = j f j (y j ) + j<k f j,k (y j , y k ).
Function f is submodular if and only if the following inequality holds for each pairwise term f j,k [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF]:

f j,k (0, 0) + f j,k (1, 1) ≤ f j,k (0, 1) + f j,k (1, 0).
Energy E m is non-submodular and optimizing it is a difficult problem, which constrains us to use heuristics and approximation algorithms. The QPBO method [START_REF] Rother | Optimizing binary mrfs via extended roof duality[END_REF] transforms the original problem in a max-flow/min-cut formulation and yields a full optimal labeling for submodular energies. For nonsubmodular energies the method is guaranteed to return a partial labeling with the property that the set of labeled variables is part of an optimal solution. That property is called partial optimality.

In practice, QPBO can leave many pixels unlabeled. There exist two extensions to QPBO that alleviate this limitation: QPBOI (improve) and QPBOP (probe). The first is an approximation method that is guaranteed to not increase the energy, but loses the property of partial optimality. The second is an exact method which is reported to label more variables than QPBO. The percentage of unlabeled pixels by QPBOP for E 1 is quite high, but the percentage decreases to zero as we set m equal to r. Therefore, we are more confident in taking the solution for values of m close to r. However, the way it varies across m values differs from shape to shape, as is illustrated in figure 9. We also noticed that, for m = r, all the pixels were labeled, which may indicate that E r is an an easy instance of the general non-submodular energy E m , but this remains to be proved. The number of pairwise terms in E r is roughly half of those in E 1 (see figure 8).

We have used QPBOI to solve E m . Naturally, in the case where all pixels are labeled by QPBOP, QPBOI returns the same labeling as QPBOP.

Application to image segmentation

We present an application of our digital curve evolution algorithm to supervised image segmentation. The DCE acts as a contour correction method. Here we use a data fidelity term in order to characterize the object of interest. The highest ring has rougly half the number of pairwise terms than the lowest ring.

Given foreground and background seeds selected by the user, we derive mixed gaussian distributions of color intensities G f and G b , and we define the data fidelity term as the cross-entropy, i.e.

g(y) = -(1 -y) log G f (I(y)) -y log G b (I(y)). (8)
We use the DCE algorithm to regularize an initial contour outputed by some segmentation algorithm or delineated by the user. In this application, the similarity term of the DCE is set to the data fidelity term (8). 9 end Algorithm 3: Contour correction algorithm. The algorithm can be initialized by a collection of compact sets, or with the result of a third-party segmentation algorithm, as GrabCut [START_REF] Rother | grabcut": Interactive foreground extraction using iterated graph cuts[END_REF]. We include an additional parameter d that dilates the initial sets using a square of side one before executing the flow. Fig. 12 tations and demonstrates the completion property of curvature. Moreover, it does not suffer from over segmentation, it is much faster than Schoenemann's, and in several opportunities produces better segmentations than GrabCut.

Conclusion

We have studied in depth several digital curve evolution models based on a digital version of the Elastica energy and we have presented an application to image segmentation. The processes we have described are completely digital, and do not suffer from issues that typically arise in models that passes through a discretization stage, such as rounding. Moreover, the model can handle changes in topology and its results are competitive with similar approaches while achieving reasonable running times. Future developments of this work will include extending the DCE algorithm to 3d, as the integral invariant estimator is available in higher dimensions. As a result, a denoising application may be derived by lifting the flow to 3d with a data fidelity defined by image intensities or by distance to initial position.

Fig. 1 :

 1 Fig. 1: The blue pixels in figure (a) illustrates a 2-ring set, while the green pixels its inner-pixel boundary. In figure (b) we highlight an element of set G 11 (C 1 , C 2 ). The curve segments of C 1 and C 2 are colored in red and green, while the junction linels in blue.

5 Fig. 2 :

 52 Fig. 2: Local combinatorial optimization process results for the square and flower shapes. Shapes displayed at every 2 iterations for the square and 7 iterations for the flower.

Fig. 3 :

 3 Fig. 3: Elastica computation along iterations of algorithm 11

r

  (y) := O (i) ∩ B r (y) and an analogous definition holds for F (i) r (y). Expanding (5), we get κ2 (y) = c 1 c 2xj ∈Br(y)

15 Fig. 4 :

 154 Fig.4: The algorithm is very sensitive to the little variations of the estimator, which are particularly important in regions of low squared curvature. Artifacts are somewhat reduced with a length penalization but increases if we use a higher ball radius.

(a) r = 5 , 3 Fig. 5 :

 535 Fig. 5: By positioning the estimation ball on farther rings, we minimize artifacts creation.

Fig. 6 :

 6 Fig.6: The choice of radius impacts the flow. In the figures, the flow ceases to evolve for all shapes when r = 3 (a). In figure (b), for r = 5, the triangle evolves to a single point, while the others stop in an intermediate shape, as in (a). In figure (c), we observe that for a given choice of radius, the elastica may increase after a certain number of iterations.

Fig. 7 :

 7 Fig. 7: High curvature regions evolves faster than lower ones (a). The flow can handle topological changes (b).

Fig. 8 :

 8 Fig. 8: We plot the ratio of pairwise terms among all |X (i) | 2 combinations.The highest ring has rougly half the number of pairwise terms than the lowest ring.

Fig. 9 : 6 S 8 i

 968 Fig.9: For each plot, we first produce shapes S (i) executing DCE with m = r. Then, for each shape in S (i) , we execute one iteration of DCE for different values of m and we count the unlabeled pixels. The number of unlabeled pixels by QPBOP remains high for lower values of m, and goes to zero when m = r. We observe the same behaviour for varying radius values.

Fig. 10 :

 10 Fig. 10: Comparison of squared curvature regularization (first row) and length regularization (second row).

  5, α = 0.1, β = 1.0, γ = 3.0) (α = 0.1, β = 1.0)

Table 1 :

 1

Running times for the local combinatorial optimization algorithm with k = 1, N = 20. Four threads were used.

https://www.github.com/danoan/BTools

This work has been partially funded by CoMeDiC ANR-15-CE40-0006 research grant.

GrabCut

Contour Correction Schoenemann (r = 5, α = 0.1, β = 1.0, γ = 3.0) (α = 0.1, β = 1.0) Fig. 11: The proposed method regularizes GrabCut contours and returns meaningful results. We can observe the completion feature of curvature in the second row, and we don't suffer from oversegmentation issues as Schoenemann's method. However, our flow may stop in a local optima as in the fourth row, while Schoenemann's is able to extrapolate such solutions.

We evaluate our method using the BSD300 database [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. All images contains the same number of pixels, the resolution being 321x481 in portrait mode. We compare the results of our method with segmentations given by GrabCut and Schoenemanns's method [START_REF] Schoenemann | Curvature regularity for region-based image segmentation and inpainting: A linear programming relaxation[END_REF]. We report an average of 3s per flow iteration, and an average of 30 iterations per image. While GrabCut executes in less than one second, Schoenemann's method may take several hours to complete.

In figure 10 we can observe the results of a curvature regularization in comparison with a pure length regularization. The curvature can fill gaps and is smoother than the one produced by length only, resulting in more pleasant segmentations. In figures 11 and 12 we list several results of our method and we compare them with segmentation produced by GrabCut and Schoenemann's method. Our method is successful in producing curvature regularized segmen-


Minh Trinh-Hoang

Mohammed Nabil

El Korso

Marius Pesavento

A PARTIALLY-RELAXED ROBUST DOA ESTIMATOR UNDER NON-GAUSSIAN LOW-RANK INTERFERENCE AND NOISE

Keywords: DOA Estimation, EM Algorithm, Robust Estimation, Partial Relaxation

In practical applications, non-Gaussianity of the signal at the sensor array is detrimental to the performance of conventional Direction-of-Arrival (DOA) estimators developed under the Gaussian model. In this paper, we propose a novel robust DOA estimator from the data collected at the sensor array under the corruption of non-Gaussian interference and noise. Additionally, the Cramér-Rao bound for DOA parameters under the considered signal model is derived. Simulation results show that the proposed estimator exhibits near-optimal estimation performance under the assumed model while being robust to model mismatch and/or the presence of outliers.

INTRODUCTION

Recently, a class of Direction-of-Arrival (DOA) estimators have been proposed under the Partial Relaxation (PR) framework [START_REF] Trinh-Hoang | Partial Relaxation Approach: An Eigenvalue-Based DOA Estimator Framework[END_REF][START_REF] Trinh-Hoang | Cramér-Rao Bound for DOA Estimators Under the Partial Relaxation Framework: Derivation and Comparison[END_REF]. Similar to the MUSIC algorithm [START_REF] Schmidt | Multiple Emitter Location and Signal Parameter Estimation[END_REF], PR estimators employ a one-dimensional spectral search to estimate the DOAs. Simulation results show that, under the Gaussian signal model, the threshold performance of PR estimators is comparable to that of conventional maximum likelihood estimators [START_REF] Böhme | Estimation of Source Parameters by Maximum Likelihood and Nonlinear Regression[END_REF][START_REF] Viberg | Sensor Array Processing Based on Subspace Fitting[END_REF][START_REF] Stoica | Performance Study of Conditional and Unconditional Direction-of-Arrival Estimation[END_REF] while the computational complexity is in the same order of magnitude as the MUSIC algorithm.

As the derivation of PR estimators in [START_REF] Trinh-Hoang | Partial Relaxation Approach: An Eigenvalue-Based DOA Estimator Framework[END_REF] relies on the assumption of Gaussian source signals and spatially white Gaussian noise, the DOA estimation performance of the PR methods degrades in practical scenarios where the receive signal at the sensor array is not Gaussian-distributed. For example, in MIMO radar systems, the receive signal is corrupted by radar clutter, which is generally modeled as realizations of some stochastic processes with heavy-tailed distributions [START_REF] Ollila | Complex Elliptically Symmetric Distributions: Survey, New Results and Applications[END_REF]. In the context of radio astronomical imaging, aside from thermal noise, low-rank terrestrial interferers [START_REF] Abdallah | Bayesian Signal Subspace Estimation with Compound Gaussian Sources[END_REF] are imminent, which adversely impact the imaging quality, especially for subspace-based methods in general and estimators under the PR framework in particular [START_REF] Leshem | Radio-Astronomical Imaging in the Presence of Strong Radio Interference[END_REF][START_REF] Besson | Directionof-Arrival Estimation in a Mixture of K-distributed and Gaussian Noise[END_REF]. Therefore, deriving robust DOA estimators under consideration of the non-Gaussianity of both interference and noise [START_REF] Rangaswamy | Robust Adaptive Signal Processing Methods for Heterogeneous Radar Clutter Scenarios[END_REF][START_REF] Ginolhac | Performance of Two Low-Rank STAP Filters in a Heterogeneous Noise[END_REF][START_REF] Ollier | Robust Distributed Calibration of Radio Interferometers with Direction Dependent Distortions[END_REF] is of great interest. In this paper, we propose a robust DOA estimator under the consideration of non-Gaussian observations which does not rely on the robust estimation of the covariance matrix as in [START_REF] Abdallah | Bayesian Signal Subspace Estimation with Compound Gaussian Sources[END_REF][START_REF] Huber | Robust Statistics, Wiley Series in Probability and Statistics[END_REF][START_REF] Williams | Robust Estimation of Structured Covariance Matrices[END_REF][START_REF] Couillet | Robust Spiked Random Matrices and a Robust G-MUSIC Estimator[END_REF][START_REF] Breloy | Robust Covariance Matrix Estimation in Heterogeneous Low Rank Context[END_REF]. The receive signal at the sensor array is assumed to be corrupted by Student t-distributed interference and noise. The proposed algorithm, which is based on the Expectation Maximization (EM) algorithm [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF] and the PR framework, jointly estimates the non-Gaussianity of the receive signal and the DOA parameters in each iteration. From the simulation results, the proposed algorithm exhibits favorable performance even in the case of model mismatch. The paper is organized as follows. The signal model is presented in Section 2. In Section 3, the Cramér-Rao bound (CRB) for the DOA parameters is derived in closed-form. Details on the derivation of the proposed DOA estimator is provided in Section 4. Numerical results in Section 5 illustrates the performance of the proposed algorithm under different scenarios. Lastly in Section 6, remarks and extensions for future research are discussed.

SIGNAL MODEL

Consider an array of M sensors receiving N narrowband signals emitted from sources with the corresponding unknown DOAs θ = [θ 1 , . . . , θ N ]

T and under corruption of R interferers. The baseband receive signal x(t) = [x 1 (t), . . . , x M (t)] T ∈ C M at the time instant t is modeled as

x(t) = A(θ)s(t) + Bj(t) + n(t) for t = 1, . . . , T. (1
) In (1), the steering matrix A(θ) = [a(θ 1 ), . . . , a(θ N )] is a full rank matrix whose n-th column a(θ n ) represents the sensor array response for the DOA θ n . B ∈ C M ×R denotes the unstructured loading matrix of R interferers. In addition, s(t) ∈ C N , j(t) ∈ C R and n(t) ∈ C M denote the baseband source signal vector, the baseband interference vector, and the noise vector at the sensor array, respectively. While the source signal s(t) is assumed to be deterministic but unknown, the interference vector j(t) and the noise vector n(t) admit the stochastic representations:

j(t) d = 1 τ (t) d(t) and n(t) d = 1 τ (t) e(t). (2)
In (2), the speckles d(t) and e(t) are i.i.d. vectors drawn from d(t) ∼ CN (0, I) and e(t) ∼ CN 0, σ 2 I , respectively.

The noise power at the sensor array is denoted by σ2 . The textures τ (t) ∼ Gamma (ν, ν)1 for t = 1, . . . , T are i.i.d. random scalars. The positive number of degrees of freedom (DOF) is denoted by ν, which models the non-Gaussianity of the inteferers and the noise. Equation (1) can be rewritten for multiple snapshots t = 1, . . . , T in a compact notation as

X = AS + BJ + N , (3)
where X = [x(1), . . . , x(T )] ∈ C M ×T is the baseband receive signal matrix. In a similar manner, we define the source signal matrix S ∈ C N ×T , the interfering signal matrix J ∈ C R×T and the sensor noise matrix N ∈ C M ×T as S = [s(1), . . . , s(T )], J = [j(1), . . . , j(T )] and N = [n(1), . . . , n(T )], respectively.

CRAM ÉR-RAO BOUND

The derivation of the CRB for the DOA parameters θ under the signal model in ( 1) is outlined in this section. By marginalizing w.r.t. the texture τ (t), we can show that the observation vectors x(t) for t = 1, . . . , T are i.i.d. drawn from a circularly complex Student's t-distribution x(t) ∼ Ct ν (A(θ)s(t), Σ) with the scale matrix Σ = BB H + σ 2 I. As Σ does not depend on the DOA parameters θ, the CRB for θ under the Student's t-Distribution case is proportional to that under the Gaussian case with the mean A(θ)s(t) and the covariance matrix Σ [START_REF] Besson | On the Fisher Information Matrix for Multivariate Elliptically Contoured Distributions[END_REF]. Therefore, combining [19, Eq. ( 22)] with [20, Eq. ( 31)], the CRB for the DOA parameters θ is given by

C(θ) = ν + M + 1 2(ν + M )T Re DH Π ⊥ A D P s -1 , (4)
with

D = Σ -1/2 ∂a(θ1) ∂θ , . . . , ∂a(θ N ) ∂θ , Π ⊥ A = I -AA † and P s = T t=1 s(t)s(t) H .

PROPOSAL OF THE EM-PR-DML ALGORITHM

In this section, we propose a novel maximum likelihoodbased DOA estimator under the signal model in Section 2.

The parameter vector φ which characterizes the density function of the observation vector x(t) is given by

φ = θ T , s(1) T , . . . , s(T ) T , ζ T , ν, σ 2 T , (5)
where ζ contains the real and imaginary part of non-redundant entries of Σ. The log-likelihood L = log p(X|φ) is given by

L = -T log det Σ -(M + ν) T t=1 log 1 + 1 ν r(t) H Σ -1 r(t) (6
) with the residuals r(t) = x(t) -A(θ)s(t). As the direct maximization of the log-likelihood function in (6) w.r.t. the parameter vector φ in (5) appears to be challenging, we adopt the EM algorithm that estimates the DOA parameters in an iterative manner: first, in the E-step (expectation step) of the k-th iteration, the lower-bound function Q(φ|φ (k) ) of the loglikelihood log p(X|φ) at the current iteration point φ (k) is constructed as

Q(φ|φ (k) ) = E Y |X,φ (k) {log p(Y |φ)} .
(7) In [START_REF] Ollila | Complex Elliptically Symmetric Distributions: Survey, New Results and Applications[END_REF], Y denotes the so-called (unobserved) complete data matrix. Then, in the M-step (maximization step) in the k-th iteration, the parameter vector φ is updated by maximizing the lower-bound function as follows:

φ (k+1) = arg max φ Q(φ|φ (k) ). (8)
The two steps in [START_REF] Ollila | Complex Elliptically Symmetric Distributions: Survey, New Results and Applications[END_REF] and ( 8) are alternatively performed until convergence. Details on the E-steps and M-steps are provided in the following subsections.

Expectation

Step By choosing the complete data as Y = [y(1) . . . , y(T )] with y(t) = x(t) T , j(t) T , τ (t) T , the lower-bound function 7) is rewritten as:

Q(φ|φ (k) ) in (
Q(φ|φ (k) ) = T t=1 Q t (φ|φ (k) ) (9) with Q t (φ|φ (k) ) = E j(t),τ (t)|x(t),φ (k) {log p( τ (t)|φ)} (10) + E j(t),τ (t)|x(t),φ (k) {log p( j(t)|τ (t), φ)} + E j(t),τ (t)|x(t),φ (k) {log p( x(t)|j(t), τ (t), φ)} .
For notational simplicity, we drop the time instant t and the DOA parameters θ in the remainder of this subsection. From ( 1) and ( 2), the expression in ( 10) is computed as

Q t (φ|φ (k) ) = -(M + R) log π + ν log ν -log Γ (ν) -M log σ 2 + (M + R + ν -1)log τ -||r|| 2 /σ 2 + ν τ + 2 σ 2 Re r H Bτ j - 1 σ 2 tr B H Bτ jj H -tr τ jj H , (11)
where (•) = E j,τ |x,φ (k) {•} denotes the expectation of the argument conditioned to the observed data x and the previously estimated parameters φ (k) . Using Bayes theorem and after some algebraic manipulation, the conditional expectation quantities in [START_REF] Rangaswamy | Robust Adaptive Signal Processing Methods for Heterogeneous Radar Clutter Scenarios[END_REF] are computed as follows

τ = α (k) β (k) , log τ = ψ α (k) -log β (k) , (12)
τ j = τ j, τ jj H = τ j j H + σ (k) 2 F (k) , (13)
with

r (k) = x-A (k) s (k) , Σ (k) = B (k) B (k) H + σ (k) 2 I, F (k) = B (k) H B (k) + σ (k) 2 I -1 , α (k) = M + ν (k) , β (k) = ν (k) + r (k) H Σ (k) -1 r (k) , j = F (k) B (k) H r (k)
and ψ(•) denotes the digamma function.

Maximization Step

In the M-step, the parameters θ, S, B, ν and σ 2 are updated by maximizing the expression in [START_REF] Rangaswamy | Robust Adaptive Signal Processing Methods for Heterogeneous Radar Clutter Scenarios[END_REF]. Neglecting constant terms, the function Q(φ|φ (k) ) in ( 9) is rewritten as

Q φ|φ (k) ∝ -M T log σ 2 - f A, S, B|φ (k) σ 2 -g(ν|φ (k) ) (14)
where two auxiliary functions g = g(ν|φ (k) ) and f = f A, S, B|φ (k) are defined by:

g = ν T t=1 τ (t) -log τ (t) + T log Γ (ν) -νT log ν, (15) f = ||(X -AS) D|| 2 F + tr B H BV -2Re tr BU (X -AS) H , (16)
respectively. Regarding the definition of the function f in ( 16), the auxiliary matrices

D = diag τ (1), . . . , τ (T ) , (17)
U = τ (1)j(1), . . . , τ (T )j(T ) ,

V = T t=1 τ (t)j(t)j(t) H , (18)
are constant. We remark from ( 14) that the number of DOF ν are decoupled from the remaining parameters. In the following, we provide the update of ν from the minimization of g = g(ν|φ (k) ) in [START_REF] Williams | Robust Estimation of Structured Covariance Matrices[END_REF].

Updating the number of degrees of freedom ν: A necessary condition for the updated number of DOF ν (k+1) is that the derivative g (ν|φ (k) ) = h(ν) in ( 15) vanishes at ν = ν (k+1) . The function h(ν) for ν > 0 is given by

h(ν) = T t=1 τ (t) -log τ (t) + T ψ(ν) -T log ν -T (20)
and its first derivative of h (ν) is given by h

(ν) = T ψ 1 (ν) -T /ν, (21)
where ψ 1 (•) denotes the trigamma function. By noting that τ -log τ ≥ 1 for all τ > 0, we obtain

T t=1 τ (t) -log τ (t) > T. (22)
From ( 22) and [21, Sec. 2], h(ν) is monotonically increasing from negative to positive. Therefore, the root of h(ν), which corresponds to the minimum of g(ν|φ (k) ), is unique [START_REF] Liu | ML Estimation of the t-Distribution using EM and its extensions, ECM and ECME[END_REF]. In order to obtain the minimum numerically, we adopt the rational approximation method [START_REF] Bunch | Rank-One Modification of the Symmetric Eigenproblem[END_REF]: At a given iteration point ν q , the function h(ν) is approximated by a rational function of type u(ν) = a -b ν with two parameters a and b such that h(ν q ) = u(ν q ) and h (ν q ) = u (ν q ).

(23) Then, the next iteration point ν q+1 is computed as the root of the rational approximant, i.e., u(ν q+1 ) = 0. The expression to compute the next iteration point ν q+1 from the current iteration point ν q is given in Step 4 of Algorithm 1. Using arguments similar to the proof in [START_REF] Bunch | Rank-One Modification of the Symmetric Eigenproblem[END_REF]Sec. 3], the convergence rate of Algorithm 1 is quadratic, which is identical to the Newton's method. In practice, the number of iterations of Algorithm 1 is approximately one third of that using the Newton's method for a given accuracy. This phenomenon can be explained as follows: Newton's method, which is based on linear approximation of the derivative, requires several restarts if the constraint ν q+1 > 0 is violated. In contrast, if the current point ν q is positive, it can be proven that the pole ν = 0 of the rational approximant u(ν) prevents the next iteration point ν q+1 from being negative. As a result, restart is not required in Algorithm 1 for any starting point ν 0 > 0. Updating the remaining parameters {θ, S, B, σ 2 } Taking the derivative of the cost function in ( 14) w.r.t. σ 2 while keeping A = A(θ), B and S fixed, the optimal value for σ 2 is given by

σ (k+1) 2 = f A, S, B|φ (k) M T , ( 24
)
Algorithm 1 Updating the number of DOF ν 1: Initialization: Inner iteration index q = 0, starting point ν q = ν (k) > 0 from the k-th EM iteration. 2: while not converged do 3:

Compute h(ν q ) in (20) and h (ν q ) in ( 21).

4:

Compute ν q+1 = ν 2 q h (νq) h(νq)+νqh (νq) .

5:

q ← q + 1. 6: end while 7: return Number of DOF ν (k+1) = ν q .

where f A, S, B|φ (k) is defined in [START_REF] Couillet | Robust Spiked Random Matrices and a Robust G-MUSIC Estimator[END_REF]. Substituting [START_REF] Mennad | Slepian-Bangs-Type Formulas and the Related Misspecified Cramér-Rao Bounds for Complex Elliptically Symmetric distributions[END_REF] back to [START_REF] Huber | Robust Statistics, Wiley Series in Probability and Statistics[END_REF], the maximization step in the EM algorithm is equivalent to the minimization of the function f = f A, S, B|φ (k) . As the function f is quadratic w.r.t. B, the minimizer B (k+1) is given by

B (k+1) = (X -AS) U H V -1 .
(25) Substitute B = B (k+1) in (25) back to ( 16), the concentrated cost function f w.r.t. A and S reads

f = ||(X -AS) D|| 2 F -(X -AS) U H V -1/2 2 F . (26)
From ( 13) and ( 16), it is not difficult to show that 27) and direct algebraic manipulation, the expression in ( 26) is reformulated as follows:

W = DD H -U H V -1 U DΠ ⊥ D -1 U H D H 0. (27) From (
f = (X -AS) W 1/2 2 F . (28)
The minimizer S (k+1) of (28) in dependence of A is given by S (k+1) = A † X.

(29) Substitute S = S (k+1) back to (28), the DOA vector θ (k+1) is obtained from the following optimization problem θ (k+1) = arg min

θ tr Π ⊥ A(θ) XW X H . (30)
We remark that the DOA estimation step in (30) resembles the conventional DML estimator under the spatially uniform Gaussian noise [START_REF] Böhme | Estimation of Source Parameters by Maximum Likelihood and Nonlinear Regression[END_REF]. Nevertheless, a key difference is the introduction of the weighting matrix W in (27), which accounts for the expected texture realization in each iteration. As the DOA estimation problem in (30) requires an N -dimensional search, we propose the use of the PR-DML estimator in [START_REF] Trinh-Hoang | Partial Relaxation Approach: An Eigenvalue-Based DOA Estimator Framework[END_REF], which is a relaxed variant of the conventional DML estimator. Specifically, the DOA parameters are estimated from the location of the N -deepest local minima of the PR-DML onedimensional null-spectrum:

F PR-DML (θ) = M m=N λ m (Π ⊥ a(θ) XW X H ). (31)
In (31), λ m (•) denotes the m-th largest eigenvalue of the matrix argument. We remark that the null-spectrum in (31) can be efficiently implemented as in [1, Sec. V.] without the full eigenvalue decomposition. Furthermore, the estimated DOAs from (31) results in a decrease of the Q function in (9) in practice, which implies the convergence of the EM Algorithm. To summarize, the proposed EM algorithm for the DOA estimation under non-Gaussian, spatially-colored interference-plusnoise is presented in Algorithm 2.

Algorithm 2 EM-PR-DML Algorithm 1: Initialization: Outer iteration index k = 0, initial quantities θ (0) , S (0) , B (0) , σ (0) 2 and ν (0) 2: while not converged do 3:

k ← k + 1.

4:

E-step: Compute necessary expectation quantities in ( 12)-( 13) and update the matrices D, U , V and W in ( 17)-( 19) and (27), respectively.

5:

M-step:

6:
Update ν (k+1) as presented in Algorithm 1.

7:

Update θ (k+1) from the local minima of the nullspectrum in (31) and

A (k+1) = A θ (k+1) . 8: Update S (k+1) = A (k+1) † Y . 9: Update B (k+1) = Y -A (k+1) S (k+1) U H V -1 . 10: Update σ (k+1) 2 = tr Π ⊥ A (k+1) XW X H M T .
11: end while 12: return The estimated DOAs θ (k)

SIMULATION RESULTS

In this section, simulation results regarding the estimation performance of different DOA estimators are presented and compared with the CRB in [START_REF] Böhme | Estimation of Source Parameters by Maximum Likelihood and Nonlinear Regression[END_REF]. The number of Monte-Carlo trials is N R = 5000. The Root-Mean-Squared-Error (RMSE) is computed as To illustrate the robustness of the proposed EM-PR-DML algorithm, in the first setup, we assume two uncorrelated source signals with unit power located at θ = [10 T .

N R r=1 N n=1 θ(r) n -θ n 2 /(N R N )
The SNR is varied between -10 and 10 dB. The receive signal at the sensor array is additionally corrupted by R = 2 interfering signals with ν = 2. However, we assume an incorrect number of interferers R = 4 in all considered algorithms. In Fig. 1, we observe that the AML estimator [START_REF] Besson | Directionof-Arrival Estimation in a Mixture of K-distributed and Gaussian Noise[END_REF] cannot detect two source signals due to the corruption of spatially-correlated low-rank interference. The estimation performances of the conventional MUSIC [START_REF] Schmidt | Multiple Emitter Location and Signal Parameter Estimation[END_REF] and PR-DML [START_REF] Trinh-Hoang | Partial Relaxation Approach: An Eigenvalue-Based DOA Estimator Framework[END_REF] slightly degrades due to the non-Gaussianity of the interference and noise. The RG-MUSIC [START_REF] Couillet | Robust Spiked Random Matrices and a Robust G-MUSIC Estimator[END_REF] and our proposed EM-PR-DML algorithm exhibit similar performances. Nevertheless, only our proposed estimator achieves the CRB in the high SNR regime. In Fig. 2, we assume no interference corrupting the receive signal. The sources are closely-spaced at θ = [10 • , 15 • ] T .

However, we assume the presence of a Gaussian noise We assume the number of interferers in our proposed EM-PR-DML estimator to be R = 2 to account for the model mismatch [START_REF] Mennad | Slepian-Bangs-Type Formulas and the Related Misspecified Cramér-Rao Bounds for Complex Elliptically Symmetric distributions[END_REF] of the noise when compared with that under the signal model in Section 2. We observe that the non-robust MUSIC and PR-DML estimator cannot resolve the DOAs of two sources in all Monte-Carlo trials. The proposed EM-PR-DML estimator outperforms RG-MUSIC and AML methods in both threshold and asymptotic regimes.

CONCLUSION AND OUTLOOK

In this paper, we propose a novel robust DOA estimator based on the EM algorithm and the PR framework. The proposed estimator jointly estimates the DOAs and the non-Gaussianity of the receive signal in each EM iteration. In order to maintain a low computational complexity, the PR framework is applied to substitute for the conventional deterministic ML. Simulation results show that the proposed DOA estimator possesses near-optimal performances under the assumed model without sacrificing robustness against model mismatch.

For future work, the theoretical explanation of the robustness of the EM-PR-DML estimator requires further comprehensive investigation. In addition, extensions to consider more realistic scenarios, e.g., decentralized DOA estimation among local subarrays or missing data, are of great interest.

  This work was supported by the DFG PRIDE Project PE 2080/2-1.

2 F

 2 N ] T and the true DOAs θ = [θ 1 , . . . , θ N ] T are sorted in ascending order. In our simulations we assume a ULA of M = 10 antennas with the spacing equal to half of the wavelength. The number of snapshots is set to T = 100. The matrix B associated with the interfering signals is chosen such that ||B|| = M R. The Signal-to-Noise-Ratio (SNR) is calculated as SNR = 1/σ 2 .

Fig. 1 :Fig. 2 :

 12 Fig. 1: RMSE vs SNR, overestimate the interference rank

Gamma (κ 1 , κ

) denotes the Gamma distribution with the shape parameter κ 1 and the inverse-scale parameter κ 2


P Landois

S Nanot

B Jouault

D Cakiroglu

S Contreras

S Juillaguet

H Peyre

J-M Decams

O Drachenko

J Leotin

C Coillot

GRAPHENE: FROM MATERIAL TO APPLICATIONS

come

Magnetotransport mode Graphene Hall cross design by electron beam lithography

-Size = 600µmx80µm -Graphene etching -Contact pads Ti/Au (10nm/100nm) deposition -Wire bonding -Check Raman

Magnetometer mode Doping-meter mode Conclusion & perspectives

The spinning current method combined to current modulation allows to remove LF noise giving rise to equivalent magnetic noise at low frequency (i.e. <1Hz) significantly impoved with respect to previous work [7]. Improving factors will permit to make graphene suitable candidate for magnetic field sensing.

The doping-meter mode is a first step towards the gas sensing ability [8] which will need to fonctionalize the graphene.

Acknowledgments

This project is supported by the LabEx NUMEV incorporated into the I-Site MUSE. Special thanks to L2C-Team « Structures Carbonées de Basses Dimensionnalités ». [START_REF]Special Report on Climate Change and Land[END_REF]. In France, agriculture is responsible for 19 % of the national GHG emissions caused by 3 gases : N 2 O (50%), CH 4 (40%) et CO 2 (10%). Agriculture is doubly concerned since it already faces dramatic consequences, due to climate change, like decrease of crop yields. That motivates the transition toward a new agriculture paradigm, known as agroecology. Agroforestry is one of the agroecologic practice, which combines trees and cereal crops. Growing trees and roots will sequester carbon in soils allowing to mitigate the carbon impact of the cultivated parcell [START_REF] Seghieri | Agroforesterie et services écosystémiques en zone tropicale[END_REF]. These practices require careful studies of the GHG of the cultivated parcel which would be sustained by new gaz sensor technologies : cheap, robust and sensitive.

n=4×10 12 cm -2 µ=1110 cm 2 / (V . s)

From graphene to sensor : development of a dedicated electronic board (@300K)

Graphene which has been early identified as very promising gaz sensor [START_REF] Schedin | [END_REF] still focuses interest in this direction [4,5].

Is graphene sensitive enough to detect greenhouse gaz at a level allowing to sustain agroecologic studies ?

⇒V Hall = IB cos (ω B t

) nq

Biasing current modulation/demodulation + Spinning current [6] => cancellation of the offset

⇒ n(t )= IBcos( ω B t ) V Hall ( t) q
Magnetic field modulation + spinning current + current modulation (option)

→ temporal tracking of doping ⇒ σ=qn μ

I =10μ A-Δ f =25 mHz
B detectable ≈25 µ T / √ Hz @ 1 Hz

  a) Conductivity measurement -Constant current biasing → 2 wires V(I) (Hall cross) or 4 wires Van der Pauw (Square sample) b) Doping measurement@300K -Magnetic field modulation -Constant current biasing → Hall voltage measurement 00by 1: graphene on SiC -2: magnetic coil -3: analog board -4: digital board Context : agriculture is one of the major cause of the climate change

  Icos(ω I t) B nq+ R offset Icos( ω I t )


Emmanuel Moulay

Wilfrid Perruquetti

Lyapunov-based approach for finite time stability and stabilization

Finite time stability is investigated for continuous system ẋ = f (x) which satisfies uniqueness of solutions in forward time. A necessary and sufficient condition is given for this class of systems using Lyapunov functions. Then, a necessary and sufficient condition is developed for finite time stabilization of class CL k -affine systems involving a class CL 0 -settling time for the closed-loop system. Finally an explicit feedback control is given using a control Lyapunov function verifying a certain inequality.

I. INTRODUCTION

This paper deals with the finite time stability of systems with the uniqueness of solutions in forward time and the finite time stabilization problem of affine systems. The aim is to provide a necessary and sufficient condition for finite time stability and finite time stabilization. Lyapunovlike techniques have been successfully used to solve these problems. Finite time stability and stabilization have often been a subject of research. Thus, Haimo gives a sufficient condition for finite time stability of continuous systems

ẋ = f (x), x ∈ R n (1)
in [START_REF] Haimo | Finite time controllers[END_REF]. Bhat and Bernstein provided an important contribution in [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF] by proving that there is a necessary and sufficient condition for finite time stability involving the continuity of the settling-time function at the origin. A part of our results is based on [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF] (the reader can find some additional and useful results in this paper). In general (without the continuity of the settling-time function at the origin) such necessary and sufficient conditions have not been discussed in the literature. The paper is organized as follows. After some notations in section II, section III states a general necessary and sufficient condition for the finite time stability using the properties of the settling-time function. Then a necessary and sufficient condition of finite time stabilization using the control Lyapunov function is addressed in section IV. Based on this result, an explicit feedback control is given under some less restrictive sufficient conditions.

II. NOTATIONS

Let us introduce some notations and definitions that will be useful later. W. Perruquetti is with the LAGIS (UMR CNRS 8021), Ecole Centrale de Lille, (e-mail: wilfrid.perruquetti@ec-lille.fr) Notation 1: Let E and F be two vector spaces and k ≥ 1, we denote by CL 0 (E, F) (respectively CL k (E, F)) the set of continuous functions on E, locally Lipschitz on E \ {0} with value in F (respectively the set of continuous functions on E, C k on E \ {0} with value in F).

Throughout this paper, V will be a non empty neighborhood of the origin in 1) possesses unique solutions in forward time on U ⊂ R n if for all x 0 ∈ U and two right maximally defined solutions of (1):

R n , B n the open unit ball in R n . As usually, a function V : V → R is proper if for every compact set K ⊂ R, V -1 (K) is compact. System (
φ x 0 : [0, T φ [ → R n and ψ x 0 : [0, T ψ [ → R n , there exists 0 < T x 0 ≤ min {T φ , T ψ } such that φ x 0 (t) = ψ x 0 (t) for all t ∈ [0, T x 0 [
. We may assume that for each x 0 ∈ U , T x 0 is chosen to be the largest in R ≥0 . In the following, φ x0 (t) denotes a solution of system (1) starting from x 0 ∈ R n at t = 0. Various sufficient conditions for forward uniqueness can be found in [START_REF] Agarwal | Uniqueness and Nonuniqueness Criteria for Ordinary Differential Equations[END_REF], [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF]Chapter 10] or [START_REF] Kawski | Stabilization of nonlinear systems in the plane[END_REF].

Now, let us recall some concepts of non-smooth analysis.

Let [a, b] ⊂ R, the upper Dini derivative of a function f : [a, b] ⊂ R → R is the function D + f : [a, b] → R defined by: D + f (x) = lim sup h→0 + f (x + h) -f (x) h . If V : V → R is a continuous function, V is the upper- right Dini derivative of V along the solutions of (1), that is V (x) = D + (V • φ x ) (0). If V is locally Lipschitz at x ∈ V, then V (x) = lim sup h→0 + V (x + hf (x)) -V (x) h and if V is continuously differentiable at x ∈ V, then V (x) = ∇V (x), f (x) .
As it is customary in control theory, a Lyapunov function V for the system (1) is a continuous positive definite function such that V is negative definite. The Lie derivative of V : R n → R along f : R n → R n is defined by:

L f V : R n → R, L f V (x) = ∇V (x), f (x) .
In section IV, we need some concept on set-valued functions.

A set-valued function Φ from X to Y is a function that maps x ∈ X to a set Φ (x) ⊂ Y. Let X and Y be two vector spaces and Φ :

X → Y a set-valued function, Φ is lower semi-continuous if {x ∈ X : Φ (x) ∩ O = ∅} is open in X for every open set O ⊂ Y.
Φ is locally Lipschitz if for any x 0 ∈ X , there exists a neighborhood N (x 0 ) ⊂ X and a constant l ≥ 0 such that for all x, x ∈ N (x 0 ),

Φ(x) ⊂ Φ(x ) + l x -x X B Y
where B Y is the unit ball in Y.

III. FINITE TIME STABILITY OF SYSTEMS WITH THE UNIQUENESS OF SOLUTIONS IN FORWARD TIME

Let us consider the system (1) where f ∈ C 0 (R n ), f (0) = 0 and when f possesses the property of uniqueness of solutions in forward time. Let us recall the asymptotic stability which is a local property. The origin is asymptotically stable for the system (1) if:

1) the origin is stable for the system (1), 2) the origin of the system (1) is attractive: for all > 0, there exists δ ( ) > 0 such that each solution starting from x 0 ∈ δ ( ) B n tends to the origin as t tends to infinity. Now, we may recall the notion of finite time stability involving the settling-time function (see [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF]Definition 2.2])

Definition 1: The origin is finite time stable for the system (1) if there exists a non empty neighborhood of the origin V in R n such that:

1) there exists a function T :

V \ {0} → R ≥0 such that if x 0 ∈ V \ {0} then φ x 0 (t) is defined (and particularly unique) on [0, T (x 0 )[, φ x 0 (t) ∈ V \ {0} for all t ∈ [0, T (x 0 )[ and lim t→T (x0) φ x 0 (t) = 0.
T is called the settling-time of the system (1). 2) for all > 0, there exists δ ( ) > 0, for every

x 0 ∈ (δ ( ) B n \ {0}) ∩ V, φ x 0 (t) ∈ B n for all t ∈ [0, T (x 0 )[. Remark 1:
First, note that if the origin of system (1) is finite time stable, then f cannot possess uniqueness in backward time at the origin, in particular f cannot be locally Lipschitz at the origin. Then, if system (1) is finite time stable, Lyapunov asymptotic stability implies that φ 0 ≡ 0 is the unique solution starting from x 0 = 0. So, the settling-time function T may be extended at the origin by T (0) = 0. We will also call this extension the settling-time of the system [START_REF] Agarwal | Uniqueness and Nonuniqueness Criteria for Ordinary Differential Equations[END_REF].

The following result is given in [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF]Proposition 2.3]. Lemma 1: Suppose that the origin is finite time stable for the system (1) with the settling-time function T : V → R ≥0 , then for all x ∈ V the flow Φ(t, x) = φ x (t) of the system (1) is defined and continuous on R ≥0 × V and Φ(t, x) = 0 for all t ≥ T (x).

This result shows that the finite time stability of system (1) implies:

• the uniqueness in forward time of solutions starting from V, • the asymptotic stability,

• the existence of a continuous flow. Finally, T (x) is the time for the solution φ x to reach the origin, and as the system is autonomous, the Lyapunov stability ensures that the solution stays at the origin for any time longer than T (x). Then, the equality is as follows

T (x) = inf {t ∈ R ≥0 : Φ(t, x) = 0} . ( 2
)
Let us recall the fundamental theorem of Kurzweil which is in [START_REF] Kurzweil | On the inversion of Liapunov's second theorem on stability of motion[END_REF]Theorem 7].

Theorem 1 (of Kurzweil): Let us consider the system (1) such that f is continuous, the system (1) is asymptotically stable if and only if there exists a smooth Lyapunov function for the system (1).

Here, the theorem 1 of Kurzweil is of importance because we cannot have the Lipschitz continuity at the origin of the right-hand side of the system (1). Let us give the main result of this section which is a general necessary and sufficient condition for finite time stability.

Theorem 2: Let us consider the system (1) with uniqueness of solutions in forward time outside the origin, the following properties are equivalent:

(i) the origin of the system ( 1) is finite time stable on V, (ii) there exists a smooth Lyapunov function V : V → R ≥0 for the system (1) satisfying for all

x ∈ V 0 V (x) ds V (Φ (θ x (s) , x)) < +∞ (3)
where

θ x is the inverse of t → V (Φ (t, x)), (iii) all smooth Lyapunov functions V : V → R ≥0 for the system (1) satisfy for all x ∈ V 0 V (x) ds V (Φ (θ x (s) , x)) < +∞.
Moreover, if V is such a Lyapunov function then

T (x) = 0 V (x) ds V (Φ (θ x (s) , x))
.

Proof: (i) ⇒ (ii) If the system (1) is finite time stable with the settling-time function T : V → R ≥0 , then, there exists a smooth Lyapunov function V : V → R n for the system (1) given by the theorem 1 of Kurzweil. So, the well defined application [0,

T (x)[ → ]0, V (x)], t → V (Φ (t, x)) is strictly decreasing and differentiable, so its inverse ]0, V (x)] → [0, T (x)[, s → θ x (s) is differentiable and satisfies for all s ∈ ]0, V (x)], θ x (s) = 1 V (Φ (θ x (s) , x))
.

The use of the change of variables s = V (Φ (t, x)) leads to the following equalities

T (x) = T (x) 0 dt = 0 V (x) θ x (s) ds = V (x) 0 ds -V (Φ (θ x (s) , x)) < +∞. (4)
(ii) ⇒ (iii) The set of smooth Lyapunov function for the system (1) S L is non empty. Let V ∈ S L , then by using the same argument as before, we deduce that V satisfies (3). (iii) ⇒ (i) As there exists a Lyapunov function for the system (1), the theorem of Lyapunov (see [START_REF] Hahn | Theory and Application of Liapunov's Direct Method[END_REF]) ensures that it is asymptotically stable. The equalities (4) imply the finite time convergence.

Remark 2: If the Lyapunov function V is defined on R n , proper, and if the condition (ii) is globally held, then, the origin of the system (1) is globally finite time stable.

Even if theorem 2 is a theoretical result, we may give a simple example.

Example 1: Let us consider the Cauchy problem

ẋ = -x x (0) = 1
and the smooth Lyapunov function

V (x) = x 2 2 . Then, V (φ (t)) = e -2t , θ (s) = -1 2 ln (s) and φ (θ (s)) = √ s lead to V (φ (θ (s)
)) = -s where s > 0. We have

T (1) = 1 2 0 ds s = +∞.
Theorem (2) ensures that the system ẋ = -x is not finite time stable.

In general, the settling-time function is not continuous at the origin. A fundamental example is given in [4, Example 2.2] which shows that the settling-time function of a finite time stable system is generally non continuous at the origin. Now, we may give a corollary to theorem 2 which gives a sufficient condition for a continuous (or class CL 0 ) settlingtime function.

Corollary 1: Let us consider the system (1) with the uniqueness of solutions in forward time outside the origin. Let us assume that there exists a smooth Lyapunov function V : V → R ≥0 for the system (1) and g ∈ L 1 ([0, sup x∈V V (x)]) such that for all x ∈ V \ {0}, and all s ∈ [0, V (x)] -1 V (Φ (θ x (s) , x)) ≤ g (s)

then the system (1) is finite time stable with a continuous settling-time function.

Proof: If there exists a smooth Lyapunov function V : V → R ≥0 for the system and a function

g ∈ L 1 ([0, sup x∈V V (x)]) such that for all x ∈ V \ {0}, and all t ∈ [0, V (x)] , -1 V (Φ(θx(s),x)) ≤ g (s) then V (x) 0 -ds V (Φ(θ x (s),x)) ≤ V (x) 0 g (s) ds < +∞ for all x ∈ V \ {0}.
As T (0) = 0, we may deduce that the system (1) is finite time stable. Moreover, lim

x →0 T (x) ≤ lim x →0 V (x) 0
g (t) dt = 0, so the settling-time function is continuous at the origin. To conclude, we may invoke the following result [4, Proposition 2.4.] which shows that T is continuous at the origin if and only if T is continuous on its domain of definition V.

Theorem 2 is quite general. Nevertheless, its application is not easy because the flow is generally unknown. In order to study the stabilization problem, we prefer to restrict the problem to the case of a continuous settling-time function at the origin. We could refer to a result given in [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF] in order to use a necessary and sufficient condition involving a Lyapunov function only. We may be inclined to use a more regular settling-time function in order to find a more regular Lyapunov function. This is important for the problem of stabilization in section IV when using the control Lyapunov functions. For this, we might need the following lemma which is similar to a result given in [5, Lemma 5.1], and [START_REF] Teel | A smooth Lyapunov function from a class-KL estimate involving two positive semidefinite functions[END_REF]Lemma 16].

Lemma 2: Let O be a non empty open subset of R n , and

β : O → R, µ, ν : O → R >0 three continuous functions. Suppose V : O → R is locally Lipschitz on O, if ∀x ∈ O, ∇V (x), f (x) ≤ β(x)
then there exists a smooth function V : O → R such that, for all x ∈ O,

V (x) -V (x) ≤ µ(x) ∇ V (x), f (x) ≤ β(x) + ν(x).
Let us recall a result which can be found in [4, Theorem 4.2] for the sufficient condition of finite time stability and [4, Theorem 4.3] for the necessary one.

Proposition 1: Consider the system (1) with the uniqueness of solutions in forward time outside the origin, the following properties are equivalent:

1) the origin of the system (1) is finite time stable with a continuous settling-time function at the origin, 2) there exists a real number c > 0, α ∈ ]0, 1[ and a

Lyapunov function

V : V → R ≥0 satisfying V (x) ≤ -c(V (x)) α (5)
for all x ∈ V. The construction of the Lyapunov function, in the proof given in [4, Theorem 4.3], involved the settling-time function in the following sense: V (x) = T (x)

1 1-α with α ∈ ]0, 1[ ( 1 1-α > 1)
. As a Lyapunov function is at least continuous, it involves the continuity of the settling-time function at the origin which is equivalent to the continuity of the settlingtime function on its domain of definition (see [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF]Proposition 2.4.]). Now, we may give a variant of this result dedicated to the class CL k -systems by using the fact that the class of systems with the uniqueness of solutions in forward time is included in the class CL k -systems for all k ≥ 0.

Proposition 2: Let k ≥ 0, if f belongs to the class CL k then the following properties are equivalent:

1) the origin of the system (1) is finite time stable with a class CL 0 settling-time function, 2) there exists a real number c > 0, α ∈ ]0, 1[ and a class CL ∞ -Lyapunov function V : V → R ≥0 satisfying the condition (5). Moreover, if V is a Lyapunov function satisfying the condition (5) then for all x ∈ V,

T (x) ≤ V (x) 1-α c(1 -α) .
Proof: Suppose that 2) is verified. We can find in [4, Theorem 4.2] that the existence of the class CL ∞ -Lyapunov function V satisfying condition (5) implies the finite time stability. Moreover, as x →

V (x) 1-α c(1-α) is a class CL ∞ -function,
we deduce that T belongs to the class CL 0 . Conversely, let us suppose that the origin of the system (1) is finite time stable with a class CL 0 -settling-time function. By using the proof of proposition 1 which can be found in [4, Theorem 4.3], we know that V (x) = T (x) 1 1-α with 0 < α < 1 belongs to the class CL 0 (V) and is a Lyapunov function for the system (1) satisfying condition [START_REF] Clarke | Asymptotic stability and smooth Lyapunov function[END_REF]. Let 0 < k < c, we apply the lemma 2 with the open set V \ {0}, and

β(x) = -cV (x) α , µ(x) = 1 2 V (x), ν(x) = kV (x) α to obtain a class CL ∞ -Lyapunov function V such that, 1 2 V (x) ≤ V (x) ≤ 3 2 V (x), x ∈ V, ∇ V (x), f (x) ≤ -c V (x) α , x ∈ V \ {0} , with c = c -k > 0.
As it is shown in the next two examples, the Lyapunov function satisfying condition (5) may be smooth everywhere.

Example 2 (scalar system): Let α ∈ ]0, 1[ and k > 0, it is easy to see that the basic system

ẋ = -k |x| α sgn(x), x ∈ R
is finite time stable using the smooth Lyapunov function

V (x) = x 2 2 with the well known class CL ∞ -settling-time function T (x) = |x| 1-α k(1-α) . Indeed, we have for all x ∈ R V (x) = -k |x| 1-α = -2 1+α 2 k V (x) 1+α 2
with 1+α 2 ∈ ]0, 1[. For this basic example, the solutions are explicit

φ x 0 (t) = g(t, x 0 ) if 0 ≤ t ≤ |x0| 1-α k(1-α) , 0 if t > |x 0 | 1-α k(1-α) , with g(t, x 0 ) = |x 0 | 1-α -k(1 -α)t 1 1-α sgn(x 0
), so we do not need a Lyapunov function. As it is recalled in the introduction, the settling-time function is given by T (x) = 0 x dy f (y) for finite time stable scalar systems (see [START_REF] Haimo | Finite time controllers[END_REF] and [START_REF] Moulay | Finite time stability of non linear systems[END_REF] for a proof of this basic result).

Example 3 (two dimensional system): Let us consider the system:

ẋ1 = -|x 1 | α sgn(x 1 ) -x 3 1 + x 2 ẋ2 = -|x 2 | α sgn(x 2 ) -x 3 2 -x 1 . Taking V (x) = x 2 2 , we obtain V (x 1 , x 2 ) = - 2 i=1 (x 4 i + |x i | α+1 ) ≤ 0. V is a Lyapunov function for the sys- tem satisfying V (x 1 , x 2 ) ≤ -2 α+1 2 V (x 1 , x 2 ) α+1 2 . Indeed, 2 i=1 (x 4 i + |x i | α+1 ) ≥ x 2 1 + x 2 2 α+1 2 = x α+1
. Thus the origin is finite time stable with a continuous settling-time function verifying T (x) ≤ 2 x 1-α 1-α .

IV. FINITE TIME STABILIZATION OF THE CLASS CL k -AFFINE SYSTEMS

Let k ≥ 0, and consider the following affine system

ẋ = f 0 (x) + m i=1 f i (x)u i , x ∈ R n and u ∈ U (6)
where U is a non empty open set of R m containing the origin, f i ∈ CL k (R n , R) for all 0 ≤ i ≤ m and f 0 (0) = 0 and the closed-loop system

ẋ = f 0 (x) + m i=1 f i (x)u i (x), x ∈ R n . ( 7
)
Let us recall the definitions of the stabilization and the finite time stabilization. We will restrict our study to the case of a class CL 0 -settling-time function for the finite time stabilization. The control system ( 6) is stabilizable (respectively finite time stabilizable) if there exists a non empty neighborhood of the origin V in R n and a feedback control law u ∈ C 0 (V \ {0} , U) such that:

1) u (0) = 0, 2) the origin of the system ( 7) is asymptotically stable (respectively finite time stable with a class CL 0 -settlingtime function). Here, we give a necessary and sufficient condition for the finite time stabilization of the system (6) involving the continuity of the settling-time function at the origin for the closed-loop system [START_REF] Hahn | Theory and Application of Liapunov's Direct Method[END_REF]. We add a condition to the concept of control Lyapunov function first defined in [START_REF] Artstein | Stabilization with relaxed controls[END_REF], which leads to the finite time stabilization.

We are going to recall some usual definitions. A positive definite function V ∈ CL ∞ (V, R ≥0 ) is a control Lyapunov function for the system (6) if for all x ∈ V \ {0},

inf u∈U (a (x) + B (x) , u ) < 0. where a (x) = L f 0 V (x), B (x) = (b 1 (x) , ..., b m (x)) with b i (x) = L f i V (x) for 1 ≤ i ≤ m.
To obtain the finite time stabilization, we have to bring in the control Lyapunov function V the following condition which holds for all x ∈ V \ {0} and for a real number

α ∈ ]0, 1[ inf u∈U (a (x) + B (x) , u ) ≤ -c (V (x)) α . ( 8
)
As usual, such a control Lyapunov function satisfies the small control property if for each > 0, there exists δ > 0 such that, if x ∈ δB n , then there exists some u ∈ B m such that a (x) + B (x) , u < 0.

Remark 3: If U ⊂ R, the small control property is equivalent to lim sup

x n →0 a(x) |B (x)| ≤ 0.
The limit may very well be -∞. We set b (x) = B (x) 2 . The theorem 3 of Mickael will help us to show our main result on finite time stabilization.

Theorem 3 (of Mickael): Let X and Y be two metric spaces, for every lower semi-continuous (respectively locally Lipschitz) set-valued function Φ : X → 2 Y , x → Φ(x) where 2 Y will denote the family of non-empty, closed, convex subsets of Y it is possible to extract a continuous (respectively locally Lipschitz) function f such that f (x) ∈ Φ(x) for all x ∈ X .

Theorem 4: The system ( 6) is finite time stabilizable under a class CL 0 -feedback control if and only if there exists a control Lyapunov function for the system (6) which satisfies the condition (8) for a real number α ∈ ]0, 1[ and the small control property.

Proof: If the control system ( 6) is finite time stabilizable, then the closed-loop system ( 7) is finite time stable with a class CL 0 -settling-time function. By using proposition 2, there exists a class CL ∞ -Lyapunov function V for the closed-loop system [START_REF] Hahn | Theory and Application of Liapunov's Direct Method[END_REF] satisfying the condition (5) which implies that ( 8) is valid. Moreover, by using the feedback control u(x) and its continuity, it is easy to see that the control Lyapunov function satisfies the small control property. Conversely, if there exists a control Lyapunov function V : V → R ≥0 for the system (6) satisfying condition ( 8), then we introduce the set valued function Φ defined for x ∈ V\ {0} by

Φ(x) = {v ∈ U : a (x) + B (x) , v ≤ -c(V (x)) α } .
As v → a(x) + B (x) , v is affine, it implies that for all x ∈ V\ {0}, Φ(x) belongs to the family of non-empty closed convex subsets of U for the subspace topology. As f i belongs to the class CL k for all 0 ≤ i ≤ m and V ∈ CL ∞ (V, R ≥0 ), a (x) + B (x) , v + c(V (x)) α is locally Lipschitz for all x ∈ V \{0}. Thus, we may deduce that Φ is locally Lipschitz on V \ {0}. As V satisfies the small control property, it is shown in [START_REF] Artstein | Stabilization with relaxed controls[END_REF]Theorem 4.3] that we may extend Φ on V by Φ(0) = {0} such that Φ now is lower semi-continuous on V. We may apply the theorem 3 of Mickael to find a selection u ∈ CL 0 (V, U ). Then V is a class CL ∞ -Lyapunov function for the closed loop system (7) satisfying the condition 5. Thus, by using proposition 2 we deduce that the system ( 6) is finite time stabilizable.

Theorem 4 provides a tool for the finite time stabilization with a class CL 0 -settling-time function.

In practical terms, the resolution of the finite time stabilization is a delicate task which has generally been studied for homogeneous systems of negative degree with respect to a flow of a complete vector field. Indeed, for this kind of systems, finite time stability is equivalent to asymptotic stability (see [START_REF] Bhat | Continuous, bounded, finite-time stabilization of the translational and rotational double integrator[END_REF], [START_REF] Hong | Finite-time stabilization and stabilizability of a class of controllable systems[END_REF] for more details). Nevertheless, if we want to use a control Lyapunov function to obtain a constructive feedback control for finite time stabilization, we can use a modified version of the Sontag feedback control given in [START_REF] Sontag | A universal construction of Arststein's theorem on nonlinear stabilization[END_REF].

Lemma 3: If there exists a continuously differentiable control Lyapunov function V : V → R ≥0 for the control system (6), then it is stabilizable under the feedback control u (x) = (u 1 (x) , . . . , u m (x)) defined by

u i (x) = -b i (x) a(x)+ p √ a(x) p +b(x) q b(x) if x ∈ V\ {0} 0 if x = 0 (
9) where p, q ≥ 2 are even integers. If furthermore V satisfies the small control property, then the feedback control ( 9) is also continuous at the origin.

Proof: Suppose there exists a smooth control Lyapunov function

V : V → R ≥0 . Let E = (x, y) ∈ R 2 : x < 0 or y > 0 and ϕ a function defined on E by ϕ (x, y) =    x + p √ x p + y q y if y = 0 0 if y = 0 . As lim y→0 x + p √ x p + y q y = lim y→0 x + |x| p 1 + y q x p y = lim y→0 -y q-1 px p-1 = 0, ϕ is continuous on E.
As V is a control Lyapunov function, then we know that (a(x), b (x)) ∈ E for all x ∈ V\ {0}. Thus, we define the feedback control by

u i (x) = -b i (x) ϕ (a(x), b (x)). u (x) is continuous on V\ {0} and we obtain for all x ∈ V\ {0} ∇V (x) , f 0 (x) + m i=1 f i (x)u i (x) = -p a(x) p + b (x) q < 0.
So, V is a Lyapunov function for the closed-loop system [START_REF] Hahn | Theory and Application of Liapunov's Direct Method[END_REF], and by using the Lyapunov theorem we know that the origin of the closed loop system ( 7) is asymptotically stable. The proof concerning the stabilization under the small control property is similar to the one given in [START_REF] Sontag | A universal construction of Arststein's theorem on nonlinear stabilization[END_REF]Theorem 1].

Proposition 3: If there exists a continuously differentiable control Lyapunov function V : V → R ≥0 for the control system (6) verifying the small control property and for all

x ∈ V, p a(x) p + b (x) q ≥ cV (x) α
where p, q ≥ 2 are even integers, and where c > 0 and 0 < α < 1 then the system (6) is finite time stabilizable under the continuous feedback control [START_REF] Hong | Finite-time stabilization and stabilizability of a class of controllable systems[END_REF]. Proof: The asymptotic stability is proved using Lemma 3. It is shown in [13, Theorem 1] that if V satisfies the small control property, the feedback control ( 9) is continuous at the origin. The inequality

∇V (x), f 0 (x) + m i=1 f i (x)u i (x) = -p a(x) p + b (x) q ≤ -cV (x) α
ensures the finite time convergence by using proposition 1.

Remark 4: Proposition 3 is also true for an only continuous control system [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF]. The proof of the finite time convergence can be given by using a result of Haimo in [8, Proposition 1].

Example 4: Let 0 < β, γ < 1 such that 4 (β + γ) < β +1 (for example β = γ = V. CONCLUDING REMARKS The problem of finite time stability of systems with the uniqueness of solutions in forward time is solved for differential equations by giving a necessary and sufficient condition for the finite time stability involving a Lyapunov function. Our results bring an answer to the question asked by Bhat and Bernstein in the conclusion of their paper [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF] concerning a stronger converse result for finite time stability. Moreover, by using their results on finite time stability involving continuity of the settling-time function at the origin, we have succeeded to solve the problem of the finite time stabilization of class CL k -affine systems involving a class CL 0 -settling-time function. The universal controller given by Sontag in [START_REF] Sontag | A universal construction of Arststein's theorem on nonlinear stabilization[END_REF] is extended to design a feedback control for the finite time stabilization. Nevertheless, our paper raises certain questions that are important from the point of view of the stabilization theory, in particular the construction of a universal finite time feedback control using a control Lyapunov function satisfying condition [START_REF] Haimo | Finite time controllers[END_REF].

  This work was supported in part by programs Maths-STIC CNRS 19 and ARCIR ROBOCOOP E. Moulay is with Laboratoire Paul Painlevé, UFR de Mathématiques Pures et Appliquées, Université des Sciences et Technologies de Lille and LAGIS (UMR-CNRS 8021), Ecole centrale de Lille, Cité Scientifique, 59651 Villeneuve d'Ascq Cedex, France (e-mail: emmanuel.moulay@eclille.fr)

1 8 )< 1 .

 81 and consider the systemẋ1 = -|x 1 | γ sgn(x 1 ) -x 2 ẋ2 = |x 1 | β sgn(x 1 ) |x 2 | 1-β + |x 2 | γ u . Using the class C 1 -function V (x) = |x 1 | β+1 + |x 2 | β+1 , we obtain a (x) = -(β + 1) |x 1 | β+γ B (x) = (β + 1) |x 2 | β+γ sgn(x 2 ) b (x) = (β + 1) 2 |x 2 | 2(β+γ) . As inf u∈R (a (x) + B (x) u) < 0 for x = 0, V is a controlLyapunov function for the system. The fact thata (x) |B (x)| = -|x 1 | β+γ |x 2 | β+γ ≤ 0and remark 3 implies that V satisfies the small control property. Now, by using the feedback control (9) with (p, q) = (4, 2), we obtaina (x) 4 + b (x) 2 = (β + 1) 4 |x 1 | 4(β+γ) + |x 2 | Thus 4 a(x) 4 + b (x) 2 ≥ (β + 1) V (x)α with α = β+γ β+1 < 1, by using proposition (3) we know that the control system is finite time stabilizable under the continuous feedback controlu (x) = |x 1 | β+γ -4 |x 1 | 4(β+γ) + |x 2 | 4(β+γ)|x 2 | β+γ sgn(x 2 ) .