Nelly Pustelnik
email: nelly.pustelnik@univ-paris-est.fr

Caroline Chaux
email: caroline.chaux@univ-paris-est.fr

Jean-Christophe Pesquet
email: jean-christophe.pesquet@univ-paris-est.fr

† N Pustelnik

Parallel Proximal Algorithm for Image Restoration Using Hybrid Regularization -Extended Version *

Regularization approaches have demonstrated their effectiveness for solving ill-posed problems. However, in the context of variational restoration methods, a challenging question remains, namely how to find a good regularizer. While total variation introduces staircase effects, wavelet domain regularization brings other artefacts, e.g. ringing. However, a trade-off can be made by introducing a hybrid regularization including several terms non necessarily acting in the same domain (e.g. spatial and wavelet transform domains). While this approach was shown to provide good results for solving deconvolution problems in the presence of additive Gaussian noise, an important issue is to efficiently deal with this hybrid regularization for more general noise models. To solve this problem, we adopt a convex optimization framework where the criterion to be minimized is split in the sum of more than two terms. For spatial domain regularization, isotropic or anisotropic total variation definitions using various gradient filters are considered. An accelerated version of the Parallel Proximal Algorithm is proposed to perform the minimization. Some difficulties in the computation of the proximity operators involved in this algorithm are also addressed in this paper. Numerical experiments performed in the context of Poisson data recovery, show the good behaviour of the algorithm as well as promising results concerning the use of hybrid regularization techniques.

Introduction

During the last decades, convex optimization methods have been shown to be very effective for solving inverse problems. On the one hand, algorithms such as Projection Onto Convex Sets (POCS) [START_REF] Bregman | The method of successive projection for a common point of convex sets[END_REF][START_REF] Gurin | Projection methods for finding a common point of convex sets[END_REF][START_REF] Youla | Image restoration by the method of convex projections. Part Itheory[END_REF][START_REF] Combettes | The foundations of set theoretic estimation[END_REF][START_REF] Combettes | The Convex Feasibility Problem in Image Recovery[END_REF] have become popular for finding a solution in the intersection of convex sets. POCS was used in data recovery problems [START_REF] Trussell | The feasible solution in signal restoration[END_REF] in order to incorporate prior information on the target image (e.g. smoothness constraints). Some variants of POCS such as ART (Algebraic Reconstruction Technique) [START_REF] Gordon | Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography[END_REF] or PPM (Parallel Projection Method) [START_REF] Iusem | Convergence results for an accelerated nonlinear Cimmino algorithm[END_REF][START_REF] Combettes | Inconsistent signal feasibility problems : least-squares solutions in a product space[END_REF] were also proposed to achieve iteration parallelization. Additional variants of POCS can be found in [START_REF] Combettes | Convex set theoretic image recovery by extrapolated iterations of parallel subgradient projections[END_REF]. Other parallel approaches such as block-iterative surrogate constraint splitting methods were considered to solve a quadratic minimization problem under convex constraints [START_REF] Combettes | A block-iterative surrogate constraint splitting method for quadratic signal recovery[END_REF] which may include a total variation constraint (see also [START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF]). However the method in [START_REF] Combettes | A block-iterative surrogate constraint splitting method for quadratic signal recovery[END_REF] based on subgradient projections is not applicable to non-differentiable objective functions.

On the other hand, some denoising approaches were based on wavelet transforms [START_REF] Mallat | A wavelet tour of signal processing[END_REF], and more generally on frame representations [START_REF] Candès | Recovering edges in ill-posed inverse problems: Optimality of curvelet frames[END_REF][START_REF] Pennec | Sparse geometric image representations with bandelets[END_REF][START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF][START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF]. In [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Figueiredo | An EM algorithm for wavelet-based image restoration[END_REF][START_REF] Bect | A l 1 -unified variational framework for image restoration[END_REF][START_REF] Harmany | This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms -theory and practice[END_REF], algorithms which belong to the class of forward-backward algorithms were proposed in order to restore images degraded by a linear operator and a noise perturbation. Forward-backward iterations allow us to minimize a sum of two functions assumed to be in the class Γ 0 (H) of lower semicontinuous convex functions defined on a Hilbert space H, and taking their values in ] -∞, +∞], which are proper (i.e. not identically equal to +∞). In addition, one of these functions must be Lipschitz differentiable on H. In [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], this algorithm was investigated by making use of proximity operator tools [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] firstly proposed by Moreau in [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]. In [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF], applications to frame representations were developed and a list of closed form expressions of several proximity operators was provided. Typically, forwardbackward methods are appropriate when dealing with a smooth data fidelity term e.g. a quadratic function and a non-smooth penalty term such as an ℓ 1 -norm promoting sparsity in the considered frame [START_REF] Bioucas-Dias | A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. The computation of the proximity operator associated with the ℓ 1 -norm indeed reduces to a componentwise soft-thresholding [START_REF] Donoho | De-noising by soft-thresholding[END_REF][START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]. Another optimization method known as the Douglas-Rachford algorithm [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF][START_REF] Douglas | On the numerical solution of the heat conduction problem in two and three space variables[END_REF][START_REF] Eckstein | On the Douglas-Rachford splitting methods and the proximal point algorithm for maximal monotone operators[END_REF][START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] was then proposed for signal/image recovery problems [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] to relax the Lipschitz differentiablity condition required in forward-backward iterations. In turn, the latter algorithm requires the knowledge of the proximity operators of both functions. This algorithm was then extended to the minimization of a sum of a finite number of convex functions [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF], the proximity operator of each function still being assumed to be known. One of the main advantages of this algorithm called Parallel ProXimal Algorithm (PPXA) is its parallel structure which makes it easily implementable on multicore architectures. PPXA is well suited to deconvolution problems in the presence of additive Gaussian noise, where the proximity operator associated with the fidelity term takes a closed form [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]. To minimize a sum of two functions one of which is quadratic, another interesting class of parallel convex optimization algorithms was proposed by Fornasier et al. in [START_REF] Fornasier | Domain decomposition methods for linear inverse problems with sparsity constraints[END_REF][START_REF] Fornasier | Subspace correction methods for total variation and ℓ 1 -minimization[END_REF]. However, in a more general context, particularly when the noise is not additive Gaussian and a wider class of degradation operators is considered, the proximity operator associated with data fidelity term does not have a closed form, which prevents the direct use of PPXA and the algorithms in [START_REF] Fornasier | Domain decomposition methods for linear inverse problems with sparsity constraints[END_REF] and [START_REF] Fornasier | Subspace correction methods for total variation and ℓ 1 -minimization[END_REF]. Therefore, other solutions have to be looked for. For Poisson noise, a first solution is to resort to the Anscombe transform [START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF], while a second one consists of approximating the Poisson data fidelity term with a gradient Lipschitz function [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF]. Both approaches require the use of a nested iterative algorithm [START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF][START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF], combining forward-backward and Douglas-Rachford iterations. Nested algorithms may however appear limited for two main reasons: the parallelization of the related iterations is difficult, and the number of functions to be minimized is in practice limited to three. More recently, approaches related to augmented Lagrangian techniques [START_REF] Hestenes | Multiplier and gradient methods[END_REF][START_REF] Fortin | Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems[END_REF] have been considered in [START_REF] Goldstein | The split Bregman method for ℓ 1 regularized problems[END_REF][START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF][START_REF] Figueiredo | Restoration of Poissonian images using alternating direction optimization[END_REF][START_REF] He | A splitting method for separate convex programming with linking linear constraints[END_REF]. These methods are well-adapted when the linear operator is a convolution and Fourier diagonalization techniques can thus be used, but for more general linear degradation operators, a large-size linear system of equations has to be solved numerically at each iteration of the algorithm.

