Emma Hoarau

Claire David

Pierre Sagaut

Thiên-Hiêp Lê

LIE GROUP STABILITY OF FINITE DIFFERENCE SCHEMES

Differential equations arising in fluid mechanics are usually derived from the intrinsic properties of mechanical systems, in the form of conservation laws, and bear symmetries, which are not generally preserved by a finite difference approximation, and leading to inaccurate numerical results. This paper develops a method that enables us to build a scheme that preserves those symmetries. The method is based on the concept of the differential approximation. A comparison of numerical performance of the invariant schemes, standard ones and higher order one has been realised for the Burgers equation.

1. Introduction. Lie groups were introduced by Sophus Lie in 1870 in order to study the symmetries of differential equations, yielding thus analytical solutions. Literature provides substantial works and applications, [START_REF] Olver | Applications of Lie Groups to Differential Equations[END_REF], [START_REF] Ibragimov | Symmetries exact solutions and conservation laws[END_REF]. Symmetry groups can be determined by an automatic procedure, but it often turn out to be tedious and induce errors. A large amount of packages using symbolic manipulations of mathematical expressions have been written. We mention here some of those works: Schwartz [START_REF] Schwarz | A REDUCE package for determining Lie symmetries of ordinary and partial differential equations[END_REF], Vu and Carminati [START_REF] Vu | Symbolic computation and differential equations: Lie symmetries[END_REF], Herod [START_REF] Herod | MathSym: a Mathematica program for computing Lie symmetries[END_REF], Baumann [START_REF] Baumann | Lie symmetries of differential equations: A mathematica program to determine Lie symmetries[END_REF], Cantwell [START_REF] Cantwell | Introduction to symmetry analysis[END_REF]. In this paper we are interesting in the application of the theory of Lie group to numerical analysis. Finite difference equations used to approximate the solutions of a differential equation generally do not respect the symmetries of the original equation, and can lead to inaccurate numerical results. Various techniques, that enable us to build a scheme preserving the symmetries of the original differential equation, have been studied. One of these techniques consists in constructing an invariant scheme from a given one by applying the method of the moving frame in [START_REF] Olver | Geometric foundations of numerical algorithms and symmetry[END_REF], [START_REF] Kim | Invariantization of numerical schemes using moving frames[END_REF]. Another one consists in constructing an invariant scheme with the help of the discret invariants of its symmetry group [START_REF] Dorodnitsyn | Finite difference models entirely inheriting continuous symmetry of original differential equations[END_REF], [START_REF] Dorodnitsyn | Lie group classification of second order difference equations[END_REF], [START_REF] Dorodnitsyn | A heat transfert with a source: the complete set of invariant difference schemes[END_REF], [START_REF] Bakirova | Symmetry-preserving difference schemes for some heat transfert equations[END_REF], [START_REF] Valiquette | Discretization of partial differential equations preserving their physical symmetries[END_REF] and provides the building of symmetry-adapted meshes, in preserving the differential equation symmetries. This technique is based on a direct study of the symmetries of difference equations and lattices. Yanenko [START_REF] Yanenko | Group classification of difference schemes for a system of one-dimensional equations of gas dynamics[END_REF] and Shokin [START_REF] Yu | The method of differential approximation[END_REF], proposed to apply the Lie group theory to finite difference equations by means of the differential approximation. Thus, they have set down conditions under which the differential representation of a finite difference scheme preserves the group of continuous symmetries of the original differential equation. They provide a dissipative scheme, which is called invariant scheme. The resulting scheme is independent of any change of the reference frame, and its differential representation is invariant under the symmetries of the original equation. Ames, Postell and Adams [START_REF] Ames | Optimal numerical algorithms[END_REF] have already used the approach of Yanenko and Shokin to present invariant schemes in which terms are added to the original 1 Ê difference scheme. They showed that, in specific cases, the invariant scheme is as accurate as high order numerical methods. In this paper, we focus on the last approach. The method is implemented on some standard schemes for the Burgers equation. A comparison is made between the numerical solutions of these schemes and the invariant scheme. The paper is organized as follows. Definitions and invariance condition for differential equations are provided in section 2. Section 3 recalls the approach of Yanenko and Shokin. Section 4 concentrates on classical schemes. In section 5, we present a method that enables us to build an invariant scheme with respect to an otherwise lost symmetry.

2. Definitions and invariance condition for differential equations. A rparameter Lie group G r of point transformations in the Euclidean space E(x, u) can be written under the form:

Gr = {x * i = φ i (x, u, a); u * j = ϕ j (x, u, a), i = 1, . . . , m; j = 1, . . . , n} (1)
Consider a system of l th -order differential equations:

F λ x, u, u (k 1 ) , u (k 1 ,k 2 ) , . . . , u (k 1 ...k l ) = 0, λ = 1, . . . , q (2)
Denote by u (k1...kp) the vector, the components of which are partial derivatives of order p, namely, u

(k1...kp) j = ∂ p uj ∂x k 1 .
..∂x kp j = 1, . . . , n and k 1 , . . . , k p ∈ {1, . . . , m}. Denote by x = (x 1 , . . . , x m ) the independent variables, u = (u 1 , . . . , u n ) the dependent variables, and (x k1 . . . x kp ) a set of elements of the independent variables. Equation ( 2) is a subset of the Euclidean space E x, u, u (k1) , . . . , u (k1...k l ) . In order to take into account the derivative terms involved in the differential equation, the action of the group G r of transformations in the space E x, u) needs to be extended to the space of the derivatives of the dependent variables. Denote by G (l) r a r-parameter Lie group of point transformation in the space E x, u, u (k1) , . . . , u (k1...k l ) of the independent variables, dependent variables and the derivative of the dependent variables with respect to the independent ones. The l th -prolongation operator of G r is:

L (l) α = ξ α i (x, u) ∂ ∂x i + η α j (x, u) ∂ ∂u j + σ α,(k 1 ) j ∂ ∂u j (k 1 ) + • • • + σ α,(k 1 ...k l ) j ∂ ∂u j (k 1 ...k l ) , (3)
i = 1, . . . , m; j = 1, . . . , n; α = 1, . . . , r.

ξ α i , η α j , σ α,(k1) j and σ α,(k1...ko) j are given by:

ξ α i = ∂φ i ∂aα a=0 , η α j = ∂ϕ j ∂aα a=0 , σ α,(k 1 ) j = Dη α j Dx k 1 - m i=1 ∂u j ∂x i Dξ α i Dx k 1 σ α,(k 1 ...ko) j = Dσ α,(k 1 ...k o-1 ) j Dx ko - m i=1 ∂ o u j ∂x i ∂x k 1 . . . ∂x k o-1 Dξ α i Dx ko , o = 2, . . . , l
where:

D Dx k = ∂ ∂x k + n j=1 ∂u j ∂x k ∂ ∂u j
The system of l th -order differential equations is invariant under the group G (l) r if and only if:

L (l) α F λ F λ =0 = 0, α = 1, . . . , r; λ = 1, . . . , q (4)
3. Lie group for the differential approximation. The finite difference scheme, which approximates the differential system (2), can be written as:

Λ λ (x, u, h, T u) = 0, λ = 1, . . . , q (5)
where h = (h 1 , h 2 , . . . , h m ) denotes the space step vector, and T = (T 1 , T 2 , . . . , T m ) the shift-operator along the axis of the independent variables, defined by:

T i [u](x 1 , x 2 , . . . , x i-1 , x i , x i+1 , . . . , xm) = u(x 1 , x 2 , . . . , x i-1 , x i + h i , x i+1 , . . . , xm). (6)
Definition 1. The differential equation:

