Pierre Dérian
email: pierre.derian@inria.fr

Patrick Héas
email: patrick.heas@inria.fr

Cédric Herzet
email: cedric.herzet@inria.fr

Étienne Mémin
email: etienne.memin@inria.fr

Wavelet-Based Fluid Motion Estimation

Keywords: Optical flow, continuous fluid motion, wavelet multi-resolution analysis, high-order regularization, polynomial approximation

Based on a wavelet expansion of the velocity field, we present a novel optical flow algorithm dedicated to the estimation of continuous motion fields such as fluid flows. This scale-space representation, associated to a simple gradient-based optimization algorithm, naturally sets up a well-defined multi-resolution analysis framework for the optical flow estimation problem, thus avoiding the common drawbacks of standard multi-resolution schemes. Moreover, wavelet properties enable the design of simple yet efficient high-order regularizers or polynomial approximations associated to a low computational complexity. Accuracy of proposed methods is assessed on challenging sequences of turbulent fluids flows.

Introduction

Recent years have seen significant progress in signal processing techniques for fluid motion estimation. The wider availability of image-like data, whether coming from experimental facilities (e.g. particle image velocimetry) or from largerscale geophysical study systems such as lidars or meteorological and oceanographical satellites, strongly motivates the development of computer-vision methods dedicated to their analysis. Correlation-based and variational methods have proven to be efficient in this context. However, the specific nature of fluid motion highly complicates the process. Indeed, one has to deal with continuous fields showing complex structures evolving at high velocities. This is particularly problematic with optical flow methods, where the problem non-linearity requires to resort to an ad-hoc multi-resolution strategy. Although leading to good empirical results, this technique is known to have a number of drawbacks. Moreover, the underdetermined nature of the optical flow estimation problem imposes to add some prior information about the sought motion field. In many contributions dealing with rigid-motion estimation, first-order regularization is considered with success. However, when tackling more challenging problems such as motion estimation of turbulent fluids, this simple prior turns out to be inadequate. Second-order regularizers allowing to enforce physically-sound properties of the flow are considered [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Heas | Bayesian selection of scaling laws for motion modeling in images[END_REF][START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF][START_REF] Yuan | Simultaneous higher-order optical flow estimation and decomposition[END_REF], but their implementation raises up several issues.

In this paper, we propose an optical-flow estimation procedure based on a wavelet expansion of the velocity field. This approach turns out to offer a nice mathematical framework for multi-resolution estimation algorithms, which avoids some of the drawbacks of the usual approach. Note that algorithms based on wavelet expansion of the data [START_REF] Bernard | Wavelets and ill posed problems: optic flow and scattered data interpolation[END_REF] or the velocity field [START_REF] Wu | Image registration using wavelet-based motion model[END_REF] have been previously proposed. However, unlike the algorithm presented hereafter, the computational complexity of the later seriously limits its application to small images and/or the estimation of the coarsest motion scales. Moreover, we consider the effective implementation of high-order regularization schemes, based upon very simple constraints on the wavelet coefficients at small scales. We finally assess the relevance of proposed methods on challenging image sequences of turbulent fluid motions. Simulation results prove that the proposed approach outperforms the most effective state-of-the-art algorithms.

Optical Flow Background

Optical flow estimation is an ill-posed inverse problem. It consists in estimating the apparent motion of a 3D scene through image brightness I(x, t) variations in space x = (x 1 , x 2 ) T ∈ Ω ⊂ R 2 and time t ∈ R. Optical flow, identified by a 2D velocity field v(x, t) : Ω × R + → R 2 is the projection on the image plane of the 3D scene velocity. Its estimation involves two main aspects: a data model that links image data to the velocity field and a regularization scheme to overcome the ill-posedness.

Non-Linear Data Model

Data models are commonly built upon assumptions about the temporal variations of the image brightness. The integration of a conservation assumption leads to the well-known Displaced Frame Difference (DFD) equation, which is studied in the following. However, the approach remains valid for any other integrated data model. Let us denote by I 0 (x) and I 1 (x) two consecutive image samples of the continuous sequence I(x, t) which has been discretized in time with a unit interval. Under rigid motion and stable lighting conditions, v = (v 1 , v 2 ) T satisfies the standard DFD equation, which reads:

∀x ∈ Ω, f DFD (I, v) = I 1 (x + v(x)) -I 0 (x) = 0 . (1)
The estimated motion field v is obtained by minimizing a cost function, which we chose quadratic in the following to clarify the presentation:

v = arg min v J DFD (I, v) , with J DFD (I, v) = 1 2 Ω |f DFD (I, v)| 2 dx . (2)
The data model being non-linear w.r.t. the velocity field v, estimation of optical flow therefore requires a specific optimization approach.

Classical Multi-resolution Strategy

Indeed, the equations for the inversion are only valid if the solution remains in the linearity region of the image intensity function. A standard approach for tackling non-linearity is to rely on an incremental multi-resolution strategy. This approach consists in choosing some sufficiently coarse low-pass-filtered version of the images at which the linearity assumption is valid, and to estimate a first displacement field assumed to correspond to a coarse representation of the motion. Then, a so-called Gauss-Newton strategy is used by applying successive linearizations around the current estimate and warping accordingly a representation of the images of increasing resolution. More explicitly, let us introduce the following incremental decomposition of the displacement field at resolution1 2 j :

v j = ṽj + v ′ j ( 3
)
where v ′ j represents the unknown incremental displacement field at resolution 2 j and ṽj i<j P j (v ′ i ) is a coarse motion estimate computed at the previous scales; P j (v ′ i ) denotes a projection operator which projects v ′ i onto the grid considered at resolution 2 j . In order to respect the Shannon sampling theorem, the coarse scale data term is derived by a low-pass filtering of the original images with a kernel2 G j , followed by a subsampling at period 2 j . Using (3), at coarse scale, image I j (x) and the motion-compensated image Ĩj (x) are then defined as:

I j (x) =↓ 2 j • (G j ⋆ I 0 (x)) Ĩj (x) =↓ 2 j • (G j ⋆ I 1 (x + ṽj (x))) , (4)
where ↓ 2 j denotes a 2 j -periodic subsampling operator. It yields a functional J j OBS defined as a linearized version of (1) around ṽj (x):

J j OBS (I j , v ′ j ) = 1 2 Ωj Ĩj (x)-I j (x) +v ′ j (x) • ∇ Ĩj (x) 2 dx . (5)
Finally, the sought motion estimate v is given by solving a system of coupled equations associated to resolutions increasing from 2 C to 2 F :

         v = v ′ F + ṽF = v ′ F + F -1 i=C P F (v ′ i ) , v ′ j = arg min v ′ J j OBS (I j , v ′ ) , ∀j ∈ {C, • • • , F } . (6)
where the finest scale s = 2 -F corresponds to the pixel whereas the coarsest scale is noted s = 2 -C . In practice, equations in (6) are usually solved independently, starting from the coarsest to the finest scale. This coarse-to-fine approach has the drawback of freezing (i.e. leaving unchanged), at a given scale, all the previous coarser estimates. Moreover, the major weakness of this strategy is the arbitrary approximation of the original functional in (2) by a set of coarse scale data terms [START_REF] Heitz | Final report on the evaluation of the tasks of the workpackage 2, FLUID project deliverable 5.4[END_REF], which are defined at different resolutions by a modification of original input images with (4) and by a linearization of model ( 1) around the previous motion estimate. In the next section, we will see that this multi-resolution strategy has a mathematically-sound formulation within the framework of wavelet representations.

The Aperture Problem and Usual Regularization Schemes

Previously introduced data models remain under-constrained, as they provide for each time t a single equation for two unknowns (v 1 , v 2 ) at each spatial location

x = (x 1 , x 2 ) T .
To deal with this under-constrained estimation problem, the most common setting consists in enforcing some spatial coherence to the solution.

Implicit Regularization

The motion field is constrained to be of the form v = Φ(Θ), where Φ is a function parametrized by Θ (piece-wise polynomial functions are often used). Implicit regularization schemes penalize discrepancies from model ( 1) by minimizing J DFD with respect to Θ, i.e.

v = Φ arg min Θ J DFD (I, Φ(Θ)) . (7)
Associated to a low-order parametric representation, this simple approach reduces drastically the dimension of the problem, hence addressing its underconstrained nature. However, when spatiotemporal gradients of the images vanish, it is impossible to guarantee the existence of an unique solution: this is the aperture problem.

Explicit Regularization Global regularization schemes in their simplest form define the estimation problem through the minimization of a functional composed of two terms balanced by a regularization coefficient µ > 0:

J(I, v, µ) = J DFD (I, v) + µJ reg (v) . (8)
Thus, motion estimate v satisfies v = arg min v J(v, I, µ). The data term J DFD is still defined by [START_REF] Corpetti | Dense estimation of fluid flows[END_REF] . The second term, J reg (the "regularization term"), encourages the solution to follow some prior smoothness model formalized with function f reg :

J reg (v) = 1 2 Ω f reg (v, x)dx . (9)
An n-order regularization writes in its simplest form:

f reg (v, x) = i=1,2 j=1,2 ∂ n v i ∂x n j (x) 2 . (10)
A first-order regularizer (i.e. n = 1) enforcing weak spatial gradients of the two components v 1 and v 2 of the velocity field v is very often used [START_REF] Horn | Determining optical flow[END_REF]. Secondorder regularizers (i.e. n > 1) have been proposed in the literature in the case of fluid flows [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF][START_REF] Yuan | Simultaneous higher-order optical flow estimation and decomposition[END_REF]. However, since motion variables are considered on the pixel grid, an approximation of continuous spatial derivatives by discrete operators is required. For regular pixel grids, it is usually done using finite difference schemes. Nevertheless, it is well known that ensuring stability of the discretization schemes of high-order regularizer may constitute a difficult problem.

Wavelet Formulation

As shown in Sect. 2, the common optical flow estimation approach suffers from two main drawbacks: the necessary "empirical" multi-resolution approach and the implementation of efficient regularizations terms. The use of wavelet bases is a simple answer to both problems. Moreover, it has been shown that a wavelet expansion is appropriate for representing turbulent flows [START_REF] Farge | Wavelet transforms and their applications to turbulence[END_REF].

Wavelet Decomposition

In order to avoid the limitations of the classical multi-resolution strategy, we consider the projection of each scalar component v 1 , v 2 of the velocity field v onto multi-resolution approximation spaces exhibited by the wavelet formalism. Let us introduce briefly this context for real 1D scalar signals. We consider a multiresolution approximation of L 2 (R) as a sequence {V j } j∈Z of closed subspaces, so-called approximation spaces, notably verifying3

V j ⊂ V j+1 ; lim j→-∞ V j = +∞ j=-∞ V j = {0} ; lim j→+∞ V j = Closure   +∞ j=-∞ V j   = L 2 (R) .
Since approximation spaces are sequentially included within each other, they can be decomposed: V j+1 = V j ⊕ W j . Those W j are the orthogonal complements of approximation spaces, they are called detail spaces.

Practically, scalar 1D signals being finite, they belong to a given approximation space according to their resolution, i.e. number of samples. Let w be a 1D signal of 2 F +1 samples, then w ∈

V F +1 = V C ⊕W C ⊕W C+1 ⊕• • •⊕W F ⊂ L 2 ([0, 1]),
where 0 ≤ C ≤ F . The projection of w on this multiscale basis writes:

w(x) = 2 C -1 k=0 w, φ C,k L 2 φ C,k (x) + F j=C 2 j -1 k=0 w, ψ j,k L 2 ψ j,k (k) . (11)
Here, {φ C,k } k and {ψ j,k } k are orthonormal bases of V C and W j , respectively. They are defined by dilatations and translations 4 of the so-called scale function φ and its associated wavelet function ψ. Functions φ and ψ verify the following two-scale relations:

φ(x) = √ 2 k∈Z h[k]φ(2x -k) ; ψ(x) = √ 2 k∈Z g[k]φ(2x -k) , (12)
where sequences h

[k] = φ(x), √ 2φ(2x -k) and g[k] = ψ(x), √ 2φ
(2x -k) are called conjugate mirror filters. Those filters play an important role in the fast implementation with filter banks of forward and inverse wavelet transform, i.e. projection on the wavelet basis and reconstruction, from [START_REF] Yuan | Simultaneous higher-order optical flow estimation and decomposition[END_REF] [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]. Finally, the representation of a signal projected onto the multiscale wavelet basis is given by the set of coefficients appearing in [START_REF] Yuan | Simultaneous higher-order optical flow estimation and decomposition[END_REF]: a C,k w, φ C,k L 2 and d j,k w, ψ j,k L 2 are approximation and detail coefficients, respectively. Those results are extended to the case of 2D signals, in order to obtain separable multiscale orthonormal bases of L 2 ([0, 1] 2 ).

Wavelet Data Term

In this work, the representation of the velocity field v is obtained by the wavelet decomposition [START_REF] Yuan | Simultaneous higher-order optical flow estimation and decomposition[END_REF] of each component. We denote by Θ 1 and Θ 2 the sets of coefficients respectively associated to v 1 and v

2 ; Θ = (Θ 1 , Θ 2 )
T is the set of all coefficients. Denoting the linear reconstruction operator by Φ for convenience, we may write

∀x ∈ Ω, v(x) = Φ(x)Θ . (13)
Here the constant coefficients vector Θ is the unknown of our optical flow estimation problem. Replacing v(x) by ( 13) in DFD data term (1), we obtain

J DFD (Θ) = 1 2 Ω [I 1 (x + Φ(x)Θ) -I 0 (x)] 2 dx ( 14
)
and the estimation problem becomes

v = Φ Θ ∈ V F +1 , where Θ = arg min Θ J DFD (Θ) . ( 15
)

Multiscale Estimation

Unknown coefficients are estimated sequentially from coarsest scale C to a chosen finest one L (with C ≤ L ≤ F ) using a gradient-descent algorithm. At each scale j, all coefficients from scales C to j are estimated. Coefficients previously estimated at coarser approximation spaces are used to initialize the gradient descent; this strategy enables the update of the latter coarser coefficients while estimating "new" details at current scale j. In other words, solution v is sequentially sought within higher resolution spaces:

V C ⊂ V C+1 ⊂ • • • ⊂ V L .
This way, the projection of the current solution v ∈ V j onto every coarser space V p with C ≤ p < j is constantly updated, contrary to the standard incremental approach (Sect. 2.2). The use of wavelet bases thus leads to a "natural" and well-defined multi-resolution framework. At each refinement level, minimization of functional J DFD is efficiently achieved with a gradient-based quasi-Newton algorithm (L-BFGS) [START_REF] Nocedal | Numerical Optimization[END_REF], to seek the optimum Θ. For any coefficient θ

i,p ∈ Θ i ⊂ Θ, ∂J DFD ∂θ i,p (Θ) = ∂I 1 ∂x i (• + Φ(•)Θ) [I 1 (• + Φ(•)Θ) -I 0 (•)] , Φ p L 2 ([0,1] 2 ) ( 16
)
where Φ p is the wavelet basis atom related to θ i,p . As a consequence, components of the spatial gradient of the data-term functional (14) are simply given by the coefficients of the wavelet decomposition of the two terms

∂I 1 ∂x i (x + Φ(x)Θ) [I 1 (x + Φ(x)Θ) -I 0 (x)] , i = 1, 2 ,
on the considered wavelet basis. It is easy to see that the proposed coarse-tofine estimation strategy enables to capture large displacements: at large scales, the decomposition of (3.3) is obtained by convolutions with the atoms of the wavelet basis having the largest support. Note that conversely to the algorithm proposed in [START_REF] Wu | Image registration using wavelet-based motion model[END_REF], the low-complexity of gradient computation via fast wavelet transform does not restrict motion estimation to large scales and/or images of small size.

Regularizations

Wavelet Properties

Wavelet-based regularizers which are described in the following are based upon wavelet properties such as polynomial reproduction, differentiation and interpolation. Those aspects are linked to the notion of vanishing moments (VM). A wavelet ψ(x) ∈ L 2 (R) has n VM if : R x ℓ ψ(x)dx = 0, for 0 ≤ ℓ < n .

Wavelets as polynomial approximations From (17), a wavelet with n VM is hence orthogonal to any polynomial of degree n -1. Consequently, piece-wise5 polynomials of degree n -1 belonging to V F +1 are exactly described in V F , since the atoms of the basis that belong to its orthogonal complement W F have vanishing coefficients.

Wavelets as Differential Operators Given a signal w ∈ C n , the behavior of its small scales coefficients resulting from an n-VM wavelet decomposition can be related to its n th derivative [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]:

lim j→∞ w(x), ψ j,k (x) 2 -j(n+ 1 2 ) ∝ ∂ n w(x) ∂x n . ( 18
)
This result can be extended to the case of 2D signals.

Wavelet-Based Interpolation A multiscale interpolation is the orthogonal projection of a signal w estimated at a given resolution 2 j onto the next finer approximation space V j+1 :

P Vj+1 w 2(p + 1 2 ) = +∞ k=-∞ w(2k)ϕ j+1 (p -k + 1/2) . ( 19
)
The interpolation function ϕ is defined as the autocorrelation of the scaling function: ϕ = φ ⋆ φ, where ⋆ and • denote respectively convolution and timereverse6 operators. It can be shown that ϕ interpolates exactly polynomials of order n if and only if the wavelet associated to scaling function φ has n + 1 VM [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]. This linear interpolation operator P Vj+1 is also implemented with filter banks using filter h i , where

h i [n] = h ⋆ ȟ [2n + 1]
and h is defined in (12).

Polynomial Approximation on a Truncated Basis

As seen in Sect. 2.3, a first way to overcome the under-constrained nature of the optical flow estimation problem consists in reducing the number of unknowns through a parametric formulation of the velocity field. Using the proposed wavelet formulation (15), this can be easily achieved by estimating the velocity field on a truncated wavelet basis. This means that the solution v belongs to a lower-resolution space V L ⊂ V F +1 and therefore is a piece-wise polynomial of order n -1 in V L+1 . Details coefficients associated to non-estimated small details scales (W j with L ≤ j ≤ F ) are thus not estimated, but set to zero.

v = Φ Θ ∈ V L , L < F + 1 , where Θ = arg min Θ J DFD (Θ) . ( 20
)
Since the basis truncation reduces the number of unknowns, it is theoretically possible to estimate detail coefficients up to penultimate scale F -1, i.e. v ∈ V F . Practically, it is impossible due to the aperture problem.

High-Order Regularization

It has been previously mentioned that smallest scale coefficients might be interpreted as the signal's n th derivative (Sect. 4.1), with n number of VM of the considered wavelet. The penalization of small scale coefficients' amplitude thus enables to control the amplitude of the derivative of the estimated signal. However, due to the dyadic structure of the discrete wavelet decomposition, only a "piecewise" control is possible. In order to control the derivative at junctions of those dyadic blocks, interpolated signal ṽ of the velocity field v on a shifted 2D grid is considered. Small scale coefficients { ΘF } and {Θ F } of both ṽ and v are penalized. "Interpolated coefficients" Θ are expressed as a linear combination of Θ through wavelet inverse and forward transformations (Φ, Φ -1 = Φ T , resp.)

and interpolation: Θ = Φ T • P V F +1 • Φ Θ. We finally get the regularization term

J reg (Θ) = 1 2 Θ F 2 + 1 2 ΘF 2 and ∇J reg (Θ) = Θ F + Φ T • P T V F +1 • Φ ΘF (21)
The gradient in (21) is a linear form which can be efficiently computed using the recursive filter banks presented in Sect. 3.1 and 4.1. The addition of the regularization term (21) therefore does not increase significantly the computational burden. Supplementing (15), the estimation problem becomes:

v = Φ Θ ∈ V F +1 , where Θ = arg min Θ J DFD (Θ) + µJ reg (Θ) . ( 22
)
5 Results

Daubechies wavelets have been chosen since they have a minimum support size for a given number of VM [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]. Daubechies wavelet with n VM will be referred to as D n hereafter. Wavelet transform is implemented with periodic boundary conditions. 1. Sample synthetic PIV image (1a) with below the vorticity of the underlying reference velocity field (1f). End-point error maps on velocity field estimations for a polynomial approximation (upper row ) and high-order regularization (lower row ) with Dn wavelets are presented, i.e. polynomial (resp. derivative) order of n -1 (resp. n), for D1 (1b, 1g), D2 (1c, 1h), D3 (1d, 1i) and D10 (1e, 1j).

Synthetic PIV Sequence

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) Fig.
The first data set used for evaluation is a synthetic sequence of Particle Imagery Velocimetry (PIV) images of size 256 × 256 pixels, representing small particles (of radius below 4 pixels) advected by a 2D periodic forced turbulent flow. The dynamic of the fluid flow is given by numerical simulation of 2D Navier-Stokes equations at Re = 3000, using the vorticity conservation equation and the Lagrangian equation for non-heavy particles transported by the flow (simulation details can be found in [START_REF] Heitz | Final report on the evaluation of the tasks of the workpackage 2, FLUID project deliverable 5.4[END_REF]). This simulated flow has a null-divergence by construction. An image of the PIV sequence is displayed in Fig. 1a together with its associated ground truth motion vorticity (1f). Estimated velocity field evaluation is based on the Root Mean Square end-point Error (RMSE) 7 .