The objective of this paper is to propose an adaptation of PPXA to minimize criteria used in a wide panel of restoration problems such as those involving a convolution or decimated convolution operator using a finite-support kernel and non-necessarily additive Gaussian noise. Decimated convolutions are important in practice since they are often encountered in super-resolution problems. To apply proximal methods, it seems that we should be able to compute the proximity operator associated with the fidelity term for a large class of noise distributions. When the proximity operator cannot be easily computed, we will show that a splitting approach may often be employed to circumvent this difficulty. This is the main contribution of this paper.

Moreover, based on similar splitting techniques and in the spirit of existing works on deconvolution in the presence of Gaussian noise [START_REF] Bect | A l 1 -unified variational framework for image restoration[END_REF][START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF][START_REF] Bioucas-Dias | An iterative algorithm for linear inverse problems with compound regularizers[END_REF][START_REF] Wen | Iterative algorithms based on decoupling of deblurring and denoising for image restoration[END_REF], a twofold regularization composed of a sparsity term and a total variation term is performed in order to benefit from each regularization. We will consider this type of hybrid regularization by investigating different discrete forms of the total variation. The paper is organized as follows: first, in Section 2, we present the considered restoration problem and the general form of the associated criterion to be minimized. Then, in Section 3, the definition and some properties of proximity operators as well as explicit forms related to the data fidelity term in a restoration context and to a discretization of the total variation are provided. Section 4 introduces an accelerated version of PPXA which allows us to efficiently solve framebased image recovery problems. Section 5.1 shows how the results obtained in the two previous sections can be used for solving restoration problems where a regularization is performed both in the spatial and in the wavelet domains. Finally, in Section 5.2, the effectiveness of the proposed approach is demonstrated by experiments for the restoration of images degraded by a blur (or a decimated blur) with finite-support kernel and by a Poisson noise. Some conclusions are drawn in Section 6.

Background

Image restoration

The degradation model considered throughout this paper is the following:

z = D α (T y) (1)
where y denotes the original image of global size N degraded by a non-negative valued convolutive operator T : R N → R M and contaminated by a noise non necessarily additive, the effect of which is denoted by D α . Here, α is a positive parameter which characterizes the noise intensity. The vector z ∈ R M represents the observed data of size M . For example, D α may denote the addition of a zero-mean Laplacian noise with standard-deviation α, or the corruption by an independent Poisson noise with scaling parameter α. T represents a convolution or a decimated convolution operator using a finite-support kernel.

Our objective is to recover the image y from the observation z by using some prior information on its frame coefficients and its spatial properties.

Frame representation

In inverse problems, certain physical properties of the target solution y are most suitably expressed in terms of the coefficients x = (x (k) ) 1≤k≤K ∈ R K of its representation y = K k=1 x (k) e k with respect to a family of vectors (e k ) 1≤k≤K in the Euclidean space R N . Recall that a family of vectors (e k ) 1≤k≤K in R N constitutes a frame if there exist two constants ν and ν in ]0, +∞[ such that1

(∀y ∈ R N ) ν y 2 ≤ K k=1 |e k ⊤ y| 2 ≤ ν y 2 . ( 2
)
The associated frame operator is the injective linear operator

F : R N → R K : y → (e k ⊤ y) 1≤k≤K ,
the adjoint of which is the surjective linear operator [START_REF] Bregman | The method of successive projection for a common point of convex sets[END_REF], (e k ) 1≤k≤K is said to be a tight frame. In this case, we have

F ⊤ : R K → R N : (x (k) ) 1≤k≤K → K k=1 x (k) e k . When ν = ν = ν in
F ⊤ F = νId, (3)
where Id denotes the identity matrix. A simple example of a tight frame is the union of ν orthonormal bases, in which case ν = ν = ν. For instance, a 2D real (resp. complex) dual-tree wavelet decomposition is the union of two (resp. four) orthonormal wavelet bases [START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF]. Curvelets [START_REF] Candès | Recovering edges in ill-posed inverse problems: Optimality of curvelet frames[END_REF] constitute another example of tight frame. Historically, Gabor frames [START_REF] Daubechies | Ten lectures on wavelets[END_REF] have played an important role in many inverse problems. Under some conditions, contourlets [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF] also constitute tight frames. When F -1 = F ⊤ , an orthonormal basis is obtained. Further constructions as well as a detailed account of frame theory in Hilbert spaces can be found in [START_REF] Han | Frames, bases, and group representations[END_REF].

In such a framework, the observation model becomes

z = D α (T F ⊤ x) (4)
where x represents the frame coefficients of the original data (y = F ⊤ x ∈ R N is the target data of size N ). Our objective is now to recover x from the observation z.

Minimization problem

In the context of inverse problems, the original image can be restored by solving a convex optimization problem of the form:

Find x ∈ Argmin x∈R K J j=1 f j (x) (5)
where (f j ) 1≤j≤J are functions of Γ 0 (R K ) (see [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] and references therein) and the restored image is y = F ⊤ x.

A particular popular case is when J = 2; the minimization problem thus reduces to the minimization of the sum of two functions which, under a Bayesian framework, can be interpreted as a fidelity term f 1 linked to noise and an a priori term f 2 related to some prior probabilistic model put on the frame coefficients (some examples will be given in Section 5).

In this paper, we are especially interested in the case when J > 2, which may be fruitful for imposing additional constraints on the target solution. At the same time, when considering a frame representation (which, as already mentioned, often allows us to better express some properties of the target solution), the convex optimization problem (5) can be re-expressed as:

Find x ∈ Argmin x∈R K S j=1 g j (F ⊤ x) + J j=S+1 f j (x) (6)
where (g j ) 1≤j≤S are functions of Γ 0 (R N ) and (f j ) S+1≤j≤J are functions of Γ 0 (R K ), related to the image or to the frame coefficients, respectively. The terms for j ∈ {1, . . . , S} related directly to the pixel values may be the data fidelity term, or a pixel range constraint term, whereas, the functions of indices j ∈ {S + 1, . . . , J} defined on frame coefficients are often chosen from some classical prior probabilistic model. For example, they may correspond to the minus log-likelihood of independent variables following generalized Gaussian distributions [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF].

We will now present convex analysis tools which are useful to deal with such minimization problems.

3 Proximal tools

Definition and examples

A fundamental tool which has been widely employed in the recent convex optimization literature is the proximity operator [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] first introduced by Moreau in 1962 [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]. The proximity operator of ϕ ∈ Γ 0 (R X ) is defined as

prox ϕ : R X → R X : u → arg min v∈R X 1 2 v -u 2 + ϕ(v). (7)
Thus, if C is a nonempty closed convex set of R X , and ι C denotes the indicator function of C, i.e., ∀u ∈ R X , ι C (u) = 0 if u ∈ C, +∞ otherwise, then, prox ι C reduces to the projection P C onto C. Other examples of proximity operators corresponding to the potential functions of standard log-concave univariate probability densities have been listed in [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF][START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF][START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]. Some of them will be used in the paper and we will thus recall the proximity operators of the potentials associated with a Gamma distribution (which is closely related to the Kullback-Leibler divergence [START_REF] Titterington | On the iterative image space reconstruction algorithm for ECT[END_REF]) and with a generalized Gaussian distribution, before dealing with the Euclidean norm in dimension 2.

Example 3.1 [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] Let α > 0 and set

ϕ : R → ]-∞, +∞] η →      -χ ln(η) + αη, if χ > 0 and η > 0; αη, if χ = 0 and η ≥ 0; +∞, otherwise. (8)
Then, for every η ∈ R,

prox ϕ η = η -α + |η -α| 2 + 4χ 2 . ( 9
) Example 3.2 [26] Let χ > 0, p ∈ [1, +∞[, and set ϕ : R → ]-∞, +∞] : η → χ|η| p . (10)
Then, for every η ∈ R, prox ϕ η is given by

                       sign(η) max{|η| -χ, 0} if p = 1 η + 4χ 3 . 2 1/3 (ǫ -η) 1/3 -(ǫ + η) 1/3 where ǫ = η 2 + 256χ 3 /729 if p = 4 3 η + 9χ 2 sign(η) 8 1 -1 + 16|η| 9χ 2 if p = 3 2 η 1+2χ if p = 2 sign(η) √ 1+12χ|η|-1 6χ if p = 3 ( 11
)
where sign denotes the signum function. In Example 3.2, it can be noticed that the proximity operator associated with p = 1 reduces to a soft thresholding.