P λ x, u, u (k 1 ) , . . . , u (k 1 ...k l ′ ) = F λ x, u, u (k 1 ) , . . . , u (k 1 ...k l ) + s β=1 m i=1 (h i ) l β R λ i (x, u, u (k 1 ) , . . . , u (k 1 ...k l ′ λ,i ) ), λ = 1, . . . , q; l ′ = max (λ,i) l ′ λ,i (7)
is called the s th -order differential approximation of the finite difference scheme (5). In the specific case s = 1, the above equation is called the first differential approximation.

Equation ( 7) is obtained from equation ( 5) by applying Taylor series expansion to the components of T u about the point x = (x 1 , . . . , x m ) and truncating the expansion to a given finite order. Denote by G ′ r a group of transformations in the space E(x, u, h):

G ′ r = {x * i = φ i (x, u, a); u * j = ϕ j (x, u, a); h * i = ψ i (x, u, h, a), i = 1, . . . , m; j = 1, . . . , n} (8)
by L α ′ the basis infinitesimal operator of G ′ r :

Lα ′ = ξ α i (x, u) ∂ ∂x i + η α j (x, u) ∂ ∂u j + ζ α i (x, u, h) ∂ ∂h i , α = 1, . . . , r (9)
where

ζ α i = ∂ψ i ∂aα a=0 , α = 1, . . . , r (10)
and by G (l ′ ) r a group of transformation in the space E(x, u, h, u (k1) , . . . , u (k1...k l ′ ) ). The l ′ th -prolongation operator of G ′ r , L

(l ′ )
α can be written as:

L (l ′ ) α = Lα ′ + n j=1 l ′ p=1 σ α,(k 1 ...kp) j ∂ ∂u (k 1 ...kp) j (11)
Theorem 1. The differential approximation ( 7) is invariant under the group G

(l ′ ) r
if and only if