When the true velocity field is decomposed on a Daubechies wavelet with a number of VM higher than 3, the velocity field reconstructed with the p = 6 coarser scales (out of 8) carries out more than 99,95% of the total kinetic energy. Those 6 scales are represented with only 12.5% of atoms from the full wavelet basis. Moreover, when n is chosen high enough (> 7), the reconstruction error stabilizes around 0, 013.

Motion Estimation

From the previous analysis of the true velocity fields, it seems that 6 detail scales out of 8 should give an accurate representation of the motion in terms of kinetic energy, as long as the chosen number of VM is high enough. Figure 2a (black curve) shows RMS errors computed on an estimated velocity field with truncated Daubechies wavelet bases (Sect. 4.2) having different VM n, i.e. with a polynomial approximation of order n -1. As expected, RMSE converges rapidly towards a median value of 0.0613 when n increases. Figure 1 (upper row) shows corresponding end-point error maps for motion estimated with wavelets bases D 1 , D 2 , D 3 and D 10 . Although errors effectively lower when higher VM wavelets are employed, artifacts due to high-amplitude errors on 6 th scale coefficients (small white "dots") and on coarse coefficients (white straight "lines") remain clearly visible. With the proposed high-order regularization scheme (Sect. 4.3), all scales are estimated, which should highly improve results with n ≤ 3 VM, i.e. for penalization of derivatives of order lower than 3. This is confirmed on Fig. 2a (red dashed curve), with a reduction of 35% and 30% of the RMS obtained with D 1 and D 2 wavelets bases, respectively, whereas the diminution observed using D 10 wavelet basis is of 10% at best. At the same time, derivative penalization eliminates most of the artifacts observed on estimates with truncated bases, which is displayed on the lower row of Fig. 1. Note also that there are less differences between estimations with different VM, in comparison to the previous case. Fig. 2. Left: Comparison of RMS errors on velocity fields estimated from a pair of the synthetic PIV sequence, with the proposed methods and Daubechies wavelets with 1 to 10 VM, i.e. for polynomial approximation (resp. derivative regularization) of order 0 to 9 (resp. 1 to 10). RMSE obtained with the polynomial approximation (6 scales out of 8, dashed line) and with derivatives regularization (best case, solid line). Right: Comparison of RMS errors on a sequence of velocity fields estimated with proposed regularization (thick solid ) and with state-of-the-art methods: correlation (thin solid ), first order regularization [START_REF] Horn | Determining optical flow[END_REF] (long-dashed ), div-curl regularization [START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF] (dashed ), selfsimilar regularization [START_REF] Heas | Bayesian selection of scaling laws for motion modeling in images[END_REF] (dotted ).