Example 3.3 [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] Let µ > 0 and set

ϕ : R 2 → R : (η 1 , η 2 ) → µ |η 1 | 2 + |η 2 | 2 . ( 12
)
Then, for every (η

1 , η 2 ) ∈ R 2 , prox ϕ (η 1 , η 2 ) =    1 - µ √ |η 1 | 2 +|η 2 | 2 (η 1 , η 2 ), if |η 1 | 2 + |η 2 | 2 > µ;
(0, 0), otherwise.

(13)

Proximity operators involving a linear operator

We will now study the problem of determining the proximity operator of a function g = Ψ • T where T : R N → R M is a linear operator,

Ψ: R M → ]-∞, +∞] : (u (m) ) 1≤m≤M → M m=1 ψ m (u (m) ) (14)
and, for every m ∈ {1, . . . , M }, ψ m ∈ Γ 0 (R). As will be shown next, the proximity operator of this function can be determined in a closed form for specific cases only. However, g can be decomposed as a sum of functions for which the proximity operators can be calculated explicitly. Firstly, we introduce a property concerning the determination of the proximity operator of the composition of a convex function and a linear operator, which constitutes a generalization of [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF]Proposition 11] for separable convex functions. The proof of the following proposition is provided in Appendix 7.

Proposition 3.4 Let X ∈ N * , Y ∈ N * , and let (o m ) 1≤m≤Y be an orthonormal basis of R Y . Let Υ be a function such that

(∀u ∈ R Y ) Υ(u) = Y m=1 ψ m (o m ⊤ u) (15)
where

(ψ m ) 1≤m≤Y are functions in Γ 0 (R). Let L be a matrix in R Y ×X such that LL ⊤ D = Y m=1 ∆ m o m o ⊤ m ( 16
)
where (∆ m ) 1≤m≤Y is a sequence of positive reals.

Then Υ • L ∈ Γ 0 (R X ) and, for every v ∈ R X prox Υ•L v = v + L ⊤ D -1 prox DΥ (Lv) -Lv ( 17
)
where DΥ is the function defined by

(∀u ∈ R Y ) DΥ(u) = Y m=1 ∆ m ψ m o m ⊤ u . (18)
The function Ψ defined in ( 14) is separable in the canonical basis of R M . However, for an abitrary convolutive (or decimated convolutive) operator L = T , ( 16) is generally not satisfied. Nevertheless, assume that (I i ) 1≤i≤I is a partition of {1, . . . , M } in nonempty sets. For every i ∈ {1, . . . , I}, let M i be the number of elements in