L (l ′ ) α P λ (x, u, u (k 1 ) , . . . , u (k 1 ...k l ′ ) P λ =0 = 0, α = 1, . . . , r; λ = 1, . . . , q (12)
or

L (l) α F λ + L (l ′ ) α s β=1 m i=1 (h i ) l β R λ i P λ =0 = 0, α = 1, . . . , r; λ = 1, . . . , q (13)
Theorem 1 provides the equations which enable us to obtain the symmetry groups of the differential approximation. The unknowns are the infinitesimal functions η α j , ξ α i and ζ α i , i = 1, . . . , m; j = 1, . . . , n. The infinitesimals σ α,(k1,...,ko) j

, j = 1, . . . , n, are functions of the partial derivatives of η α j and ξ α i . Equation ( 13) is simplified by means of the condition [START_REF] Olver | Geometric foundations of numerical algorithms and symmetry[END_REF]. They lead to an overdetermined system of differential equations, the unknowns of which are the infinitesimal functions.

4. The specific case of the Burgers equation.

4.1.

Symmetries of the Burgers equation. The Burgers equation can be written as:

F (x, t, u, ν, ux, ut, uxx) = ut + u ux -ν uxx = 0 ( 14
)
where ν ≥ 0 is the dynamic viscosity. Denote by G a group of transformations of the Burgers equation in the space E(x, t, u, ν) of the independent variables (x, t), the dependent variable u, and the viscosity ν. The viscosity is taken as a symmetry variable in order to enable us to take into account variations of the Reynolds number.

G is a set of transformations acting smoothly on the space E(x, t, u, ν).

The six-dimensional Lie algebra of the group G is generated by the following operators:

L 1 = ∂ ∂x , L 2 = ∂ ∂t , L 3 = x ∂ ∂x + 2t ∂ ∂t -u ∂ ∂u L 4 = xt ∂ ∂x + t 2 ∂ ∂t + (-ut + x) ∂ ∂u , L 5 = t ∂ ∂x + ∂ ∂u , L 6 = -t ∂ ∂t + u ∂ ∂u + ν ∂ ∂ν (15)
which respectively correspond to:

• the space translation : (x, t, u, ν) -→ (x + ǫ 1 , t, u, ν);

• the time translation

: (x, t, u, ν) -→ (x, t + ǫ 2 , u, ν); • the dilatation : (x, t, u, ν) -→ (ǫ 3 x, ǫ 2 3 t, ǫ -1 3 u, ν); • the projective transformation : (x, t, u, ν) -→ x 1-ǫ 4 t , t 1-ǫ 4 t , xǫ 4 + u(1 -ǫ 4 t), ν ; • the Galilean transformation : (x, t, u, ν) -→ (x + ǫ 5 t, t, u + ǫ 5 , ν); • the dilatation : (x, t, u, ν) -→ (x, ǫ -1 6 t, ǫ 6 u, ǫ 6 ν).
(ǫ i ) i=1,...,6 are constants.

Symmetries of first differential approximations.

Denote by h the mesh size, τ the time step, N x the number of mesh points, N t the number of time steps, and u n i , i ∈ {0, . . . , N t }, n ∈ {0, . . . , N x } the discrete approximation of u(ih, nτ ). In order to shorten the size of the finite difference scheme expressions, we use the following notations introduced by Hildebrand in [START_REF] Hildebrand | Introduction to Numerical Analysis[END_REF]:

δ(u n i ) = u n i+ 1 2 -u n i-1 2 h , µ(u n i ) = u n i+ 1 2 + u n i-1 2 2 δ + (u n i ) = u n i+1 -u n i h , δ -(u n i ) = u n i -u n i-1 h , E α u n i = u n i+α
The Burgers equation can be discretized by means of:

• the FTCS (forward-time and centered-space) scheme:

u n+1 i -u n i τ + µδ h u 2 2 n i -ν δ 2 h 2 u n i = 0
• the Lax-Wendroff scheme:

u n+1 i -u n i τ + µδ h u 2 2 n i -ν δ 2 h 2 u n i + A n i = 0
where:

A n i = - τ 2h 2 E 1 2 u n i δ + u 2 2 n i -E -1 2 u n i δ -u 2 2 n i - ν 2 τ 2 δ 4 h 4 u n i + ντ 2h 3 E 1 2 u n i δ 2 (E 1 2 u n i ) -E -1 2 u n i δ 2 (E -1 2 u n i ) + ντ 2 µδ 3 h 3 u 2 2 n i
• the Crank-Nicolson scheme:

u n+1 i -u n i τ + µδ h u 2 2 n+1 i + u 2 2 n i -ν δ 2 h 2 [u n+1 i + u n i ] = 0
Linear stability properties and the related orders of approximation are displayed in Table 1 (where CF L = aτ h , S = ντ h 2 and S * = ν + ahCF L 2 τ h 2 ).

Scheme Stablility condition Error FTCS

S ≤ 1 2 , CF L ≤ 1 O(τ, h 2 )
Lax-Wendroff

S * ≤ 1 2 , CF L ≤ 1 O(τ 2 , h 2 )
Crank-Nicolson

unconditional stability O(τ 2 , h 2 )
Table 1.

Table of finite difference schemes

Consider u i n as a function of the time step τ , and of the mesh size h, expand it at a given order by means of its Taylor series, and neglect the o(τ α ) and o(h β ) terms, where α and β depend on the order of the schemes. This yields the differential representation of the finite difference equation. The following differential representations are obtained:

• for the FTCS scheme:

ut + 1 2 (u 2 )x -ν uxx + τ 2 g 2 + h 2 12 (u 2 )xxx - νh 2 12 uxxxx = 0
• for the Lax-Wendroff scheme:

ut + 1 2 (u 2 )x -ν uxx + τ 2 6 g 3 + h 2 12 (u 2 )xxx - νh 2 12 uxxxx = 0
• for the Crank-Nicolson scheme:

ut + 1 2 (u 2 )x -νuxx + τ 2 g 3 6 + 1 4 (g 2 1 + ug 2 )x - ν 4 (g 2 )xx + h 2 12 (u 2 )xxx - νh 2 12 uxxxx = 0
where

g 1 = -u 2 2 x + νuxx, g 2 = -g 1 u x + ν g 1 xx , g 3 = -g 2 u -g 2 1 x + ν g 2 xx
Denote by G ′ the group of transformations of a first differential approximation in the space E(x, t, u, h, τ, ν) of the independent variables (x, t) and the dependent variable u, the step size variables (h, τ ) and the viscosity ν.

The l ′ th -prolongation of G ′ can be written as:

L ′(l ′ ) α = ξ α 1 ∂ ∂x + ξ α 2 ∂ ∂t + η α ∂ ∂u + l ′ p=1 σ α,(k 1 ...kp) j ∂ ∂u (k 1 ...kp) j + ζ α 1 ∂ ∂h + ζ α 2 ∂ ∂τ + θ α ∂ ∂ν (16)
where l ′ has been defined in definition 1. Theorem 1 enables us to obtain the necessary and sufficient condition of invariance of the first differential approximation P:

L ′(l ′ ) α P P=0 = 0 ( 17
)
Theorem 1 is applied to the differential representations of the above schemes. The resolution of the determining equations of each first differential approximation yields the 4-parameter group:

ξ α 1 = a + b x, ξ α 2 = c + (2b -d) t, η α = (-b + d) u (18)
ζ α 1 = b h, ζ α 2 = (2b -d) τ , θ α = eν
The 4-dimensional Lie algebra of G ′ is generated by:

L 1 = ∂ ∂x , L 2 = ∂ ∂t , L ′ 3 = x ∂ ∂x + 2t ∂ ∂t -u ∂ ∂u + h ∂ ∂h + 2τ ∂ ∂τ L ′ 4 = -t ∂ ∂t + u ∂ ∂u -τ ∂ ∂τ + ν ∂ ∂ν (19)
These operators are respectively related to:

where (ǫ i ) i=1,...,4 are constants.

The above finite difference equations are preserved by the space translation, the time translation and both dilatations.

Approximating the Burgers equation by the above finite difference equations results in the loss of the projective and Galilean transformations.

5. The invariant scheme.

5.1. Invariant scheme construction. An invariant scheme is constructed in such a way that the related differential approximation preserves the symmetries of the Burgers equation. We propose to approximate the Burgers equation by the following finite difference scheme:

u n+1 i -u n i τ + 1 h µδ - µδ 3 6 u 2 2 n i -ν 1 h 2 δ 2 - δ 4 12 (u n i ) -Ω n i+ 1 2 δ + -Ω n i-1 2 δ -u n i = 0 ( 20
)
where Ω n i = Ω(x i , t n , u n i ) is defined next so that the related differential representation is preserved by the symmetries of the Burgers equation. The scheme has second-order accuracy in space and first-order accuracy in time. The derivatives u 2 2 x and u xx are approximated by fourth order accuracy difference expressions:

µδ h - µδ 3 6h (u n i ) = ux - h 4 30 u 5x n i + O(h 6 ), δ 2 h 2 - δ 4 12h 2 (u n i ) = uxx - h 4 90 u 6x n i + O(h 6 ) (21)
The truncation error of the difference scheme (20) can be written as:

ǫ = τ 2 utt -h 2 Ωux x + O(τ 2 ) + O(h 4 )
u tt is replaced by an expression involving partial derivatives with respect to x, by using the Burgers equation:

utt = (u 2 ux)x -ν(uuxx)x -ν u 2 2 xxx + ν 2 uxxxx ( 22
)
Replacing the previous expression in the truncation error leads to:

ǫ = Cux x - ντ 2 uuxx x - ντ 2
u 2 2 xxx + ν 2 τ 2 uxxxx + O(τ 2 ) + O(h 4 )
where C = τ 2 u 2h 2 Ω. It is convenient for the calculation of C that the truncation error is reduced to:

ǫ = Cux x + O(τ 2 ) + O(h 4 )
The related finite difference scheme is the following first order accuracy in time and second order accuracy in space:

u n+1 i -u n i τ + 1 h µδ - µδ 3 6 u 2 2 n i -ν 1 h 2 δ 2 - δ 4 12 (u n i ) -Ω n i+ 1 2 ∆ 1 -Ω n i-1 2 ∆ -1 u n i + ντ 2 u n i+ 1 2 µδ 2 h 2 (u n i+ 1 2 ) -u n i-1 2 µδ 2 h 2 (u n i-1 2 ) - ν 2 τ 2 δ 4 h 4 u n i + ντ 2 µδ 3 h 3 u 2 2 n i = 0 (23)
and the differential approximation can be written as:

P(x, t, u, ν, ux, ut, uxx) = ut + u ux -ν uxx + (Cux)x = 0 (24)
The von Neumann stability analysis of scheme (23) under a linearized form provides the following necessary conditions for S, CF L and Ω τ = Ωτ :

CF L 2 -2S -2Ωτ ≤ 0, 0 ≤ 4S 3 -2S 2 + Ωτ ≤ 1 2 (25)
If Ω takes is sufficiently close to zero, these conditions become then sufficient for the linear formulation.

5.2. Calculation of the artificial viscosity term. Here we describe the method for determining the artificial viscosity term (Cu x ) x , which is constructed in such a way that the differential approximation (24) is preserved by the symmetries of the Burgers equation. C is a function of the variables (x, t, u, τ, h), and also depends on the partial derivatives of u with respect to x: u x and u xx . C = C(x, t, h, τ, u, u x , u xx ). The necessary and sufficient condition for the differential approximation to be an invariant of the Burgers equation symmetry group is:

L (2) α (ut + u ux -ν uxx) P=0 + L (3) α ((Cux)x) P=0 = 0 (26)
Equation ( 26) provides the determining equations of the symmetry group of equation (24). The determining equations involve partial derivatives of the unknown function C and partial derivatives of the infinitesimal functions of G ′ , which is the symmetry group of the differential representation of the invariant scheme. The infinitesimal functions of G ′ have the following expressions:

ξ α 1 = a + b x + c t + d tx, ξ α 2 = e + d t 2 + (2b -f ) t, ζ α 1 = b h, (27)
ζ α 2 = (2b -f ) τ , η α = c + d x + (-b -d t + f ) u, θ α = f ν
The determining equation with respect to the unknown function C is simplified in using the infinitesimal functions of each subgroup of G ′ . The determining equations of each subgroup of G ′ provides the following linear partial differential equations and the expressions for C:

• the space translation

∂ ∂x C = 0 ⇒ C = C 1 (t, h, τ, u, ux, uxx); • the time translation ∂ ∂t C = 0 ⇒ C = C 2 (x, h, τ, u, ux, uxx); • the dilatation x ∂ ∂x C+2t ∂ ∂t C-u ∂ ∂u C+h ∂ ∂h C+2τ ∂ ∂τ C = 0 ⇒ C = C 3 ( t x 2 , h x , ux, τ x 2 , ux, uxx);
• the projective transformation

∂ ∂x C = 0, ∂ ∂u C = 0, ∂ ∂uxx C = 0, t 2 ∂ ∂t C + 2 ∂ ∂ux C = 0 ⇒ C = C 4 (h, τ, 2+tux t );
• the Galilean transformation ∂ ∂u C + t ∂ ∂x C ⇒ C = C 5 ( ut-x t , t, h, τ, u, ux, uxx);

• the dilatation -t ∂ ∂t C + u ∂ ∂u C + ν ∂ ∂ν C -τ ∂ ∂τ C = 0 ⇒ C = 1 t C 6 (x, h, τ t ).
5.3. Numerical application. The numerical resolution of the Burgers equation has been implemented for scheme (23), the standard schemes (cf. section 4.2) and a scheme with second-order accuracy in time and fourth-order accuracy in space, which is obtained from the invariant scheme when C = 0. The solutions are calculated in the reference frame (F 1) and in the one (F 2) resulting from the Galilean transformations (x, t, u, ν) -→ (x + t, t, u + 1, ν). The artificial viscosity has the following expression:

C = -0.01t(tu -x) 2 (ux) 2 (28) Ω = 1 h 2 ( τ 2 u 2 -C
) is in a sufficiently small neighborhood of zero that we have the sufficiency of conditions (25) for the linear formulation. The problem consists in solving the following differential system:

ut + uux -νuxx = 0, x ∈ [0, 40], t ∈ [0, 20] u(x, 0) = f (x), u(0, t) = g(t), u(40, t) = h(t)
The initial and boundary conditions, f , h, and g are provided by an exact solution of the Burgers equation: In each frame, the numerical solutions are compared to the exact one. The error analysis of the invariant scheme in the reference frame through the features of the truncation error and the graphical representation of the norms of the error (cf. Figures 1, 3 and5) allows to say that the invariant scheme is dissipative and slightly dispersive. The presence of the dissipative term (Cu x ) x in the differential representation of the invariant scheme and the presence of the higher order error terms involving the even-order derivative u 6x (cf. Equation ( 21)) show that the scheme produces numerical damping. Particularly, the amplitudes are not correctly represented for high frequencies, since the solution is subjected to rather rough variation during the first iterations. The dissipation is stronger for Re h = 2, CF L = 0.08 in the reference frame (see Figure 4). Moreover, the presence of higher order error terms involving the odd-order derivative u 5x corresponds to a phase error. The non-invariant schemes are more altered by the change of the frame than the invariant one. Moreover, the invariant scheme appears to be as accurate as the higher order one in the frame (F2).

Figure 4 .Figure 5 .

 45 Figure 4. Space variation of the numerical solutions of the schemes and the exact solution in (F 1) and (F 2). Re h = 2, CF L = 0.08

Figure 6 .

 6 Figure 6. Space variation of the schemes numerical and of the exact solution in (F 1) and (F 2). Re h = 3, CF L = 0.08

  Figures 1, 3 and 5 show the time evolution of the L 2 -norm of the error for the considered schemes, for specific values of the CF L number and the mesh Reynolds number Re h . Figures 2, 4 and 6 display the variations, as functions of the space variable, of the numerical solutions of the considered schemes for the specific value t = 5.

	u(x, t) =	1 + ν 2 √	(x -2t)/(t + 0.1) t + 0.1 exp ((x -2t) 2 /(4ν(t + 0.1)))	+ 2	(29)

  Space variation of the numerical solutions of the schemes and the exact solution in (F 1) and (F 2). Re h = 2, CF L = 0.04 Evolution of the error L 2 -norm in (F1) and (F2). Re h = 2, CF L = 0.08

	L 2 error norm	L 2 error norm
	1	FTCS F1		FTCS F2
		L W F1		L W F2
	0.8	CN F1	1	CN F2
		Invariant F1	0.8	Invariant F2
	0.6	high. ord. F1		high. ord. F2
			0.6
	0.4
			0.4
	0.2		0.2
		Time	Time
	0.2 0.4 0.6 0.8 1 L 2 error norm 1.8 2 2.2 2.4 u x,t 1.6 1.8 2 2.2 2.4 u x,t	2 Figure 1. Evolution of the error L 2 -norm in (F1) and (F2). Re h = 2, CF L = 0.04 4 6 8 10 Time high. ord. F1 Invariant F1 CN F1 L W F1 FTCS F1 2 4 6 8 10 Time 0.2 0.4 0.6 0.8 1 L 2 error norm high. ord. F2 Invariant F2 CN F2 L W F2 FTCS F2 4 6 8 10 12 14 16 18 x high. ord. F1 Invariant F1 Crank Nicol. F1 Lax Wend. F1 FTCS F1 Exact F1 4 6 8 10 12 14 16 18 x 0.6 0.8 1 1.2 1.4 u x,t high. ord. F2 Invariant F2 Crank Nicol. F2 Lax Wend. F2 FTCS F2 Exact F2 4 6 8 10 2 4 6 8 10 Figure 2. 2 6 8 10 12 14 16 18 x high. ord. F1 Invariant F1 Crank Nicol. F1 Lax Wend. F1 FTCS F1 4 6 8 10 12 14 16 18 x 0.6 0.8 1 1.2 1.4 1.6 u x,t high. ord. F3 Invariant F2 Crank Nicol. F2 Lax Wend. F2 FTCS F2 Figure 3. 4 Exact F1 Exact F2


Adel Benamara

Serge Miguet

Mihaela Scuturici

Real-time Multi-Object Tracking with Occlusion and Stationary Objects Handling for Conveying Systems

Multiple object tracking has a broad range of applications ranging from video surveillance to robotics. In this work, we extend the application field to automated conveying systems. Inspired by tracking methods applied to video surveillance, we follow an on-line trackingby-detection approach based on background subtraction. The logistics applications turn out to be a challenging scenario for existing methods. This challenge is twofold: First, conveyed objects tend to have a similar appearance, which makes the occlusion handling difficult. Second, they are often stationary, which make them hard to detect with background subtraction techniques. This work aims to improve the occlusion handling by using the order of the conveyed objects. Besides, to handle stationary objects, we propose a feedback loop from tracking to detection. Finally, we provide an evaluation of the proposed method on a real-world video.

Introduction

Multiple object tracking is a fundamental problem in computer vision; It plays a significant role in many applications, such as video surveillance, robot navigation, human-machine interaction, self-driving cars. Many studies in the recent years focus on tracking pedestrian, cars, bicycles. In this work, we investigate a novel application of multiple object tracking, which is automated conveyor systems. Automated conveyor systems play an important role in many industries, including postal services, packaging, automobile, aeronautics. The purpose of the conveying system is to transport material from one location to another. The controller of the conveying system needs to know the position of the conveyed objects to route them to the right location. Therefore, a robust real-time tracking system is required. Several tracking technologies have been used in the logistics industry, such as photo-electrical sensor network and RFID.

In this work, we follow an on-line tracking-by-detection approach that uses background subtraction techniques for object detection. We address the following problems:

Complex occlusions management: Some methods [START_REF] Rogez | A 3d tracker for groundmoving objects[END_REF] [5] [START_REF] Dziri | Real-time multiple objects tracking on raspberry-pi-based smart embedded camera[END_REF] introduce the notion of group in order to track collectively occluded and occluding objects during occlusions. The model of each member of the group is used at the end of occlusion to achieve correct re-association. In these approaches object reidentification relies on appearance model (color and shape). However, in the logistics scenario, conveyed objects tend to have similar appearance both in color and size (e.g. parcels in the postal application). As a result, these methods fail to recover objects identity. Furthermore, conveyed objects are close to each other, which results in complex group aggregation as depicted in Fig 1 .a. We propose to overcome this issue, an association scheme based on the ordering relation defined by conveyor pathways. In other words, we attempt to exploit the fact that conveyed objects travel serially along the conveyor routes.

Stationary object detection: This is a known issue of background subtraction methods [3] [4]. The absorption of moving objects that become motionless is caused by the adaptation process of background subtraction techniques, where pixels of the stationary objects are progressively included to the background model as shown in Fig 1 .b. In particular, in the logistics scenario, this issue has a dramatic effect on the tracking performance since conveyed objects may be stopped for certain periods of time for a routing purpose. In this paper, we propose a feedback loop from tracking to the detection module; the idea is to stop the adaptation process for regions that correspond to stationary objects. An idea was first presented in [START_REF] Pnevmatikakis | Kalman tracking with target feedback on adaptive background learning[END_REF], we extend the original idea to several state-of-the-art background subtraction techniques [START_REF] Barnich | Vibe: A universal background subtraction algorithm for video sequences[END_REF] [12], we show significant improvements of tracking results in real video sequence. Fig. 2 presents the overview of the pipeline. The proposed algorithm is a trackingby-detection approach. Therefore, the pipeline is decomposed into two main modules: detection and tracking. The acquisition module grabs frames from the camera sensor and feeds sequentially the tracking pipeline with input frames. The first block of the detection module is background subtraction. This block extracts foreground mask by comparing the input frame to the background model. The morphology block removes noise and fills the holes present in the foreground mask. Erosion is first performed to remove noisy pixels, followed by dilation to fill the holes. Connected pixels are labeled using connected component analysis [START_REF] Suzuki | Topological structural analysis of digitized binary images by border following[END_REF], the output of this block are called blobs (a blob is a set of 8-connected pixels). The blobs are filtered, blobs with small size are discarded. Features are extracted from the set of selected blobs, each blob is represented by its center coordinates, bounding box, and color histogram. The tracking module is decomposed into two blocks. The data association block detects associations between the tracked objects and detected blobs. Complex associations are first detected based on the overlapping criteria between objects and blobs. Blobs and objects who are not involved in a complex association are then matched using the Hungarian algorithm based on their centroid. Depending on the type of detected associations between blobs and objects, the tracking management block creates, suppress, associates, recovers objects after splits or merges objects into groups. The model of associated objects is also updated.

The paper is structured as follows: Section 3 presents the tracking algorithm that deals with occlusions. Section 4 explains the tracking to detection feedback mechanism used to stop static objects integration to the background. Section 5 presents extensive results on a real video sequence.

Tracking

Data association

The aim of tracking is to maintain objects identity by resolving at each step the association between tracked and detected objects. In the ideal case, where individual detections are provided for each object in the scene, a greedy one-toone association scheme would be sufficient to link detections reliably. However, in real situations, individual detection cannot be achieved for partial or total occlusions. Especially when a background subtraction technique is used, where a single pixel may connect two distinct blobs into a merged blob. Therefore, we perform data association in two steps. The first step detects complex situations that may result from merging and splitting blobs or in the presence of occlusions between objects. The second step considers only simple cases where individual objects are visible and can be matched with at most one blob. The above observations are valid when the detection step is fast enough to observe a slow evolution of the objects bounding box.

Complexe association

Simple association

For objects (denoted O) and blobs (denoted B) that are not involved in a complex assignment, detected blobs should be assigned to separate and visible objects based on spatial proximity. Hence, we define the cost φ ij of assigning blob j to object i using the Euclidean distance as

φ ij = b c j -o c i if b c j -o c i < τ max , ∞ otherwise (1)
where b c j and o c i are respectively the centroid of blob j and object i, τ max is a distance threshold used to define the maximum allowed displacement of objects between two successive frames. In order to obtain the optimal assignment, we use Hungarian algorithm [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] for computing the assignment matrix A * = [a ij ],a ij ∈ {0, 1}, which minimizes the total assignment cost:

A * = argmin A |O| i=1 |B| j=1 a ij φ ij , s.t. |O| i=1 a ij = 1 , ∀j ∈ {1, . . . , |B|} , |B| j=1 a ij = 1 , ∀i ∈ {1, . . . , |O|} , (2)
Although the Hungarian algorithm is designed for square matrix (|O| = |B|), the algorithm can easily be extended for rectangular cost matrix by padding with impossible cost ∞.

Occlusion handling with ordering

We follow a merge-split approach to deal with occlusions. In this approach, as soon as the beginning of an occlusion is detected, individual objects are frozen, and a group containing these objects is created. The group is tracked as any other objects ( i.e. with its model and attributes ). When a split is detected, the problem is to identify the object that is splitting from the group. Several methods rely on appearance model to re-establish identities at each split. In the context of logistics, objects of interest have a similar appearance, making object re-identification hard. On the other hand, objects are conveyed in procession. Their trajectory follows the conveyor pathways. The ordering of the objects is then preserved during their routing. We propose to use this order to re-identify objects at the end of occlusion. The ordering relation is derived from the conveyor model. A simple polyline that fit the pathways of the conveyor is sufficient to order conveyed objects. This model can be either defined manually during a configuration step or by clustering objects trajectories in a training phase. In our method when a merge occurs, we maintain the object ordering inside the group using the conveyor model. When the group split into several blobs, we use Algorithm-1 to re-establish objects identity. The algorithm finds the group partition that maximizes the likelihood to the blobs involved in the split such as illustrated in Fig 4-a. Blobs are first sorted using the ordering relation. The algorithm processes the ordered objects and blobs sequentially. Let be i the index of the current object and j the index of the current blob. For each blob b j , the algorithm tries to determine if the blob corresponds to a single object or several grouped objects using the area ratio as a clue. The algorithm constructs a new group G j initialized with a single object o i , if the area ratio match with the blob b j , the algorithm moves to the next object and blob. Otherwise, the next objects are added to the candidate group until a match is found. Aside from the sort performed on the blobs, the complexity of the algorithm is O(max(N, M )), where N is the number of objects member of the group and M is the number of blobs detected in the split.

Procedure 1 Association based on order

Input: A group of ordered N tracked objects G = {o1, o2, . . . , oN } that split into M blobs B = {b1, b2, . . . , bM }. Output: Assignments pairs < G | , bj > where {G1, . . . , GM } is M partition of G.

1: procedure Resolve-Split 2: Sort(B) sort the blobs using the ordering relation. 4 Feedback Framework

Feedback of the parameter map

In this section, we will describe the proposed feedback mechanism. We will first discuss background subtraction technique regarding the strategy they follow to update their background model. The update strategy can follow two schemes blind and conservative [START_REF] Barnich | Vibe: A universal background subtraction algorithm for video sequences[END_REF]. Conservative methods update their background model only with pixels classified as background. These methods can indefinitely detect stationary objects. However, pure conservative update scheme leads to everlasting foreground pixel in the case of misclassification. Conversely, blind methods update the background model with pixels whether they have been classified as background or foreground. These methods are not subject to the deadlock situation. However, blind methods are sensitive to slow moving objects.

We have modified the original methods in order to be able to selectively control the update parameters at pixel level rather than a global level. As shown in Table .1, the Gaussian mixture model (GMM) incorporates pixel samples with the learning rate α. We zero out α for stationary objects pixels. ViBe is a conservative method uses a spatial update scheme to incorporate background information. This mechanism can cause the absorption of stationary objects. Therefore, we disable the neighborhood update process for stationary objects pixels. We also adapt the parameter of the method for stationary objects pixels. min is the minimum number of samples of the background model close to the current pixel (the distance in the color space is under R), φ controls the frequency of the update process.

Stationary objects are detected using the speed estimation derived by the Kalman filter associated to each tracked object. An object is considered as stationary when its speed is less than a fixed threshold v sta for c sta successive frames.

Evaluation

In this section, we will discuss the evaluation protocol. We describe first the collected sequence and the methodology used to generate the ground truth dataset. We use a smart embedded camera to perform the video acquisition and the performance evaluation. The embedded platform includes a quad processor ARM Cortex-A9 running at 1GHz with 1 GB of memory. The video was taken with an Omnivision OV5640 at a frame rate of 25 frames per second at VGA resolution (640x480). The video is 4 minutes and 9 seconds long (approximately 6231 frames). The annotated sequence contains 5419 frames. We use a dedicated tool for this task. Objects are manually annotated only on some keyframes by defining object's bounding box, the annotation between key frame is obtained by linear interpolation. The camera is mounted above a conveyor with a bird view angle. The camera capture only a portion of the conveyor system. We use 5 parcels in the experiments: 4 yellow boxes with similar color and size and one white rectangular parcel. The objects are circulating on the conveyor system in a loop. 22 objects are annotated in the whole sequence, including a human operator that manipulates a blocked parcel.

We use the CLEARMOT [START_REF] Bernardin | Evaluating multiple object tracking performance: the clear mot metrics[END_REF] metrics, particularly : the Multiple Object Tracking Accuracy(MOTA), Multiple Object Tracking Precision (MOTP), False Positive (FP), Missed Object (Miss), ID Switch (IDs). To measure the improvement of each contribution, we run the tracking pipeline with different settings. We execute the pipeline with the proposed tracking-to-detection feedback enabled and disabled for Vibe and GMM. In conjunction with the two possibles occlusion handling methods : ours with order and appearance as a baseline. For occlusion handling with appearance, we use a 2D histogram and quantize the H,V components of HSV color space using 3 and 2 bits respectively.

Quantitative results of our experiments are listed in Table 2. As expected GMM is more sensitive to stationary objects, the feedback loop demonstrates a significant improvement to lower missed objects (i.e. stationary objects). ViBe, on the other hand, is less sensitive to stationary objects due to the conservative update scheme. However, the inhibition of the spatial update mechanism with the feedback loop has also a significant effect in lowering missed objects. We achieve significant improvements with the proposed occlusion handling with ordering relation, we lower both id switches and missed objects in comparison to the appearance method. Indeed, the lowering of id switches is explained by the presence of identical yellow boxes, which makes the appearance method fail to recover their identity reliably. On the other hand, the lowering of missed objects is explained by the ability of our method to handle complex occlusion situations, where a group splits into several smaller groups. In this situation, the appear- ance method can only handle blobs that correspond to individual objects, the blobs that correspond to more than one objects will not be matched and then generate misses.

Conclusion

We have proposed a multi-object tracker adapted for conveying systems. We have proposed an occlusion handling method that exploits the ordering of the object to reliably re-establish objects identity in complex occlusion situations. We have also addressed stationary objects detection with by introducing a feedback loop from tracking to detection. Our evaluations on real data demonstrate significant improvements of tracking results in logistics scenario compared to state-of-theart approach for occlusion handling. In future works, we will extract the ordering relation by learning conveyed objects trajectories rather than relying on manual configuration.

Fig. 1 .

 1 Fig. 1. Examples of the addressed issues.

Fig. 2 .

 2 Fig. 2. Tracking-by-detection pipeline with tracking feedback to detection.

Fig. 3 .

 3 Fig. 3. Split and merge detection

  Fig 4b illustrates the conveyor model used during the experiments.

Table 1 .

 1 Background subtraction methods and parameters map generation. The (x, y) subscript is left for ViBe in the first row, Mstationary is a binary mask generated based on the union of the bounding boxes of objects detected as Stationary.

	3:	j = 1
	4:	i = 1
	5:	while i ≤ |G| and j ≤ |B| do
	6:	Gj = {oi}
	7:	while M atch(Gj, bj) = true and i + 1 ≤ |G| do
	8:	Gj = Gj ∪ {oi+1}
	9:	i = i + 1
	10:	end while
	11:	j = j + 1
	12:	i = i + 1
	13:	end while
	14: end procedure
	15: function M atch(Gj, bj)
	16:	rG j = area(Gj)/area(G)	compute the contribution of Gj to group G.
	17:	r b j = area(bj)/area(B)	compute the contribution of bj to blobs set B.
	18:	if rG j ≈ r b j then
	19:	return true
	20:	else
	21:	return false
	22:	end if
	23: end function

static , static , φ static ), if (x, y) ∈ Mstatic (R bg , bg , φ bg ), otherwise GMM [12] [10] Blind α(x, y) = 0, if (x, y) ∈ Mstatic α bg , otherwise