Comparison with State-of-the-Art Estimators

Real PIV Sequence

This data set consists in 128 PIV pictures of a transversal view of a planar concomitant jet flow, of size 1024 × 1024 pixels. The flow has a "top-hat" velocity profile and is poorly turbulent, but shows two high-shear regions featuring development of Kelvin-Helmholtz instabilities. Motion is estimated with proposed wavelet-based estimator (22), using the following settings: 2 VM and derivatives penalization with factor µ = 10 7 . Figure 3 presents a PIV image of the sequence and streamlines of an estimated velocity field along with two consecutive vorticity maps. A qualitative evaluation of the presented motion field shows a remarkably good agreement with the physics of concomitant jets. A very good temporal coherence is also observed, although no prior dynamic model is considered (i.e. successive pairs of images are processed independently).

Conclusion

An optical flow estimation algorithm dedicated to continuous motion has been introduced. The choice of the wavelet formalism sets-up a well-defined multiresolution framework that avoids most drawbacks of such usual approaches. Being associated to a gradient-based quasi-Newton optimization method, its low complexity makes possible the estimation of the full range of scales composing the motion. Moreover, high numbers of vanishing moments enable to truncate the wavelet basis without increasing the error of the polynomial reconstruction, thus significantly reducing the number of unknowns and the problem complexity. A high-order regularization scheme, involving small scale coefficients penalization, highly enhances estimation results and generally helps reducing errors by removing noise of the solution, as emphasized by experiments on a synthetic PIV sequence. Application to a real PIV sequence shows the capability of the estimation method to reconstruct accurately vortices of large amplitude.

Figure 2

 2 is a comparison of RMS errors obtained on the synthetic PIV sequence with the proposed highorder regularizer and various state-of-the-art estimators, after a null-divergence projection. Our wavelet-based estimator clearly outperforms other methods.

Fig. 3 .

 3 Fig. 3. Sample estimated motion fields from 2D planar jet PIV dataset: detail of input PIV image (3a), streamlines (3b) and vorticity (3c). Figure 3d is the vorticity field corresponding to motion estimated at the next time step. Three different areas are visible: at the output of the jet (top of the field), shear regions begin to oscillate slowly. The middle region clearly shows the development of vortices characteristic of the Kelvin-Helmholtz instability. Finally, in the lower part of the field, structure of vortices collapse due to their tri-dimensionalization.

In this paper, we shall use the following convention: indices j ≥ 0 represent the resolution

j -contrary to[START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]. Corresponding scale is 2 -j . 2 A Gaussian kernel of variance proportional to 2 j is commonly used.

See[START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF] for a complete presentation of wavelet bases.

Written in a general form f j,k (x) = 2 j/2 f (2 j x -k).

On the support of {ψ F,k }.

More explicitly, f : t → f (t) = f (-t).

Ground-truth velocity fields being given on a shifted grid (by 1/2 pixel) by the numerical simulation, they have been interpolated in order to compute accurate RMSE on the pixel grid.

Acknowledgments

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant MSDAG (ANR-08-SYSC-014) "Multiscale Data Assimilation for Geophysics".


Pierre Derian
email: pierre.derian@inria.fr

Patrick Heas
email: patrick.heas@inria.fr

Etienne Memin
email: etienne.memin@inria.fr

WAVELETS TO RECONSTRUCT TURBULENCE MULTIFRACTALS FROM EXPERIMENTAL IMAGE SEQUENCES

In the context of turbulent fluid motion measurement from image sequences, we propose in this paper to reverse the traditional point of view of wavelets perceived as an analyzing tool: wavelets and their properties are now considered as prior regularization models for the motion estimation problem, in order to exhibit some well-known turbulence regularities and multifractal behaviors on the reconstructed motion field.

INTRODUCTION AND MOTIVATION ON WAVELETS AND TURBULENT MOTIONS

The fractal nature of turbulent motion fields may be derived from A. Kolmogorov statistical theory. More recently, Frisch (1995) showed how multifractal models are adapted to the description of intermittency phenomena in turbulence. These notions may be characterized by the Lipschitz regularity α or the singularity spectrum D(α) of a signal, resp. In this context, wavelets are known to constitute a simple yet efficient tool for regularity measurement. Given a signal, the asymptotic decay across scales of the amplitude of its wavelet coefficients can be linked to its Lipschitz regularity, under simple constraints on the analyzing wavelet. In the following, we consider a wavelet with n vanishing moments 1 and n fast-decay derivatives. Then, considering wavelet coefficients 2 {d j,p } j,p of given function f ∈ L 2 (R) uniformly Lipschitz α ≤ n over the support [0, 1], there exist A such that (Mallat, 2008):

|d j,p | = | f , ψ j,p L 2 | ∼ A(2 j ) α+1/2
(1) * Authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant MSDAG (ANR-08-SYSC-014) "Multiscale Data Assimilation for Geophysics".

1 The notion of vanishing moments (VM) is simply orthogonality to polynomials up to a certain order.

2 See Sect. 2.2 for an introduction to wavelet formalism.