I i ( I i=1 M i = M ) and let Υ i : R M i → ]0, +∞[ : (u (m) ) m∈I i → m∈I i ψ m (u (m)
). If, for every i ∈ {1, . . . , M }, t i is the vector of R N corresponding to the i-th row vector of T , we have then

g = I i=1 Υ i • T i where T i is a linear operator from R N to R M i associated with a matrix t m 1 . . . t m M i ⊤ ( 19
)
and

I i = {m 1 , . . . , m M i }.
The following assumption will play a prominent role in the rest of the paper: Assumption 3.5 For every i ∈ {1, . . . , I}, (t m ) m∈I i is a family of non zero orthogonal vectors.

Then, g can be decomposed as a sum of I functions (Υ i • T i ) 1≤i≤I where, for every i ∈ {1, . . . , I}, D i = T i T ⊤ i is associated with an invertible diagonal matrix Diag(∆ i,1 , . . . , ∆ i,M i ). According to Proposition 3.4, we have then, for every

y ∈ R N , prox Υ i •T i y = y + T ⊤ i D -1 i prox D i Υ i (T i y) -T i y . ( 20
)
Remark 3.6 1. Note that Assumption 3.5 is obviously satisfied when I = M , that is when, for every i ∈ {1, . . . , I}, I i reduces to a singleton.

It can be noticed that the application of T i or T ⊤

i reduces to standard operations in signal processing. For example, when T corresponds to a convolutive operator, the application of T i consists of two steps: a convolution with the impulse response of the degradation filter and a decimation for selected locations (m ∈ I i ). The application of T ⊤ i also consists of two steps: an interpolation step (by inserting zeros between data values of indices m ∈ I i ) followed by a convolution with the filter with conjugate frequency response.

The fundamental idea behind the previously introduced partition (I i ) 1≤i≤I , is to form groups of non-overlapping -and thus orthogonal -shifts of the convolution kernel so as to be able to compute the corresponding proximity operators. To reduce the number of proximity operators to be computed, one usually wants to find the smallest integer I such that, for every i ∈ {1, . . . , I}, (t m ) m∈I i is an orthogonal family. For the sake of simplicity, we will consider the case of a 1D deconvolution problem, where N represents the original signal size whereas M corresponds to the degraded signal size, the extension to 2D deconvolution problems being straightforward. Different configurations concerning the impact of boundary effects on the convolution operator will be studied: first, we will consider the case when no boundary effect occurs. Then, boundary effects introduced by zero padding and by a periodic convolution will be taken into account. Finally, the special case of decimated convolution will be considered. Q designates in the sequel the length of the kernel and (θ q ) 0≤q<Q its values.

1. One-dimensional convolutive models without boundary effect.

We typically have the following Toeplitz structure:

   t ⊤ 1 . . . t ⊤ M    =       θ Q-1 . . . θ 1 θ 0 0 . . . 0 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0 . . . 0 θ Q-1 . . . θ 1 θ 0       (21)
where

M = N -Q + 1 ≥ Q.
In order to satisfy Assumption 3.5, we can choose I = Q and, for every i ∈ {1, . . . , I},

I i = m ∈ {1, . . . , M } (m -i) mod I = 0 . (22)
Hence, we have for all i ∈ {1, . . . , I},

∆ i,1 = . . . = ∆ i,M i = Q-1 q=0 |θ q | 2 . ( 23
)
In this case, g can be decomposed as a sum of Q functions, whose proximity operators can be easily calculated.

2. One-dimensional zero-padded convolutive models.

The following Toeplitz matrix is considered:

   t ⊤ 1 . . . t ⊤ M    =                θ 0 0 . . . .
. . . 0 θ Q-1 . . . θ 1 θ 0                (24)
where M = N ≥ 2Q. In this case, I can be chosen equal to Q and the index sets (I i ) 1≤i≤I are still given by [START_REF] Harmany | This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms -theory and practice[END_REF]. However, the diagonal parameters are not all equal as in the previous example. We have indeed, for every i ∈ {1, . . . , I},

∆ i,1 = i-1 q=0 |θ q | 2 ∆ i,2 = . . . = ∆ i,M i = Q-1 q=0 |θ q | 2 . (25)
3. One-dimensional periodic convolutive models.

In this case, a matrix having a circulant structure [START_REF] Golub | Matrix computations[END_REF] is involved:

   t ⊤ 1 . . . t ⊤ M    =                θ 0 0 . . . 0 θ Q-1 . . . θ 1 θ 1 θ 0 . . . . . . . . . . . . . . . . . . . . . . . . θ Q-1 θ Q-1 . . . . . . 0 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0 . . . 0 θ Q-1 . . . θ 1 θ 0                (26)
where M = N ≥ Q. In order to satisfy Assumption 3.5, we subsequently set

I = min{i ≥ Q | (M -i) mod Q = 0}
and, for every i ∈ {1, . . . , I},

I i =      {i} if i ≤ Q -1 m ∈ {Q, . . . , M } (m -i) mod Q = 0 otherwise. ( 27
)
The diagonal parameters are then given by [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. Another choice which was made in [START_REF] Pustelnik | Hybrid regularization for data restoration in the presence of Poisson noise[END_REF] is to set I = min{i ≥ Q | M mod i = 0} and to proceed as in [START_REF] Harmany | This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms -theory and practice[END_REF] and [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. This solution may be preferred due to its simplicity, when the resulting value of I is small. 4. One-dimensional d-decimated zero-padded convolutive models.

We get the following matrix of R M ×N where N = M d ≥ 2Q.

   t ⊤ 1 . . . t ⊤ M    =                  θ d-
0 θ Q-1 . . . θ 2d-1 . . . θ d-1 . . . θ 0                  . ( 28
)
In order to satisfy Assumption 3.5, we subsequently set I = Q d and the index sets (I i ) 1≤i≤I are still given by [START_REF] Harmany | This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms -theory and practice[END_REF]. We have indeed, for every i ∈ {1, . . . , I},

∆ i,1 = min(id,Q)-1 q=0 |θ q | 2 ∆ i,2 = . . . = ∆ i,M i = Q-1 q=0 |θ q | 2 . (29)
Note that, when d ≥ Q, (t m ) m∈{1,...,M } is an orthogonal family, and thus I = 1.

Remark 3.7 In the previous example (the non-decimated example being a special case when d = 1), the computational complexity of applying each operator T i or T ⊤ i with i ∈ {1, . . . , I} is O(M d) and we have about Q/d proximity operators prox Υ i •T i to compute. Assuming a complexity O(M i ) for computing prox D i Υ i , the overall computational complexity is O(M (2Q + 1)). In turn, if we choose I = M , the complexity of computation of T i or T ⊤ i is O(Q), but we have about M proximity operators prox Υ i •T i to compute. Thus, the overall computational complexity remains of the same order as previously. This means that limiting the number of proximity operators to be computed has no clear advantage in terms of computational complexity, but it allows us to reduce the memory requirement (gain of a factor M d/Q for the storage of the results of the proximity operators).

Discrete forms of total variation and associated proximity operator

Total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] represents a powerful regularity measure in image restoration for recovering piecewise homogeneous areas with sharp edges [START_REF] Rudin | Total variation based image restoration with free local constraints[END_REF][START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelet for image restoration[END_REF][START_REF] Aujol | Structure-texture image decompositionmodeling, algorithms, and parameter selection[END_REF][START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF][START_REF] Bresson | Fast dual minimization of the vectorial total variation norm and applications to color image processing[END_REF]. Different versions of discretized total variation can be found in the literature [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]. Our objective here is to consider discrete versions for which the proximity operators can be easily computed. The main idea will be to split the total variation term in a sum of functions the proximity operators of which have a closed form. The considered form of the total variation of a digital image y = (y

n 1 ,n 2 ) 0≤n 1 <N 1 ,0≤n 2 <N 2 ∈ R N 1 ×N 2 is tv(y) = N 1 -P 1 n 1 =0 N 2 -P 2 n 2 =0 ρ tv (y h ) n 1 ,n 2 , (y v ) n 1 ,n 2 , (30)
where ρ tv ∈ Γ 0 (R 2 ), and y h and y v are two discrete gradients computed in orthogonal directions through FIR filters with impulse responses of size P 1 × P 2 . More precisely, in the above expression, we have

(y h ) n 1 ,n 2 = tr(H ⊤ Y n 1 ,n 2 ) (y v ) n 1 ,n 2 = tr(V ⊤ Y n 1 ,n 2 ) ( 31
)
where H ∈ R P 1 ×P 2 and V ∈ R P 1 ×P 2 are the filter kernel matrices here assumed to have unit Frobenius norm, and for every (n 1 , n 2 ) ∈ {0, . . . , N 1 -

P 1 } × {0, . . . , N 2 -P 2 }, Y n 1 ,n 2 = (y n 1 +p 1 ,n 2 +p 2 ) 0≤p 1 <P 1 ,0≤p 2
<P 2 denotes a block of P 1 × P 2 neighbouring pixels. Since the proximity operator associated with the so-defined total variation does not take a simple expression in general, (30) can be split in "block terms" by following an approach similar to that in Section 3.2:

(∀y ∈ R N 1 ×N 2 ) tv(y) = P 1 -1 p 1 =0 P 2 -1 p 2 =0 tv p 1 ,p 2 (y) (32)
where, for every p 1 ∈ {0, . . . , P 1 -1} and p 2 ∈ {0, . . . , P 2 -1},

tv p 1 ,p 2 (y) = ⌊ N 1 -p 1 P 1 ⌋-1 n 1 =0 ⌊ N 2 -p 2 P 2 ⌋-1 n 2 =0 ρ tv (y h ) p 1 ,p 2 n 1 ,n 2 , (y v ) p 1 ,p 2 n 1 ,n 2 (33)
and the notation (•) p 1 ,p 2 n 1 ,n 2 = (•) P 1 n 1 +p 1 ,P 2 n 2 +p 2 has been used. A closed form expression for the proximity operator of the latter function can be derived as shown below (the proof is provided in Appendix 8). Proposition 3.8 Under the assumption that tr(HV ⊤ ) = 0, for every

y = (y n 1 ,n 2 ) 0≤n 1 <N 1 ,0≤n 2 <N 2 ∈ R N 1 ×N 2
and µ > 0, we have

(∀(p 1 , p 2 ) ∈ {0, . . . , P 1 -1} × {0, . . . , P 2 -1}) prox µtvp 1 ,p 2 y = (π n 1 ,n 2 ) 0≤n 1 <N 1 ,0≤n 2 <N 2 (34)
where, for every (n 1 , n 2 ) ∈ {0, . . . ,

⌊ N 1 -p 1 P 1 ⌋ -1} × {0, . . . , ⌊ N 2 -p 2 P 2 ⌋ -1}, (π P 1 n 1 +p 1 +p ′ 1 ,P 2 n 2 +p 2 +p ′ 2 ) 0≤p ′ 1 <P 1 ,0≤p ′ 2 <P 2 = (β p 1 ,p 2 n 1 ,n 2 -h p 1 ,p 2 n 1 ,n 2 )H + (κ p 1 ,p 2 n 1 ,n 2 -v p 1 ,p 2 n 1 ,n 2 )V + Y p 1 ,p 2 n 1 ,n 2 (35) with h p 1 ,p 2 n 1 ,n 2 = tr(H ⊤ Y p 1 ,p 2 n 1 ,n 2 ), v p 1 ,p 2 n 1 ,n 2 = tr(V ⊤ Y p 1 ,p 2 n 1 ,n 2 ) ( 36
) (β p 1 ,p 2 n 1 ,n 2 , κ p 1 ,p 2 n 1 ,n 2 ) = prox µ ρtv (h p 1 ,p 2 n 1 ,n 2 , v p 1 ,p 2 n 1 ,n 2 ), ( 37
)
and (∀(n 1 , n 2 ) ∈ {0, . . . , N 1 -1} × {0, . . . , N 2 -1}) π n 1 ,n 2 = y n 1 ,n 2 if            n 1 < p 1 or n 2 < p 2 or n 1 ≥ P 1 ⌊ N 1 -p 1 P 1 ⌋ or n 2 ≥ P 2 ⌊ N 2 -p 2 P 2 ⌋. (38)
The result in Proposition 3.8 basically means that, for a given value of (p 1 , p 2 ) ∈ {0, . . . , P 1 -1}× {0, . . . , P 2 -1}, the image is decomposed into non-overlapping blocks

Y p 1 ,p 2 n 1 ,n 2 = (y P 1 n 1 +p 1 +p ′ 1 ,P 2 n 2 +p 2 +p ′ 2 ) 0≤p ′ 1 <P 1 ,0≤p ′
2 <P 2 of P 1 × P 2 pixels. Eq. ( 35) then provides the expression of the proximity operator associated with each one of these blocks, whereas [START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF] deals with boundary effects. Remark 3.9 The above result offers some degrees of freedom in the definition of the discretized total variation for the choices of the function ρ tv and of the gradient filters.

• Two classical choices for the function ρ tv [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] are the following:

1. If ρ tv : (η 1 , η 2 ) → |η 1 |+|η 2 |
then, an anisotropic form is obtained. According to Example 3.2, (37) reduces to

β p 1 ,p 2 n 1 ,n 2 = sign(h p 1 ,p 2 n 1 ,n 2 ) max(|h p 1 ,p 2 n 1 ,n 2 | -µ, 0) κ p 1 ,p 2 n 1 ,n 2 = sign(v p 1 ,p 2 n 1 ,n 2 ) max(|v p 1 ,p 2 n 1 ,n 2 | -µ, 0). ( 39
)
2. If ρ tv : (η 1 , η 2 ) → (η 1 ) 2 + (η 2 ) 2 , then the standard isotropic form is found. The proximity operator involved in ( 37) is given in Example 3.3.

• Some examples of kernel matrices H and V satisfying the assumptions of Proposition 3.8 are as follows:

1. Roberts filters such that H =

-1/ √ 2 0 0 1/ √ 2 and V = 0 -1/ √ 2 1/ √ 2 0
were investigated in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF].

2. Finite difference filters can be used, which are such that

H = V ⊤ =   0 0 0 -1/ √ 2 0 1/ √ 2 0 0 0   .
3. Prewitt filters also satisfy the required assumptions. They are defined by

H = V ⊤ =   -1/ √ 6 0 1/ √ 6 -1/ √ 6 0 1/ √ 6 -1/ √ 6 0 1/ √ 6   .

Sobel filters such that

H = V ⊤ =   -1/ √ 12 0 1/ √ 12 -2/ √ 12 0 2/ √ 12 -1/ √ 12 0 1/ √ 12
  are possible choices too.

Proposed algorithm

In the class of convex optimization methods, an algorithm recently proposed in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] appears wellsuited to solve the class of the minimization problems formulated as in Problem (5). However, when synthesis frame representations are considered (Problem (6)) and when the function number S is large, the frame analysis and synthesis operators have to be applied several times in the algorithm which induces a long computation time. In this section, we briefly recall the Parallel ProXimal Algorithm and its convergence properties. Then, we propose an improved version of PPXA to efficiently solve Problem (6).

Parallel ProXimal Algorithm (PPXA)

An equivalent formulation of the convex optimization problem ( 5) is:

Find x ∈ Argmin x 1 ∈R K ,...,x J ∈R K x=x 1 =...=x J J j=1 f j (x j ). ( 40
)
This formulation was used in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] to derive Algorithm 1.