Table 2 .

 2 Tracking performance. MOTA (higher is better), MOTP (lower is better), FP false positive (lower is better), Miss (lower is better), IDS id switch (lower is better).

	Occlusion Handling	Methods	MOTA [%] MOTP [px] FP MISS IDS
		ViBe	80.4237	65.09	515 2223 34
	Appearance	ViBe + Feedback 83.5805	60.98	566 1731 28
		GMM	61.9845	81.18	607 4737 39
		GMM + Feedback 84.7811	64.86	772 1357 26
		ViBe	86.815	67.13	742 1098 27
	Order	ViBe + Feedback 88.6017	56.14	1577 23 14
		GMM	66.363	81.43	111 4043 37
		GMM + Feedback 82.3517	70.31	577 1897 25


Ameni Sassi
email: ameni.sessi.tn@ieee.org

Ben Chokri
email: chokri.benamar@ieee.org

Amar

Serge Miguet
email: serge.miguet@univ-lyon2.fr

Skyline-based Approach for Natural Scene Identification

Keywords: skyline, geometric descriptor, landscape, curvature

The skyline, defined as the line separating the sky from other objects on the ground, could provide unique and useful information for a variety of applications. This line was used as a key data, especially, for geo-localization and aerial robotic applications. The particular shape and the geometric features of a skyline may be the identity of the landscape itself. The skyline, once well extracted, could show the silhouette of a famous tower, the mountain peaks, or the landscape topography. In this paper, we proposed a geometric description of the extracted skyline from landscapes. Based on some geometric descriptors, we tried to pick up practical measurements for each skyline. The first proposed approach was the straight lines' classification to differentiate between urban and natural landscapes from their horizon line. The second one is the curvature analysis using a Curvature Scale Space descriptor. This descriptor was used to enhance the first one and to distinguish between natural part and buildings in the same skyline. The results obtained from these geometric description tools were very competitive and they will be the inputs for a classification process.