Conversely, a bounded signal whose coefficients verify (1) for α / ∈ N is uniformly Lipschitz α over ]0, 1[. Similarly, the singularity spectrum may be estimated from the computation of the partition functions of the coefficients and an inverse Legendre transform (Mallat, 2008). It has been used to exhibit the multifractal nature of turbulent velocity fields [START_REF] Bacry | Singularity spectrum of fractal signals from wavelet analysis: Exact results[END_REF].

Fractality of signals is also known to reveal a power-law behavior in the power spectrum. Applied to turbulent motion, it rises up the famous -5/3 decay predicted by Kolmogorov theory. Using continuous wavelet transform (CWT), Perrier et al. (1995) define a (continuous) wavelet power spectrum (CWPS) from the energy conservation property and the multiscale nature of the CWT. At a given scale, this CWPS is proportionnal to the L 2 norm of the wavelet coefficients. Moreover, it is related to the usual Fourier power spectrum (FPS) for a finite-energy signal. In particular, considering a power-law behavior at small scales (wavenumber k → ∞) of the Fourier spectrum E(k):

E(k) ∼ k -γ
then a sufficient condition for the CWPS to exhibit the same behavior as the FPS is to use a wavelet with n > γ-1 2 VM.

Aside from their use in turbulence analysis, wavelets have been employed to design motion estimators from fluid flow image sequences. Indeed, methodological advances in correlation-based or model-based measurement techniques have enhanced the accuracy of the computer-vision approach to fluid motion measurement (particle image velocimetry [START_REF] Adrian | Particle imaging techniques for experimental fluid mechanics[END_REF], scalar imagery, meteorological or oceanographical satellite imagery (Heitz et al., 2010), . . . ). In order to measure motion from an image sequence, any motion estimation method needs to impose either explicitly (e.g. optical flow) or implicitly (e.g. correlation) smoothness assumption on the solution through smoothness functions (so-called regularizers) or implicit parametric spatial constraints (constant, affine, etc.). Contrary to correlation-based methods, model-based techniques known as optical flow approaches present the advantage of offering a wide range of observation models (Liu & Shen, 2008) and regularizers, which may be chosen accordingly to the specificities of the observed flow. Authors in Heas et al. (2009) have designed an efficient regularization model based on statistical properties of turbulence, using the velocity increment function. More recently, a wavelet-based observation model, taking advantage of the intrinsic multiscale nature of wavelet bases, was proposed in [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF]. In the continuation of those works, the idea developed in this paper is to design wavelet-based regularizers to reconstruct the motion field with some desired physical properties such as fractality. Since [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF] wavelet-based observation model relates directly wavelet coefficients of the sought motion field to image intensity variations, fractal properties on the solution can be imposed simply by constraining its wavelet coefficients decay trough scales. The wavelet framework therefore enables a much simpler and more efficient implementation of fractal regularization, and is in the same time well-adapted to optical flow estimation.

WAVELET-BASED OPTICAL FLOW

When addressing fluid motion estimation, one has to deal with continuous complex motion, involving a wide range of scales and velocity magnitudes. Correlation-based methods have been used successfully, especially in the context of experiments featuring flows seeded with particles (PIV). Those methods use a single data model, namely Digital Image Correlations (DIC), which is notably adapted to PIV. However, they have the drawbacks of producing sparse velocity fields, while their regularization scheme (i.e. blockwise constant fields) imposes an erroneous regularity to the reconstructed motion field. On the opposite, optical flow methods reconstruct dense motion field (i.e. one velocity vector per pixel) and offer a wide range of data models and regularizers, which can be chosen in adequacy with the specific nature of the problem. Fluid-motiondedicated regularization schemes preserving rotational and/or divergence quantities, or built upon statistical properties of the turbulence, have been proposed. They have similar accuracy to correlations methods on PIV and outperform them on scalar transport images.

Optical Flow Problem

Optical flow aims at estimating the apparent motion field v(x) = (v 1 (x), v 2 (x)) T transforming an image I 0 = I(t) of the sequence into the next one I 1 = I(t + 1), ∀x ∈ Ω, where Ω ⊂ R 2 denotes image domain. This problem is classically solved through the minimization of a functional J(v). This functional J includes a data model linking observed images to the sought motion fields. The simplest assumption used is the luminance consistency along a point trajectory:

dI(x(t),t) dt = 0,
which constitutes an exact physical model in the case of non-diffusive transport by bi-dimensional turbulent flows (Liu & Shen, 2008). The data model used hereafter, so-called the Displaced Frame Difference (DFD), is obtained by time-integration along the motion trajectories of the previous assumption:

∀x ∈ Ω, I 0 (x) -I 1 (x + v(x)) = 0 (2)
The velocity field estimate v is finally obtained by minimizing quadratic discrepancies to the data model ( 2):

     v = arg min v J DFD (v) J DFD (v) = 1 2 Ω [I 0 (x) -I 1 (x + v(x))] 2 dx.
(

) 3
This problem is ill-posed, as there are twice as much unknowns as equations. Moreover, this data model does not provide any information when image spatial or temporal gradients vanish. Regularization schemes are needed to circumvent those limitations.

Projection on Orthogonal Wavelet Bases

We consider the projection of each scalar component v 1 , v 2 of the velocity field v onto multi-resolution approximation spaces exhibited by the wavelet formalism. Let us introduce briefly this context for real 1D scalar signals -see Mallat ( 2008) for more details. We consider a multi-resolution approximation of L 2 (R) as a sequence {V j } j∈Z of closed subspaces, so-called approximation spaces, notably verifying:

V j+1 ⊂ V j ; lim j→+∞ V j = +∞ j=-∞ V j = {0} lim j→-∞ V j = Closure   +∞ j=-∞ V j   = L 2 (R)
Since approximation spaces are sequentially included within each other, they can be decomposed as

V j-1 = V j ⊕ W j .
The W j are the orthogonal complements of approximation spaces, they are called detail spaces.

Practically, scalar 1D signals being finite, they belong to a given approximation space according to their resolution, i.e. number of samples. Writing

N j = 2 -j ∀ j ≤ 0 hereafter, let w be a 1D signal of N F -1 sam- ples. Then, w ∈ V F -1 = V C ⊕ W C ⊕ W C-1 ⊕ • • • ⊕ W F ⊂ L 2 ([0, 1]), where F ≤ C ≤ 0.
The projection of w on this multiscale basis writes:

w(x) = N C -1 k=0 w, φ C,k L 2 φ C,k (x) + F j=C Nj -1 k=0 w, ψ j,k L 2 ψ j,k (k) (4)
Here, {φ C,k } k and {ψ j,k } k are orthonormal bases of V C and W j , respectively. They are defined by dilatations and translations 3 of the so-called scaling function φ and its associated wavelet function ψ. Finally, the representation of a signal projected onto the multiscale orthonormal wavelet basis is given by the set of coefficients appearing in (4): a C,k w, φ C,k L 2 and d j,k w, ψ j,k L 2 are approximation and detail coefficients, respectively. Those results are extended to the case of 2D signals, in order to obtain separable multiscale orthonormal bases of L 2 ([0, 1] 2 ).

Wavelet-Based Problem Formulation

In this work, we consider the wavelet expansion of each component v i of the sought motion field v; Θ i is the set of coefficients representing v i on the considered wavelet basis. Writing Θ = (Θ 1 , Θ 2 ) T the superset of all coefficients, we have:

v(x) = v 1 (x) v 2 (x) = Φ(x)Θ 1 Φ(x)Θ 2 Φ(x)Θ (5)
Supplementing (3), the optical flow problem becomes

       v = Φ arg min Θ J DFD (Θ) J DFD (Θ) = 1 2 Ω [I 0 (x) -I 1 (x + Φ(x)Θ)] 2 dx (6)
This wavelet-based motion estimator has proven to be efficient in the context of fluid motion estimation [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF], using simple regularizers (basis truncation) or more elaborate high-order regularizers.

RECONSTRUCTION OF TURBULENCE MULTIFRACTALS 3.1 Lipschitz Regularity

Singularity spectrum D(α), with α the Lipschitz regularity of a turbulent process, is related to the decay across scales of the p th order partition function (PF) of the wavelet coefficients, which is given by the p-norm of coefficients vector at a given scale. Because of the correspondence with wavelet spectrum, we consider the 2 nd order PF defined as:

Z(2 j ) = Nj -1 p=0 d j,p 2 (7)
For a 2D multifractal signal, Legendre transform relates the scaling exponent of the 2 nd order PF (power-law behavior) to its singularity spectrum -see Mallat (2008):

ζ 2 = lim j→-∞ log 2 Z(2 j ) log 2 2 j = min α∈[αmin,αmax]
(2(α + 1) -D(α))

(8)
3 Written in a general form ψ j,p (x) = 2 -j/2 ψ(2 -j xp).

In the case of a monofractal signal, ζ 2 = 2α.

Discrete Wavelet Spectrum

Similarly to relations involving continuous Fourier and wavelet spectra (Sect. 1), it is also possible to link discrete wavelet power spectrum (DWPS) to the Fourier one. For 1D signal f ∈ L 2 ([0, 1]), we define the Fourier (9) and mean wavelet (10) spectra:

E(k) = 1 N F f [k] 2 , ∀k ∈ [0 • • • N F -1] (9) Ẽ(N j ) = 1 2 Nj -1 p=0 d j,p 2 , ∀j ∈ [-(F -1) • • • 0] (10)
where { f [k]} k and {d j,p } j,p are Fourier and wavelet coefficients, resp. Derivation of the following relation between the two spectra can be found in App. A.

Ẽ(N j ) = 1 (N F ) 2 k E(k) Nj -1 p=0 ψ * [(k + p)/N j ] 2 (11)
Considering a power-law behavior k -γ at small scales of FPS, it is shown in App. B that, under simple conditions on the VM of chosen wavelet, we obtain a bound on DWPS decay :

E(k) ∼ k -γ and ψ with m > γ + 1 2 VM ⇒ Ẽ(k) ∼ k -(γ-1) (12)
In the 2D case, the bias is of 2:

E(k) ∼ k -γ ⇒ Ẽ(k) ∼ k -(γ-2) (13)
A notable difference with the CWPS relation (Sect. 1) is the apparition, in the last sum of (11), of wavelet spectrum foldings. A consequence in ( 12) is the bias between decays of FPS and DWPS.

Regularization Scheme

Using properties (8) and ( 13), we now define adequate smoothing functions for the optical flow problem exhibiting particular fractal behavior or FPS decay. Since expressions of the 2 nd order PF (7) and the DWPS (10) are very similar, the following derivation of the regularization scheme is presented for the spectrum only. In 2D, from the hypothesis of a power-law behavior at small scales:

Ẽ(N j ) = A(2 -j ) -(γ-2) , ∀ F ≤ j ≤ L ≪ 0 We obtain d log 2 Ẽ dj (N j ) = -(γ -2)
, which is approached by a finite difference. At scale 2 j , the powerlaw behavior finally corresponds to:

-1 γ -2 log 2 Ẽ(N j ) -log 2 Ẽ(N j+1 ) -1 = 0 (14)
Then the constrained estimation problem writes:

Θ = arg min Θ J DFD (Θ) s.t. c j (Θ) = (14) ∀F ≤ j ≤ L (15)
In order to exhibit a power-law of slope -γ to the FPS, one has to pick constraints with -(γ -2) in 2D ( 13). Considering a Lipschitz regularity α for the reconstructed motion, since in 2D we have ζ 2 = γ -2, constraints ( 14) write for a monofractal with γ -2 = 2α whereas for a multifractal γ -2 = minα (2(α + 1) -D(α)) (Sect. 3.1). Note that in the latter case, we known that the minimizer is lower than arg maxα D(α) -see (Mallat, 2008). Therefore, constraints apply on a point located on the left part of singularity spectrum.

Optimization

The proposed wavelet-based estimation is optimized using fast filter-banks [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF]. Daubechies wavelets are chosen since they have the shortest support for a given number of VM, hence minimizing computations. The algorithm processes as follows:

1. coarser (unconstrained) scales are estimated first, scale by scale, using L-BFGS (Nocedal & Wright, 1999) to achieve minimization of DFD functional; 2. to exhibit the desired property (regularity or spectrum decay), involved fine scales are also estimated, enabling constraints. Minimization of the DFD functional under constraints (15) is achieved using augmented Lagrangian method (App. C) which is appropriate here since it does not require the initialization point to respect constraints.

RESULTS

Synthetic PIV Dataset

Following results were obtained from a synthetic PIV image sequence generated from a 2D turbulence simulation at Re = 3000, using vorticity conservation equation and Lagrangian equation for non-heavy particle transported by the flow -details can be found in [START_REF] Carlier | Report 1 on production and diffusion of fluid mechanics images and data[END_REF]. Ground truth velocity fields were thus available for evaluating the reconstruction error reduction brought by the method. Figure 1 presents an example of synthetic PIV picture, along with vorticity of the underlying simulated and estimated velocity fields. Figure 2 shows the evolution of the DWPS decay slope and dominant Lipschitz regularity, as functions of the wavelet number of VM. Slope converges towards the FPS slope of -5.05 when the analyzing wavelet has enough VM. Similarly, the estimated dominant Lipschitz regularity α 0 converges towards a median value of 1.85. Figures 3 and4 compare power and singularity spectra, resp., of one ground-truth velocity field, of an estimate with proposed power-law regularization and of another estimate using 2 nd order regularization as introduced in [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF]. Note that because of the bias (13), all DWPS graphs are compensated by k -2 to facilitate comparisons with FPS. DWPS at scale range 2 -5 to 2 -7 has been constrained to decay according to a power-law of (compensated) slope -ζ = -5.05, in adequacy with the slope observed on ground-truth spectra. The reconstructed FPS oscillate at finer scales, inducing a slightly lower slope decay for sufficient VM -as expected from the bound given in App. B. Root Mean Square (RMS) error on estimated motion obtained with power-law regularization is 0.086, which is to be compared with the RMS of 0.083 obtained with 2 nd order regularization, keeping in mind that the objective of this work is the reconstruction of statistical properties rather than RMS reduction.

Scalar Diffusion

The second dataset consists in pictures of a passive scalar diffused by an electromagnetically-forced 2D turbulent motion field, investigated in Jullien et al. (2000). The flow features a direct enstrophy cascade, with a power-law energy spectrum E(k) = k -3 . Figure 5 shows an input image, along with the vorticity field obtained by motion estimation with a prior on the spectrum slope. Figure 6 superimposes FPS and DWPS of the estimated motion file, where PF of scales 2 -4 to 2 -8 are constrained to decay according to a power-law of slope -ζ = -3.

CONCLUSION

An extension of wavelet-based optical flow algorithm dedicated to turbulent fluid motion estimation has been presented. Using wavelet properties, it enables the reconstruction of turbulence multifractals or power-law spectrum by controlling decay across scales of coefficients amplitude. Results on synthetic PIV as well as on experimental scalar diffusion images show enhancement brought by the proposed methodology, in terms of spectral reconstruction, in comparison to previously introduced high-order regularizers.

tional Methods in Computer Vision (SSVM2011

A Discrete Wavelet Spectrum Derivation

Let f be a signal of Ω F , i.e. with N F = 2 -F samples. We consider its Fourier and Wavelet transforms:

f [n] = 1 N F N F -1 k=0 f [k] exp ikn = -F <j≤0 Nj -1 p=0 d j,p ψ j,p [n]
For simplicity we here considered f with a zero mean, i.e. its coarser approximation in the wavelet transform vanishes. The wavelet coefficients can be linked to the signal and wavelet Fourier coefficients:

d j,p = 1 √ N F f ⋆ ψj [2 j-F p] = 1 N F √ N F N F -1 k=0 f [k] ψj [k] exp i2πkp/Nj (16)
Then it is possible to link DWPS (10) to FPS (9). Wavelet coefficients are replaced by ( 16), then Fubini theorem allows to swap sums order.

Ẽ(N j ) = 1 2 Nj -1 p=0 |d j,p | 2 = k1,k2 f [k 1 ] f * [k 2 ] ψj [k 1 ] ψj * [k 2 ] N j δ[k 1 -k 2 ] 2(N F √ N F ) 2 = N j (N F ) 2 k E(k) ψj [k] 2 (17)
Spectrum of the scaled wavelet ψj can be expressed in terms of mother wavelet spectrum. Since ψ j [n] 1 √ 2 j ψ[n/2 j ], j ≤ 0, using the following leads to (11). ψj

[k] = 1 N j Nj -1 p=0 ψ * [(k + p)/N j ] (18)

B Power-Law Behavior

Similarly to what is done in Perrier et al. (1995) with continuous transforms, we consider a power-law behavior of the Fourier spectrum at large wavenumbers:

E(k) = k -γ , ∀ k > k 0 ( 19
)
and a wavelet with m VM, which implies the existence of a continuous and bounded function φ such that:

ψ[k] = k m φ[k] and φ[0] = 1 (20)
Using ( 19), the wavelet spectrum ( 17) is split in two terms at k = k 0 , and ( 20) is used to replace E(k) where k > k 0 . In the limit of infinite signals (F → -∞), using Cauchy-Schwarz inequality, then Fubini theorem to swap sums and finally bounding one of the sums by its analogous continuous integral, we obtain:

k ≤ k 0 : Ẽ(N j ) ≤ (N j ) -2m Nj -1 p=0 k0+p k=p k 2m E(k -p) φ( k N j ) 2 ∼ (N j ) -2m+1 k > k 0 : Ẽ(N j ) ≤ (N j ) -γ

C Augmented Lagrangian

Augmented Lagrangian (AL) is related to quadratic penalty, but avoid ill-conditioning by introducing explicit Lagrange multipliers (Nocedal & Wright, 1999). From (15), AL writes at iteration k:

L A (Θ k , λ k , µ k ) =J DFD (Θ k ) (21) - F j=L λ k j c j (Θ k ) + F j=L 1 2µ k j (c j (Θ)) 2
Here {λ k j } j and {µ k j } j are respectively Lagrange multipliers and penalization factors associated to constrained scales F ≤ j ≤ L. At iteration k, minimizer Θk of functional ( 21) is obtained using L-BFGS method. After convergence, multipliers are updated according to

λ k+1 j = λ k j - cj ( Θk )
µ k , whereas penalization factors decay according to µ k+1 j ∈ [0; µ k j ]; here we used µ k+1 j = 1 3 µ k j .

Figure 1 :

 1 Figure 1: Synthetic PIV picture (left) used for motion estimation, along with vorticity of the underlying ground-truth motion field (middle) and of the reconstructed motion with proposed spectrum prior (right).

Figure 2 :

 2 Figure 2: Analyze of ground-truth velocity field: evolution of the WPS slope (continuous blue) and dominant Lipschitz regularity α 0 (continuous red), as functions of the wavelet VM. WPS (compensated) slope of groundtruth converges towards the FPS slope -γ = -5.05, α 0 converges towards a median value of 1.85. Results from motion estimation with different VM and the proposed power-law/Lipschitz regularizations, using ground-truth values as priors: Lipschitz regularity rapidly converge towards the expected value (dashed red), whereas the convergence slows down with the FPS slope (dashed blue).

Figure 3 :

 3 Figure 3: Comparison of FPS (continuous) and DWPS (dashed lines) obtained with a 9 VM wavelet. Graphs have been shifted vertically by 3 decades relatively so as to enhance visualization. Red graphs are groundtruth spectra exhibiting a power-law behavior of slope -α = -5.05; FPS and DWPS are in good agreement.Blue spectra are obtained by proposed estimator with spectrum slope regularization. Pink graphs are obtained with 2 nd order regularization as presented in[START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF]. Both WPS and FPS are better reconstructed by the power-law regularization (blue) compared to the 2 nd order one (pink), although FPS slightly oscillate at finest scales.

Figure 4 :

 4 Figure 4: Comparison of singularity spectra. Red is obtained from ground-truth field, analyzed with a 7 VM wavelet. Blue results from motion estimation using ground-truth prior. Pink gives for comparison the spectrum of the estimate with 2 nd order regularization.

Figure 5 :

 5 Figure 5: Passive tracer diffusion picture (left), along with vorticity of estimated motion by optical flow with spectrum slope regularization (right).

Figure 6 :

 6 Figure6: Comparison of FPS (continuous line) and DWPS (dashed line) of the estimated motion field obtained with a 9 VM wavelet, with a power-law behavior of slope -α = -3 (blue). The vertical shift of 3 decades between blue and pink graphs has been introduced for visualization purpose. DWPS (compensated) slope is exactly -3, while FPS slope slightly deviates by ∼ 3% In pink are given, for comparison, spectra corresponding to motion estimated using 2 nd order regularization, showing clearly an erroneous behavior at medium and fine scales.

y

  2m |φ(y)| 2 dy ∼ (N j ) -γ+1 if γ < 2m + 1Thus, the second term dominates when N j → ∞ if m > γ 2 . Finally, a signal whose FPS is a power-law of slope -γ will give a DWPS of slope steeper or equal to -(γ -1), if previous condition on the wavelet VM is verified. Extension to the case of 2D signals gives a DWPS decay bound of -(γ -2) for a -γ FPS slope.

  ).

	Ein-Gedi, Israel.
	Frisch, U. 1995 Turbulence : the legacy of A.N. Kol-
	mogorov. Cambridge university press.
	Heas, P., Memin, E., Heitz, D. & Mininni, P.D. 2009
	Bayesian selection of scaling laws for motion model-
	ing in images. In International Conference on Com-
	puter Vision (ICCV'09). Kyoto, Japan.
	Heitz, D., Mémin, E. & Schnörr, C. 2010 Variational
	fluid flow measurements from image sequences: syn-
	opsis and perspectives. Experiments in fluids 48 (3),
	369-393.
	Jullien, M.C., Castiglione, P. & Tabeling, P. 2000 Ex-
	perimental Observation of Batchelor Dispersion of
	Passive Tracers. Physical Review Letters 85 (17),
	3636-3639.
	Liu, T. & Shen, L. 2008 Fluid flow and optical flow.
	Journal of Fluid Mechanics 614, 253.
	Mallat, S. 2008 A Wavelet Tour of Signal Processing:
	The Sparse Way. Academic Press.
	Nocedal, Jorge & Wright, Stephen J. 1999 Numerical
	Optimization. New York, NY: Springer-Verlag.
	Perrier, V., Philipovitch, T. & Basdevant, C. 1995
	Wavelet spectra compared to Fourier spectra. Jour-
	nal of Mathematical Physics 36, 1506.


Fawaz Paraiso

Gabriel Hermosillo

Romain Rouvoy

Philippe Merle

Lionel Seinturier

Fawaz Paraïso

Distributed Complex Event Processing Engine

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

Introduction Experimental Deployment

 DiCEPE (Distributed Complex Event Processing) is a platform that focuses on the integration of CEP engines in distributed systems, and which is capable of using various communication protocols.

 It was built using a component-based approach, and inherits the flexibility and adaptation facilities provided by FraSCAti.

 The DiCEPE platform is used to federate complex event processing.

 We developed a solution for a crisis management scenario using DiCEPE .

 This scenario was deployed in a cloud environment using CloudBees, a public Platform as a Service (PaaS) provider.

 Each actor of the scenario runs its own instance of DiCEPE.

Platform Overview

Scalability

 To evaluate the scalability, we used two different datasets of event generators: one with 10,000 and a second one with 15,000 which generated around 500,000 and 750,000 events respectively.

 As shown in the graphs, the processing time for each event remained stable and very low during both benchmarks (around a 10th of a millisecond), despite the fact that the average number of simultaneous sessions had a significant increase of about 50% (from 89 with the first dataset, to 135 with the second).

Integration

 The DiCEPE platform facilitates the integration of complex event processing engines.



  The architecture of DiCEPE is composed of four parts:  Engine: This component acts as the engine instance, by which Statement components, events, and outputs (Listener component) are registered.



  Statement: This component is used for quering the inbound event streams. It is registered within the Engine component, which at the same time is connected to one or many Statement components.  Listener: This component generates a new complex event when an action is detected.  Context: This component collects information of the executing environment, like the number of statement rules deployed in the engine at run-time.



  We integrated DiCEPE with the Esper and Etalis CEP engines. Etalis Conclusion  DiCEPE is a platform that offers interoperability for Distributed Complex Event Processing engines, via federation. It focuses on providing a very flexible component architecture, which supports the interaction of different complex event processing engines simultaneously, while enabling communication among them with a distributed execution and deployment system.

Acknowledgments. This work is partially funded by the French Ministry of Higher Education and Research, Nord-Pas de Calais Regional Council and FEDER through the Contrat de Projets Etat Region Campus Intelligence Ambiante (CPER-CIA) 2007-2013, and the ANR (French National Research Agency) ARPEGE SocEDA project.


Anas El

Husseini • Abdallah M'hamed

Bachar El

Mounir Mokhtari

Trust-Based Authentication Scheme with User Rating for Low-Resource Devices in Smart Environments

Keywords: Smart environments, privacy preservation, trust evaluation, authentication, low-resource devices

In smart environments, pervasive computing contributes in improving daily life activities for dependent people by providing personalized services. Nevertheless, those environments do not guarantee a satisfactory level for protecting the user privacy and ensuring the trust between communicating entities. In this paper, we propose a trust evaluation model based on user past and present behavior. This model is associated to a lightweight authentication key agreement protocol (EC-SAKA). The aim is to enable the communicating entities to establish a level of trust and then succeed in a mutual authentication using a scheme suitable for low-resource devices in smart environments. An innovation in our trust model is that it uses an accurate approach to calculate trust in different situations, and includes a human-based feature for trust feedback, which is user rating. Finally, we tested and implemented our scheme on Android mobile phones in a smart environment dedicated for handicapped people.

Introduction

Trust, Anonymity and Privacy preservation are known as major factors to the acceptance and the success of Anas EL HUSSEINI • Abdallah M'HAMED • Mounir MOKHTARI Telecom SudParis, Handicom Lab, Evry, France Anas EL HUSSEINI • Bachar EL HASSAN Lebanese University, LaSTRe Lab, Azm Center for Scientific Research, Tripoli, Lebanon pervasive computing systems. Service Platforms tend to collect and manage a large amount of personal information about individuals in order to authenticate users and/or provide personalized services, therefore threating privacy and causing a conflict between service providers and personal information owners. People dislike automatic spread of personal and identifiable data, especially when it is transferred to other parties beyond control [START_REF] Langheinrich | Privacy by Design-Principles of Privacy-Aware Ubiquitous Systems[END_REF] [START_REF] Stajano | Security for Ubiquitous Computing[END_REF]. The trade-off between privacy and collecting private data for authentication poses nowadays a great challenge to security designers in Smart Environments.

For disabled and aging users, smart environments are deployed to facilitate every day life activities and adapt to their needs. A ubiquitous computing environment entails an extensive and complex computer architecture deployment, sophisticated data control, and a judicious external interface facilitating user interaction with the system [START_REF] Modahl | UbiqStack: a taxonomy for a ubiquitous computing software stack[END_REF]. The characteristics of ubiquitous systems amplify the concern of security problems, by promoting spontaneous interactions between diverse heterogeneous entities [START_REF] English | Towards self-protecting ubiquitous systems: monitoring trust-based interactions[END_REF]. However, within those environments, the security and confidentiality of sensitive data has to be guaranteed. Beside that, the anonymity and privacy of the users should be protected, therefore personal information concerning users like names, addresses, financial data, and medical profile should not be allowed to flow freely without protection. Therefore, one of the foundations of the security of users in smart environments is to protect individuals privacy.

Several researchers [START_REF] Langheinrich | Privacy by Design-Principles of Privacy-Aware Ubiquitous Systems[END_REF] [START_REF] Stajano | Security for Ubiquitous Computing[END_REF] have admitted that smart environments are vulnerable to many security and privacy threats, and that securing pervasive computing present critical challenges at many levels [START_REF] Shuxin | A Trust Model for Pervasive Computing Environments[END_REF][6] [START_REF] Heng | Protection of Privacy in Pervasive Computing Environments[END_REF]. Below, some of the challenges addressed in [START_REF] Campbell | Towards Security and Privacy for Pervasive Computing[END_REF] are outlined:

1. Privacy Issues: Sensors and actuators distributed in space expose a great danger to user privacy since the information collected can be disclosed to intruders, malicious insiders and tracking systems.

User Interaction Issues:

The access control mechanisms in pervasive environments should allow users to interact easily with devices while ensuring an appropriate authentication.

3. Security Policies: Smart environments should have a convenient method to define and manage security policies with dynamicity and flexibility, with respect to the behavior of entities in their systems. [START_REF] English | Towards self-protecting ubiquitous systems: monitoring trust-based interactions[END_REF]. Two new security challenges, introduced by [START_REF] Tentori | Supporting Quality of Privacy (QOP) in Pervasive Computing[END_REF], are to be added: Quality of Privacy (QoP) and Trustworthy Authentication. The second is combination of trust and authentication in one scheme.

In this paper, we design a trust model and implement it associated to an authentication protocol. The aim is to introduce an intermediate phase that evaluates the trustworthiness of communicating entities before the phase of service provision. Our scheme will preserve and protect user privacy since it will use non-sensitive information in the evaluation process of the trust. What makes trust evaluation models really needed in smart environments is that they create a secure yet more flexible environments that what security policies alone cannot do. That is because the strictness of static policies may limit the freedom of normal users while it is potentially vulnerable to malicious users; on the other hand, trust systems are very adaptive to the user needs, actions and behaviors. In our enhanced trust model, we also integrated two new ideas. The first is Questions of Trust that aims to enhance the model immunity against malicious manipulation of trust values, as well as dealing more efficiently with situations where there are only new neighbors. The other innovation is User Rating that allows users to give feedback about the quality of service they received, which will accordingly affect how much they trust the service provider in the future.

The remainder of this paper is organized as follows.

In section 2, we present several researches in literature concerning trust and repuation models. In section 3, we describe our proposed trust model featuring its advantages, pointing out the two new features: Questions of Trust and User Rating. Section 4 describes the features of the Elliptic Curve-based Simple Authentica-tion Key Agreement (EC-SAKA) protocol, designed for low-resource mobile devices. In section 5, we show some simulations concerning the trust model metrics and discuss the results. Section 6 shows the implementation of our trust based authentication scheme. Finally, section 7 concludes our work.

Related Work

Several researches and studies were done regarding trust models for different purposes and in different environments. Among those, we choose several prominent ones to discuss their strong and weak points, and build our model accordingly.

In their work, Mihaela et al [START_REF] Ion | A Peer-to-Peer Multidimensional Trust Model for Digital Ecosystems[END_REF] propose a new trust model for DEs which has several innovative features. The model is based on the concept of social networks and addresses trust at different levels: user, data, service and node. The model allows fast bootstrapping of trust by importing existing trust relationships from outside DE systems and by relying on certificates issued by trusted authorities external to the DE. Furthermore, trust can be measured in a variety of contexts by using user-defined tags â folksonomy . The model abstracts from specific reputation algorithms by providing necessary interfaces for plugging-in those on one's own choice.

Yan Lindsay et al present in their paper [START_REF] Lindsay Sun | Information Theoretic Framework of Trust Modeling and Evaluation for Ad Hoc Networks[END_REF] an information theoretic framework to quantitatively measure trust and model trust propagation in ad hoc networks. In the proposed framework, trust is a measure of uncertainty with its value represented by entropy. They develop four Axioms that address the basic understanding of trust and the rules for trust propagation. Based on these Axioms, they present two trust models: entropy-based model and probability-based model, which satisfy all the Axioms. Techniques of trust establishment and trust update are presented to obtain trust values from observation. The proposed trust evaluation method and trust models are employed in ad hoc networks for secure ad hoc routing and malicious node detection. A distributed scheme is designed to acquire, maintain, and update trust records associated with the behaviors of nodes forwarding packets and the behaviors of making recommendations about other nodes. Simulations show that the proposed trust evaluation system can significantly improve the network throughput as well as effectively detect malicious behaviors in ad hoc networks.

With the growing popularity of wireless mobile ad hoc networks (MANETs), many security concerns have arisen from MANETs especially in that misbehaving nodes pose a major threat during the construction of a trusted network. A reputation-based trust system can track the behavior of nodes and thereby proceed by rewarding well-behaving nodes and punishing misbehaving ones. However, existing techniques are usually either energy-consuming or complicated since the relevant reputation information is propagated throughout the network. In their paper, Yonglin et al propose [START_REF] Ren | Modeling and Managing the Trust for Wireless and Mobile Ad hoc Networks[END_REF] a novel trust computation and management system, called TOMS, which not only establishes the new concepts of trust and community but also includes both the trust computation model and trust management mechanism.

In order to obtain more suitable trust evaluation in MANETs, Junhai et al [START_REF] Luo | A trust model based on fuzzy recommendation for mobile ad-hoc networks[END_REF] suggest that the measurement and computation of trust to secure interactions between mobile nodes is crucial for the development of trust mechanisms. The calculation and measurement of trust in unsupervised ad-hoc environment involves complex aspects such as credibility rating for opinions delivered by a node, the honesty of recommendations provided by a mobile node, or the assessment of past experiences with the node one wishes to interact with. The deployment of suitable algorithms and models imitating fuzzy logic can help to solve these problems. In this paper, RFSTrust, a trust model based on fuzzy recommendation similarity, is proposed to quantify and to evaluate the trustworthiness of nodes, which includes five types of fuzzy trust recommendation relation ships based on the fuzzy relation theory and a mathematical description for MANETs. Fuzzy logic provides a natural framework to deal with uncertainty and the tolerance of imprecise data inputs for the subjective tasks of trust evaluation, packet forwarding review and credibility adjustment. Theoretical analysis and experimental results show that RFSTrust is still robust under more general conditions where selfish nodes cooperate in an attempt to deliberately subvert the system, end-to-end packet delivery ratio more quickly, and decreases the average energy consumes more effectively. The effect of node rating data's sparsity can be greatly reduced and show the excellent performance on typical data set.

According to Masthoff [START_REF] Masthoff | Computationally Modelling Trust: An Exploration[END_REF], Trust is a popular and much disputed topic in various research communities. In his paper, he attempts to integrate existing knowledge on trust into a simple computational model. The model incorporates the impact of direct experiences, reputation, stereotypes, empathy and user characteris-tics on trust. He also presents the results of two exploratory experiments testing and improving aspects of the model. Felix Gomez Marmol et al [START_REF] Gomez Marmol | Security threats scenarios in trust and reputation models for distributed systems[END_REF] describe in their review the different scenarios where trust evaluation is threatened by different attacks of malicious users. Examples of those attacks are attacks of individual malicious users, attacks of collective malicious users, attacks of malicious collective with camouflage, malicious spies, sybil attacks, man in the middle attacks, malicious pretrusted peers, partially malicious users, etc. The paper also explains how affecting with peer's reputation can also become a security threat, such as driving down the reputation of a good user. Finally, the paper discusses how some of the known trust models deal with those trust evaluation threats. The trust models mentioned in that section were: EigenTrust, PeerTrust, BTRM-WSN , and PowerTrust. Finally, it concludes with a table summarizing the vulnerabilities and resiliencies of those trust and reputation models against the mentioned attacks.

Our Proposed Trust Model

Before an authentication takes place, each of the communicating parties needs to trust the others, hence the need of a trust evaluation model. Many trust models for smart environments have been proposed in the last decade, some of them presented context awarenessbased security [START_REF] Moloney | A Context-aware Trust-based Security System for Ad Hoc Networks[END_REF], while others concentrated on the quick convergence of updated trust values [START_REF] Ajay | Rapid Trust Establishment for Pervasive Personal Computing[END_REF][19] [START_REF] Sheikh | A Risk-aware Trust Based Secure Resource Discovery (RTSRD) Model for Pervasive Computing[END_REF], addressing the different situations with similar approaches. In our model, we tried to work on the points that are missed or lacking in the previously mentioned trust models, such as the guarantee of the accuracy of trust values, the rapid convergence of trust values during estimation, the compatibility with low-resource devices, and the adaptivity to the needs and abilities of dependent people.

Our proposed model has the following properties:

1. an enhanced compatibility with low-resource devices, since it uses the memory-efficient elliptic curve cryptography and a bandwidth-efficient trust evaluation model;

2. useful features like service-dependent trust, mentioned in [START_REF] Ghorbel | Secured and Trusted Service Provision in Pervasive Environment[END_REF];

3. a lightweight security monitoring to prevent malicious attackers that attempt to forge trust data or alter it in any way;

4. imitation of human rational thinking in evaluating the judgment ability of recommending entities. That means whenever an entity is asked to recommend another, we take a look on how much that entity know the others, which will give a metric representing the judgment ability of this entity, before taking its recommendations as granted;

5. ability to deal with situations which it has no prior experience dealing with, using Questions of Trust;

6. a user rating system that allows to take feedback from users and integrate it in the trust evaluation process.

Figure 1 illustrates the architecture of our proposed trust model with its different modules, showing the steps of the trust evaluation process. The trust updates take into consideration the feedback obtained by user ratings. Judgment is calculated based on the reports of experience messages. Indirect trust is obtained by a multiplicative relation between judgments and recommendation values given by recommendation messages. Direct trust is calculated through a risk assessment based on number of positive and negative actions of the node in question (also provided by experience messages). Finally, the net trust is a linear combination of the direct and indirect trusts. The trust updates occur only on demand or when the trust values have expired. We will explain those modules in more details in the following paragraphs. Next, we are going to present the properties of our trust model and discuss the novel attributes introduced in it.

Minimization of Resource Usage

Our proposed trust model aims to serve smart environments where the hardware equipment have small processors and limited memories. Some trust models such as in [START_REF] Ghorbel | Secured and Trusted Service Provision in Pervasive Environment[END_REF] used a mesh-like approach in calculating trust values. That is each node of the network has always an updated trust table of all nodes in the network. That means that resources are wasted on calculating trust values that may expire before being used. Our proposed model evades this issue by evaluating the trust values only on demand and expiry.

Our proposed model tries also to use as minimal memory resources as possible. Essentially, each node in the network, independently of its nature, will need two matrices. The node will store trust values of nodes communicating with it only. Because of the memory limitation, the node will get rid of trust data related to nodes that had left the neighborhood or went dead in the network.

Trustworthiness

Trustworthiness is a factor that recommenders have to consider in the selection of reliable peers for collaboration. Most approaches in this regard estimates trust base on global user profile similarity or history of exchanged opinions [START_REF] Godoy | Enabling topic-level trust for collaborative information sharing[END_REF]. Trustworthiness is used in this paper to refer to the level of trust of an entity B in respect to a separate entity A. The net trustworthiness is obtained by calculation of two values: direct trust and indirect trust. Direct trust is what is commonly called "Risk Assessment". It is used for dealing with newcomers which the entity has not yet any records of trust evaluation. In case where trust is service-dependent, we added a multiplicative factor, called the Security Action Coefficient (SAC), to the number of negative actions. This coefficient refers to the security level of a service. If the application possesses a high security level, the associated SAC should be high (e.g. 10). On the other hand, if the service didn't need any security, such as a weather broadcasting service, SAC can be as low as 1. Direct trust is obtained using the following equation:

DT = ΣP A i ΣP A i + SAC × ΣN A i . ( 1
)
where P A i represents the number of positive actions done by the node in question and noticed by node i. N A i refers to the number of negative actions, and SAC is the Security Action Coefficient related to the security level of the service.

The indirect trust, representing the recommendations of other nodes, is:

IT = ΣT w i × J i n . (2)
where T w i and J i are the trustworthiness and judgment values corresponding to the node i.

The value of the net trustworthiness is a combination of direct and indirect trust:

T w = α DT × DT + α IT × IT . ( 3
)
where α IT the indirect trust coefficient is:

α IT = T S self T S self + Σ T Si nrecomm × ΣJ i n tot . ( 4
)
and α DT the direct trust coefficient is:

α DT = 1 -α IT . ( 5
)
where T S self refers to the timestamp of the trust value of the node itself, while T S i denotes the timestamp of the trust value of the node i. n tot is the total number of nodes in the subnetwork, whereas n recomm is the number of nodes that responded with recommendations.

Judgment

Judgment is one of the new features introduced that aims to imitate the human behavior in a technical approach. The judgment ability is represented by the overall experience of dealing with the node in question. That experience includes both the total number of control messages exchanged and the total number of actions whether positive or negative.

The judgment related to the number of actions is equal to the total number of actions ΣA i over the maximum number of actions M aximum A, as follows:

J A = ΣA i M ax A , if J A > 1 then J A = 1 . (6)
Similarly, the judgment related to the number of messages exchanged is:

J M = Σ messages i M ax messages , if J A > 1 then J A = 1 . (7)
At last, the overall judgment value is:

J = J A × J M . (8)

Control Messages

In order to control the aspect of trust evaluation and share the trust data, short control messages are used for that purpose. Most of these messages are trigger-based type, except for consistency and hello messages. Those messages are:

1. Recommendation messages: used as request for trust recommendations. The addressed nodes will reply, if possible, with a recommendation reply containing the trustworthiness value requested.

2. Experience messages: used to retrieve information about statistical behavior. The reply contains information about messages exchanged and positive and negative actions, used later in the calculation of Judgment and Direct Trust.

3. Hello messages: periodic messages that are issued to inform the neighbors about self existence.

4. Consistency messages: periodic messages that aim to test the consistency behavior of a certain node in order to prevent any suspicious behavior trying to affect the trust evaluation.

5. Knowledge Migration messages: issued only when a node is about to pass out or leave the network. The message is a notification of the availability of trust data that is going to be lost. Interested nodes will respond by asking for recommendations and experience data.

Questions of Trust

Our basic trust model described in [START_REF] Husseini | A Novel Trust-Based Authentication Scheme for Low-Resource Devices in Smart Environments[END_REF] is suitable to deal with calculating indirect trust based on recommendations of neighbors. However, it may not be efficient enough in case there are several malicious nodes giving false recommendations in order to earn the user's trust. To prevent that, the user can carry out a procedure called "Questions of Trust" before asking for recommendations. The procedure is about sending questions to the targeted nodes and comparing the answers upon receiving them, either with each other or with a set of previously known answers. The questions can be service-related. For example, if the service provided is weather-casting, the questions can be regarding the weather table at a certain geographic location at a certain period of time. Alternatively, the questions can be about the trust tables themselves; and by comparison with trust tables of other neighbors, a contradictory or non-consistent trust values can be detected. With this, the user can be wary of nodes trying to deceive him and thus he will not ask for their services nor their recommendations.

This procedure can be repeated several times to increase the odds of detecting suspicious behaviors. The questions can be also sent directly to targets, or indirectly through other nodes as proxies for masking the identity of the original sender. Questions of Trust are also very handy if the user is entering the network for the first time or is dealing with neighbors he has no prior experience with. This metric, like all other metrics mentioned in this paper, can generally target a large audience of users, independently of their individual capabilities. When we are targeting dependent people in particular for our system, the same mathematical model still applies, but what differs is the type of feedback taken from the users which will become an input for our model. Depending on the type of disability of the user, this feedback can be either vocal or physical (movement of fingers on touchpad or keypad).

User Rating

So far, all the indicators in our trust model are machinebased, but since trust is a concept initially borrowed from humans we decided to add a human-based indicator in our model. User ratings were used before, as in [START_REF] Borzymek | Enriching Trust Prediction Model in Social Network with User Rating Similarity[END_REF][25], for evaluating trust in social networks, eshopping websites, etc. They indicate how much customers are satisfied by the quality of service provided by sellers, or how much a social network user trust another based on real-life activities. However, user rating were not used for trust models in smart environments before, although the user satisfaction/dissatisfaction after using a service from a certain provider may be quite a good indicator of the accuracy of trustworthiness calculated by machine-based indicators. In fact, the quality of service value estimated by the user rating is closely related to the trustworthiness of the service provider, since untrustworthy service providers are not expected to provide a good level of service all the time. Many advantages can be gained from adding user rating to the trust model. First, the machine-based trust models do not take into consideration the human satisfaction/dissatisfaction of a certain service, which is generally a good indicator whether the service is real or fake, and whether there is an encouraging reason for a new user to try this service with this particular service provider. Second, user ratings, unlike trust recommendations, last much more longer and therefore presents a good alternative in case of the absence of the latter. Therefore, we can consider the user rating as a human recommendation and integrates it in our formula to calculate trust. Similarly to how we defined the weights of direct and indirect trusts, the amount of participation of human recommendation in the final value of trust (also referred as Human-Machine Trust) will depend on the numbers of users that have rated the target, and how many times they used the service in question provided by the target. This can be elaborated by the following formulas:

α H = 1 N R × ΣQ i if Σi ≤ N R (9) or α H = ΣQ i if Σi > N R ( 10
)
Q i = U i U max but if Q i > 1 then Q i = 1 (11) α M = 1 -αH (12)
T HM = α H × T H + α M × T M ( 13
)
where T HM is Human-Machine Trust. T H and T M are respectively Human Trust and Machine Trust, and α H and α M are their respective weights. Σi is the total of nodes that rated the service. N R is a constant for normalizing the number of nodes that have rated the target (typically N R = 10). U i is the number representing how many times the node i used the service from that provider, and U max is the normalizing value (typically U max = 100). Q i is the normalized ratio of the number of utilizations of the service done by node i.

Our user rating mechanism is designed in a way that ensures the fair usage of this feature. Since we cannot expect all users to be rational in their ratings, we tried at least to integrate our own approach to fight against abuse of rating, coming from either users or service providers. In our scheme, the interacting parties are the user, the service provider, the Rating Authority and the Abuse Control Authority. The role of the Rating Authority is only to update and sign rating tables and verify the signatures. So it only stores public keys used for the verification of digital signatures, and it does not store any rating data, as to preserve the decentralized nature of the network. Each service provider will store his own rating table(s) (one rating table per service provided), and since the table(s) is digitally signed by the Rating Authority, it will prevent any party from modifying the ratings. Another positive point about the User Rating feature is that they are permanent even if the user who did the rating leaves the network, since the service provider always carry the rating tables along; whereas trust tables depend on the presence of the node, who is going to recommend, in the neighborhood. The Abuse Control Authority is an optional entity that aims to prevent abuse coming from either users or service providers. If the service provider abstain from giving the user his rating table after service provision (i.e. preventing the user from rating him), the user can report him to the Abuse Control Authority using the service receipt he obtains at the beginning of service provision. On the other hand, if the user abuses his right of rating repetitively, he'll be detected by the rating investigation carried out by the Abuse Control Authority. The abuse control mechanism, like spam filtering, can be an automated task based on pattern recognition in abuse behaviors.

To begin with, we assume that we already have a key management scheme enabled and working in the network. The process of the user rating goes as it is shown in Figure 2. First, the user get a copy of the service provider's public key from the Rating Authority. Next, the user will request the service from the service provider. The service provider will send him a service receipt digitally signed and the service will be executed. After the service provision is complete, the service provider will take the Rating Authority's digitally signed copy of his rating table, sign it with his private key, and send the double-signed copy to the user. The user will verify the signatures, add his signature to the rating table, rate and sign his own rating, and send his rating and the rating table to the Rating Authority. The authority will verify all signatures, and if they are valid it will update the rating table based on the new user rating, and send the updated rating table back to the service provider after signing it. To defend against replay attacks, all digitally-signed exchanged data implicitly include a timestamp and/or a nonce.

Figure 3 shows two possible cases of rating abuse. In the first case the abuse is coming from the service provider, where the server abstains from giving the rating table to the user, thus blocking away the user's right of rating. As a deterrent for such situations, we allow the user to file a complaint by reporting to the Abuse Control Authority what happened, using his service receipt as proof of his right to rate. In the second case in the illustration, the abuse is coming from the user where he is giving untruthful ratings to increase/decrease the reputation of a certain service provider. In both abuses, the Abuse Control Authority registers the abuse either by received a justified report or by detecting using abuse pattern recognition. The algorithm describing the mechanism of the Abuse Control Authority is out of the scope of this article.

EC-SAKA Authentication Protocol

Elliptic Curve-based Secure Authenticated Key Agreement protocol (EC-SAKA) was proposed by [START_REF] Pierre | A Fast and Secure Elliptic Curve Based Authenticated Key Agreement Protocol For Low Power Mobile Communications[END_REF] in order to suit the needs for low-resource mobile devices that face difficulties when dealing with large-sized cryptographic keys. This protocol benefits from the power of elliptic curve discrete logarithm problem that is even more powerful than traditional discrete logarithm problems. It is based on an asymmetric Diffie-Hellman scheme to generate a common secret key without exposing it to eavesdroppers. Indeed, the asymmetry in this approach Fig. 3 Two cases involving rating abuse and Abuse Control Authority is to prevent malicious attacks that try to impersonate both of the entities communicating and forward one's data to the other aiming to expose their generated key in the process. EC-SAKA scheme also prevents impersonation attacks. In addition to key agreement, this protocol provides identity verification through El Gamal Signature scheme (ECEGS). This way, it provides identity verification and common key generation at the same time. Figure 4 illustrates the three passes of EC-SAKA protocol.

In Fig. 4, P and Q represent two public points on the elliptic curve chosen. n is large public number chosen by Bob and Alice, used to define the elliptic curve. a and b represent two secret values, each known only by its respective owner, Alice and Bob. A and B are the respective public keys for Alice and Bob generated from a and b. h() is a hashing function, and n is sufficiently large number. K is the common secret key that Alice and Bob will be able to generate at the end of this 3-way authentication.

The EC-SAKA protocol will be used in our scheme in order to authenticate entities and make them confirm each other identities. The authentication phase comes after a level of trust is already established between the Fig. 4 The three passes of EC-SAKA protocol entities. If the result of the trust phase is higher than a specific threshold, which is either dependent of the service or globally set by the administrator, the application will proceed to the authentication phase. Otherwise, the application won't pass to the authentication phase and the service won't be provided, because one of the communicating party appears untrustworthy to the other. The process of trust evaluation is performed using the trust model that was already detailed in the previous section.

Simulation and Discussion

The object of the judgment value is to increase the accuracy of the trust calculation. Not only that, the judgment value helps making the net trust evaluation converge quickly. Unlike the trust evaluation techniques [START_REF] Luo | A trust model based on fuzzy recommendation for mobile ad-hoc networks[END_REF][17] [START_REF] Ajay | Rapid Trust Establishment for Pervasive Personal Computing[END_REF][19] [20][21] where the trust metrics oscillate before reaching a stable value, we made a simulation that demonstrates how the net trust in our model instantaneously reflects the variation in the trust metrics.

The reason that the trust values in our trust model converges rapidly and do not oscillate like in the others is the fact that we are already depending on an alwaysready judgement database. This database consists of two elements: the experience gathered by monitoring and reporting the negative actions performed by other nodes during any communication, and the user feedback through their rating which will enrich this database and makes the trust evaluation easier and faster.

Fig. 5 The variation of net trust and judgment with respect to positive actions Fig. 6 The variation of net trust and judgment with respect to negative actions Figures 5 and6 show the graphs of two simulations: the first one represents the variations of net trust and judgment when the positive actions are increasing, while the second shows the trust variations when the negative actions are increasing. The judgment value increase linearly with the number of actions, positive or negative. That's because the judgment represents the experience which is directly related to the number of actions. On the other hand, the net trust increases/decreases faster when the number of actions increases. The reason is that the grown experience increases the weight of indirect trust (recommendations) and accelerates the variations of net trust. The trustworthiness in our model directly reflects the behavior of the nodes, while in other trust models it takes some time oscillating before converging to an accurate value [START_REF] Ghorbel | Secured and Trusted Service Provision in Pervasive Environment[END_REF].

In comparison with other trust models, our proposed model is using metrics that directly reflect the present and the past line of actions committed by an entity in the network. In addition to that, our trust model tends to decrease the overall power consumption since it triggers trust updates only on demand or expiration. Furthermore, the security monitoring part of the model protects the entities from potential malicious attacks. Questions of Trust also assist in preventing malicious manipulations that affect trust values. Finally, the user rating adds a human-based feature that allows human users to assist in trust decisions based on the satisfaction obtained from using services. Further simulations on real scenarios are taking place by testing the model on dependent users in a special residence in the aim of validation and confirmation of the proprieties of the model.

Implementation

In order to validate the proposed scheme composed of an authentication module and a trust module, we have implemented it using Java language and Eclipse IDE platform. A server in our architecture is the device that provides the service to other nodes. It can be a computer, an RFID reader, or even a sensor. The server part of the implementation contains no graphical interface, since it only calculates different data and communicates them to other nodes. On the other hand, the client part includes a graphical interface and can run on almost any mobile device that supports Java. We chose to test our implementation on Android mobile phones using Android SDK tools in Java. The code itself is implementable on simpler platforms than Android phones, such as active RFID tags and sensors. Android systems remain easier to work on since the implementation will only take place software-wisely.

For the authentication module, it has the responsibility to establish secret key generation and identity confirmation using the EC-SAKA protocol. The implementation of this module was done through 4 Java classes. The first two ECCPoint and ECCurve take care of all mathematical definitions and calculations related to Elliptic Curves, which are the base of the EC-SAKA protocol. The other two classes, Requester and Provider are used by the authenticating parties to exchange the messages needed for the establishment of the authentication. Since EC-SAKA uses a 3-way asymmetric scheme, the messages sent received by one authenticating party are not alike to those sent and received by the other party, thus the need of two different classes for the two parties. Figure 7 shows the class diagram of the authen-tication module. The other module implemented, the trust module, uses two Java classes. The first class is for exchanging control messages. The other class uses the information provided by recommendation and experience messages to calculate direct trust, indirect trust and judgment. This class also calculates the trust weights, described in the previous section, in order to evaluate the net trust.

Our implementation of this scheme takes part in a project called "Cohabit", a smart environment project for dependent people. The project takes place in a particular residence for disabled people called ADEP. The services provided in that smart environment are daily services needed by dependent people, such as opening the door, turning on/off the light, closing the curtains, etc. The security modules we have developed evaluate the trust between the users and their environment before giving them access to use those services. When the user of the Android phone chooses a service, he will be directed to another window that lists the security steps needed to activate that service. As it appears in Fig. 10, when clicking on 'Evaluate Trust' button, the trustworthiness value of the service provider will be calculated, as described in section 3, and the threshold value set by the administrator will be displayed. For more flexibility, the procedure was divided to several steps, to allow the user to use services on his own risk -if he wishes to -when the trustworthiness value is less than the threshold. The next button uses the EC-SAKA protocol to generate a secret shared key and enable the two nodes to verify each other's identities. If the user sees two checkmarks next to the two buttons, he can now execute the service he has chosen. Some services might not need trust evaluation and authentications, such as date & time and weather-forecast services. In this case, the corresponding buttons for trust and authentication will be disabled and only the service execution button will be enabled. Nevertheless, the actions of non-authenticated services will be recorded in the action matrices.

Conclusion

After emphasizing on the necessity and importance of security, privacy and trust in smart environments, we have demonstrated the effectiveness of Elliptic Curve Cryptography as new candidate in publickey cryptosystems. We have adopted a lightweight authentication protocol called Elliptic Curve Secure Authenticated Key Agreement (EC-SAKA) protocol in our scheme. That authentication protocol uses a 3-pass scheme to generate a common secret key, in addition to an elliptic curve-based digital signature.

Next, we have proposed a new trust model that respects the limitation of tiny mobile devices in terms of resources and bandwidth. Our trust model contain two new features that enhance the process of trust evaluation. The first one is a new trust metric introduced called the ability of judgment. This value tends to imitate the human rational thinking in trust and recommendation acceptance. The second feature is a lightweight security monitoring ability within the trust model to defend against security threats. Furthermore, we introduced two new features in our enhanced trust model. One is Questions of Trust that enable users to ask service-related questions to the service provider prior to service execution, in order to detect any fraud or fake services. The other is User Rating that gives the model a more humanistic side by integrating human feedback in the trust evaluation process.

Finally, we have implemented our scheme into independent modules using Java language. Our trust based authentication scheme will be embedded on Android smartphones to be used by dependent people in their residential areas.

Fig. 1

 1 Fig. 1 Our proposed Trust model with its basic parts

Fig. 2

 2 Fig. 2 Diagram showing actions performed between user, service provider and rating authority

Fig. 7

 7 Fig. 7 The different classes used in the implementation of the authentication module

  Figures 8,[START_REF] Tentori | Supporting Quality of Privacy (QOP) in Pervasive Computing[END_REF] and 10 shows several snapshots the client application on the Android phone listing the existing services spotted in the environment.

Fig. 8

 8 Fig. 8 Snapshot of the service provision in the Android app

Fig. 9

 9 Fig. 9 Snapshot of the Trust action matrix in the Android app

Fig. 10

 10 Fig. 10 Snapshot of the trust-authentication scheme in the Android app