Algorithm 1 General form of PPXA Set γ ∈ ]0, +∞[. For every j ∈ {1, . . . , J}, set (ω j ) 1≤j≤J ∈]0, 1] J such that J j=1 ω j = 1. Set (u j,0 ) 1≤j≤J ∈ (R K ) J and x 0 = J j=1 ω j u j,0 . For ℓ = 0, 1, . . .

             For j = 1, . . . , J ⌊ p j,ℓ = prox γf j /ω j u j,ℓ + a j,ℓ p ℓ = J j=1 ω j p j,ℓ Set λ ℓ ∈ ]0, 2[ For j = 1, . . . , J ⌊ u j,ℓ+1 = u j,ℓ + λ ℓ (2 p ℓ -x ℓ -p j,ℓ ) x ℓ+1 = x ℓ + λ ℓ (p ℓ -x ℓ )
PPXA involves real constants γ and (ω j ) 1≤j≤J , and, at each iteration ℓ ∈ N, a relaxation parameter λ ℓ . It also includes possible error terms (a j,ℓ ) 1≤j≤J in the computation of the proximity operators, which shows the numerical stability of the algorithm. The sequence (x ℓ ) ℓ≥1 generated by Algorithm 1 can be shown to converge to a solution to Problem [START_REF] Hestenes | Multiplier and gradient methods[END_REF] (or equivalently to Problem (5)) under the following assumption [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF].

Assumption 4.1 1. lim x →+∞ f 1 (x) + . . . + f J (x) = +∞. 2. ∩ J j=1 rint dom f j = ∅. 2 3. (∀j ∈ {1, . . . , J}) ℓ∈N λ ℓ a j,ℓ < +∞. 4. ℓ∈N λ ℓ (2 -λ ℓ ) = +∞. Remark 4.2
The fact that the algorithm involves several parameters should not be viewed as a weakness since the convergence is guaranteed for any choice of these parameters under the previous assumption. These parameters bring out flexibility in PPXA in the sense that an appropriate choice of them (typical values will be indicated in Section 5.2) may be beneficial to the convergence speed.

Consider now Problem [START_REF] Combettes | The Convex Feasibility Problem in Image Recovery[END_REF] where a tight frame is employed (F ⊤ F = ν Id). By setting (∀j ∈ {1, . . . , S}) f j = g j • F ⊤ and by invoking Proposition 3.4 with L = F ⊤ and D = ν Id, the iterations of Algorithm 1 become as described in Algorithm 2.

Algorithm 2 PPXA iterations for Problem [START_REF] Combettes | The Convex Feasibility Problem in Image Recovery[END_REF] For ℓ = 0, 1, . . .

                   For j = 1, . . . , S ⌊ p j,ℓ = u j,ℓ + 1 ν F prox νγg j /ω j (F ⊤ u j,ℓ ) -F ⊤ u j,ℓ + a j,ℓ For j = S + 1, . . . , J ⌊ p j,ℓ = prox γf j /ω j u j,ℓ + a j,ℓ p ℓ = J j=1 ω j p j,ℓ Set λ ℓ ∈ ]0, 2[ For j = 1, . . . , J ⌊ u j,ℓ+1 = u j,ℓ + λ ℓ (2 p ℓ -x ℓ -p j,ℓ ) x ℓ+1 = x ℓ + λ ℓ (p ℓ -x ℓ )
However, the first loop can be costly in terms of computational complexity because it requires to apply S times the operators F and F ⊤ at each iteration. We will now see how it is possible to speed up these iterations.

Accelerated version of PPXA

In Algorithm 3, we propose an adaptation of PPXA in order to reduce its computational load by limiting the number of times the operators F and F ⊤ are applied. Details concerning the derivation of this algorithm can be found in Appendix 9.

Let us make the following assumption:

Assumption 4.3 1. lim x →+∞ g 1 (F ⊤ x) + . . . + g S (F ⊤ x) + f S+1 (x) + . . . + f J (x) = +∞. 2. ∩ S j=1 rint dom (g j • F ⊤ ) ∩ J j=S+1 rint dom f j = ∅.
3. (∀j ∈ {1, . . . , S}) ℓ∈N λ ℓ a j,ℓ < +∞ and (∀j ∈ {S + 1, . . . , J}) ℓ∈N λ ℓ a j,ℓ < +∞.

4.

ℓ∈N λ ℓ (2 -λ ℓ ) = +∞.

Then, Algorithm 3 converges to a solution to Problem [START_REF] Combettes | The Convex Feasibility Problem in Image Recovery[END_REF]. In addition, this algorithm requires only 3 applications of F or F ⊤ at each iteration. Hence, a gain w.r.t. Algorithm 2 is obtained as soon as S ≥ 2. This fact will be illustrated by our simulation results in Section 5.2.1.

Algorithm 3 Accelerated PPXA Let γ ∈ ]0, +∞[. For every j ∈ {1, . . . , J}, set (ω j ) 1≤j≤J ∈]0, 1] J such that J j=1 ω j = 1. Set (u j,0 ) 1≤j≤J ∈ (R K ) J and x 0 = J j=1 ω j u j,0 . For every j ∈ {1, . . . , S}, set v j,0 = F ⊤ u j,0 and u ⊥ j,0 = u j,0 -1 ν F v j,0 . For ℓ = 0, 1, . . .

                            For j = 1, . . . , S ⌊ q j,ℓ = 1 ν prox νγg j /ω j v j,ℓ + a j,ℓ For j = S + 1, . . . , J ⌊ p j,ℓ = prox γf j /ω j u j,ℓ + a j,ℓ p ℓ = S j=1 ω j u ⊥ j,ℓ +F S j=1 ω j q j,ℓ + J j=S+1 ω j p j,ℓ r ℓ = 2 p ℓ -x ℓ ; r ℓ = F ⊤ r ℓ ; r ⊥ ℓ = r ℓ -1 ν F r ℓ Set λ ℓ ∈ ]0, 2[ For j = 1, . . . , S u ⊥ j,ℓ+1 = u ⊥ j,ℓ + λ ℓ (r ⊥ ℓ -u ⊥ j,ℓ ) v j,ℓ+1 = v j,ℓ + λ ℓ ( r ℓ -νq j,ℓ ) For j = S + 1, . . . , J ⌊ u j,ℓ+1 = u j,ℓ + λ ℓ (r ℓ -p j,ℓ ) x ℓ+1 = x ℓ + λ ℓ (p ℓ -x ℓ )
5 Application to restoration

Hybrid regularization

In restoration problems, one of the terms in the criterion to be minimized usually is a fidelity term measuring some distance between the image degraded by the operator T and the observed data z. We will assume that this function takes the form g = Ψ • T where Ψ ∈ Γ 0 (R M ). In the case of data corrupted by a additive zero-mean white Gaussian noise with variance α, a standard choice for Ψ is a quadratic function such that Ψ = 1 2α • -z 2 . Then, the associated proximity operator of g can be computed explicitly (see [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]). In the case of data contaminated by an independent Poisson noise with scaling parameter α, a standard choice is Ψ = D KL (z, α •) where D KL is the generalized Kullback-Leibler divergence [START_REF] Titterington | On the iterative image space reconstruction algorithm for ECT[END_REF][START_REF] Fessler | Hybrid Poisson/polynomial objective functions for tomographic image reconstruction from transmission scans[END_REF][START_REF] Zheng | Parallelizable bayesian tomography algorithms with rapid, guaranteed convergence[END_REF][START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF][START_REF] Figueiredo | Restoration of Poissonian images using alternating direction optimization[END_REF][START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF][START_REF] Chouzenoux | Majorize-minimize linesearch for inversion methods involving barrier function optimization[END_REF] such that,

∀u = (u (m) ) 1≤m≤M ∈ R M , Ψ(u) = M m=1 ψ m (u (m) ) (41)
and

ψ m (u (m) ) =            αu (m) -z (m) + z (m) ln z (m)
αu (m) if z (m) > 0 and u (m) > 0, αu (m) if z (m) = 0 and u (m) ≥ 0, +∞ otherwise.

(

) 42
The proximity operator of Ψ can then be derived from Example 3.1.

Concerning regularization functions, a standard choice of penalty function in the wavelet domain is:

(∀x = (x (k) ) 1≤k≤K ∈ R K ), Φ(x) = K k=1 φ k (x (k) )
where, for every k ∈ {1, . . . , K}, φ k is a finite function of Γ 0 (R) such that lim |x (k) |→+∞ φ k (x (k) ) = +∞. Power functions as in Example 3.2 are often chosen for (φ k ) 1≤k≤K (see e.g. [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]). The main problem with wavelet regularization is the occurence of some visual artefacts (e.g. ringing artefacts), some of which can be reduced by increasing the redundancy of the representation. Another popular type of regularization that can be envisaged consists of employing a total variation measure [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. Its major drawback is the generation of staircase-like effects in the recovered images. To combine the advantages of both regularizations, we propose to:

Find x ∈ Argmin x∈R K Ψ(T F ⊤ x) + µ tv(F ⊤ x) + ι C (F ⊤ x) + ϑ Φ(x). ( 43
)
As already mentioned, Φ corresponds to the regularization term operating in the wavelet domain. tv represents a discrete total variation term as defined by [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]. Finally, ι C is the indicator function of a nonempty closed convex set C of R N (for example, related to support or value range contraints). This kind of objective function was also recently investigated in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] but the approach was restricted to the use of a quadratic data fidelity term and of a specific form of the total variation term.

The non-negative real parameters ϑ and µ control the degree of smoothness in the wavelet and in the space domains, respectively.

The main difficulty in applying Algorithm 1 to our restoration problem is that it requires to compute the proximity operators associated with each of the four terms in [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF]. In general, closed forms of the proximity operators are known only for the indicator function ι C and for Φ [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]. However, as explained in Section 3.2, provided that the function Ψ is separable, the data fidelity term can be decomposed as a sum of I functions (Υ i • T i ) 1≤i≤I for which the proximity operators can be calculated according to [START_REF] Figueiredo | An EM algorithm for wavelet-based image restoration[END_REF]. Similarly, by using the results in Section 3.3, the tv function can be split in P 1 P 2 functions (tv p 1 ,p 2 ) 0≤p 1 <P 1 ,0≤p 2 <P 2 , the proximity operators of which are given by Proposition 3.8. Algorithm 3 can then be applied with S = I + P 1 P 2 + 1 and J = I +P 1 P 2 +2. In the present case, it can be noticed that if ϑ > 0, Assumption 4.3 1) is satisfied.