I. INTRODUCTION

Watching a cityscape from high rise, we can notice that besides the indentation created by tall buildings, the "nature" elements such as vegetation and topographical context may increase the complexity of urban landscape representations. The appearance of these natural elements in a cityscape could affect the reviews of citizens or visitors about that landscape. Our distant goal is to objectify the influence of natural elements on the representations of urban landscapes using a photo-questionnaire method. This was in the context of a sustainable city project that aims to identify and systematically analyze the skyline perceptions of practitioners and the general public by matching the characteristics taken from skyline photographs with the perceptions of several audiences in some European cities (Lyon and London). The first characteristic identified and extracted from landscapes was the skyline, defined as the line that marks the boundary between the sky and the ground elements (natural or artificial). Based on the obtained skyline, we will try to extract other useful measurements to describe that line and use it for our purposes. Our first goal is the identification of natural elements within the skyline. In fact, this line can give a preliminary idea about the existing elements in a landscape. We can, for example, know that there are tall buildings or large mountains by scanning the skyline of a specific photography. The Fig. 1 shows the skyline of Chicago where we can easily discriminate the towers and the tall buildings. Because of the unique aspects that could contain a skyline, we have chosen it to be the first tool to interpret the images of landscapes. This paper is mainly divided into four parts. The first part reveals some related works about how the skyline was used for different applications. Then, in the second part, we will describe the geometric description tools used to analyze the skyline of landscapes. After that we present some first results in the third section. The last section of this paper contains the conclusions and some perspectives that will be the goal of a future work.

