Liyan Ma
email: maliyan@ime.ac.cn

Lionel Moisan
email: lionel.moisan@parisdescartes.fr

Jian Yu
email: jianyu@bjtu.edu.cn

Tieyong Zeng
email: zeng@hkbu.edu.hk

A STABLE METHOD SOLVING THE TOTAL VARIATION DICTIONARY MODEL WITH L ∞ CONSTRAINTS

Keywords: 2000 Mathematics Subject Classification: Primary: 68U10; Secondary: 49N45, 94A08, 65T60 Image restoration, total variation, nonlocal total variation, proximal gradient method, wavelet packet decomposition

Image restoration plays an important role in image processing, and numerous approaches have been proposed to tackle this problem. This paper presents a modified model for image restoration, that is based on a combination of Total Variation and Dictionary approaches. Since the well-known TV regularization is non-differentiable, the proposed method utilizes its dual formulation instead of its approximation in order to exactly preserve its properties. The data-fidelity term combines the one commonly used in image restoration and a wavelet thresholding based term. Then, the resulting optimization problem is solved via a first-order primal-dual algorithm. Numerical experiments demonstrate the good performance of the proposed model. In a last variant, we replace the classical TV by the nonlocal TV regularization, which results in a much higher quality of restoration.

Introduction

A classical inverse problem in modern imaging sciences is image restoration, the aim of which is to recover an ideal image u ∈ L 2 (Ω) from a degraded observation g ∈ L 2 (Ω), knowing that (1)

g = Au + b,
where Ω is a connected bounded open subset of R 2 with compact Lipschitz boundary, and b ∈ L 2 (Ω) is a realization of a white Gaussian noise with standard deviation σ.

Here A is a known linear operator from L 2 (Ω) to L 2 (Ω), that depends on the considered inverse problem: A = I (identity) for denoising, A is convolution operator for deblurring, and A is the canonical projection on known pixels in the case of inpainting.

Image restoration is an ill-posed problem. Directly solving (1) may not yield meaningful solutions. To overcome this difficulty, many regularization methods have been proposed in the literature, e.g. Tikhonov's regularization [START_REF] Tikhonov | Solution of Ill-Posed Problems[END_REF], the Total Variation (TV) criterion [START_REF] Osher | Nonlinear total variation based noise removal algorithms[END_REF], non-local graph regularization [START_REF] Peyre | Non-local regularization of inverse problems[END_REF], and wavelet-based sparsity assumptions [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF]. In this paper, we mainly focus on the TV approach, first proposed by Rudin, Osher and Fatemi [START_REF] Osher | Nonlinear total variation based noise removal algorithms[END_REF]:

(2) min

u T V (u) + λ 2 Au -g 2 2 .
Here the Total Variation of an image u ∈ L 2 (Ω) is defined by

T V (u) = Ω |Du|,
the gradient Du being understood in the distribution sense.

TV regularization has been extensively used in image processing applications, for it is a simple term that manages to preserve sharp features in images. However, due to the non-differentiability of the TV operator, the functional (2) is difficult to minimize. In the past decades, many methods have been developed to address this problem. To name a few, let us mention variable splitting methods [START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF][START_REF] Goldstein | The split Bregman method for l1 regularized problems[END_REF][START_REF] Setzer | Operator splittings, Bregman methods and frame shrinkage in image processing[END_REF][START_REF] Wu | Augmented Lagrangian method for total variation restoration with non-quadratic fidelity[END_REF], the alternating direction method of multipliers (ADMM) [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF][START_REF] Ng | Inexact alternating direction methods for image recovery[END_REF], dual methods [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] and primal-dual methods [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF][START_REF] Chan | A nonlinear primal-dual method for total variationbased image restoration[END_REF][START_REF] Zhu | An efficient primal-dual hybrid gradient algorithm for total variation image restoration[END_REF]. Some of these methods have strong connections, which are studied in [START_REF] Setzer | Operator splittings, Bregman methods and frame shrinkage in image processing[END_REF][START_REF] Tai | Augmented Lagrangian method, dual methods and split Bregman iteration for ROF model[END_REF].