In addition, Assumption 4.3 2) is fulfilled ∩ I i=1 y ∈ R N T i y ∈ rint dom Υ i rint C = ∅ (since dom Φ = R K and (∀(p 1 , p 2 ) ∈ {0, . . . , P 1 -1} × {0, . . . , P 2 -1}) dom tv p 1 ,p 2 = R N ). This condition is verified if ]0, +∞[ M ⊂ dom Ψ and C = [0, 255]
N since for every i ∈ {1, . . . , I}, T i has been assumed non-negative real valued in Section 2.1, and with non-zero lines (see Assumption 3.5).

Experimental results for convolved data in the presence of Poisson noise

In our simulations, we will be first interested in studying the performance in terms of convergence rate of the accelerated version of PPXA. Algorithms 2 and 3 are implemented by setting γ = 50, λ ℓ ≡ 1.6 and, for every j ∈ {1, . . . , J},

ω j =      1 4I if 1 ≤ j ≤ I 1 4P 1 P 2 if I + 1 ≤ j ≤ I + P 1 P 2 1 4
otherwise [START_REF] Figueiredo | Restoration of Poissonian images using alternating direction optimization[END_REF] if (g j ) 1≤j≤I are the functions corresponding to the decomposition of the data fidelity term and (g j ) I+1≤j≤I+P 1 P 2 correspond to the decomposition of tv. The weights are thus chosen to provide equal contributions to the four functions in Criterion [START_REF] Goldstein | The split Bregman method for ℓ 1 regularized problems[END_REF]. For the first and second functions which are splitted, the corresponding 1/4 weight is further subdivided in a uniform manner. Note however that the behaviour of the algorithm did not appear to be very sensitive to an accurate choice of these parameters. A comparison between the different total variation regularization terms defined in Section 4.2 will also be made. Another discussion will be held concerning the boundary effects. Two cases will be considered: the use a periodic convolution and then, of a convolution with zero-padding. Results for a decimated convolution will also be presented. Finally, the interest in combining total variation and wavelet regularization terms will be shown with respect to classical regularizations. A tight frame version of the dual-tree transform (DTT) proposed in [START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF] (ν = 2) using Symlets of length 6 over 3 resolution levels is employed. We choose potential functions of the form: for every k ∈ {1, . . . , K}, φ k = χ k | • | p k where χ k > 0 and p k ∈ {1, 4/3, 3/2, 2}, the proximity operators of which are given by Example 3.2.

Convergence rate comparison between PPXA and its accelerated version

Table 1 gives the iteration numbers and the CPU times for the original PPXA algorithm and the proposed accelerated one in order to reach convergence when considering different image sizes ("Sebal": N = 128 × 128, "Peppers": N = 256 × 256 and "Marseille": N = 512 × 512) and various kernel blur sizes. The stopping criterion is based on the relative error between the objective function computed at the current iteration and at the previous one. 3 The stopping tolerance has been set to 10 -3 . These results have been obtained with an Intel Core2 6700, 2.66 GHz. The last line of Table 1 illustrates the gain in CPU-time when using Algorithm 3. Moreover, in Figure 1, the mean square error on the image iterates F ⊤ (x nx) It can be noticed that the larger the kernel blur size is, the higher the gain is. This is due to the fact that the number of proximity operators to compute increases with the kernel size.

Comparison between different forms of total variation

In Section 3.3, we have introduced the proximity operator associated with discretized total variation functions and the possibility of choosing various filters has been mentioned. In Figure 2, tests have been carried out on "Peppers" degraded by a 3 × 3 uniform blur and corrupted by Poisson noise with scaling parameter α = 0.1. We compare the restored images for different kinds of total variation, in terms of Signal to Noise Ratio -SNR and structural similarity measure -SSIM [START_REF] Wang | Mean squared error: love it or leave it?[END_REF].

The SSIM takes a value between -1 to 1, the maximum value being obtained for two identical images. Each curve represents the resulting SNR and SSIM versus µ (the regularization parameter related to the total variation), for a given form of tv (i.e. a given filter associated with either an isotropic or anisotropic function ρ tv ). A small wavelet regularization parameter ϑ = 10 -3 has been chosen in order to better illustrate the influence of the different tv forms on restoration quality. It can be concluded from Figure 2 that the choice of the gradient filters and of the form (isotropic/anisotropic) of ρ tv has a significant influence on the restoration quality when the wavelet regularization is small. However, we also noticed in our numerical experiments that when the wavelet regularization parameter ϑ becomes larger, the choice of the tv form has a lower influence on the restoration quality provided that the regularization parameters are appropriately chosen.

Boundary effects on restored images

This section illustrates the influence of boundary effect processing. More precisely, we degraded an extended version of "Boat" image by a 7 × 7 uniform blur, and the resulting blurred image was cropped to create an image of size 256×256. As a consequence, the boundary values are functions of pixel locations which are no longer present in the blurred image. The scaling parameter associated with Poisson noise is α = 0.5. The objective is then to restore the image (which was centered) by using one of the convolution models discussed in Section 3.2, namely either a periodic convolution or a convolution with zero-padding. Visual and quantitative results are given in Figure 3. As it can be noticed from this figure, the periodic convolution model introduces significant boundary artefacts unlike the convolution with zero-padding. The results obtained when considering "Peppers" led to the same conclusion. For "Sebal", zero-padding or periodic models provided similar results.

Decimated convolution

We now present experimental results for a 256×256 SPOT image degraded by a uniform decimated blur with a kernel size Q = 3 × 3 and a decimated factor d = 2. The scaling parameter of the Poisson noise is equal to α = 1. Due to the structure of the degradation operator, the data fidelity is splitted in a sum of I = 4 functions. The results are presented in Figure 4 where the good behaviour of the model can be observed.

Influence of each regularization term

We now present numerical and visual results for the different kinds of regularization when a generalized Kullback-Leibler divergence is used as a data fidelity term. This experiment allows us to compare the hybrid regularization with existing approaches based on a wavelet-frame [START_REF] Figueiredo | Restoration of Poissonian images using alternating direction optimization[END_REF] regularization or a total variation regularization [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF][START_REF] Figueiredo | Restoration of Poissonian images using alternating direction optimization[END_REF]. The latter regularized solutions can be computed either by using augmented Lagrangian techniques [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF][START_REF] Figueiredo | Restoration of Poissonian images using alternating direction optimization[END_REF] or with our splitting approach (by setting ϑ = 0 or µ = 0). In our experiments, the computation time of the two approaches was observed to be similar. Note that comparisons performed in [START_REF] Figueiredo | Restoration of Poissonian images using alternating direction optimization[END_REF] led to the conclusion that the wavelet-frame regularization is quite competitive with respect to other existing restoration methods [START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF].

In the images displayed in Figures 5,6, and 7, one can observe the artefacts related to the wavelet regularization, the staircase effects which are typical of the total variation penalization, some checkerboard patterns resulting of the chosen gradient discretization, and also the benefits which can be drawn from the use of a hybrid regularization.

Similarly to [START_REF] Figueiredo | Restoration of Poissonian images using alternating direction optimization[END_REF], the values of µ and ϑ have been adjusted so as to maximize the SNR. Optimizing the hyperparameters manually as we did is a common practice in imaging applications, especially when a data set of test images having similar characteristics as the one to be restored (medical images, satellite images,...) is available. Automatic methods for the optimization of the hyperparameters can also be found in the literature such as cross-validation [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF], Stochastic EM [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF], MCMC [START_REF] Robert | Monte Carlo statistical methods[END_REF][START_REF] Chaâri | A hierarchical bayesian model for frame representation[END_REF] or Stein-based methods [START_REF] Ramani | Monte-Carlo SURE: A black-box optimization of regularization parameters for general denoising algorithms[END_REF]. These automatic procedures often are relatively intensive. They will not be addressed in this paper due to the lack of space.

Conclusion

A new convex regularization approach to restore data degraded by a (possibly decimated) convolution operator and a non necessarily additive noise has been proposed. The main advantages of the method are (i) to deal directly with the "true" noise likelihood (i.e. the Kullback-Leibler divergence in the case of Poisson noise) without requiring any approximation of it; (ii) to permit the use of sophisticated regularization functions, e.g. one promoting sparsity in a wavelet frame domain and a total variation penalization. In addition, the proposed algorithm has a parallel structure which makes it easily implementable on multicore architectures. Numerical and visual results demonstrate the effectiveness of the proposed approach. One can note that, even if the paper is devoted to the case of convolutive operators, this approach could be generalized to more general linear operators.

Note that the primal-dual approaches [START_REF] Zhu | An efficient primal-dual hybrid gradient algorithm for total variation image restoration[END_REF][START_REF] Zhang | A unified primal-dual algorithm framework based on Bregman iteration[END_REF][START_REF] Esser | A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science[END_REF][START_REF] Setzer | Infimal convolution regularizations with discrete ℓ 1type functionals[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] can offer alternative solutions to the ones developed in this paper. However, one of the advantages of PPXA is that it easily leads to efficient parallel implementations 7 Proof of Proposition 3.4

Since D = LL ⊤ is the matrix associated with a bijective operator, L is associated with a surjective one and dom Υ = ∅ ⇒ dom (Υ • L) = ∅. This allows us to conclude that h = Υ • L is a function of Γ 0 (R X ). To calculate the proximity operator of h, we now come back to the definition of this operator. We have thus, for every w ∈ R X , prox h w = arg min

v∈R X 1 2 v -w 2 + Υ(Lv). ( 45
)
We can write any vector v ∈ R X as a sum of an element L ⊤ t ∈ ran L ⊤ and v ⊥ ∈ (ran L ⊤ ) ⊥ = ker L.

We have then Lv = LL ⊤ t = Dt. Similarly, we can write w = L ⊤ u + w ⊥ where u ∈ R Y and w ⊥ ∈ ker L. So, prox h w can be determined by finding min

(t,v ⊥ )∈R Y ×R X 1 2 L ⊤ t + v ⊥ -L ⊤ u -w ⊥ 2 + Υ(Dt) = min (t,v ⊥ )∈R Y ×R X 1 2 L ⊤ (t -u) 2 + 1 2 v ⊥ -w ⊥ 2 + Υ(Dt). ( 46
)
This yields v ⊥ = w ⊥ = w -L ⊤ u [START_REF] Wen | Iterative algorithms based on decoupling of deblurring and denoising for image restoration[END_REF] and it remains to find min t∈R Y

It can be deduced from [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]Lemma 2.6] that, for every m ∈ {1, . . . , Y },

