R Alexandre
email: radja.alexandre@univ-evry.fr

Thanh Nguyen

Long
email: longnt@hcmc.netnam.vn

A Pham

Ngoc Dinh

A MATHEMATICAL MODEL FOR THE EVAPORATION OF A LIQUID FUEL DROPLET, SUBJECT TO NONLINEAR CONSTRAINTS

Keywords: Hyperbolic system, droplet radius, nonlinear parabolic equation AMS Subject Classification: 35Q35, 76T10

We study the mathematical evolution of a liquid fuel droplet inside a vessel. In particular, we analyze the evolution of the droplet radius on a finite time interval. The model problem involves an hyperbolic system coupled with the pressure and velocity of the surrounding gas. Existence of bounded solutions for the mass fraction of the liquid, submitted to nonlinear constraints, is shown. Numerical simulations are given, in agreement with known physical experiments.

Introduction

In this paper, we are mainly concerned with the mathematical analysis of the evaporation of a single droplet in a gas, in the continuation of our previous work [START_REF] Pham | Long A mathematical model for the evaporation of a liquid fuel droplet inside an infinite vessel[END_REF] . Experimental studies around this subject are of course important for industrial purposes. Let us refer for instance to the works [START_REF] Bneson | A Modeling of droplet evaporation and coalescence for direct injection into an argon inductively coupled plasma[END_REF][START_REF] Consolini | S A molecular dynamics simulation of droplet evaporation[END_REF][START_REF] -W | A numerical study of high-pressure droplet vaporization[END_REF][START_REF] Has | B Transient heat transfer of deforming droplets at high Reynolds numbers[END_REF][START_REF] Jia | High pressure droplet vaporization. II. Laminar natural convection heat transfer over a sphere[END_REF][START_REF] Lafon | M Numerical Analysis of droplet vaporization and nurning under high-pressure conditions[END_REF][START_REF] Morin | M Vaporization and oxidation of liquid fuel droplets at high temperature and high pressure: application to n-alkanes and vegetable oil methyl esters[END_REF][START_REF] Odeide | A Simulation numérique de l'évaporation d'une goutte de combustible liquide[END_REF] . One of our main interest in this paper is to analyze the time evolution of the droplet radius, a study that we began in our previous work [START_REF] Pham | Long A mathematical model for the evaporation of a liquid fuel droplet inside an infinite vessel[END_REF] . Let us just mention here that the experimental evolution of this radius is well known as the d 2 law [START_REF] Odeide | A Simulation numérique de l'évaporation d'une goutte de combustible liquide[END_REF][START_REF] Williams | A Combustion theory[END_REF] , where d denotes the diameter of the droplet, see below for more details.

Let us recall the standard physical framework for this evolution. The evaporation of a single droplet in a gas involves simultaneous heat and mass transfer processes. In particular, heat from evaporation is transferred to the droplet boundary by conduction and convection, while vapor is carried by convection and diffusion back into the gas stream. Evaporation rate depends on the pressure, temperature and physical properties of the gas, the temperature, volatility and diameter of the drop in the spray.

To fix the ideas, in the experimental study [START_REF] Odeide | A Simulation numérique de l'évaporation d'une goutte de combustible liquide[END_REF] of a single droplet evaporation performed by the LCSR (Combustion Laboratory of the University of campus from Orleans, France), the droplet is suspended from a silicate tube. The elliptic shape of the droplet is assimilated to a sphere of equal volume. Important quantities of interest for these experiments are time evolution of the droplet radius, as well as classical quantities such as mass fractions or temperatures of the liquid and gas. In the experimental studies performed above, the so-called d 2 law is used to simplify two-phase fluid models and then propose adequate numerical schemes. This law simply states that the time evolution of the radius behaves as d 2 d 2 0 in time flow, and is purely phenomenological.

Our purpose in this paper is exactly in the opposite sense. We start from phenomenological fluid (mixtures) PDE modeling the drop evaporation process, compute the time-evolution of the drop radius, and then deduce other quantities of interest such as mass fractions of the liquid and gas. In particular, our numerical experiments are in good agreement with this phenomenological d 2 law, at least for small time evolution.

Our framework is therefore as follows: we consider a droplet initially represented as a single component mixture (liquid chemical specie 1) while the surrounding gas at time t = 0 is made of only one (gas) chemical specie, say 2.

During the evaporation process, the liquid vapor is transferred into the gas, while by condensation at the droplet surface and then by diffusion, gas chemical specie 2 appears inside the droplet. We make the important simplification that the moving interface between the droplet and the surrounding gas (i.e. between the two species) is spherical, with radius R = R(t) evolving in time.

Let ρ G (resp. ρ L ) denote gas density (resp. liquid density), and v G (resp. v L ) denote the gas velocity (resp. liquid velocity). Then, one has the classical overall continuity and momentum conservation laws

∂ t ρ k + div(ρ k v k ) = 0 (1.1) ρ k ∂ t v k + ρ k v k .∇v k = -∇p. (1.2)
Above, subscript k refers to the gas G or to the liquid L, depending on whether one considers the gas or liquid.p the state equation of the gas. Let Y L1 , Y L2 (resp. Y G1 , Y G2 ) the mass fractions of the liquid (resp. gas) obtained after diffusion of species in the surrounding gas. Therefore for two species, one has

Y G1 + Y G2 = Y L1 + Y L2 = 1.
Along with equation (1.1), we have to add the equation giving species conservation. So for the liquid, we have

ρ L ∂ t Y Lk + ρ L v L .∇Y Lk + div(ρ L Y Lk v Lk ) = -ρ L f (Y Lk ), k = 1, 2, (1.3)
Y Lk denoting mass fraction of the liquid, and f being a continuous function modeling a friction or a resistance for the drop. We assume that the liquid speed is so small that is can be settled to 0. Equations (1.1) and (1.3) can then be written under conservative form as

∂ t (ρ g) + div(ρ gv) = F ( g) (1.4) or ∂ t u + ∂ ∂x (f (u)) = F (u), (1.5)
in a system of particular coordinates.

If Γ is a curve of discontinuity of u, then one has

[f (u)] = [u] dx dt , (1.6)
where [.] denotes the jump of the inner quantity, s = dx dt is the speed of discontinuity along Γ. The jump relation (1.6) is known as Rankine-Hugoniot condition. It merely means that discontinuities cannot be completely arbitrary. The above considerations are all classical facts [START_REF] Landau | M Fluid mechanics[END_REF][START_REF] Smoller | J Shock waves end Reaction-diffusion equation[END_REF][START_REF] Ta-Tsien | L Global classical solutions for quasilinear hyperbolic ssytems[END_REF][START_REF] Williams | A Combustion theory[END_REF] . In the case of our droplet, in order to find interface condition at the droplet surface, i.e. for r = R(t), it is sufficient to use (1.5) and (1.4) in polar coordinates, getting

[ρ gv] = [ρ g] dR dt . (1.7)
Thus taking g = 1, one has

(ρ G -ρ L ) dR dt = ρ G v G -ρ L v L , that is also with v L = 0 ρ G v G - dR dt = -ρ L dR dt . (1.8)
Taking g = Y in (1.7), Y denoting the mass fraction of the liquid or the gas after diffusion, we get

(ρ G Y Gk -ρ L Y Lk ) dR dt = ρ G Y Gk (v G + v Gk ) -ρ L Y Lk v Lk
and this is equivalent to

ρ G Y Gk (v G -R ′ ) + ρ G Y Gk v Gk = -ρ L Y Lk R ′ + ρ L Y Lk v Lk . (1.9)
Above v Gk (resp.v Lk ) is the speed of specie Gk (resp.Lk), k = 1, 2.

Combining relation (1.9) with Fick's law [START_REF] Jia | High pressure droplet vaporization. II. Laminar natural convection heat transfer over a sphere[END_REF][START_REF] Landau | M Fluid mechanics[END_REF][START_REF] Williams | A Combustion theory[END_REF] , that is

Y G1 v G1 = -D 12 ∇Y G1 , Y G2 v G2 = -D 21 ∇Y G2 , D 12
and D 21 being diffusion coefficients, and with equations relating the thermodynamic state at the interface r = R(t)

Y Gk = K k Y Lk , k = 1, 2,
we obtain, for the mass fraction of the liquid Y L1 , the boundary condition Once ρ G (r, t) and v G (r, t) determined, radius R(t) of the drop suspended in the gas will be computed through the ordinary differential equation (1.8). Then we shall determine the mass fraction Y L1 of liquid after evaporation process, through the PDE (1.3) along with boundary condition (1.10), for a given suitable function f . For this last purpose, within the framework of weighted Sobolev spaces on initial data and for some continuous function f subject to increasing condition, we shall provide an unique local solution for the mass fraction Y L1 of the liquid. In addition, we shall show that if the initial condition is bounded, then so is our solution.

