P Riedinger

J.-C Vivalda

An LQ sub-optimal stabilizing feedback law for switched linear systems

The aim of this paper is the design of a stabilizing feedback law for continuous time linear switched system based on the optimization of a quadratic criterion. The main result provides a control Lyapunov function and a feedback switching law leading to sub optimal solutions. As the Lyapunov function defines a tight upper bound on the value function of the optimization problem, the sub optimality is guaranteed. Practically, the switching law is easy to apply and the design procedure is effective if there exists at least a controllable convex combination of the subsystems.

Introduction

Over the past decade, the design of stabilizing laws for switched systems (in continuous and discrete time) has been the focus of considerable research attention. Several approaches have been used to tackle this problem, one can cite for example [START_REF] Lincoln | Relaxing dynamic programming. Automatic Control[END_REF][START_REF] Seatzu | Optimal control of continuous-time switched affine systems[END_REF][START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF][START_REF] Zhang | On the value functions of the discrete-time switched lqr problem[END_REF] for dynamic programming approaches, [START_REF] Axelsson | Gradient descent approach to optimal mode scheduling in hybrid dynamical systems[END_REF][START_REF] Sorin | Optimal control of switching systems[END_REF][START_REF] Shaikh | On the hybrid optimal control problem: Theory and algorithms[END_REF] for variational approaches, or [START_REF] Colaneri | Stabilization of continuous-time switched nonlinear systems[END_REF][START_REF] Deaecto | Switched state-feedback control for continuous time-varying polytopic systems[END_REF][START_REF] José | Stability and stabilization of discrete time switched systems[END_REF] for Lyapunov based approaches. This problem is not easy, even numerically [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF][START_REF] Shaikh | On the hybrid optimal control problem: Theory and algorithms[END_REF] and the design of a stabilizing feedback law based on the optimization of a criterion is a challenging task.

LQ regulators are widely used for the control of linear systems because of their simple design and their robustness properties. These regulators can also been used for the 1 design of stabilizing feedback laws for linear switched systems but, as yet, one cannot obtain the exact solution of a switched LQ problem. Moreover, the main drawback of this method is the difficulty to get a good numerical approximation of the solution of the optimal problem; this approximation is difficult to obtain even for small dimensional systems. Another possibility is to use open loop control law, which can be achieved through direct or indirect methods [START_REF] Shaikh | On the hybrid optimal control problem: Theory and algorithms[END_REF][START_REF] Axelsson | Gradient descent approach to optimal mode scheduling in hybrid dynamical systems[END_REF], nevertheless, in this case singular solutions [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF][START_REF] Sorin | Optimal control of switching systems[END_REF] cause numerical complications [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF].

In [START_REF] José | Stability and stabilization of discrete time switched systems[END_REF], the author address two Lyapunov based strategies for stabilization of discrete time linear switched systems. The first one is of open loop nature while the second one is of closed-loop nature and is designed from the solution of the Lyapunov-Metzler inequalities. Their approach uses a family of quadratic Lyapunov functions and an upper bound on the cost is provided but the distance from the optimality of the stabilizing feedback law is not estimated.

In this paper, we consider a linear and controlled switched linear system. Together with this system, we consider a quadratic cost function. Our aim is to design a stabilizing feedback law that approaches the solution of the optimal problem related to this cost function. To the best of our knowledge, the problem of finding a continuous stabilizing law which satisfies the optimal criterion has not yet been solved.

In Section 2, the problem statement is given as well as a relaxed version that takes into account all the convex combinations of the subsystems. We explain why this relaxation is useful to solve the problem. Then, the necessary condition of the Pontryagin Maximum Principle are recalled. We also discuss the numerical difficulties encountered when singular controls enter in the solution. To circumvent this, a numerical framework is proposed to solve properly the optimization problem.

In Section 3, assuming that there exists at least a globally asymptotically stable (GAS) convex combination of the subsystems and as the positive definite solution of an algebraic Riccati equation is a continuous function with respect to the convex combination, we are able to build a parametrized family of positive definite function whose parameters belong to a a compact set. A control Lyapunov function is then defined as the point-wise infimum of this familly.

We prove that the Lyapunov function is locally Lipschizt and homogenous of degree two. Then, we show that its directional Dini derivative is well defined along trajectories and we deduce a state feedback leading to a cost value less than the value of the Lyapunov function. In addition, it is proved that the feedback makes the system globally exponentially stable. In section 4, we show that the sampled time version of the state feedback law is also globally exponentially stable. Finally, in Section 5, numerical examples are given which show actually that the optimal cost is finely approached.

Problem statement and necessary conditions

We consider the class of continuous time linear switched systems:

ẋ(t) = A σ(t) x(t) + B σ(t) u σ(t) (t) x(0) = x 0 (1)
where σ : [0, +∞) → S = {1, • • • , s} denotes the switching law that selects the active mode at time t by choosing among a finite collection of linear systems defined by the pairs

(A i , B i ) ∈ R n×n × R n×m i , i ∈ S.
Each subsystem is also govern by a control

u i (t) ∈ R m i , 0 ≤ m i ≤ n.
Our aim is to design a state feedback switching law (i.e.

x → (σ(x), u σ(x) (x))) for system (1) that approaches the optimal solution of the following optimization problem:

Problem 1: Minimize the switched quadratic criterion:

min σ(•),u σ(•) (•) 1 2 ∞ 0 x T (t)Q σ(t) x(t) + u T σ(t) (t)R σ(t) u σ(t) (t)dt (2)
where