o m ⊤ t = prox 1 ∆m ψm(∆m•) (o m ⊤ u) = 1 ∆ m prox ∆mψm (∆ m o m ⊤ u), (50)
which, according to [26, Proposition 2.10], leads to

t = D -1 Y m=1 prox ∆mψm (∆ m o m ⊤ u) o m = D -1 prox DΥ (Du). (51)
Altogether, ( 47) and ( 51) yield v = w + L ⊤ D -1 prox DΥ (Du) -u .

In addition, since L ⊤ u is the projection of w onto ran L ⊤ , u = (LL ⊤ ) -1 Lw = D -1 Lw and ( 17) follows.

8 Proof of Proposition 3.8

By using the proximity operator definition [START_REF] Trussell | The feasible solution in signal restoration[END_REF],

π = prox
) 68
This allows us to transpose the convergence results concerning Algorithm 2 to Algorithm 3.

Figure 1 :

 1 Figure 1: Convergence profiles of the Algorithm 2 (dotted line) and Algorithm 3 (solid line) versus computation time in seconds for a 3 × 3 uniform blur (left) and a 7 × 7 uniform blur (right) and a 128 × 128-image.

Figure 2 :

 2 Figure 2: SNR (left) and SSIM (right) for different total variation terms with respect to µ. Roberts: thin-black line, Finite difference: thick-black, Prewitt: thick-gray, ρ tv = | • | 2 + | • | 2 : solid line and ρ tv = | • | + | • |: dashed line