∂ r Y L1 + R ′ (t)(K 1 -1) K 2 ρ G (R(t), t) -K 3 Y L1 = 0 at r = R(t), (1.10
In the numerical applications (2nd example) we have chosen the experimental conditions made by the LCSR in the study of single drop evaporation and in this case the study of the radius of the drop shows us that the graphic associated to our mathematical model presents the same features as in the experimental curves.

Plan of the paper: In Section 2, by using Riemann invariants, we determine the droplet radius. This enables to get, in Section 3, the liquid mass fraction, using a variational method. Finally, we have presented some numerical simulations in the last Section, which shows that our model is is good agreement with experimental simulations, at least for short time.

Hyperbolic system and droplet radius

The gas velocity v G (r, t) and its density ρ G (r, t) satisfy the following system, using polar coordinates

∂ρ G ∂t + 1 r 2 ∂ ∂r r 2 ρ G v G = 0 ∂ ∂t (ρ G v G ) + 1 r 2 ∂ ∂r r 2 ρ G v 2 G = - ∂p ∂r .
(2.1)

Setting ρ(r, t) = r 2 ρ G (r, t) and v(r, t) = v G (r, t), we have ∂ρ ∂t + v ∂ρ ∂r + ρ ∂v ∂r = 0 ∂v ∂t + v ∂v ∂r = - 1 ρ ∂p 1 ∂r , (2.2)
where p 1 (r, t) is an auxiliary function connected to the state equation of the gas p(r, t), by

∂p 1 ∂r = r 2 ∂p ∂r .
In (2.2), according to the discussion in the Introduction, we choose the auxiliary function p 1 (r, t) as

p 1 =ρ γ , γ > 1.
With this choice, we get the following system ∂ρ ∂t

+ v ∂ρ ∂r + ρ ∂v ∂r = 0 ∂v ∂t + v ∂v ∂r + γρ γ-2 ∂p ∂r = 0. (2.3)
We note that (2.3) is equivalent to matrix form

∂ ∂t ρ v + A. ∂ ∂r ρ v = 0, (2.4)
A being the (2, 2) matrix

A = v ρ γρ γ-2 v . Eigenvalues of A (characteristics speeds) are given λ = v -c and µ = v + c, c = γρ γ-1 . Since λ < µ, system (2.4
) is therefore hyperbolic. Thus there exists two functions W (ρ, v) and Z(ρ, v) (Riemann invariants) such that

W (ρ, v) = constant on dX 1 dt = λ, (2.5)
Z(ρ, v) = constant on dX 2 dt = µ. (2.6)
W (ρ, v) is determined by the system dv c/ρ = dρ 1 , vector R 1 = (1, c/ρ) being and eigenvector associated to the eigenvalue µ. Thus

W (ρ, v) = v - 2c γ -1 . (2.7)
Similarly, the Riemann invariant Z(ρ, v) corresponding to λ is given by

Z(ρ, v) = v + 2c γ -1 . (2.8)
Functions W (ρ, v) = W (t, r) and Z(ρ, v) = Z(t, r) satisfy the following system, equivalent to system (2.4)

∂W ∂t + λ(W, Z) ∂W ∂r = 0 ∂Z ∂t + µ(W, Z) ∂Z ∂r = 0, (2.9)
where λ(W, Z) and µ(W, Z) are given by

λ = -γ-3 4 Z + γ+1 4 W µ = γ+1 4 Z -γ-3 4 W, (2.10)
as follows from (2.7) and (2.8).

It is well known that a sufficient condition in order that (2.9) is authentically nonlinear is that ∂λ ∂W > 0 and ∂µ ∂Z > 0, which is the case here according to (2.10). Integration along the characteristics defined by

dX 1 dt = λ(W, Z), X 1 (0) = β gives X 1 (0,β) (t) = β + t 0 λ(W (s, X 1 (s)), Z(s, X 1 (s)))ds. (2.11)
Therefore, the solution of the initial value problem

W t + λ(W, Z)W r = 0, W (0, r) = W 0 (r) (2.12)
can be written as

W (t, r) = W 0 X 1 (0,β) (0) = W 0 (β), (2.13)
where

β = r - t 0 λ(W (s, X 1 (s)), Z(s, X 1 (s)))ds. Similarly we have Z(t, r) = Z 0 X 2 (0,α) (0) = Z 0 (α), (2.14)
where

X 2 (0,α) (t) = α + t 0 µ(W (s, X 2 (s)), Z(s, X 2 (s)))ds.
The above considerations lead to the following Proposition 2.1. Assume that W ′ 0 (β) < 0 or Z ′ 0 (α) < 0. Then solution of system (2.3) is defined on a finite interval [0, T [. Proof: Differentiation of (2.5) and (2.11) with respect to β gives

dX 1 β dt = λ β (W, Z) with X 1 β (t = 0) = 1, (2.15)
and in the same way

dX 2 α dt = µ α (W, Z) with X 2 α (t = 0) = 1. (2.16) Since λ β = λ W W β + λ Z Z β = λ W W ′ 0 (β) and µ α = µ Z Z ′ 0 (α), from (2.15
) and (2.16), integrating w.r.t. t along the characteristics yields

X 1 β (t) = 1 + t 0 λ W W ′ 0 (β)dt = 1 + γ + 1 4 W ′ 0 (β)t, (2.17) X 2 α (t) = 1 + t 0 µ Z Z ′ 0 (α)dt = 1 + γ + 1 4 Z ′ 0 (α)t. (2.18) From (2.17), it follows that X 1 β (t 1 ) = 0 for t 1 = -4 (γ + 1)W ′ 0 (β) > 0. Similarly X 2 α (t 2 ) = 0 for t 2 = -4 (γ + 1)Z ′ 0 (α) > 0.
Hence ∂W ∂r (t, r) becomes infinite for

T = inf {t 1 , t 2 }, since ∂W ∂r = W β . dβ dX = W ′ 0 (β) X β .
From Proposition 2.1, it follows that Proposition 2.2. System (2.9) admits an unique C 1 solution on [0, T [, for all r ∈ R + and for initial data

ρ G (0, r) = ρ 0 (r) and v G (0, r) = v 0 (r) belonging to C 1 (R + ).
Concerning the droplet radius, it follows from (1.8), that we have the following ode for this radius

dR(t) dt = v G (t, R(t))ρ G (t, R(t)) ρ G (t, R(t)) -ρ L , R(0) = R 0 . ( 2
(r) such that W ′ 0 (r) < 0 or Z ′ 0 (r) < 0.

Liquid Mass Fraction

The liquid mass fraction Y L1 satisfies the conservation equation of specie (1.3), which can be rewritten as

∂ t Y L1 - 1 r 2 ∂ ∂r r 2 ∂ ∂r Y L1 + f (Y L1 ) = 0. (3.1)
We have used polar coordinates, and taken the diffusion constant D 12 as being equal to 1. Of course, (3.1) is equivalent to

∂ t Y L1 -∆Y L1 - 2 r ∂ ∂r Y L1 + f (Y L1 ) = 0 , for 0 < r < s(t), (3.2)
where s(t) = R(t) denotes the droplet radius determined in section 2.

The boundary condition at the surface s(t) is given by the Rankine-Hugoniot condition connected to the thermodynamic equilibrium, i.e. formula (1.10).

Performing the change of variable r = R(t)x, function Y L1 (t, r) turns to function Y L1 (t, R(t)x) = u(t, x), which satisfies the following initial boundary value (i.b.v.) problem

∂ t u -a(t) ∆u + 2 x ∂ x u -x R ′ (t) R(t) ∂ x u + f (u) = 0, 0 < x < 1, t > 0, (3.3
)

|lim x→0+ xu x (t, x)| < ∞, u x (t, 1) + k(t)u(t, 1) = 0, (3.4) u(0, x) = u 0 (x), (3.5)
where we used the following notations We shall do so by setting this problem in a variational framework, using weighted Sobolev spaces.

a(t) = 1 R 2 (t) , k(t) = R(t)R ′ (t)(K 1 -1) K 2 ρ G (t, R(t)) -K 3 . ( 3
Let Ω =]0, 1[ and define H as the Hilbert space given by H = {v : Ω → R, measurable and such that

1 0 x 2 v 2 (x)dx < +∞}. Note that H is the closure of C 0 ( Ω) w.r.t. the norm v H = 1 0 x 2 v 2 (x)dx 1/2
. We also introduce the real Hilbert space V = v ∈ H |v ′ ∈ H . In the following, we shall often use the fact that V

is the closure of C 1 ( Ω) w.r.t. the norm v V = v 2 H + v ′ 2 H 1/2
. V is continuously embedded in H. Identifying H with his dual H ′ , one has V ⊂ H ⊂ V ′ with continuous injections. Note also that the norms . H and . V can be defined, respectively, from the inner products

< u, v >= 1 0 x 2 u(x)v(x)dx and < u, v > + < u ′ , v ′ >.
We then have the following results, the proofs of which can be found in the paper 12 , Lemma 3.1.

For every v ∈ C 1 ([0, 1]), ǫ > 0 and x ∈ [0, 1] we have v 2 0 ≤ 1 2 v ′ 2 0 + v 2 (1), v 2 (1) ≤ ǫ v ′ 2 0 + C ǫ v 2 0 , v(1) ≤ 2 v 1 , xv(x) ≤ √ 5 v 1
where

C ǫ = 3 + 1 ǫ and . 0 = . H , . 1 = . V . Lemma 3.2. The embedding V ⊂ H is compact. Remark 3.1. From Lemma 3.1, it follows that v ′ 2 0 + v 2 (1) 1/2
and v 1 are two equivalent

norms on V since 2 3 v 2 1 ≤ v 2 (1) + v ′ 2 0 ≤ 5 v 2 1 , for all v ∈ V . Remark 3.2. We have xv(x) ∈ C 0 ([0, 1]
), for all v ∈ V . Indeed, on one hand, lim x→0+ xv(x) = 0, ∀v ∈ V (see the book 1 , p.128), and on the other hand

v |[ǫ,1] ∈ C 0 ([ǫ, 1]), ∀ǫ, 0 < ǫ < 1, since we have H 1 (ǫ, 1) ⊂ C 0 ([ǫ, 1]) and ǫ v H1(ǫ,1) ≤ v 1 ∀v ∈ V ∀ǫ, 0 < ǫ < 1.
If X is any Banach space, we denote by . X its norm, and by X ′ the dual space of X. We denote by L p (0, T ; X), 1 ≤ p ≤ ∞, the standard Banach space of real functions u : (0, T ) → X, measurable, such that

u L p (0,T ;X) = T 0 u(t) p X dt 1/p < +∞, for 1 ≤ p < ∞ and u L ∞ (0,T ;X) = ess sup 0<t<T u(t) X , for p = ∞. Let u(t), u ′ (t), u x (t), u xx (t) denote u(t, x), ∂u ∂t (t, x), ∂u ∂x (t, x), ∂ 2 u ∂x 2 (t, x) respectively.
We shall make the following set of assumptions:

• (H 1 ) u 0 ∈ H; • (H 2 ) a, k ∈ W 1,∞ (0, T ), a(t) ≥ a 0 > 0; • (F 1 ) f ∈ C(R, R); • (F 2 ) There exists positive constants C 1 , C ′ 1 , C 2 and p, 1 < p < 3, such that (i) uf (u) ≥ C 1 |u| p -C ′ 1 , (ii) |f (u)| ≤ C 2 (1 + |u| p-1 ). Let u ∈ C 2 ([0, T ] × [0, 1]
) be a solution of problem (3.3)-(3.5). Then, after multiplying equation (3.3) by x 2 v, v ∈ V w.r.t. the scalar product of H, integrating by parts and taking into account boundary condition given by (3.4), we get

d dt < u(t), v > + a(t) 1 0 x 2 u x v x dx + a(t)k(t)u(1)v(1) - R ′ (t) R(t) 1 0 x 3 u x vdx+ < f (u), v >= 0
The weak formulation of the ibv problem (3.3)-(3.5) can then be given in the following way: Find u(t), defined on the open set (0, T ), such that u(t) satisfies the following variational problem

d dt < u(t), v > + ã(t; u(t), v)+ < f (u(t)), v >= 0, ∀v ∈ V, (3.7)
together with the initial condition

u(0) = u 0 . (3.8)
Above, we have used the following bilinear form

ã(t; u, v) = a(t) 1 0 x 2 u x v x dx + a(t)k(t)u(1)v(1) - R ′ (t) R(t) 1 0 x 3 u x vdx, u, v ∈ V . (3.9)
We first note the following lemma, the proof of which can be found in our previous paper 2 Lemma 3.3. There exists constants K T , α T and β T depending on T , such that

| ã(t; u, v)| ≤ K T u 1 v 1 , for all u, v ∈ V, (3.10) ã(t; u, u) ≥ α T u 2 1 -β T u 2 0 , u, v ∈ V. (3.11)
Integrating (3.17), using (3.15), it follows that

S m (t) ≤ C 0 + 2 3 T C ′ 1 + 2β T t 0 S m (s)ds, (3.18)
where

S m (t) = u m (t) 2 0 + 2α T t 0 u m (s) 2 1 ds + 2C 1 t 0 ds 1 0 x 2 |u m (s, x)| p dx, (3.19)
and C 0 is a constant depending only on u 0 with u 0m 2 0 ≤ C 0 ∀m. Applying Gronwall's lemma, we obtain from (3.18)

S m (t) ≤ C 0 + 2 3 T C ′ 1 exp(2β T t) ≤M T , ∀m, ∀t, 0 ≤ t ≤ T m ≤ T, (3.20) that is T m = T .
In the following, we denote by M T any generic constant depending only on T .

(b) Second estimate. Replacing w j by t 2 u m in (3.8) gives

tu ′ m 2 0 + 1 2 d dt a(t) tu m 2 0 + a(t)k(t)t 2 u 2 m (1) + 1 2 d dt t 2 1 0 x 2 f (u m )dx = u mx 2 0 d dt t 2 a(t) + 1 2 u 2 m (1) d dt (t 2 a(t)k(t)) + 2t 1 0 x 2 f (u m )dx + R ′ (t)t 2 R(t)
2 0 + a(t)t 2 u 2 m (t, 1) = a(t)(1 -k(t))t 2 u 2 m (t, 1) + t 0 s 2 a(s) ′ u mx (s) 2 0 ds + t 0 s 2 a(s)k(s) ′ u 2 m (s, 1)ds + 2 t 0 R ′ (s) R(s) s 2 < xu mx (s), u ′ m (s) > ds +4 t 0 sds 1 0 x 2 f (u m (s, x))dx -2t 2 1 0 x 2 f (u m )dx. (3.23)
By means of assumption (H 2 ) and Remark 3.1, we get

a(t) tu m (t) 2 0 + a(t)t 2 u 2 m (t, 1) ≥ 2 3 a 0 tu m (t) 2 1 ∀t ∈ [0, T ], ∀m. (3.24)
We fix ǫ > 0 such that

a ∞ 1 -k ∞ ǫ < a 0 3 , (3.25)
where . ∞ = . L∞(0,T ) .

Using again Lemma 3.1, Remark 3.1 with ǫ > 0 as in (3.25) and first estimate (3.20), the terms on the r.h.s. of (3.23) can be estimated as follows From assumptions (F 1 ) and (F 2 ), we note also that

a(t)(1 -k(t))t 2 u 2 m (t, 1) ≤ a ∞ 1 -k ∞ ǫ tu m (t) 2 1 + C ǫ tu m (t) 2 0 ≤ a 0 3 tu m (t) 2 1 + M T , (3.26) t 0 s 2 a(s) ′ u mx (s) 2 0 ds + t 0 s 2 a(s)k(s) ′ u 2 m (s, 1)ds ≤ (t 2 a) ′ ∞ + (t 2 ak) ′ ∞ t 0 u mx (s) 2 0 + u 2 m (s, 1) ds ≤ 5 (t 2 a) ′ ∞ + (t 2 ak) ′ ∞ t 0 u m (s) 2 1 ds ≤ M T , (3.27) 2 t 0 s 2 R ′ (t) R(t) < xu mx (s), u ′ m (s) > ds ≤ t 0 su ′ m (s) 2 0 ds + R ′ R 2 ∞ t 0 su m (s)
-m0 = - z0 -z0 |f (y)|dy ≤ f (z) = z 0 f (y)dy ≤ C 2 |z| + |z| p p , ∀z ∈ R, (3.29)
where

z 0 = (C ′ 1 /C 1 )
x 2 f (u m (s, x))dx -2t 2 1 0 x 2 f (u m (t, x)dx ≤ 4C 2 t 0 sds 1 0 x 2 |u m (s, x)| + 1 p |u m (s, x)| p dx + 2t 2 1 0 x 2 m0 dx ≤ 4C 2 t 0 s u m (s) 0 ds + 4 p C 2 t t 0 ds 1 0 x 2 |u m (s, x)| p dx + 2 3 T 2 m0 ≤ 2C 2 T M T + 2C 2 pC 1 T M T + 2 3 T 2 m0 ≤ M T . (3.30)
Hence, we deduce from (3.23), (3.24), (3.26)-(3.28) and (3.30) that

t 0 su ′ m (s) 2 0 ds + a 0 3 tu m (t) 2 1 ≤ M T + R ′ R 2 ∞ t 0 su m (s) 2 1 ds. (3.31)
By Gronwall's lemma, we get

t 0 su ′ m (s) 2 0 ds + a 0 3 tu m (t) 2 1 ≤ M T exp R ′ R 2 ∞ . 3T a 0 ≤ M T , ∀t ∈ [0, T ]. (3.32)
Finally, using (3.20) and assumption (F 2 ,ii) we have also

t 0 ds 1 0 x 2/p ′ f (u m (s, x)) p ′ dx ≤ (2C 2 ) p ′ t 0 ds 1 0 x 2 |u m (s, x)| p dx ≤ M T , (3.33)
with p ′ = p p -1 .

• Step 3, the limiting process.

From (3.20), (3.32) and (3.33), we deduce that there exists a subsequence of {u m }, still denoted {u m } such that

u m → u weakly * in L ∞ (0, T ; H) , u m → u weakly in L 2 (0, T ; V ) , x 2/p u m → x 2/p u weakly in L p (Q T ) , tu m → tu weakly * in L ∞ (0, T ; V ) , (tu m ) ′ → (tu) ′ weakly in L 2 (0, T ; H). (3.34)
Using a standard compactness lemma 10 (p.57) together with (3.34), we can extract from the sequence {u m }, a subsequence still denoted by {u m } such that

tu m → tu strongly in L 2 (0, T ; H). (3.35)
Continuity of f also implies (up to a sub-sequence)

f (u m (t, x)) → f (u(t, x)) a.e. (t, x) ∈ Q T = (0, T ) × (0, 1). (3.36)
Applying a standard weak convergence lemma 10 , we have also

x 2/p ′ f (u m ) → x 2/p ′ f (u) weakly in L p ′ (Q T ).
Passing to the limit in (3.13) and (3. First of all, we note the following slight extension of a lemma used in our previous paper 2 (see also the book 10 ) Lemma 3.4. Let w be the weak solution of the following i.b.v. problem We now turn to the boundness of the above solutions.

w t -a(t)(w xx + 2 x w x ) = f (t, x), 0 < t < T, 0 < x < 1, lim x→0+ xw x (t, x) < +∞, w x (t, 1) + k(t)w(t, 1) = 0, w(0, x) = 0, w ∈ L 2 (0, T ; V ) ∩ L ∞ (0, T ; H), x 2/p w ∈ L p (Q T ), tw ∈ L ∞ (0, T ; V ), tw t ∈ L 2 (0, T ; H). Then 1 2 w(t) 2 0 + t 0 a(s) w x (s) 2 0 + k(s)w 2 (s, 1) ds - t 0 < f (s), w (
For this purpose, we shall make use of the following assumptions

• (H' 1 ) u 0 ∈ L ∞ (0, 1), |u 0 (x)| ≤ M, a.e.
x ∈ (0, 1)

• (H' 2 ) a, k ∈ W 1,∞ (0, T ), a(t) ≥ a 0 > 0, k(t) ≥ k 0 > 0 • (F' 1 ) uf (u) ≥ 0 ∀u ∈ R such that |u| ≥ u 0 ∞ , for a.e.
x ∈ (0, 1).

We then have the following result Theorem 3.2. Let (H' 1 ), (H' 2 ), (F 1 )-(F 3 ) and (F' 1 ) hold. Then the unique weak solution of the ibv problem (3.7)-(3.9), as given by theorem 1, belongs to L ∞ (Q T ).

Proof of Theorem 3.2. Firstly, we note that Z = u -M satisfies the i.b.v. problem x 2 Z t vdx + a(t)

∂ t Z -a(t) ∆Z + 2 x ∂ x Z -x R ′ (t) R(t) ∂ x Z + f (Z + M ) = 0, 0 < x < 1, t ∈ (0, T ), ( 3
1 0 x 2 Z x v x dx + a(t)k(t)Z(t, 1)v(1) - R ′ (t) R(t) 1 0 x 3 Z x vdx + 1 0 x 2 f (Z + M )vdx = -M a(t)k(t)v(1), ∀v ∈ V , (3.41)
hence for v = Z + = 1 2 Z + |Z| , since u 0 ∈ L ∞ (0, 1). It follows that 1 2 d dt 1 0 x 2 |Z + | 2 dx + a(t) 1 0 x 2 |(Z + ) x | 2 dx + a(t)k(t)|Z + (t, 1)| 2 - R ′ (t) R(t) 1 0 x 3 Z + x Z + dx + 1 0 x 2 f (Z + + M )Z + dx = -M a(t)k(t)Z + (t, 1) ≤ 0, since 1 0 x 2 Z t Z + dx = 1 0,Z>0 x 2 (Z + ) t Z + dx = 1 2 d dt 1 0 x 2 |Z + | 2 dx.
On the other hand, by assumption (H' 2 ) and Remark 3.1, one has a(t)

1 0 x 2 Z + x 2 dx + a(t)k(t) Z + (t, 1) 2 ≥ C0 Z + (t) 2 1 , (3.42)
where C0 = 2 3 a 0 min {1, k 0 }. Using the monotonicity of f (u) + δu and (F' 1 ), we have

1 0 x 2 f (Z + + M )Z + dx = 1 0 x 2 f (Z + + M ) -f (M ) Z + dx + 1 0 f (M )x 2 Z + dx ≥ -δ 1 0 x 2 Z + 2 dx + 1 0 f (M )x 2 Z + dx ≥ -δ Z + 2 0 . ( 3
x 3 Z + x Z + dx yields d dt Z + (t) 2 0 + C0 Z + (t) 2 1 ≤ 1 C0 R ′ R 2 ∞ + 2|δ| Z + (t) 2 0 . (3.44)
Integrating (3.44), we get

Z + (t) 2 0 ≤ Z + (0) 2 0 + 1 C0 R ′ R 2 ∞ + 2|δ| t 0 Z + (s) 2 0 ds. (3.45) Since Z + (0) = u(0, x) -M + = u 0 (x) -M + = 0, Gronwall's lemma yields Z + (t) 0 = 0. Thus u(t, x) ≤ M a.e. (t, x) ∈ Q T .
The case u 0 (x) ≥ -M is similar, by considering Z = u + M and Z -= 1 2 |Z| -Z . Thus we get Z -= 0 and hence u(t, x) ≥ -M a.e. (t, x) ∈ Q T . All in all, one obtains |u(t, x)| ≤ M a.e. (t, x) ∈ Q T and this ends the proof of Theorem 3.2.

Numerical applications

For the numerical applications, we have taken in (2.3) γ = 3, so that equation (2.9) reduces to Burger's equation

W t + W W r = 0, W (0, r) = W 0 (r) Z t + ZZ r = 0, Z(0, r) = Z 0 (r). (4.1)
It is well known that classical Burger's equation

u t + uu r = 0, u(0, r) = u 0 (r)
admits the solution u(t, r) = u 0 (ξ(t, r)), ξ(t, r) being defined by the parametrization r = u 0 (ξ)t + ξ.

Having in mind (4.1), we have considered two examples.

• First example.

For the first example, we have chosen the initial conditions W 0 (r) = 1, r > 0; W 0 (0) = 0 and Z 0 (r) = 2, r > 0; Z 0 (0) = 0. The continuous solutions of (4.1) are then given by

W (t, r) = 1 if 0 ≤ t ≤ r, r t if 0 ≤ r ≤ t, Z(t, r) = 2 if 0 ≤ 2t ≤ r, r t if 0 ≤ r ≤ 2t.
According to Section 2, the droplet radius is given by formula (2.19) which we consider here with an initial condition taken equal to be 1

dR(t) dt = v G (t, R(t))ρ G (t, R(t)) ρ G (t, R(t)) -ρ L , R(0) = 1, (4.2)
and where

   v G (t, r) = 1 2 (W (t, r) + Z(t, r)) ρ G (t, r) = 1 2 √ 3r 2 (Z(t, r) -W (t, r)). ( 4

.3)

In figure 1 below, we have drawn the curve t -→ R(t) on the time interval [0, 1] with a step h = 0.05 and ρ L = 0.9.

Fig. 1 • Second example.

For the second example, we have chosen the truly experimental conditions made by the LCSR in the study of single drop evaporation, the drop being suspended from a silicate tube. Drops are made up of n-heptane fuel (ρ L = 683 kg/mm 3 )) in air at normalized atmospheric pressure and with an initial speed v G (0, r) = C 1 = 35mm/s. The initial density ρ G (0, r) of the gas is taken as C 2 = 348 T 0 , T 0 = 373K. In this case, the solution of (4.1) are given by

W (t, r) = C 1 - √ 3C 2 ξ 2 , ξ = 1 + 1 -4(r -C 1 t) √ 3C 2 t 2 √ 3C 2 t , Z(t, r) = C 1 + √ 3C 2 η 2 , η = -1 + 1 + 4(r -C 1 t) √ 3C 2 t 2 √ 3C 2 t .
We then compute v G and ρ G by formula (4.3), and then solve the ode for R given by (4.2). We note that R ′ (0) = C1C2 C2-ρL < 0 where C 2 = 348 372 < 1 < ρ L = 683. On the other hand, one has

v G (t, r) = 1 4 √ 3C 2 t [C 1 + √ 3 12t 2 C 2 (α + √ 1 -α - √ 1 + α)],
where 2 and3 represent resp. the velocity v G (t, r) and the pressure ρ G (t, r) given by (4.3) for (t, r) ∈ (0, 1) × (0, 1). The curve of the radius t -→ R(t) for this case is drawn in figure 4. Since W ′ 0 (ξ) < 0 the maximal existence interval is finite (Proposition 2.1) as can be seen in our graphic. Let us remark that looking on the experimental curves 11 made by the LCSR (figure 5) at the beginning, the function t -→ R(t) is increasing around the vicinity of t = 0. This fact is confirmed by our model which represents a good improvement of our previous model 2 in which the velocity v G (t) was a given function of t.

α = α(t, r) = 4(r -C 1 t) √ 3C 2 t,

. 6 )

 6 Our purpose in this Section is to analyze the boundary value problem (3.3)-(3.5).

  14), it follows from (3.15), (3.34) and (3.36), that function u(t) satisfies the i.b.v. problem (3.7), (3.8). • Step 4. Uniqueness of the solutions.

  s) > ds = 0, a.e. t ∈ (0, T ) Uniqueness of solutions for our initial i.b.v problem will then be deduced as follows. Let u and v be two weak solutions of (3.3)-(3.5). Then w = u -v is a weak solution of problem mentioned in Lemma 3.4, with r.h.s. given by f(t, x) = xR ′ (t) R(t) w x -f (u) + f (v). Therefore, Lemma 3; w(s), w(s))ds + 2 t 0 < f (u(s)) -f (v(s)), w(s) > ds = 0.Using Lemma 3.3 and assumption (F 3 ) we obtain δ + β T ≥ 0 we have w(t) 0 = 0 by applying Gronwall's lemma. In the case where δ + β T < 0 the result is clearly still true. This ends the proof of Theorem 3.1.

  .38)|lim x→0+ xZ x (t, x)| < ∞, Z x (t, 1) + k(t)[Z(t, 1) + M ] = 0, (3.39) Z(0, x) = u 0 (x) -M(3.40) Multiplying equation (3.38) by x 2 v, for v ∈ V , integrating by parts w.r.t. variable x and taking into account boundary condition (3.39), one has 1 0

  Fig.2

  Fig. 3

Fig. 4

 4 Fig. 4

Fig. 5

 5 Fig. 5

  1/p .

	Using first estimate (3.20), (3.29) and Lemma 3.1, we obtain
	t	1
	4	sds
	0	0

  |α| ≤ 1. One has sgnv G (t, r) = sgn[C 1 12t 2 C 2 + G (t, r) > 0 if C 1 12t 2 C 2 ≥ √ 6 i.e. t ≥ t 2 = 1 G (t, r) < 0 on the interval [t 1 , t 2 [ . Similarly 0 ≤ ρ G (t, r) ≤ ≤ 1 3t 2 r 2 . Thus, if t ≥ t m , r ≥ r m , ρ G (t, r) ≤ 1 L ≃ 683, that is for t m r m ≥ 1 3.683 ≃ 1 45 . Since R ′ (0) < 0 it follows that R(t)is decreasing on [0, t 1 [, increasing on [t 1 , t 2 [ and then from the starting point t 2 always non increasing. Figures

	√ Thus v C12 3(α + √ 1 -α -√ 3(α + √ 1 -α -√ √ 6. 1 + α) ≥ -1 4 ≃ 1 35 and v 1 √ 1 + α)]. Since |α| ≤ 1, it follows that C 2 2 t 2 6.2.r 2 [2 + √ 1 -α -√ 1 + α] 3t 2 m r 2

m ≤ ρ

We then have the following existence theorem Theorem 3.1. Let T > 0 and assumptions (H 1 ),(H 2 ),(F 1 ),(F 2 ) hold true. Then, there exists a solution u of the variational problem (3.7), (3.8)

Furthermore, if f satisfies the additional condition

for all u, v ∈ R, for some δ ∈ R, then the above solution u is unique.

Proof of Theorem 3.1. We divide it in several steps.

• Step 1, Galerkin method. Denote by {w j }, j = 1, 2, ..., an orthonormal basis of the separable Hilbert space V . We wish to find u m (t) of the form

where c mj (t) satisfy the following system of nonlinear differential equations

together with the initial condition

and

Clearly, for each m, there exists an unique local solution u m (t) of the form (3.12), which satisfies (3.13) and (3.14) almost everywhere on 0 ≤ t ≤ T m , for some T m , 0 < T m ≤ T. The following estimates allow us to take T m = T for all m.


Dorota Duda
email: d.duda@pb.edu.pl

Marek Kretowski
email: m.kretowski@pb.edu.pl

Johanne Bezy-Wendling

Effect of Slice Thickness on Texture-Based Classification of Liver Dynamic CT Scans

Keywords: texture analysis, classification, dynamic CT, liver, slice thickness, stability of the parameters, tissue characterization, CAD

This paper assesses the impact of slice thickness on texture parameters. Experiments are performed on liver dynamic CT scans, with two slice thicknesses. Three acquisition moments are considered: without contrast, in arterial and in portal phase. In total, 155 texture parameters, extracted with 9 methods, are tested. Classification of normal and cirrhotic liver is performed using a boosting algorithm. Experiments reveal that slice thickness does not considerably influence the stability of the parameters. They also enable to assess the rate of parameter dependency on slice thickness. Finally, they show that applying different slice thicknesses for training and testing the CAD system requires slice thickness-independent parameters.

Introduction

Various medical imaging techniques are available for acquiring a diagnostic information, e.g., Computed Tomography (CT), Positron Emission Tomography (PET), ultrasound, Medical Resonance Imaging (MRI), Single Photon Emission Computed Tomography (SPECT). With their rapid development, the image quality has been significantly improved. The images acquired over the years are characterized by increasing spatial and grey level resolution, and decreasing slice thickness. The number of images obtained within a single examination increases. The appropriate interpretation of an image information is thus a very complex task, but decisive for a correct diagnosis and treatment proposal.

Due to the fact, that the men are not able to identify such huge amount of information stored in the images, the (semi)automatic Computer Aided Diagnosis (CAD) systems become of a rapidly growing interest. Combining different image analysis techniques (like texture analysis [START_REF] Haralick | Statistical and structural approaches to texture[END_REF] or segmentation) and classification algorithms, they appear to be a very promising diagnostic tool.

In most of the CAD systems based on imaging data two main stages of work can be distinguished. The first one, called training, is a preparation of the system for the recognition of several classes of tissue, representing different pathologies. The second stage is the applicationof the system for the (semi)automatic recognition of new cases, not yet diagnosed. Due to the constant changes of image acquisition protocols, the question arises whether the protocols for new recognized images can be different from those that were used during the system training. The question seems all the more important given the fact that there is still no universal consensus on image acquisition protocols, so the images obtained from different machines could be of different properties and qualities.

The aim of our study is to examine the effect of slice thickness on texture parameters characterizing hepatic tissue on dynamic (contrast-enhanced) CT images. Our choice was dictated by the fact that the image database that we have been creating for about 10 years includes images of several slice thicknesses: from the oldest ones, of 10 mm, to the most recent ones, of 1.3 mm. So far, we have not found any research concerning the impact of slice thickness on texture-based classification of liver CT images with different contrast product concentrations.

In this study, we first assess the influence of slice thickness on parameter stability. Secondly, we study the possibility of tissue differentiation, with the most known parameters and different combinations of slice thicknesses used for system training and testing. Then, the parameter dependency on slice thickness is evaluated. Finally, the classification of two types of liver tissue, characterized by parameters which are least dependent on slice thickness is performed.

The next section includes a short description of related works. In Sect. 3 the system for classification of multiphasic textures is presented. Then, the methods for assessing the effect of slice thickness on the system performance are proposed. An experimental validation is described in Sect. 4. Conclusions and future work are presented in the last section.

Related Work

The effect of image acquisition protocols on CAD system performances has already been investigated in several studies.

In [START_REF] Guggenbuhl | Reproducibility of CT-based bone texture parameters of cancellous calf bone samples: Influence of slice thickness[END_REF] the influence of slice thickness on CT-based texture parameters for cancellous calf bone was studied. Four different slice thicknesses were considered.

The work [START_REF] Mayerhoefer | Effects of MRI acquisition parameter variations and protocol heterogeneity on the results of texture analysis and pattern discrimination: an application-oriented study[END_REF] analyzed the sensitivity of texture features of five different categories to variations in the number of acquisitions, repetition time, echo time, and sampling bandwidth at different spatial resolutions of T2-weighted MR images.

In [START_REF] Miles | Colorectal cancer: texture analysis of portal phase hepatic CT images as a potential marker of survival[END_REF] the effect of acquisition parameters, such as tube currents and voltages (10 different combinations), on texture was assessed on CT images of a cylindric phantom filled with water. Eight different texture parameters were analyzed.

The effect of slice thickness on brain MRI texture classification was studied in [START_REF] Savio | Effect of slice thickness on brain magnetic resonance image texture analysis[END_REF]. In the work, thick slices were simulated on the basis of thin ones.

The first three of the cited studies showed that changes in all analyzed acquisition parameters could influence the texture parameter values, and thus -the texture discrimination. Only in the last work, the two considered tissue classes were separable even if slice thicknesses differed between training and testing sets.

Texture-Based Classification of Hepatic Tissues

Two Stages of Work: Training and Aiding a Diagnosis

The system we develop for texture-based classification of multiphasic liver CT scans also works according to previously described two-stage model. Since typical CT exams of abdominal organs often consider three acquisition moments, which are related to the contrast product propagation in hepatic vessels, we proposed to analyze triplets of images [START_REF] Duda | Texture-based classification of hepatic primary tumors in multiphase CT[END_REF], [START_REF] Duda | Texture characterization for Hepatic Tumor Recognition in Multiphase CT[END_REF]. The first of the three simultaneously analyzed images corresponds to acquisition without a contrast. The second and the third are taken after its injection, at arterial and portal phases of contrast propagation.

The main steps of system training are shown in Fig. 1. First, a database of image triplets (preprocessed, if needed) is created. Then a Region of Interest (ROI) is drawn at the same position on each of the three corresponding images. Each ROI is characterized by the same vector of texture parameters. Three parameter vectors are thus created, each of them characterizes a texture at different acquisition moment. In the next step, parameters from those three vectors are concatenated, in order to make one "complex" parameter vector, describing a triplet of textures. The label representing a pathology is associated with each complex vector. At this moment, a set of labeled complex vectors could be subjected to a feature selection. Finally, it is used for the construction of the classifiers. Then the second stage of the system work can take place.

In order to recognize a new observation (see Fig. 2), a triplet of images is necessary. The three simultaneously analyzed images (in no contrast, arterial and portal phase) are subjected to the same preprocessing as it was in the training stage. Then three ROIs are drawn -one ROI on each of the three images. The triplet of ROIs is characterized in the same way as in the training stage. Finally, a complex vector of concatenated parameters (corresponding simultaneously to three acquisition moments) is classified. The tissue class that is attributed to this vector, is one of the classes considered in the training stage.

During the system design, many aspects must be investigated to ensure the best possible tissue recognition. One of them is the choice of the most relevant texture parameters. Such a choice must consider the variability of image acquisition settings that could result in different image properties, like those depending on slice thickness. Some ideas for adapting the system for working with images of different slice thicknesses are presented in the next part of our study.

Parameter Stability

The estimation of parameter stability could help to decide if the parameter is reliable for proper tissue characterization. Parameters sensitive to small changes in ROI size, or small ROI displacements, should be excluded from further analyses.

We adopt the approach considered in [START_REF] Lefebvre | Computerized ultrasound B-scan characterization of breast nodules[END_REF] in order to assess how the parameter changes over the different ROI locations or sizes. As a measure of its changeability, the classical coefficient of variation (CV, the ratio of the standard deviation to the mean) is used. Stable parameters are characterized by low CV values, and the more unstable is the parameter, the greater is its CV. In our work, the CV of each texture parameter will be calculated for different slice thicknesses.

We propose to apply the following approaches: Displace and Size Changing. In the Displace approach, the initial ROI dimensions are first slightly decreased, then the reduced ROI is displaced in order to take all the possible positions inside its initial boundaries. On the basis of each new ROI location, the value of a parameter is calculated. With the Size Changing approach, the initial ROI size is successively reduced, pixel by pixel, by moving each time one of the subsequent ROI vertices. For each of thus obtained ROIs a parameter value is calculated. In both cases, a set of several parameter values, obtained for different ROI locations (Displace) or sizes (Size Changing) serves to calculate a CV.

Effect of Slice Thickness on Classification Accuracy

In order to evaluate the effect of slice thickness on the classification accuracy, we propose to perform several experiments, each time considering a different combination of slice thicknesses for training of classifier and for its test. If the image of a particular slice thickness is included in a training set, its counterpart (of the same slice position in a patient's body) of another slice thickness should not be included in the test set.

Effect of Slice Thickness on Parameter Values

To access the parameter dependency of the slice thickness, we propose a similar approach that we apply for the assessment of the parameter stability. For a given parameter, its dependency on slice thickness will be measured by a "variation" between several parameter values measured by the classical coefficient of variation. These values will be obtained for the same ROI position on images of the same acquisition moment, but characterized by different slice thicknesses.

In order to validate the usefulness of the slice thickness-independent parameters, the experiments on images with different slice thicknesses (different for a training and for a testing stage) will be performed.

Experiments

Database Description

The images, from 29 patients, were gathered at the Department of Radiology of the Pontchaillou University Hospital in Rennes, France. They were acquired with LightSpeed16 device (GE Medical Systems). For each patient, three scan series were performed: without contrast, at arterial and at portal phase. The contrast material, of 100 ml, was injected at 4 ml/s, in an arm vein. The arterial phase acquisitions started about 20 seconds after the contrast product injection, the portal phase acquisitions started from 30 to 40 seconds later. For each of the three scan series, two "versions", corresponding to slice thicknesses of 5 mm and of 1.3 mm, were available. Each thick-slice image had its equivalent in thin-slice one, with the same slice position in the patient's body. All images had the size of 512×512 pixels. They were recorded in DICOM format, with 4096 gray levels. Since only the range of 248 gray levels sufficed for describing the pixels in the considered ROIs, the images were converted to a 8-bit BMP format.

The 303 pairs of images (of a thick slice and of a corresponding thin slice) were taken into account for each of the three acquisition moments. Two classes of liver tissue were represented: cirrhotic and healthy liver (171 and 132 pairs of image triplets, respectively). A square ROI of 60×60 pixels was drawn at the same location on each of the 2•3=6 considered simultaneously images. An example of six corresponding ROIs is given in Table 1.

Texture Parameters Chosen for Evaluation

In total, 155 texture parameters, extracted with 9 methods, were tested (see Table 2). They are based on: First Order Statistics (FO), Gradients (GB), Co-Occurrence Matrices (COM) [START_REF] Haralick | Textural features for image classification[END_REF], Run Length Matrices (RLM) [START_REF] Galloway | Texture analysis using gray level run lengths[END_REF], Gray Level Difference Matrices (GLDM) [START_REF] Weszka | A comparative study of texture measures for terrain classification[END_REF], Laws Texture Energy (LTE) [START_REF] Laws | Textured image segmentation[END_REF], Fractals (FB) [START_REF] Chen | Fractal feature analysis and classification in medical imaging[END_REF], Texture Feature Numbers (TFN) [START_REF] Horng | Texture feature coding method for classification of liver sonography[END_REF], and Autocorrelation (AC) [START_REF] Gonzalez | Digital Image Processing[END_REF].

When applying the COM, GLDM, and RLM methods, the number of gray levels was reduced from 256, used initially, to 64. The CO Matrices and the GLD Matrices were constructed separately for 4 standard directions (0 • , 45 • , 90 The normalized autocorrelation coefficients (AC) and the two TFN parameters (among 7) were calculated separately for 5 different pixel distances, from 1 to 5. They were included, respectively, in the sets AC 5 , and TFN 15 (with 5 other distance-independent TFN parameters). In the set TFN 7 , the values of 2 parameters calculated separately for the 5 distances were averaged.

The LTE method provided two parameter sets: LTE 14 and LTE 5 . The first one, composed of 14 parameters, was obtained by the application of 24 filtering masks of size 5×5: 4 symmetric, and 10 pairs of asymmetric ones, each pair consisted of a mask and its transposition. The second set was composed of 5 parameters, corresponding to the application of 3×3 masks, 2 symmetric and 3 pairs of asymmetric ones. The sum of elements of each convolution matrix was zero. For each pair of asymmetric masks, the resulting images were added. Images obtained by the application of symmetric masks were multiplied by two. Finally, the entropies of thus obtained images served as the texture parameters.

The FB method is based on the fractional Brownian motion model [START_REF] Chen | An automatic diagnostic system for CT liver image classification[END_REF] and considers 4 pixel distances [START_REF] Haralick | Statistical and structural approaches to texture[END_REF][START_REF] Guggenbuhl | Reproducibility of CT-based bone texture parameters of cancellous calf bone samples: Influence of slice thickness[END_REF][START_REF] Mayerhoefer | Effects of MRI acquisition parameter variations and protocol heterogeneity on the results of texture analysis and pattern discrimination: an application-oriented study[END_REF][START_REF] Miles | Colorectal cancer: texture analysis of portal phase hepatic CT images as a potential marker of survival[END_REF]. It provides a set of two parameters: FB 2 .

Effect of Slice Thickness on Parameter Stability

In this experiment, two approaches were applied for the CV calculation: Displace and Size Changing. The coefficient of variation was calculated on the basis of 9 parameter values. In the Displace approach, the ROI was reduced to a 58×58 square in order to take the 9 possible positions inside its initial boundaries. With the Size Changing one, the reduced ROI sizes were going from 60×60 to 52×52.

The experiment was performed separately for 12 texture sets, corresponding to all combinations of: 2 slice thicknesses, 2 tissue classes, 3 acquisition moments. For each parameter, the CV values obtained from the same texture set were averaged. We observed that regardless of the approach, the three corresponding averaged CV values, obtained for three acquisition moments, did not differ sig-nificantly. Nor did they differ between 2 tissue classes. The Displace approach almost always produced slightly higher CVs than the Size Changing approach. Table 3 shows the averaged CV values obtained for the two different slice thicknesses. Each value is the maximum one of the 2•2•3=12 values, obtained for each combination of the approaches, classes and acquisition moments. Since the presentation of all the results occupies too much space, we present the results for 16 selected parameters. However, the CVs calculated on all 155 parameters yield the same conclusions as those obtained from the presented subset.

We can conclude, that the parameter stability, expressed by its coefficient of variation, is not considerably influenced by the slice thickness. Regardless of the extraction method, the CVs corresponding to thick and thin slices were similar. None of the thicknesses had proven to ensure better parameter stabilities.

For further experiments we decided to use the most stable parameters, for which the coefficient of variation does not exceed 0.01. The sets of parameters satisfying this condition were nearly identical for the two slice thicknesses. Only 4 parameters turned out relatively stable with only one of the thicknesses: CodeVar (rejected for thin slices, with CV=0.0102), Contrast, DiffVar, and Hom (rejected for thick slices, with CV equal to 0.0144, 0.0195, and 0.1721, respectively). In total, 93 parameters were accepted for both slice thicknesses (see bolded parameters in Table 2), and the 62 parameters were rejected.

Effect of Slice Thickness on Classification Accuracy

In this step, we considered the following combinations of slice thicknesses for training and for testing of the classifier:

-"T/T": training and testing performed on only thick slices (of 5 mm), -"t/t": training and testing performed on only thin slices (of 1.3 mm), -"T/t": training on only thick slices, testing on only thin slices, -"t/T": training on only thin slices, testing on only thick slices, -"mix": both thick and thin slices in a training and a testing set.

For each combination, the whole set of 303 pairs of image triplets (thick and thin version) was randomly divided into two subsets: a training set (202 pairs) and a testing set (101 pairs). For each subset, the original proportion between tissue classes was preserved. The experiment was repeated 10 times.

Table 4 presents the classification results obtained for the 5 combinations of slice thicknesses used for training and testing. Only the most efficient stable parameter sets are taken into account. Classification was performed with the Weka software [START_REF] Hall | The WEKA data mining software: an update[END_REF], using Ensemble of Classifiers with adaptive boosting voting scheme [START_REF] Freund | A decision-theoretic generalization of online learning and an application to boosting[END_REF] (AdaBoostM1 ) and a C4.5 tree [START_REF] Quinlan | C4.5: Programs for Machine Learning[END_REF] as the underlying algorithm.

We can observe that using the same slice thickness, both for training and for testing the classifier, results in a quite high classification accuracy. The best results are more frequent for thick slice thicknesses ("T/T" possibility). In this case, the maximal classification accuracy (86.73%) is obtained for the parameter sets COM 55 and RLM 8 . Just below the best results are those obtained when only Unsatisfactory results are obtained when the slice thicknesses used for training and testing are different. Most of them do not exceed 60%. The best results for "T/t" and "t/T" combinations are: 64.16% and 58.82%, respectively.

We can conclude that, one should be careful using the system to aid the diagnosis, basing on images of a slice thickness different from those used for training.

However, we suppose that, in this case, the use of the texture parameters independent of the slice thickness could improve the system performance. The search for such parameters will be the subject of our next experiment.

Effect of Slice Thickness on Parameter Values

In this experiment, we considered only those parameters that have been found to be stable for both slice thicknesses: 5 mm and 1.3 mm (see bolded parameters in Table 2). The parameter dependency from the slice thickness was expressed by the coefficient of variation of its two values, obtained from two corresponding ROIs (drawn on the images of a 5 mm and of a 1.3 mm slice thickness).

In total, 303 pairs of images were considered for each of 3 acquisition moments. For each parameter, the 3•303 values of the coefficient of variation were averaged. The ranking of parameters, according to the average coefficient of variation, was performed. Parameters with the lowest average coefficient are considered as least dependent on slice thickness. Table 5 presents the ranking of the first 40 parameters which are least dependent on the slice thickness.

Since we do not know how many parameters are acceptable and sufficient for the best possible tissue recognition, when different slice thicknesses are considered, we will test several sets of the first parameters from the ranking.

Classification of Textures Characterized by Parameters the Least Dependent on Slice Thickness

In our final experiment, 7 parameter sets were tested: First-8, First-13, First-16, First-23, First-30, and First-39. They were composed, respectively, of the first 8, 13, 16, 23, 30, and 39 parameters from the ranking presented in Table 5. Table 6 presents the classification accuracy obtained for those 7 sets by Ensemble of Classifiers (with the same settings as in the previous classification experiment). We can conclude that the classification accuracy, obtained for the 5 considered combinations of slice thicknesses are similar only when the First-8 parameter set is used. This set provides the classification accuracy ranging from 84.06%, for the "t/T" combination, to 87.13%, for the "mix" one. We can also notice that the application of the First-8 set guaranties the best classification results for the combinations, for which the slice thickness used for testing is different from that used for training (85.35% for "T/t" and 84.06% for "t/T"). For the first of these combinations, the best result is now 21.39% better than the best one obtained in the previous classification experiment, which did not consider the slice thickness-independent parameters (see Tab. 4). For the second onethe best classification accuracy augmented of 25.24%. Nevertheless, this accuracy still remain of about 6%-8% worse than the best results obtained now for the one-thickness cases ("T/T", "t/t").

In general, for the combinations "T/t" and "t/T", the more numerous is the set of parameters, the worse is the classification accuracy. Contrarily, for the one-thickness cases, the classification is generally better for more numerous sets.

With the "mix" combination, the application of the First-8 set leads to obtaining quite good classification accuracy, but not the best. In this case, increasing the number of parameters first results in lowering the classification accuracy (82.58% for the First-13 parameter set), then leads to its successive increase. Finally, the best classification result, for "mix" case (88.32%) is obtained by the most numerous parameter set, First-39, and is 9.60% better than the best accuracy observed, for "mix" case, in the previous classification experiment.

Conclusions and Future Work

The experiments allowed us to: (i ) analyze the influence of the slice thickness on parameter stability, (ii ) assess the parameter dependency on slice thickness, (iii ) find parameters which are the least dependent on slice thickness, (iv ) evaluate the classification accuracy obtained by parameters which are less and more dependent on slice thickness, when different slice thicknesses were simultaneously considered. Several conclusions can be formulated to sum up our study.

First, the parameter stability does not considerably depend on slice thickness. The sets of parameters recognized as stable were nearly identical for the two slice thicknesses. Second, one should be particularly careful when applying a CAD system, when recognized images are of the slice thickness different from the one used for the system training. With quite popular parameters, obtained by the COM or RLM methods, a satisfactory tissue recognition is not possible. In this case, it would be safer to use the slice thickness-independent parameters.

In the future, we plane to perform the experiments with more slice thicknesses. Due to the fact that acquiring the series of images of many different slice thicknesses is practically impossible within a single patient study, we plane to use a phantom. It will enable to assess not only the effect of the slice thickness in the texture-based classification, but also the effect of several other parameters, like image resolution, or the time elapsed from the contrast product injection.

Fig. 1 .Fig. 2 .

 12 Fig. 1. System for texture-based classification of liver tissues: Training

  135 • ) and for 5 different distances between the pixel pairs, going from 1 to 5. From each of 20 thus obtained matrices, the same parameters were calculated, 11 parameters by the COM method and 5 parameters by GLDM. The RL Matrices were also constructed for 4 standard directions, each of them served to calculate 8 parameters. For the three aforementioned methods, an averaging of the same parameter corresponding to 4 different directions was done. As the result, the following parameter sets were obtained: COM 55 (11•5 parameters), GLDM 25 (5•5 parameters), and RLM 8 . The sets COM 11 and GLDM 5 are the result of averaging of 20 parameter values, calculated for 4 directions and for 5 distances.

5 GLDM5 5 TFN7

 55 , (d )AngSecMom, (d )SumVar, (d )DiffVar, (d )DiffAvg, (d )Contrast, (d )Corr, where d = 1DAvg, (d )DEntr, (d )DAngSecMom, (d )DInvDiffMom, (d )DContrast, where d = 1, 2, 3, 4, CodeEntr, Coarse, Hom, CodeVar, ResSim, (d )CodeSim, where d = 1, 2, 3, 4, MeanConv, CodeEntr, Coarse, Hom, CodeVar, ResSim, CodeSim

Table 1 .

 1 Six ROIs at the same slice position

	Acquisition:	Without Contrast	Arterial Phase	Portal Phase
	ROI example:

• ,

Table 2 .

 2 Texture parameters chosen for evaluation. In bold -parameters which later turned out to be stable for both slice thicknesses

Table 3 .

 3 Maximum

		parameter CV value (among 12 values, corresponding to: 2 ap-
	proaches, 2 classes and 3 acquisition moments), obtained for 2 different slice thicknesses
	Parameter 5 mm 1.3 mm	Parameter	5 mm 1.3 mm
	(1)Autocorr 0.0004 0.0007	HighGLREmp 0.0035 0.0037
	FractalDim 0.0012 0.0009	(4)CodeEntr 0.0036 0.0027
	(5)SumAvg 0.0016 0.0017	S3E3	0.0038 0.0032
	(1)SumAvg 0.0017 0.0018	(1)InvDiffMom 0.0039 0.0048
	Avg	0.0018 0.0018	(3)AngSecMom 0.0047 0.0046
	E3L3	0.0023 0.0021	GradAvg	0.0054 0.0052
	DiffEntr	0.0026 0.0019	MeanConv	0.0058 0.0062
	S5E5	0.0034 0.0031	DAvg	0.0063 0.0057

Table 4 .

 4 Classification accuracy obtained by Ensemble of Classifiers for 5 combinations of slice thicknesses, used for training and testing. Each line corresponds to a different set of texture parameters. Only stable parameters are considered

	Set	T/T	t/t	T/t	t/T	mix
	COM55 86.73±2.48 84.75±2.73 63.07±2.24 56.64±5.52 78.72±4.81
	COM11 84.16±3.33 82.97±3.29 64.16±2.08 53.77±5.57 77.13±4.26
	GLDM25 67.03±2.47 69.81±6.40 56.84±1.76 48.02±2.61 57.93±1.94
	GLDM4 63.47±3.51 61.09±4.55 56.94±0.70 49.21±6.32 56.84±0.96
	LTE14 78.02±3.45 78.62±2.65 58.42±2.95 44.75±8.87 69.81±3.62
	LTE5	74.46±3.67 74.46±4.89 59.71±7.61 49.50±4.48 63.27±5.25
	RLM8 86.73±4.43 85.15±1.81 63.96±3.65 58.82±6.65 78.62±3.65

Table 5 .

 5 Ranking of parameters by their dependency on the slice thickness 15% with the RLM 8 set. Slightly inferior results are observed when both slice thicknesses are considered for training and for testing ("mix"). The best of them is 78.72%, with the set COM 55 .

	Rank	Parameter	Variation	Rank Parameter Variation
	1	Avg	0.0024	21 (2)CodeEntr 0.0515
	2	(5)SumAvg	0.0024	22 (3)CodeEntr 0.0515
	3	(4)SumAvg	0.0024	23	RLEntr	0.0658
	4	SumAvg	0.0024	24	E3E3	0.0750
	5	(3)SumAvg	0.0024	25	E3L3	0.0787
	6	(2)SumAvg	0.0024	26	E5E5	0.0803
	7	(1)SumAvg	0.0025	27	S5S5	0.0811
	8	FractalDim	0.0027	28	S5L5	0.0826
	9	(1)Autocorr	0.0104	29	E5L5	0.0847
	10	(4)Autocorr	0.0144	30	S5E5	0.0884
	11	(5)Autocorr	0.0145	31	(1)SumEntr 0.0905
	12	(3)Autocorr	0.0147	32	SumEntr	0.0931
	13	(2)Autocorr	0.0155	33	(2)SumEntr 0.0932
	14 HighGLREmp 0.0157	34	S3L3	0.0935
	15	ShortEmp	0.0296	35	(3)SumEntr 0.0937
	16	Fraction	0.0394	36	(4)SumEntr 0.0939
	17	(1)CodeEntr	0.0511	37	W5E5	0.0941
	18	CodeEntr	0.0514	38	(5)SumEntr 0.0945
	19	(5)CodeEntr	0.0514	39	W5L5	0.0953
	20	(4)CodeEntr	0.0515	40	S3E3	0.1016
	thin slices are considered ("t/t"): 85.

Table 6 .

 6 Classification accuracy obtained by Ensemble of Classifiers for 5 combinations of slice thicknesses, used for training and testing. Each line corresponds to a different set of stable texture parameters, the least dependent of the slice thickness

	Set	T/T	t/t	T/t	t/T	mix
	First-8 85.94±4.12 86.64±2.77 85.35±4.78 84.06±4.08 87.13±3.90
	First-13 87.43±2.64 87.53±1.76 67.03±3.59 75.94±3.62 82.58±3.10
	First-16 87.82±2.84 87.63±3.17 63.47±2.11 73.96±6.42 85.05±3.21
	First-23 89.51±3.50 88.42±2.76 62.38±3.10 74.46±4.82 86.24±3.38
	First-30 91.39±3.10 88.52±6.00 64.06±2.92 71.09±5.00 88.22±2.49
	First-39 90.79±3.93 92.48±2.10 62.78±1.88 66.74±6.78 88.32±2.71

Acknowledgments

We thank Dr D. Olivie for his precious help. This work was supported by the grant S/WI/2/2013 from Bialystok University of Technology.


Izabela Rejer
email: irejer@wi.zut.edu.pl

Pawel Górski
email: pagorski@wi.zut.edu.pl

Paweł Górski

Independent Component Analysis for EEG Data Preprocessing -Algorithms Comparison

Keywords: Brain Computer Interface, BCI, EEG, preprocessing, ICA, Independent Component Analysis

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

Introduction

A BCI (Brain Computer Interface) is defined as a communication system in which messages or commands that a user sends to the external world do not pass through the brain's normal output pathways of peripheral nerves and muscles [START_REF] Wolpaw | Braincomputer interfaces for communication and control[END_REF]. Although, at first brain computer interfaces were dedicated mainly to people in locked-in state or in completely locked-in state, nowadays, the range of their potential recipients is much wider. They are tested and used by the army, they are used in the entertainment industry and there are also some attempts to use them to control limb prosthesis.

A BCI system is composed of seven parts, followed one after the other in a closedloop. Succeeding parts of the BCI system are responsible for: measuring brain activity, preprocessing of the acquired signals, describing the signals by a few relevant features (feature extraction), selecting the most relevant features (feature selection), assigning a class to a set of selected features (classification), executing a command assigned to the chosen class and providing feedback to the user informing him about the mental state recognized by the BCI system. Although, all steps of the loop have to be carefully designed in order to build a successful BCI system, four of them play the major role: preprocessing [START_REF] Wang | Common Spatial Pattern Method for Channel Selection in Motor Imagery Based Brain-Computer Interface[END_REF][START_REF] Oja | The FastICA Algorithm Revisited: Convergence Analysis[END_REF][START_REF] Augustyniak | Adaptive Wavelet Discrimination of Muscular Noise in the ECG[END_REF], feature extraction [START_REF] Yang | Feature Extraction of Motor Imagery EEG Based on Wavelet Transform and Higher-Order Statistics[END_REF][START_REF] Müller-Putz | Fast set-up asynchronous brain-switch based on detection of foot motor imagery in 1-channel EEG[END_REF], feature selection [START_REF] Koprinska | Feature Selection for Brain-Computer Interfaces[END_REF][START_REF] Rejer | Genetic Algorithms in EEG Feature Selection for the Classification of Movements of the Left and Right Hand[END_REF][START_REF] Burduk | Recognition Task with Feature Selection and Weighted Majority Voting Based on Interval-Valued Fuzzy Sets in Computational Collective Intelligence[END_REF] and classification [START_REF] Pfurtscheller | Separability of EEG Singals Recorded During Right and Left Motor Imagery Using Adaptive Autoregressive Parameters[END_REF][START_REF] Bobrowski | Separable Linearization of Learning Sets by Ranked Layer of Radial Binary Classifiers[END_REF][START_REF] Wozniak | Combined classifier based on feature space partitioning[END_REF]. This paper deals with the first of these four main steps, it is with the preprocessing step. The aim of this step is to transform the acquired brain signals to its purer form by eliminating artifacts, reducing the ongoing brain activity, enhancing the primary components. One of the methods used in the preprocessing step is Independent Component Analysis (ICA).

Independent Component Analysis is a method for transforming a set of mixed signals into a set of independent components. It has been reported in some publications that this method allows to detect artifacts in originally recorded signals [START_REF] Chen | A symmetric orthogonal FastICA algorithm and applications in EEG[END_REF]. Due to this, ICA seems to be a good choice for preprocessing of electroencephalographic data (EEG) used as control signals in brain computer interfaces (BCI). After applying ICA on a set of EEG data, some components should reflect original data sources and one or more components should reflect artifacts. In the process of feature selection, which is a further step of signal processing leading to classification, features calculated for the "artifact components" should be discarded from the set of features while features calculated for the components essential for the classification precision should be preserved.

One can ask why features calculated for the "artifact components" will be discarded from the feature set. The answer is quite straightforward. There are two main types of artifacts activity. Some artifacts, such as eye artifact, appear rhythmically during each trial, while others, such as unexpected body movements, appear in unexpected moments along the whole experiment. Since artifacts from both groups do not dependent on the class coded in the recorded signals, they do not enhance the classification precision. Hence, when the classification precision is a measure used in the feature selection process, features calculated for the components reflecting artifacts should be discarded from the feature set.

There are a lot of scientific papers describing applications of ICA in BCI research [START_REF] Chen | A symmetric orthogonal FastICA algorithm and applications in EEG[END_REF][START_REF] Peterson | Feature Selection and Blind Source Separation in an EEG-Based Brain-Computer Interface[END_REF][START_REF] Boord | Discrimination of left and right leg motor imagery for brain-computer interfaces[END_REF]. In most of them there is underlined that by applying ICA for data preprocessing, the classification precision increases. However, these papers do not specify the exact profits (in terms of classification accuracy) gained by using different algorithms for calculating independent components. Does it mean that the choice of the ICA algorithm is of no importance for the classification precision?

The aim of this paper is to examine whether ICA transformation, performed according to three most popular algorithms, has a major or only minor influence on the classification precision in different experimental settings. Hence, three questions are posed in the paper:

 Is the classification precision essentially higher when ICA transformation is performed?  Has a choice of a specific algorithm for performing ICA an influence on the classification precision?  Is there any correlation between the classification precision (after performing ICA) and the number of input features used in the classification process?

In order to answer these questions, a set of experiments was planned and performed. In all experiments a data set from the second BCI Competition was used. The results of the experiments, together with a short discussion, are presented in the paper.

Independent Component Analysis

The problem of a blind source separation (BSS) consist in finding a matrix W such that the linear transformation will allow to recover the source signals from a set of mixed signals [START_REF] Park | A modified infomax algorithm for blind signal separation[END_REF][START_REF] Stone | Independent component analysis: an introduction[END_REF]. The term 'blind' means that no prior information about the source signals or the mixing process is available [START_REF] Park | A modified infomax algorithm for blind signal separation[END_REF].

Independent Component Analysis (ICA) is one of the most popular BSS method. ICA problem can be stated as follows. Let's assume that there are n linear mixtures of n independent components. Vector x (observed signals) can be written as:

x = As (1)
where A represents a mixing matrix with the size of , and s is the vector of independent components. The aim of ICA is to find a matrix W (i.e. an inverse of the matrix A) to reverse the mixing effect. Then, after computing the matrix W, we can obtain the independent components by [START_REF] Langlois | An Introduction to Independent Component Analysis: InfoMax and FastICA algorithms[END_REF][START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF]:

(2) Most of the popular ICA algorithms put some constraints on the mixed signals. First of them is a statistical independence between source signals ; second, a non-Gaussian distribution of the source signals and the third -the equality of the number of source signals and the number of mixture signals. While two first constrains are main assumptions utilized by many algorithms, the third one is introduced only to decrease the algorithm complexity (it causes that the mixing matrix is square). Furthermore, it is assumed that each source signal has the unit variance { } . To hold this assumption, the matrix of the source signals is whitened before the ICA calculation [START_REF] Langlois | An Introduction to Independent Component Analysis: InfoMax and FastICA algorithms[END_REF][START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF]. One more assumption, introduced only to simplify the algorithm, is that all mixture signals are centered.

As was mentioned earlier, ICA does not require any prior information about the source signals. Instead, ICA algorithms utilize the concept of statistical independency of the mixed signals. According to the formal definition, the variables a and b are said to be independent if information about the value a does not give any information about the value b and vice versa [START_REF] Stone | Independent component analysis: an introduction[END_REF], [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF]. Technically, independence can be defined in terms of the probability density function (pdf) [START_REF] Langlois | An Introduction to Independent Component Analysis: InfoMax and FastICA algorithms[END_REF]:

( ) ( ) ( ) ( ) (3)
where are random variables. There are two main approaches to measuring independence: maximization of non-Gaussianity and minimization of mutual information. Most of the existing ICA algorithms are based on one of them. When the first approach is applied, the task for the algorithm is to modify the components in such a way to obtain the source signals of strong non-Gaussian distribution (the assumption is: the stronger non-Gaussianity, the stronger independence [START_REF] Langlois | An Introduction to Independent Component Analysis: InfoMax and FastICA algorithms[END_REF]). In other words, the distributions of the mixture signals have to be more Gaussian than the source signals. This approach utilizes various measures of non-Gaussianity, like: kurtosis, negentropy, approximations of negentropy and others [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF].

Mutual information, utilized in the second approach, informs how much information about the variable a can be gained from the information about the variable b. Since smaller value of mutual information means that more information about a given system is stored in the variables [START_REF] Langlois | An Introduction to Independent Component Analysis: InfoMax and FastICA algorithms[END_REF], ICA algorithms based on mutual information approach minimize the mutual information of the system outputs [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF].

FastICA -Deflation Approach

The FastICA algorithm, proposed by Hyvärinen and Oja, is an iterative method to find local maxima of a defined cost function [START_REF] Langlois | An Introduction to Independent Component Analysis: InfoMax and FastICA algorithms[END_REF][START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF], [START_REF] Oja | The FastICA Algorithm Revisited: Convergence Analysis[END_REF]. The purpose of this algorithm is to find the matrix of weights such that the projection ( ) maximizes non-Gaussianity [START_REF] Oja | The FastICA Algorithm Revisited: Convergence Analysis[END_REF], [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF]. As a measure for non-Gaussianity, simple estimation of negentropy based on the maximum entropy principle is used [START_REF] Langlois | An Introduction to Independent Component Analysis: InfoMax and FastICA algorithms[END_REF][START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF]:

( ) { ( )} { ( )} (4)
where:standardized non-Gaussian random variable,standardized random variable with Gaussian distribution, ( ) -any non-quadratic function.

There are two classes of FastICA algorithms, the deflation algorithms (called also one-unit algorithms) and the symmetric algorithms [START_REF] Tichavský | Performance Analysis of the FastICA Algorithm and Cramér-Rao Bounds for Linear Independent Component Analysis[END_REF]. In the deflation approach, the independent components (ICs) are extracted sequentially, one by one. The algorithm can be summarized as follows [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF], [START_REF] Naik | An Overview of Independent Component Analysis and Its Applications[END_REF]:

1. Choose an initial vector w (e.g. random) 2. Do steps 3-6 3.

{ ( )} { ( )}

4.

|| ||

5. Do the Gram-Schmidt orthogonalization:

∑ √ 6. Stop if not converged
Gram-Schmidt procedure, used in the algorithm, prevents different vectors from matrix w from converging to the same maxima [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF]. The order, in which the independent components are extracted, depends on the initial value of .

Independent Component Analysis for EEG Data Preprocessing -Algorithms Comparison

FastICA -Symmetric Approach

The only difference between deflation approach and symmetric approach is the procedure of weights calculation. While in deflation approach vectors of weights are calculated one by one, in symmetric approach the estimation of all components (all weights vectors) proceeds in parallel [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF][START_REF] Tichavský | Performance Analysis of the FastICA Algorithm and Cramér-Rao Bounds for Linear Independent Component Analysis[END_REF]. Instead of Gram-Schmidt procedure, the following formula is used in the orthogonalization step:

( ) (5)
where is the matrix of weights vectors ( ) . The square root of ww T is obtained from the eigenvalue decomposition of as [START_REF] Fan | A New FastICA Algorithm with Symmetric Orthogonalization[END_REF]:

( ) (6)
where is the matrix of eigenvectors and is the diagonal matrix of eigenvalues.

The algorithm is performed until the stop condition (e.g. given by 7 [START_REF] Tichavský | Performance Analysis of the FastICA Algorithm and Cramér-Rao Bounds for Linear Independent Component Analysis[END_REF]) is met:

( ( ( ))) (7)
where is a chosen constant.

Infomax

Infomax algorithm is based on the general optimization principle for neural networks and other processing systems described by Linsker in 1987 [START_REF] Linsker | Self-Organization in a Perceptual Network[END_REF]. In general this principle says that a function that maps a set of input values a to a set of output values b should be chosen or learned so as to maximize the average Shannon mutual information between a and b. The ICA algorithm utilizing this principle was first proposed in 1995 by Bell and Sejnowski [START_REF] Bell | An information-maximisation approach to blind separation and blind deconvolution[END_REF] and then in 1997 optimized by Amari [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF], [START_REF] Naik | An Overview of Independent Component Analysis and Its Applications[END_REF].

Infomax algorithm for calculating independent components is based on the maximization of the output entropy of a neural network with non-linear outputs [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF]. The most essential parameter of this algorithm is a learning rate which does not need to be constant over time and which should give a good compromise between speed of learning and estimation precision [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF], [START_REF] Raimondo | CUDAICA: GPU Optimization of Infomax-ICA EEG Analysis[END_REF]. The weights of this neural network are updated according to the following formula [START_REF] Langlois | An Introduction to Independent Component Analysis: InfoMax and FastICA algorithms[END_REF], [START_REF] Naik | An Overview of Independent Component Analysis and Its Applications[END_REF], [START_REF] Karhunen | Neural approaches to independent component analysis and source separation[END_REF]:

( ) (8)
where: ymatrix of source estimation (y=Wx);number of iteration;the identity matrix; learning rate which may depend on k; ( )a nonlinear function.

Mostly a classic logistic function is used as a nonlinear function g [START_REF] Karhunen | Neural approaches to independent component analysis and source separation[END_REF]:

( ) , (9)
however, sometimes also its extended version is applied:

( ) ( ) (10)
Using ( 9), the Infomax algorithm can be summarized as follows [START_REF] Naik | An Overview of Independent Component Analysis and Its Applications[END_REF]:

1. x = perm(sources); 2. 3. 4. 5. = 6.
where is random permutation.

Experimental Settings

The comparison of ICA algorithms described in Section 2 was carried out with a data set submitted to the second BCI Competition (data set IIImotor imaginary) by Department of Medical Informatics, Institute for Biomedical Engineering, Graz University of Technology [START_REF]Data set III, II BCI Competition, motor imaginary[END_REF]. The data set was recorded from a normal subject (female, 25y) whose task was to control the movements of a feedback bar by means of imagery movements of the left and right hand. Cues informing about the direction in which the feedback bar should be moved were displayed on a screen in the form of the left and right arrows. The order of left and right cues was random. The experiment consisted of 280 trials, each trial lasted 9 seconds. The first 2s was quiet, at t=2s an acoustic stimulus was generated and a cross "+" was displayed for 1s; then at t=3s, an arrow (left or right) was displayed as a cue. The EEG signals were measured over three bipolar EEG channels (C3, Cz and C4), sampled with 128Hz and preliminary filtered between 0.5 and 30Hz. The whole data set, containing data from 280 trials, was then divided into two equal subsetsthe first one was intended for classifier training and the second intended for external classifier test. Since only data from the first subset was published with target values (1 -left hand, 2-right hand), only this subset could be used in the research.

In the preprocessing step, the data from the original data set was transformed according to the algorithms described in Section 2. After performing this step, three different sets of signals were obtained:

1. Set of components obtained with FastICA -deflation approach algorithm. 2. Set of components obtained with FastICA -symmetric approach algorithm.

Set of components obtained with Infomax algorithm.

All these sets of components, together with the fourth set, composed of original signals from the channels C3, Cz and C4 were used in the experiments.

The data from each set of components, was transformed to a set of frequency band power features. The signal power was calculated separately for:

1. 12 frequency bands: alpha band (8-13Hz) and five sub-bands of alpha band (8-9Hz; 9-10Hz; 10-11Hz; 11-12Hz; 12-13Hz); beta band (13-30Hz) and also five sub-bands of beta band (13-17Hz; 17-20Hz; 20-23Hz; 23-26Hz; 26-30Hz), 2. each of 7 seconds of the trial (data from the first and second seconds of the recordings of each trial were discarded because they covered the period before the clue presentation), 3. each of 3 canals (C3, Cz, C4).

In this way 252 band power features were obtained per each of four sets of components. Taking into account a very small number of trials equal to 140, the number of features had to be significantly reduced before the classification step. In fact with 140 trials, no more than several features should be used without the threat of overfitting. According to Raudys and Jain, at least 10 times more training data per class than the features should be gathered to train the classifier correctly [START_REF] Raudys | Small sample size effects in statistical pattern recognition: Recommendations for practitioners[END_REF].

In order to reduce the number of features, a genetic algorithm, described in details in [START_REF] Rejer | Genetic Algorithms in EEG Feature Selection for the Classification of Movements of the Left and Right Hand[END_REF] was used. Some basic features of this algorithm are as follows:

1. An individual is composed of the number of genes equal to the critical number of features, given by the user or calculated automatically in terms of number of observations, number of classes and classifier type. 2. Each gene can take an integer value from the interval {0,1...F}, where F denotes the dimension of the feature set. 3. The basic genetic operation in the algorithm is a very aggressive mutation. The aggressive mutation means that not only each individual in the population is mutated, but also each gene of each individual. 4. Since, after the mutation a lot of new individuals is born, the selection step is performed after the reproduction step. The selection is made from the population composed of parent individuals and their mutated children. 5. The fitness function is pure classifier accuracy.

As it was stated above, a fitness function of the algorithm used for feature selection was a classification accuracy. A classifier was built per each individual of each generation. Input features introduced to each classifier were encoded in succeeding genes of the evaluated individual.

A linear SVM method was used in the classification process. The classification threshold was set to 0.5 and hence, all classifier results greater than 0.5 were classified as class "2" (right hand) and results smaller or equal to 0.5 were classified as class "1" (left hand). The classifiers accuracy was tested with 10-fold crossvalidation. The final accuracy measure of a given feature set was the mean value calculated on the basis of classification accuracy obtained for all validation sets. The accuracy of one validation set was calculated according to the following equation: [START_REF] Bobrowski | Separable Linearization of Learning Sets by Ranked Layer of Radial Binary Classifiers[END_REF] where: A k -accuracy of k validation subset (k=1...10), R k -number of properly classified cases from k validation subset, U k -number of all cases in k validation subset.

Results

The experiments were performed with the application prepared for the Matlab 7.12.0 environment. All analyzed preprocessing algorithms, it is: FastICA-deflation, FastI-CA-symmetric and Infomax, were implemented according to the general schemes given in Section 2. The learning rate in Infomax algorithm was variable in time. In three first iterations it was equal to: 0.01, 0.001 and 0.0001, respectively and in the remaining 500 iterations it was equal to 0.005. The main aim of the experiments described in the paper was to find out whether the preprocessing with ICA has an essential influence on the classification accuracy.

In order to answer this question, the classification accuracy calculated over raw signals and signals preprocessed with the analyzed ICA algorithms had to be compared. The comparison was carried out for different number of input features introduced to the classifier (from one to six). To perform the analysis, six genetic algorithms were prepared. The first one processed individuals composed of only one gene, second processed individuals composed of two genes and so forth up to the last one which processed individuals composed of six genes. Since a genetic algorithm is a heuristic optimization method which gives only sub-optimal solutions, the algorithm was run five times for each from the given six settings. Each from these 120 algorithms (six features, 4 sets of input signals, 5 runs) processed data by 30 generations. After each generation of each genetic algorithm, the individual of the highest value of the fitness function was stored in a table. At the end of the experiments, the results of the best individuals were averaged separately for each set of signals and for each number of input features. The average values of fitness function (classification accuracy) of the best individuals, together with standard deviations are presented in Table 1.

Discussion

In Section 1 three questions regarding the influence of ICA transformation on the classification precision were posed. To answer the first question, about the overall increase in the classification accuracy, the average value of the classification accuracy obtained with all three ICA algorithms was calculated. The result was surprising because the average classification precision calculated for signals preprocessed with ICA algorithms was equal to 87.2% which was slightly lower than the average classification precision calculated for raw signals (87.5%).

In order to answer the second question about the influence of the choice of a specific ICA algorithm on the classification precision, the comparison between algo-rithms had to be done. This time the results were more consistent with the theory and other analysis. The highest average classification precision for all six genetic algorithms was obtained with Infomax algorithm (87.6%), a slightly smaller precision was obtained with FastICA-symmetric algorithm (87.1) and the smallest precision was obtained with FastICA-deflation (86.9%) (Fig. 2). Looking closer at the average precisions of all three algorithms, one can easily notice that the differences between the performance of ICA algorithms are very small, in fact they are so small that they are not regarded as significant by any reasonable statistical test. Therefore, it is difficult to state that one algorithm is better than the others only on the basis of the classification precision. However, when standard deviations of results are taken into consideration, the conclusion is quite different. The standard deviation of the results gathered in Tab. 1 was on average by 33% lower in case of Infomax (0.8%) than in case of the raw signals (1.2%). This means that the Infomax algorithm gave more stable results and that these results are more reliable than results calculated over raw signals. Also the uncertainty of the classification precision calculated over signals preprocessed with two other ICA algorithms was significantly smaller than in case of raw data (std of FastICA-deflation -0.9% and std of FastICA-symmetric -0.9%). Hence, taking into account the standard deviation of results a conclusion of a practical usefulness of ICA transformation should be drawn. And finally, the third question posed in Section 1 was about a correlation between difference in the classification precision obtained with signals preprocessed with analysed ICA algorithms and the number of input features used in the classification process. This question is difficult to address because of very small, in fact non-significant differences in classification accuracy obtained over signals preprocessed with all three algorithms at each step of the experiment. On average Infomax gave better results but for different number of input features different algorithms exhibits slightly better performance: FastICA-deflation gave the best result for two features, FastICAsymmetric for five features and Infomax for one, three, four and six features.

Conclusion

The classification accuracy is a very important factor of a successful BCI because wrong classification results in executing improper action/command. Misclassification is not allowed, especially when the BCI is used as a system controlling actions of real devices.

The overall goal of this paper was to determine whether the preprocessing with ICA results in increasing the classification accuracy. Since the classification results were very similar with preprocessing and without it, it is very difficult to answer this question decidedly. Undoubtedly, the application of ICA increases the complexity of the whole BCI system so, on the one hand, taking into account only the lack of improvements in the accuracy, the answer should be "no". However, on the other hand, after applying ICA the uncertainty of classification results decreased rapidly which is a very important fact pro ICA application. Hence, before the question posed in the paper will be answer definitely, more experiments have to be done, experiments which will allow to find out which components survive the selection process, whether these components are the same in all ICA algorithms and lastly what is their distribution over the head model.

  signals preprocessed by one of the analyzed ICA algorithms (Fig. 1a -FastICAdeflation, Fig. 1b -FastICA-symmetric, Fig. 1c -Infomax).

Fig. 1 .

 1 Fig. 1. The average classification accuracy calculated over five runs for each of the six genetic algorithms. Each sub-figure presents the results obtained with raw signals and the results obtained with signals preprocessed by one of the analyzed ICA algorithms (Fig. 1a -FastICAdeflation, Fig. 1b -FastICA-symmetric, Fig. 1c -Infomax).

Fig. 2 .

 2 Fig. 2. The comparison of classification accuracy calculated over signals preprocessed with the analyzed ICA algorithms for different numbers of input features.

Table 1 .

 1 The average values of fitness function (classification accuracy) and standard deviation calculated over the best individuals obtained for different sets of signals and different numbers of input features. A much more compact comparison of results obtained with ICA algorithms and results obtained with raw signals is given in Fig. 1. Each sub-figure of Fig. 1 presents the comparison of the results obtained with raw signals with the results obtained with

		Raw data	Deflation	Symmetric	Infomax
	No. of features	Accuracy [%]	std	Accuracy [%]	std	Accuracy [%]	std	Accuracy [%]	std
	1	76.71	0.64	75.14	0.53	75.00	0.45	75.86	0.83
	2	82.57	2.7	82.71	1.31	82.43	1.32	82.29	0.29
	3	89.00	0.81	88.43	0.53	87.86	1.36	88.57	1.43
	4	91.00	1.08	90.86	0.53	90.71	0.9	91.71	0.73
	5	92.71	1.06	92.14	1.01	92.57	0.35	92.43	0.57
	6	93.14	1.08	92.29	1.53	93.86	0.73	94.57	0.73