Second order decomposition model for image processing : numerical experimentation

Keywords: Second order total variation, image decomposition, variational method, inf-convolution, texture extraction Classification: 65D18, 68U10, 65K10

 where a second order image decomposition model to perform denoising and texture extraction has been studied. Here we perform some numerical experimentation to make the behavior of the model as clear as possible. For highly textured images the model gives a two-scale texture decomposition.

Introduction

The Rudin-Osher-Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is the most famous denoising variational model. It involves a total variation regularization term that preserves discontinuities. The image to recover/denoise is split in two parts namely the noise and a smooth part (here a function of bounded variation [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF][START_REF] Attouch | Variational analysis in Sobolev and BV spaces[END_REF][START_REF] Evans | Measure theory and fine properties of functions[END_REF]). A lot of people have investigated such decomposition models based on variational formulation, considering that an image can be decomposed into many components, each component describing a particular property of the image ( [START_REF] Aubert | Modeling very oscillating signals. Application to image processing[END_REF][START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Osher | Image decomposition and restoration using total variation minimization and the h 1 norm[END_REF][START_REF] Osher | Modeling textures with total variation minimization and oscillating patterns in image processing[END_REF][START_REF] Osher | Image denoising and decomposition with total variation minimization and oscillatory functions. special issue on mathematics and image analysis[END_REF][START_REF] Yin | A comparison of three total variation based texture extraction models[END_REF]] and references therein for example).

In [START_REF] Bergounioux | A full second order variational model for multiscale texture analysis[END_REF] we have presented a second order model where the first order classical total variation term and a second order total variation one are involved. The suitable functional framework is the space of functions with bounded hessian introduced in [START_REF] Demengel | Fonctions à hessien borné[END_REF] (see also [START_REF] Bergounioux | A second-order model for image denoising[END_REF][START_REF] Bergounioux | A full second order variational model for multiscale texture analysis[END_REF][START_REF] Bergounioux | On poincaré-wirtinger inequalities in bv -spaces[END_REF]). To achieve this goal K. Bredies and al. have recently introduced a second order generalized total variation definition [START_REF] Bredies | Total generalized variation[END_REF][START_REF] Bredies | Properties of l 1 -tgv 2 : the one-dimensional case[END_REF][START_REF] Bredies | Regularization of linear inverse problems with total generalized variation[END_REF] which is a nice compromise/mixture between the first and second order derivatives. In [START_REF] Bergounioux | Mathematical analysis of a inf-convolution model for image processing[END_REF], we have investigated an inf-convolution type model for texture extraction, mixing first and second order regularization terms. This is not more e cient than the Total Generalized Variation approach of [START_REF] Bredies | Total generalized variation[END_REF][START_REF] Bredies | Properties of l 1 -tgv 2 : the one-dimensional case[END_REF][START_REF] Bredies | Regularization of linear inverse problems with total generalized variation[END_REF] for denoising. However, it provides a decomposition of the image at di erent scales what the other model does not a priori. We focused on this decomposition that provides a multiscale description of textured images and gave a rough mathematical analysis of this model in [START_REF] Bergounioux | Mathematical analysis of a inf-convolution model for image processing[END_REF].

More precisely, we assumed that an image (in L 2 ( )) can be split in three components: a smooth (continuous) part v, a cartoon (piecewise constant) part u and an oscillating part w that should involve noise and/or fine textures. The oscillating part of the image is included in the remainder term w = u d ≠ u ≠ v, while v is the smooth part (in BH( )) and u belongs to BV ( ): we hope u to be piecewise constant so that its jump set gives the image contours. For highly textured images, the model provides a two-scale texture decomposition: u can be viewed as a macro-texture (large scale) whose oscillations are not too large and w is the micro-texture ( much more oscillating) that contains the noise. The same model has been considered in [START_REF] Benning | Higher-order tv methods enhancement via bregman iteration[END_REF] from a di erent point of view and the numerical realization has been performed with Bregman-type algorithms.

In this paper, we perform numerical experimentations to check the behavior of the method that we do not completely understand. Preliminary results et conjectures have been set in [START_REF] Bergounioux | Mathematical analysis of a inf-convolution model for image processing[END_REF]. The aim of these numerical tests is to validate some theoretical results (as uniqueness in particular cases) and reinforce conjectures. From that point of view, the numerical method to compute the solutions is far to be optimal. In particular, the methods we used to compute the di erent projections are not the most e cient with respect to the Chambolle-Pock algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] or Bregman-type algorithms.

Presentation of the model

We now assume that u d belongs to L 2 ( ) and that the image we want to recover can be decomposed as u d = w + u + v where u, v and w are functions that characterize di erent parts of u d . Components belong to di erent functional spaces: v is the (smooth) second order part and belongs to BH( ), u is a BV ( ) component and w oe L 2 ( ) is the remainder term. We consider the following cost functional defined on BV ( ) ◊ BH( ):

F ⁄,µ (u, v) = 1 2 Îu d ≠ u ≠ vÎ 2 L 2 ( ) + ⁄T V (u) + µT V 2 (v), (1)
where ⁄, µ > 0. We are looking for a solution to the optimization problem

inf{ F ⁄,µ (u, v) | (u, v) oe BV ( ) ◊ BH 0 ( ) } (P ⁄,µ )
We refer to [START_REF] Bergounioux | Mathematical analysis of a inf-convolution model for image processing[END_REF] and the references therein for the di erent spaces definition.

We expect v to be the smooth part of the image, u to be a BV ( )\BH( ) function which derivative is a measure supported by the contours and w := u d ≠ u ≠ v oe L 2 is the noise and/or small textures (we shall detail this point later). We have proved in [START_REF] Bergounioux | Mathematical analysis of a inf-convolution model for image processing[END_REF] that problem (P ⁄,µ ) has at least an optimal solution (u ú , v ú ) in BV ( ) ◊ BH 0 ( ). Moreover the dual problem to (P ⁄,µ ) writes

inf woe⁄K1flµK2 1 2 Îu d ≠ wÎ 2 2 . ( 2
)
where where

K 1 = K 1 is the L 2 -closure of K 1 := ) › = div Ï | Ï oe C 1 c ( ), ÎÏÎ OE AE 1 * . ( 3
)
and and

K 2 ∏ K 2 with K 2 is the L 2 -closure of K 2 := ) › = div 2 Â |Â oe C 2 c ( , R d◊d ), ÎÂÎ OE AE 1 * . ( 4
)
The unique solution w ú is the L 2 -projection of u d on the closed convex set ⁄K 1 fl µK 2 :

w ú = ⁄K1flµK2 (u d ) .
Next we have a relation between the solutions to (P ⁄,µ ) and the (unique) solution of the dual problem.

Theorem 2.1. 1. Let w ú be the (unique) solution to the dual problem (P ⁄,µ ) ú :

w ú = ⁄K1flµK2 (u d ) .
Then there exists (ū, v) oe BV ( ) ◊ BH 0 ( ) an optimal solution to (P ⁄,µ ) such that

w ú = u d ≠ ū ≠ v and w ú oe ˆ 2 µ (v) fl ˆ 1 ⁄ (ū) . 2. Conversely, if (ū, v) oe BV ( ) ◊ BH 0 ( ) is any solution to (P ⁄,µ ) then w = u d ≠ ū ≠ v = ⁄K1flµK2 (u d ) . ( 5
)
Here

1 ⁄ (u) = ; ⁄T V (u) if u oe BV ( ) +OE else. and 2 µ (v) = ; µT V 2 (v) if v oe BH 0 ( ) +OE else.
This theorem is useful from many points of view. First, we have a necessary and su cient condition for a pair (u, v) to be a solution. As we explain in the sequel, this theorem is the key theorem for a fixed point method to compute a solution.

Moreover, solutions (u, v) are not unique but we get uniqueness for the remainder part w = u d ≠ u ≠ v. In addition, we have (partial) results in [START_REF] Bergounioux | Mathematical analysis of a inf-convolution model for image processing[END_REF]. We expect numerical results to be consistent and give hints for uniqueness open cases. Let us recall these results thereafter. Theorem 2.2 ([8] Section 4.2). 1. Assume (u 1 , v 1 ) and (u 2 , v 2 ) are two optimal solutions of (P ⁄,µ ). Then, there exists Ï oe BV ( ) fl BH 0 ( ) such that u 2 = u 1 ≠ Ï and v 2 = v 1 + Ï.

2. Let (⁄, µ) be nonnegative real numbers such that ⁄ Ø Îu d Î G and µ >> ⁄, where ΕΠis the Meyer-norm [START_REF] Meyer | Oscillating Patterns in Image Processing and Nonlinear Evolution Equations[END_REF]. Then the C( ) functions are the only solutions to (P ⁄,µ ),where

C( ) := {(u, v) oe BV ( ) ◊ BH 0 ( ) : ÷c oe R u = c and v = ≠c a.e on }.
3 Numerical aspects

Discretized problem and algorithm

We assume that the image is rectangular with size N ◊ M . We note X := R N ◊M ƒ R NM endowed with the usual (normalized) inner product and the associated Euclidean norm

Èu, vÍ X := 1 NM ÿ 1AEiAEN ÿ 1AEjAEM u i,j v i,j , ÎuÎ X := Û 1 NM ÿ 1AEiAEN ÿ 1AEjAEM u 2 i,j . (6) We set Y = X ◊ X.
It is classical to define the discrete total variation with finite di erence schemes as following (see for example [START_REF] Aubert | Mathematical Problems in Image Processing, Partial Di erential Equations and the Calculus of Variations[END_REF]): the discrete gradient of the numerical image u oe X is Òu oe Y and may be computed by the following forward scheme for instance:

(Òu) i,j = 1 (Òu) 1 i,j , (Òu) 2 i,j 2 , ( 7
)
where

(Òu) 1 i,j = ; u i+1,j ≠ u i,j if 1 < i < N 0 if i = 1, N, and (Òu) 2 i,j = ; u i,j+1 ≠ u i,j if 1 < j < M 0 if j = 1, M.
Note that the constraint ˆu ˆn = 0 is involved in the discretization process of the gradient. Therefore, in a discrete setting, the sets K 2 and K 2 coincide. The (discrete) total variation corresponding to 1 (u) is given by

J 1 (u) = 1 NM ÿ 1AEiAEN ÿ 1AEjAEM . . (Òu) i,j . . R 2 , ( 8
)
where .

.

(Òu) i,j . . R 2 = . . ! Òu 1 i,j , Òu 2 i,j ". . R 2 = Ú 1 Òu 1 i,j 2 2 + 1 Òu 2 i,j 2 2 .
The discrete divergence operator -div is the adjoint operator of the gradient operator:

'(p, u) oe Y ◊ X, È≠div p, uÍ X = Èp, ÒuÍ Y .
To define a discrete version of the second order total variation 2 we have to introduce the discrete Hessian operator. For any v oe X, the Hessian matrix of v, denoted Hv is identified to a X 4 vector: " .

(Hv) i,j = ! (Hv)
We refer to [START_REF] Bergounioux | A full second order variational model for multiscale texture analysis[END_REF][START_REF] Bergounioux | A second-order model for image denoising[END_REF] for the detailed expressions of these quantities. The discrete second order total variation corresponding to 2 (v) writes

