Keywords: compressive sensing, image reconstruction, regularization, total variation, augmented Lagrangian, nonlocal self-similarity, wavelet denoising

since the beginning and for being by my side at the difficult moments of this thesis with a lot of love and understanding.

Résumé

Dans les applications de télédétection et d'imagerie médicale, l'un des points clés est l'acquisition, le prétraitement en temps réel et le stockage des informations. En raison de la grande quantité d'informations présentes sous forme d'images ou de vidéos, la compression de ces données est nécessaire. La détection comprimée est une technique efficace pour relever ce défi. Elle consiste à acquérir un signal, en supposant qu'il puisse avoir une représentation éparse, en utilisant un nombre minimal de mesures linéaires non adaptatives. Après ce processus de détection comprimée, une reconstruction du signal original doit être effectuée au niveau du récepteur.

Les techniques de reconstruction sont souvent incapables de préserver les textures de l'image et ont tendance à en lisser les détails. Pour surmonter ce problème, nous proposons, dans ce travail, une méthode de reconstruction par acquisition compréssée qui combine la régularisation par variation totale et la contrainte d'autosimilarité non locale. L'optimisation de cette méthode est réalisée à l'aide du Lagrangien augmenté qui évite le problème de la non-linéarité et de la non-différenciabilité des termes de régularisation. L'algorithme proposé, appelé denoising-compressed sensing by regularization (DCSR) terms, permet non seulement de reconstruire l'image mais aussi de la débruiter. Pour évaluer les performances de l'algorithme proposé, nous le comparons à quelques méthodes , telles que l'algorithme de Nesterov, la représentation clairsemée par groupes et les méthodes basées sur les ondelettes, en termes de débruitage et de préservation des bords, de la texture et des détails de l'image, ainsi que du point de vue de la complexité de calcul.

Notre approche permet un gain jusqu'à 25% en termes d'efficacité de débruitage et de qualité visuelle en utilisant deux métriques : le rapport signal/bruit maximal et la similarité structurelle.

Mots clés: acquisition comprimée; reconstruction d'image; régularisation; variation totale; Lagrangien augmenté; autosimilarité non locale; débruitage par ondelettes.

Publications

Journal Papers

A. El Mahdaoui, A.Ouahabi, M.S.Moulay, . "Image Denoising Using a Compressive Sensing Approach Based on Regularization Constraints," Sensors 2022, https://doi.org/10.3390/s22062199

International Conferences

A. El Mahdaoui, M.S.Moulay, A.Ouahabi. " Image Recovery Using Total Variation Minimization on Compressive Sensing", " international conference on image and signal processing and their applications (ISPA 2019). IEEE. [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]

Introduction

Image processing is undergoing considerable expansion with the emergence of more advanced technology nowadays. Indeed, while on the one hand computers, televisions, mobile phones and other consumer products, whose commercial and economic stakes are immense, bring new challenges, more "fundamental" areas also raise many problems. For example, we can mention medicine and the use of X-rays [START_REF] Zehani | Trabecular bone microarchitecture characterization based on fractal model in spatial frequency domain imaging[END_REF], MRI [START_REF] Atia | Particle swarm optimization and two-way fixed-effects analysis of variance for efficient brain tumor segmentation[END_REF] or ultrasound [START_REF] Girault | Estimation of the blood Doppler frequency shift by a time-varying parametric approach[END_REF], [START_REF] Girault | Analytical formulation of the fractal dimension of filtered stochastic signals[END_REF] or ultrasound, the observation of the infinitely small,of the infinitely large ...

In each case mentioned above, the challenge is to extract the maximum amount of information from the image. This generally requires special processing: restoration, decomposition, segmentation, etc.

Compressed sensing (CS) has already attracted great interest in various fields. Examples include medical imaging [START_REF] Yang | Fast and accurate compressed sensing model in magnetic resonance imaging with median filter and split Bregman method[END_REF], [START_REF] Labat | Harmonic propagation of finite amplitude sound beams: Experimental determination of the non linearity parameter B/A[END_REF], communication systems [START_REF] He | Compressive multi-attribute data gathering using hankel matrix in wireless sensor networks[END_REF] - [START_REF] Haneche | Compressed sensing-speech coding scheme for mobile communications[END_REF]), remote sensing [START_REF] Li | Compressed sensing imaging with compensation of motion errors for MIMO Radar[END_REF], radar imaging [START_REF] El Mahdaoui | Multilevel fast multipole acceleration for Fast ISAR imaging based on compressive sensing[END_REF],

reconstruction algorithm design [START_REF] Andras | A time domain reconstruction method of randomly sampled frequency sparse signal[END_REF], image storage in databases [START_REF] Mimouna | A heterogeneous multimodal data set for advanced environment perception[END_REF], etc. Compressed sensing provides an alternative approach to Shannon's vision to reduce the number of samples and/or reduce transmission/storage costs. There are other approaches that also address this issue, such as random sampling [START_REF] Ouahabi | Spectrum estimation from randomly sampled velocity data LDV[END_REF]. Compressed sensing recovery is a linear optimization problem. The most common CS retrieval algorithms explore the prior knowledge that a natural image is sparse in certain domains, such as in the wavelet domain, where simple and efficient noise reduction is possible [START_REF] Ouahabi | A review of wavelet denoising in medical imaging[END_REF]- [START_REF] Femmam | Perception and characterization of materials using signal processing techniques[END_REF], or in the discrete gradient domain, which we will develop in this work.

In real life, images are usually noisy, but noise is random and is, therefore, unknown.This noise can have many origins: It can be due to poor weather conditions (wind, haze, fog, mist,...), light fluctuations, the electronic image sensor of a digital camera, the conditions under which the image was acquired or simply the manner the image was stored and the techniques used to compress it.

Therefore, denoising is an essential pre-processing step to recover improved image quality. We can consider the problem of image recovery as an inverse problem because the goal is to reconstruct an image close to the original image from a degraded image while respecting two constraints: an optimized image quality and the speed of execution. Our strategy is part of this framework, which consists in recovering a quality image a debatable notion that is imprecise and depends not only on objective criteria but also on the "eye" of the observer in a relatively short time. The regularization of an inverse problem corresponds to the idea that data alone do not allow obtaining an acceptable solution and that it is, thus, necessary to introduce a priori information on the regularity of the image to be estimated, reconstructed or re-covered. The regularization of an inverse problem corresponds to the idea that data alone cannot make obtaining an acceptable solution possible and that it is, therefore, necessary to introduce a priori information, likely allowing an estimation, reconstruction or recovery of the image of interest.

Many inverse problem optimization approaches for image denoising have been proposed in the literature. Some are based on deep learning and more precisely on the deep generative network [START_REF] Chen | Image denoising using a novel deep generative network with multiple target images and adaptive termination condition[END_REF], and others are based on models [START_REF] Zha | Group sparsity residual constraint with non-local priors for image restoration[END_REF].

In particular, the total variation (TV)-based approach has been one of the most popular and successfully applied approach, where, for example, Chambolle [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] introduced the dual approach to the unconstrained real case. Subsequently, Beck and Teboulle [START_REF] Beck | Fast Gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] presented a fast computational method based on a gradient optimization approach to solve the TV-regularized problem. Recently, other denoising approaches have been proposed, such as a remote sensing image denoising method via a low-rank tensor approximation [START_REF] Ma | Remote sensing image denoising via low-rank tensor approximation and robust noise modeling[END_REF], which can be formulated as a generative Bayesian model. Another method for regularizing nonuniformly sampled data based on least-squares spectral analysis [START_REF] Ghaderpour | Multichannel antileakage least-squares spectral analysis for seismic data regularization beyond aliasing[END_REF] can be applied to non stationary signals by introducing classical sliding windowing, as was the case for processing marine seismic data. Other approaches to denoising hyperspectral images [START_REF] Zhang | Hyperspectral image denoising with total variation regularization and nonlocal low-rank tensor decomposition[END_REF], [START_REF] Wang | Tensor low-rank constraint and l0 total variation for hyperspectral image mixed noise removal[END_REF] combine total variation regularization and low-rank tensor decomposition.

However, since the total variation model favors piecewise constant image structures, the total variation models tend to oversmooth image details; it tends to smooth out the fine details of an image. To overcome these intrinsic drawbacks of the total variation model, we introduce a nonlocal self-similarity constraint as a complementary regularization. Nonlocal self-similarity can restore high-quality images. To make our algorithm robust, an augmented Lagrangian method was used efficiently to solve the above inverse problem.

Extensive numerical results are presented, obtained from experiments. The results establish the feasibility and robustness, to various extent, of the proposed data processing schemes, models and algorithms.

There remain many challenges to be further resolved in each area. But, still hope the progress made in this thesis will represent a useful first step toward meeting these challenges in the future.

Thesis Outline

The remainder of thesis is organized as follows:

Chapter 1: "Notation and Preliminaries"

This chapter is devoted to preliminaries and reminders, an overview of image characteristics, some types of deterioration image noise, and image denoising.

Chapter 2: "Compressed Sensing"

Chapter 2, provides the basic concepts of CS. The mathematical formulation is given, in addition to an overview of its applications in many areas, including those in communications.

Chapter 3: "Inverse problems and solving methods"

Introduction

We present the main problem of interest, inverse problems and illustrate the ill-posedness nature of some of them. Next, we will explain how, starting from an inverse problem, we led to an optimization problem. We also review some of the most common reconstruction methods encountered in the literature.

Chapter 4: "Image denoising using a compressive sensing approach based on regularization constraints"

Describes our algorithm's denoising-compressed sensing by regularization (DCSR) and under certain assumptions, a theorem about the uniqueness of the problem's solution is founded and approved.

Also, it contains descriptions and results of numerical experiments to demonstrate the proposed methods or framework's efficiency, effectiveness, and limitations.

Finally, we conclude the whole thesis, points out the future work and a bibliography.

Chapter 1

Notation and Preliminaries

This chapter briefly recalls some of the definitions and analytical tools that will be used frequently in the following chapters and also, introduces some generalities on image processing.

Notations

f is a function defined on Ω an open of R n .
• ∇ is the gradient operator, ∇f = ( ∂f ∂x i ) i=1,..,n

• D u is the derivative operator in the direction u ∈ R n , D u f =< ∇f, u >, where < •, • > is the inner product.

• div(u) is the divergence operator defined for u ∈ C 1 (Ω, R n ) as:

div(u) = n i=1 ( ∂u i ∂x i ) • 1 ≤ p ≤ ∞, L p (Ω) = {f : Ω → R, f measurable and Ω |f (x)| p < ∞} • f p is the norm in L p (Ω), f p = ( Ω |f (x)| p dx) 1 p
• Let p ≥ 1 and Ω denote the set of all real-valued sequences. The space l p is the set of all real sequences x = (x i ) ∈ Ω such that