SNR = 16

 16 .9 dB -SSIM = 0.62 SNR = 17.7 dB -SSIM = 0.64

Figure 3 :

 3 Figure 3: Periodic (left) and zero-padded (right) restoration.

Figure 4 :

 4 Figure 4: Restoration results for "Spot" image.

1 2 L

 12 ⊤ (t -u) 2 + Υ(Dt) = min

  u) ⊤ D(t -u) + Υ(Dt).[START_REF] Daubechies | Ten lectures on wavelets[END_REF] By using the separability of Υ, this is equivalent to findingmin o m ⊤ t -o m ⊤ u) 2 + ψ m (∆ m o m ⊤ t).

Table 1 :

 1 2 is plotted as a function of computation time, where (x n ) n>0 denotes the sequence generated by Algorithm 2 or Algorithm 3. Comparisons between PPXA and its accelerated version.

	Image size (uniform) blur size Iteration numbers	128 × 128 3 × 3 7 × 7 3 × 3 7 × 7 3 × 3 7 × 7 256 × 256 512 × 512 30 50 41 50 50 50
	CPU time (in second)	117.2 633.0 411.7 1298	1458 4514
	CPU time -accelerated version (in second) 13.53 29.82 60.59 89.48 263.6 405.0
	Gain	8.67	21.2	6.79	14.5	5.53	11.1

  µtvp 1 ,p 2 (y)(n 1 ,n 2 )∈B (π n 1 ,n 2 -y n 1 ,n 2 ) 2 ,p 2 n 1 ,n 2 -Y p 1 ,p 2 n 1 ,n 2 2 F + µ ρ tv tr(H ⊤ Π p 1 ,p 2 n 1 ,n 2 ), tr(V ⊤ Π p 1 ,p 2 n 1 ,n 2 ) (52)where • F denotes the Frobenius norm, for every (n 1 , n 2 ) ∈ {0, . . . ,⌊ N 1 -p 1 P 1 ⌋-1}×{0, . . . , ⌊ N 2 -p 2 P 2 ⌋-1}, Π p 1 ,p 2 n 1 ,n 2 = (π P 1 n 1 +p 1 +p ′ 1 ,P 2 n 2 +p 2 +p ′ 2 ) 0≤p ′ 1 <P 1 ,0≤p ′ 2 <P 2(53)which leads to Algorithm 3. We finally note that Assumption 4.3 3) implies that Assumption 4.1 4) is satisfied since, for every j ∈ {1, . . . , S},

	minimizes
	1 2	π-y 2 + µtv p 1 ,p 2 (π)
		=	1 2
			+	⌊	N 1 -p 1 P 1 n 1 =0 ⌋-1	⌊	N 2 -p 2 P 2 n 2 =0 ⌋-1	1 2	Π p 1

ℓ∈N λ ℓ a j,ℓ = ℓ∈N λ ℓ F a j,ℓ ≤ F ℓ∈N λ ℓ a j,ℓ < +∞.

(

In finite dimension, the upper bound condition is always satisfied.

The relative interior of a set S of R X is designated by rint S and the domain of a function f :R X →] -∞, +∞] is dom f = {x ∈ R X |f (x) < +∞}.

The relative error was evaluated based on Criterion[START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF] where the indicator function was discarded.

* Part of this work appeared in the conference proceedings of EUSIPCO 2009 [1]. This work was supported by the Agence Nationale de la Recherche under grant ANR-09-EMER-004-03.

Total variation regularization

and

It is then clear that [START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF] holds since the variables

In addition, since it has been assumed that tr(H ⊤ V ) = 0 and

,n 2 can be decomposed in an orthogonal manner as follows:

where

and (h p 1 ,p 2 n 1 ,n 2 , v p 1 ,p 2 n 1 ,n 2 ) is given by [START_REF] Fornasier | Domain decomposition methods for linear inverse problems with sparsity constraints[END_REF]. After some simplications, we have thus to minimize, for every

This shows that ( 37) is satisfied and that (Π

Eq. ( 35) straightforwardly follows.

Derivation of Algorithm 3

Let Π F denote the projector on ran F . For every u ∈ R K , we have

where u ⊥ ∈ (ran F ) ⊥ is the projection error and there exists q ∈ R N such that

By combining this with the fact that F ⊤ u ⊥ = 0, we obtain the relation,

which allows us to deduce from (59) that

Now, consider the first step of Algorithm 2: (∀j ∈ {1, . . . , S})

where a j,ℓ is assumed to belong to ran F , i.e. a j,ℓ = F a j,ℓ with a j,ℓ ∈ R N . Defining q j,ℓ ∈ R N similarly to (60) yields Π F p j,ℓ = F q j,ℓ . According to (61), q j,ℓ is such that

By combining ( 62) and ( 63), q j,ℓ = 1 ν prox νγg j /ω j (v j,ℓ ) + a j,ℓ where v j,ℓ = F ⊤ u j,ℓ .

Moreover, since p j,ℓ = F q j,ℓ + p ⊥ j,ℓ , the computation of the variable p ℓ = J j=1 ω j p j,ℓ in Algorithm 2 can be rewritten as

where, according to (59), ( 60), ( 62) and ( 64),

In the new formulation, the last steps of the algorithm consist of updating u ⊥ j,ℓ and v j,ℓ , for all j ∈ {1, . . . , S}. We propose to define r ℓ = 2p ℓ -x ℓ , r ℓ = F ⊤ r ℓ and r ⊥ ℓ = r ℓ -1 ν F rℓ , which yields v j,ℓ+1 = v j,ℓ + λ ℓ r ℓ -F ⊤ p j,ℓ and u ⊥ j,ℓ+1 = u ⊥ j,ℓ + λ ℓ r ⊥ ℓ -p ⊥ j,ℓ . By using ( 63) and (66), these relations can be simplified as      v j,ℓ+1 = v j,ℓ + λ ℓ r ℓ -νq j,ℓ and u ⊥ j,ℓ+1 = u ⊥ j,ℓ + λ ℓ r ⊥ ℓ -u ⊥ j,ℓ ,