II. RELATED WORKS

The fact that the horizon line is highly informative was proved in the literature by using the information extracted from skylines for different purposes, such as navigation in airborne [START_REF] Bazin | Dynamic programming and skyline extraction in catadioptric infrared images[END_REF][2][4] and marine vehicles [START_REF] Gershikov | Horizon line detection in marine images: which method to choose?[END_REF], visual geo-localization [START_REF] Ramalingam | Geo-localization using skylines from omni-images[END_REF][6] [START_REF] Baatz | Large scale visual geo-localization of images in mountainous terrain[END_REF] and also annotation of landscapes [START_REF] Zhu | Adaptive rendering for large-scale skyline characterization and matching[END_REF] [START_REF] Liu | Automatic peak recognition for mountain images[END_REF].

For vision-based navigation, many researchers [START_REF] Bazin | Dynamic programming and skyline extraction in catadioptric infrared images[END_REF][2] [START_REF] Gershikov | Horizon line detection in marine images: which method to choose?[END_REF] suggested to use the skyline as a key information, instead of the GPS and the INS, to estimate the position and the orientation of the Unmanned Arial Vehicles UAV. Using vision-based methods to extract the skyline and some useful points on it, and then matching the extracted data with a Digital Elevation Map (DEM) could be the only required steps to estimate accurately the position and the attitude of the UAV camera. The challenge in these methods is to extract the skyline and the appropriate features under noisy images, bad weather and variable luminosity conditions. Reference [START_REF] Woo | Robust horizon and peak extraction for vision-based navigation[END_REF] proposed a vision-based navigation system to estimate the localization of UAV for a mountain-region use. The purpose was to accurately extract the skyline and the position of some key peaks on the mountain from infrared images. The use of IR (infrared) images was a solution for some challenges, but the narrow field of view of traditional cameras restricts the viewed part of the skyline. To deal with this problem, the authors in [START_REF] Bazin | Dynamic programming and skyline extraction in catadioptric infrared images[END_REF] proposed a challenging algorithm to extract the skyline in catadioptric infrared images. Compared to traditional infrared cameras, the catadioptric IR systems could cover a much larger part of the horizon scene which improve the matching with the Digital Elevation Map. Another type of images was used in [START_REF] Stone | Sky segmentation with ultraviolet images can be used for navigation[END_REF], that proposed a sky segmentation using panoramic ultraviolet UV images.