Q i = Q T i ≥ 0, R i = R T i > 0, i ∈ S subject to ẋ(t) = A σ(t) x(t) + B σ(t) u σ(t) (t), x(0) = x 0 .
A usual framework [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF][START_REF] Sorin | Optimal control of switching systems[END_REF] to solve optimal control problem for switched systems ( ẋ = f i (x), i ∈ S) is to solve its relaxed version, replacing the vector field set (f i (x)) by its convex hull ( ẋ = co{f i (x)}). At least, three reasons justify the convexification of the problem: (i) the solutions are well defined (Fillipov; [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides: Control Systems[END_REF]); (ii) the density of the switched system trajectories into the trajectories of its relaxed version [START_REF] Ingalls | An infinite-time relaxation theorem for differential inclusions[END_REF]; (iii) the existence of singular optimal solutions are taking into account [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF][START_REF] Sorin | Optimal control of switching systems[END_REF].

The relaxed version of Problem 1 is then given as a pure continuous time optimal control problem given by :

ẋ(t) = s i=1 λ i (t)(A i x(t) + B i u i (t)) x(0) = x 0 (3) λ(t) ∈ Λ = λ ∈ R s : s i=1 λ i = 1 λ i ≥ 0 . ( 4
)
and by the convexified cost:

min λ(•), u i (•) 1 2 ∞ 0 s i=1 λ i (t)(x T (t)Q i x(t) + u T i (t)R i u i (t))dt (5)
In the sequel we denote by u = (u 1 , u 2 , • • • , u s ) and we use the following notation:

A(λ) = i∈S λ i A i , B(λ) = [λ 1 B 1 , λ 2 B 2 , • • • , λ s B s ] Q(λ) = i∈S λ i Q i , R(λ) = diag([λ 1 R 1 , λ 2 R 2 , • • • , λ s R s ])
Then the dynamics of the relaxed system can be redefined as:

ẋ = A(λ)x + B(λ)u
and the cost as

min λ(•),u(•) 1 2 ∞ 0 x T (t)Q(λ(t))x(t) + u T (t)R(λ)u(t)dt
To apply Pontryagin Maximum Principle (PMP) for Problem 1 or its relaxed version, the Hamiltonian function is defined as follow:

H(x, λ, u, p) = s i=1 λ i H i (x, u i , p) (6)
with

H i (x, u i , p) = p T (A i x + B i u i ) + 1 2 (x T Q i x + u T i R i u i )
and where p defines the costate. This leads to the following classical necessary conditions for optimality [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]: Theorem 1. Suppose that (λ * , u * ) is optimal with the corresponding state x * . Then, there exists an absolutely continuous function p * , called co-state, such that:

1. p * = 0, 2. ṗ * = s i=1 λ * i (t)(-A T i p * -Q i x * ) for almost all t ∈ R + , 3. (λ * (t), u * (t)) ∈ arg min (λ∈Λ,u) H(x * (t), λ, u, p * (t)), 4. H(x * (t), λ * (t), u * , p * (t)) = 0.
Theorem 1 can be simplified thanks to the following lemma:

Lemma 2. The optimal value of the u i 's are given by u * i (t) = -R -1 i B T i p * (t) and λ * satisfies:

λ * (t) ∈ arg min λ∈Λ s i=1 λ i H i (x * , -R -1 i B T i p * , p * ). ( 7
)
Proof. From Equation 6, the minimum of H with respect to the u i 's is clearly independent of the value of λ and the result follows.

From [START_REF] Cortes | Discontinuous dynamical systems[END_REF], it is clear that if there exists i ∈ S at time t such that

H i (x * (t), u * i , p * (t)) < H j (x * (t), u * j , p * (t)), ∀j ∈ S \ {i},
then the optimal control has to satisfy λ * i (t) = 1 and λ * j (t) = 0, ∀j ∈ S \ {i}.

A switching instant can occur at time t if there exists at least a pair (i, j) ∈ S 2 such that H i = H j = 0. At this time, the value of λ cannot be determined directly. Actually, if we suppose that 0 = H i = H j < H k , ∀k ∈ S \ {i, j} then the values that satisfy the relation, λ i + λ j = 1, are potential candidate for optimality. Moreover, a so called singular control λ can exist, for which 0 = H i = H j < H k , ∀k ∈ S \ {i, j} on a non empty time interval (a, b). This is a well known situation in the literature [START_REF] Robbins | A generalized Legendre-Clebsch condition for the singular cases of optimal control[END_REF][START_REF] Bryson | Applied Optimal Control: Optimization, Estimation, and Control[END_REF][START_REF] Bonnard | Singular Trajectories and Their Role in Control Theory[END_REF] and second order necessary conditions given by the Generalized Legendre-Clebsh Condition [START_REF] Volker | Singular optimal control -the state of the art[END_REF][START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF] can be necessary to solve the optimal control problem.

Definition 3. We call singular control, a control λ(.) such that there exist at least two indices i, j, for which

H i = H j = 0 on a non zero measure time interval (a, b) and satisfying λ(t) ∈ Λ, λ k (t) = 1 ∀k ∈ S, ∀t ∈ (a, b).
The corresponding part of the trajectory is called a singular arc.

A singular control defines a Fillipov solution [START_REF] Cortes | Discontinuous dynamical systems[END_REF] for the original switched system (1). Hence, it allows to extend properly the notion of optimal solution for switched systems.

Roughly speaking when an optimal solution of the relaxed problem possesses singular arcs, these arcs define sliding surfaces for the switched system (1) which lead to chattering if the surface is attractive. It is noteworthy that only suboptimal solutions can be achieved for the switched systems due to the limited switching frequency; see for example [START_REF] Sorin | Optimal control of switching systems[END_REF].

Numerical resolution

If one attempts to solve numerically an optimal control problem in which singular arcs appear, numerical difficulties will be encountered due to the insensitive of the Hamiltonian with respect to the control. On the one hand, standard indirect numerical methods such as multiple shooting methods are not appropriate to deal with singular arcs without a priori information on the structure of the trajectories (see [START_REF] Fraser-Andrews | Finding candidate singular optimal controls: A state of the art survey[END_REF], [START_REF] Maurer | Numerical solution of singular control problems using multiple shooting techniques[END_REF] and [START_REF] Bonnans | Singular arcs in the generalized Goddard's problem[END_REF] to apply multiple shooting methods in this context). This information can be achieved using regularization techniques such as the continuation method used in [START_REF] Moylan | Generalizations of singular optimal control theory[END_REF], [START_REF] David | Computation of optimal singular controls. Automatic Control[END_REF], [START_REF] Martinon | Numerical resolution of optimal control problems by a simplicial method[END_REF]. On the other hand, when direct methods such as nonlinear programming (NLP) are used to solve (1)-( 2), a very bad results are generally obtained on the optimal value for the control due to ∂H ∂u = 0, ∀u in an open subset of U . In [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF], we have proposed to use a mix direct-indirect method.

The idea consists to take implicitly into account the singular arcs using the necessary condition of the PMP and the Hamiltonian systems and then to solve directly an augmented constraint optimization problem.

For the LQ switched problem, this leads to the following constraint optimization problem (here z = (x, p)):

Problem 2: Minimize (using NLP):

min λ(•) 1 2 ∞ 0 s i=1 λ i (t)(x T (t)Q i x(t) + p T (t)B i R -1 i B T i p(t))dt (8) subject to ż(t) = s i=1 λ i (t) A i -B i R -1 i B T i -Q i -A T i z(t) (9) 0 ≤ λ i ⊥ H i (x, -R -1 i B T i p, p) ≥ 0, i ∈ S (10) λ(t) ∈ Λ, x(0) = x 0 (11)
where the sign x ⊥ y means xy = 0.

The complementarity constraints (10) allow λ to be multivalued function when the H i 's vanish for at least two subscripts. This is a key point since the necessary condition of PMP does not imply the unicity of solution of system (9) from an initial condition.

When the admissible values for λ are multiple, this formulation by the minimization of the cost (8) yields the optimal value.

Practically, the complementarity constraint are taken into account with penalization terms ρλ i H i in the cost with weight ρ. Hight order constraints can be also added to improve the result see [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF] for more details.

Lyapunov based switching law

The aim of this part is to define a Lyapunov function as a tight upper bound on the value function. We means tight in the sense that the two functions may coincide at some points.

We denote by B(

√ λ) the matrix B( √ λ) = √ λ 1 B 1 | √ λ 2 B 2 | . . . | √ λ s B s
and we assume there exists a λ 0 ∈ Λ such that the pair (A(λ 0 ), B( √ λ 0 )) is controllable. Then, the Riccati equation:

A(λ) T P λ + P λ A(λ) -P λ B( √ λ)R -1 B( √ λ) T P λ + Q(λ) = 0. (12) R = diag([R 1 , R 2 , • • • , R s ]
), admits a (unique) symmetric positive definite solution P λ 0 for λ = λ 0 . Notice that the pack-writing term

P λ B( √ λ)R -1 B( √ λ) T P λ can be expanded as P λ B( √ λ)R -1 B( √ λ) T P λ = i∈S λ i (P λ B i R -1 i B T i P λ ).
In fact to ensure the existence of a positive definite solution to the Riccati equation ( 12), as the Q i 's are positive definite, it is enough to assume that the pair (A(λ 0 ), B( √ λ 0 )) is stabilizable (cf. Lemma 4). Now, if the pair (A(λ 0 ), B( √ λ 0 )) is controllable, the same is true if λ belongs to a sufficiently small neighborhood of λ 0 ; so the Riccati equation ( 12) admits a unique (positive definite) solution P λ for every λ in some neighborhood of λ 0 . Moreover, it is well known that the positive definite solution of a Riccati equation is a continuous function of (A(λ), B( √ λ), R, Q(λ)) [START_REF] Lancaster | Algebraic Riccati Equations[END_REF] and so P λ is a continuous function of λ. Notice also that the condition for the existence of a symmetric nonnegative solution of the Riccati equation ( 12) can be weakened: in [START_REF] Kučera | A review of the matrix Riccati equation[END_REF], V. Kučera proved that if the pair (A(λ), B( √ λ)) is stabilizable and if the matrix

M =   A(λ) -B( √ λ)R -1 B( √ λ) T -Q(λ) -A(λ) T  
has no purely imaginary eigenvalues, then there exists a symmetric nonnegative solution to equation [START_REF] José | Stability and stabilization of discrete time switched systems[END_REF]. This result allow us to prove the following lemma about Eq. ( 12).

Lemma 4. If the pair (

A(λ), B( √ λ)
) is stabilizable and Q(λ) is positive definite, then there exists a positive definite solution to Eq. [START_REF] José | Stability and stabilization of discrete time switched systems[END_REF].

Proof. We denote by S(λ) the matrix

B( √ λ)R -1 B( √ λ)
T and we take x = (x T 1 , x T 2 ) T a vector of C 2n such that M x = iαx (with α ∈ R); we shall see that x = 0. We have

A(λ)x 1 -S(λ)x 2 = iαx 1 (13) -Q(λ)x 1 -A(λ) T x 2 = iαx 2 . ( 14
)
Multiplying on the left the members of equation ( 13) (resp. Eq. ( 14)) by xT 2 (resp. by xT 1 ) (the bar denotes the conjugate), we get

xT 2 A(λ)x 1 -xT 2 S(λ)x 2 = iαx 1 (15) -x T 1 Q(λ)x 1 -xT 1 A(λ) T x 2 = iαx 2 . ( 16
)
by adding the conjugate of the members of equation ( 15) to the members of equation ( 16), we get

-x T 1 Q(λ)x 1 -xT 2 S(λ)x 2 = 0. If x 1 = 0, as Q(λ) is positive definite,

this last equality implies xT

2 S(λ)x 2 < 0, but this inequality cannot occur because S(λ) is nonnegative, so we must have x 1 = 0 and reporting this equality in ( 13) and ( 14), we get S(λ)x 2 = 0 and A(λ

) T x 2 = -iα x 2 , which in turn implies A(λ) T x 2 = -iα x 2 and B( √ λ) T x 2 = 0. As the pair (A(λ) T , B( √ λ) T ) is
detectable, the Hautus lemma implies that x 2 = 0.

By applicantion of the above-mentionned result from Kučera, we deduce that there exists a symmetryc nonnegative solution P λ to Eq. [START_REF] José | Stability and stabilization of discrete time switched systems[END_REF]. Now this solution is necessarily definite, assume indeed that v is a vector such that P v = 0, left-multiply both sides of ( 12) by vT and right-multiply by v, we get vT Qv = 0 which implies v = 0 since Q is assumed to be positive definite.

This lemma proves that for every λ ∈ Λ such that the pair (A(λ), B( √ λ)) is stabilizable, there exists a positive definite solution, denoted by P λ , to the Ricatti equation [START_REF] José | Stability and stabilization of discrete time switched systems[END_REF].

We denote by Λ + the set

Λ + = λ ∈ Λ | the pair(A(λ), B( √ λ)
) is stabilizable and max spec(P λ ) ≤ ν max where spec(P λ ) denotes the spectrum of P λ ; this set Λ + satisfies the following property.

Lemma 5. The matrices Q i being positive definite, if there exists λ 0 ∈ Λ such that (A(λ 0 ), B( √ λ 0 )) is controllable, then, for every ν max large enough, set Λ + is compact and its interior is not empty in Λ. Moreover, there exist positive real numbers α m and α M defined as

α m = min λ∈Λ + min(spec(P λ )) > 0 α M = max λ∈Λ + max(spec(P λ )) > 0.
Proof. As noticed above, we can find a compact neighborhood U of λ 0 such that the Riccati equation ( 12) admits a positive definite solution P λ for every λ ∈ U . The mapping λ → P λ being continuous and U being compact, we have sup λ∈U max spec(P λ ) < ∞; this implies that if ν max is chosen large enough, the interior of Λ + is non empty. Now, set Λ + is included in Λ, therefore it is bounded; we shall show that it is also closed. Suppose that there exists a sequence (λ k ) k≥1 ∈ Λ + such that lim k→∞ λ k = λ. As Λ is a compact set, λ ∈ Λ. Moreover, the sequence (P λ k ) k≥1 is bounded, so we can assume that it converges to a symmetric matrix P . As a limit of a sequence positive definite matrices, this matrix is positive (semi) definite; moreover it is a solution of [START_REF] José | Stability and stabilization of discrete time switched systems[END_REF] with λ = λ. We claim first that P is definite, assume indeed that v is a vector such that P v = 0, left-multiply both sides of ( 12) by vT and right-multiply by v, we get vT Qv = 0 which implies v = 0 since Q is assumed to be positive definite. Moreover the pair (A(λ), B( √ λ) is stabilizable, to see this let µ be an eigenvalue of A(λ) T such that ℜ(µ) ≥ 0 (ℜ(•) stands for the real part) and let v ∈ R n be a vector such that A(λ) T v = µ v and B T ( √ λ)v = 0, we shall see that v = 0 which implies that the rank of the matrix (A(λ

) T -µ Id, B T ( √ λ))
T is equal to n for every µ in the closed right half plane and so the result will follow from the Hautus lemma. Matrix P being definite, there exists x ∈ R n such that P x = v, left-multiply both sides of ( 12) by xT and right-multiply by x, we get

µ xT v + μ vT x + xT Q x = 0.
If x = 0, as matrix Q is definite positive, this equality implies that ℜ(µ xT v) < 0 but xT v = xT P x so xT v > 0 since P is positive definite, therefore we must have ℜ(µ) < 0. This contradicts the fact that ℜ(µ) ≥ 0, so, we must have x = 0, which implies that v = 0. Finally, the existence of α m and α M follows from the compactness of Λ + and the continuity of the mapping λ → P λ .

For the sake of readability, let us introduce the following notations. We denote by:

M (λ) := i∈S λ i M i (λ) (17)
where

M i (λ) := A i -B i K i (λ) (18)
K i (λ) := R -1 i B T i P λ ( 19
)
and by

N (λ) := i∈S λ i N i (λ) (20)
where

N i (λ) := Q i + K i (λ) T R i K i (λ).
The Riccati equation ( 12) can then be rewritten as :

M (λ) T P λ + P λ M (λ) + N (λ) = 0. ( 21
)
Lemma 6. For every (x, λ) ∈ R n {0} × Λ + , we have

min i∈S 2x T M T i (λ)P λ x + x T N i (λ)x ≤ 0.
Proof. Take x ∈ R n {0} and λ ∈ Λ + , then equation ( 21) admits a solution and we can write

i∈S λ i 2x T M T i (λ)P λ x + x T N i (λ)x = 0 (22) so we cannot have 2x T M T i (λ)P λ x + x T N i (λ)x > 0 (23)
for every i ∈ S because in this case the left-hand member of equality ( 22) would be positive. Thus, for every pair (x, λ), there exists an index i which is such that the left-hand member in ( 23) is non positive.

Let us now introduce the following Lyapunov function

V m (x) := inf λ∈Λ + x T P λ x ( 24
)
where P λ denotes the solution of equation ( 21). Clearly, as every P λ is positive definite when λ belongs to Λ + and as the set {P λ | λ ∈ Λ + } is compact, V m is a positive definite function; notice also that V m is homogeneous of degree 2 and locally Lipschitz.

Proposition 7. The function defined by ( 24) is locally lipschitzian.

Proof. We have z T P λ z ≤ α M z 2 for every (z, λ) ∈ R n × Λ + . So, if we take x and y in the ball B(0, R), we obtain easily that, for every λ ∈ Λ + ,

|x T P λ x -y T P λ y| = (x -y) T P λ (x + y)| ≤ K x -y
where K = 2α M R. Thus, the family of functions indexed by λ ∈ Λ + and defined by x → x T P λ x is uniformly locally Liptchitz with a Lipschitz constant equals to K := 2M R on the ball B(0, R).

Now, as the function λ → x T P λ x is continuous and set Λ + is compact, there exists a pair (λ 1 , λ 2 ) ∈ (Λ + ) 2 such that:

V m (x) = x T P λ 1 x V m (y) = y T P λ 2 y.
From the definition of V m , we deduce easily that

x T P λ 1 x -x T P λ 2 x ≤ 0 y T P λ 1 y -y T P λ 2 y ≥ 0
therefore, by continuity, there exists z on the line segment [x y] such that z T P λ 1 z = z T P λ 2 z and it follows:

|V m (x) -V m (y)| ≤ |x T P λ 1 x -z T P λ 1 z| + |z T P λ 2 z -y T P λ 2 y| (25) ≤ K x -z + K z -y = K x -y . ( 26
)
Finally function V min is proper because it is continuous and

α m z 2 ≤ V m (z) for every z ∈ R n .
Let f be a function defined on R n and d be a vector of R n , as in [START_REF] Furukawa | Optimality conditions in nondifferentiable programming and their applications to best approximations[END_REF] we shall denote by f ′ (x; d) the following limit (if it exists)

f ′ (x; d) := lim h→0 h>0 f (x + h d) -f (x) h . ( 27
)
For the sake of readability, we let v λ (x) := x T P λ x. In order to compute V ′ m (x; d), we use Theorem 6.1 in [START_REF] Furukawa | Optimality conditions in nondifferentiable programming and their applications to best approximations[END_REF] whose conditions of application are clearly met, thus we have

V ′ m (x; d) = inf λ∈L(x) v ′ λ (x; d). ( 28
)
Here L(x) denotes the set of λ ∈ Λ + such that V m (x) = v λ (x) (this set is clearly nonempty and compact because Λ + is compact and the function λ → v λ (x) is continuous). As the function x → v λ (x) is smooth, from ( 28), we infer that

V ′ m (x; d) = 2 inf λ∈L(x) d T P λ x.
From this formula and from Lemma 6, we get the following properties.

Lemma 8. For every (x, λ 0 ) ∈ R n × L(x), there exist i(x, λ 0 ) such that

V ′ m (x; M i(x,λ 0 ) (λ 0 )x) ≤ -x T N i(x,λ 0 ) (λ 0 )x ( 29
)
Proof. For every (x, λ 0 ) ∈ R n × L(x), Lemma 6 implies that there exists i(x, λ 0 ) such that 2x T M T i(x,λ 0 ) (λ 0 )P λ 0 x + x T N i(x,λ 0 ) (λ 0 )x ≤ 0 It follows directly that

V ′ m (x; M i(x,λ 0 ) (λ 0 )x) = inf λ∈L(x) 2x T M T i(x,λ 0 ) (λ 0 )P λ x ≤ 2x T M T i(x,λ 0 ) (λ 0 )P λ 0 x ≤ -x T N i(x,λ 0 ) (λ 0 )x
In the following theorem, we shall consider mappings from R n to S × Λ + of the form x → (i(x), λ(x)) such that λ(x) ∈ L(x). To such a mapping, we relate the following feedback law for system (1): the mode σ(t) is equal to i(x(t)) for every t ≥ 0 and u σ(t) is equal to -K i(x) (λ(x)) x. Theorem 9. We assume that the matrices Q i are positive definite and there exists at least a λ ∈ Λ such that the pair (A(λ), B( √ λ)) is controllable. For every x ∈ R n , we choose (i(x), λ(x)) ∈ arg min

(i,λ)∈S×L(x) 2x T M i (λ)P λ x + x T N i (λ)x .
Then the feedback related to (i(x), λ(x)) stabilizes the switched system (1) with a cost smaller than 1 2 V min (x 0 ). Moreover the convergence is exponential with a rate β = η 0 α 1 where η 0 and α 1 are given by:

η 0 = min i∈S inf x∈S n-1 inf λ∈L(x) x T N i (λ)x, α 1 = max x∈S n-1 V m (x)
Remark. From a practical point of view, Th. 9 remains still valid if the feedback switching law is simplified as follows: for a given x, choose λ(x) in L(x) and take i(x) as i(x) ∈ arg min

(i)∈S 2x T M i (λ)P λ x + x T N i (λ)x .
Proof of Th. 9. We shall compute the derivative of V m along the trajectories of system (1) in closed-loop with the feedback introduced in the theorem. Hereafter, for the sake of readability, we denote by

M i(x) the matrix M i(x) (λ(x)) and d dt V m (x(t)) denotes the Dini derivative (cf. ( 27
)) of V m . d dt V m (x(t)) = V ′ m (x; M i(x) x) = 2 min λ∈L(x) x T M T i(x) P λ x ≤ 2x T M T i(x) P λ(x) x ≤ -x T N i(x) (λ(x))x from the definition of (i(x), λ(x)) ≤ - η 0 α 1 V m (x).
This inequality implies that

V m (x(t)) ≤ e -βt V m (x 0 )
with β = η 0 α 1 . As V m is homogeneous of degree 2, this last inequality, implies the global exponential stability.

The upper bound on the cost (1/2V m (x 0 )) comes from the fact that

x T Q i(x) x + x T K T i(x) R i(x) K i(x) x = x T N i(x) x ≤ d dt V m (x(t)).
Why do we claim that the Lyapunov function can be a tight upper bound on the value function ? Observe first that in the case where all subsystems (related to the pairs (A i , B i )) are stabilizable, then the solution P i of the Algebraic Riccati Equation exists for each mode i and the Lyapunov function satisfies always the following inequality:

1 2 V m (x) ≤ min i∈S 1 2 x T P i x
One can also observe that for a given state x 0 , the value 1 2 V m (x 0 ) is the best cost related to every constant convex combination that stabilizes the relaxed system. The corresponding control is of the form :

λ(t) = λ(x 0 ), ∀t ≥ 0, u i (t) = -R -1 i B T i P λ x(t).
In the general case, when can we say that 1 2 V m (x 0 ) is optimal ? The answer is :'Along the part of trajectories where the optimal control λ * is constant to reach the origin". At least two cases can be mentioned: if the number of switchings is finite which means that a same mode is used after a time t or if the trajectory is steered to the origin by a constant singular control λ (cf. Definition 3) for which P λ > 0. Note that singular control in dimension n = 2 can be algebraically determined [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF] and are constant.

Sampled time switched control law

Practically, a sampled time version of the continuous time algorithm is applied. So, we shall show that the sampled time version of the above algorithm stabilizes the system for appropriate choice of sampled period.

Let τ be a given sampling period . The control is now piecewise constant and updated every times t k = kτ , k ∈ N following the state feedback provided by Theorem 9. To be more precise, at time t 0 = 0, we start with the initial condition x 0 and we choose the mode i(x 0 ) and λ(x 0 ) as in Theorem 9. Thus we have

2x T 0 M T i 0 (x 0 ) P λ(x 0 ) x 0 ≤ -x T 0 N i(x 0 ) (λ(x 0 ))x 0 ≤ -η 0 x 0 2 ,
and, by the way, notice that v λ(x 0 ) (x 0 ) = V m (x 0 ). We apply the feedback law related to i(x 0 ), λ(x 0 ) to system [START_REF] Axelsson | Gradient descent approach to optimal mode scheduling in hybrid dynamical systems[END_REF], that is to say, we choose the mode i(x 0 ) (σ(t) = i(x 0 )) and u σ(t) = -K i(x 0 ) (λ(x 0 )), we do so during a time τ > 0. At time, t 1 = t 0 + τ , we arrive at a point x 1 and we choose a new index i(x 1 ) and a new λ(x 1 ) as in theorem 9. In this way, we build a sequence of points (x k ) k≥0 together with a sequence of pairs (i(x k ), λ(x k )) k≥0 chosen as in theorem 9. On each interval [t k , t k+1 ) (t k = t 0 + k τ ), we choose the mode i(x k ) and we apply the feedback u = -K i(x k ) (λ(x k )). The switched system (1) in closed loop with this feedback writes

   ẋ = A i(x k ) -B i(x k ) K i(x k ) (λ(x k )) x, t ∈ [t k , t k+1 ) x k = x(t k ) (30)
We state the following theorem about this algorithm; notice that we still assume that the Q i 's are positive definite.

Theorem 10. If τ is chosen sufficiently small, the sampled time switching law described above stabilizes the switched system (30) globally exponentially

For the proof of the theorem, we shall need the following lemma.

Lemma 11. Consider the solution t → x(t) of (30) starting from x 0 at some time t 0 = 0 where the mode i is chosen such that

x T 0 P λ(x 0 ) x 0 = V m (x 0 ) 2x T 0 M T i (λ(x 0 ))P λ(x 0 ) x 0 ≤ -x T 0 N i (λ(x 0 ))x 0 ≤ -η 0 x 0 2
where η 0 is defined in Theorem 9. We define the time T as

T = inf t ≥ 0 | 2x(t) T M T i (λ(x 0 ))P λ(x 0 ) x(t) ≥ -η 0 γ x(t) 2
where γ is a parameter chosen in the interval (0, 1). Then there exists τ 0 > 0 independent from x 0 such that T ≥ τ 0 .

Proof. For the sake of readability, we shall denote by i and λ the terms i(x 0 ) and λ(x 0 ) respectively. We have

2x(T ) T M T i (λ)P λ x(T ) -2x T 0 M T i (λ)P λ x 0 = 2 T 0 d dt x(t) T M T i (λ)P λ x(t) dt = 2 T 0 x(t) T (M T i (λ)) 2 P λ x(t) + x(t) T M i (λ)P λ M i (λ)x(t) dt. (31)
Now, for every t ∈ [0, T ), we have x(t) = e t M i (λ) x 0 and so x(t) ≤ e t M i (λ) x 0 .

Notice that the matrix M i (λ(x 0 )) is bounded because λ(x 0 ) evolves on the compact set Λ + and the mapping λ → M i k (λ) is continuous. So there exists µ > 0 such that M i (λ) ≤ µ for every (i, λ) ∈ S × Λ + ; also, from Lemma 5, we know that P λ ≤ α M for every λ ∈ Λ + . So, from equality (31), we deduce that

2x(T ) T M T i (λ)P λ x(T ) -2x T 0 M T i (λ)P λ x 0 ≤ 4µ 2 α M x 0 2 T 0 e 2t µ dt = 2µα M (e 2(T µ -1) x 0 2 . (32) Now, since 2x T 0 M T i (λ)P λ x 0 ≤ -η 0 x 0 2 and 2x(T ) T M T i (λ)P λ x(T ) = -η 0 γ x(T ) 2 ≥ -η 0 γ e 2T µ x 0 2 , we have 2x(T ) T M T i (λ)P λ x(T ) -2x T 0 M T i (λ)P λ x 0 ≥ η 0 x 0 2 (1 -γ e 2T µ ).
Substituting this inequality in (32), we get 2µα M (e 2(T µ -1)

x 0 2 ≥ η 0 x 0 2 (1 -γ e 2T µ )
and so (2µ

α M + η 0 γ)e 2T µ ) ≥ 2µ α M + η 0 which implies that T ≥ τ 0 := 1 2µ ln 2µ α M + η 0 2µ α M + η 0 γ .
Proof of Theorem 10. We take a sampling time τ > 0 no greater than τ 0 , and we shall prove that system (30) is globally exponentially stable about the origin. Recall that, due to the choice of the pair (i k , λ k ), we have v λ k (x k ) = V m (x k ) for every index k. Recall also that from Lemma 5, we have the inequality x T P λ x ≥ α m x 2 for every λ ∈ Λ + and every x ∈ R n . The expression 2x(t) T M T i k (λ k )P λ k x represents the derivative of v λ k along the trajectories of system (30) on the time interval [t k , t k+1 ). So due to the choice of our feedback, this derivative is less or equal to -η 0 γ x(t) 2 for every t ∈ [t k , t k+1 ), so we have

vλ k (x) ≤ -η 0 γ x 2 ≤ - η 0 γ λ M (P λ k ) v λ k (x)
where λ M (P λ k ) denotes the greatest eigenvalue of matrix P λ k . From Lemma 5, we have λ M (P λ k ) ≤ α M and so we obtain

vλ k (x) ≤ - η 0 γ α M v λ k (x) which implies that v λ k (x k+1 ) ≤ v λ k (x k )e - η 0 γ α M τ ( 33
)
for every index k. Now, we have

V m (x k+1 ) ≤ v λ k (x k+1 ) and V m (x k ) = v λ k (x k ), so from (33), we deduce V m (x k+1 ) ≤ θ V m (x k ) (34)
where θ = e -η 0 γ α M τ ∈ (0, 1) and so

V m (x k ) ≤ θ k V m (x 0 ) which implies , as V m (x) ≥ α m x 2 , x k 2 ≤ 1 α m θ k x 0 2 . ( 35
)
Now, on the interval [t k , t k+1 ), we have v λ k (x(t)) ≤ v λ k (x k ), this inequality together with (35), proves the result.

Illustratives examples

Before presenting some examples, it is important to mention that it is not necessary to ensure a stabilizing switched law to determine all the possible values of the set Λ + . Only one value is sufficient to guarantee the stability. So, a reasonable finite number of values ensures performances. Practically, a finite number have been used using a discretization of the set Λ.

Example 1 : a regular case

Consider a two mode switched system with the following design parameters:

A 1 = -2.7 3.9 4.4 -12.6 , A 2 = -9.5 -5.1 -7.5 -3.3 , B 1 = 0.1 0 , B 2 = 4.6 0 , Q 1 = Q 2 = Id, R 1 = 1 and R 2 = 2.
For each subsystem, an LQ design can be be performed separately. Thus, in order to make some comparisons, Figure 1 shows the state space trajectories for the switching law given by ( 30) and the optimal one. The later is obtained by NL programming in a suitable formulation taking into account singular arcs [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF]. If not, numerical difficulties in the control determination are often encountered. We can see that the two solutions match well together. In the example, it can be observe that two singular arcs (defined by two lines in the state space) occur in the solutions.

Figure 2 compares the optimal cost with the costs obtained by using the switching law, only mode 1 and only mode 2, respectively. This comparison is made for initial states taken on the unit ball, the x-axis represent the angle θ. We have also added the guarantee on the cost provided by upper bound i.e. 1 2 V min (x). It can be observed that the cost associated to the switching law coincides with the cost of the optimal numerical solution. Of course, the essential difference is that the numerical solution is an open loop control while the switching law defines a closed loop control. It is also clear that the used of a single mode with no switching leads to lower performances.

Example 2 : non stabilizable subsystems

For this second example, we have chosen two non stabilizable subsystems:

A 1 = 1 0 1 -1 , A 2 = 2 1 0 1 , B 1 = 0 1 , B 2 = 1 0 .
Therefore, there is no LQ design that can be defined separately for each subsystem. However, in a switched framework and taking Q 1 = Q 2 = Id, R 1 = 1 and R 2 = 2, as the set Λ + is non empty, the switching law presented in this paper can be applied. Once again the optimal solution and the one provided by the switching law are very closed as showed by Figure 3. Figure 4 compares the optimal cost with the costs obtained by using the switching law. The upper bound V min is also plotted.

Figure 1 :Figure 2 :

 12 Figure 1: Ex. 1: State space trajectories: (red) optimal solution (NLP); (blue) switching law

Figure 3 :Figure 4 :

 34 Figure 3: Ex. 2: State space trajectories: (red) optimal solution (NLP); (blue) switching law


Yang Chen

Xindao Yin

Luyao Shi

Huazhong Shu

Limin Luo

Jean-Louis Coatrieux

Christine Toumoulin

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 1

Keywords: Low-dose CT (LDCT), abdomen tumor, dictionary learning, noise, artifacts

In abdomen computed tomography (CT), repeated radiation exposures are often inevitable for cancer patients who receive guided surgery or radiotherapy. Low-dose scans should thus be considered in order to avoid too high accumulative harm of radiation. This work is aimed at improving abdomen tumor CT images from low-dose scans by using a fast dictionary learning (DL) based processing. Stemming from sparse representation theory, the patch-based DL approach proposed in this paper allows effective suppression of both mottled noise and streak artifacts. The experiments carried out on clinical data show that the proposed method brings encouraging improvements in abdomen low-dose CT images with tumors.

1.Introduction

Since the introduction of computed tomography (CT) in the 1970s, a wide-spread concern on CT is the increasing risk of cancer induced by X-ray radiation [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF][START_REF] Brenner | Computed tomography-an increasing source of radiation exposure[END_REF]. The radiation doses delivered to patients during X-ray CT examinations are relatively high when compared to other radiological examinations [START_REF]Patient safety: Radiation Exposure in X-ray and CT Examinations[END_REF]. Based on a recent report, the overall averaged radiation dose associated with a routine abdomen/pelvis CT scan is around 10 mSv, which is roughly 5 times of head CTs and 100 times of chest X-ray Radiography [START_REF]Patient safety: Radiation Exposure in X-ray and CT Examinations[END_REF]. Additionally, the dose in CT is cumulative in lifetime, and successive CT scanning can significantly increase the lifetime radiation risk of fatal cancers [START_REF]Patient safety: Radiation Exposure in X-ray and CT Examinations[END_REF][START_REF] Smith-Bindman | Radiation dose associated with common computed tomography examinations and the associated lifetime attributable risk of cancer[END_REF].

CT is also frequently used to guide the surgery or radiotherapy by providing localized contrast information between tumors, organs and other surrounding human tissues. Patients with diagnosed or suspicious abdomen tumors would be submitted to repeated CT scans over a long observing period before or after surgery or therapy. Low dose CT (LDCT) is therefore of major importance in order to alleviate the harm caused by cumulated radiations for the patients with abdomen tumors. Among all the methods proposed so far to obtain LDCT images, the most practical and widely used method is lowering the X-ray tube current by modulating the mA or mAs setting, but at a cost of degraded CT image quality due to increased quantum noise and artifacts [START_REF] Mannudeep | Strategies for CT Radiation Dose Optimization[END_REF][START_REF] Yazdi | Artifacts in Spiral X-ray CT Scanners: Problems and Solutions[END_REF][START_REF] Watanabe | Improvement of image quality of low radiation dose abdominal CT by increasing contrast enhancement[END_REF]. In the past ten years, other approaches have been explored to improve the quality of LDCT images. They can be divided into three categories: pre-processing approaches, reconstruction approaches and post-processing approaches.

The first one refers to techniques that improve the CT imaging by restoring the projected raw data before filtered-backprojection (FBP) reconstruction. Adaptive filtering, multiscale penalized weighted least-squares and bilateral filtering have been reported to suppress the excessive quantum noise in projected raw data [START_REF] Hsieh | Adaptive streak artifact reduction in computed tomography resulting from excessive x-ray photon noise[END_REF][START_REF] Wang | Multiscale penalized weighted least-squares sinogram restoration for low-dose X-Ray computed tomography[END_REF][START_REF] Ehman | Noise reduction to decrease radiation dose and improve conspicuity of hepatic lesions at contrast-enhanced 80-kV hepatic CT using projection space denoising[END_REF]. Reconstruction approaches treat the LDCT imaging as an ill-posed inverse problem, and solve the problem via maximizing a prior-regularized cost function using iterative optimizations [START_REF] Kamphuis | Accelerated iterative transmission CT reconstruction using an ordered subset convex algorithm[END_REF][START_REF] Nuyts | Iterative reconstruction for helical CT: A simulation study[END_REF][START_REF] Thibault | A three-dimensional statistical approach to improved image quality for multislice helical CT[END_REF][START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF][START_REF] Chen | Nonlocal Prior Bayesian Tomographic Reconstruction[END_REF][START_REF] Chen | Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting nonlocal prior[END_REF][START_REF] Lubner | Reduced image noise at low-dose multidetector CT of the abdomen with prior image constrained compressed sensing algorithm[END_REF][START_REF] Prakash | Reducing abdominal CT radiation dose with adaptive statistical iterative reconstruction technique[END_REF][START_REF] Silva | Innovations in CT dose reduction strategy: application of the adaptive statistical iterative reconstruction algorithm[END_REF]. Many prior options have been proposed in the past decade, for example the total-variation minimization in [START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF], the nonlocal prior reconstruction in [START_REF] Chen | Nonlocal Prior Bayesian Tomographic Reconstruction[END_REF][START_REF] Chen | Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting nonlocal prior[END_REF], and the prior image constrained compressed sensing (PICCS) algorithm in [START_REF] Lubner | Reduced image noise at low-dose multidetector CT of the abdomen with prior image constrained compressed sensing algorithm[END_REF]. It

Improving Abdomen Tumor Low-dose CT Images Using A Fast Dictionary Learning Based Processing

should be noted that results of clinical value in abdomen LDCT have been achieved by using the PICCS algorithm [START_REF] Lubner | Reduced image noise at low-dose multidetector CT of the abdomen with prior image constrained compressed sensing algorithm[END_REF] and the adaptive statistical iterative reconstruction (ASIR) [START_REF] Prakash | Reducing abdominal CT radiation dose with adaptive statistical iterative reconstruction technique[END_REF][START_REF] Silva | Innovations in CT dose reduction strategy: application of the adaptive statistical iterative reconstruction algorithm[END_REF]. However, due to the difficult access to well-formatted projection data of the main CT vendors, researches on pre-processing and reconstruction approaches are often limited in practice. Another well-known concern for iterative reconstructions is the intensive computation cost required for reconstruction, which may delay clinical workflow and diagnosis. Additionally, the CT scanners equipped in most current hospitals are based on FBP algorithms and upgrading to the latest CT scanners with iterative algorithm is often too expensive.

The approach described in this paper falls into the third category, i.e post-processing methods. They can be directly applied on LDCT images to suppress noise and artifacts. The main merits of post-processing methods include their low computation cost and their easy implementation in most current CT systems with no built-in low dose solutions. The objective when applying post-processing on LDCT images is to obtain images with visual appearances close to the corresponding SDCT (standard dose CT) images without introducing structure lost and new artifacts (or false structures). However, due to the fact that the back-projection process in the FBP algorithm distributes the noise and artifacts non-uniformly over the whole image, CT noise and artifacts with mottled or streak effects are often difficult to be removed from already reconstructed images. The noise and artifacts in abdomen LDCT images cannot be modeled by Gaussian or Poisson distributions and restorations based on such priors lead to poor performance. This explains why post-processing techniques rely on more heuristic approaches. Furthermore, compared to other abdominal pathological features such as hepatic cyst, tumor tissues are often of lower contrast, and it is rather challenging to preserve their characteristics through image processing [START_REF] Watanabe | Improvement of image quality of low radiation dose abdominal CT by increasing contrast enhancement[END_REF]. In [START_REF] Kalra | Low-dose CT of the abdomen: evaluationof image improvement with use of noise reduction filters pilot study[END_REF][START_REF] Kalra | Detection and characterization of lesions on low-radiation-dose abdominal CT images postprocessed with noise reduction filters[END_REF], several noise reduction filters were proposed to enhance the conspicuity of lesions in abdomen LDCT images. An adaptive noise reduction filtering which combines smoothing and edge-enhancing was reported to gain 50% dose reduction (tube current reduced from 160mAs to 80mAs) without losing low-contrast detectability [START_REF] Funama | Improvement of low-contrast detectability in low-dose hepatic multidetector computed tomography using a novel adaptive filter: evaluation with a computer-simulated liver including tumors[END_REF]. In [START_REF] Chen | Improving Low-dose Abdominal CT Images by Weighted Intensity Averaging over Large-scale Neighborhoods[END_REF], a filtering technique, named large-scale nonlocal means (LNLM), was proposed to improve the quality of abdomen LDCT images with hepatic cysts. This LNLM method was further combined with a multiscale directional diffusion to suppress the streak artifacts in thoracic CT images [START_REF] Chen | Thoracic low-dose CT image processing using an artifact-suppressed largescale nonlocal means[END_REF].

Recent years have witnessed a growing interest in the study of sparse and redundant representations over dictionary learning (DL) [START_REF] Lewic | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF][START_REF] Lewicki | Learning overcomplete representations[END_REF][START_REF] Delgado | Dictionary learning algorithms for sparse representation[END_REF][START_REF] Donoho | Maximal sparsity representation via l1 minimization[END_REF][START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Mairal | Sparse representation for color image restoration[END_REF][START_REF] Mairal | Learning multiscale sparse representations for image and video restoration[END_REF][START_REF] Wright | Robust face recognition via sparse representation[END_REF]. Some successful applications have been explored in medical imaging [START_REF] Ravishankar | MR Image Reconstruction From Highly Undersampled k-space Data by Dictionary Learning[END_REF][START_REF] Xu | Low-Dose X-ray CT Reconstruction via Dictionary Learning[END_REF][START_REF] Li | An efficient dictionary learning algorithm and its ication to 3-D medical image denoising[END_REF]. Instead of being based on a pixel-wise intensity update, patch-wise DL processing inherently enables a more effective representation of patch-shaped features such as tumors or organs. In the paper, we propose to apply patch-based DL processing to improve abdomen tumor LDCT images. This method is referred to as DL algorithm. In the experiments presented below, clinical LDCT images of cancer patients were used. The corresponding SDCT images were also acquired to have a ground truth reference. The proposed DL algorithm can be efficiently implemented by using a global dictionary and a parallelization technique. The relevance of this DL approach is shown through experiments conducted on a large set of patients. The structure of this paper is as follows: in Section 2, the acquisition protocol, the patient set and the proposed algorithm are described. Results are given and discussed in Section 3. Section 4 concludes this paper with a brief description of its contributions and some open questions for future work.

Material and methods

Clinical data and processing environment

The protocol of this study (data collection and post-processing) was approved by our institutional ethical review board. 25 patients were involved in the experiments. All these patients have given their written consent to the participation and received remuneration for it. A non-conflict of interest for this work was declared. CT images were acquired on a multi-detector row Siemens Somatom Sensation 16 CT scanner.

The patient cohort includes 11 women and 14 men with an average age of 64 years (age range: 54-73 years). All the patients suffer from cancers confirmed by biopsy examinations. Both LDCT and SDCT images were collected using a reduced tube current 50mAs and the routine tube current 260mAs under abdomen scan mode. 40 slices were collected for each patient scan. Other scan parameters include: kVp, 120; slice thickness, 5 mm; Gantry rotation time, 0.5s; detector configuration (detector rows section thickness), 16 mm 1.5 mm ; table feed per gantry rotation, 24 mm; pitch, 1:1; reconstruction method: FBP algorithm with convolution kernel "B20f" ("B20f" is the routine smoothing kernel used in reconstructing abdomen images on Siemens CT). Here we consider the abdomen window (center, 50HU; width, 350HU). The volume CT dose index (CTDI vol ) is a linear function of the tube currents [START_REF] Yazdi | Artifacts in Spiral X-ray CT Scanners: Problems and Solutions[END_REF]. We recorded the accumulated doses from the workstation for each scan with 40 slices. The recorded doses are 16.38 mGy for the routine 260 mAs protocol, and 3.09 mGy for the low dose 50 mAs protocol. All the CT images were exported as DICOM files and then processed offline under a PC workstation (Intel Core™ with multiple hepatic metastases and Fig. 1 (b) is the corresponding SDCT image for reference. The LDCT image and SDCT image were respectively acquired with 50mAs and 260mAs, as specified above, and other scanning parameters were both kept as routine setting. We can see that (i) the LDCT image is significantly degraded by mottled noise with strong intensity variation and (ii) the patch-shaped tumor tissues (pointed by arrows) present a low contrast with the surrounding tissues. Therefore, processing LDCT tumor images in abdomen scans should be targeted at improving LDCT image quality and obtaining an enhanced detectability of abdomen tumors.

Method and parameter setting

The proposed DL-based patch processing was performed on 2-D slices because of the large slice thickness setting in routine abdomen CT scanning (>=1mm, 5mm in this work).

Assuming the patches in the target LDCT images are sparsely representable, DL-based patch processing can be carried out by coding each patch as a linear combination of only a few patches in the dictionary [START_REF] Lewicki | Learning overcomplete representations[END_REF][START_REF] Delgado | Dictionary learning algorithms for sparse representation[END_REF]. This method finds the best global over-complete dictionary and represents each image patch as a linear combination of only a few dictionary vectors (atoms). The coefficients of the linear combination are estimated through a sparse coding process [START_REF] Donoho | Maximal sparsity representation via l1 minimization[END_REF]. Based on [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] and [START_REF] Li | An efficient dictionary learning algorithm and its ication to 3-D medical image denoising[END_REF], the DL-based patch processing aims at solving the following problem:

2 2 2 2 0 ,, min s.t. , ij ij ij D ij x R D i j x y x T (1)
where, x and y denote the processed and the original LDCT images respectively, and the subscript ij denotes the pixel index (i, j) in the image. ij R represents the operator that extracts the patch ij x of size n n (centered at (i, j)) from image x. The patch-based dictionary D is a n K matrix, which is composed of K n-vector atoms (columns). Each n-vector column corresponds to one n n patch. denotes the coefficient set

ij ij ij D ij R D i j x T (2) 2 2 2 2 min ij ij ij x RD x y x (3)
The objective of Eq. ( 2) is to train the coefficients and dictionary D from a set of image patches. It can be efficiently solved by k-means singular value decomposition (K-SVD) with the replacement of x by the known observed image y [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]. Starting from an initial dictionary (e.g. the Discrete Cosine Transform dictionary, DCT dictionary), this K-SVD operation estimates and D by alternatively applying the orthogonal marching pursuit (OMP) algorithm and the SVD decomposition [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]. The columns of the dictionary D are constrained to be of unit norm to avoid scaling ambiguity in practical calculation [START_REF] Li | An efficient dictionary learning algorithm and its ication to 3-D medical image denoising[END_REF]. Then, when dictionary D and are obtained, we can obtain the output image x by solving the first order derivative of (3) with respect to each pixel in x:

1 ij ij ij ij ij ij T T R R R x I y D (4)
It was pointed out in [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] that the dictionary trained from the original noisy image itself can lead to better noise suppression than using a global dictionary trained from the other available images in the database. However, in this study, we found that the dictionary trained from a typical SDCT abdomen image always leads to visually close LDCT images when compared to the dictionary trained from the LDCT image itself. The reason might be that most abdomen CT images are composed of similar human tissues in rather low number, and that the dictionary discrepancy caused by differences between different CT abdomen images only leads to small differences in the sparsified features. For one specific LDCT abdomen image, organs or other human tissues can be efficiently represented by the dictionary atoms learned from some other abdomen CT images. In this study, we used the pre-calculated dictionary p D (the left image in Fig. 2) trained from one typical SDCT abdomen image (the right image in Fig. 2) to process all the LDCT images. One important merit of the approach lies in that the intensive computation required in training dictionary is avoided by using this pre-calculated global dictionary. Here, the size of overlapping patches is set to 88 to give an effective representation of local tumor and organ tissues, and the atom number K is set to 64 for it is found large enough to represent the structures in abdomen CT images. Choosing larger K cannot yield result with obvious visual improvement but will significantly increase the computation load. An analysis of the parameter setting in dictionary learning is provided in Section 3.4. With the pre-calculated dictionary p D (obtained via (2)) used as the global dictionary, the whole DL processing can be defined by the following two steps:

2 2 0 (S1), min s.t. , ij ij ij ij p R x i j D (5) 2 2 2 2 (S2), min ij ij ij x p R x y x D (6)
Here, the sparse coefficient set and the image x can be calculated by solving ( 5) and ( 6) using OMP and the solution given in (4). In [START_REF] Mannudeep | Strategies for CT Radiation Dose Optimization[END_REF], denotes the tolerance parameter used in calculating by OMP method. For evaluation purpose, we compared the proposed method with the LNLM method, which is a pixel-wise weighted intensity summation with weights determined by patch similarities [START_REF] Chen | Improving Low-dose Abdominal CT Images by Weighted Intensity Averaging over Large-scale Neighborhoods[END_REF]. The intensive pixel-wise operations in LNLM method was accelerated using GPU (Graphics Processing Unit) techniques based on [START_REF] Chen | Improving Low-dose Abdominal CT Images by Weighted Intensity Averaging over Large-scale Neighborhoods[END_REF]. The parameters involved in the LNLM and proposed DL methods were specified under the guidance of one radiological doctor (X.D.Y. with 15 years of experience) to provide the best visual results. Practically, we found that the same parameter setting can be used to process the LDCT images with the same scan protocol for the two methods. TABLE I lists this parameter setting. For the LNLM method, the parameter h, the patch size Np and the neighborhood size Nn were set to 2, 77 and 81 81, respectively. For the proposed DL method, as shown in Fig. 3, the dictionary size was set to 64 64 for 88 patch (n = 64) and 64 atom numbers (K = 64). Sparsity level 0 L was set to 3 atoms, and 20 iterations (Itern=20) were used in the K-SVD calculation to update the dictionary D and the sparse coefficient . The initial dictionary for K-SVD is the DCT dictionary obtained by sampling the cosine wave functions in different frequencies. The tolerance parameter and in ( 5) and ( 6) were set to 21 and 0.8.

Results

Visual assessment

Here, the complete images are displayed using the standard abdomen window and in full size in order to facilitate an overall quality evaluation. Fig. 2 shows the trained dictionaries (D) used in the implementation of the proposed DL method.

The low dose CT data in Fig. 3-Fig. 5 are respectively from a 61 years female patient, a 53 years female patient, and a 56 years male patient. To be more precise, Fig. 3(a) and Fig. 4(a) are respectively two LDCT images with live tumors, and Fig. 5(a) is a LDCT image with hepatic metastases. (b) and (d) in Fig. 3-Fig. 5 show the corresponding SDCT images and the DL processed LDCT images. Fig. 3(c), Fig. 4(c) and Fig. 5(c) show the LNLM processed images. In Fig. 3-Fig. 5, (e)-(h) illustrate the zoomed regions of interest (ROI) identified by squares in (a). We should note here that, for the SDCT images, the ROIs with the closest visual match with the LDCT images are selected because of the unavoidable displacements between the two different scans. From all (a) in Fig. 3-Fig. 5, we observe that, under low dose CT scanning condition, mottled noise and streak artifacts severely degrade the image quality and lead to unclear tumor boundaries. With the SDCT images as reference, we can see in (c) in Fig. 3-Fig. 5 that the LNLM method can suppress noise and artifacts but tends to introduce some strange striped artifacts and so leads to lowered perceptive detectability of small lesions (arrows in Fig. 5(g)). In particular, referring to the case of multiple hepatic metastases in Fig. 5, we can note that the proposed DL processed image allows a better discrimination of small lesions (arrows in Fig. 5(h)) when compared to the original LDCT image and the LNLM processed image.

Fig. 6 provides the processing results of another five patients with abdomen tumors. The first, second and third columns correspond to the original LDCT images (a-l, b-l, c-l, d-l, e-l), the original SDCT images (a-s, b-s, c-s, d-s, e-s) and the DL processed LDCT images (a-p, b-p, c-p, d-p, e-p). To be specific, in Fig. 6, the first row depicts images with liver metastases (a-l, a-s, a-p); the second row displays images of uterine malignancies (b-l, b-s, b-p); the third row illustrates images with gastric cancer and liver metastases (c-l, c-s, c-p); the fourth row shows images with kidney cancer and liver metastases (d-l, d-s, d-p); the fifth row, images with pelvic malignancies (e-l, e-s, e-p) are given. We can clearly see that a significant improvement of image quality can be obtained by processing the original LDCT images using the proposed DL algorithm. In the processed LDCT images, both noise and artifacts are effectively suppressed, which leads to better visibility of tumor tissues.

By comparing the processed LDCT images with the original SDCT images in Fig. 3-Fig. 6, we can observe that the proposed DL method can produce LDCT image with textures visually closer to those of the original SDCT images. However, in the DL processed LDCT images we can still notice that some original high contrast artifacts still remain (see arrows in Fig. 6 a-p, b-p ande-p). The DL algorithm might be not effective in suppressing such high contrast singular artifacts.

Quantitative assessment

The LDCT and the corresponding SDCT images always have no exact spatial correspondence to each other because of the displacements caused by patient breath and movements in scans. This makes it impossible to quantitatively evaluate image quality using some Euclidean distance metrics (e.g. the mean squared error (MSE)). We thus chose to compare the standard deviation (STD) of the regions of interest (ROI) in both the original and processed LDCT images with respect to those of the SDCT images. The ROIs (denoted by ) for Figures 3 to 5 is the pixel number in .

TABLE II depicts the calculated STDs of tumor and background regions specified in the images in Fig. 7. Tumor-1 and background-1 correspond to the specified regions in the CT images of the first column in Fig. 7; tumor-2 and background-2 to the regions indicated in the second column; tumor-3 and background-3 to the regions defined in the third column of Fig. 7. We can note in TABLE II that the DL processed LDCT images have closer STDs to those of the SDCT images than the LDCT images and also the LNLM processed LDCT images. 8 (a)) for the original LDCT image, the LNLM processed and the DL processed LDCT images. The images were chosen from the case reported in Fig. 4. We also provide the reference histogram map (in red) for the corresponding tumor region in the corresponding SDCT image (selected from the set of slices as the best anatomical match, the regions being of equal size). A larger difference in the histogram maps for the original LDCT and the SDCT images can be observed in Fig. 8 We can also see that, in comparison with the LNLM processing, the DL processed image histogram is rather similar to the original SDCT image one.

Qualitative assessment

For qualitative assessment, 50 original images (including 25 LDCT and 25 SDCT images), 50 processed images (including 25 LNLM processed LDCT and 25 DL processed LDCT images) were considered. 5 subjective features -noise suppression, artifact suppression, contrast preservation, tumor discrimination and overall image quality -were evaluated using a five-point subjective scale (1=unacceptable, 2=substandard, 3=acceptable, 4=above average, 5=excellent). Here, we define artifacts as any tissue region having a passive effect on subjective diagnosis. Three radiological readers (X.D.Y. with 15 years of experience, X.H.Y. with 8 years of experience, Y.M.D. with 5 years of experience.) independently evaluated the randomized set of LDCT images, SDCT images, the LNLM processed LDCT images and the DL processed LDCT images in a digital DICOM archiving/assessing workstation (ViewDEX 2.0 [START_REF] Hakansson | ViewDEX 2.0: a Java-based DICOM-compatible software for observer performance studies[END_REF]). So, the 5 subjective features were assessed for all the 100 images (50 original CT images, and 50 processed CT images). This results in 1500 parameter ratings in total (5 image quality parameters together and the 3 readers, i.e. 100 5 3 1500 ). For each subset of images, the 5 image scores were reported as SDs means

(averaged scores of the 3 radiologists ± standard deviations ). The subjective quality parameters of the original LDCT images and the processed LDCT images were compared with those of the original SDCT images and the differences between each two groups were evaluated by the Student t test (Excel; Microsoft) with P<0.05 considered as statistically significant difference. As illustrated in TABLE III, for all the 5 scores, the original LDCT images are of lower quality than the original SDCT images and the processed LDCT images. The DL processed LDCT images obtained quality scores substantially higher than the LNLM processed LDCT images. Statistically significant differences (P<0.05) with respect to the original SDCT are noticed in all the subjective scores for the LDCT images, and in all the subjective scores except noise suppression for the LNLM processed LDCT images. The differences between the DL processed LDCT images and the original SDCT images for the 5 subjective scores are not statistically significant (P>0.05).

Sensitivity analysis for dictionary training

In this section, we analyze the sensitivity of the proposed DL approach to different dictionary options and also to the atom number K of dictionary. The LDCT image in Fig. 3 is used for this analysis. Fig. 9 (a) illustrates the DL processed image using the dictionary trained from the LDCT image itself. The same parameters, as shown in TABLE I, are used in DL processing. Fig. 9 (b) displays the difference image between Fig. 9 (a) and Fig. 3 (d) (the DL processed image using the global dictionary). The visual as well as the quantitative difference is very small (less than 5 HU). We can thus validate that a global dictionary can lead to almost the same result as the dictionary trained from the LDCT image itself. With respect to the processed LDCT images using global dictionary with larger values of K (from 4 to 441) the root-mean-square deviation (RMSD) is calculated using (8): x denotes each pixel intensity in the DL processed LDCT images using the global dictionary with number K . We can see in Fig. 10 that the RMSD decreases sharply up to K = 50 before going stable. Though not obeying the redundant constraint 64 = n < K given in [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], the moderate atom size K (=64) works well in this study. The reason can be ascribed to the fact that abdomen CT images can always be decomposed into a limited number of tissues with quite similar properties. Fig. 10 The RMSD with respect to the largest atom numbers K for DL processed images using global dictionary. IV lists the computation costs (in seconds) required in implementing different methods for the studies in Fig. 3-Fig. 5. For brevity, we use D-step, O-step, and I-step to represent the K-SVD step ((2) dictionary training), OMP step ((5) sparse coefficient estimation) and the image update step (6) in the DL processing. We can see that the dictionary learning in D-step is much more computationally intensive than the O-step and I-step, and the O-step cost varies a little for different data. TABLE IV shows that, with the D-step avoided by using a pre-trained dictionary (such as the one shown in Fig. 2), the proposed method is computationally much more efficient than the GPU accelerated LNLM method.

Computation costs

Conclusion

This paper proposes to improve abdomen tumor LDCT images by applying a patch-based DL processing. The experiments have shown that the proposed method is clinically effective (5 mm slice thickness) and preserves patch-shaped tumor features with less than 1/5 routine tube current values (or radiation doses). Also, this approach can be of interest for other LDCT abdomen applications.

With no access to well-formatted raw data, the proposed processing can be easily applied to almost all the existing CT systems. It can also be efficiently implemented by using GPU parallelization and a pre-calculated global dictionary. In addition, the inter-patient similarity of abdomen CT images allows a higher acceleration by using dictionary with moderate atom numbers such as 64 in this study. It has been found in our experiments that the same parameter setting is applicable to LDCT images with the same scan protocols, thus offering a tractable strategy for parameter determination. It is expected that the proposed DL approach can also be well used in the processing of other LDCT abdomen data.

Nonetheless, we also found that some artifacts with strong contrast were hard to be suppressed without blurring tumor structures. This effect is expected to be more severe in the case of thinner slices. The reason might be that only limited similarity information within a small neighboring region is used in the proposed DL algorithm. Also, the whole computation cost of the DL processing still needs a further acceleration to fulfill practical clinical requirements (less than 0.5 second per 2-D slice). Currently, some parameters (e.g. the sparsity level and tolerance parameter) still need to be empirically set. Future work will thus be devoted to improvements by incorporating some artifact-suppressing constraints trained from available SDCT images, extending the application to 3D in order to deal with cases with thin slice thicknesses, further accelerating the OMP computation in DL processing by using parallelization techniques, and also developing more robust estimation for parameters.

2

  Quad CPU and 4096 Mb RAM, GPU (NVIDIA GTX465)) with Visual C++ as the developing language (Visual Studio 2008 software; Microsoft).

Fig. 1 A

 1 Fig.1 A typical LDCT image (50mAs) (a) and the corresponding SDCT image (260mAs) (b).

Fig. 1 (

 1 Fig.1 (a) illustrates one typical abdomen LDCT tumor image

  ij ij for all the sparse representations of patches, and each patch ijx can be approximated by a linear combination l norm that counts the nonzero entries of vector ij , and T is the preset parameter of sparsity level that limits the maximum nonzero entry number in ij . Based on[START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], solving (1) includes the following two steps (2) and (3

Fig. 2

 2 Fig.2The left picture is the trained, normalized, global dictionary used in this study with n = 64 ( 88 patch size) and K = 64. The right image corresponds to the abdomen SDCT image from which the global dictionary is trained.

Fig. 3

 3 Fig.3 Results for a 61 years female patient with live tumor. (a), the original LDCT image; (b), the SDCT image; (c), the LNLM-processed LDCT image; (d), the DL processed LDCT image; (e)-(h) show the zoomed regions specified in (a).

Fig. 4

 4 Fig.4 Results for a 53 years female patient with live tumor. (a), the original LDCT image; (b), the SDCT image; (c), the LNLM-processed LDCT image; (d), the DL processed LDCT image; (e)-(h) show the zoomed regions specified in (a).

Fig. 5

 5 Fig.5 Results for a 56 years male patient with multiple hepatic metastases. (a), the original LDCT image; (b), the SDCT image; (c), the LNLM-processed LDCT image; (d), the DL processed LDCT image; (e)-(h) show the zoomed regions specified in (a).

Fig. 6

 6 Fig.6 Five patients (in rows) with abdomen tumors. The first, second and third columns correspond to the original LDCT images (a-l, b-l, c-l, d-l, e-l), the original SDCT images (a-s, b-s, c-s, d-s, e-s) and the DL processed LDCT images (a-p, b-p, c-p, d-p, e-p). The first row illustrates the images with liver metastases (circles); the second row, textures of uterine malignancies; the third row, images corresponding to gastric cancer (circles) and liver metastases (arrows); the fourth row, CT images with liver metastases (arrows) from kidney cancer; the fifth row, CT images with pelvic malignancies (circles) .Note that the arrows in the first, second and third rows point out some remaining artifacts in the processed images.

  include both tumor regions and backgrounds specified by a radiological doctor (X.D.Y. with 15 years of experience). Fig.7 illustrates the ROIs for tumor regions (surrounded by red circles) and backgrounds (dark circles) for the original SDCT and LDCT images. The standard deviation Ω pixel intensity and the averaged intensity within , respectively.

Fig. 7

 7 Fig.7 The three columns from left to right illustrate the selected ROI ( ) for tumor regions (tumor region-1, tumor region-2, tumor region-3, red circles) and background regions (background-1, background-2, background-3, black circles) in the abdomen CT images corresponding to Fig.3, Fig.4 and Fig.5, respectively. In each column, the upper and lower images denote the SDCT image and the corresponding LDCT image, respectively.

Fig. 8

 8 Fig.8The histogram maps (black) of the tumor region (specified in Fig.4 (a)) for the original LDCT image, the LNLM processed and the DL processed LDCT images. The histogram map of the corresponding tumor region in the original SDCT image is overlaid in red color.

Fig. 8 (

 8 Fig.8 (b), (c) and (d) plot the histogram maps (in black) of the specified tumor region (the red frame in Fig.8(a)) for the original LDCT image, the LNLM processed and the DL processed LDCT images. The images were chosen from the case reported in Fig.4. We also provide the reference histogram map (in red) for the corresponding tumor region in the corresponding SDCT image (selected from the set of slices as the best anatomical match, the regions being of equal size). A larger difference in the histogram maps for the original LDCT and the SDCT images can be observed in Fig.8We can also see that, in comparison with the LNLM processing, the DL processed image histogram is rather similar to the original SDCT image one.

Fig. 9 (

 9 Fig.9 (a). The DL processed LDCT image using the dictionary trained from the LDCT image itself; (b). The difference image between Fig.9 (a) and Fig.3 (d).

  Fig.3 study

  (a) Processed LDCT image using self-trained dictionary (b) Difference image between Fig.9 (a) and Fig.3(d) RMSD SD KS

TABLE II CALCULATED

 II STDS OF ROI FOR THE ORIGINAL LDCT IMAGES, THE ORIGINAL SDCT IMAGES , AND THE PROCESSED LDCT IMAGES

				LNLM	DL
	ROI	LDCT	SDCT	processed	processed
				LDCT	LDCT
	Tumor-1	235.02	46.06	95.50	62.31
	Background-1	272.99	78.03	102.27	82.86
	Tumor-2	108.18	23..80	47.20	25.04
	Background-2	106.63	21.94	55.80	25.32
	Tumor-3	219.09	29.41	84.45	26.28
	Background-3	314.98	51.74	95.56	65.84

TABLE

TABLE IV THE

 IV COMPUTATION COST (IN CPU SECONDS) FOR THE LNLM METHOD IN[START_REF] Chen | Improving Low-dose Abdominal CT Images by Weighted Intensity Averaging over Large-scale Neighborhoods[END_REF]

Acknowledgments

This research was supported by National Basic Research Program of China under grant (2010CB732503), National Natural Science Foundation under grants (81000636), and the Project supported by Natural Science Foundations of Jiangsu Province (BK2011593).


Chen Yang

Fei Yu

Limin Luo

Christine Toumoulin

Yang Chen

Yu Fei

Improving abdomen tumor low-dose CT images using dictionary learning based patch processing and unsharp filtering

Keywords: Low-dose CT (LDCT), abdomen tumor, Gaussian kernel, preprocessing, learning dictionary

Reducing patient radiation dose, while maintaining a highquality image, is a major challenge in Computed Tomography. The purpose of this work is to improve abdomen tumor low-dose CT (LDCT) image quality by using a two-step strategy: a first patch-wise non linear processing is first applied to remove the noise, that is based on a sparsity prior in term of a learned dictionary, then an unsharp filtering aims to enhance the contrast of tissues and compensate the contrast loss caused by the DL processing. Preliminary results show that the proposed method is effective in suppressing mottled noise as well as improving tumor detectability.

INTRODUCTION

CT imaging is increasingly incorporated into clinical decision making and despite rapid progresses in CT technology over the past decade, one major concern appears today related to the associated radiation rate rising [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF][START_REF] Brenner | Computed tomography-an increasing source of radiation exposure[END_REF]. A large number of researches in CT have been motivating by the need to reduce patient radiation dose. Among the possible solutions, one is to consider lowering the X-ray tube current. Nevertheless, low dose CT provides degraded images by mottled noise and different kinds of non stationary artifacts [START_REF] Mannudeep | Strategies for CT Radiation Dose Optimization[END_REF][START_REF] Yazdi | Artifacts in Spiral X-ray CT Scanners: Problems and Solutions[END_REF], which render the interpretation of these images particularly difficult. Tumor tissues often thus appear under the form of mosaic shapes within a low contrasted environment [START_REF] Watanabe | Improvement of image quality of low radiation dose abdominal CT by increasing contrast enhancement[END_REF][START_REF] Funama | Improvement of low-contrast detectability in lowdose hepatic multidetector computed tomography using a novel adaptive filter: evaluation with a computer-simulated liver including tumors[END_REF]. Two kinds of methods are applied to enhance image quality. They, either, directly proceed in the reconstruction domain or within a postprocessing denoising stage. In both cases, efficient noise suppression and tumor tissue preservation remains challenging. Neighborhood filters have shown interested properties for the restoration of noisy low dose CT images.

Let cite for instance, adaptive filters [START_REF] Funama | Improvement of low-contrast detectability in lowdose hepatic multidetector computed tomography using a novel adaptive filter: evaluation with a computer-simulated liver including tumors[END_REF] that allow to reduce the X-ray dose by 50% for the same image quality and without loss in low contrast detectability. Other filters such as multiscale penalized weighted least-squares [START_REF] Wang | Multiscale penalized weighted least-squares sinogram restoration for low-dose X-Ray computed tomography[END_REF], bilateral filters [START_REF] Ehman | Noise reduction to decrease radiation dose and improve conspicuity of hepatic lesions at contrast-enhanced 80-kV hepatic CT using projection space denoising[END_REF] and Non Local Mean (LNLM) [START_REF] Chen | Improving Low-dose Abdominal CT Images by Weighted Intensity Averaging over Large-scale Neighborhoods[END_REF] have also shown some efficiency in enhancing anatomical/ pathological features in Low dose CT images.

Recent years have reported a growing interest in the study of sparse representation based dictionary learning (DL) [START_REF] Lewic | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF][START_REF] Lewicki | Learning overcomplete representations[END_REF][START_REF] Delgado | Dictionary learning algorithms for sparse representation[END_REF][START_REF] Donoho | Maximal sparsity representation via l1 minimization[END_REF][START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]. Compared to pixel-wise intensity update-based restoration methods, patch-wise DL processing are considered as being more robust to mottled noise and generally provides a more efficient representation of patchshaped features such as tumors or organs. We describe in this paper, a new patch-wise processing, based on a sparsity prior in terms of a learned dictionary to suppress mottled noises in abdomen tumor LDCT images and a contrast-enhancing unsharp filter whose role is to compensate the contrast loss induced by the DL process. This method referred as DL-unsharp algorithm, is described in section 2. The flowchart of the method is given in Fig. 1. Section 3 provides a comparative study between our algorithm and a LNLM restoration filter [START_REF] Chen | Improving Low-dose Abdominal CT Images by Weighted Intensity Averaging over Large-scale Neighborhoods[END_REF]. Preliminary results show that the proposed DL-unsharp algorithm provides a good restoration of structures in LDCT images with an image quality that is comparable with the original SDCT images.

METHOD

The core idea is to impose a sparsity prior patch-wise on the LDCT images in terms of a dictionary D. Assuming the patches in the LDCT image are sparsely representable, DL based patch processing is carried out by coding each patch as a linear combination of just a few patches in the dictionary i.e. each patch of the image can be approximated by a linear combination of just a few columns from D [START_REF] Lewicki | Learning overcomplete representations[END_REF][START_REF] Delgado | Dictionary learning algorithms for sparse representation[END_REF]. This way to proceed leads to find the best global overcomplete dictionary. The coefficients of the linear combination can be estimated through the sparse coding process described in [START_REF] Donoho | Maximal sparsity representation via l1 minimization[END_REF].

The DL based patch processing aims to solve the following optimization problem [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]:

2 0 2 0 2 2, 0 ,, min s.t. , ij ij ij D ij x R D L i j x y x           (1)
where, x and y denote the treated and original LDCT images, respectively. ij R is the operator that extracts the square The numerical solution of the optimization problem ( 1) is obtained by a weighted version of the K-SVD algorithm [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]. It consists of two main steps:

2 0 2, 0 , min s.t. , ij ij ij D ij R D T i j x          (2) 2 2 2 2, min ij ij ij x RD x y x       (3)
(2) aims to train the dictionary D and  from a set of image patches. It is solved with the K-means Singular Value Decomposition (K-SVD) after replacing x by the known observed image y. This operation is iteratively performed in two steps: (1) sparse coding of  (including all   ij ij  ) using the orthogonal marching pursuit (OMP) algorithm; (2) dictionary update by minimizing (2) with D being a matrix with unit-norm columns in order to avoid scaling ambiguity [START_REF] Li | An efficient dictionary learning algorithm and its application to 3-D medical image denoisin g[END_REF]. Given the dictionary D, we compute for each training image, the patch-wise representations  using the OMP algorithm. Finally, given D and  , we compute the output image x by solving (3) according to a simple least squares approach:

1 ij ij ij ij ij ij ij ij T T R R R x I y Dx                         (4)
Contrary to what was pointed out in the literature, we found that the dictionary trained from a SDCT abdomen image always provided results that were visually close to LDCT images (when compared with the dictionary trained from the LDCT image itself). The reason is that most abdomen CT images have similar tissue compositions and the dictionary discrepancy often expresses a very few difference in the final sparsified features. In consequence, we decided to use a pre-computed general dictionary p D (cf. Fig. 2.a for illustration) that was preliminary trained from a high quality SDCT abdomen reference image (Fig. 2.b). One interesting advantage of using this general dictionary is that the intensive computations, involved in the dictionary construction, are avoided. We also improved the original K-SVD algorithm by associating a weighting to each pixel in the extracted patches, which is function of its distance to the center of the patch (a large weight is given to points that are close to the centre of the window and decreases as the distance of the point increases). We finally performed the optimization process using the dictionary p D obtained from ( 2) and ( 3) and considering the following three steps:

2 2 2, 2, (S1), min s.t. , ij ij ij ij ij p R D R x Dp i j x             (5) 2 2 2 2, (S2), min ij ij ij x p RD x y x       (6)   (S3), unsharpFilter , = p x x  (7)
The sparse coefficients  and DL treated image x are calculated using ( 5) and ( 6). The  in (5) denotes the tolerance parameter used in the computation of the  coefficients.

Step S3 (7) characterizes the final contrast enhancing unsharp filtering with the  -weight kernel [16]:

- - - - 1+8  - - - -

EXPERIMENT

Approval of this study was granted by our institutional review board. A non-conflict of interest for this work was declared. SDCT images of abdomen were acquired on a multi-detector row Siemens Somatom Sensation 16 CT scanner with a current of 260mAs while LDCT images were obtained from a reduced current of 50 mAs. The scanning parameters were the following: kVp, 120; slice thickness, 5 mm; Gantry rotation time, 0.5s; detector configuration (detector rows  section thickness), 16mm 1.5mm



; table feed per gantry rotation, 24 mm; pitch, 1:1; reconstruction method: FBP algorithm with convolution kernel "B31f" ("B31f" is the routine smoothing kernel used in abdomen scans for Siemens CT). The windows and level setting were chosen to optimize the visualization of these data (center, 50HU; width, 350HU). For evaluation purpose, we compared the proposed method with the LNLM method in [START_REF] Chen | Improving Low-dose Abdominal CT Images by Weighted Intensity Averaging over Large-scale Neighborhoods[END_REF]. The LNLM method was accelerated using GPU (Graphics Unit) techniques based on [START_REF] Chen | Improving Low-dose Abdominal CT Images by Weighted Intensity Averaging over Large-scale Neighborhoods[END_REF]. All the CT images were exported as DICOM files and then processed offline on a PC workstation (Intel Core™ 2 Quad CPU and 4096 Mb RAM, GPU (NVIDIA GTX465)) using Visual C++ as programming (Visual Studio 2008 software; Microsoft).

For both algorithms, the parameter setting was completed applying a greedy algorithm to find the optimal parameter that provided the best qualitative results. This qualitative evaluation was carried out in collaboration with a radiologist (X. D.Y, 15 years clinical experience). These optimal parameters are listed in TABLE I with the computation time costs for each method. To specify the computation cost at each different stages of the proposed DL-unsharp processing, we used K-step, Ostep, I-step and F-step techniques to represent the K-SVD stage (2)-(3) (dictionary training), OMP step (5) (sparse coefficient estimation), the image update stage [START_REF] Funama | Improvement of low-contrast detectability in lowdose hepatic multidetector computed tomography using a novel adaptive filter: evaluation with a computer-simulated liver including tumors[END_REF], and the unsharp filtering step [START_REF] Wang | Multiscale penalized weighted least-squares sinogram restoration for low-dose X-Ray computed tomography[END_REF], respectively. We see in Table I that the dictionary learning in the K-step method is computationally more time consuming than the O-step and Istep techniques. Si if we remove the K-step and replace it by a pre-trained global dictionary (Fig. 3(a)), the proposed implementation appears thus more efficient in computation time (2.28+ 0.96 +0.12 = 3.36 seconds) than the LNLM method (8.07 seconds). . Their restoration provides a visual appearance of the texture that is close to those of the SDCT images. In the one including multiple hepatic metastases in Fig. 4, the DLunsharp algorithm allows enhancing the small structures as the small lesions which appears better discriminated than in the LNLM processed images (see arrows).

TABLE I. PARAMETER SETTING AND COMPUTATION COST (IN

CONCLUSION

The algorithm described here and named DL-unsharp is employed for improving abdomen LDCT image quality, the objective being both to suppress the mottled noise and streak artifacts while enhancing the structure edges especially on tumors or lesions. This method makes use of a patch based DL processing followed by a contrast restoration unsharp filtering. Furthermore, the dictionary training has been built from a set of abdomen SDCT images to optimize the algorithm performance. We demonstrated the potential of the proposed approach on abdomen tumor LDCT datasets. Experiment results showed the proposed approach can greatly improve the quality of images acquired with a reduced X-Ray dose by 80%.

However, some improvements are needed: First, the whole computation cost of the DL-unsharp processing still need to be accelerated to meet the clinical requirement (less than 0.5 second per image). Second some parameters are currently set empirically and need more experiments to validate their value. Thirdly, extensive experiments have to be led to confirm these preliminary results. In conclusion, future work will be devoted to all these points: parallelization of pair-wise operations, automatic estimation of the best parameters to optimize the tradeoff between data fidelity and priors and extensive evaluation.

Fig 1 .

 1 Fig 1. Outline of the proposed DL-unsharp algorithm.

  patch ij x of size n n  (centered at point (i, j)) from the image x. This patch is encoded by D ij . D is a n K  matrix, which is composed by K columns of n-vectors. Each nvector column corresponds to a patch of size n n  .  denotes the coefficient set   ij  for the sparse representation of all patches.  denotes the standard deviation of the Gaussian kernel which assigns a large weight to points that are closer to the center of the patch. It is defined by ij Rx and counts the nonzero entries of vector ij  , and 0 L is the sparsity level that limits the maximum nonzero entry number in ij  .

  a) Dictionary example; (b) Abdomen SDCT image from which the dictionary has been trained.

Fig. 3

 3 Fig.3 and 4illustrate the results for two patient datasets. Fig.3(a) and Fig.4(a) depict two abdomen LDCT images including tumors (pointed by red circles or arrows) of a 61 years old female and 56 years old male patient respectively. Fig.3and Fig.4 (b), (c) and (d) show the corresponding SDCT, LNLM processed LDCT and DL-unsharp processed LDCT images respectively. We observe that, under low dose scanning condition, mottled noise severely degrades the images and the tumor boundary appears fuzzy. Considering the SDCT images as references, we observe that the LNLM processed images (Fig.3 (c) and Fig.4 (c)) have been smoothed but still contain noise and stripe artifacts. The result appears more convincing with the DL-unsharp method since we can observe a more efficient noise reduction with a good preservation of tumor structures (Fig.3 (d) and Fig.4(d)). Their restoration provides a visual appearance of the texture that is close to those of the SDCT images. In the one including multiple hepatic metastases in Fig.4, the DLunsharp algorithm allows enhancing the small structures as the small lesions which appears better discriminated than in the LNLM processed images (see arrows).

Fig. 3

 3 Fig.3 Results for a dataset of a 61 years female patient having a liver tumor (red circles). (a): original LDCT image; (b): SDCT image; (c): LNLM processed LDCT image; (d): DL-unsharp processed LDCT image.

Fig. 4

 4 Fig.4 Results for a dataset of one 56 years male patient having multiple hepatic metastases (red arrows) in the abdomen. (a): original LDCT image; (b): SDCT image; (c): LNLM processed LDCT image; (d): DL-unsharp processed LDCT image.

ACKNOWLEDGEMENTS

This research was jointly supported by National Basic

Research Program of China under grant (2010CB732503), National Natural Science Foundation under grant 81000636, and the Natural Science Foundations of Jiangsu Province (BK2011593) and the French National Institute for Health (INSERM).