J 2 (v) = 1 NM ÿ 1AEiAEN ÿ 1AEjAEM Î(Hv) i,j Î R 4 , (9) with Î
F ⁄,µ := 1 2 Îu d ≠ u ≠ vÎ 2 X + ⁄J 1 (u) + µJ 2 (v). (10)
Problem [START_REF] Bergounioux | A full second order variational model for multiscale texture analysis[END_REF] has obviously a solution ũ and ṽ that satisfies the following necessary and su cient optimality conditions

ũ = u d ≠ ṽ ≠ ⁄K1 (u d ≠ ṽ) , ( 11a
) ṽ = u d ≠ ũ ≠ µK2 (u d ≠ ũ) , ( 11b
)
where K 1 and K 2 are the following convex closed subsets :

K 1 = {div p | p oe X 2 , Îp i,j Î R 2 AE 1 'i = 1, . . . , N, j = 1, . . . , M}, ( 12a
)
K 2 = {H ú p | p oe X 4 , Îp i,j Î R 4 AE 1, 'i = 1, . . . , N, j = 1, . . . , M}, (12b)
and Ki denotes the orthogonal projection on K i . These projections are computed with a Nesterov-type scheme as in [START_REF] Weiss | E cient schemes for total variation minimization under constraints in image processing[END_REF] where the maximal number of (inner) iterations has been set to 30 and the accuracy to 10 ≠7 . We refer to [START_REF] Bergounioux | A second-order model for image denoising[END_REF] for more details. As already mentioned our main concern is to check the relevance of the model and highlight some results and/or conjectures about uniqueness or the structure of solutions. Thus, our numerical method is far to be optimal.

The numerical implementation could be widely improved using methods as the ones described in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Combettes | Variable metric forward-backward splitting with applications to monotone inclusions in duality[END_REF][START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF] for example. This leads to the following fixed-point algorithm:

Algorithm -Initialization step. Choose u 0 , v 0 , set 0 < -< 1/2 and n = 1. -Iteration. Define the sequences ((u n , v n )) n as I u n+1 = u n + -(u d ≠ u n ≠ v n ≠ ⁄K1 (u d ≠ v n )) v n+1 = v n + -(u d ≠ u n ≠ v n ≠ µK2 (u d ≠ u n )) . -Stopping test. If max(Îu n+1 ≠ u n Î L 2 , Îv n+1 ≠ v n Î L 2 ) AE Á ( 13
)
where Á > 0 is a prescribed tolerance, or if the iterations number is larger than a prescribed maximum number itmax, then STOP.

We have proved in [START_REF] Bergounioux | A full second order variational model for multiscale texture analysis[END_REF] that for anyoe (0, 1/2), the sequence generated by the algorithm converges to a stationary point, solution of (11) that we generically denote (u ú , v ú ) in the sequel. The tolerance was set to Á = 10 ≠2 so that the stopping criterion is de facto the maximum number of iterations itmax. In the sequel, we have set itmax=10 000 for the 1D case and itmax = 400 for the 2D case.

We do not report on CPU time since all tests have been done with MATLAB © and the code is not optimized. A parallelled C ++ is version is written that reduces the computational time significantly. Moreover, as explained before, the projections can be computed more e ciently with recent algorithms to spare computation time and improve accuracy; nevertheless, we are interested in the behavior of the model and the optimization of projections computation is not our main concern here.

Examples

We use 1D and 2D examples.

For the first (1D) example we set s = s 0 + s 1 + s 2 on [0,1] with s 0 a white gaussian noise with standard deviation ‡ = 0.02 and The second example is a 2D picture of a butterfly and the third one an highly textured image (old wall). We used geometrical images as well but we do not report on them. 1 We present some results and comments in the next subsections.

s 1 = ; 0.4 on [ 3 10 , 6 10 ] 0 elsewhere , s 2 (x) = ; 0.8 x + 0.2 on [0, 1 2 ] ≠1.2 (x ≠ 1) elsewhere.

Initialization process

We have tested many initialization choices for algorithm. Indeed, we have not proved uniqueness (though we conjecture it). So the computed solution is only a stationary point. As we may have many, we may think that the initialization process has a significant influence on the generated sequence.

More precisely, we used u 0 = 0, v 0 = u d , that we call initialization (a) in the sequel, -u 0 = u d , v 0 = 0 that we call initialization (a') in the sequel, -u 0 = 0, v 0 = 0 : initialization (b), -randomized initializations around u d mean value.

Initialization (a) (resp. (a')) provides a stationary pair (u ú , v ú ) such that u ú (resp. v ú ) has null mean value. given by the algorithm satisfies

⁄ u ú = 0. Similarly, if u 0 = u d and v 0 = 0, the pair (u ú , v ú )
given by the algorithm satisfies

⁄ v ú = 0.
Proof. Though we consider a discrete setting we use a continuous setting notation (using for example a piecewise a ne approximation). We first note that

w oe K 1 fi K 2 =∆ ⁄ w = 0 .
We prove the first assertion. Assume that u 0 = 0 and v 0 = u d . It is easy to see by induction that

'n oe N ⁄ u n = 0 and ⁄ (v n ≠ u d ) = 0 . ( 14
) using I u n+1 = u n + -(u d ≠ u n ≠ v n ≠ ⁄K1 (u d ≠ v n )) v n+1 = v n + -(u d ≠ u n ≠ v n ≠ µK2 (u d ≠ u n )) .
Passing to the limit we get ⁄ u ú = 0 and

⁄ (v ú ≠ u d ) = 0 .
The second assertion is proved similarly.

Proposition 3.1 yields that the BV -part u ú (or the BHpart v ú ) belongs to the discrete Meyer space G (see [START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF]) if we perform the appropriate initialization step. This means it is an oscillating function. More precisely, choosing u 0 = u d , v 0 = 0 gives a BHpart that belongs to G. This is not what we want, since the BHpart should not be oscillating. Therefore, we shall never use such an initialization. Initializations (a) and (a') seem to give di erent results from initialization (b). We shall see in the sequel that the di erence is small if the iteration number is large enough. Therefore, we think that the initial guess has no influence on the result, but only on the convergence speed. , µ = 1 u 0 = 0, v 0 = u d 1.716468e-04 6.577422e-03 6.381764e-04 8.619797e-05

Initialization F ,µ (u ⇤ , v ⇤ ) kw ⇤ k L 2 T V (u ⇤ ) T V 2 (v ⇤ ) = 10 3 , µ = 10 2 u 0 = 0, v 0 =
u 0 = 0, v 0 = 0 1.144536e-
u 0 = u d , v 0 = 1.
716467e-04 6.577421e-03 6.381750e-04 8.619794e-05

u 0 = 0, v 0 = 0 1.
ú = u d ≠ u ú ≠ v ú .
-The number of inner iterations (in the projections) has been set to 100, to achieve the best accuracy as possible.

We can see on Figure 3 (1D) the oscillating e ect of initialization u 0 = 0, v 0 = u d : 2 give the computed pairs with initializations (a) (a') (b) and a randomized initialization around the mean value of u d for a 2D example.

(a) Original signal (without noise) - F ⁄,µ (u ú , v ú ) = 8.19e ≠ 07 (b) Initialization u0 = 0, v0 = u d (a) - F ⁄,µ (u ú , v ú ) = 1.65e ≠ 06 (c) Initialization u0 = u d , v0 = 0 (a') - F ⁄,µ (u ú , v ú ) = 1.43e ≠ 06 (d) Initialization u0 = 0, v0 = 0 (b) - F ⁄,µ (u ú , v ú ) = 8.95 ≠ 07
Initialization F ,µ (u ⇤ , v ⇤ ) kw ⇤ k L 2 T V (u ⇤ ) T V 2 (v ⇤ ) Error # it. = 1, µ = 10 u 0 = 0, v 0 = u d 23.
ú = u d ≠ u ú ≠ v ú .
The error reported in column 5 is defined as in (13) (and should be less that 10 ≠2 ). The blue (grey) lines of Table 2 show what we call the optimal solution, that is the computed pair whose cost functional value is the lowest. We observed that -the randomized initialization gives the same result as initialization (b), -the component w ú = u d ≠ u ú ≠ v ú is always the same, which is consistent with the theoretical result of uniqueness, -the values of the cost functional may be quite close and the computed pairs quite di erent: see for example ⁄ = 5, µ = 7 (and figure 4), -initialization (b) gives a pair (u b , v b ) such that neither u b nor v b has null mean value.

In the sequel, (u a , v a ) denotes the pair given by the algorithm with initialization (a) and (u b , v b ) the one given by the algorithm with initialization (b). Moreover, we set the signed relative error as

"F ⁄,µ = F ⁄,µ (u a , v a ) ≠ F ⁄,µ (u b , v b ) min(F ⁄,µ (u a , v a ), F ⁄,µ (u b , v b )) . ( 15
)
Figures 5 and6 show the behavior of "F ⁄,µ with respect to ⁄ and µ. and (b) for ⁄ = 7, µ = 9, as the number of iterations increases. We observe convergence : the cost functional is exponentially decreasing. 4. First and second order total variations T V and T V 2 of pairs given by initializations (a) and (b) for ⁄ = 7, µ = 9, as the number of iterations increases. Here Ï = u b ≠ ua = v b ≠ va and the error is given by the stopping criterion (13) of Algorithm.

# it. T V (ua) T V (u b ) T V (') T V 2 (va) T V 2 (v b ) T V 2 (')
Though F ⁄,µ (u a , v a ) ƒ F ⁄,µ (u b , v b ) the pairs (u a , v a )
and (u b , v b ) may be very di erent. More precisely, we have u b = u a ≠ Ï and v b = v a + Ï. Though the computed function Ï k at iteration k is not a constant function (see Figure 7), we infer that Ï k converges to a constant function as the iteration number increases. Indeed, we have numerically observed (see Table 4) that both T V (Ï) and T V 2 (Ï) decreases to 0 as the iteration number increases. Nevertheless, we can perform only a limited number of iterations. So the computed solutions di er from a (small) piecewise constant function (see Figure 7). In addition, it is numerically confirmed that w a = u d ≠ u a ≠ v a = w b (what was theoretically proved).

Convergence

We chose -= 0.25 in the fixed point algorithm and we always observed convergence. We set the maximal number of iterations quite large but we noticed that the solution is satisfactory with less iterations (400 for 2D case and 1000 for 1D case).

# it. F ,µ (ua, va) F ,µ (u b , v b ) | F ,µ | = 1, µ =
# it. T V (ua ) T V (u b ) T V (') T V 2 (va ) T V 2 (v b ) T V 2 (') Error (a) Error (b) = 1, µ =

Sensitivity with respect to sampling and quantification

We first investigate the sensitivity with respect to sampling. Table 7, figures 13 and 14 show that the model is robust with respect to sampling. Here, we have discretized the analogical signal of example 1D with 10 3 , 10 4 and 10 5 points respectively.

µ F ,µ (u b , v b ) kw b k L 2 T V (u b ) T V 2 (v b ) 10
⁄,µ (u ⁄,µ , v ⁄,µ ) = 1 2 Îu d ≠ u ⁄,µ ≠ v ⁄,µ Î 2 + ⁄T V (u ⁄,µ ) + µT V 2 (v ⁄,µ ) = 1 2-2 Î-u d ≠ -u ⁄,µ ≠ -v ⁄,µ Î 2 + ⁄ - T V (-u ⁄,µ ) + µ - T V 2 (-v ⁄,µ ) = 1 -2 3 1 2 Î-u d ≠ u ⁄,μ ≠ v ⁄,μ Î 2 + -⁄T V (u ⁄,μ ) + -µT V 2 (v ⁄,μ ) 4 = 1 -2 F ⁄,μ (u ⁄,μ , v ⁄,μ
) with ⁄ = -⁄ and μ = -µ. Table 8. Sensitivity with respect to quantification-Initialization (b) -itmax = 400. The model is robust with respect to quantifications as expected.

Sensitivity with respect to parameters

As mentioned before the computed stationary pair depends on the initialization guess via the convergence speed. We consider three cases and we illustrate them on test 2D (Butterfly).

-If µ << ⁄, then initialization (a) : u 0 = 0 and v 0 = u d is the best choice to make the algorithm converge quickly. So we use this initialization to get the solution (u ú , v ú ). In this case, the BV part is close to 0. However, we note that if we fix µ then T V (u ú (⁄, µ)) decreases to 0 and T V 2 (v ú (⁄, µ)) increases to become constant (see Figure 15) when ⁄ ae +OE. This means that if ⁄ is large then u ú is constant. As we know that u ú has a null mean value, then u ú = 0. [START_REF] Bredies | Properties of l 1 -tgv 2 : the one-dimensional case[END_REF]. In this case the cartoon part tends to a constant function as ⁄ = µ increases (the total variation tends to 0) while the second order total variation of the smooth part converges to a limit (not equal to 0, a priori) (Butterfly test) -If ⁄ << µ, then we choose initialization (b) : u 0 = 0 and v 0 = 0 to get the solution. The behavior is similar to the case µ < ⁄: if we fix ⁄, then T V (u ú (⁄, µ)) increases to a constant value and T V 2 (v ú (⁄, µ)) converges to 0 as µ ae +OE (see figure 19). This means that if µ is large enough then solution is always the same : v ú is an a ne function. µ 10. Wall-example, initialization u0 = 0, v0 = 0 -800 iterations -The larger ⁄ is, the more the texture information is involved in the L 2 part, w ú .

(a) Îw ú Î L 2 (b) T V (u ú ) (c) T V 2 (v ú )
= µ = 0.5, 1, 2, 3 • • • 25-init (a) - 800 iterations -Butterfly (a) Cost functional (b) Îw ú Î L 2 (c) T V (u ú ) (d) T V 2 (v ú )
(a) T V (u ú ) (b) T V 2 (v ú )
F ,µ (u ⇤ , v ⇤ ) kw ⇤ k L 2 T V (u ⇤ ) T V 2 (v

Conclusion

The model is well adapted to texture extraction. In the case, where the data is noiseless and/or is not too much textured, the decomposition given par ⁄µ and initialization u 0 = v 0 = 0, gives a cartoon part which is piecewise constant as expected. This means that u = q i u i 1 i where t i i is the contour set. In this case, the remainder L 2 term is the texture and/or noise. The decomposition is robust with respect to quantification, sampling and is always the same for any µ >> ⁄, once ⁄ has been chosen.

In the case where the image is highly textured the model provides a twoscale decomposition. The T V part represents the macro-texture and the L 2 part the micro-texture and/or noise. The scaling is tuned via the ratio fl = ⁄ µ .

The notion of highly textured may be quantified par the G-norm. In our 2D examples, the butterfly G norm was ƒ 7.71 and the wall one was ƒ 4.92.

Figure 21 shows the behavior of the di erent components with respect to ⁄ and µ. We have chosen the 1D noiseless case, to see the multi-scale e ect on components u and w when µ < ⁄. Moreover, the initialization process has no influence on the solution (up to a constant function) but rather on the algorithm speed. The choice has to be made with respect to the parameters: roughly speaking, if ⁄ < µ we choose u 0 = 0, v 0 = 0 and if ⁄ Ø µ we choose u 0 = 0, v 0 = u d . Finally, we have observed (numerically) that the L 2 -component w is unique.

Next issue is to speed up the algorithm (using more performant algorithms) and set an automatic parameter tuning with respect to data properties (G norm, Signal to Noise Ratio , and so on.) From the theoretical point of view, we infer that problem (P ⁄,µ ) has a unique solution (up to constant functions) but the question is still open.

Fig. 1 .

 1 Fig. 1. 1D example -1000 points. First line : signal s, second line : BH-part s2, third line : BV -part s1 and last line ; L 2 noise s0.

Fig. 2 .

 2 Fig. 2. 2D examples

Proposition 3 . 1 .

 31 Assume u 0 = 0 and v 0 = u d . Then any solution (u ú , v ú )

  u d 6.401836e-06 1.424060e-03 5.134659e-03 2.532040e-05 u 0 = u d , v 0 = 6.242638e-06 1.434416e-03 5.076678e-03 1.371843e-05 u 0 = 0, v 0 = 0 6.235598e-06 1.437834e-03 5.022263e-03 1.796523e-05 Random 6.226676e-06 1.441939e-03 5.007722e-03 1.793599e-05 = 10 2 , µ = 10 2 u 0 = 0, v 0 = u d 2.230749e-05 4.897618e-03 8.834750e-04 1.479412e-04 u 0 = u d , v 0 = 2.201165e-05 4.972975e-03 8.022872e-04 1.623540e-04 u 0 = 0, v 0 = 0 2.200473e-05 4.977910e-03 7.968618e-04 1.646320e-04 Random 2.200287e-05 4.978969e-03 7.959566e-04 1.648234e-04 = 1, µ = 10 1 u 0 = 0, v 0 = u d 1.144536e-04 1.027027e-02 8.406575e-06 5.330783e-04 u 0 = u d , v 0 = 1.144536e-04 1.027027e-02 8.406575e-06 5.330783e-04

Fig. 3 .

 3 Fig.3. Example 1D without noise (s1 + s2). The number of inner iterations (to compute projections) is 100 and the number of outer iterations is 50 000. We present the case where ⁄ = 10 ≠3 , µ = 10 ≠2 with di erent initializations. Note that both u ú and w ú have null mean value for initialization (a) (u0 = 0, v0 = u d ). We recover the original decomposition with initialization (b) (u0 = 0, v0 = 0). Initialization (a') (u0 = u d , v0 = 0 ) gives a di erent solution : we shall comment in the sequel.

Figure 4

 4 Figure 4 and Table2give the computed pairs with initializations (a) (a') (b) and a randomized initialization around the mean value of u d for a 2D example.

Fig. 4 .

 4 Fig. 4. BH-part v (first line) and BV -part u (second line) given by initializations (a), (a'), (b) and random for ⁄ = 5, µ = 7 -Butterfly example with 400 iterations.

Fig. 5 .

 5 Fig. 5. Behavior of "F ⁄,µ for 400 iterations (Butterfly example) with respect to ⁄ and µ. If ⁄ and µ are large enough (⁄ > 0.1 and µ > 0.1 for example), both optimal values are very close.

Fig. 6 .

 6 Fig. 6. Behavior of "F ⁄,µ for 400 iterations (Butterfly example). Slices ⁄ = 0.1 (red dotted line) and µ = 0.1 (blue solid line) .

Fig. 7 .

 7 Fig. 7. Di erence between the solutions given by initializations (a) and (b) for ⁄ = 7, µ = 9 and 5000 iterations . ÎÏÎ2 = 0.1518, TV (Ï) = 0.1484, TV 2 (Ï) = 0.1803. The function Ï is close to be piecewise constant as we see it on the gradient norm. This is consistent with Theorem 2.2

Fig. 8 .Fig. 9 .Fig. 10 .

 8910 Fig. 8. Behavior of the cost functional for ⁄ = 7, µ = 9, 100 iterations and initialization (a) (u0 = 0, v0 = u d )

2 . 2 Fig. 11 .Fig. 12 .

 221112 Fig. 11. BV component ua and Ï corresponding to Table 6 -800 iterations. Functionn Ï turns to be constant. The choice ⁄ = 10, µ = 15 gives a satisfactory cartoon part.

Fig. 13 .

 13 Fig. 13. Test 1D (with noise) -Pair given by initialization (b) for ⁄ = 0.1, µ = 1, 10 000 iterations and 10 3 points sampling.

Fig. 14 .

 14 Fig. 14. Test 1D (with noise) -Pair given by initialization (b) for ⁄ = 0.1, µ = 1, 10 000 iterations and 10 4 points sampling.

F

= 7 , 13 = 10 , µ = 2

 713102 µ = 9 Initialization (a) ku ↵ ,↵µ ↵u ,µ k1/↵ 3.0291e-01 3.0291e-01 kv ↵ ,↵µ ↵v ,µ k1/↵ 3.1291e-01 3.1291e-01 = 7, µ = 9 Initialization (b) ku ↵ ,↵µ ↵u ,µ k1/↵ 1.8006e-01 1.8006e-01 kv ↵ ,↵µ ↵v ,µ k1/↵ 1.7924e-01 1.7924e-01 = 10, µ = 2 Initialization (a) ku ↵ ,↵µ ↵u ,µ k1/↵ 8.0280e-15 8.4421e-15 kv ↵ ,↵µ ↵v ,µ k1/↵ 1.1324e-13 1.4552e-Initialization (b) ku ↵ ,↵µ ↵u ,µ k1/↵ 3.9216e-03 4.5475e-14 kv ↵ ,↵µ ↵v ,µ k1/↵ 1.1324e-13 1.0914e-13

Fig. 15 .

 15 Fig. 15. Generic L 2 -norm, T V and T V 2 behavior (µ fixed ) 400 iterations -Example 2D (Butterfly).

(a) L 2 Fig. 16 .

 216 Fig. 16. ⁄ = 7, µ = 5 -initialization u0 = 0 and v0 = u d , 400 iterations

Fig. 17 .

 17 Fig. 17. Cost functional, L 2 -norm, T V and T V 2 for ⁄ = µ = 0.5, 1, 2, 3 • • • 25 -Dotted (blue) line is initialization (a) and solid (red) line is initialization (b) -800 iterations -Error is given by[START_REF] Bredies | Properties of l 1 -tgv 2 : the one-dimensional case[END_REF]. In this case the cartoon part tends to a constant function as ⁄ = µ increases (the total variation tends to 0) while the second order total variation of the smooth part converges to a limit (not equal to 0, a priori) (Butterfly test)

Fig. 18 .

 18 Fig. 18. T V and T V 2 for ⁄ = µ = 0.5, 1, 2, 3 • • • 25 -Dotted (blue) line is initialization (a) and solid (red) line is initialization (b) -800 iterations -Wall test -The behavior of T V and T V 2 is similar to the previous example.

( a )Fig. 19 . 20 Fig. 20 .

 a192020 Fig. 19. Generic L 2 -norm, T V and T V 2 behavior -400 iterations -Example 2D (Butterfly). Surfaces giving L 2 -norm, T V and T V 2 respectively with respect to ⁄ and µ.

Fig. 21 .

 21 Fig. 21. Test 1D without noise (1000 points)

Table 1 .

 1 Comparison

		716468e-04 6.577422e-03 6.381758e-04 : 8.619796e-05
	Random	1.716468e-04 6.577422e-03 6.381759e-04	8.619796e-05

of di erent initializations -1D case with noise -itmax=50 000 -The stationary pair is denoted (u ú , v ú ) and w

Table 2 .

 2 Comparison of di erent initializations (Butterfly )-The number of outer iterations is limited to itmax=400 -The stationary pair is denoted (u ú , v ú ) and w

		68	1.04	12.70	1.04	2.35
	u 0 = u d , v 0 = 0	18.29	1	13.43	0.43	0.73
	u 0 = 0, v 0 = 0	20.36	1.03	12.87	0.69	0.89
	Random	20.39	1.03	12.88	0.69	0.87
			= 2, µ = 0.1
	u 0 = 0, v 0 = u d	1.5414	2.24 e-01 3.64 e-04	15.15	8.48 e-03	22
	u 0 = u d , v 0 = 0	8.0239	2.76 e-01	3.31	13.52	4.35
	u 0 = 0, v 0 = 0	3.3335	2.45 e-01	0.92	14.65	3.12
	Random	3.5384	3.18 e-01	1.02	14.62	3.30
			= 5, µ = 7
	u 0 = 0, v 0 = u d	61.7005	4.22	5.71	3.45	1.67
	u 0 = u d , v 0 = 0	62.7803	4.02	7.29	2.60	3.25
	u 0 = 0, v 0 = 0	61.6248	4.15	6.34	3.04	1.50
	Random	61.6331	4.15	6.35	3.04	1.55
			= 7, µ = 7
	u 0 = 0, v 0 = u d	69.6775	5.23	2.29	5.69	9.76 e-01
	u 0 = u d , v 0 = 0	72.6262	4.96	4.09	4.51	4.74
	u 0 = 0, v 0 = 0	70.2957	5.13	2.97	5.19	2.45
	Random	70.3114	5.12	2.98	5.18	2.52
			= 7, µ = 9
	u 0 = 0, v 0 = u d	79.7064	5.42	4.10	4.03	1.33
	u 0 = u d , v 0 = 0	80.1229	5.18	5.40	3.33	4.09
	u 0 = 0, v 0 = 0	79.8224	5.33	4.58	3.72	1.86
	Random	79.8297	5.33	4.58	3.72	1.89
			= 10, µ = 15
	u 0 = 0, v 0 = u d	116.2130	7.04	3.59	3.69	1.39
	u 0 = u d , v 0 = 0	116.9598	6.79	4.33	3.36	4.67
	u 0 = 0, v 0 = 0	116.0822	6.95	3.83	3.56	2.02
	Random	116.0918	6.95	3.84	3.56	2.10

Table 3 .

 3 Cost functional and relative error (log scale) for pairs given by initializations (a)

	50	82.38439	81.69328	2.072
	100	80.9713	80.9679	4.379
	200	80.18443	80.35509	2.672
	400	79.83481	79.94497	2.860
	600	79.73564	79.80224	3.078
	800	79.68948	79.73411	3.252
	1000	79.66213	79.69571	3.375
	1200	79.64396	79.67121	3.466
	1500	79.62567	79.64659	3.580
	5000	79.5738	79.5718	4.618

) . # it. F ,µ (ua, va) F ,µ (u b , v b )log 10 | F ,µ |

Table 5 .

 5 Sensitivity with respect to number of iterations : cost functional value. On can refer to Table3as well. Figures 8-10 illustrate the generic behavior of the cost-functional F ⁄,µ .

			10
	50	39.132	25.513	5 e-01
	100	31.727	23.113	3.7 e-01
	200	26.907	21.440	2.5 e-01
	400	23.711	20.377	1.6 e-01
	600	22.410	19.978	1.2 e-01
	800	21.688	19.774	9.6 e-02
		= 10, µ = 15
	50	119.448	117.102	2 e-02
	100	117.578	116.601	8.3 e-03
	200	116.612	116.257	3 e-03
	400	116.215	116.083	1.1 e-03
	600	116.106	116.031	6.5 e-04
	800	116.052	116.006	4 e-04
		= 10, µ = 2
	50	25.90989	39.586	5.2 e-01
	100	25.91003	33.501	2.9 e-01
	200	25.91008	29.512	1.4 e-01
	400	25.91009	27.558	6.3 e-02
	600	25.91009	26.986	4.1 e-02
	800	25.91009	26.699	3 e-02

2

2 One can look at http://maitinebergounioux.net/PagePro/Movies.html to see the convergence process.

Table 6 .

 6 Sensitivity with respect to number of iterations. Here Ï = u b ≠ ua = v b ≠ va and the error is given by the stopping criterion (13) of Algorithm. On can refer to Table4as well.

						10
	50	10.73	12.12	3.98	2.74	1.27	2.19	14.58	6.51
		11.81	12.51	3.13	1.92	1	1.50	8.30	3.33
		12.39	12.74	2.47	1.39	0.81	1.09	4.39	1.72
		12.70	12.87	1.96	1.04	0.69	0.82	2.37	0. 90
		12.80	12.92	1.70	0.9	0.65	0.70	1.53	0.58
		12.85	12.94	1.53	0.83	0.63	0.63	1.09	0.42
					= 10, µ = 15
	50	2.59	4.83	2.74	4.42	3.09	3.14	9.04	9.33
		3.19	4.44	1.85	3.98	3.27	2.14	5.29	6.17
		3.48	4.07	1.18	3.77	3.45	1.39	2.49	3.77
		3.59	3.83	0.73	3.69	3.56	0.87	1.39	2.03
		3.61	3.76	0.55	3.67	3.60	0.67	1.08	1.40
		3.62	3.72	0.45	3.67	3.62	0.54	0.87	1.07
					= 10, µ = 2
	50	8.93 e-03	1.61	1.61	11.21	9.88	1.80	1.1 e-02	18.18
		8.95 e-03	0.89	0.89	11.21	10.56	9.6 e-01	4.7 e-03	12.92
		8.95 e-03	0.41	0.41	11.21	11.01	3.3 e-01	1.9 e-03	8.92
		8.95 e-03	0.18	0.18	11.21	11.18	6.0 e-02	9.3 e-04	4.48
		8.95 e-03	0.11	0.11	11.21	11.20	1.5 e-02	6.2 e-04	2.54
		8.95 e-03	0.08	0.08	11.21	11.21	6.4 e-03	4.2 e-04	1.77

Table 7 .

 7 3 points 10 3 points 10 3 points 10 3 points 10 4 points 10 4 points 10 4 points 10 4 points 10 5 points 10 5 points 10 5 points 10 Test

	5 points

Table 9 .

 9 Cost functional, L 2 -norm, T V and T V 2 for ⁄

	= µ 0.5	F , 6.3577	kwak 2 1.459 e-03	T V (ua) 4.685	T V 2(va) 7.646	Error 4.13 e-01
	1	12.2448	2.655 e-03	4.756	6.853	5.19 e-01
	5	52.4043	9.654 e-03	3.020	5.782	5.83 e-01
	10	93.2718	1.562 e-02	1.733	5.395	5.59 e-01
	13	114.9198 1.835 e-02	1.268	5.237	5.11 e-01
	17	141.3794 2.128 e-02 8.506 e-01	5.066	4.78 e-01
	21	165.9126 2.357 e-02 5.768 e-01	4.940	4.36 e-01
	25	188.8569 2.537 e-02 3.941 e-01	4.840	3.86 e-01

Complete results (text files, movies, other examples) and MATLAB © code, are available at http://maitinebergounioux.net/PagePro/Movies.html.

We now investigate the sensitivity of the model with respect to quantification. Let u d a data (with values in [0, 255] for example). Let (⁄, µ) be chosen parameters and (u ⁄,µ , v ⁄,µ ) the corresponding computed pair (with the appropriate initialization). Let -> 0 and consider the new data -u d . This is the case, for example, if we get 16 bits images and convert them to 8 bits : in this case -= (2 8 ≠ 1)/(2 16 ≠ 1). We may want to normalize the data as well: in this case -= 1/ max(u d ). The question is to know what new parameters ( ⁄, μ) must be chosen to get u ⁄,μ = -u ⁄,µ and v ⁄,μ = -v ⁄,µ . For any (u ⁄,µ , v ⁄,µ ) solution to (P ⁄,µ ), we get


M Bentoumi
email: bentoumi@esstin.uhp-nancy.fr

P Aknin
email: aknin@inrets.fr

G Bloch
email: bloch@esstin.uhp-nancy.fr

On-line rail defect diagnosis with differential eddy current probes and specific detection processing

Keywords: detection -differential probe -inverse filtering -Wiener filter -eddy currents -rail defect -wavelet PACS. 06.90.+v Other topics in metrology, measurements, and laboratory procedures

This paper presents a non destructive evaluation system for rail inspection with a device embedded in a subway train. The difficult operating conditions lead to the design of a specific double-coils and double-frequencies sensor that delivers 8 complementary signals. After the presentation of the technological solution, the paper deals with the detection of minor defects and presents three approaches: the first one is a time heuristic approach, the second one is based on wavelet basis projections, and the last one implements an inverse filtering which allows the transformation of the differential probe outputs into estimated absolute ones. The methods are validated statistically on a labeled defect data set with the help of ROC curves.

The context of the application

The introduction of automatic driverless trains involves the redefinition of the whole traffic control system. Currently, the basic device of traffic control is the track circuit whose principal function is the train detection along a track section and secondary function is the detection of broken rail. This device will disappear in an automatic control context, and must be replaced by new means of broken rail detection. The frequency of rail breakage is weak, but its occurrence degrades immediately the quality of the complete transport system, of mechanical and electrical points of view (traction currents return to the electric substation through the rails). This paper introduces a suitable differential eddy current sensor planned to replace the track circuit in its function of broken rail detection. An eddy current sensor is sensitive to any modification of the geometry and/or electromagnetic characteristics of the target. Obviously, transverse splits of the rail are detected, but minor defects as well (shelling, welded joints or corrugation for example). The detection and monitoring of these kinds of defect allow to set up predictive maintenance policy. This paper is focused on the detection procedures for minor defects.

The instrumentation and signal processing chain

Previous works led to the realization of a suitable doublecoils and double-frequencies sensor [START_REF] Oukhellou | Parametrisation et classification de signaux en contr le non destructif. Application la reconnaissance des defauts de rails par courants de Foucault[END_REF]. The sensor structure was designed and optimized according to the following specifications: positioning at 40 mm height, vertical and horizontal displacements of the sensor due to the bogie dynamics and 100 km/h maximum speed of the train. Moreover a particular attention must be focused on strong acceleration levels (until 10g) and electromagnetic compatibility problem caused mainly by the traction currents that circulate in the rails.

Two differential coils are used to ensure a greater immunity to the target-sensor distance variations for the detection of localized cracks as well as large shellings. The two supply frequencies make the sensor sensitive to the rail electromagnetic characteristics for different skin depths. After preprocessing, four complex channels (active and reactive parts) are available, which are equivalent to eight real signals d 1 (t), d 2 (t), ..., d 8 (t) (figure 1). The device thus generates complementary information that will be aggregated during the phase of defect recognition.

The digitalization is carried out with respect to the distance covered by the train with a fixed step dx of 5 mm. It must therefore be noted that, in all that follows, processed signals are discrete signals depending on distance After perusal of records, track visits are required, in order to precisely label each particular detected event. A series of tests led to gather about 600 defects distributed into four classes.

For on-line diagnosis, a detection (or segmentation) stage is required to isolate particular data segments on which a pattern recognition procedure can be applied.

The complexity of the processing lies in the fact that the track presents a great number of singularities which must be distinguished from the real defects. Furthermore, processing must be compatible with real time calculation capabilities of the system. Three approaches for detection of minor defects (shellings, Sh, and welded joints, WJ) are presented. The first one is a time heuristic approach, the second one is based on wavelet basis projections, and the last one implements an inverse filtering which allows the transformation of the differential probe outputs into estimated absolute ones. The methods are validated statistically on the labeled defect data set with the help of ROC curves.

3 Heuristic time detection

Isochronism of differential signatures

As indicated previously, the sensor includes differential coils. A defect signature in the signals is directly related to an imbalance between the reluctances of the two associated windings. In the case of localized defect, the complete differential signature takes a symmetrical form with two successive swings, positive and negative (figure 3b). The equivalent frequency is proportional to the inner-distance of the differential coils.

Due to the bogie dynamics, the main source of errors lies in the rotating movements of the sensor relatively to the target (pitching). In this case, the reluctances of the elementary windings vary simultaneously (figure 3a). The signature is also symmetrical, but the equivalent frequencies are the same for the two differential coils, depending only on the natural frequency of the bogie and the speed of the train.

One solution to visualize the isochronism of the two signatures is to plot each signature with respect to one of the others into a pseudo-Lissajous plane (figure 3). In the case of pitching, the curve is roughly linear, whereas for a real defect the curve has two loops indicating the sequential run of the defect under the different windings. A detector can be defined in such a plane, using linear decision boundaries, with threshold q and slope p (figure 3). As the sensor has four active channels with two different inner-distances, two distinct detectors can be built in planes {d 1 , d 5 } and {d 3 , d 7 }. The fusion of these detectors is operated in a conjunctive way:

[(d 5 < pd 1 -q) ∩ (d 5 > 0)] ∩[(d 7 < pd 3 -q) ∩ (d 7 > 0)] = 1 =⇒ minor def ect (1)

The ROC Curve

This detection requires to tune the two parameters q and p of the boundary. The Receiver Operating Characteristic (ROC) was chosen to tune the detector. Well detection probabilities (Pw) are drawn according to the false detection probabilities (Pf) and the (0,1) point is the optimal operating point. These probabilities are estimated with the help of the experimental database built during in-site testing (435 minor defects). The nearest to the (0,1) point the ROC curve is, the better the tuning is. The optimal tuning is obtained with p=1 and q=0.008. However, the minimal distance to the (0,1) point remains relatively large (dist≈0.39). Note that this detector is implemented after a bandpass filtering between 0.05 m and 2.5 m.

4 Detection in the time-scale plane

Choice of the mother wavelet and optimum scales

The second detector transposes the decision procedure in a transformed space. The transformation used is the continuous wavelet transform [START_REF] Strang | Wavelets and filter banks[END_REF], applied to a complex spacedependent signal f (x):

W [f ](x, s) = +∞ -∞ f (t) 1 √ s Ψ * ( t -x s )dt (2)
where

1 √ s Ψ ( t-x s
) is a version dilated with a scale parameter s and translated by a space-shift x of a function Ψ (t) called "mother wavelet", such as +∞ -∞ Ψ (t)dt = 0, and Ψ * stands for the complex conjugate of Ψ . Four continuous wavelet transforms are applied to the four eddy current channel outputs d 2i-1 (x) + jd 2i (x), i = 1, . . . , 4:

W i (x, s) = W [d 2i-1 (x) + jd 2i (x)](x, s), i = 1, . . . 4 (3)
The choice of the mother wavelet is important, particularly when the wavelet coefficients W i (x, s) are used as detectors. From a certain viewpoint, this kind of detection is close to adapted filtering [START_REF] Joubert | [END_REF]. The statistical defect base allows us to define a typical form of defect signature. Eight typical signatures (corresponding to the eight real signals) for each defect have been built by filtering and averaging the real defect signatures. Figure 5 shows two typical signatures of minor defects corresponding to the first channel (d 1 (t), d 2 (t)). A spline biorthogonal wavelet [START_REF] Strang | Wavelets and filter banks[END_REF] has been chosen, its form being close to the differential minor defect signatures (figure 5). Once the mother wavelet is chosen, the optimum scales s i , i = 1, . . . , 4, must be determined. The modulus of the first wavelet transform, applied on minor defect typical signature, W 1 (x, s), corresponding to the first channel, is presented in figure 5. The higher the modulus is, the better the detection is. In this example, the highest modulus is obtained for scale value s 1 = 27. the optimal scale values for the four channels have been found equal to [27, 27, 39, 39]. The optimum scale mainly depends on the dimension of the differential coil arrangement (inner-distance between elementary coil).

Fusion of detectors

Each detector is dedicated to one coil and one frequency. Their redundancy can improve the robustness of the detection by aggregate them. For this purpose, the following total detector was built :

D(x) = 4 |W i (x, s i )| > T h ⇒ minor def ect (4)
where T h is the detection threshold. Figure 6 presents an example of wavelet transform for the four channels and their multiplicative fusion.

-0 100 200 300 400 500

Th Π|W i | |W 4 | |W 3 | |W 2 | |W 1 |
x (m) Fig. 6. fusion of wavelet detectors.

This multiplicative fusion increases the contrast between defect and no defect decisions. It can be referred to the evidence theory where belief masses are aggregated using the Dempster-Shaffer rule [START_REF] Aknin | Systemes multi-capteurs: diagnostic et fusion[END_REF]. In the case of two exclusive classes (defect or no defect), nonstandard probabilistic fusion makes a normalized multiplication of the belief masses; here the normalization coefficient is equal to 1, indicating that no conflict exists between the detectors.

Inverse filtering

Differential probes allow relevant measurements in noncooperative contexts such as modifications of target-sensor distance, variations of temperature, electromagnetic perturbations... However, the differential mode modifies the signal shape and its interpretation becomes more difficult.

The third detection procedure is built from estimated absolute signals, outcome of a deconvolution phase. First, the direct problem is introduced for differential measurements.

Direct model of single differential measurements

The principle of differential measurement is to subtract two close absolute measurements in order to constitute a local approximation of the absolute signal derivative. A compromise must be made between the quality of the derivative and the amplitude of the differential signal, both related to the inner-distance of the absolute coils (the smaller this distance is, the better the derivative quality is and the lower the differential signal amplitude is). The differential signal d is given by :

d(x) = a(x) -a(x -L) ( 5
)
where a is the absolute signal and x corresponds to a spatial variable. L is the inner-distance between the two absolute coils. An example of differential response is shown figure 8. The discrete version of equation ( 5) is:

d(i) = a(i) + a(i -) (6)
where = L/dx and dx is the sampling width. Equation (6) corresponds to a linear filtering with an impulse response h containing two opposite Dirac impulses separated by points (cf. figure 8).

In the general case of additive noise, the estimated absolute signal â is given by : â = g * (h * a) + g * b [START_REF] Ragot | Validation de donnees et diagnostic[END_REF] where * denotes the convolution operation, a and b are respectively the absolute signal and the noise, h the impulse response of the differentiation and g the impulse response of the inverse filter. The corresponding block diagram is given figure 9. The case of double-differential measurements can be taken up with the same approach [START_REF] Bentoumi | Advances in Signal Processing for Non Destructive Evaluation of Materials[END_REF].

Wiener filtering

Using for G(jω) a direct inverse filter 1/H(jω) induces some problems. Indeed, a differential measurement cancels the harmonic components 1/L, 2/L, 3/L... So the pure inverse filter has an infinite impulse response and becomes unusable in practice.

The inverse Wiener filter G minimizes the mean-square error between a(x) and â(x), for a stationary context and uncorrelated signals. Its expression is based on Wiener-Hopf equation [START_REF] Hayes | Statistical digital signal processing and modeling[END_REF].

In the frequency domain, the signals a(x) (absolute signal), b(x) (noise) and d(x) (differential signal) (cf. figure 9), of Fourier transforms A(jω), B(jω), D(jω), have respectively P a (jω), P b (jω), P d (jω) as power spectral densities, and the exact formulation of G is:

G(jω) = P ad (jω)/P d (jω) (8)
where P ad (jω) is the cross-power spectral density between absolute and differential signals. According to the variable definitions,

D(jω) = H(jω)A(jω) + B(jω) (9)
the power spectral density of the output is:

P d (jω) = P a (jω) HH * (jω) + P b (jω) (10)
Assuming that b(x) and a(x) are uncorrelated, the crosspower spectral density is:

P ad (jω) = P a (jω)H * (jω) (11)
By using equations ( 8), (10) and (11), the inverse filter becomes:

G(jω) = 1 H(jω) HH * (jω) |H(jω)| 2 + P b (jω)/P a (jω) (12)
Therefore, the optimal inverse filter is the multiplication of the direct inverse filter by a regularization function depending on the noise-to-signal ratio N S 2 r = P b (jω)/P a (jω).

If N S 2 r tends to 0, the inverse filter tends to the direct inverse filter; if N S 2 r tends to +∞, the inverse filter tends to 0. Obviously the noise-to-signal ratio can be adjusted along the different frequencies, on a priori noise characteristics.

Examples of inverse filter realizations

The inverse Wiener filter is an infinite impulse response (IIR) filter, tuned with the N S r ratio. Figure 10 The inverse filter tuning requires the knowledge of the noise characteristics (white noise or colored noise, signalto-noise ratio). Once this adjustment is carried out, the output signal can be considered as an estimation of the absolute signal. Figure 11 presents processed signatures. Note that the estimated signal is similar to an absolute one.

0.1 0 -0.1 -0.2 0 1 2 3 4 5
d 1 (x) â 1 (x)
x (m) Fig. 11. two minor defects and their inversions.

Figure 12 presents the frequential behavior of direct and inverse filters for three values of N S r ratio. The inverse filter G(jω) tends to the direct inverse filter 1/H(jω) when the N S r ratio becomes low.

Detector structure

On-line pattern recognition procedure requires two successive operations: detection and classification. The first solution consists in a detector deliberately permissively tuned [START_REF] Ragot | Validation de donnees et diagnostic[END_REF], followed by a complete classifier including a rejection class. Following this approach, previous works have been developed with the help of neural network techniques [START_REF] Oukhellou | [END_REF].

The present results are based on a single processing which aggregates detection and classification phases and is dedicated to the minor defects (shellings and welded joints). The detector structure associates four elementary detectors. Each of them processing a particular inverted signal (figure 13). These detectors are paired on in order to increase the detection contrast. Then, the final decision phase is carried out in a particular plane whose coordinates are the enhanced detectors, according to the existing expert knowledge.

The detection is validated when the curve comes into the gray zones of the pseudo-Lissajous plane (figure 14): The adjustment of the thresholds λ and µ was made with the help of ROC curve (figure 15). The optimal setting stands at the point where the distance between the ROC curve and the optimal point (0,1) is the lowest (dist ≈ 0.3). The results are similar to the wavelet detector ones. Nevertheless, contrary to wavelet detector, the enhanced Wiener detector allows the separation of shellings from welded joints: the first quarter of the plane is dedicated to the shelling detection whereas the third is dedicated to the welded joint detection. Moreover, the distributions of false detection probability and well detection probability are homogenous on the whole minor defect class (union of shelling class and welded joint class). This is due to the fact that each part of the detection area has its own threshold µ or λ.

[(â 1 â2 > λ) ∩ (â 5 â6 > 0)] ∪[(â 1 â2 > 0) ∩ (â 5 â6 > 4λ)] ∪[(â 1 â2 < -µ) ∩ (â 5 â6 < 0)] ∪[(â 1 â2 < 0) ∩ (â 5 â6 < -4µ)] = 1 =⇒ minor def ect (13)

Conclusion

This paper presented an original non destructive evaluation system dedicated to the rail inspection. After the presentation of the technological solution, the paper dealt with the detection of minor defects (shellings and welded joints) and presented three approaches: the first one is a time heuristic approach, using directly the differential measurements, the second one is based on wavelet decomposition, and the last one implements an inverse filtering which allows the transformation of the differential probe outputs into estimated absolute ones.

The wavelet approach gives best results, thanks to the fusion of four detectors. But this detection procedure penalizes the shelling class to the detriment of the welded joint class. On the contrary, the inverse Wiener filtering allows to define a more flexible detector structure which is able to separate shellings and welded joints. However, its practical realization requires some a priori knowledge about noise model. This knowledge is generally not easy to collect, but once it was carried out, the detection results are improved.

Fig. 1 .FJFig. 2 .

 12 Fig. 1. data processing chain.

Figure 2

 2 Figure 2 gives an example of the active part d 1 (x) on 500 m of rail. Particular points can be located, like switch (Sw), undulatory wear (UW) of the rail, fishplated joints (FJ), welded joints (WJ) and shellings (Sh).

Fig. 3 .

 3 Fig. 3. effects of pitching (a) and defect (b).

FalseFig. 4 .

 4 Fig. 4. ROC curves for heuristic detector.

Figure 4

 4 Figure 4 presents ROC curves for 0.3 <p< 10 and 0.001<q<0.1. The optimal tuning is obtained with p=1 and q=0.008. However, the minimal distance to the (0,1) point remains relatively large (dist≈0.39). Note that this detector is implemented after a bandpass filtering between 0.05 m and 2.5 m.

7 |W 1 Fig. 5 .

 715 Fig. 5. mother wavelet, minor defect typical signature and wavelet transform.

Fig. 7 .

 7 Fig. 7. ROC curve of the time-scale detector.

Figure 7

 7 Figure 7 illustrates the tuning of the detector threshold with the help of ROC curve. The optimal value T h=1.5e-4 corresponds to the minimal distance dist ≈ 0.29.

Fig. 8 .

 8 Fig. 8. impulse response of single differential measurements with =6.

Fig. 9 .

 9 Fig. 9. principle of inversion.

Fig. 10 .

 10 Fig. 10. example of impulse response of one inverse Wiener filter.

Fig. 12 .

 12 Fig. 12. frequency responses in dB for different tunings of differential measurement inversion.

6 Fig. 13 .

 613 Fig. 13. enhanced detectors.

6 Fig. 14 .

 614 Fig. 14. detector representation.

Fig. 15 .

 15 Fig. 15. ROC with enhanced Wiener inversion.

This work was supported by the French PREDIT organization (Interministerial land transport research and innovation program) with the coordination of RATP Company.


Tieyong Zeng
email: zeng@math.univ-paris13.fr

INCORPORATE T V -L ∞ MODEL WITH SPARSE REPRESENTATIONS FOR IMAGE DENOISING, A POST-PROCESSING APPROACH

Keywords: Sparse representation, total variation, image denoising, dictionary, K-SVD, post-processing

Sparse representations of images have revoked remarkable interest recently. The assumption that natural images admit a sparse decomposition over a redundant dictionary leads to efficient algorithm for image processing. In particular, the K-SVD method has been recently proposed and shown to perform very well for gray-scale and color image denoising task ([1],[2]). Meanwhile, the T V -l ∞ model with special choice of dictionary has been proved to be very effective for image restoration ([3],[4]). In this paper, we propose a hybrid model which combines these two methods and may be regarded as a post-processing procedure for K-SVD. Due to the excellent work of K-SVD and the fact that T V -l ∞ can reconstruct lost information quickly, this hybrid model lead to a new state-ofart denoising performance.

INTRODUCTION

In this paper, we consider the classic image denoising problem: an ideal image x is observed in the presence of an additive zero-mean Gaussian white noise of standard deviation σ. Thus the observed image y is:

y = x + b. (1)
Our task is to design an algorithm to remove the noise from y, getting as close as possible to the original image x and in the mean time, getting a better visual quality if that is possible.

Recently, Michael Elad and Michal Aharon proposed an image denoising method via sparse and redundant representations over learned dictionaries ([1]). This leads to state-of-the art denoising performance. The importance of this method is that it can recover most of the information in the noisy image while there is very little wash-out effect(especially for the face region of image Barbara). That's to say, it is able to avoid the shortcoming of most total variation (T V ) based denoising method.

The little drawback of this approach is that it has some so called "checkboard" effect and it sometimes still lost some texture information. It is well known that total variation model can avoid the checkboard effect and TV-Gabor( [START_REF] Zeng | Using gabor dictionaries in a tv -l ∞ model, for denoising[END_REF]) has proved to be very effective for texture restoration. So in this paper, we try to incorporate total variation concept into K-SVD model. This will lead to a new state-of-art denosing performance.

K-SVD Denoising Algorithm

In this section, we will briefly review the main mathematical framework of the approach of Elad([1]), as this is one of starting point of our work. First let the clean image x be written as a column vector of length N . Given fixed-size pathes √ n × √ n, we assume that all such patches in the clean image x admit a sparse representation. Addressing the denoising problem as a sparse decomposition technique for each patch leads to the following energy minimization problem:

{ α i,j , D, x} = arg min D,αij ,x λ x -y 2 2 + i,j µ i,j α i,j 0 + i,j Dα i,j -R i,j x 2 2 . (2)
In this equation, x is the estimator of x, and the dictionary D ∈ R n×K is an estimator of the best dictionary which gives the sparsest representation of the patches in the restoration image. The indices [i, j] mark the position of the patch in the image(representing it's top-left corner). The binary matrix R i,j of n × N extracts the square √ n × √ n patch of coordinates [i, j] from the image represented by a column vector of size N . µ i,j is hidden parameter which is decided implicitly by the method of solving the above equation. The first term of the above equation demands a proximity between x and y. The second and the third term both pose the image prior. This regularization term assumes that goodbehaved natural images are to have a sparse representation for every patch, and from every position in the image, over the learned dictionary D. The second term ensures the sparsest representation, and the third term forces the consistency of the decomposition.

T V -l ∞ Algorithm for Denoising

Our second starting point is the works of ( [START_REF] Zeng | Using gabor dictionaries in a tv -l ∞ model, for denoising[END_REF]). They investigate the solution provided by the following model:

minimize T V (x) under the constraint x -y D0,∞ ≤ τ (3)
where . D0,∞ is defined by

x D0,∞ = sup ψ∈D0 | x, ψ |,
for a finite dictionary D 0 ⊂ R N and the choice of τ depends on the noise level σ. For an explicitly definition of discrete definition of T V of image x, please see [START_REF] Zeng | Using gabor dictionaries in a tv -l ∞ model, for denoising[END_REF] and references cited therein. Typical choice of D 0 is but not limited to wavelet basis, wavelet packets, translate-invariant Gabor dictionary or translate of small special patches. Using the idea of penalization, ( 3) is solved by minimizing the unconstrained energy

x = arg min x T V (x) + γ ψ∈D0 ϕ τ ( x -y, ψ ), (4)
with

ϕ τ (t) = (sup(|t| -τ, 0)) 2 ,
and for a large number γ. This model assumes the original image x has small total variation and forces the residual x -y to be smaller than τ in every special direction ψ ∈ D 0 . This leads to an efficient restoration method which preserves edges and texture very well and has little checkboard effect. Now before presentation of our new hybrid model, let's do some comparisons for these two models. Comparing to T V -l ∞ model, K-SVD model gives higher P SN R and has very few washout effect, especially for the face region of Barbara; comparing to K-SVD model, T V -l ∞ model has more chance of avoiding the checkboard effect. Both methods can recover most part of the texture and sometime they can be regarded as complementary. And Do NOT forget that both methods use energy-minimization idea. So why not combine these two?! This will lead to our new approach. Meanwhile, we should point out that incorporate T V into sparse representation is not a totally new idea. In fact, the authors of [START_REF] Starck | Image decomposition via the combination of sparse representation and a variational approach[END_REF] have considered similar model together with special choosing dictionaries which leads to an efficient image decomposition method.

TV-SPARSEST HYBRID MODEL

From the above discussion, we now propose our new approach. In order to solve the denoising task, we try to use the following new energy minimization problem:

{ α i,j , D, x} = arg min D,αij ,x λ x -y 2 2 + i,j µ i,j α i,j 0 + i,j Dα i,j -R i,j x 2 2 + βT V (x) + γ ψ∈D0 ϕ τ ( x -y, ψ ).
(5) Now in the new Eq.5, the first three terms are just the same as the Eq.2 and the last two are from Eq.3. The term T V (x) forces the estimator image x to be total variation limited and ψ∈D0 ϕ τ ( x -y, ψ ) makes the energy of residual x -y on the special interesting direction no higher than τ . With all these constraints, better result is expected.

NUMERICAL ASPECTS

To approximate a solution for Eq.5, we first use the K-SVD procedure of ( [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]) to minimize (2) and later we use another T V iteration to minimize the last two terms. Basically, this algorithm can also be regarded as initialization with K-SVD result and then T V -iterations as post-processing procedure. We present the details in Table 1.

In the above denoising procedure, follow with ([1]), we assume that σ is known and then we set J = 10, λ = 30/σ, C = 1.15 (another choice is 0.93 which is based on Rayleigh law, see ( [START_REF] Mairal | Sparse representation for color image restoration[END_REF])). As our algorithm only depends on γ/β,we can fix β as 1.0. Being a penalization parameter,the typical choice of γ is between 9999 and 999999. τ = 3.5σ is empirically from ( [START_REF] Zeng | Using gabor dictionaries in a tv -l ∞ model, for denoising[END_REF]). We discuss a little about the choice of k.

After K-SVD procedure, we can already get a fairly good estimator x. Now if we aim to obtain an estimator x which is as near as possible to the original image, i.e. we need a higher P SN R,we can stop the T V iteration when:

• T V (x) reaches its top This is not strange, when the estimator is near the original image, the higher the T V , the higher possibility of presence of recovered structures, then the higher the P SN R. Typical choice is k = 2 or 3.

On the other side, if we want better visual quality (less artifact, better texture, better structure), we can iterate more (but not too much else this will cause wash out effect)to get a smaller T V . But in the same time, we should note that this will surely make P SN R a slightly lower. Our tests show that this time k = 15 is a good choice.

EXPERIMENTS RESULTS

We report our experiments for σ = 20, 30 on Barbara image. As the same in [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], we assume that σ is already known or could be estimated from elsewhere. We use Gabor dictionary (Gabor II of ( [START_REF] Zeng | Using gabor dictionaries in a tv -l ∞ model, for denoising[END_REF]), large size), as we think that this one can recover better high-frequency information which has higher possibility than low-frequency to have been lost during the K-SVD procedure. The sum of FFT of the 145 filters in this dictionary is shown as Fig. 1.

Task: Denoise a given image y from Gaussian noise of standard deviation σ via solve Equation [START_REF] Starck | Image decomposition via the combination of sparse representation and a variational approach[END_REF]. Parameters: n-block size, K-size of first dictionary, J-number of K-SVD iterations, λ-Lagrange multiplier, C-noise gain, β-penalization parameter for T V , γpenalization parameter for special direction of residual x -y, k-number of T V iterations.

1. Set x = y, D=overcomplete DCT dictionary.

Repeat J times:

• Sparse Coding Stage:

Use OMP pursuit algorithm to compute the representation vectors α i,j for every fixed patch R i,j x, by approximating the solution of

min αi,j α i,j 0 subject to R ij x-Dα i,j 2 2 ≤ (Cσ) 2 .
• Dictionary Update Stage: For each column l = 1, 2, . . . , k in D, update it by -Find the set of patches that use this atom, ω l = {(i, j)α i,j (l) = 0} -Find every index (i, j) ∈ ω l , compute its representation error

e l i,j = R i,j xi,j - m =l d m α i,j (m)
-Set E l as the matrix whose columns are {e l i,j } (i,j)∈ω l -Apply SVD decomposition E l = U △V T .

Choose the updated dictionary column dl to be the first column of U . Update the coefficient values {α i,j (l)} i,j ∈ ω l to be entries of V multiplied by △(1, 1).

3. Set:

x = (λI + i,j R T i,j R i,j ) -1 (λy + i,j R T i,j Dα i,j )
4. Repeat k times:

• Calculate direction of gradient:

w = div( ∇x |∇x| ) - γ β ψ∈D0 ϕ ′ τ ( x -y, ψ )ψ
• Find the optimal step :

t = arg min s∈R βT V (x)+γ ψ∈D0 ϕ τ ( x -y, ψ ) • Update x: x = x + tw.
Table 1. General form of the algorithms. σ is known.

Step 1 to 4 is typical K-SVD procedure( [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]) to minimize (2). This algorithms can be regarded as initialization with K-SVD result and then a T V -iteration as post-processing. where k is the T V iteration number of Table 1. Note that k = 0 is result of K-SVD of Elad. .

Noise level of σ = 20.

Our first experiment is to denoise the Barbara image with noise level σ = 20, the P SN R of noisy image is 22.0977.

Fig. 2 shows the relationship between T V (x) and P SN R(x) with the iteration number k in the main algorithm (see step 4 of Table 1). This Fig. tells us that after about 2 or 3 T V iteration, T V (x) and P SN R(x) reach their tops at the same time. So if we aim to get a higher P SN R denoising result, we should stop at the iteration where T V (x) reaches its top.

For out experiment, the highest P SN R is 30.9376 of our new approach when k = 2, this is a slightly higher than K-SVD of Elad (30.8113) which claimed state-of-art denoising performance and much higher than the classical Rudin-Osher-Fatemi method (24.6759).

Meanwhile, if we want to a better visual quality, we can continue the T V iteration and k = 10 to 15 is a typical choice. Fig. 3 displays a slice of left-bottom part of the Barbara with k = 15 from which we clearly know our approach is more effective.

Noise level of σ = 30.

In our experiment, when σ = 30 the P SN R of the slice is 18.5448. Fig. 4 shows the result for the same slice as Fig. 3. From this Fig. we obviously see that our new approach performs best. Both Rudin-Osher-Fatemi method and K-SVD of Elad fail to recover the texture of tablecloths on the desk, while our approach still can recover most of the information. And for the left part this slice, Rudin-Osher-Fatemi lost the texture information and K-SVD can recover some of this information but the texture is a little disordered. But for our new approach, this is a different story, we have recovered more information and the texture is still arranged.

Globally,for this level of noise,the P SN R of noisy image and the result of Rudin-Osher-Fatemi,K-SVD,and our approach are respectively 18.5867, 24.0429, 28.5947,and 28.8376.

CONCLUSION

In this paper, we have proposed a hybrid model which combines the sparse representation and total variation concept. In our current methodology of solving the model, this can also be regarded as initialization T V -l ∞ model with K-SVD or as a post-processing for K-SVD. Solving the hybrid model directly or choosing D 0 adaptively to noisy image may lead to better result and thus could be future work. Due to the leading performance of K-SVD and the fact that T V -l ∞ with Gabor dictionary procedure can reconstruct the lost information quickly, our hybrid model attains a new state-of-art denoising performance, equivalent or surpassing recently published leading alternative methods including K-VSD itself.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Sum of Fourier transforms of the 145 filters in the Gabor II dictionary(big size). This image is larger size version of Bottom-Right of Fig.1 in [3].

Fig. 3 .

 3 Fig. 3. Denoising a 128 × 128 slice of Barbara. From Left to right and from top to bottom: Noisy image( σ = 20), PSNR 22.0896; Rudin-Osher-Fatemi, PSNR 24.2663); K-SVD of Elad, PSNR 28.9013; Our new approach,PSNR 29.1148.

Fig. 4 .

 4 Fig. 4. Denoising the same slice of Barbara as Fig.3. From Left to right and from top to bottom: Noisy(σ = 30), PSNR 18.5448;Rudin-Osher-Fatemi, PSNR 23.4331; K-SVD of Elad, PSNR 26.4467;Our new approach,PSNR 27.032.

Acknowledgement

The author would like to thank for Professeur M. Elad for his great help.


Professor Juan R Pimentel

Problem C: Calculation of P merr for FlexCAN

Introduction

The occurrence of certain incidents in the normal operation of in-vehicle networks, such as bit errors caused by electromagnetic interference (EMI), produces a special form of network partitioning (i.e., a node is temporarily isolated from the rest of the network).

Recovering from these situations takes some time, leading to an increase of the message latency (also called response time). This may induce failures of expected hard real-time properties or services of the network. Timeliness of applications may be compromised and may lead to application failures. Thus, methods to analyze the impact of such faults are important.

Several approaches have been proposed to analyze the impact of externally induced faults. The approach proposed in the CANELy architecture [], where this type of network partitioning is called network inaccessibility, requires: I1 -Study the accessibility constraints, ensuring that the number of inaccessibility periods and their duration have a bound. I2 -Show that such a bound is suitably low for service requirements. I3 -Accommodate the effects of inaccessibility events in the timeliness model and in protocol operation, at all the relevant levels of the system.

We classify the CANELy approach as a bottom-up approach because it starts at the protocol level (using CAN). In this study, we propose a top-down analysis approach starting at the application and regressing to lower levels of the application and details of the communication system. Nevertheless there are common issues in both approaches. The characterization of network partitioning including techniques to overcome them is crucial in simplex networks and cannot be ignored in networks with node and channel redundancy. Indeed, network partitioning affecting individual network replicas would lower their fault coverage in the time domain.

In (Broster at al., 2004), a numerical method to evaluate the impact transient random errors (e.g., EMI) on the real-time delivery capability of CAN and TTCAN has been developed. This method has proven particularly useful to evaluate error recovery mechanisms inherent in some protocols such as CAN. The method can be also applied to TDMA-based networks (e.g., TTP/C and FlexRay). Preliminary evaluations in [] has shown that because of the error recovery inherent in CAN, its performance against transient random errors is much better than equivalent time-triggered protocols (e.g., TTCAN) that lack explicit error recovery mechanisms.

The FlexCAN communication architecture is both time-triggered and event triggered and its performance under transient random errors can be also analyzed using this method. In this section, we use the numerical evaluation method of [] to make a detailed study of message latency in the FlexCAN communication architecture under transient random errors. Such studies are important because they are the basis for application oriented safety calculations.

Problem formulation

Given the FlexCAN communication architecture, find the probability of failure delivery of message m i , P merr defined as the probability that the maximum message latency R i corresponding to message m i is greater than the sub-cycle interval T sc , (i.e., R i > T sc ).

The solution to this problem depends on three elements: the communication protocol, the structure of the time-triggered communication cycle, and the error model assumed for the external random disturbances. Details of the communication protocol (e.g., the frame structure, the data rate, etc.) are needed to numerically evaluate R i . However the probability of failure delivery of a message is particularly sensitive to the error recovery policy or mechanism inherent in the protocol. CAN uses the feature of automatic message re-transmission as the error recovery policy; whenever any error is detected, the protocol will automatically retransmit the message in error. The structure of the communication cycle is important because it imposes a time reference, not only for the calculation of R i but also for the application, with important simplifications and improved accuracy. Finally, it should be intuitive obvious that the error model used for the external disturbances is important because it will dictate precisely how R i is affected by errors and the impact of errors in the probability of failure delivery of a message.

Evaluating P merr is one important step in the overall safety evaluation of communication architectures against transient random disturbances. Other steps involve the nature and duration of the external disturbance and the detailed characteristics of the application in question (e.g., steer-by-wire).

Fault Model

The error model is used to characterize external disturbances of abnormal potentially dangerous environments rather that background noise typically found in normal operating environments. The bit error rates in normal operating environments are rather low and the CAN protocol is excellent in dealing with them. The undetected error rate for CAN is [],

Prob. {undetected error rate, normal environments} < 4.7x10 -11 For this reason, most studies consider that the probability of undetected errors in CAN is zero. An example of an abnormal potentially dangerous environment is an section of road crossed by an vehicle near powerful radio frequency transmitters or high voltage lines emitting unusual levels of electromagnetic interference. Several models have been suggested for the probabilistic characterization of these error source ranging from constant [Wilwert, EFTA2005], to Poisson [Broster-Rod, Navet], to generalized distributions [Simonot, FET2005]. Following references [], in this study we also assume the external disturbances follow a random Poisson model that states that the probability of having n errors in a time interval t is given by the following Poisson probability law,

p(n,t) = Prob{n errors in an interval t} = ! ) ( n t e n t λ λ -
Where λ is the error rate.

The FlexCAN Communication Model

The FlexCAN architecture has been developed to support safety oriented applications with features to help meet safety integrity requirements. It is based on a time-triggered communication cycle that is composed of a number of sub-cycles of equal or different length. All messages are scheduled on a sub-cycle basis and for each sub-cycle, messages access the bus according to the CAN protocol. Thus, the error recovery policy of FlexCAN is the same as that of CAN (i.e., it automatically retransmit frames in error) with one simplifying assumption, any message not transmitted at the end of their transmission sub-cycle are flushed from the system (i.e., removed from the CAN transmit buffer). The main reason for this is to respect the schedule of the remaining messages and thus to achieve time-domain composability. An advantage of this policy is that messages from one sub-cycle do not interfere with messages in any-other subcycle. It is important to note that for FlexCAN applications, the number of messages per sub-cycle typically range from 2 to 6 [Padova Lift Truck] The FlexCAN protocol also has node and channel replication features to deal with permanent failures as well as a bus guardian to deal with babbling idiot failures but these are not detailed here.

Summary of CAN message latency calculations

Assume a CAN network with a data rate of R bits/sec, a one bit propagation delay τ b = 1/R, with a certain number of messages m 1 , m 2 , to be transmitted. Each message m i in a multi-rate system is characterized by a tuple (C i , T i , D i , J i , P i ) where:

• C i is the time to transmit the message of size b bytes and including overhead bits, • T i is the transmission period, • D i is the relative deadline defining the maximum time interval tolerated between the transmission request and the reception of the frame at all nodes, • J i is the queueing jitter • P i is the CAN message identifier (ID) defining its priority

The so-called Tindell Equations [] involve a procedure for calculating the worst case message latency (also called message response time) R i , defined as the time interval from when a transmitter node enqueues a message for transmission until such message is successfully read by any receiver node. The procedure summarized below is for the general case in which a message m i is subjected to external noise interferences and undergoes K errors before it is eventually transmitted successfully.

Intuitively, the worst case message response time (R i ) is equal to the queueing jitter (J i ) plus the transmission time of the message (C i ) plus the interference caused by three sources: lower priority messages that may be pending at the time message m i is submitted for transmission (B i ), higher priority messages that must be transmitted before m i , (I i ), and periods of network inaccessibility due to errors (E i ). Thus,

R i = J i + t i t i = B i + C i + I i (t i ) + E i (t i ) Where, b i b b C τ ) 4 1 8 34 8 44 ( ⎥ ⎦ ⎥ ⎢ ⎣ ⎢ - + + + = S C B k i lp k i + = ∈ ∀ ) ( max ) ( E i (t i ) = K M i M i = E + j j hep j C ) ( max ∈ ∀
It is difficult to know ahead of time the interference due to lower priority messages thus it is set to the worst case which is the maximum of all messages with equal or lower priority than m i plus S (the inter-frame space in a CAN frame) as given by Eq. (). In a multi-rate system, such as the one considered, we need to calculate the number of higher priority messages that must be transmitted before m i , and this includes an estimation of the multiple times the same message will be sent. This estimation is the term in square brackets in Eq. (), thus the interference due to higher priority messages is simply the sum of the estimated values a message is sent times the total time taken by each message (C i +S). Finally, the interference due to K errors is simply K times the extra time involved in sending a message in error which is estimated as the maximum length of an error frame (E) plus the maximum value of all messages of higher or equal priority than m i . For additional details the reader is referred to [].

Equation () must be solved iteratively by forming a recurrence relation with 0 i t =C i which terminates when

n i n i t t = +1 or fails when 1 + n i t > D i -J i , and D i ≤ T i .
If there is a solution R i ∀ i and R i ≤ D i then the analysis guarantees that all messages will always meet their deadlines, provided that there are no faults.

Difficulties with Tindell Equations:

1. For application engineers, the equations are not straightforward to evaluate as the formula for t i is in closed-form and must be solved iterately as outlined above. 2. The formulas assume that one has accurate estimations of the queueing jitter J i which is difficult to have in a system that is not time-triggered. 3. The numerical values for R i are not accurate because there is uncertainty in the value of B i (i.e., it is not known in advance all the messages in the set lp(i), thus the approximation of Eq. ().

) ( ) ( ) ( S C T J C t t I j i hp j j b j i i + ⎥ ⎥ ⎥ ⎤ ⎢ ⎢ ⎢ ⎡ + + - = ∑ ∈ ∀ τ
Evaluating Message Latencies in the FlexCAN Architecture As noted, there are a number of difficulties and inaccuracies in the evaluation of Tindell equations for CAN. The main sources of these difficulties and inaccuracies are: 1. Widely varying message transmission periods T i in the message set.

Lack of a time reference to base the calculations

The first difficulty stems from usage of CAN in early applications where a CAN network was the only network in the system and thus had to support all messages with widely varying message transmission periods. In fact, the application used in [Broster] has 6 messages with the smallest and largest period of 2 and 240 msec respectively, a factor of 120. With such widely differing factors in the transmission periods, the term in brackets in Eq. () is needed. The reason for having equation in closed form is also due to this difficulty. Current and future systems have several communication networks each serving an application type (e.g., entertainment, dash-board, powertrain, steer-by-wire, etc.). These networks are typically configured as a backbone network and several sub-networks each dedicated to one or few functional units of a vehicle. Messages in the sub-networks do not have widely varying message transmission periods, in fact, factors of 4 to 8 are sufficient. The second difficulty stems from the event triggered nature of the CAN protocol that does not use the notion of global time. Because of this, there is uncertainty in the assumptions of values for the queueing jitter and the interference due to lower priority messages.

FlexCAN overcomes these previous two limitations (and others which are not discussed here) by defining a sub-network for messages with closely related transmission periods and also defining a time-triggered time reference made of communication cycles divided into a number of sub-cycles. Just as is the case with other TT architectures, FlexCAN requires an off-line global message schedule to be configured. Furthermore, all messages scheduled in a sub-cycle are submitted for transmission at exactly the same time (at the beginning of the sub-cycle). For example, the following Table shows a communication cycle with four sub-cycles where the rate (reciprocal of communication period) of message m1 and m2 are 4 and 2 times respectively the rate of the remaining messages.

Table . . Example of FlexCAN message scheduling supporting multi-rates.

Sub-cycle

I II III IV messages m1, m5, m6 m1, m2, m3 m1, m4 m1, m2, m7

The above described features greatly simplifies the equations for calculating message latencies in the FlexCAN architecture. To begin with, message latency calculations are done on a sub-cycle basis for all sub-cycles. For each sub-cycle, the messages are relabeled as m1, m2, etc. with m1 the highest priority message, m2 the next highest priority message, and so on. Because all messages are submitted for transmission at exactly the same time, the jitter J i = 0 and B i = S. Thus R i = t i , resulting in the following simplifications,

R i = S + C i + I i + E i b i b b C τ ) 4 1 8 34 8 44 ( ⎥ ⎦ ⎥ ⎢ ⎣ ⎢ - + + + = E i = K ( E + j j hep j C ) ( max ∈ ∀ )
Table , shows the results of WCRT calculations for CAN and FlexCAN, using the same set of messages as that in [] for a data rate R of 250 Kbps and assuming no errors (i.e., K = 0). It can be noticed that the bound of the FlexCAN calculations are much more tighter than that for the case of CAN. 5. All messages scheduled in a sub-cycle are submitted for transmission at exactly the same time (at the beginning of the sub-cycle). 6. All messages in a sub-cycle have the same deadlines D i = T sc .

Solution of the Problem

Going back to the problem of finding the probability of failure delivery of message m i , P merr = Prob. (R i > T sc ), we start by using the notation p(R i/K ) to represent the upper bound on the probability that a frame is affected by exactly K faults and hence may arrive no later than R i/K , that is,

p(R i/K ) = Prob { R i/K < T sc } Broster et. al.
[] have provided a solution for finding P merr which is slightly modified here for the case of FlexCAN, where D i = T sc , and the periods T i are all the same in the sub-cycle in question1 . Further FlexCAN Properties 1. Successful transmission of message mi implies successful transmission of messages mk, where k < i. 2. The worst case response time R k < R i , for k < i. 3. R k/K < R i/K , for k < i, ∀ K 4. For a fixed response time R, p(R k/n ) < p(R i/j ) , for k < i, such that R k/n = R i/j = R, ∀ n,j. 5. For a fixed probability level P, R k/n > R i/j , for k < i, such that p(R k/n ) = p(R i/j ) = P, ∀ n,j.

P merr = Prob {R i/K > T sc } = 1 -∑ < ∀ i K i D R K K i R p / / / ) ( Where p(R i/n ) is computed recursively as follows, p(R i/n ) = p(n, R i/n ) -∑ - = 1 0 n j p(R i/j ) p(n-j, R i/n -R i/j )

Numerical Results

Table .

Conclusions

• Message latency calculations for FlexCAN are much simpler than those for CAN.

• Several graphs were generated that are useful for sizing, dimensioning, and configuring FlexCAN networks for specific applications. • Just like CAN, FlexCAN outperforms TT-CAN and TDMA networks for transient random failures.

• FlexCAN retains most of the error recovery properties of CAN at the expense of a smaller bandwidth utilization.

  Worst case response time (WCRT) for a FlexCAN schedule with T sc = 2.5 msec, a set of 6 messages of varying lengths, R = 1 Mbps, and K = 0the FlexCAN probabilistic analysis for R= 1 Mbps with λ = 30 faults/second appears in Figure A as a cumulative probability distribution. The graph shows the probabilities of each message exceeding a given response time. Fig. B shows the probabilities as a function of message priorities with the message response time as a parameter. Figures C and D are likewise but for a data rate R of 250 Kbps.

Figure A .

 A Figure A. Worst case probability of exceeding a response time for six messages arranged in priority order (R = 1 Mbps).

Figure B .

 B Figure B. Worst case probability of exceeding a response time as a function of message priority. Priority 1 is the highest priority (R = 1 Mbps).

Figure C .

 C Figure C. Worst case probability of exceeding a response time for six messages arranged in priority order (R = 250 Kbps).

Figure D .

 D Figure D. Worst case probability of exceeding a response time as a function of message priority. Priority 1 is the highest priority (R = 250 Kbps).

Table .

 . Comparison of worst case response time (WCRT) results for CAN and FlexCAN for a typical set of 6 messages of varying lengths, R = 250 Kbps. The calculations are easy to evaluate as the equations are not in closed-form and furthermore both I i (t i ) and E i (t i ) do not depend on t i . 3. The formulae assume that the jitter J i is zero which is enforced by the FlexCAN message schedule. 4. The formulae are accurate because the term B m = S, (i.e., there is no uncertainty regarding B m ). 5. The multi-rate case is taken into account ahead of time by the FlexCAN message global schedule thus the various periods T 1 , T 2 , … do not appear explicitly in the equations.

	I	i	( = ∑ hp j ∈ ∀	i	)	(	C	j	+	S	)
	Priority						Length	Period	Deadline	WCRT R i	WCRT R i
									C i ,(μs)	T i , (μs)	Di, (μs)	CAN	FlexCAN
	1								288	2000	2000	828	299
	2								328	4000	4000	1168	638
	3								328	4000	4000	1508	977
	4								528	8000	8000	2048	1516
	5								248	12000	12000	2608	1775
	6								528	240000	240000	2320	2314
	Calculating message response times in FlexCAN has the following advantages:
	1. The calculations are more accurate since they are done on a sub-cycle basis relative to
			a time-triggered communication cycle with each sub-cycle being independent from
			the next (i.e., all messages are sent in their respective sub-cycles)
	2. In summary, note the following properties of FlexCAN,
	1. It is a TT system with communication cycles and sub-cycles.
	2. It is a sub-network where messages do not have widely varying transmission periods
			T i .
	3. An off-line global message schedule takes care of the multi-rate case.
	4. Message latency calculations are done per sub-cycle

The deadlines and transmission periods correspond to messages already scheduled on a sub-cycle. They may be different to the deadlines and transmission periods of the entire set of messages. It should be noted that a multi-rate case has already taken into account by the scheduling process.