Since the skyline could provide robust and unique features, we can recognize the area where a photo was taken by looking only at the skyline. So, extracting the skylines and matching them to a database of known skyline shapes or digital elevation maps could be an effective tool for automated geo-localization of images. The geo-localization can be applied for urban landscapes [START_REF] Ramalingam | Geo-localization using skylines from omni-images[END_REF] or also for natural ones [6] [7]. Reference [START_REF] Ramalingam | Geo-localization using skylines from omni-images[END_REF] relies only on a skyline extraction algorithm from omnidirectional images for an accurate geo-localization within the streets of a city. The extracted skyline is matched to the ones extracted from 3D city models to estimate the position of the omnidirectional camera. Another research work on urban geolocalization [START_REF] Bansal | Geometric urban geolocalization[END_REF], proposed to use the building skyline, some useful points on it and the skylines pre-rendered from digital elevation maps, to estimate the camera position. The Digital Elevation Models and 3D models have shown promising results for the geo-localization in natural landscapes as well. In this work [START_REF] Baatz | Large scale visual geo-localization of images in mountainous terrain[END_REF], the authors have proposed a work for geolocating queries in a large scale mountainous terrain by matching against skylines pre-rendered from digital elevation models. These works of vision-based geo-localization are almost useful when the GPS information is not available and they could be helpful for geo-tagging applications. In this geotagging application [START_REF] Zhu | Adaptive rendering for large-scale skyline characterization and matching[END_REF], they have proposed an automatic extraction algorithm of the skyline to match it to a database of reference skylines in order to geo-tag photos and images in natural terrain. The extracted skyline was used here [START_REF] Liu | Automatic peak recognition for mountain images[END_REF] to annotate the positions and names of mountains peaks.

As proved in the previously-mentioned research works, the skyline could be the key data for a variety of applications. The way of operating this line depends on the requirements of each application. In our work, the main aspects used to interpret the skyline were the specific geometric features that differ from an horizon line to another.

III. SKYLINE'S GEOMETRIC DESCRIPTION

Based on the extracted skyline, we tried to take out useful measurements to distinguish between urban landscape photos and natural landscape photos and thereafter the natural part and the buildings one in the same landscape. As we can notice in the Fig. 2, buildings and towers are generally represented by regular geometric forms and straight segments, unlike the natural elements that tend to occur more irregular. The buildings are mostly represented by long segments, whereas the natural skyline part is often composed by many short segments. If we define the corner as the intersection between straight segments, then the number of corners is not the same in the two parts. Moreover, the curvature seems to have dissimilar values all along the skyline; which could be a way to detect some key points on the skyline.

The skyline was seen as an open curve or a polyline. Thus, some geometric descriptors could be applied to reach our goal. The first geometric description tool was based on the length of straight segments. The second geometric descriptor was the Curvature Scale Space representation to analyze the horizon line point per point.

A. The straight lines classification

As we have mentioned before, the buildings and towers in landscapes are characterized by long straight lines, especially when there is a high tower. So, the idea was to break the skyline into segments through a polygonal approximation algorithm. The skyline is, then, approximated into a fewer number of vertices and straight segments depending on a precision parameter. The value given to this parameter must preserve the small details of the approximated skyline; it is supposed to eliminate the noise and obtain truly straight segments. The next step was to classify the straight lines, obtained from the polygonal approximation, according to their length. The result of this process is a histogram that presents the number of segments on each category of segments length. To avoid getting categories with zero segments and to have an appropriate histogram to each skyline and landscape photography, we thought about making variable categories' intervals. The limits of intervals follow a pattern of geometric sequence. The width of the interval of each category depends on three elements: (i) the number of categories itself categ, (ii) the shortest segment l in the approximated skyline and (iii) the height of the landscape image. The equations (14) explain better the choice of categories limits.

[ ] The comparison between the number of straight lines in each category for different skylines, could be a way to differ urban landscapes from natural ones. Nevertheless, this manner of comparison would not be adequate since the width of each category is not the same for all landscapes. Therefore, the idea of classifying straight lines using histograms was manipulated to obtain an indicative value of segments' length for each skyline. This step will be better explained in section IV-A while showing the results.

1 0, interval l = [ ] 2 , interval l l lx = + 2 3 , interval l lx l lx lx   = + + +   ........ 1 1 1 , 1 1 categ categ categ x x interval l l x x -       - - =       - -      

B. The curvature analysis

Looking at a skyline, we can affirm that the degree of curvature goes through a lot of changes along this line. It is obvious that the indentation created by a tree, the mountain peaks and the structure of a building seem to have different values of curvature. The formula to calculate the curvature, we admitted, is [START_REF] Ramalingam | Geo-localization using skylines from omni-images[END_REF] where the skyline is defined as the curve [START_REF] Bansal | Geometric urban geolocalization[END_REF]. The representation of the curvature along the skyline gives a graphic as shown in the Fig. 3; the curvature value could be positive, negative or null. The sudden changes in these values label the key points on the skyline.

( )

3/2 2 2 ( ) ( ) ( ) ( ) ( ) ( ) ( )
x t y t x t y t k t

x t y t

- = +       (5)
( )

( ) ( ), ( ) t x t y t τ =  (6)
To better discover the geometric aspects along the skyline and select the persistent key points, we have proceeded to the curvature scale-space description. The CSS (Curvature Scale Space) descriptor was introduced by F. Mokhtarian and A.K. Mackworth [START_REF] Mokhtarian | A Scale-based description and recognition of planar curves and two-dimensional shapes[END_REF] and it can be used for shape similarity retrieval [START_REF] Mokhtarian | Shape similarity under affine transforms[END_REF] and also for corner detection [START_REF] Mokhtarian | Robust image corner detection through curvature scale space[END_REF]. The process consists on computing the curvature values for a curve that have to be smoothed gradually using Gaussian kernels that precise the scale. So, at each scale we have different values of curvatures corresponding to the smoothed curve. As we have said before, the curvature goes through negative values and positive ones. The zero-crossings of curvature at different scales are the way to find the concavities/convexities of a curve which are helpful to find some perceptual features within this curve. Fig. 4 explains the output of the CSS descriptor for a closed contour. The CSS map illustrates the lifetime of the curve's points all over the increasing scales. The zero-crossings of the curvature could be picked from the CSS image for each scale. The application of the CSS descriptor to the skyline gives important information that could be practical to take out the geometric features of each part on the skyline. The manner how we analyze the curvature and the results obtained are illustrated in the section IV-B.

IV. EXPERIMENTAL RESULTS

A. Lengths histograms and segments 'mean size'