x p = ( ∞ i=1 (|x i | p ) 1 p < ∞ and x ∞ = sup i (|x i |)
Chapter 1

Preliminaries and Generality

Analysis notion

We recall hereafter some notion of functional and convex analysis, valid for the continuation. We invite a reader to references [START_REF] Brezis | Opérateurs maximaux monotones et semi groupes de contractions dans les espaces de Hilbert[END_REF]- [START_REF] Ekeland | Convex analysis and variational problems[END_REF] for more details.

Definition 1.1. (Hilbert space) A Hilbert space H is a vector space over R or C with inner product noted < •, • >, such that it is a Banach space (i.e. a complete metric space) with the norm inducted by the inner product, that is:

x = √ < x, x >, ∀x ∈ H (1.1)
Properties 1.1. Let H be Hilbert space.

1. If H is a real space, then

∀x ∈ H, < x, x >≥ 0 and x = √ < x, x > (1.2)
2. Cauchy-Schwarz inequality:

∀x, y ∈ H, |< x, y >|≤ x y (1.3) 3. Parallelogram identity ∀x, y ∈ H : x 2 + y 2 = 1 2 ( x + y 2 + x -y 2 ) (1.4)
4. The orthogonal complement of a set S is the closed subspace

S ⊥ = {x ∈ H : x ⊥ s, ∀s ∈ S} (1.5)
5. Let x, y ∈ H. We say that the two vectors x and y are orthogonal if the inner product:

< x, y >= 0 and we put x ⊥ y.

6. If x ⊥ y, then we have the Pythagorean theorem

x + y 2 = x 2 + y 2 (1.6) 7.
Let (e n ) n∈N an orthogonal basis in H.

• For any x ∈ H : x = n≥0 a n e n , a n =< x, e n > Chapter 1

Preliminaries and Generality

• For any x ∈ H :

x 2 = n≥0 |a n | 2 ,
• For any x, y ∈ H with x = ∞ n=1 a n e n , y = ∞ n=1 b n e n we have:

< x, y >= ∞ n=1 a n • b n
As an example, the space L 2 ([a, b]), (a set of measurable functions from an interval [a, b] ⊂ R to R or C, whose squares are Lebesgue integrable) is a Hilbert space with the following inner product

< f, g >= b a f (t)g(t)dt, ∀ f, g ∈ L 2 ([a, b]) (1.
A = sup x∈H 1 ,x =0 Ax x (1.8) is finite. Exemple 1.1. The bounded linear operator is an integral operator A : L 2 ([a, b]) -→ L 2 ([c, d]) as: Af := b a K(x, s)f (s)ds (1.9)
Definition 1.3 (Adjoint operator). et H 1 and H 2 be Hilbert spaces and A : H 1 -→ H 2 be a bounded linear operator. Then, the adjoint operator A + of A is an operator

A + : H 2 -→ H 1 satisfies: ∀x ∈ H 1 , ∀y ∈ H 2 :< Ax, y >=< x, A + y > (1.10)
In the following, we recall some notion of convex analysis [START_REF] Brezis | Functional analysis, sobolev spaces and partial differential equations[END_REF]. Let X an Euclidian space equipped with inner product < •, • > and euclidian norm • .

The following notions could generalized to Hilbert spaces.

Let f : X -→ R = [-∞, +∞] with the effective domain:

dom(f ) = {x ∈ X : f (x) < +∞} (1.11) Definition 1.4. (Coercivity) • f is called coercive on a subset K ⊂ X if: For (x n ) n a sequence in K, x n -→ +∞ then f (x n ) -→ +∞.
Chapter 1

Preliminaries and Generality

• Furthermore, it is strongly convex if :

For (x n ) n a sequence in K, x n -→ +∞ then f (xn) xn -→ +∞. Definition 1.5. (Convexity)
• f is convex if and only if:

∀x, y ∈ X; t ∈ [0, 1] : f (tx + (1 -t)y) ≤ tf (x) + (1 -t)f (y).
(1.12)

• f is strictly convex if: ∀x, y ∈ X, x = y, t ∈]0, 1[: f (tx + (1 -t)y) < tf (x) + (1 -t)f (y). (1.13)
• f is proper if f is not identically equal to +∞ and has no value equal to -∞.

• f is lower-semicontinuous on X if: for any sequence (x n ) n converging to x in X, we have:

f (x) ≤ lim inf n-→∞ f (x n ). (1.14)
We will state several interesting properties of convex functions. These properties make them particularly well suited to minimization problems.

Proposition 1.1 (Properties of convex function). Let X and Y be euclidian sets,f : X -→ R and g : Y -→ R convex functions. We have the following:

1. αf is convex for any α > 0.

2. f + g is convex.

for φ

: Y -→ X affine linear, f • φ is convex on Y . (i.e.: Φ(λu + (1 -λ)v) = λΦ(u) + (1 -λ)Φ(v), ∀λ ∈ [0, 1], ∀u, v ∈ Y ) 4. let (h i ) i∈I a family of convex functions, h i : X -→ R then: h(x) = sup i∈I h i (x) for any x ∈ X, is convex. (h is a pointwise supremum) Definition 1.6. (Sub-differential) Let f : X -→ R = [-∞, +∞], be a convex function. The sub-differential of f is defined as: ∀x ∈ X, ∂f (x) = {y ∈ X : f (z) ≥ f (x)+ < y, z -x >, ∀z ∈ dom(f )}. (1.15) Remark 1.1. • If f is differentiable at x ∈ X, then ∂f (x) = ∇f (x)
Chapter 1

Preliminaries and Generality

• For x ∈ X, ∂f (x) is a closed convex subset possibly empty.

Denote by: intC = interior of C, intrC = relative interior of C where C is a convex subset.

Theorem 1.1. (Properties of sub-differentials) Let: X, Y be real normed spaces and f, g two functions: X -→ R proper and convex, A : Y -→ X linear and continuous operator.

Then, the sub-differential obeys the following rules:

• ∂(λf )(x) = λ∂f (x) , ∀λ > 0. • ∂(f + g)(x) ⊃ ∂f (x) + ∂g(x) on ∀x ∈ X. • ∂(f + g)(x) = ∂f (x) + ∂g(x), if intr(domf ) ∩ intr(domg) = ∅. • ∂(f • A)(x) ⊃ A + ∂f (Ax), ∀x ∈ X. • ∂(f • A)(x) = A + ∂f (Ax) if f is continuous at some point x 0 ∈ rg(A) ∩ dom(f ) , such that rg(A) = A(Y ) and A + the adjoint operator of A.
Proof For the complete proof, see [START_REF] Gilbert | Fragments d'optimisation différentiable-théorie et algorithmes[END_REF] Chapter 3, Proposition 3.70 and Proposition 3.72 .

An important notion of convex analysis, is:

Definition 1.7. (Convex conjugate of a function) Let f : X -→ R = [-∞, +∞].
The convex conjugate noted f * of f is the function f * : X -→ R, defined as:

f * (x ) = sup x∈X {< x , x > -f (x)}, x ∈ X (1.16) Remark 1.2.
1. The convex conjugate of a function f is always lower semi continues (l.s.c).

2. The biconjugate f * * (the convex conjugate of the convex conjugate f * ) is the closed convex hull (convex envelope) with: f * * ≤ f .

3. For a proper function f , we have:

f * * = f ⇐⇒ f is convex and l.s.c.
We have the following relation ship between the sub-differential and the convex conjugate of some function.

Theorem 1.2. Let f : X -→ [-∞, +∞] be proper, then:

∂f (u) = {v ∈ X : f (u) + f * (v) =< v, u >} (1.17)
If f is also convex and l.s.c, then:

∂f * (v) = {u ∈ X : f (u) + f * (v) =< v, u >} (1.18)
Chapter 1

Preliminaries and Generality

Proof From the definition of the sub-differential, we have:

v ∈ ∂f (u) ⇐⇒ f (z) ≥ f (u)+ < v, z -u > ∀z ∈ X, ⇐⇒< v, u > -f (u) ≥< v, z > -f (z) ∀z ∈ X, ⇐⇒< v, u > -f (u) ≥ sup z∈X < v, z > -f (z)
, which shows that the sup in z is reached at point u, hence:

v ∈ ∂f (u) ⇐⇒< v, u > -f (u) = sup z∈X < v, z > -f (z) = f * (v),
hence the first result. The second result is a direct consequence of the first result applied to f * and the equality f * * = f in the case where f is proper, convex and l.s.c.

One of the important characteristics of the notion of the sub-differential [START_REF] Aubin | L'analyse non linéaire et ses motivations économiques[END_REF] is, it's a relationship with optimization. Among the minimization problems, we recall the following results:

Theorem 1.3. Let f : X -→ R,
0 ∈ ∂f (x * ) Proof 0 ∈ ∂f (x * ) ⇐⇒ f (z) ≥ f (x * )+ < 0, z -x >= f (x * ), ∀z ∈ dom(f ) ⇐⇒ f (x * ) = min z∈X f (z).
For the same reason of optimization, we generally introduce their interesting notions to treat some complex problems, like the following :

Definition 1.8. (Moreau envelope)For a proper, lower semi-continuous function f , the Moreau envelope or Moreau-Yoshida regularization of f is given by the infimal convolution:

f γ (x) = inf y∈X { 1 2γ x -y 2 2 + f (y)} , ∀y ∈ X, γ > 0. (1.19)
Definition 1.9. (Proximal operator) For a real-valued convex function f , the proximity operator (proximal mapping) of f is defined as:

prox f (x) = arg min y∈X { 1 2γ x -y 2 2 + f (y)}. (1.20)
This operator has been widely studied for usual convex functions in signal processing.

Chapter 1

Preliminaries and Generality

Functions of bounded variation

In this part, recall the space of functions with bounded variation (BV) and some primary results.

The books [START_REF] Attouch | Variational analysis in Sobolev and BV spaces : applications to PDEs and optimization[END_REF], [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF] propose precise studies of this space. We refer the reader to them for the proof of the different results given below.

The interest of the BV space in image processing is that it allows discontinuous along the curves and surfaces.

Before, recall some notions related to measure theory.

Let Ω be an open subset if X = R n . We denote:

C c (Ω) = {f : Ω -→ R continuous, with compact support in Ω}. C k c (Ω) = {f : Ω -→ R differentiable to order k ∈ N,
L(f ) = Ω f dµ, f ∈ C c (Ω) (1.21)
Proof See [START_REF] Arveson | Notes on measure and integration in locally compact spaces[END_REF], Theorem 3.3 for the complete proof. i.e. µ ∈ M(Ω) ⇐⇒ µ is (or corresponds to) a linear form on C 0 (Ω).

Theorem 1.5. Let L ∈ M(Ω). There exist two positive Radon measure µ + , µ -on Ω such that: µ = µ + -µ -and the corresponding linear form L on C 0 (Ω) is given by:

L(f ) = Ω f dµ + - Ω f dµ -, (1.22)
Put: |µ| = µ + + µ -, then the norm of L on the space M(Ω) is given by:

L M(Ω) = |µ|(Ω) (1.23)
|µ| is said the variation measure of µ.

Chapter 1

Preliminaries and Generality

Proof We refer the reader to see [START_REF] Simonnet | Measures and probabilities[END_REF], Proposition 20.3.2 for detailed proof. Now, to see the relationship of Radon measure with BV space, let us consider:

X an Euclidean space, Ω a sub set open in X and f ∈ L 1 (Ω, R n ) , ϕ ∈ C ∞ c (Ω).
In the sense of distribution D (Ω), ∇f is well defined and we have:

< ∇f, ϕ > D ×C ∞ c (Ω) = - Ω f divϕ(x)dx (1.24)
Hence, the functional:

L : ϕ -→< ∇f, ϕ > verifies: | L(ϕ) |≤ f L 1 • ϕ C 1 c (Ω)
Definition 1.12. Let u ∈ L 1 (Ω). We say that u is a function of bounded variation, if the deviation of u in the sense of the distribution is a finite Radon measure, denoted Du,that is if:

Ω u divϕ(x)dx = - Ω ϕ d(Du) (1.25) Definition 1.13. (Total variation) A function u of L (Ω) has bounded variation if J(u) < ∞,
where:

J(u) = sup{ Ω u(x) divϕ(x)dx/ϕ ∈ C 1 c (Ω), ϕ L ∞ (Ω) ≤ 1}. (1.26)
Definition 1.14. Let u ∈ L 1 (Ω). We say that u is a function of bounded variation in Ω if the distributional derivative of u is representable by a finite Radon measure in Ω,

Ω u ∂φ ∂x i dx = - Ω φ dD i u, ∀φ ∈ C ∞ c (Ω), i = 1, ..., n (1.27)
for some R n valued Radon measure Du = (D 1 u, D 2 u, ..., D n u) in Ω.

The vector space of all functions of bounded variation is denoted BV (Ω).

Remark 1.1. The formula (1.27) remains valid for ϕ ∈ C 1 c (Ω). It can still be written as:

Ω u divφdx = - n i=1 Ω φ i d D i u, ∀φ ∈ C 1 c (Ω), ∀φ ∈ C 1 c (Ω, R n ) (1.28)
In what follows, Ω, a boundary open of R 2 of Lipschitz boundary and

C 1 c (Ω, R 2 ) is the space of functions C 1 with compact support in Ω has values in R 2 .
Definition 1.15. A function u of L 1 (Ω) (with values in R) has Bounded Variation (BV) in Ω if J(u) < +∞ where:

J(u) = sup{ Ω u(x)divφ(x)dx | φ ∈ C 1 c (Ω, R 2 ), φ ∞ ≤ 1}. (1.29)
Chapter 1

Preliminaries and Generality we note,

BV (Ω) = {u ∈ L 1 (Ω) | J(u) < +∞}
the space of such functions.

Theorem 1.6. (Optimality for Convex Minimization Problems)

Let f : X -→ [-∞, +∞] be a convex functional on a real normed space. Then

u * ∈ X solves min u∈X f (u) ⇐⇒ 0 ∈ ∂f (u * ). Proof An element u ∈ X is a solution of min u∈X f (u) if and only if f (u * ) ≤ f (u) for all u ∈ X.But f (u * ) = f (u * )+ < 0, u -u * > is, by the definition of subgradients, equivalent to 0 ∈ ∂f (u * ).
We invite readers to refer to Rockafellar [START_REF] Rockafellar | Extension of Fenchel duality theorem for convex functions[END_REF] for background information found on the theory of Fenchel duality.

Definition 1.16. The Fenchel conjugate f * of a function f : X -→ [-∞, +∞] is defined as:

f * (v) = max u∈X < u, v > -f (u), ∀v ∈ X (1.30)
where < •, • > denotes the inner product defined on X.

This function is also referred to as the convex conjugate or LegendreFenchel transformation of f . Theorem 1.7. (Fenchel-Rockafellar Duality) Let f, g : X -→ [-∞, +∞] be proper, convex and lower semi-continuous on the real Banach space X. Further, let A : X -→ Y be linear and continuous mapping. Suppose that the minimization problem:

min u∈X {f (u) + g(Au)} (1.31)
has a solution u * ∈ X. If there exists some u 0 ∈ X such that f (u 0 ) < ∞, g(u 0 ) < ∞ and g is continuous at (Au 0 ), then

max v∈Y {-f * (-A + v) -g * (v)} = min u∈X {f (u) + g(Au)} (1.32)
In particular, the maximum is realized at some v * ∈ X * .

Recall that, we use A + to denote the adjoint linear operator of A and X * denotes the domain of f * Theorem 1.8. (Fenchel-Rockafellar Optimality System) Let f, g : X -→ [-∞, +∞] be proper, convex and lower semi-continuous functionals, and assume strong duality is fulfilled and that the primal problem has a solution, that means, we have:

max v {-f * (-A + v) -g * (v)} = min u {f (u) + g(Au)} (1.33)
Chapter 1

Preliminaries and Generality then a pair (u * , v * ) is a solution of the primal-dual problem if and only if

-A + v * ∈ ∂f (u * ), v * ∈ ∂g(Au * ) (1.34)

Image characteristics

A digital image is acquired from an imaging device (digital cameras, medical scanner, satellite, etc.)

and represents a certain scene in the continuous universe.

We observe an image on a finite grid Ω included in R d . We have thus had access only to a discretized version of the image. That is if one is working in R 2 .

For classical images, the dimension d of the grid is 2, but the same model can be used for films, in which case d = 3. We then speak of voxels rather than pixels. Our formalism applies to any dimension and thus allows us to deal with possibly larger sizes.

Several characteristics detail the content of an image. An image is represented as a vector in a color space. Most classical representation is done by grayscale in 8-bits: the possible values for the intensity of the image pixels are integers between 0 and 255. It is possible to manage color images in the same way: their most standard presentation is called RGB for Red-Green-Blue, so the intensity is a vector of R 3 , indicating the level of each color component. For example, noise in an image is considered a phenomenon of abrupt pixel intensity variation compared to its neighbors. It comes from the lighting of the optical and electronic devices of the sensor. The dimension is the size of the image, presented in the form of a matrix whose elements are numerical values representative of the light intensities (pixel); the number of rows of this matrix multiplied by the number of columns gives us the total number of pixels in an image. The marked opposition between two regions of an image, more precisely between the dark and light regions of this image, is said to be the image's contrast.

The word resolution is the size of an image (number of pixels), which refers to the total number of horizontally or vertically displayable pixels; where this number is bigger, the resolution is better.

Mathematical representation of an image

A continuous gray-scale image is considered that is a two-dimensional function:

u : R 2 -→ [0, 1] (x, y) -→ u(x, y)
The variable (x, y) is the pixel of the image, and it belongs to N 2 for digital images, and u(x, y)

represents the luminous intensity or the grey level of the image at pixel (x, y).

Grey-scale images are often presented as a matrix of scalars, and color images are represented as vectors.

Chapter 1

Preliminaries and Generality

An image consists of grey values of a rectangular point grid as

{u ij : i = 1, • • •, N and j = 1, ••, M }
with u ij denotes the grey value at pixel (i, j); N and M are the width and height of the image in pixels, respectively.

For a color image, it is not enough to know its gray level for each pixel: it is necessary to know the intensity of each of the three channels of the fundamental colors, the red R, the green G, and the blue B. A color image is a vector function as follows:

u : R 3 -→ [0, 1] 3 (x, y) -→ u(x, y) = (R(x, y); G(x, y); B(x, y))
Thus, even if we assume that there is a true "beautiful" image u, the observer sees only a degraded version u ε ,with ε : R 2 -→ R, the noise added to the image u. The difference between the real image and the perceived image will be called noise (a term that goes back to the genesis of signal processing and thus to sound signals).

Degradation of image

Several types of noise can degrade the quality of the digital image during acquisition and transmission.

Thus it is important to know the nature of the noise that contaminates the image. First the noise may be produced by the sensor electronics, which generate an additive Gaussian noise signal superimposed on the useful signal. Second, it may be due to environmental conditions during acquisition, such as insufficient lighting during shooting and the low light limit of the sensor. A second noise source may be due to a possible excessive temperature of the environment and/or the sensor; electronic noise increases linearly with temperature. Images can also be corrupt during transmission, mainly due to interference in the channel used for transmission. For example, an image transmitted using a wireless network may be damaged due to lighting or other atmospheric disturbances.

the most important forms of noises, these are Gaussian noise, salt and pepper noise, poison noise, impulse noise, and speckle noise. This thesis used two types of noise (Gaussian noise and salt and pepper noise) to affect the initial image. First, for simplicity, we make the even more restrictive assumption that the noise is additive, i.i.d and Gaussian. White Gaussian noise is already relevant in many cases and is standard in the image processing community. The Central Limit Theorem can legitimize this Gaussian model in the case where many minor errors accumulate. It is also a simple and classical model in statistics. Moreover, the acquisition noise can be modelled Gaussian when the light intensity is sufficiently large.

Chapter 1

Preliminaries and Generality

Gaussian noise

Gaussian noise is the type of degradation produced by the sensor electronics (linear relationship with sensor temperature) and is related to the low-light limit of the sensor. It is, therefore, the noise that mainly taints digital images and that we will consider later. To synthetically create this uncorrelated additive Gaussian noise (white), a Gaussian random variable is added to the system at the "ideal" image.

Mathematically the additive Gaussian noise model is then expressed as follows:

u ε (i, j) = u(i, j) + ε(i, j), ∀(i, j) ∈ Ω (1.35)
where, for all (i, j) ∈ Ω , ε(i) is a Gaussian vector, moreover, the vectors ε(i) i∈Ω are independent and identically distributed (i.i.d). We will also assume that we know the intensity of the noise, measured in terms of variance. That is, it is assumed that V ar(ε(i)) = σ 2 .I d where the matrix I d is the identity

matrix of size d × d.
This assumption amounts to supposing that the noise affects the various color components independently and with the same independently and with the same intensity.

The Gaussian distribution or normal distribution in probability theory and statistics is usually denoted by N (µ, σ 2 ), where the parameters µ and σ 2 are the mean and the variance of the distribution, respectively. The probability density function describes the Gaussian distribution.

P (x) = 1 √ 2πσ 2 e -(x-µ) 2 2σ 2
(1.36)

The figure(1.1) below shows an illustration of Gaussian noise, the Gaussian distribution function (probability distribution function) of Gaussian noise. White noise is a realization of a random process in which the power spectral density is the same for all frequencies in the baseband. Additive white Gaussian noise is a white noise that follows a normal distribution of a given mean and variance.

Salt and pepper noise

The saturation of several sensors can create this noise. Some pixels are saturated. This translates visually into white and black dots distributed with a specific density in the image. Salt and pepper noise is also known as bipolar impulse noise and expressed as a percentage of the total number of pixels in the picture, with a value of 0 or 255. The noise generated is evenly distributed throughout the whole image. The corrupted pixels are alternately set to the minimum and highest value, giving Chapter 1

Preliminaries and Generality

Speckle Noise

Speckle noise is multiplicative noise, unlike Gaussian or Salt and Pepper noise; this reduces the quality of the image by giving images a backscattered wave appearance caused by many microscopic dispersed reflections flowing through internal organs and makes it more difficult for the observer to distinguish fine details in the pictures. This noise is found in many systems, including synthetic aperture radar (SAR) images, ultrasound imaging, and many more. Speckle noise can exist similarly in an image as Gaussian noise. Its probability density function follows gamma distribution given as:

P (x) = x k-1 e -x θ (k -1)!θ k (1.38)
The distribution representation of this noise is shown below (1.3):

Poisson Noise

Poisson noise is produced by the image detector and nonlinear recorder responses. The image data determine this type of noise. Because detecting and recording procedures incorporate arbitrary electron emission with a Poisson distribution and a mean response value, this expression is utilized.

Because a Poisson distributions mean and variance are the same, the image-dependent term is considered to have a standard deviation if it is assumed that the noise has one variance.

The number of electrons received x can then be modelled by a Poisson law P λ (x) of parameter λ:

P λ (x) = λ x e -λ x! (1.39)
The plot of probability distribution function of Poisson noise is shown in figure(1.4) bellow:

Denoising of image

Denoising is a vast and delicate subject in image processing. There are many approaches to image denoising from different mathematical theories. The main objective of image denoising is to attenuate Chapter 1

Preliminaries and Generality

• Textures, a visual rendering of a material whose details can be of very fine dimensions.

• The contours, also called edges, are discontinuities of the signal in intensity.

In recent years, digital acquisition technologies have continued to improve. These different developments have allowed the appearance of miniaturized sensors, energy-efficient and especially capable of acquiring data with a quality increasingly fine. However, despite these various improvements, the sensors remain subject to disturbances (electronic noise, magnetic resonance, etc.).

These electronic problems added to the constraints of acquisition (distortions, masking effects, etc.) come to disturb the capture, leading to defects such as noise, parasitic points, holes in the acquired data, etc. These imperfections are inherent to any capture technology and cannot be corrected by hardware. They, therefore, require special software processing. Thus, at the same time as the appearance of the first digital sensors at the beginning of the 90 s, the first restoration methods, such as denoising, reconstruction, oversampling, etc., appeared.

In this section, we will give an overview of some methods for image denoising.

The variational approach

The dominant image reconstruction approach consists of a broad class of methods based on optimization, known as the variational approach [START_REF] Aubert | Mathematical problems in image processing, Partial differential equations and the calculus of variations[END_REF]. An objective function is constructed that penalizes both a data fidelity term and a regularization term. The latter term penalizes implausible images with large values, thus favoring some form of the regularity of the solution. An iterative minimization procedure is applied to find a solution.

The first attempt to recover the target image u from the data f is to calculate

û = arg min u Au -f 2 (1.40)
which corresponds to a least squares solution. This type of method is not appropriate when dealing with problems with much missing information. To overcome this problem, Tikhonov proposed a way to weigh the solution space and, consequently, favor one type of solution over another. Here is a general variational formulation [START_REF] Ekeland | Convex analysis and variational problems[END_REF] for solving linear inverse problems:

û = arg min u Au -f 2 + λϕ(u) (1.41)
This objective function is the sum of a data fidelity term and a regularization term λϕ(u), both positive. For the former, we can choose the Euclidean distance as in (1.40), but other options are also possible. Now, we list many regularization methods [START_REF] Brezis | Opérateurs maximaux monotones et semi groupes de contractions dans les espaces de Hilbert[END_REF] encountered in the literature:

Chapter 1

Preliminaries and Generality

• Total Variation: The regularization term is

ϕ(u) = N i,j=1 ∇ i,j u p , ∇ i,j u = (u i+1,j -u i,j , u i,j+1 -u i,j ) (1.42)
where u i,j is the pixel value of u at position (i, j) and p = 1 or 2

• Sparsity: ϕ(u) = Lu 1 for an appropriate matrix L, so promoting sparsity of the coefficients of the image computed in some representation space.

• Combined regularization: ϕ(u) can be the sum of various regularization terms, combining the goodness of different assumptions into the same solution

1.4.3.

Bayesian approach

Through statistical modeling, we can derive our prior belief about the unknown signal in the form of a probability distribution, called the prior distribution. When the probability distribution is combined with the observed data gives rise to the so-called Bayesian methods. The resulting posterior distribution contains a piece of important information about the unknown signal. In particular, we can compute different point estimates such as maximum a posteriori estimator but also perform uncertainty estimation as confidence intervals.

Principle

In this part, we will focus on the choice of the attachment function to the data h to solve the inverse problem modeled by equation 1.40. One approach to choosing this function is to rely on a Bayesian formulation of the problem.

Suppose that u and f are realizations of random vectors U and F continuous. Using a maximum a posteriori approach, we can define an estimate û of u by :

û ∈ arg max u P U/F =f (u) (1.43)
where,P U/F =f (u) the a posteriori probability density of U such that F = f . Thanks to Bayes' theorem, which allows us to invert the conditional probabilities, we have:

P U/F =f (u) = P (F/U =u) (f )P U (u) P F (f ) (1.44)
and we obtain:

û = arg max u P (F/U = u)P U (u) P F (f ) = arg max u P (F/U = u)P U (u) (1.45) Chapter 1
Preliminaries and Generality which allows us to invert the conditional probabilities, which can be rewritten as the minimizer of a sum of logarithms

û = arg min u {-log P (F/U =u) (f ) a -log P U (u) b } (1.46)
Interpretation The term (a): represents the probability of obtaining f from the image u: its attachment term to the data.

The term (b): represents the probability of obtaining the image u, independently of the observation.

f : it is the regularization term. It is interpreted as a penalty on the image u to verify satisfactory intrinsic properties decided a priori. This term is sometimes called a priori.

Chapter 2

Compressed Sensing

Compressed Sensing Theory (CS)

Since the founding works concerning sampling, it has been known that it is possible to accurately reconstruct a continuous band-limited signal with a countable number of samples (Nyquist, 1928; For the proof of this theorem, we refer the reader to see [START_REF] Shannon | A mathematical theory of communication[END_REF].

Moreover, there is an interpolation formula that is exact if the sampling frequency is at least twice the maximum frequency of the signal. These results have motivated the study and use of digital signals rather than analog ones. However, these data's storage, transfer, and processing remain limited by computer means. Therefore, it has become necessary to compress this information. The compression of information can be interpreted as a change of representation. We hope that the signal can be expressed or approximated with fewer coefficients in this new representation and will therefore be less expensive to store and process.

Compressed sensing theory is related to other fields like optimization theory, random matrix theory, linear algebra, etc. The process of compressed sensing is composed of three essential parts:

• The first concerns signal sparsity, basis representation to obtain a vector with a fewer non zero values;

• The second is the measurement matrix incoherence in obtaining the compressed sensing measurements;

• The last one is the implementation of an efficient reconstruction algorithm's robustness.

Sparsity

Definition 2.1. We say that a signal u ∈ R N is sparse of order K if it contains at most K non-zero coefficients with u 0 ≤ K. Let us note N K , the subspace of R N formed by the K-order sparse signals.

Definition 2.2. Let u ∈ R N such that, the support of a vector u is the set

supp(u) = {i ∈ 1, ..., N /u i = 0}, (2.1)
with, 1 ≤ i ≤ N , u j is the i-th coordinate of u and the size of the support of u is denoted:

u 0 = card(supp(u)) (2.2)
it's the number of non-zero coefficients of u.

where card(.) denotes the cardinality of the set. Mathematically, the l 0 is not a norm, but it is commonly (and abusively) called the l 0 norm. It counts the total number of nonzero elements in a vector.

The three essential parts of compressed sensing as described in figure (2.1): Compressed sensing is based on the fact that most real-life signals, such as images, music, and ultrasonic signals, are linear combinations of only a small number of relevant coefficients from an appropriate base or dictionary. This signal property is called sparsity, and it is prevalent for sparse recovery.

Chapter 2

Compressed Sensing

As we have said, the notion of sparsity is linked to the notion of economy. Thus, a representation of a real vector is said to be parsimonious if there is a more economical way to describe the vector than to give the value of all its elements.

Mathematically, a signal x with a sparse coefficients vector u contain only a few nonzero values (K entries).

Let us assume a real signal x ∈ R N possessing a sparse representation in some transform domain, with K sparsity, K N . The signal x can be represented by a linear combination as:

x = N i=1 u i ψ i (2.3)
where Ψ = (ψ 1 , ψ 2 , .., ψ N ) T is the basis matrix (also called representation basis)and u ∈ R N is a weighted N dimensional vector in the basis matrix Ψ, u i =< x, ψ i >; where, < •, • > denotes the scalar product. Then, for simplicity a signal can be written as follows:

x = Ψu (2.4)
Ψ is a N × N matrix with ψ i as column i.

Other measures describe the sparsity of a vector besides the l 0 norm. Commonly used sparsity measures are described in the following section.

Sparsity Measures

Spasticity is usually defined from the point of view of the l 0 norm as the number of non-zero inputs.

Nevertheless, this definition may not be practical in many real situations, such as in the noisy case.

In [START_REF] Karvanen | Measuring sparseness of noisy signals[END_REF], several alternative sparse measures are described. The kurtosis, log and tanh functions are examples of sparse measures.

Karvanen and Cichocki, in 2003 [START_REF] Karvanen | Measuring sparseness of noisy signals[END_REF] defined the kurtosis κ 4 which measures the peakedness of a distribution and also be used as sparsity measure. It is defined as:

κ 4 (u) = N i=1 u 4 i ( N i=1 u 2 i ) 2 (2.5)
with u is a sparse vector.

Rickard and Fallon, 2004 [START_REF] Rickard | The Gini index of speech[END_REF] measure sparsity using the log measure, which enforces sparsity outside some range. It is defined as:

log(1 + u 2 ) = N i=1 (1 + u 2 i ) (2.6) Chapter 2
Compressed Sensing tanh a,b is defined as follows by Karvanen and Cichocki, 2003 [65]:

tanh a,b (u) = N i=1 tanh(|au i | b ) (2.7)
with a and b are positive constants and the difference between tanh a,b , and the l p norm is that the former is limited to the range [0, 1].

In [122], the authors propose using the Gini index (GI) to measure sparsity, computed as:

GI(u) = 1 -2 N i=1 | u (i) | u 1 ( N -i + 1 2 N ) (2.8)
such that: u is a vector in R N with its element reorder in ascended order as:

| u (1) |≤| u (2) |≤ • • • ≤| u (N ) |
, and u 1 is the l 1 norm of u.

Acquisition

We have seen that a signal's sparsity could help reduce the dimension of a signal while keeping the information necessary for its reconstruction or approximation. This process of dimension reduction in the K-best term approximation is adaptive because the indices of the kept coefficients depend on

x. The main objective of the compressed acquisition is twofold. The first is to reduce the dimension of x in a non-adaptive way. This property is known as the universality of acquisition. The second is that in transformation compression, to retain the best terms, it is necessary to the best terms, it is necessary to calculate all the coefficients of the signal in the new base. Compressed acquisition directly calculates the coefficients necessary for the reconstruction. There is no need to throw away coefficients after their calculation.

Let us assume a real signal x ∈ R N possessing a sparse representation in some transform domain, with K sparsity, K N . Then, the signal can be written as follows:

x = Ψu (2.9)
where Ψ = (ψ 1 , ψ 2 , .., ψ N ) is the basis matrix, u ∈ R N ×1 is the column vector of weighting coefficients with

u i =< x, ψ i >= Ψ T x using Φ ∈ R M ×N sensing matrix the dimension measure reducing into K ∼ M . < •, • > denotes the inner product.
Using Φ as a sensing matrix of dimension M × N with K ∼ M N . The measurement vector is acquired according to the following linear model:

f = Φx = Au (2.10)
with, A = ΦΨ and A is the (M × N ) matrix which must obey two conditions known as incoherence condition and restricted isometry property (Candes, Romberg, and Tao, 2006) [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF].

Incoherence

The incoherence is the low value of mutual coherence between the sensing matrix Φ and the sparse basis Ψ and the mutual coherence measures the largest correlation between any two elements of Φ and Ψ. The incoherence is defined by (Candes and Tao, 2006) [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] as:

µ(Φ, Ψ) = α max 1≤l,j≤N |< Φ l , Ψ j >| (2.11)
where Φ l and Ψ j are the l -th and j -th columns of Φ and Ψ, respectively. α can be take the values N or √ N according to various authors.

Restricted Isometry Property

The Restricted Isometry Property (RIP) condition is studied in a compressed acquisition because it is sufficient to get good reconstruction guarantees in case the measurements are disturbed. The RIP requires that all subsets of K columns taken from the product A are, in fact, almost orthogonal; they can't be precisely orthogonal since we have more columns than rows. It is defined as:

1 -δ K ≤ Au lp u lp ≤ 1 + δ K (2.12)
where δ K is the RIP constant of order K and say matrix A satisfy the restricted isometry condition of order K.

The restricted isometry condition imposes a minimum number of measurements on the measurement matrix. This number depends on the choice of the elements of the matrix. We can quote as examples the following matrices:

• Gaussian random matrix (Candes and Tao, 2006) [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] means that all elements of matrix A are i.i.d. according to a centered normal distribution of variance 1/n : n is of the order of K log( N K ).

• Bernoulli matrix, i.e. when the elements of A are i.i.d according to a Bernoulli law of parameter 1 2 : n is of the order of K log( N K ).

Reconstruction Method

The compressed sensing recovery of x from f means find the sparsest u subject to formula (2.10).

The matrix A = ΦΨ is the size (M × N ). It is therefore of rank less than N and is not invertible;

this implies that there are infinitely many solutions û such that: Aû = Au with û = u. The decoding problem is to determine which of these solutions corresponds to u.

The Moore-Penrose pseudo-inverse [START_REF] Albert | Regression and the Moore-Penrose pseudoinverse[END_REF] is a conventional approach which allows to obtain the solution of minimal l 2 norm, defined as:

û = arg min u u 2 subject to f = Au (2.13)
If the matrix A is of rank M , the solution is given by the equation below,

û = A + f = A T (AA T ) -1 f (2.14)
where A + is the Moore-Penrose pseudo-inverse, and A T is the transpose of A.

Such that the Moore-Penrose pseudo-inverse matrix is defined below:

Definition 2.

(Moore-Penrose inverse)

A matrix A + satisfying the following conditions is called the Moore-Penrose inverse of A:

1. AA + A = A, 2. A + AA + = A + , 3. (A + A) T = A + A, 4. (AA + ) T = AA + .
A + is unique.

Indeed, the minimization of the l 2 norm does not invite the sparsity solution, and the hypothesis of the existence of the signal's sparsity is not taken into account.

In order to exploit this a priori information, we must use an approach that penalizes non-sparse solutions. The natural approach is to try to minimize the norm l 0 of the solution, which results in the following problem of minimization:

û = arg min u∈R N ×1 u 0 subject to f = Au (2.15)
where, u 0 = N i=1 |u i | 0 the cardinal of u. However, searching for such a solution cannot be done in polynomial time because it is an NPcomplete problem and cannot be solved in practice for large enough values of N . Let us define an NP-complete problem below:

Chapter 2 Compressed Sensing Definition 2.4. NP-complete A problem X is said to be NP-complete if: • X ∈ N P , • Y ≤ X, ∀Y ∈ N P .
Y ≤ X means that Y is not much harder than (easier than) X.

To overcome the feasibility problem of finding a solution as defined in problem (2.15), l 0 constraint is replaced by l 1 norm as the following optimization problem [START_REF] El Mahdaoui | Image recovery using total variation minimization on compressive sensing[END_REF]:

û = arg min u∈R N u 1 subject to f = Au (2.16)
• 1 denotes the l 1 norm of a vector,i.e., u 1 = N i=1 |u i |. It is then a question of replacing the l 0 norm with a norm that favors parsimony and makes the problem solvable (polynomial complexity algorithm). The problem l 1 norm is the most used sparsity measure, probably due to its advantages; it can be interpreted as a convex optimization problem and can be solved by linear programming if the coefficients are real [START_REF] Candes | An introduction to compressive sampling[END_REF].

After estimating the sparse coefficient vector, one can recover the estimated image as x = Ψû.

In the framework of the compressed sensing, the measurement f ∈ R M is carried out via a multitude of linear projections of the signal x, element of R N . These projections are non-adaptive but must preserve the structure of the signal. M measurements are performed via a known real sensing matrix φ of size M × N with ideally M N .

This measurement may be pertubed by an additive noise δ. Therefore, this measurement operation is written mathematically as:

f = Φx + δ (2.17)
where δ ∈ R M is an unknown noise such that δ 2 ≤ ε, the minimization problem (2.16) needs to be changed to the following problem known as basis pursuit denoising [START_REF] Candes | An introduction to compressive sampling[END_REF].

û = arg min u∈R N u 1 subject to Au -f 2 ≤ ε (2.18)
where ε bounds the amount of noise in the data.

In the literature, the penalty method is the most famous way to transform constrained optimization (2.18) into an unconstrained version, by which the problem becomes:

min u∈R N 1 2 Au -f 2 +α u 1 (2.19)
Problems (2.18) and (2.19) are equivalent because solving one will determine the parameter in the other such that both provide the same solution [START_REF] Zhu | First-order optimality condition of basis pursuit denoise problem[END_REF]. The signal x can be reconstructed by solving the l 1 norm minimization problem under the condition that x is sufficiently sparse, and the measurement matrix ψ is inconsistent with the orthogonal basis φ.

The objective is to reconstruct the x signal from the f measurement via l 1 minimization. It is an inverse problem.

Several iterative reconstruction algorithms solve the problem of denoising algorithms for CS recovery and obtain high performances for nature images other than l1 minimization, the greedy methods could also handle compressive sensing problems by iteratively computing the support of the signal.Such as matching pursuit algorithm [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF], orthogonal matching pursuit [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] and approximate message passing (AMP) [START_REF] Donoho | Message-passing algorithms for compressed sensing[END_REF], an extension of AMP and denoising based AMP [START_REF] Metzler | From denoising to compressed sensing[END_REF].

Mallat and Zheng proposed the Matching Pursuit algorithm (MP) [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF], which proceeds as follows:

From a residue, which has an initial value of the compressed observation, the algorithm selects the atom of the dictionary that is closest to the observation, then removes the contribution of this atom from the residue; and this is repeated until the residue is sufficiently small relative to the initial observation. The solution is, by construction, parsimonious since after n iterations, the obtained signal has at most n non-zero components.

Orthogonal Matching Pursuit (OMP) [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] is an evolution of Matching Pursuit. The problem with

Matching Pursuit is that it can select the same atom several times. OMP overcomes this problem.

The beginning of the iteration is identical to MP:

First, we select the dictionary atom most correlated with the residue. Then, we ensure that the set of selected atoms is orthogonal to the residue by recalculating the decomposition on the set of atoms selected by the algorithm.

However, l 1 minimization methods usually require fewer measurements and better stability than greedy algorithms. When noise exists, l 1 minimization methods provide a much more stable solution and apply the methods to real-world problems.

Chapter 3

Inverse problems and solving methods

Inverse problem

According to J.B. Keller [START_REF] Keller | Inverse problems[END_REF], two problems are said to be inverses if the formulation of one of them implicates the other. This definition is somewhat arbitrary, making the two problems considered an asymmetrical role. A more operational definition is that an inverse problem determines causes knowing effects. Thus, this problem is the inverse of the direct problem, consisting of deducing the effects, the causes being known.

The prediction of the future state of a physical system, knowing its current state, is a typical example of a direct problem. One can consider various inverse problems: for example, reconstructing the past state of the system knowing its current state (if this system is irreversible), or determining the parameters of the system, knowing (part of) its evolution.

The study of inverse problems is complicated and this is due to the possibility of having several solutions because different causes sometimes lead to the same effects. Another practical difficulty of studying inverse problems is that it often requires a good knowledge of the direct problem, so success in solving an inverse problem is usually based on elements specific to this problem. There are however some techniques that have a wide range of applications.

Success in solving an inverse problem usually depends on elements specific to that problem. However, some techniques have a wide range of applicability, and this course is an introduction to the main ones: the regularization of ill-posed problems and the method of least squares. The most important one is the reformulation of an inverse problem to minimize an error functional between the real measurements and the synthetic measurements (i.e., the solution of the direct problem). It will be convenient to distinguish between linear and nonlinear problems. Let us specify that nonlinearity refers to the inverse problem itself and not to the direct problem (considering the parameters as known).

Among the fields in which inverse problems play an important role, we can mention :

• medical imaging (ultrasound, scanners, X-rays ...) ;

• petroleum engineering (seismic and magnetic prospecting, identification of probabilities in a reservoir ...) ;

• radar (determination of the shape of an obstacle);

• quantum mechanics (determination of the potential) ;

• image-processing (restoration of blurred images). This last area will be our main interest in the following. From a mathematical point of view, these problems fall into two main groups : linear problems (ultrasound, image processing), which are reduced to the resolution of an integral equation of the first kind ; non-linear problems, which are most often questions of parameter estimation in differential or partial differential equations.

Examples of an inverse problem

We are interested in the estimation of parameters in a partial differential equation:

         ∂y ∂t -n i=1 ∂ ∂x i (a ∂y ∂x i ) = f inΩ×]0, T [ y(x, 0) = y 0 (x) inΩ ∂y ∂n = g inΓ×]0, T [ (3.1)
It is the heat equation, y is the temperature, f is a source term, a is the thermal conductivity, and g is the heat flux (incoming or outgoing). We can use the same equation to model a single-phase flow (such as oil): y is the pressure, f represents the pumping wells, a is the permeability of the medium, and g = 0 for the closed space. The problem is the following: from measurements of y at some points and times, we must find a. The direct problem is trivial, but the inverse problem can be more complicated.

In signal processing, from medical imaging to analog and digital conversions, from seismic exploration to high definition video display, recovering high resolution and high-quality signals from partial and noisy measurements are the dream behind inverse problems. The measurements are modeled by a linear operator applied to the input signal, but this operator is generally not invertible. Therefore, computing an accurate estimate of the signal is not possible without strong a priori information about the signal.

In signal and image processing, many applications are related to the resolution of inverse problems.

This resolution estimate an unknown original signal from a degraded (or incomplete) version. Mathematically, the original signal can be represented by a vector u ∈ R N and the degraded observation of this signal by a vector f ∈ R M linked to u by the relation:

f = Au + ε (3.2)
where A ∈ R M ×N is a matrix modelling a degradation (in image processing, this operator can, for example, represent a blur appearing during the capture of the image), and ε ∈ R M is a realization of a random variable ξ, representing noise and/or a model error. The objective is then to produce an estimate û from f and A.

• Existence of Solution: ∀f ∈ Y, ∃u ∈ X Au = f • Unicity: ∀f ∈ Y, ∃!u ∈ X.
• Stability: The solution u depends continuously on the data f .

If at least one of these three conditions is not verified, then the problem is said to be ill-posed.

Ill-posed problems

In the sense that the previous definition, a problem that is not well-posed is said to be ill-posed.

The non-existence and non-uniqueness of the solution to an ill-posed problem are undoubtedly serious difficulties, but they are restored. However, the lack of continuity is more problematic, especially given an approximate or numerical solution. That is to say, it will not be possible (independently of the numerical method) to approach a satisfactory solution to the inverse problem because the available data will be noisy and close but different from the real data.

Exemples

Numerical Derivation The problem of numerical calculation of the derivative of a function we measure is an ill-posed problem. To see this, we consider the Banach space X of C ∞ functions defined on [0, 1] with values in R equipped with the uniform convergence norm. We define the linear operator I : X -→ X by:

If (t) = t 0 f (s)ds, ∀f ∈ X. (3.3)
The immediate problem is, therefore, to compute the primitive of f ∈ X which is 0 in 0. This problem is well-posed. To solve the inverse problem, one must, knowing F ∈ X compute its derivative f ∈ X.

The inverse problem is ill-posed. Indeed, we know there is a unique derivative of F ∈ X, but the derivation operator is not continuous from X into itself.

This instability can be observed in Figure (3.2). We have calculated the derivative of the function t → sin(3t) + p(t) where the function p is a small perturbation (white noise). We observe a significant deviation between the derivative of the function F and that of the function perturbed. We define it in R 2 . Note that in the case of R 2 the Radon transform is also called the "X-ray" transform. We consider a function f : R 2 -→ R has a compact support Ω ⊂ R 2 . For an angle Chapter 3

Inverse problems and solving methods θ ∈ [0, 2π[ and s > 0, we define the Radon transform Rf of f by:

Rf (θ, s) = {x∈R 2 |<x,ν θ >=s} f (x)dx, ν θ = cos θ sin θ (3.4)
This corresponds to integrating the function f on all lines δ defined by the angle θ and the number s > 0 (see Figure 3.3).

The inverse problem is determining the function f from its Radon transform Rf . Radon [START_REF] Radon | On the determination of functions from their integral values along certain manifolds[END_REF], in 1917, was the first to give an explicit inversion formula. This inversion formula is written as:

f (x) = 2π 0 R + ∂ s Rf (θ, s) < x, ν θ > R 2 -s dsdθ (3.5)
Where the integral on R + is defined in the sense of principal values, and ∂ s is the derivation operator to the variable s. The problem is badly posed because the inversion formula provides an unstable solution. As in the previous example, the derivation operator used in the inversion formula is the source of this instability. For complete details on the Radon transform, we refer to [START_REF] Natterer | The mathematics of computerized tomography[END_REF].

Inverse problem solving methods

We recall that in the case of linear inverse problems, the discrete problem is written as:

Au = f (3.6)
where, A is a rectangular matrix.

If we want to solve the equation Au = f , we have two cases:

• If A is square and invertible, the solution is unique and simple; it is expressed by u = A -1 f .

• If A is not invertible and we suppose that the problem is ill-posed, several methods of solution exist.

Regularization methods

In this part of the chapter, we discuss the principle of some regularization methods. This part is an introduction to the most common regularization methods: Tikhonov's method and Morozov's method, and we end with an iterative method that is Landweber's method. For a more detailed reading, see the book by Kirsch [START_REF] Kirsch | An introduction to the mathematical theory of inverse problems[END_REF].

To regularize an ill-posed problem is to replace it with another one, well-posed so that the error committed is compensated by the gain in stability.

The main difficulty in applying an accrual method is determining the accrual parameter itself.

In this part, we will consider A : X -→ Y , a bounded linear operator defined between two Hilbert spaces X and Y .

Tikhonov method

The principle of Tikhonov's method for solving the inverse problem posed by Au = f is to choose as a solution the element u µ ∈ X, which minimizes the quantity

Au -f 2 Y + µ u 2 X (3.7)
The existence and uniqueness of the minimum are immediate by coercivity and strict convexity of u -→ u 2 X . The parameter µ is called the regularization parameter. This parameter must be chosen small enough for the element u µ that realizes the minimum to have a small error with the data f . It must also be large enough so that the strict convexity of the term u X corrects the instability of the problem. Theorem 3.1. Let: µ > 0 and A : X -→ Y a bounded linear operator from the Hilbert space X to the Hilbert space Y . Then the Tikhonov functional admits a unique minimum in u µ ∈ X. The element u µ is the solution of the normal equation

αu µ + A + Au µ = A + f (3.8)
For the proof of this theorem, see [START_REF] Kirsch | An introduction to the mathematical theory of inverse problems[END_REF] page37.

The solution of equation (3.8) can be written as u µ = R µ f , such that:

R µ := (µI + A + A) -1 A + : Y -→ X. (3.9)
By choosing a singular system (α j , u j , f j ) for the compact operator A, we see that R µ f admits the following representation:

R µ f = ∞ n=0 α j µ + α j 2 (f, f j )u j = ∞ n=0 q(µ, α j ) α j (f, f j )u j (3.10)
where, q(µ, α) = α 2 µ+α 2 . This function q is called the filter function. It remains to show that this operator is indeed a regularization operator and under which conditions the choice of µ as a function of the noise level is admissible. It is the object of the following theorem. Theorem 3.2. Let A : X -→ Y be a compact linear operator and µ > 0.

The operator (µI + A + A) is inversible and the operator R µ defined in (3.9)is a regularization strategy with R µ ≤ 1 2 √ µ .It is called the method of regularization of Tikhonov. R µ f is determined as the unique solution u µ, of the equation of the second degree.

µu µ, + A + Au µ, = A + f , (3.11) All choice of µ( ) -→ 0 ( -→ 0) with 2 µ( ) -→ 0 ( -→ 0) is admissible.
For the proof of this theorem and the results on the speed of convergence when -→ 0, see [START_REF] Tikhonov | Numerical methods for the solution of ill-posed problems[END_REF].

The Morozov's principle

We give here an example of a posteriori method of choice of the regularization parameter. We expose the most classical of these, "the discrepancy principle" of Morozov [START_REF] Morozov | Methods for solving incorrectly posed problems[END_REF], or the Morozov's discrepancy principle. According to Kirsch [START_REF] Kirsch | An introduction to the mathematical theory of inverse problems[END_REF], we present a principle based on the Tikhonov regularization method.

We assume that, A : X -→ Y is a compact injective operator defined on two Hilbert spaces X, Y with Im(A) ⊂ Y .

Again we study the equation Au = f, f ∈ Y . We now compute the regularization parameter µ = µ( ) > 0, such that the corresponding Tikhonov solution is the solution of the equation (3.11) and it is the minimum of (3.7)which satisfies the equation

Au µ, -f = , (3.12)
note that the choice of µ by "the discrepancy principle" guarantees, on the one hand, that the error is , on the other hand, µ very small. Definition 3.2. A family of bounded linear operators R µ : Y -→ X is a "regularization strategy" if:

∀u ∈ X, lim µ-→0 R µ Au = u (3.13)
i.e., the operator R µ A simply converges to the identity.

Definition 3.3. A regularization strategy -→ µ( ) is admissible if for any u ∈ X lim -→0
µ( ) = 0 and lim

-→0 sup f ∈Y { R µ( ) f -u X ; such that Au -f Y ≤ } = 0 (3.14)
The initial data f ∈ Y is never known exactly: some noise always disturbs it. Let us note f the perturbed data where the number > 0 is the noise level, i.e.

f -f Y ≤
Let u µ, := R µ f be the approximation of the solution of the inverse problem Au = f obtained with the regularization operator and the perturbed data.

Theorem 3.3. Let A : X -→ Y linear compact operator with dense image on Y .

Let Au = f, u ∈ X, f ∈ Y, f ∈ Y such that: f -f Y ≤ < f . Let u µ( ) the Tikhonov solution which satisfies
Au µ( ), -f = , for all ∈ (0, 0 ).

Then:

1. u µ( ), -→ u when -→ 0. Then "the discrepancy principle" is admissible.

Let

u = A + f ∈ A + (Y ) with f ≤ K then, u µ( ), -u ≤ 2 √ K
For this purpose, "the discrepancy principle" is an optimal regularization strategy under the condition below:

(A + ) -1 u ≤ K.
The proof of this theorem is in [START_REF] Kirsch | An introduction to the mathematical theory of inverse problems[END_REF], page 48.

The determination of µ( ) is equivalent to the problem of finding the root of the monotonic function:

φ(µ) = Au µ( ), -f 2 -2
for fixed .

It is not necessary to satisfy the equation Au µ, -f = , exactly, the inclusion of form

c 1 ≤ Au µ( ), -f ≤ c 2
is sufficient to prove the assertions of the previous theorem.

The following theorem proves that the convergence order O( √ ) is better for the Morozov shift principe .

It's easy to see the explicit form u m = R m y, such that, the operator R m : Y → X is defined as:

R m = a m-1 k=0 (I -aA + A) k A + , m = 1, 2, .. ( 3
u 0, = 0 and u m, = (I -aA + A)u m-1, + aA + f for m = 1, 2, .. all strategy m( ) -→ ∞( -→ 0) with 2 m( ) -→ 0( -→ 0) is admissible. 2. Let u = A + f ∈ Im(A + ) with f ≤ K and 0 < c 1 < c 2 , for each choice m( ) with c 1 K ≤ m( ) ≤ c 2 K
, the following estimate is verified:

u m, -u ≤ c 3 √ K
For c 3 , which depends on c 1 , c 2 , and a.

Then the Landweber iteration is optimal for (A + ) -1 u ≤ K.

Now, let u

= A + Af ∈ Im(A + A), f ≤ K and 0 < c 1 < c 2 , for each choice m( ) with c 1 ( K ) 2 3 ≤ m( ) ≤ c 2 ( K ) 2 3
, we have:

u m, -u ≤ c 3 K 1 3 2 3
For c 3 , which depends on c 1 , c 2 , and a, also the Landweber iteration is optimal for

(A + A) -1 u ≤ K.
For this method, we notice that a solution is more accurate when the number of iterations m is large, but the stability forces us to keep m as small as possible.

Non linear inverse problem

In the case of nonlinear inverse problems, the problem: u = F (x), where F is the operator (or the set of operations) to be performed on x to obtain u, is transformed into an optimization problem that consists in searching for the minimum of the function defined by:

J(x) = 1 2 u -f 2 (3.19)
then we apply one of the optimization methods to solve this type of problem. Difficulties arise at different levels of existence, uniqueness, and stability. So we can use a classical technique to find the solution by adding a constraint or constraints to the cost function and this by using Lagrange multipliers.

Definition 3.4. The Lagrange multiplier is a method to find the stationary points (maximum, minimum) of a derivable function of one or more variables.

This technique helps solve optimization problems under linear constraints.

Solver of inverse problem

As indicated in the above section above, the solution by the pseudo-inverse is an ill-posed problem because the uniqueness of the solution is not guaranteed. Additional information must be incorporated into the image to be restored into the model to ensure the stability of estimate u.

During the last ten years, the methods developed to solve form problem (2.18) have seen a significant interest in connection with the robustness to outliers and specific optimality of this formulation l 1 .

The classical optimization algorithms, such as interior-point methods, were quickly discarded because they were not adapted to the large dimensions of image processing applications.

Augmented Lagrangian method (ALM/MM)

A fundamental class of methods for constrained optimization is to seek the minimizer or maximizer by solving a sequence of unconstrained subproblems iteratively. The solutions to subproblems should converge to a minimizer or maximizer eventually. In 1943, Courant [START_REF] Courant | Variational methods for the solution of problems with equilibrium and vibration[END_REF] proposed the quadratic penalty method, which could be viewed as the precursor to the augmented Lagrangian method. This method penalizes equality constraint violation by adding a multiple of the square of the constraint violation into the objective function and turns the constrained optimization problems to be unconstrained. This approach has been used and studied comprehensively due to its simplicity and intuitive appeal. However, it requires the penalty parameter to go to infinity to guarantee convergence, which May cause a deterioration in the numerical conditioning of the method. In 1969, Hestenes [START_REF] Hestenes | Multiplier and gradient methods[END_REF] and

Powell [START_REF] Powell | A Method for Nonlinear Constraints in Minimization Problems, Optimization[END_REF] independently proposed the augmented Lagrangian method, which, by introducing and adjusting Lagrangian multiplier estimates, no longer requires the penalty parameter to go to infinity for the process to converge.

Let us consider the following equality-constrained convex optimization problem:

min u f (u) subject to Hu = g (3.20) Here, u ∈ R N , g ∈ R M , H is a matrix of dimension M × N and f ∈ R N is a convex function.
The augmented Lagrangian method, originally known as the multipliers method (MM) [START_REF] Hestenes | Multiplier and gradient methods[END_REF], combines the Lagrangian function and a quadratic penalty term. It is applied to solve a constrained optimization problem iteratively, applied to equation (3.20), the augmented Lagrangian can be expressed as follows:

L(u, λ, ρ) = f (u) -λ T (Hu -g) + ρ 2 Hu -g 2 2 . (3.21)
where, λ ∈ R M is a vector of the Lagrange multiplier and ρ > 0 is the augmented Lagrangian penalty parameter for the quadratic infeasibility term [START_REF] Nocedal | Numerical Optimization[END_REF]. We recall in Algorithm 1 the so-called Augmented Lagrangian method: this consists of searching for the optimal solution giving the signal u by alternately updating u and then the Lagrange multiplier λ. The alternate iteration method keeps one vector fixed and updates the other successfully until the stopping criterion is satisfied.

These different steps are summarized in the algorithm.

Algorithm 1 Augmented Lagrangian method Fixed: k = 0, ρ > 0, λ 0 , u 0 Iteration: for k = k + 1 until the stopping criterion is satisfied, repeat steps (1 -2)

Step 1 Keep λ fixed and update u:

u k+1 = min u (f (u) -λ T k (Hu -g) + ρ 2 Hu -g 2 2 )
Step 2 Keep u fixed and update λ:

λ k+1 = λ k + ρ(Hu k+1 -g)
Output: u the optimal solution. Remark 3.1. A stopping condition must be defined whatever the iterative optimization algorithm is considered. It can be determined as a function of the total number of iterations, or function of the error between the restored and original images (when the latter is known) or the convergence of any criterion (PSNR, SSIM, etc.), or as a function of the term MSE of the restored image.

Alternatif direction method of multipliers (ADMM)

Extending the classic augmented Lagrangian method as described in section above,the common advantage of both methods includes the capability of handling the non-differentiability and side constraints. Rather than attacking the problem (3.20) head-on, the ADMM methods Eckstein and Bertsekas use a strategy of alternating minimizations with respect to u and v, keeping the other variable fixed each time. This technique is known in the literature as "non-linear block-Gauss-Seidel". Eckstein and Bertsekas [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF] demonstrated proof of convergence of this method. ADMM algorithm solves problems in the following form:

min u,v f (u) + g(v) subject to Au = v (3.22)
The augmented Lagrangian of (3.22) is in the following form:

L(u, v, λ) = f (u) + g ( v) -λ T (Au -v) + β 2 Au -v 2 2 (3.23)
The steps of the ADMM method are described in algorithm 2.

Algorithm 2 A Summary of direction alternatif method (ADMM)

Initialization: k = 0, β, λ 0 , v 0 repeat minimization of AL with respect to u u (k+1) ∈ arg min L(u, v k , λ k ) minimization of AL with respect to v v (k+1) ∈ arg min L(u k+1 , v, λ k ) Update multipliers λ k+1 = λ k + β(Au k+1 -v k+1 ) k = k + 1 until maximum iteration number is reached.
Thus, the ADMM method inherits the algorithmic framework of ALM methods with the advantage of minimizing u and v sequentially via the solution of two subproblems [START_REF] Carpentier | Décomposition-coordination en optimisation déterministe et stochastique[END_REF] of different complexity. Minimizing u and v sequentially via the resolution of two sub-problems of lower complexity. In this sense, the ADMM approach is a variant of ALM techniques particularly adapted to the solution of structured problems like (3.22). Compared to classical ALM approaches, it is more adapted to the case of non-quadratic, non-smooth functionals whose gradient is not Lipchitz differentiable. Finally, this method allows internal approximations without affecting the stability of the convergence [START_REF] He | New inexact alternating directions method for monotone variational inequalities[END_REF], which makes it a powerful and robust tool compared to classical techniques.

ADMM implementation

The ADMM formalism is adapted to our problem with a new formulation. The steps described in the previous section can be efficiently implemented using point-to-point operations and proximal operators with analytical solutions.

Chapter 3

Inverse problems and solving methods

Solving step 1: diagonalizable least squares

Step 1 of the ADMM algorithm is written:

u k+1 ∈ arg min u∈R n L(u, v k , λ k ) (3.24) = arg min u∈R n < λ k , Au -v k > + β 2 Au -v k 2 2 (3.25)
u k+1 is thus the solution of a linear system in the least-squares sense which is written, omitting the independent terms of u

u k+1 = arg min u∈R n (λ k ) T Au + β 2 (u T A T Au -2(v k ) T Au). (3.26)
The normal equation that the functional of this problem must satisfy is:

A T λ k + β 2 (2A T Au k+1 -2A T v k ) = 0 (3.27)
and leads to the characteristic equation

A T Au k+1 = A T (v k - λ k β ) (3.28)
The solution is then given by:

u k+1 = (A T A) -1 A T (v k - λ k β ) (3.29)

Solving step 2: Proximonial operator

The second step of ADMM algorithm is given:

v k+1 ∈ arg min v L(u k+1 , v, λ k ) (3.30) = arg min v g(v)+ < λ k , Au k+1 -v > + β 2 Au k+1 -v 2 2 (3.31) = arg min v g(v) + β 2 v -(Au k+1 + λ k β ) 2 2 (3.32)
To solve this last equation, we can use proximal operators [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF], sometimes called Moreau operators [START_REF] Moreau | Proximité et dualité dans un espace Hilbertien[END_REF], which have been popularized by works such as those of Combettes and Wajs [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. The proximal operator associated with the convex function f is defined for

x 0 ∈ R n as prox f (x 0 ) = arg min x∈R n f (x) + 1 2 x -x 0 2 2 (3.33)
These proximal operators are potent tools for solving optimization problems (3.32) because, in some cases, the solution is unique and can be obtained analytically [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. Equation (3.32) can then be solved as:

v k+1 = prox(Au k+1 + λ k β ) (3.34)
Chapter 4

Image denoising using a compressive sensing approach based on regularization constraints

Technical framework for image denoising

Several methods exist for image denoising, including total variation image regularization [START_REF] Blomgren | Color tv: Total variation methods for restoration of vector-valued images[END_REF], wavelet thresholding [START_REF] Mallat | A wavelet tour of signal processing[END_REF], [START_REF] Mallat | A wavelet tour of signal processing: The sparse way[END_REF], non-local means [START_REF] Buades | Nonlocal image and movie denoising[END_REF], basis pursuit denoising [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], block-matching, and 3D

filtering [?], among others. Moreover, these methods can perform, to a certain extent, image smoothing and preserve edges. Therefore, we can reconstruct the image well by CS theory to obtain more precise measurements of the original images than the corresponding noisy images.

In particular, we can get them into two classes, as exceptional cases of the proposed framework.

Regularization Functions

Image recovery in these application domains can be formulated as a linear inverse problem, which can be modeled as follows:

f = Au + ε (4.1)
where f ∈ R M is the noisy image observation, u ∈ R N is the cleanest image unknown, ε is an additive noise and A ∈ R M ×N is a linear operator. Given A, image reconstruction extracts û from f , making classical least-squares approximation alone unsuitable. To stabilize recovery, regularization techniques are frequently used, producing a general reconstruction model of the following form:

u * = arg min u∈R N ( 1 2 Au -f 2 2 + λφ reg (u)) (4.2)
where λ > 0 is a regularization parameter, and • 2 denotes the l 2 norm. The fidelity term, Au-f 2 2 , forces the reconstructed image close to the original image, and the regularization function, φ reg (u), performs noise reduction.

The choice of the regularization function is very important for reconstructing an image that reflects, as accurately as possible, the original image of interest. In our work, we have combined total variation and nonlocal self-similarities. The interest of such a choice is linked to the fact that the total variation (TV) model shows high efficiency in preserving the contours and recovering smooth regions. However, this operator is local.It,therefore, does not take into account nonlocal features of the data, such as repetitive structures (such as texture, for example). However, nonlocal selfsimilarity describes the repetitiveness of textures [START_REF] Ouahabi | Multifractal analysis for texture characterization: A new approach based on DWT[END_REF], [START_REF] Djeddi | Discrete wavelet for multifractal texture classification: Application to medical ultra sound imaging[END_REF] or embodied structures in realistic images, which allows the preservation of sharp edges.

Total Variation

The history of total variation is rich and long. In 1881, the total variation was introduced firstly by Jordan for real-valued functions while studying the convergence of the Fourier series.

Decades after research, total variation is widely used to compute discontinuous solutions to inverse problems after being thoroughly investigated. In 1992, Rudin, Osher and Fatemi [START_REF] Rudin | Non linear total variation based noise removal algorithms[END_REF] introduced first the total variation as a regularizing criterion for solving inverse problems [START_REF] Acar | Analysis of bounded variation penalty methods for ill-posed problems[END_REF]- [START_REF] Teboul | Variational approach for edge-preserving regularization using coupled PDE's[END_REF]. Its major interest lies in the fact that it preserves the discontinuities, which can be seen as contours, and edges present in natural images. For this characteristic, TV regularization is widely used in many image processing applications, such as blind deconvolution, inpainting and super-resolution.

Geman and Yang [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] proposed a joint minimization method to solve half-quadratic models. Grounded proposed a fast half quadratic method to solve deconvolution and denoising problems with TV regularization [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF] in their work. Further extended this method to image reconstruction [START_REF] Yang | A fast TVL1-L2 minimization algorithm for signal reconstruction from partial fourier data[END_REF] and multichannel image deconvolution problems [START_REF] Yang | An efficient TVL1 algorithm for deblurring of multichannel images corrupted by impulsive noise[END_REF]. The two central ideas in this approach are splitting and alternating. Total variation (abbreviated TV) regularization can be a generalized l 1 regularization in compressive sensing problems. Instead of assuming the signal is sparse, the premise of TV regularization is that the gradient of the underlying signal or image is sparse. In other words, total variation measures the discontinuities and seeks the solution with the sparsest gradient.

TV minimization has attracted more research activities in the broad area of compressive sensing. TV regularization instead of the l 1 term makes the reconstructed images sharper by accurately preserving the edges or boundaries. In most cases, the edges of the underlying image are essential to characterize different properties than the smooth part. For example, in seismic imaging, detecting boundaries of distinct media plays a crucial role in identifying the geological structure. This advantage of TV minimization stems from the property that it can recover sparse signals or images and dense staircase signals or piecewise constant images.

Definition 4.1. The total variation of u denoted T V (u) is defined as the supremum of absolute differences for any finite partition as:

T V (u) = sup P n i=1 | u(x i ) -u(x i-1 ) | (4.3)
where u is a real valued function defined on [a, b] and P = (a =

x 0 < x 1 < •• < x i < •• < x n = b) of [a, b].

Total variation function

The following result shows the relation between the total variation of the function u in the sense of measures and distributions. The material in this part is mostly extracted from the textbooks [START_REF] Ziemer | Weakly differentiable functions[END_REF], [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF] which we invite the reader to consult for further details.

Definition 4.2. The space BV (Ω) of functions with bounded variation is the set of functions u ∈ L 1 (Ω) such that T V (u) < ∞, endowed with the standard norm

BV (Ω) → R + u → u BV = u L 1 +T V (u, Ω).
The total variation of an image is defined by duality. (Ω) such that, the following quantity, named Total Variation of u, is finite:

T V (u) = sup{ Ω u(x)divφ(x)dx : φ ∈ C 1 c (Ω, R 2 ), φ ∞ ≤ 1, ∀x ∈ Ω} (4.4)
A function is said to have Bounded Variation whenever T V (u) < +∞.

When u ∈ C 1 (Ω) (i.e.u is continuously differentiable over Ω), we have:

T V (u) = Ω |∇u(x)dx| (4.5)
When u is smooth, Du(x) = ∇u(x)dx, such that this measure Du is the distributional gradient of u.

Main properties of the total variation

Lower semi-continuity Definition (4.4) has some advantages. It can be introduced for any locally integrable function (without requiring regularity or derivability). But also,T V (u) is written as a sup of linear forms

L φ : u -→ Ω u(x)divφ(x)dx (4.6)
which are continuous with respect to very weak topologies [START_REF] Ziemer | Weakly differentiable functions[END_REF].

Wavelet Denoising

The introduction of wavelets by Meyer at the beginning of the 90s [START_REF] Taffard | L'analyse par ondelettes[END_REF] allowed the introduction of apriori in image denoising. Wavelets are linear operators that will enable the signal to be split into several scales. A scale defines the order of size of the objects present in the image.

These operators have been applied to image analysis, notably by Mallat [START_REF] Mallat | A wavelet tour of signal processing[END_REF] and Daubechies [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF].

The fields of application are wide: information compression, recognition and detection of objects and image denoising.

Wavelets were originally applied to a very large class of signals: the apriori on the represented objects was weak. Later, Meyer and Mallat's theory extended to a larger class of operators related to more specific cases. We thus observe in the 90 s a race to wavelet type operators.

We denote ψ the wavelet transform function. ψu is the wavelet representation of u: it is a set of images of the same dimension (2D or 3D) as u. In However, the size can differ from one image to another. We associate a scale to an image (the converse is not true). The application of a wavelet operator to an image is called a representation.

The effectiveness of wavelets in image analysis stems from an intrinsic property:

the specificity of representing the data in a parsimonious way, i.e. to code them in the form of a set of coefficients with a large number of zero intensities. This characteristic applies to a large class of potentially represented objects. In other words, the image does not need to be constrained by strong apriori for the information to be compressible.

The basic wavelet shrinkage denoising algorithm comprises three steps:

• Discrete wavelet transform (DWT) [START_REF] Donoho | Wavelet shrinkage[END_REF];

• Denoising [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF];

• Inverse DWT.

Wavelet image denoising, which started with Mallat [START_REF] Mallat | Singularity detection and processing with wavelets[END_REF] , really took off with the work of Donoho and Johnstone [START_REF] Donoho | De-noising by soft-thresholding. Information Theory[END_REF]. The latter has proposed thresholding the coefficients of a wavelet transform applied to an image. They are returned to the reconstructed image using an inverse transform.

The following is the measurement model:

f = u + ε (4.16)
where u is the original image of size M × N corrupted by additive noise ε. The goal is to estimate denoised image û from noisy observation f .

The elimination of this additive noise permits the assumption that the appropriate decomposition basis allows the discrimination of a useful signal (image) from noise. This hypothesis justifies, in part, the traditional use of denoising by thresholding.

The soft-thresholding function (also called the shrinkage function) is as follows:

d û j (k) =          d y j (k) -S, if d y j (k) > S. d y j (k) + S, if d y j (k) < -S. 0, otherwise (4.17)
The hard-thresholding function, another popular alternative, is as follows:

d û j (k) =    d y j (k), if |d y j (k)| ≥ S. 0, otherwise (4.18)
d y j (k) includes the wavelet coefficients of the measured signal at level j, and d û j (k) is the estimation of the wavelet coefficients of the useful signal with the threshold S. Donoho and Johnstone determined the optimal parameter for soft thresholding to minimize the "worst-case error", S = σ √ 2N log N , where N and σ represent the number of pixels for the test image and the standard noise deviation.

A Note that this threshold is only relevant in the case of a Gaussian type deterioration.

Our efficient image denoising scheme

The main contribution of this thesis is to propose a new minimization scheme for the total variation To solve the problem (4.20), we define the augmented Lagrangian function as follows:

L A (w, u, x) = 1 2 Au -f 2 2 + τ w 1 + δ Θ u 1 -γ T (Du -w) -ϕ T (u -x) + µ 2 Du -w 2 2 + β 2 u -x 2 2 (4.21)
where µ and β are the penalty parameters corresponding to Du -w 2 2 and u -x 2 2 , respectively. The problem (4.19) can be rewritten as:

min u F (Au) + R(Gu) (4.22)
Where:

F (Au) = 1 2 Au -f 2 2 , R(Gu) = τ w 1 + δ Θ u 1 = R 1 (Du) + R 2 (u). Such that G ∈ R M ×N , F : R N -→ R , R : R M -→ R.
We make the following assumptions for the fidelity term F (Au) and regularization term R(Gu) defined on the whole space.

• (H1): N ull(G) ∩ N ull(A) = {0}, such that: N ull(•) is the null space of (•).

• (H2): dom(R • G) ∩ dom(F • A) = ∅. The domain of F defined as: dom(F ) = {v ∈ V, F (v) < +∞}.
• (H4): F (v) is continuous over domain of F .

• (H3): F (v) is convex, proper and coercive function.

• (H5): R(.) is closed, proper and convex function.

Under these above assumptions (H1-H5), The Weierstrass theorem [START_REF] Ekeland | Convex analysis and variational problems[END_REF] implies a minimum reached.

According to generalized Weierstrass theorem and Fermats rule [START_REF] Ekeland | Convex analysis and variational problems[END_REF], we have the following results: Theorem 4.1. They exist at least one solution ũ of problem (4.22) which satisfy:

0 ∈ A + ∂F (Au) -div∂R(Gu) (4.23)
where, A + is the adjoint

L 2 operator of A, ∂F (Au) is sub-differential of F at Au and ∂R(Gu) is sub-differential of R at Gu. The minimizer is unique if A • F (u) is strictly convex.
Finding an optimal solution to problem (4. [START_REF] Brezis | Functional analysis, sobolev spaces and partial differential equations[END_REF]) is equivalent to finding a saddle point for L A .

In other words, ( w, ũ, x) is a primal optimal solution, and (γ, φ) is a dual optimal solution if and only if:

L A ( w, ũ, x, γ, ϕ) ≤ L A ( w, ũ, x, γ, φ) ≤ L A (w, u, x, γ, φ) (4.24)
Therefore:

max γ,ϕ F (Aũ) + R(Gũ)+ < γ, Dũ -w > + < ϕ, ũ -x >= L A (ũ, w, x, γ, φ) = min w,u,x F (Au) + R(Gu)+ < γ, Du -w > + < φ, u -x > (4.25)
From inequality (4.19), we can directly read the Kuhn-Tucker optimality conditions: The simulations used two types of noise: additive white Gaussian noise(AWGN) measured by its where p = (p 1 , p 2 , ..., p n ) ∈ R n pixels, ν max is the pixel intensity of salt pixels and ν min is the pixel intensity of pepper pixels.

Dũ -w = 0 (4.26) ũ -x = 0 (4.27) γ ∈ ∂R 1 (Dũ) (4.28) φ ∈ ∂R 2 (ũ) (4.
The sum, q = q 1 + q 2 , is the level of the salt and pepper noise. The Barbara image is relatively complex given its rich texture and its geometric structure: it is clear that our algorithm (see the first line of Figure (4.3)) has this ability to effectively denoise while preserving the details and texture of this particular image. Secondly, we handled the proposed ap- Chapter 4

Visual Quality Comparison

Image denoising using a compressive sensing approach based on regularization constraints Let us recall that a higher PSNR indicates superior image quality and good performance of the algorithm. The values of PSNR illustrate that DCSR yields a higher PSNR than the other methods. Furthermore, to evaluate the effect of denoising our algorithm, DCSR, the gray scale image is corrupted by salt and pepper noise with different levels.

Chapter 4

Image denoising using a compressive sensing approach based on regularization constraints The performance of our method is confirmed in Figure The next challenge is to apply this DCSR algorithm to color or multidimensional images. This is essential since most digital images used in the modern world are not grayscale but usually operate in In the second experiment, we added salt and pepper noise with different noise levels at 10%, 15%, 20%

and 30% to the Lena test image. Then, we applied our denoising algorithm to restore the noisy images and compared them with two other algorithms: Wavelet denoising and GSR algorithm. From these results of PSNR and SSIM, in all the cases, the proposed method achieves the highest scores, which fully demonstrates that the restoration results by the proposed method are the best both objectively and visually.

Algorithm robustness

The robustness of the proposed algorithm will be confirmed in this subsection. The test images are corrupted, on the one hand, by additive white Gaussian noise and, on the other hand, by salt and pepper noise at different noise levels.

Figures (4.12-4.15) show PSNR values (after reconstruction) as a function of the number of iterations for the greyscale images of Barbara and Cameraman and for the color image of Lena using different algorithms.

Note that we did not use the NESTA algorithm in the case of the color images because NESTA is not suitable for color images.

Computational Complexity

In the following paragraph, we will estimate the computational complexity of the proposed DCSR algorithm. It is clear that the main complexity of the proposed algorithm comes from total variation (TV) and the high cost of the nonlocal self-similarities (NLS).

Knowing that the computational complexity of TV is O(N ) [START_REF] Daei | Sample complexity of total variation minimization[END_REF], let us compute that of NLS.

For an image u of N pixels, the average time to compute similar blocks for each reference block is hence,the resulting complexity is O(N (n m 2 + T s )), similarly to the computational complexity of a group-based sparse representation GSR [START_REF] Zhang | Group-based sparse representation for image restoration[END_REF].

We point to the finding that the total computational complexity of our algorithm is O(N (n m 2 + T s ) + N ).

It is interesting to compare the computational complexity of DCSR with competing methods.

The computational complexity of NESTA algorithm is O(N + N log 2 N ) [START_REF] Becker | NESTA: A fast and accurate first-order method for sparse recovery[END_REF] and that of wavelet denoising is O(N log 2 N ) [START_REF] Srivastava | A New wavelet denoising method for selecting decomposition levels and noise thresholds[END_REF].

Table (4.2) summarizes the computational complexity of the four algorithms used.

The result of this comparison is as follows. Relation (4.52) clearly shows that our proposed algorithm is more expensive in terms of computational complexity by an order N . Such an increase is not excessive in view of the good performance of this algorithm and the current computational means,which permit real-time processing. GSR O(N (n m 2 + T s )) [START_REF] Zhang | Group-based sparse representation for image restoration[END_REF] Our algorithm DCSR O(N (n m 2 + T s ) + N )

Discussion

In the real world, the images acquired or measured or recorded have generally suffered degradation of various origins:

• Bad weather conditions (wind, fog, haze, etc.);

• Chain of acquisition of the image;

Conclusion and Outlooks

In recent years, compressive sensing has been intensively investigated from different perspectives and applied to applications in diverse areas.

This thesis is centred around developing an efficient algorithm mainly based on regularization functions (TV and NLS) for l 1 minimizations and extending the idea of compressive sensing to the field of image processing which traditionally requires tremendous computation and resources. We proposed in this thesis an original image denoising method based on compressed sensing that we called denoising compressed sensing by regularizations terms (DCSR), by incorporating two regularization constraints in the model: total variation and nonlocal self-similarity. The optimization of this method is performed by the augmented Lagrangian, which avoids the difficult problem of nonlinearity and non-differentiability of the regularization terms.

The effectiveness of our approach was validated using images corrupted by white Gaussian noise and impulsive salt and pepper noise. Comparing DCSR in terms of PSNR and SSIM to state-of-the-art methods such as Nesterovs algorithm, group-based sparse representation and wavelet-based methods, it turns out that, depending on the image texture and the type of noise corrupting the image, our method performs much better: We gain at least 25% in PSNR and at least 11% in SSIM.

The price to pay is a slight increase in terms of computational complexity of the order of image size, but this does not call real time processing into question. Due to the robustness and the speed of convergence of DCSR algorithm, its application is efficient in vital and sensitive domains such as medical imaging and remote sensing.

Our future contribution is a technological breakthrough that includes introducing a layer of intelligence at the acquisition level aimed at automatically determining image texture and its quality in terms of noise level, blur and shooting conditions (lighting, inpainting, registration, occlusion, low resolution, etc.) [START_REF] Adjabi | Past, present, and future of face recognition: A Review[END_REF]- [START_REF] Arbaoui | Wavelet-based multiresolution analysis coupled with deep learning to efficiently monitor cracks in concrete[END_REF] in order to automatically adjust the parameters necessary for an optimal use of the proposed DCSR algorithm.

  be convex and l.s.c. The following problem: « Find x * ∈ X , such that f (x * ) = min x∈X f (x) » has a solution, if and only if:

Theorem 1 . 4 .

 14 with compact support in Ω}. If the boundary ∂Ω of Ω is regular, we define: C 0 (Ω) = {f : Ω -→ Rcontinuous andf |∂Ω = 0} Definition 1.10. (Radon measure) µ is said a positive Radon measure on a σ-algebra B(X) of the Borel set, if it is a regular Borel measure on B(X) finite on compact sets of X. Radon measure have the following interesting characteristics: Let µ be a Radon measure on Ω ⊂ X. Then, µ defines a linear form L on C c (Ω) and conversely, to any linear form L on C c (Ω) corresponds a unique (positive) Radon measure µ related to L by:

Definition 1 .

 1 11. (Bounded Radon measure) The topological dual of the space C 0 (Ω), noted M(Ω), is by definition the space of bounded Radon measure.

Figure 1 . 1 :Figure 1 . 2 :

 1112 Figure 1.1: Distribution plot of Gaussian noise.

Figure 1 . 3 :

 13 Figure 1.3: Speckle Noise.

Figure 1 . 4 :Figure 1 . 5 :

 1415 Figure 1.4: Poisson noise.

Shannon, 1949 )

 1949 by choosing the right sampling frequency. We recall the Shannon-Nyquist theorem. Theorem 2.1. (Shannon-Nyquist theorem)Let x(t) be an analog signal and s e (t) it's signal sampled at the frequency ν s . If the sampling frequency ν s is greater than twice the maximum frequency ν max (i.e. ν s > 2ν max ) contained in the analogue signal, the spectra of the two signals coincide for any frequency present in the analogue signal.

Figure 2 . 1 :

 21 Figure 2.1: Compressed sensing modeling.

Figure ( 3 . 1 )

 31 Figure (3.1) shows an example of inverse problems. The example corresponds to an image restoration problem.

Figure 3 . 1 :

 31 Figure 3.1: (a).original image peppers, (b).degraded image by a blur, (c).estimated image

Figure 3 . 2 :

 32 Figure 3.2: The derivation is an ill-posed inverse problem in the sense of stability. In left: the function f : t → sin(3t) and the same function perturbed by white noise. in the right: the derivative of f and its numerical evaluation.

Figure 3 . 3 :

 33 Figure 3.3: Radon Transform

. 18 )

 18 Theorem 3.1. 1. Let A : X -→ Y , a compact operator and let 0 < a < 1 A 2 . We define the linear and bounded operators R m : Y -→ X by equation (3.18). These parameter control strategy operators µ = 1 m , m ∈ N and R m ≤ √ am. The following iterations calculate the sequence u m, = R m f :

Definition 4 . 3 .

 43 Denote by Ω an open connected set of R 2 , BV (Ω) is the subspace of functions u ∈ L 1

Figure 4 . 1 :

 41 Figure 4.1: The illustration of nonlocal self-similarity in 3D transform domain.

29 )

 29 Then; ( w, ũ, x, γ, φ) is a saddle point for L A only if the optimal conditions (4.26-4.29) are well conserved.Under these assumptions (H1-H5), we give the following theorem states the relationship between the saddle point of the problem (4.24) and the solution of the problem (4.20). Theorem 4.2. ũ ∈ R N is a solution of the problem (4.20), if and only it exists, ( w, x) ∈ (R M , R N ) and (γ, φ) ∈ (R M , R N ). Such that, ( w, ũ, x, γ, φ) is a solution of the problem (4.20). Proof of Theorem 4.2. We assume that ( w, ũ, x, γ, φ) is a solution to problem (4.24). From the left inequality of equation minimization by modifying the smooth approximation of the objective function, called the NESTA algorithm. The third algorithm, group-based sparse representation (GSR), simultaneously enforces image sparsity and self-similarity under a unified framework in an adaptive group domain. We used three images, Barbara and Cameraman gray-scale images with sizes of 256 × 256 and a color Lena image sized 512 × 512 in our experiments (see Figure 4.2).

Figure 4 . 2 :

 42 Figure 4.2: Test images. (a)Barbara, (b)Cameraman and (c)Lena.

First, we evaluate

  the performance of our DCSR algorithm by performing experiments on Barbara's image corrupted by white Gaussian noise with a standard deviation σ ranging from 20 to 80. Then, to verify the superiority of the proposed method, we compared it with the GSR algorithm, wavelet denoising and the NESTA algorithm. In Figure(4.3), each line represents the Barbara image reconstructed according to the following algorithms: DCSR (ours) in the first line, GSR in the second line, wavelet denoising in third line and NESTA in the fourth line. The level of the initial Gaussian white noise varies according to each column: (a)σ = 20, (b)σ = 50, (c)σ = 60 and (d)σ = 80.

Figure 4 . 3 :

 43 Figure 4.3: Visual comparison of the reconstruction quality for different noise levels. (a) σ = 20, (b)σ = 50, (c) σ = 60 and(d) σ = 80. The algorithms used are DCSR, GSR, wavelet denoising and NESTA: the denoised images are arranged in rows 1, 2, 3 and 4, respectively.

Figure ( 4 . 4 )

 44 Figure(4.4) provides visual results for different algorithms and several noise levels; it can be observed that the Cameraman image has been reconstructed very well by our DCSR algorithm. By analyzing the images in the fifth column of Figure(4.4), obtained by the NESTA algorithm, we can observe some artifacts above the camera head and in the upper right corner due to the nature of the NESTA algorithm itself, which is characterized by a loss of precision in high-frequency components.

Figure( 4 . 6 )

 46 Figure(4.6) presents the performance analysis of four denoising methods for grayscale images corrupted by additive white Gaussian noise. We can see that DCSR, our algorithm, exhibits the best performance with high PSNR (low noise level) and low PSNR (high noise level).

Figure 4 . 6 :

 46 Figure 4.6: Grayscale image corrupted by additive white Gaussian noise (AWGN): performance analysis of four denoising methods.

  (4.7), which illustrates PSNR output variations concerning the input PSNR.

Figure 4 . 7 :

 47 Figure 4.7: Grayscale image corrupted by salt and pepper noise: performance analysis of four denoising methods.

Figure ( 4 . 9 )

 49 Figure(4.9) shows the output results (after reconstruction) evaluated at different noise levels.DCSR algorithm provides the best denoising performance.

Figure 4 . 9 :

 49 Figure 4.9: Color image corrupted by AWGN: performance analysis of three denoising methods.

Figure( 4 . 10 )

 410 Figure(4.10) shows that the DCSR algorithm provides better visual quality results.

Figure 4 . 10 :

 410 Figure 4.10: Restoration results of the salt and pepper noise-corrupted Lena image for different values of noise (see column 1 located on the left side of the figure). The algorithms used, the proposed DCSR, wavelet denoising and GSR, are arranged in columns 2, 3 and 4, respectively.

  This performance is confirmed by Figure(4.11), which represents the variations of the output PSNR vs. the input PSNR. The robustness of our algorithm relative to the noise level is, thus, verified.

Figure 4 . 11 :

 411 Figure 4.11: Color image corrupted by salt and pepper noise: performance analysis of three denoising methods.

Figure 4 . 12 :

 412 Figure 4.12: PSNR values of Barbara's grayscale image recovered by four competing methods as a function of the number of iterations. The test image is corrupted by AWGN.

Figure 4 . 13 :

 413 Figure 4.13: PSNR values of the Cameraman grayscale image recovered by four competing methods as a function of the number of iterations. The test image is corrupted by salt and pepper noise.

Figure 4 . 14 :

 414 Figure 4.14: PSNR values of the recovered Lena color image by the competing methods vs. iteration number. The test image is corrupted by AWGN.

Figure 4 . 15 :

 415 Figure 4.15: PSNR values of the recovered Lena color image by the competing methods vs. iteration number. The test image is corrupted by salt and pepper noise.

  blocks u i , i = 1, 2, .., N and (m -1) is the number of similar blocks denoted S u i , all blocks of S u i are stacked in a matrix of size (n × m) with complexity O(n m 2 );

O(N log 2 N

 2 ) < O(N + N log 2 N ) < O(N (n m 2 + T s )) < O(N (n m 2 + T s ) + N ) (4.52)

  -25 November 2019. DOI : 10.1109/ISPA48434.2019.8966926 A. El Mahdaoui, M.S.Moulay, A.Ouahabi. " Multilevel Fast Multipole Acceleration for Fast ISAR Imaging based on Compressive Sensing", " international conference (CISTEM 2018). IEEE. 28-31 October 2018, Algiers. DOI: 10.1109/CISTEM.2018.8613411 Definition of acronyms.

	AWGN	:	Additive White Gaussian Noise
	CS	:	Compressed Sensing
	DCSR	:		Denoising Compressed Sensing by Regularizations terms
	DWT	:		Discrete Wavelet Transform
	GSR	:		Group-based Sparse Representation
	MSE	:		Mean Square Error
	NESTA	:	Nesterovs algorithm
	NLS	:		Nonlocal self-Similarity
	PSNR	:		Peak Signal-to-Noise Ratio
	SSIM	:

  Table (4.1) shows the quantitative assessment results of PSNR and SSIM relative to our algorithm, DCSR, for several values of noise levels, and the results were obtained using GSR, wavelet and NESTA. The best results are highlighted in bold type font. Table(4.1) validates that our method has superiority in image reconstruction compared with the three other methods.

Table 4 . 1 :

 41 Quality metric results on different algorithms for different values of noise levels.

	Method	Noise level PSNR SSIM
		20%	33.97	0.99
	Ours	50%	30.08	0.95
		60%	27.03	0.93
		70%	25.24	0.92
		20%	26.70	0.80
	GSR algorithm	50%	24.66	0.78
		60%	23.61	0.74
		70%	22.40	0.69
		20%	31.68	0.89
	Wavelet denoising	50%	28.61	0.86
		60%	23.61	0.79
		70%	22.40	0.74
		20%	18.88	0.39
	NESTA algorithm	50%	16.71	0.38
		60%	16.33	0.37
		70%	16.03	0.36

Table 4 .

 4 2: Computational complexity.

	Algorithms	Computational Complexity
	TV	O(N ) [33]
	wavelet denoising	O(N log

2 N ) [102] NESTA O(N + N log 2 N ) [16]

Acknowledgments

The work described in this manuscript would not have been possible without the scientific and human support of my thesis supervisors, Professor Mohammed Said Moulay and Professor Abdeljalil Ouahabi. Their skills, their sympathy and the patience with which they have supervised me have been essential for the realization of this work. I thank Professor Ammar Khemmoudj, Professor at the University of Science and Technology of Algiers (USTHB), who did me the honor of chairing the jury. I also thank Professor Abdelaziz Choutri, Professor at the Higher Normal School of Kouba (ENS), Professor Mohammed Hachama, Professor at the University of Blida and Dr. Amar Rebbouh, Associate-Professor at the USTHB for the honor they have given me by accepting to be examiners of this thesis.

I express my deepest gratitude to my mother for having encouraged me and accompanied me

Theorem 3.4. Let A a compact operator and µ( ) is chosen by the shift principe. We assume that, for all, u ∈ Im(A + A), f = Au = 0 and for a sequence n -→ 0, f n∈Y , such that: f -f n ≤ n and f n > n for all n. The corresponding Tikhonov solution u n = u µ( n), n , converges to u faster than √ n to zero, therefore

Tremendous efforts are being made to modify the original shift principle. See [START_REF] Engl | Discrepancy principles for Tikhonov regularization of ill-posed problems leading to optimal convergence rates[END_REF], [START_REF] Gfrerer | An a posteriori parameter choise for ordinary and iterated Tikhonov regularization of ill-posed problems leading to optimal convergence rates[END_REF], and others.

Landweber method

Iterative methods build a sequence of approximate solutions that (in the noiseless case) converge to the desired resolution. In the context of inverse problems, the situation is more complicated: in the presence of noise, the sequence constructed by the iterative method does not converge, in general, to a solution to the starting problem. Once again, it is necessary to regularize the iterative process, and the iteration index plays the role of regularization parameter. In other words, it is appropriate to stop the iterations earlier than one would do in a noiseless case.

In this paragraph, we examine the simplest of the iterative methods: the Landweber method [START_REF] Landweber | An iteration formula for Fredholm integral equations of the first kind[END_REF],

which is suitable for a simple analysis. Unfortunately, it converges too slowly to be usable in practice, mainly because more efficient methods exist. The two most important are the Brakhage method see [START_REF] Engl | Regularization of inverse problems[END_REF], and especially the conjugate gradient method and its variants. This last method is the most used. In the context of ill-posed problems, Kirch's book [START_REF] Kirsch | An introduction to the mathematical theory of inverse problems[END_REF] gives an accessible presentation.

Landweber [START_REF] Landweber | An iteration formula for Fredholm integral equations of the first kind[END_REF], Fridman [START_REF] Fridman | Method of successive approximations for Fredholm integral equations of the first kind[END_REF] proposed to rewrite the equation Au = f as: u = (I -aA + A)u + aA + f f or a > 0 (3.15)

The iterative scheme of this equation is the following:

Then, Ψ is differentiable in the sense of Frechet for all y ∈ X and Ψ (y)u = Re(Ay -f, Ay) = Re(A + (Ay -f ), y), x ∈ E (3.17)

The linear functional Ψ (y) can be identified with A + (Ay -f ) on the Hilbert space X.

Chapter 4

Image denoising using a compressive sensing approach based on regularization constraints

Convexity For any u 1 , u 2 and t ∈ [0, 1],

Homogeneity It is obvious that for each u and t > 0,

We will now consider discrete images (the practice case). A digital image, or discrete, is a vector of two dimensions of size N × N . Let X be the Euclidean space R N ×N and Y = X × X. We give X with the usual scalar product:

and the norm:

We will give a discrete formulation of what was done before and, in particular, define a discrete total variation that we will note TV. We first introduce a discrete version of the operator gradient. If u ∈ X, the gradient ∇u is a vector of Y defined by:

where

The discrete total variation is then given by norm p of the discrete gradient:

where u ∈ R N represents an image, and

with D h and D v representing the discrete gradient [START_REF] Aubin | L'analyse non linéaire et ses motivations économiques[END_REF] of the image in the horizontal and vertical direction, respectively.

The l p norm could be the l 1 norm corresponding to the anisotropic TV [START_REF] Beck | Fast Gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] or the l 2 norm corresponding to the isotropic TV [START_REF] Zuo | A Generalized accelerated proximal gradient approach for total variation based image restoration[END_REF]. By definition, the l p norm is

In this work, we consider p to be equal to 1; in another word, we used an anisotropic TV.

Total variation regularization can be a generalized l1 regularization in compressive sensing problems.

Instead of assuming the signal is sparse, the premise of TV regularization is that the gradient of the underlying signal or image is sparse. In other words, total variation measures the discontinuities and Chapter 4

Image denoising using a compressive sensing approach based on regularization constraints seeks the solution with the sparsest gradient.

Numerically, the gradients are done by discretization on a horizontally and vertically oriented grid.

Note the similarity of this approach, in its discretized version, with the fused LASSO or F-LASSO method, introduced by Tibshirani, Saunders, Rosset, Zhu, and Knight [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] (although the authors add the penalty term l1 also to obtain the sparsity of the solution). The l1 norm then forces the resolution of the minimization problem to be piecewise constant. Indeed, the notion of sparsity induced by the l1 norm is applied to the image differences and leads to one of the significant shortcomings of this type of approach, creating a staircase effect (in Staircasing Effect). Many solutions have been proposed to circumvent this defect while keeping the denoising qualities of the TV method (the possibility of keeping discontinuities). We refer to the article by Louchet and Moisan [START_REF] Louchet | Total variation denoising using posterior expectation[END_REF] for more details.

TV denoising model

The minimization problem of total variation is :

Where λ is a coefficient balancing the effects of the regularization term and the attachment term to the data, it is used to choose the allowable noise level.

The main interest of this choice of regularization consists in the sharpness of the images it reconstructs and its ability to reconstruct images from incomplete data. In addition, this regularization tends to reconstruct piecewise constant images because of the norm l 1 applied to the image gradient, reinforcing the gradient's sparsity.

The image restoration model based on total variation regularization yields piecewise constant images.

This effect is called the 'staircasing'. Smooth regions in the original image are recovered as piecewise smooth regions. To overcome this difficulty, we propose an alternative to the total variation, the nonlocal self similarity which reduces, the staircasing phenomenon.

Nonlocal self-similarity

While Local denoising methods have low time complexity, their performance is limited when the noise level is high. The reason is that the correlations of the neighborhood pixels are seriously disturbed by the high noise level.

In recent years, the exploitation of non-local self-similarity (NSS) in images has significantly improved image denoising performance, it refers to the fact that, for a given local patch in natural image, many similar patches can be found throughout the image.

Nonlocal self-similarity [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF] is another property of natural images, which means that patches tend to repeat in different locations, this notion was first used in texture synthesis by Efros and Leung [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF] Chapter 4

Image denoising using a compressive sensing approach based on regularization constraints and then by Criminisi, Pérez, and Toyama [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] while determining the value of the hole using similar image patches and filling in the gaps in the images, which also influences the image denoising task.

Based on this idea, Buades and al. [START_REF] Buades | A non-local algorithm for image denoising[END_REF] proposed an efficient denoising model called nonlocal means (NLM), which takes advantage of this image property to perform a type of weighted filtering for denoising tasks using the degree of similarity between the surrounding pixels.

The nonlocal self-similarity prior effective for image processing, such as the state-of-the-art denoising method BM 3D [START_REF] Dabov | Image denoising by sparse 3D transform-domain collaborative filtering[END_REF] is based on this prior, it is the most efficient method at the moment; it is based on optimizing a very large number of parameters (a dozen). Moreover, it isn't easy to measure the impact of the various components separately. It is, therefore, a method with a practical rather than theoretical purpose.

A nonlocal self-similarity is a significant property of natural images too. It was obtained in the specifc process is that:

First, assuming the original image u of size N is divided into many overlapping blocks

Then, find (m -1) similar blocks. Instead of using a tunable threshold to choose similar blocks in for denoising, our choice with a fixed number is simple and robust to the similarity criterion.

Thus, for simplicity, the measure for calculating the similarity between different blocks is Euclidean distance. Here, the Euclidean distance between the search block and the reference block must satisfy the following formula:

Where: S u i is the set of similar image blocks for reference block u i , T is the metric threshold.

is the Euclidean distance between the block u i and u i j, and can be calculated by:

The block u i j is similar to the block u i , if the Euclidean distance d(u i , u ij ) is less than the metric threshold T . Then, group or stack (hence the name 3D) similar patches in groups (block creation step), (m -1) similar blocks comprise set S u i in the training window with (L × L) size.

Finally, all blocks of S u i are stacked into a 3D array

where T 3D is a transform operator, and T 3D (Z u i ) includes the transform coefficients for (Z u i ). θ u is the column vector of the lexicographically stacked representation of all 3D transform coefficients.

We will not only give an estimator by means of patches similar to the patch of interest: It is rather

Chapter 4

Image denoising using a compressive sensing approach based on regularization constraints (4.24), we have:

Thus, the inequality on the right of equation (4.24) shows that:

Let Du = w equal u = x in equation (4.31), ũ is a solution of problem (4.19).

Conversely, we assume that ũ is a solution. We take, Dũ = w and ũ = x; there exist γ ∈ ∂R 1 (Du)

and φ ∈ ∂R 2 (u). Then, we can verify that (ũ, w, x, γ, φ) is a saddle point of L, and the proof is complete.

To solve (4.20), we use the augmented Lagrangian method iteratively as follows:

Here, subscript k denotes the iteration index, and γ and ϕ are the Lagrangian multipliers associated with the constraints Du = w, u = x, respectively.

The alternative direction method is introduced to solve the problem efficiently. Due to the nondifferentiability of equation (4.21), the alternative direction method is introduced to solve the problem efficiently, which alternatively minimizes one variable while fixing the other variables to split equation In what follow, we describe how to solve efficiently solve these problem.

Uptdate w Given u and x, we obtain subproblem w arg min

The optimization problem described by equation (4.33) can be solved using the shrinkage formula [START_REF] Li | An efficient augmented Lagrangian method with applications to total variation minimization[END_REF].

Then, w can be obtained as follows:

where max{•} represents the the larger number between two elements, and sgn(•), is a piece wise function that is defined as follows:

Uptdate u Fixing w and x, the optimization associated with u is:

equation (4.36) is u quadratic, so minimizing the subproblem is simplified to a linear system.

The matrix on the left-hand side of the above system is positive, definite and tridiagonal, since D T D is a positive semidefinite tridiagonal matrix. Moreover, µ and β are both positive scalars.

Uptdate x Given u, we obtain the x subproblem as follows:

arg min

By applying a completing square method and omitting all constants independent of x, the subproblem defined in equation (4.38) can be simplified as follows:

Considering r = u -c, with c = ϕ β , as a noisy observation of x, error (or noise) e = x -r follows a probability law that is not necessarily Gaussian, with zero mean and variance σ 2 . According to the central limit theorem or law of large numbers, the following equation holds:

with e, x, r ∈ R N and θ x , θ r ∈ R i for i = 1, ...., N .

Incorporating (4.40) into (4.39) results in the following:

Chapter 4

Image denoising using a compressive sensing approach based on regularization constraints Since θ x is component-wise separable and an unknown variable, according to [START_REF] Li | User's guide for TVAL3: TV minimization by augmented Lagrangian and alternating direction algorithms[END_REF], the closed-form of the minimization problem of equation (4.41) can be written as follows:

Thus, the solution of the x subproblem of equation (4.38) is as follows:

where Ω is the reconstruction operator.

Based on the discussions above, we get the algorithm for solving (4.19) shown below:

Algorithm 3 Our algorithm DCSR Input: The measurement f and the linear measurement matrix A. The comparative performance of competing methods for the recovery of noisy images is discussed in the next section.

Experimental results

In this section, we verify the efficiency and practicability of the proposed method by describing experiments with simulated and real data sets. The proposed algorithm's performance, DCSR [START_REF] Mahdaoui | Image denoising using a compressive sensing approach based on regularization constraints[END_REF],

Image Quality Metrics

Peak signal-to-noise ratio and mean square error (MSE) have long been used as fidelity metrics in the image processing community. The formulas are simple to understand and implement; they are Chapter 4

Image denoising using a compressive sensing approach based on regularization constraints easy and fast to compute, and minimizing MSE is also very well understood from a mathematical point of view.

PSNR (peak signal-to-noise ratio, unit: dB) [START_REF] Ferroukhi | Medical video coding based on 2nd-generation wavelets: Performance evaluation[END_REF] is the ratio between the maximum possible power of a signal and the power of noise. Because of its generalization, it is possible to compare various methods using this criterion. Indeed, most of the methods presented in the image processing literature are evaluated by this criterion. PSNR facilitates the comparison of performances with previous methods.

Higher PSNR value means better visual quality. PSNR is defined as follows:

MSE is the mean square error between initial image u and estimated image û of size M × N , i and j represent the image row and column pixel position,respectively, k is the number of bits of each sample value.

Structural similarity (SSIM) textures [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF] have proven to be a better error metric for comparing the image quality with better structure preservation. They are in the range of [0, 1], which is a value closer to one indicating better structure preservation:

such that l(i, j) is the luminance comparison defined as follows:

where µ i and µ j are functions of the mean intensities of signalsi and j, respectively. c(i, j), the contrast comparison, is a function of standard deviations σ i and σ j , and it is defined as the following form:

The structure comparison s(i, j) is defined as follows:

where µ i and µ j are the mean value of images u and û,respectively. σ 2 i and σ 2 j represent the variance of u and û, respectively. σ ij is the covariance of images u and û. c 1 , and c 2 and c 3 are constant. Specifically, it is possible to choose c i = K 2 i D 2 with i = 1 and K i 1 and, D as the pixel values' dynamic range.

Quantitative Assessment

In this sub-section, we evaluate the quality of image reconstruction. We compare these methods quantitatively; the peak signal-to-noise ratio and structural similarity indices are calculated for images with different algorithms. Remember that, at first, Barbara was contaminated by Gaussian Chapter 4

Image denoising using a compressive sensing approach based on regularization constraints

• Compression of the image.

This degradation can be canceled or at least reduced by proceeding to a preprocessing procedure called a denoising operation or simply denoising. Such an operation allows, on the one hand,producing a perception of a quality image and, on the other hand, the improvement inthe performance of subsequent image processing (extraction of the desired information, prediction, classification, texture analysis, segmentation, etc).

With this in mind, the authors of reference [START_REF] Zhu | A novel fast single image dehazing algorithm based on artificial multiexposure image fusion[END_REF] propose a new image denoising method called dehazing because it tends to eliminate haze due to bad weather conditions. This original method is based on the application of artificial multiexposure image fusion [START_REF] Kaur | Image fusion techniques: A survey[END_REF] involving local and global image details. Such an approach allows the recovery of quality images but, unfortunately, halos or artifacts often appear near the edges when the inputs are sparse, which makes postprocessing (linear saturation adjustment) introduced by the authors ineffective.

Wavelet image denoising [START_REF] Bnou | A wavelet denoising approach based on unsupervised learning model[END_REF] is powerful for edge detection in three preferred directions: diagonal, vertical and horizontal. For this purpose, several types of wavelets exist, such as the Haar wavelet that preserves edge information, but technically it is not continuous and is not differentiable. This wavelet can be identified with the optimization problem using an l 1 norm. In contrast, the Morlet wavelet, which is continuous, can be identified with the use of the l 2 norm [START_REF] Ouahabi | Signal and Image Multiresolution Analysis[END_REF].

Our method mixed the total variation with nonlocal self-similarities to recover the image without artifacts and preserve details and textures of the image. On a psycho-emotional level, the visual quality of an image is necessary for the observer. In this respect, it is interesting to note that the images reconstructed in Figure 4.3 by the NESTA algorithm are visually unpleasant and lose some important details.Limitations also appear during wavelet denoising: The noise contained in the image cannot be removed if the standard deviation is set too large (above 50), i.e., if the noise level is high.

Concerning GSR (second row of Figure 4.3), we can observe that this algorithm is quite efficient in removing noise; however,it has a slight loss of contrast. On the other hand, it can be observed that our algorithm, DCSR, is efficient in removing even noise with standard deviations equal to 80. It is, therefore, confirmed that the proposed method provides the most visually satisfying results for both edges and textures.

On an objective and therefore measurable level, Table ( Chapter 4

Image denoising using a compressive sensing approach based on regularization constraints

However,when we compared our method, DCSR, in terms of computational complexity, our approach is not the most advantageous; however, real-time processing is largely feasible. Hence, complexity reduction techniques such as block-compressed sensing [START_REF] Li | Overview of compressed sensing: Sensing model, reconstruction algorithm, and its applications[END_REF] or deep learning technique [START_REF] Tian | Deep learning on image denoising: An overview[END_REF] are feasible.

We applied our algorithm in an ablation study to clarify the effect of total variation (TV) and nonlocal self-similarity (NLS) on the compressed sensing (CS) recovery model. First, we cancel all regularization constraints and leave only the CS alone. Next, we analyze CS coupled with TV and canceled NLS. Then, we delete TV, replaced it with NLS and kept CS. Table (4.3) summarizes this ablation analysis using quantitative values: PSNR (dB) and SSIM. From Table (4.3), we can conclude that the presence of TV or NLS improves the quality of the reconstruction with CS, while noting that the CS+NLS coupling performs better than CS+TV. On the other hand, by comparing the four scenarios with each other, we can conclude that the reconstruction, with both regularization functions simultaneously (CS+TV+NLS), is significantly better than the other three scenarios: it effectively removes noise and improves the robustness of our approach.