Besides TV regularization, another famous approach for image restoration is the wavelet-based method. Since wavelet bases can provide a multiscale and sparse representation for images (see [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF]), they have been widely adopted in image restoration [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF][START_REF] Malgouyres | A framework for image deblurring using wavelet packet bases[END_REF][START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelets for image restoration[END_REF][START_REF] Malgouyres | Minimizing the total variation under a general convex constraint for image restoration[END_REF] after their introduction by Donoho and Johnstone [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF]. Later, Chambolle et al. [START_REF] Chambolle | Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage[END_REF] and Steidl et al. [START_REF] Steidl | On the equivalence of soft wavelet shrinkage, total variation diffusion, total variation regularization, and sides[END_REF] established the relationship between wavelet-based and TV-based methods. Eventually, by combining TV and wavelet methods, Candes and Guo [START_REF] Candes | A new multiscale transform, minimum total variation synthesis: application to edge-preserving image reconstruction[END_REF] and Malgouyres [START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelets for image restoration[END_REF] proposed a hybrid approach for image restoration. The optimization model in [START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelets for image restoration[END_REF] is

(3) min u T V (u) s.t. Au -g D,∞ ≤ τ ,
where D ⊂ L 2 (Ω) is a finite dictionary, the parameter τ > 0 depends on the noise level, and u D,∞ = sup Ψ∈D | u, Ψ |. As pointed out by Malgouyres in [START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelets for image restoration[END_REF],

taking the (infinite) dictionary D = w ∈ L 2 (Ω), w 2 = 1 would reduce (3) to the ROF model. In the same paper, he showed that when the dictionary contains wavelet/wavelet packet bases and their opposites, the new model preserves textures better than the ROF model. The choice of the dictionary is very important, since a carefully selected dictionary can improve the quality of restoration. In order to obtain optimal results through (3), Zeng and Malgouyres [START_REF] Zeng | Using Gabor dictionaries in a TV-L ∞ model for denoising[END_REF] investigated twelve Gabor dictionaries for [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF], and the experimental results clearly demonstrated that the choice of the dictionary deeply affects the performance of this model. Later, Zeng and Ng [START_REF] Zeng | Incorporating known features into a total variation dictionary model for source separation[END_REF][START_REF] Zeng | On the total variation dictionary model[END_REF] provided a theoretical study on the model [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF] and the choice of the dictionary. They proved the existence of solutions of (3), and showed that the dictionary should represent sparsely the curvatures of the ideal image in order to achieve a better denoising performance.

The minimization problem (3) is evidently a rather difficult task, and it has been rarely addressed in the literature. In order to approximate the solution of (3), Lintner and Malgouyres [START_REF] Lintner | Solving a variational image restoration model which involves contraints[END_REF] proposed the following LM model:

(4)      min u T V β (u) 2 := Ω ϕ β (|∇u|) + β u, 1 2 2 s.t. ∀Ψ ∈ D, Au -g, Ψ ≤ τ
where 1 denotes a constant vector with all coordinates equal to 1, β is a small positive number, and ϕ β (t) = t 2 + β 2 . As the objective function in ( 4) is differentiable, convex and coercive, a dual Uzawa method was used in [START_REF] Lintner | Solving a variational image restoration model which involves contraints[END_REF] to solve this optimization task. Note that in (4), as well as in all this paper, the dictionary D is assumed to be symmetric, which means that

(5) ∀ψ ∈ D, -ψ ∈ D.
Hence, the constraints in (4) could as well be written

| Au -g, Ψ | ≤ τ .
Although improving the quality of the recovered image compared to ROF Model, the approach proposed in (4) still has some drawbacks. First, the use of T V β (u) 2 to replace T V (u) causes a well-known unpleasant smear effect; second, the term β u, 1 introduced in (4) does make the objective function coercive (this is critical to make the Uzawa method work), but it also imposes an extra assumption on the mean value of the ideal image u which might affect the quality of the restored image; third, the numerical scheme adopted in [START_REF] Lintner | Solving a variational image restoration model which involves contraints[END_REF] may be unstable in some situations, and imposes a careful (small enough) choice of the gradient-ascent time step used in Uzawa method.

The goal of this paper is to fix these issues. We propose a modified model which integrates the ROF and the dictionary model (3) seamlessly. Moreover, inspired by the recent works of Chambolle and Pock [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF], we give two stable and effective algorithms to solve the associated optimization problem, one for the general case and one in the special case of image denoising (A = I).

The outline of the paper is as follows. In Section 2, we introduce the proposed modified model and study some basic properties. In Section 3, some notations and preliminaries (especially in the discrete settings) are presented, which will be used in the rest of the paper. In Section 4, we briefly review Chambolle-Pock Algorithm, then we describe how it can be applied to the proposed model (Section 5). We propose two algorithms that solve the corresponding unconstrained and constrained models respectively. In Section 6, we report various numerical results to demonstrate the good performance of the proposed approach. Finally, we briefly investigate in Section 7 an interesting variant of the proposed model, where the TV regularization term is replaced by the nonlocal TV [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF], before we conclude in Section 8.

Proposed model

This section presents the primal formulation of the proposed model and a study of its basic properties. In Section 5, we will solve the proposed model in its primaldual formulation to exactly preserve its properties.

2.1. Model. Using TV regularization rather than the T V β (u) 2 term, the minimization problem (4) becomes ( 6)

min u T V (u) s.t. ∀Ψ ∈ D, Au -g, Ψ -τ ≤ 0.
From a statistical point of view, since we are dealing with additive white Gaussian noise, the probabilistic formulation of (1) is

p (g|u) = x∈Ω 1 √ 2πσ exp - (Au(x) -g(x)) 2 2σ 2 .
The corresponding energy function,log p(g|u), is, up to a constant additive and multiplicative term, the data-fidelity term Au-g 2 of the ROF model. We propose to add this natural term to the energy used in [START_REF] Boyd | Convex Optimization[END_REF] in order to control the coercivity and to enforce the data-fidelity constraint on the restored image. Thus, we obtain a modified model defined by the primal formulation ( 7)

   min u T V (u) + λ 2 Au -g 2 2 s.t. ∀Ψ ∈ D, Au -g, Ψ -τ ≤ 0.
Both parameters λ and τ control the regularization induced by the model, but in different ways. For fixed τ , a small value of λ makes the restored image less noisy, but may provide a solution that is too smooth; on the contrary, a large value of λ tends to produce sharper but also more noisy results. We could probably use an estimation of the noise level to automatically set the value of λ (as it is done in some other works), but this is beyond the scope of this paper. In this work, we decide to tune it empirically. Model [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF] seamlessly integrates the ROF model and the original TV dictionary model. Indeed, if λ = 0, it is exactly (3); if the dictionary is empty, then it reduces to the classical ROF model. We also avoid the extra assumption needed in (4). Experiments will later demonstrate that the least-squares term is very important to stabilize the algorithm and to provide a high-quality restored image.

2.2.

Existence of a solution. We here prove the existence of a solution for the proposed model [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF] in the continuous settings. First, we recall the definition of the space of functions with bounded variation (BV), which plays an important role in image modeling. Note that a rather complete discussion on this topic can be found in [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problem[END_REF] and [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF].

Definition 2.1. The space of functions of bounded variation is

BV (Ω) = u ∈ L 1 (Ω) ; Ω |Du| < ∞ ,
where

Ω |Du| = sup Ω u div ϕ dx; ϕ = (ϕ 1 , ϕ 2 ) ∈ C ∞ 0 (Ω) 2 , ϕ L ∞ (Ω) ≤ 1 .
Endowed with the norm u BV = Ω |u|+ Ω |Du|, BV (Ω) is a Banach space, and T V (•) is a proper convex lower semi-continuous functional from L 2 (Ω) to (-∞, +∞] (see Lemma 8 in [50]).

S τ := {u ∈ L 2 (Ω); Au -g 2 < τ }.
is nonempty and A 1 = 1. If D is a finite symmetric subset of the closed unit ball of L 2 (Ω), then for any λ ≥ 0, Problem [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF] admits at least one solution.

Proof. • We first consider the (simpler) case λ > 0. The ROF functional

(9) u → T V (u) + λ 2 Au -g 2 2
is convex, closed (proper and lower semicontinuous), and coercive (because A 1 = 1, see [START_REF] Chan | Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods[END_REF]) on L 2 (Ω), so it admits a minimizer (see [START_REF] Ekeland | Convex Analysis and Variational Problems[END_REF] p.35). Now let us consider the feasible set (10)

C τ = u ∈ L 2 (Ω); ∀Ψ ∈ D, Au -g, Ψ ≤ τ .
Obviously, C τ is a nonempty closed convex set, and we can rewrite [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF] as

min u T V (u) + λ 2 Au -g 2 2 + δ Cτ (u),
where δ Cτ (u) is the indicator function of C τ (which takes the value 0 if u ∈ C τ , and the value +∞ otherwise). The functional above is, like the ROF functional (9), convex, closed, and coercive on L 2 (Ω). Therefore, (7) admits at least one solution.

• Now let us consider the case λ = 0. In [START_REF] Zeng | On the total variation dictionary model[END_REF], Zeng and Ng proved that there exists a solution of (7) when A = I. Here, we prove that this remains true when A is a linear operator such that S τ is nonempty. Denote as in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] by C τ the feasible set of [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF], and take v * ∈ S τ . As

C ∞ 0 is dense in L 2 (Ω) and A is bounded, we can choose a v 0 ∈ C ∞ 0 such that A(v 0 -v * ) 2 < τ -Av * -g 2 . Then we have Av 0 -g 2 ≤ A(v 0 -v * ) 2 + Av * -g 2 < τ
and by Cauchy-Schwarz inequality, Av 0g, Ψ < τ for any Ψ ∈ D (because Ψ 2 ≤ 1). Hence v 0 ∈ int(C τ ) and since Ω is bounded, we also have T V (v 0 ) < +∞. Hence, Slater condition [START_REF] Zalinescu | Convex Analysis in General Vector Spaces[END_REF] holds and there exists a minimizing sequence (u j ) in C τ such that lim j→∞ T V (u j ) = M , where M is the minimum value in (7) (recall that λ = 0 here). Note that we can assume without loss of generality that T V (u j ) ≤ M + 1 for all j.

Denote by m Ω (f ) the mean of f on Ω. Since Ω is a connected bounded open subset of R 2 with compact Lipschitz boundary, it is an extension domain (using Prop. 3.21 of [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problem[END_REF]) and by Poincaré inequality (see [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problem[END_REF], Remark 3.50]), there exists a constant c 0 > 0 such that [START_REF] Bugeau | A Comprehensive Framework for Image Inpainting[END_REF] ∀j,

u j -m Ω (u j ) L 2 (Ω) ≤ c 0 • T V (u j ) ≤ c 0 (M + 1).
From [START_REF] Bugeau | A Comprehensive Framework for Image Inpainting[END_REF], it is easy to deduce that u jm Ω (u j ) BV is bounded. As every bounded sequence in BV(Ω) is relatively compact in L 1 (Ω) (see [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF]), there exists a subsequence u which converges to some w * ∈ BV(Ω) both in L 1 (Ω), and weakly in L 2 (Ω) (thanks to [START_REF] Bugeau | A Comprehensive Framework for Image Inpainting[END_REF]). By lower semicontinuity of TV with respect to the L 1 (Ω) topology (see [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF]), one obtains

T V (w * ) ≤ lim inf j→∞ T V u (1) j -m Ω u (1) j = lim inf j→∞ T V u (1) j = M.
Now let us prove that for some constant c, the image w * + c is in the feasible set C τ . As u (1) j ∈ C τ , we know that for any Ψ ∈ D, we have [START_REF] Candes | A new multiscale transform, minimum total variation synthesis: application to edge-preserving image reconstruction[END_REF] u

(1)

j -m Ω u (1) j , A * Ψ + m Ω u (1) j Ω A * Ψ ≤ τ + g, Ψ .
If Ω A * Ψ = 0 for any Ψ ∈ D, then taking the limit in (12) yields

w * , A * Ψ ≤ τ + g, Ψ ,
so that w * is in the feasible set and it is a solution of [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF].

Otherwise, we have α := Ω A * Ψ 0 = 0 for some Ψ 0 ∈ D, and using ( 12) with Ψ = εΨ 0 yields, thanks to [START_REF] Bugeau | A Comprehensive Framework for Image Inpainting[END_REF],

∀j, ∀ε ∈ {-1, 1}, εα • m Ω u (1) j ≤ τ + g, εΨ 0 + c 0 (M + 1) A * 2 ,
which proves that m Ω u

(1) j is bounded. Hence, we can choose a subsequence

u (2) j such that m Ω u (2) j
converges to some constant c. Taking the limit for this subsequence in [START_REF] Candes | A new multiscale transform, minimum total variation synthesis: application to edge-preserving image reconstruction[END_REF], we obtain

∀Ψ ∈ D, w * , A * Ψ + c Ω A * Ψ ≤ τ + g, Ψ , so that w * + c belongs to C τ ( A(w * + c) -g, Ψ ≤ τ )
and it is a solution of (7).

Notations and Preliminaries

Now let us go to the discrete settings. Before giving details of the algorithms, we need to introduce some notations and preliminaries which will be used in the following sections. In order to simplify the notations, we here focus on gray-level images. Corresponding notions for color images, that will be considered in the experiments, can be found in [START_REF] Bresson | Fast dual minimization of the vectorial total variation norm and applications to color image processing[END_REF].

3.1. Notations. Discrete images will be represented by elements of the Euclidean space E = R m×n . Let •, • denote the standard inner product, and • p be the standard l p -norm. The transpose (dual operator) A * of a linear transform A is defined as usual by Au, g = u, A * g . I denotes the identity matrix.

The discrete gradient and divergence operators. Denote

U = E × E. The gradient operator is a map ∇ : E → U defined by ∇u = (∂ x u, ∂ y u)
T as indicated below. The adjoint of the gradient (∇ * ) is the opposite of the divergence operator (-div), so that for any z

= (p q) T ∈ U , div(z) = -(∂ * x p + ∂ * y q). The maximum norm of z is defined by z ∞ = max i,j p 2
i,j + q 2 i,j . We shall consider two different boundary conditions for the discretization of ∇ and div = -∇ * , namely: 1) Symmetric Boundary Conditions. The discrete gradient operators ∂ x and ∂ y are defined as

(∂ x u) i,j = u i+1,j -u i,j if 1 ≤ i < m, 0 otherwise, and
(∂ y u) i,j = u i,j+1 -u i,j if 1 ≤ j < n, 0 otherwise.
The discrete operators ∂ * x and ∂ * y are defined accordingly as

(∂ * x p) i,j =    p i,j if i = 1, p i,j -p i-1,j if 1 < i < m, -p i-1,j if i = m, and
(∂ * y q) i,j =    q i,j if j = 1, q i,j -q i,j-1 if 1 < j < n, -q i,j-1 if j = n.
2) Periodic Boundary Conditions. In that case, all operators are defined using a circular convention for the coordinates, that is, m+1 := 1 for the i coordinate, and n + 1 := 1 for the j coordinate:

(∂ x u) i,j = u i+1,j -u i,j if 1 ≤ i < m, u 1,j -u i,j if i = m, (∂ y u) i,j = u i,j+1 -u i,j if 1 ≤ j < n, u i,1 -u i,j if j = n. (∂ * x p) i,j = p i,j -p m,j if i = 1, p i,j -p i-1,j if 1 < i ≤ m, (∂ * y q) i,j = q i,j -q i,n if j = 1, q i,j -q i,j-1 if 1 < j ≤ m.

Proximal operators.

The proximal operator, introduced by Moreau [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF][START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF], is an important concept in convex optimization and has been widely used recently [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF][START_REF] Combettes | Signal recovery by proximal forward backward splitting[END_REF][START_REF] Micchelli | Proximity Algorithms for Image Models: Denoising[END_REF]. Definition 3.1 (Proximal Operator). Let ϕ be a convex, proper and closed (hence lower semi-continuous) function on R d . The proximal operator of ϕ is defined by

∀x ∈ R d , prox ϕ (x) = arg min y ϕ(y) + 1 2 x -y 2 2 .
We recall the following useful result for the calculation of prox f +δC , that plays a key role in solving the constrained formulation of the model we propose. Proposition 2. (see [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF], Proposition 2.1) Let f : R → R ∪ {+∞} be a convex, proper and closed function, and let C be a closed convex subset of R such that C ∩ domf = ∅. Then prox f +δC = P C • prox f , where P C is the projection on C.

The primal-dual problem.

Most recent papers in the literature [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Chan | A nonlinear primal-dual method for total variationbased image restoration[END_REF][START_REF] Zhu | An efficient primal-dual hybrid gradient algorithm for total variation image restoration[END_REF] use the dual or the primal-dual formulation for TV-related convex optimization problems, because they can reach better performances than the primal formulation for this kind of non-differentiable convex functionals. This section shortly reviews the primal-dual framework.

Let us consider the primal minimization problem

(13) min u F (Ku) + G(u),
where K is a continuous linear map from E = R m×n to U = E × E with norm Thus the saddle-point optimization problem or the primal-dual formulation of ( 13) can be written [START_REF] Chambolle | Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage[END_REF] min

K 2 = max { Ku 2 ; u ∈ E, u 2 ≤ 1} ,
u∈E max p∈U Ku, p + G(u) -F * (p).

Review of Chambolle-Pock algorithm

In this section, we briefly review the first-order primal-dual algorithm proposed by Chambolle and Pock [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF]. This algorithm is designed to solve the convex optimization problem [START_REF] Chambolle | Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage[END_REF]. Let (û, p) ∈ E × U be a solution of the problem ( 14), then we have the optimality conditions:

K û ∈ ∂F * (p), -K * p ∈ ∂G(û).
Assuming that prox F * and prox G are easily computed, Chambolle and Pock solve the optimization problem [START_REF] Chambolle | Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage[END_REF] with respect to u and p respectively. To accelerate the algorithm, they use a multi-step scheme adopted in many algorithms [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O 1/k 2[END_REF] to replace the classical one-step scheme. Algorithm 1 below summarizes the general scheme proposed by Chambolle and Pock [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF], when the parameter θ they introduce is set to 1.

Algorithm 1 General scheme of Chambolle-Pock Algorithm

Initialization: ū0 = u 0 , p = 0, t p > 0, t u > 0, and

t p t u K 2 2 < 1. Output: u N . Iterate for n = 0, 1, . . . N -1: Step 1. p n+1 = prox tp•F * (p n + t p • K ūn ) Step 2. u n+1 = prox tu•G u n -t u • K * p n+1 Step 3. ūn+1 = 2u n+1 -u n
Step 3 in Algorithm 1 can be seen as an approximate intergradient step. The authors proved that Algorithm 1 converges at the rate O(1/n) on the primal iterates.

Remark 1. We only presented here the basic algorithmic framework proposed by Chambolle and Pock in [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF]. Their paper also describes some variants which can achieve O(1/n 2 ) convergence under certain conditions.

Application of the Chambolle-Pock algorithm for our model

In this section, after providing the discrete form of the proposed model, we propose two algorithms for solving it. One solves the corresponding unconstrained formulation of the model, while the other solves the constrained formulation. 5.1. Discrete form. The discrete TV is defined as T V (u) = ∇u 1 , that is, T V (u) = F (Ku) with K = ∇ and F (•) = • 1 using the previous notations. Following the general scheme of Section 4, we obtain

T V (u) = sup p∈U ∇u, p -F * (p) = sup p∈U -u, div p -F * (p),
where F * = δ P is the indicator function of P = {p ∈ U ; p ∞ ≤ 1}. Rewriting the TV term of [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF] in its dual formulation, the primal-dual formulation of the proposed model becomes ( 15)

   min u max p -u, div p + λ 2 Au -g 2 2 -δ P (p) s.t. ∀Ψ ∈ D, Au -g, Ψ -τ ≤ 0
We now present two algorithms to solve this optimization problem.

5.2.

Algorithm for solving the unconstrained formulation. We first derive an unconstrained formulation of ( 15), then give an algorithm to solve it. Using the method of Lagrange multipliers, we can rewrite (15) as ( 16)

min u max p,λD L (u, p, λ D ) = -u, div p + λ 2 Au -g 2 2 -δ P (p) + ψ∈D λ Ψ ( Au -g, Ψ -τ ) -δ Q (λ Ψ ) ,
where λ D = (λ Ψ ) Ψ∈D ∈ R D is the (finite dimensional) vector of Lagrange multipliers associated to the constraints Aug, Ψτ ≤ 0, and

Q = {q ∈ R; q ≥ 0}.
Obviously, the objective function in ( 16) is convex, and this is a primal-dual problem (the primal variable is u, the dual variable is (p, λ D )). The form ( 14) can be obtained with Ku = (∇u, ( Au, Ψ ) ψ∈D ), G(u) = λ 2 Aug 2 and F * (p, λ D ) = δ P (p) + Ψ∈D δ Q (λ Ψ ) + λ Ψ ( g, Ψ + τ ). Thus, applying Algorithm 1 to (16) yields the following steps:

(17)                p n+1 = arg min p -L (ū n , p, λ n D ) + 1 2tp p -p n 2 2 , λ D n+1 = arg min λD -L (ū n , p n , λ D ) + 1 2t λ λ D -λ n D 2 2 , u n+1 = arg min u L u, p n+1 , λ n+1 D + 1 2tu u -u n 2 2 . ūn+1 = 2u n+1 -u n .
Each sub-problem is solved as follows: A. Update p: This sub-optimization problem is equivalent to

min p -∇ū n , p + δ P (p) + 1 2t p p -p n 2 2 .
Let h = δ P , then the update rule for p is

p n+1 = prox tp•h (p n + t p • ∇ū n ) = p n + t p • ∇u n max(1, |p n + t p • ∇ū n |) .
B. Update λ D : Because of the separable structure of the second optimization problem in [START_REF] Chan | Simultaneous total variation image inpainting and blind deconvolution[END_REF], we can update individually each component λ Ψ of λ D by solving the corresponding minimization problem

λ n+1 Ψ = arg min λΨ -λ Ψ ( Aū n -g, Ψ -τ ) + δ Q (λ Ψ ) + 1 2t λ (λ Ψ -λ n Ψ ) 2 .
Here h = δ Q so the update rule for λ Ψ is

λ Ψ n+1 = prox t λ •h λ n Ψ + t λ • ( Aū n -g, Ψ -τ ) = max λ n Ψ + t λ • ( Aū n -g, Ψ -τ ) , 0 . C. Update u: For this sub-problem, h = λ 2 Au -g 2 2
, so the update rule for u is ( 18)

u n+1 = prox tu•h (ũ) = (I + t u λA * A) -1 (ũ + t u λA * g), where ũ = u n -t u • -div p n+1 + A * Ψ∈D λ Ψ Ψ .
u N = 1 N N n=1 u n , p N = 1 N N n=1 p n , λ D,N = 1 N N n=1
λ n D and the restricted primal-dual gap

(19) G (u, p, λ D ) = max p ′ ,λ ′ D L (u, p ′ , λ ′ D ) -min u ′ L (u ′ , p, λ D ) . If t p t u is small enough, then there exists a constant C > 0 such that for all N , G (u N , p N , λ D,N ) ≤ C N .
One can easily prove this result by following the proof of Theorem 1 in [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF]. In the present case, the operator K that permits to write [START_REF] Chan | A nonlinear primal-dual method for total variationbased image restoration[END_REF] under the general form [START_REF] Chambolle | Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage[END_REF] satisfies

Ku 2 = ∇u 2 + Ψ∈D Au, Ψ 2 , so that K 2 ≤ ∇ 2 + A 2 D 2 with D 2 = sup u ≤1 Ψ∈D u, Ψ 2 and ∇ 2 ≤ 8
(see [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]). Thus, the announced property holds as soon as t p t u (8 + A 2 D 2 ) < 1.

Algorithm 2

The proposed algorithm for the unconstrained model

Initialization: ū0 = u 0 , p = 0, λ D = 0, t p , t λ , t u . Output: u N .
Iterate for n = 0, 1, . . . N -1:

Step 1.

p n+1 = p n + t p • ∇ū n max(1, |p n + t p • ∇ū n |) Step 2. ∀Ψ ∈ D, λ n+1 Ψ = max λ n Ψ + t λ • ( Aū n -g, Ψ -τ ) , 0 Step 3. ũ = u n -t u • -div p n+1 + A * Ψ∈D λ n+1 Ψ Ψ Step 4. u n+1 = (I + t u λA * A) -1 (ũ + t u λA * g)
Step 5. ūn+1 = 2u n+1u n 5.3. Algorithm for solving the constrained formulation. In this subsection, we propose an algorithm to solve directly the constrained optimization problem (15) instead of introducing the Lagrange multipliers (λ Ψ ) Ψ∈D . Model [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF] can be written as [START_REF] Combettes | Signal recovery by proximal forward backward splitting[END_REF] min

u∈Eτ max p -u, div p + λ 2 Au -g 2 2 -δ P (p),
where

E τ = {u ∈ E; ∀Ψ ∈ D, Au -g, Ψ -τ ≤ 0}.
Minimizing the right-hand term of ( 20) with respect to u, we get the following sub-problem:

(21) min u∈Eτ -u, div p + λ 2 Au -g 2 2 .
It is equivalent to

(22) min u∈E -u, div p + λ 2 Au -g 2 2 + δ Eτ (u). With h(u) = λ 2 Au -f 2 2
, the update rule for u writes u n+1 = prox h+δE τ u nt u • -div p n+1 , and using Proposition 2, we obtain

u n+1 = P Eτ prox h u n -t u • -div p n+1 .
Unfortunately, the projection P Eτ is not easy to compute when A is not the identity matrix. This is why we propose to use this operator only for the denoising case (A = I), which yields Algorithm 3 below. Note that in the case when D is made of orthogonal vectors (e.g. an orthogonal wavelet basis) plus their opposites, P Eτ can be computed by

P Eτ (u) = g + Ψ∈D P [0,τ ] ( u -g, Ψ ) Ψ,
where P [0,τ ] (t) = min(τ, max(0, t)).

Remark 2. There are many works about the projection onto a convex set. For example, we could compute P Eτ by the method proposed in [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF], but the computational cost would be much higher than with Algorithm 2. This is why we systematically use Algorithm 2 when A = I. When A = I, it is better to use Algorithm 3 since it exactly enforces the constraints of the model at each iteration, while showing the same convergence rate as Algorithm 2 (as will be shown on Fig. 3).

Algorithm 3

The proposed algorithm for the constrained model when A = I Initialization: ū0 = u 0 , p = 0, t p , t u . Output: u N . Iterate for n = 0, 1, . . . N -1:

Step 1.

p n+1 = p n + t p • ∇ū n max (1, |p n + t p • ∇ū n |) Step 2. ũ = u n -t u • -div p n+1
Step 3. u n+1 = P Eτ ũ + t u λg 1 + t u λ

Step 4. ūn+1 = 2u n+1u n

Applications and Experiments

In this section, we use the proposed algorithms in three image applications, namely denoising, deblurring and simultaneous deblurring and inpainting. These problems can be addressed with TV regularization (see, e.g., [START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF][START_REF] Zhu | An efficient primal-dual hybrid gradient algorithm for total variation image restoration[END_REF] for deblurring, [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF] for inpainting, [START_REF] Chan | Simultaneous total variation image inpainting and blind deconvolution[END_REF] for simultaneous deblurring and inpainting), hence we will be able in each case to compare the performance of the ROF model (solved by Algorithm 1), the LM model, and our model. For the LM model and our model, we use a dictionary made of wavelet packets (see details below). Note that the ROF model has a much lower computational cost since it does not involve wavelet packet decomposition/reconstruction as in the two other models. Most of the experiments we discuss are realized with the two gray-level images (with size 256 × 256) shown in Fig. 1. Designing a dictionary for a given class of images is still an open problem. Therefore, in each case we used the dictionary obtained by the full wavelet packet decomposition up to level four using the freeware Wavelab850 1 running under MATLAB. This choice seems to provide an acceptable representation for most images. The choice of the wavelet filter also affects the performances of the model. A large compact support can achieve better performance, because a large compact support in "time" allows more narrowness in "frequency", which is usually more efficient for denoising. However, a very large compact support would introduce a heavy computational cost, so we chose the "sym8" filter to achieve a reasonable trade-off.

The peak signal-to-noise ratio (PSNR) was used to measure the quality of the recovered images. It is defined (in dB) by

PSNR(û, u) = 10 log 10 255 2 1 256×256 ij |û ij -u ij | 2 ,
where û denotes the restored image and u is the ground-truth image. As a stopping criterion we used the relative difference between two successive iterates of the 1 This software is available on http://www-stat.stanford.edu/∼wavelab restoration scheme (RDS) which is defined by

RDS u n+1 , u n = u n+1 -u n 2 u n+1 2 .
We stop the algorithm when RDS u n+1 , u n < ε. We used ε = 10 -6 for denoising, and ε = 10 -5 in all other experiments (which do not converge as fast as the denoising case). When specified, the "number of iterations" or "CPU time" of an algorithm means the number of iterations or the CPU time that was needed to drop the RDS below the error tolerance ε.

It is not easy to tune the parameters for the LM-Uzawa algorithm. A steepest gradient descent was used to compute the minimization for u in [START_REF] Lintner | Solving a variational image restoration model which involves contraints[END_REF], with the bisection method to search a suitable time step at each iteration (it usually requires about 30 sub-iterations to find the time step). We used the same algorithm, but significantly accelerated the code to avoid unnecessarily repeated computations. We observed that with this acceleration, the CPU time for one iteration was approximately the same for the LM-Uzawa and the proposed algorithms, but as we shall see later the convergence of the LM-Uzawa algorithm remains significantly slower.

The choice of ρ 0 (the initial time step for updating λ D , see [START_REF] Lintner | Solving a variational image restoration model which involves contraints[END_REF] for details) also has a strong impact on the performance of the LM-Uzawa algorithm. Large ρ 0 results in a good control of the constraints, but may also cause bad convergence behavior. The authors pointed out that one can adjust ρ 0 by ensuring that the dual function does not decrease significantly. Unfortunately, this is also not easy to control.

The parameter in T V β (u) for the LM model was set to be β = 10 -5 . In all the following experiments, we used the same constant time step settings: t u = 1, t p = 0.1, t λ = 0.1 for the proposed model, and the experiments demonstrate that this simple choice yields good performances.

6.1. Image denoising. The proposed method can be easily applied to image denoising. In this case we have A = I, A * = I and A * A = I. We give two examples: gray-level image denoising and color image denoising. Symmetric boundary conditions are used in this section.

The first example is gray-level image denoising. The original image is corrupted by an additive white Gaussian noise with zero mean and standard deviation σ = 20. We set the threshold τ = 75 for the LM model and the proposed model, ρ 0 = 20 for the LM model, λ = 0.04 for the ROF model, and λ = 0.03 for the proposed model. Fig. 2 displays the noisy image and the restored images. We can observe that the LM model and the proposed model preserve textures much better than the ROF model. As expected, this indicates that the wavelet packet approach is very helpful to recover textures. Furthermore, the proposed model converges much faster (for the required precision) than the LM model. Note that the CPU time for one iteration of LM-Uzawa algorithm is 0.98s, and that of the proposed algorithm is 0.86s, so the large speed difference is essentially due to the faster convergence rate of the proposed method (281 iterations to meet the RDS criterion, versus 1456 for the LM-Uzawa model).

We first discuss the convergence behavior of the LM model and the proposed model. We get a solution u * for the LM model by running the LM-Uzawa algorithm for 20,000 iterations, and get a solution u * for the proposed model by running Algorithm 3 for 20,000 iterations (since Algorithm 3 can exactly achieve the constraints of the proposed model at each iteration, its solution can be seen as the true solution of the proposed model). Fig. 3 gives some comparisons with respect to the CPU time. The RDS of the LM-Uzawa algorithm is not very stable, so we stopped the algorithm when the RDS was less than ε for the first time. Table 1 gives the statistics of the restored images obtained from running the algorithms for 5000 iterations. In particular, we compute the Maximum of the violated constraints

(23) Max = max Ψ∈D [ u -g, Ψ -τ ] + ,
(with as usual [t] + = max(0, t)), the Proportion of violated constraints

(24) Pro. = |V | |D| , where V = Ψ ∈ D, u -g, Ψ -τ > 10 -5 ,
and, when V = ∅, the Mean of the violated constraints

(25) Mean = 1 |V | Ψ∈V u -g, Ψ -τ.
One can clearly see from Fig. 3 and Table 1 that the proposed model is more stable than the LM model, and that the proposed algorithms (2 and 3) manage to 23); Mean value of the violated constraints only (Equation 25). Notice that after 2500s, the set of violated constraints for Algorithm 2 is empty most of the time, which explains the dotted line.

solve it precisely and much faster than the LM-Uzawa algorithm does for the LM model. Algorithm 2, which corresponds to the unconstrained formulation, exhibits TV Max Mean Pro. LM-Uzawa 7.2722 1.9385 0.1658 0.40% Algorithm 2 7.6661 0.0000 -0% Algorithm 3 7.6661 0.0000 -0%

Table 1. Experimental statistics of the result images obtained after running the algorithms for 5000 iterations (denoising case). Max: Maximum of the constraints (Equation 23). Mean: Mean of the violated constraints only (Equation 25). Pro.: proportion of violated constraints (Equation 24).

very good performances in terms of constraint enforcement: the maximum of the constraints for Algorithm 2 decreases quickly toward 0, and after 2500s, no more constraint is violated most of the time. In comparison, the LM-Uzawa algorithm clearly does not enforce the constraints as well as the proposed algorithm. Next, we discuss the choice of the parameter λ weighting the data-fidelity term in the ROF and in the proposed model. It is well-known that the performance of the ROF model is highly dependent on this parameter. Table 2 lists the PSNR values of the restored images obtained for different values of λ in the denoising case (Fig. 2) and in the Gaussian deblurring case (Fig. 8). We also present some of the corresponding resulting images in Fig. 4 for the denoising problem. From these Table and Figures, one can clearly observe that, unlike the ROF model, the performance of the proposed model is not very sensitive to the parameter λ. This indicates that the wavelet-based term can restore much information, by contributing to a trade-off with the TV regularization. 2. PSNR values (dB) of the restored images obtained with the ROF model and the proposed model for different values of λ (the reference images are shown in Fig. 1).

To illustrate further the general behavior of the proposed model, we perform experiments on four other test images, shown in Fig. 5: Boat (size 512 × 512), Cameraman (256×256), Goldhill (512×512), and Lena (512×512). The degradation process is the same as before. Table 3 presents the resulting PSNR values. We can observe that the LM model and our model do not reach much better results than the ROF model for these relatively poorly textured images. In fact, the LM model (and consequently ours), with the standard dictionary used, are in favor of textured images. Note incidentally that the proposed model achieves slightly higher PSNR values than the LM model. To further improve the performance of the LM model and our model on images with less texture, one may consider choosing a better dictionary, as indicated in [START_REF] Zeng | On the total variation dictionary model[END_REF], or using more powerful dictionary learning method [START_REF] Liu | HuanA Weighted Dictionary Learning Model for Denoising Images Corrupted by Mixed Noise[END_REF].

For texture images, the proposed model is better than ROF model. We perform experiments on two Brodatz texture images shown in Fig. 6. The degradation The proposed model can also be easily generalized to color image denoising. We adopt the scheme proposed by Bresson and Chan [START_REF] Bresson | Fast dual minimization of the vectorial total variation norm and applications to color image processing[END_REF] to compute the dual formulation of the vectorial TV regularization. The top-left image in Fig. 7 Table 3. The resulting PSNR(dB) on the test images of Fig. 5. color image with resolution 256 × 256 × 3 (3 is the number of color channels). The noisy image is obtained by adding independent noise (with the same Gaussian distribution as before) on each color channel. On Fig. 7, we compare the proposed model with the vectorial ROF model, solved by Bresson and Chan [START_REF] Bresson | Fast dual minimization of the vectorial total variation norm and applications to color image processing[END_REF] 2 via a dual minimization algorithm. We can observe that the proposed method clearly achieves a better denoising, which is confirmed by the higher PSNR value.

6.2. Image deblurring. In this section, we apply the proposed model to a more difficult task: image deblurring. We will consider three blur kernels: a Gaussian blur kernel, a uniform blur kernel and a motion blur kernel. In these experiments, periodic boundary conditions are adopted and the FFT is used to compute the convolution operator. The initial images of the algorithms are the blurred and noisy images. We set τ = 5 for the first two experiments, τ = 2 for the last experiment since the degraded image has less noise, and ρ 0 = 5 for the LM model in all the experiments. In the first experiment, the corrupted image is obtained by convolving the original image with a Gaussian blur of size 5 × 5 with mean zero and standard deviation 2 (produced by the MATLAB function "fspecial") and then adding a Gaussian white noise with mean zero and standard deviation 2. We set λ = 5 for the ROF model and the proposed model. Fig. 8 displays the restored images. 2 The MATLAB code is available on http://www.cs.cityu.edu.hk/∼xbresson/ucla/code.html Then, we choose a uniform 7 × 7 blur to corrupt the original image, and add as before a Gaussian white noise with mean zero and standard deviation 2 to it. Again, the parameter λ is set to 5. The restored images are shown in Fig. 9.

Finally, we apply the proposed algorithm to a motion-blurred image. The kernel is generated by the MATLAB function "fspecial" with approximately 30 pixels and an angle of 45 • . Then the blurred image is corrupted with an additive white Gaussian noise with mean zero and standard deviation 1. With the parameter λ = 20, the restored images are shown in Fig. 10.

Fig. 11 shows the resulting RDS of the three experiments mentioned above. One can observe that the proposed algorithm achieves better performance than the LM-Uzawa algorithm. Table 4 gives some statistics of the restored images obtained after running the algorithms for 5000 iterations. The Max, Mean and Pro. statistics correspond to the definitions given in Equations ( 23), ( 25) and [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF], except that u is replaced with Au to take into account the blur operator A. Although the maximum of the violated constraints is a little larger than that of the LM-Uzawa algorithm in the case of a uniform blur kernel, the proportion of violated constraints is always much smaller for the proposed model. 23) to [START_REF] Ekeland | Convex Analysis and Variational Problems[END_REF] with u replaced by Au).

6.3. Simultaneous image deblurring and inpainting. In several image processing applications, some entire image regions are damaged or lost. The goal of image inpainting is to restore the missing and/or damaged areas using the valid surrounding information. Since the pioneering work of [START_REF] Masnou | Level lines based disocclusion[END_REF] and [START_REF] Bertalmio | Image inpainting[END_REF], this problem has been extensively studied [START_REF] Bugeau | A Comprehensive Framework for Image Inpainting[END_REF]. In this subsection, we consider a more complicated image degradation problem: images are blurred and pixels are missing in some regions.

Let Ω 1 denote the set of indices corresponding to missing pixels. The linear transform operator for inpainting is defined as

(A 2 u) i,j = u i,j , if u i,j / ∈ Ω 1 , 0, if u i,j ∈ Ω 1 .
Let A 1 denote the convolution operator. Then the degraded operator considered in this subsection is A = A 2 • A 1 . In this case, because the operator A * 1 • A 2 • A 1 can not be diagonalized by the FFT, the update step 4 in Algorithm 2 can not be solved as above. Fortunately, thanks to the I term, the matrix I +t u λA * A is always invertible. Thus we can use the preconditioned conjugate gradient method to solve this sub-problem.

We study two situations for missing pixels: randomly missing pixels and structured missing regions. In these experiments, we use the same parameter values as in the motion deblurring experiments, except for the randomly missing pixels case for which we take τ = 1.2 (because there are lots of missing pixels, the coefficients of the wavelet packet decomposition are smaller). To accelerate the convergence of the algorithms, we use a initial image obtained by performing the nearest neighbor interpolation on the corrupted image.

In the first experiment, the corrupted image is obtained by randomly removing 50% pixels on the motion-blurred image used in the previous section. On Fig. 12, we can see that the visual quality of the image restored by the proposed model is better. In the second experiment, we still use the same motion-blurred image but the missing pixels now correspond to a superimposed text. The results are displayed on Fig. 13. Table 5 gives some statistics of the restored images obtained by running the algorithms for 5000 iterations, and Fig. 14 shows the resulting RDS. One can observe that the proposed model still exhibits better performances on the simultaneous deblurring and inpainting problem.

Combining nonlocal total variation (NLTV)

Recently, much literature has focused on exploiting nonlocal similarities of textures and has reached higher restoration quality than methods only based on the similarity between pixels. Buades et al. [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] proposed the nonlocal means algorithm, which exhibits good performances for image denoising. Gilboa and Osher [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF] proposed the NLTV approach, based on a variational framework. Zhang et al. [START_REF] Zhang | Bregmanized Nonlocal Regularization for Deconvolution and Sparse Reconstruction[END_REF] used this regularization for image deconvolution and sparse reconstruction. We use the definitions of the NL functionals presented in [START_REF] Bresson | A short note for nonlocal TV minimization[END_REF][START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF][START_REF] Zhang | Bregmanized Nonlocal Regularization for Deconvolution and Sparse Reconstruction[END_REF]. The general NLTV framework for image restoration can be defined as [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF] min

u ∇ w u 1 + λ 2 Au -g 2 2 ,
where (∇ w u) (x, y) = (u(y)u(x)) w (x, y). The weight function w is a symmetric function that can be defined from the initial image g and h as w(x, y) = exp -1 h 2 Ω G a * g(x + z)g(y + z)

2 2 dz ,
where Ω G a * g(x + z)g(y + z)

2
2 dz is the distance between patches located at x and y on the given degraded image, G a is a Gaussian function with standard deviation a and h is a filtering parameter that has to be set according to the noise level. The divergence of the dual variable is defined by (div w p) (x) = Ω (p(x, y)p(y, x)) w(x, y)dy, and we have ∇ w u, p =u, div w p . We can easily replace TV regularization by the NLTV regularization in the proposed algorithm. In the following, we give some experimental results of the proposed algorithm and the NLTV restoration model which can be solved by Algorithm 1.

To compute the weights (the values of the all pixels of the degraded image are scaled in the range [0, 1] when computing the weights), the MATLAB code provided by Bresson3 was used. We set the patch size to 9 × 9, the searching window size to 11 × 11 and the number of neighbors to 10. The filtering parameter was set to h = 0.45 for image denoising, h = 0.01 for restoration from Gaussian blur and uniform blur, and h = 0.005 for motion deblurring. Since NLTV regularization uses more neighbors than TV regularization, the parameter λ has to be increased, so we took λ = 0.05 for denoising, λ = 10 for restoration from Gaussian blur and uniform blur, and λ = 30 for motion deblurring. We included the 4 nearest neighbors in

Conclusion

We presented a modified model for image restoration based on Total Variation and wavelet packets. In this model, the original TV regularization term is replaced by its dual formulation, and a first-order primal-dual algorithm is used to solve the corresponding optimization problems. The performance of the proposed model is demonstrated by experiments on several standard image restoration problems (denoising, deblurring and inpainting), and the results show that it manages to recover textures better than previous comparable approaches. To further improve the performance of the proposed model, we then replaced TV regularization with NLTV regularization, which resulted in a clear improvement of the quality of the restored images. Future work could focus on the choice of an adaptive dictionary to increase further the performance of the proposed model.

  and F : U → [0, +∞) and G : E → [0, +∞) are proper, convex, lower semicontinuous functions. The convex conjugate F * of F [6] is defined by ∀p ∈ U, F * (p) = sup w∈U p, w -F (w). Since F * * = F , we have ∀u ∈ E, F (Ku) = sup p∈U p, Ku -F * (p).

Figure 1 .

 1 Figure 1. Original images used for image restoration. The left image is used for denoising, and the right one is used for deblurring and inpainting. In both cases (as well as for the other images considered later in Fig. 5), the pixel values are in the range [0, 255].

Figure 2 .

 2 Figure 2. Denoising of a gray-level image. Top, left: noisy image; Top, right: ROF Model (787 iterations, CPU time = 3.9s, PSNR = 24.75dB); Bottom, left: LM Model (1456 iterations, CPU time = 1434s, PSNR = 26.09dB); Bottom, right: our model (281 iterations, CPU time = 241s, PSNR = 26.19dB).

Figure 3 .

 3 Figure 3. Evolution in function of the CPU time (in seconds) of several quantities during gray-level denoising with ROF, LM-Uzawa, and the proposed Algorithms 2 and 3 (these two curves often overlap). From left to right, then top to bottom: PSNR; distance to solution u nu * 2 (log scale); TV; RDS (log scale); Maximum value of the constraints (Equation23); Mean value of the violated constraints only (Equation25). Notice that after 2500s, the set of violated constraints for Algorithm 2 is empty most of the time, which explains the dotted line.

Figure 4 .

 4 Figure 4. Image denoising with different values of λ (top line: λ = 0.02, bottom line: λ = 0.05). Left column: ROF Model; Right column: our model.

Figure 5 .

 5 Figure 5. Test images. From left to right: Boat, Cameraman, Goldhill, Lena.

Figure 6 .

 6 Figure 6. Original Brodatz texture images.

Figure 7 .

 7 Figure 7. Color image denoising. Top, left: original color image; Top, right: noisy image; Bottom, left: vectorial ROF model (480 iterations, CPU time=10.7s, PSNR=26.54dB); Bottom, right: Proposed model (238 iterations, CPU time=626.2s, PSNR=26.92dB).

Figure 8 .

 8 Figure 8. Restoration from Gaussian blur. Top, left: blurred image; Top, right: ROF model (452 iterations, CPU time=4.0s, PSNR=25.52dB); Bottom, left: LM model (1602 iterations, CPU time=1639s, PSNR=27.16dB); Bottom, right: Proposed model (966 iterations, CPU time=861s, PSNR= 27.43dB).

Figure 9 .

 9 Figure 9. Restoration from uniform blur. Top, left: blurred image; Top, right: ROF model (522 iterations, CPU time=4.6s, PSNR=24.73dB); Bottom, left: LM model (1978 iterations, CPU time=1991s, PSNR=25.49dB); Bottom, right: Proposed model (1083 iterations, CPU time=968s, PSNR= 25.75dB).

Figure 10 .

 10 Figure 10. Restoration from motion blur. Top, left: blurred image; Top, right: ROF model (806 iterations, CPU time=8.1s, PSNR=25.58dB); Bottom, left: LM model (3348 iterations, CPU time=3402s, PSNR=26.21dB); Bottom, right: Proposed model (1546 iterations, CPU time=1382s, PSNR= 27.14dB).

Figure 11 .

 11 Figure 11. RDS in function of CPU time (in seconds) for the deblurring case. From left to right: Gaussian blur, uniform blur, motion blur.

Figure 12 .

 12 Figure 12. Image inpainting with randomly missing pixels. Top, left: corrupted image; Top, right: ROF model (993 iterations, CPU time=142s, PSNR=24.29dB); Bottom, left: LM model (2449 iterations, CPU time=2477s, PSNR=24.63dB); Bottom, right: Proposed model (1856 iterations, CPU time=2084s, PSNR=25.06dB).

Figure 13 .

 13 Figure 13. Image inpainting with structured missing pixels (superimposed text). Top, left: corrupted image; Top, right: ROF model (925 iterations, CPU time=151s, PSNR=25.06dB); Bottom, left: LM model (2721 iterations, CPU time=2803s, PSNR=25.30dB); Bottom, right: Proposed model (1388 iterations, CPU time=1546s, PSNR=26.20dB).

Figure 14 .

 14 Figure 14. RDS in function of CPU time (in seconds) for the inpainting case. Left: randomly missing pixels; Right: structured (text) missing regions.

Figure 15 .

 15 Figure 15. Image denoising based on NLTV regularization. Left: NLTV restoration model (PSNR=28.11dB); Right: Proposed model (PSNR=28.68dB).

Figure 16 .

 16 Figure 16. Results of image deblurring based on NLTV regularization. Left column: classical NLTV deblurring; Right column: proposed model combined with NLTV. From top to bottom: Gaussian blur (PSNR=26.44dB, 28.11dB), uniform blur (PSNR=25.48dB, 26.16dB), and motion blur (PSNR=25.79dB, 27.50dB).

  We use a multi-step procedure to accelerate the convergence, which results in Algorithm 2 below.The following theorem gives the convergence rate of Algorithm 2.

Theorem 5.1. Choose t u > 0, t λ > 0 and t p > 0 such that t λ = t p . Consider the sequences {u n }, {p n } and {λ n D } generated by Algorithm 2. Let

  .[START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF] 24.75 25.38 24.68 24.92 25.22 25.52 Proposed model 25.69 26.19 26.57 26.85 27.10 27.27 27.37 27.43

			denoising			deblurring
	λ	0.02	0.03	0.04	0.05	2	3	4	5
	ROF model 23.59 24Table

  is the original .[START_REF] Zhang | Bregmanized Nonlocal Regularization for Deconvolution and Sparse Reconstruction[END_REF] 28.72 28.58 28.76 28.66 28.81 Cameraman 28.21 27.99 28.36 28.05 28.41 28.12 28.46 Goldhill 29.06 29.03 29.16 29.08 29.19 29.15 29.23 Lena 31.04 31.03 31.05 31.06 31.06 31.11 31.06 Average 29.23 29.14 29.32 29.19 29.36 29.26 29.39

			τ = 75	τ = 70	τ = 65
	Images	ROF	LM Ours	LM Ours	LM Ours
	Boat	28.61 28

Table 4 .

 4 ). Experimental statistics obtained for the deblurring case after 5000 iterations: Total Variation, maximum of the violated constraints, proportion of violated constraints, mean of the violated constraints (see Equations (

		Algorithms	TV	Max	Pro.	Mean
	Gaussian blur LM-Uzawa 10.0958 7.9479 3.93% 0.0038
	τ = 5	Algorithm 2 11.2513 7.9154 0.90% 0.0027
	Uniform blur		8.8571 8.0361 2.27% 0.0031
	τ = 5	Algorithm 2 9.3452 8.0533 0.96% 0.0029
	Motion blur	LM-Uzawa 11.7072 3.8542 5.79% 0.0038
	τ = 2	Algorithm 2 12.7043 3.0598 2.50% 0.0024

Table 5 .

 5 ). Experimental statistics of the resulting images obtained after running the algorithms for 5000 iterations (inpainting case).

		Algorithms	TV	Max	Pro.	Mean
	random pixels inpainting LM-Uzawa 10.3346 1.4823 7.13% 0.0031
	τ = 1.2	Algorithm 2 10.3569 1.5168 5.17% 0.0023
	text inpainting	LM-Uzawa 10.0482 2.7645 4.57% 0.0025
	τ = 2	Algorithm 2 11.1518 2.9219 1.79% 0.0013

http://www.cs.cityu.edu.hk/∼xbresson/ucla/code.html

Acknowledgments

We would like to thank the anonymous referees for their useful comments and suggestions.

This research was supported in part by NSFC 11271049, 61033013, RGC 211710, RGC 211911, the National 863 Project (2012AA040912), the Ph.D. Programs Foundation of Ministry of Education of China (20120009110006), and the FRGs of Hong Kong Baptist University.