The classification of straight lines composing the skyline is shown as histograms with variable categories as indicated before. The obtained histograms for two different landscapes are shown in Fig. 8. We can notice that the numbers of segments in each category vary from a landscape to another. For New York skyline, the histogram shows that there is a considerable number of long straight segments unlike the histogram generated for the natural landscape skyline that depicts a very small number of long segments and a very important number of short segments. To describe the skyline with a concrete value, we decided to generate other histograms with different distribution. The idea was to make a link between the number of segments for each length and the length itself, in order to highlight the relation between the number of segments and its length. The axes and the distribution of these histograms are more explained in the Fig. 5 and Fig. 6. The cumulated histograms, obtained from the histograms showing the number of segments per each segment size multiplied by the size, give an idea about the ' mean size' of segments. We expect to have a higher value for urban landscapes than natural ones. The mean segments size, derived from the cumulated histogram, could be an index for describing the skyline segments in a general way: the skyline is almost natural or artificial. The mean segments sizes that we got for some landscapes are listed in Table 1. Observing these values, we can notice that urban landscapes are described by a high index of segments size whereas the natural ones show smaller values. This skyline segments classification was not crucial for some urban landscapes where the photo is taken from faraway or there are no tall buildings like the case of the photo of the Istanbul landscape (Fig. 7). Moreover, these histograms have given an idea about the entire skyline and could not help to specify which part represents a natural element and which one shows buildings.

B. The evolution of skyline key points based on curvature

The first step of the CSS description algorithm is computing the CSS map, starting with re-sampling the skyline into a fixed number of points, to calculate the curvature for every point at different increasing scales. After that, the values of curvature are scanned at each scale to select the points where the curvature passes through zero. These points are called the zero-crossings points and they are the dots composing the CSS map. The Fig. 9 shows an example of a CSS map generated for a skyline where the blue curves represent the zero-crossing points of curvature over different scales of smoothing. It is clear from the figure that the number of zero-crossings at high scales is fewer than the ones found in low scales since the curve is very well smoothed in coarse scales and there will be few concavities and convexities. Applying the curvature descriptor to our extracted skylines, we obtained interesting results. Fig. 11 (a)(b) depict a zoom-on parts of the colored skylines; the points where the curvature is almost constant are in yellow (straight lines) and the red or green points are the points where the skyline change the direction to a concavity point (Green) or a convexity one (Red). Observing the CSS map and the colorful skyline over different scales for two different landscapes (urban and natural), we can conclude that this descriptor was useful for tracking the evolution of key points on the skyline. The Fig. 11 shows a zoom-on the colored skyline for the New York landscape and the natural landscape at the scale 1 (a)(b) and the scale 5 (c)(d) behind their generated CSS map. What we can notice from these images is that there is an interesting number of key points (concavity and convexity) in the natural landscape at fine scale that disappear quickly by increasing the scale of smoothing (scale 5). In contrary, a visible number of inflexion points still appear at the same scale for the urban landscape. By comparing the skylines obtained in different scales for a variety of landscapes, we can notice that for the natural landscapes there was an interesting number of key points that reflect the fast variation in the curvature but these points disappeared fast from low scales. In contrast, the remarkable points in urban skylines, such as the buildings corners, resist until high scales. To translate the observations with concrete values, we chose to represent the waning of the number of key points across the different scales through scatter plots as shown in the Fig. 10. To describe these graphics using numbers, we calculated the percentage of key points that resist until a specific scale. The smoothing scale we have chosen to compute these rates was the middle of the scales' interval. To justify our observations, these percentages for cityscapes should be higher than the ones obtained for natural landscapes. The results for some landscapes of our test images are summarized in Table 2. The percentage of key points on the skyline of natural landscapes that resist until the middle scale seems to be lower than the percentage for urban landscapes; which prove that the small serration caused by natural elements on the skyline disappears from low scales of smoothing. The images collected in the Fig. 12 show the results obtained for four different landscapes: Shanghai landscape, Chicago landscape and two natural landscapes. The histograms presenting the classification of straight lines in the skyline are shown in Fig. 12 (b). The number of categories has been set to five, unlike the width of each category which is variable and depends on parameters extracted from the processed photo. The zoom on the colorful skyline resulted from the curvature analysis is shown in Fig. 12 (c) for the different landscapes. After the curvature analysis, the CSS maps are depicted in the Fig. 12 (d) to describe the evolution of skyline points across the increasing scales of smoothing. The scatter plots in Fig. 12(e), as mentioned before, are made in order to visualize the evolution in the number of the skyline key points all over the different scales. These plotted curves prove the short lifetime of the key points across the scales in the natural landscapes, in the contrary, the points in urban landscapes persist until high scales.

V. CONCLUSIONS AND FUTURE WORK

In this paper, we aimed to introduce a new approach for the processing of landscapes' photos based on the skyline. Besides its usability in different applications, the horizon line, once extracted, is the same in different weather and lighting conditions. So, if we do extract useful information from the skyline, we can minimize the energy spent for searching information otherwise and overpass some problems like: luminosity, ambiguous colors or textures. Extracting some geometric features from the skyline was our first step to analyze this curve. The straight lines' classification gave an idea about the whole skyline and especially if the background in a landscape is almost natural or artificial. To enhance this classification, we have used the curvature scale space descriptor which analyze the skyline point per point.

Our work with the curvature descriptor should be expanded by making a tracking algorithm for some key points in the skyline. Thus, we will try to answer some questions like: Until which scale the buildings' corners still appear? Which points disappear from very low scale levels? The measurements picked from landscapes using the geometrical description will be the entry data for a classification algorithm. The results obtained using our skyline-based approach to identify natural scenes will be compared with similar works in the literature through a future work.

Fig. 1 .

 1 Fig. 1. Chicago's Skyline

Fig. 2 .

 2 Fig. 2. Skyline's parts

Fig. 3 .

 3 Fig. 3. The curvature function

Fig. 5 .Fig. 4 .

 54 Fig. 5. Simple cumulated histogram

Fig. 7 .Fig. 6 .Fig. 8 .

 768 Fig. 7. The landscape of Istanbul

Fig. 9 .

 9 Fig. 9. The CSS image

Fig. 10 .

 10 Fig. 10. The number of skyline key points over smoothing scales

Fig. 11 .Fig. 12 .

 1112 Fig. 11. (a)(b) Zoom on the Colorful skyline at the scale 1(c)(d) Zoom on the Colorful skyline at the scale 5 (e)(f) The CSS MAP (g)(h) The evolution of the number of skyline key points over the scales

TABLE I .

 I THE MEAN LENGTH INDEX FOR DIFFERENT LANDSCAPES

	Landscape	Index of the Mean Segments size
	New York	63
	Boston	61
	Dubai	55
	Shanghai	35
	Chicago	53
	Madrid	58
	London	28
	Istanbul	25
	Natural Landscape2	22
	Natural Landscape4	30
	Natural Landscape5	32
	Natural Landscape6	38
	Natural Landscape8	21
	Natural Landscape9	25
	Vancouver	32

TABLE II .

 II THE PERCENTAGE OF RESISTING POINTS FOR DIFFERENT LANDSCAPES

Landscape % Percentage of resisting key points at the middle scale

		London	6.83
		Istanbul	6.66
		Natural Landscape2	3.63
		Natural Landscape4	5.71
		Natural Landscape5	3.82
		Natural Landscape6	4.73
		Natural Landscape8	5.69
		Natural Landscape9	4.11
		Vancouver	11.2
	New York	7.56
	Boston	11.68
	Dubai	6.36
	Shanghai	8.08
	Chicago	8.65
	Madrid	16.12

ACKNOWLEDGEMENTS

The authors would like to acknowledge the financial support of this work by grants from General Direction of Scientific Research (DGRST), Tunisia, under the ARUB program. This work was initiated in the project SKYLINE which is led by EVS (Lyon2), with the following participants: l'EIVP-Ville de Paris, le LIRIS (Lyon1-2) and "l'Agence d'Urbanisme de Lyon."