Caroline Chaux

Laure Blanc-Féraud

Josiane Zerubia

Wavelet-based restoration methods: application to 3D confocal microscopy images

Keywords: restoration, deconvolution, 3D images, confocal microscopy, Poisson noise, wavelets

We propose in this paper an iterative algorithm for 3D confocal microcopy image restoration. The image quality is limited by the diffraction-limited nature of the optical system which causes blur and the reduced amount of light detected by the photomultiplier leading to a noise having Poisson statistics. Wavelets have proved to be very effective in image processing and have gained much popularity. Indeed, they allow to denoise efficiently images by applying a thresholding on coefficients. Moreover, they are used in algorithms as a regularization term and seem to be well adapted to preserve textures and small objects. In this work, we propose a 3D iterative wavelet-based algorithm and make some comparisons with state-of-the-art methods for restoration.

INTRODUCTION

Confocal microscopy 1 is a powerful tool allowing to get an accurate view of cellular physiology, thus explaining its growing use in biological specimen imaging. Biological image quality is limited by two factors: a blur due to microscope optics (modeled by the Point Spread Function (PSF) denoted here by h), even if it is limited in confocal microscopy with respect to widefield microscopy, and the detection which is done with low photon flow leading to Poisson statistics. The Richardson-Lucy algorithm [START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF][START_REF] Lucy | An iterative technique for rectification of observed distributions[END_REF] is often used to restore images in case of Poisson noise because of its well fitting. A regularization term has also been added in order to stabilize the inversion. More recently, some authors proposed to use the wavelet transform [START_REF] Colicchio | Wavelet filtering applied to 3D wide field fluorescence microscopy deconvolution[END_REF][START_REF] Boutet De Monvel | Image Restoration for Confocal Microscopy: Improving the Limits of Deconvolution, with Application to the Visualization of the Mammalian Hearing Organ[END_REF] for biological image deconvolution. In the literature of Poisson noise, many people use a pre-processing like the Anscombe [START_REF] Anscombe | The transformation, of Poisson, binomial and negative-binomial data[END_REF] transform or the Fisz transform 7 in order to stabilize the noise variance and to apply more usual methods [START_REF] Besbeas | A comparative simulation study of wavelet shrinkage estimators for poisson counts[END_REF] (i.e. techniques developed for Gaussian noise). In the biological context, different approaches have been introduced: deconvolution and denoising steps are often splitted and then applied consecutively. The approach considered here combines, in a unique algorithm, denoising and deconvolution. In this paper we propose a deconvolution method based on a 3D wavelet decomposition and an iterative deconvolution algorithm operating on the wavelet coefficients. We then make a comparative study of the different approaches considered in the biological microscopy context. The organization of the paper is as follows: in Section 2, we present the relevant background in biological microscopy image restoration. Then we introduce wavelets in Section 3 and the proposed algorithm in Section 4. Finally, we will give comparative results in Section 5 both on simulated and real data.

BACKGROUND

Problem context

Image formation model The considered image formation model is the following:

i = P(o * h) ( 1
)
where i is the observed image, o the original one, P(.) stands for Poisson noise and * denotes the 3D convolution. We work with 3D objects, both for the images and the PSF (Point Spread Function) which is constructed following different models (physical PSF, Gaussian PSF...). Image restoration (including both denoising and deconvolution) consists in estimating the original image, thus obtaining ô, only knowing the observation i (we assume here that the PSF is known).

PSF model

The PSF modelizes the microscope optics. We consider here an analytical form of the PSF described by: h

(x, y, z) = |A R (x, y) * Pλem (x, y, z)| 2 .| Pλex (x, y, z)| 2 (2)
The excitation (resp. emission) wavelength is λ ex (resp. λ em ). Pλ is the 2D Fourier transform on x, y of the pupil function P λ , given by [START_REF] Born | Principles of Optics[END_REF][START_REF] Fitzgerrel | Defocus transfer function for circularly symmetric pupils[END_REF][START_REF] Goodman | Introduction to Fourier Optics[END_REF] :

P λ (u, v, z) = Π ρ u 2 + v 2 .e 2iπ λ .W (u,v,z) (3)
Here the complex term W (u, v, z) = 1 2 z (1cos 2α) is the aberration phase derived from [START_REF] Goodman | Introduction to Fourier Optics[END_REF] . ρ = N A 2λ is the lateral cut-off frequency, and N A = n o sin α the numerical aperture which is related to the amount of light entering the microscope and the immersion medium refractive index n o . The phase W depends only on the defect of focus. This theoretical model of the confocal PSF does not take into account geometrical (e.g. spherical) lens aberrations and refractive index induced diffraction.

Noise nature: Poisson statistics At the same time, we are confronted with Poisson noise. A way to face this problem and use more conventional processing (stabilizing the variance) consists in either:

• applying the Anscombe transform 6 . This transformation is defined as follows:

A(i(x, y, z)) = 2 i(x, y, z) + 3 8 . (4)
• applying the Fisz transform 7 . This transformation is defined as follows:

1. Process a non normalized Haar transform.

2. Transform the detail coefficients d j,m (k) of level j, subband m and spatial position k as

d j,m (k) = d j,m (k)/ √ a j,m (k) if a j,m (k) = 0 0 if a j,m (k) = 0
where a j,m denotes the approximation coefficients.

3. Reconstruct using the inverse non normalized Haar transform.

After such transformations, the noise can be asymptotically considered as N (0, 1). Some authors 12 make the assumptions of a mixture of Gaussian noise (with mean µ and variance σ 2 ) and Poisson noise P(.). They consider the following image formation model:

i = α P(o * h) + N (µ, σ 2 ) ( 5
)
where α > 0 and develop adapted methods for image restoration problem.

In our case, we consider the noise having Poisson statistics.

Richardson-Lucy based algorithms

Many approaches are based on the Richardson-Lucy (RL) algorithm 2, 3 which is well adapted to Poisson noise.

One iteration of this algorithm is given by:

o n+1 (x) = i(x) (o n * h) (x) * h(-x) o n (x) . (6)
Several regularizations, allowing to stabilize the inversion, have been proposed such as:

Tikhonov-Miller regularization 13 o n+1 (x) = i(x) (o n * h) (x) * h(-x) o n (x) 1 + 2λ T M △o n (x) . ( 7
)
Total variation regularization [START_REF] Dey | A deconvolution method for confocal microscopy with total variation regularization[END_REF][START_REF] Dey | Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution[END_REF] The algorithm (RL-TV) is thus yielding to:

o n+1 (x) = i(x) (o n * h) (x) * h(-x) o n (x) 1 -λ div ∇on(x) |∆on(x)| . (8)

Separation of the denoising and deconvolution steps

In the microscopy field, researchers often apply deconvolution and denoising steps consecutively, thus introducing some difficulties. We discuss in this part the advantages and drawbacks of the different approaches. The denoising step often corresponds to a wavelet coefficient thresholding using wavelets [START_REF] Colicchio | Wavelet filtering applied to 3D wide field fluorescence microscopy deconvolution[END_REF][START_REF] Boutet De Monvel | Image Restoration for Confocal Microscopy: Improving the Limits of Deconvolution, with Application to the Visualization of the Mammalian Hearing Organ[END_REF] or steerable pyramids [START_REF] Rooms | Simultaneous degradation estimation and restoration of confocal images and performances evaluation by colocalization analysis[END_REF] . Concerning the deconvolution step, RL algorithm or a MAP (Maximum A Posteriori) method is often used. When Poisson noise is considered, pre-processing stage to stabilize the noise variance is applied previously on the data. Several drawbacks can be pointed out : we wonder whether the PSF is modified during the denoising step (in [START_REF] Rooms | Simultaneous degradation estimation and restoration of confocal images and performances evaluation by colocalization analysis[END_REF] they estimate it after denoising); but also, we are not sure of the noise nature after the denoising step. Wavelets have shown to be really efficient in denoising problems, which leads naturally to method 1. The deconvolution step often corresponds to an inversion of the PSF 17 taking into account it may have zeros. A simple way to achieve it is

Deconvolution then denoising

F(ô) = F (o) F (P SF ) if |F(P SF )| > ǫ 0 otherwise (9)
where F(.) denotes the Fourier transform. ǫ can be taken equal to 0.001.

Concerning the denoising step, it often corresponds to a coefficient thresholding in a transformed domain. The noise level has to be estimated. When Poisson noise is considered, pre-processing stage to stabilize the noise variance is previously applied on the data. The approach creates an amplification of the noise during the deconvolution step and after the inversion of the PSF, the noise nature is not known. The second step (denoising) must take into account this particularity. This method is easy to implement and comes naturally to mind. Nevertheless, without an adaptation and a careful implementation, [START_REF] Kalifa | Deconvolution by thresholding in mirror wavelets bases[END_REF] the noise is too strongly amplified and this method gives rise to bad results.

In order to tackle the different problems above mentioned, we propose to combine the denoising and deconvolution steps using a single algorithm operating in the 3D transformed domain, as we will see later.

3D WAVELET ANALYSIS

This section aims at giving a brief recall of the 3D wavelet transform. We also propose here an anisotropic version of the transform.

x y z The classical 3D dyadic wavelet analysis of L 2 (R) involves one scaling function ψ 0 ∈ L 2 (R) and a mother wavelet

H * 1 H * 1 ↓ 2 ↓ 2 H * 0 H * 1 ↓ 2 ↓ 2 H * 0 H * 0 ↓ 2 H * 1 ↓ 2 ↓ 2 H * 0 H * 1 ↓ 2 ↓ 2 H * 0 H * 1 ↓ 2 H * 0 ↓ 2 H * 1 ↓ 2 H * 0 ↓ 2 ↓ 2
ψ 1 ∈ L 2 (R) such that: ∀m ∈ {0, 1}, 1 √ 2 ψ m ( t 2 ) = ∞ k=-∞ h m [k]ψ 0 (t -k),
where (h m [k]) k∈Z are real-valued square summable sequences. The analysis filter bank presented in Fig. 3 itered on the lowpass coefficients allows to realize the mutiresolution analysis. Another specification of confocal microscopy is the difference of resolution in the z direction compared to the (x, y) directions. A way to tackle this problem would be to consider an anisotropic decomposition. By anisotropic, we mean to make a different processing (i.e. better adapted) in the z direction. We propose here to use a different wavelet in the z direction (as the decomposition is separable) and also to consider a different resolution level (lower). This is equivalent in Fig. 3 to use different filters in the z direction compared to the x and y directions and stop the multiresolution analysis earlier for the z direction.

The coefficients obtained after such an analysis will be processed as presented in the next section.

PROPOSED ALGORITHM

We study here an iterative 3D wavelet-based algorithm [START_REF] Bect | A l 1 -unified variational framework for image restoration[END_REF][START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] allowing to operate conjointly denoising and deconvolution in the wavelet transform domain, considering the 3D Discrete Wavelet Transform presented in the previous section. The objective is to:

minimize o∈H 1 2 h * o -i 2 + k∈K φ k ( o | e k ). ( 10
)
This criterion is composed of two terms: a data fidelity term and a regularization term (which penalizes each wavelet coefficients) where the regularization function φ k can be chosen adaptively.

The algorithm reads:

o n+1 = o n + λ n (prox γnφ k < o n + γ n (h * (i -ho n ))|e k > -o n ) ( 11
)
where e k stands for an orthonormal basis and prox denotes the proximity operator. Some properties of this operator as well as the convergence property of this algorithm are given in [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] .

For instance, we can choose:

φ = ω | . | (12)
or φ = ω | . | 2 ( 13
)
where ω is a fixed parameter and φ a fixed regularization function. When φ is defined by ( 12), then the corresponding prox φ is a soft thresholding with threshold ω. We take γ = 1.99 (step size) and λ = 1 (relaxation parameter) in the numerical experiments. This algorithm is directly applied to the coefficients issued from a 3D wavelet transform. When Poisson noise is considered, we apply a pre-processing on the data (Anscombe or Fisz), because this algorithm has been developed for Gaussian noise.

COMPARATIVE RESULTS

Initial definitions

In order to evaluate the quality improvement, we propose to use the I-divergence measure. The I-divergence [START_REF] Csiszár | Why least squares and maximum entropy?[END_REF] between two 3D images o and ô is given by:

I o,ô = ijk o ijk . ln o ijk ôijk -(o ijk -ôijk ) (14)

Tested images

We propose here to proceed the tests on simulated, phantom and real data.

Synthetic image

The first kind of data allows us to evaluate the proposed method performances in terms of qualitative but also quantitative measures. We suppose we know the original image and we add blur and noise in order to obtain the observed data. The original and degraded data are represented in Fig. 4. We assume the microscope has a perfect pinhole and its characteristics are the following :

• The magnification of the objective is 63×.

• The refractive index of the medium is 1.33 and of the immersion lens medium (oil) is 1.518.

• The numerical aperture is 1.4 airy units.

• The excitation wavelength is 488nm whereas the emission wavelength is 520nm.

• The voxel size is 0.02µm × 0.02µm × 0.05µm. Phantom image We consider here a spherical shell (as represented in Fig. 5) of diameter 15µm and which thickness is known to lie between 0.5µm and 0.7µm.

The acquisition is made using a microscope having the following parameters:

• The magnification of the objective is 63×.

• The refractive index of the medium is 1.33 and of the immersion lens medium (oil) is 1.518.

• The numerical aperture is 1.4 airy units.

• The excitation wavelength is 488nm whereas the emission wavelength is 520nm.

• The voxel size is 89.26nm × 89.26nm × 232.75nm.

The 3D image consists in 128 stacks of size 256 × 256 and one xy slice of the obtained image is shown in Fig. 7 (left). The objective is, knowing the original object size, to recover it after restoration process. Real image We also work with real 3D confocal microscopy images.

µm

We consider here an image which consists in 30 stacks representing a regenerating mammalian skin (see Fig. 6). The acquisition parameters are the following :

• Microscope Zeiss Axiovert 200
• The magnification of the objective is 40×.

• The refractive index of the medium is 1.33 and of the immersion lens medium (oil) is 1.518.

• The numerical aperture is 1.3 airy units.

• The excitation wavelength is 488nm whereas the emission wavelength is 520nm.

• The pinhole size is equal to 71µm

• The voxel size is 0.45µm × 0.45µm × 1µm.

Tested methods

On the other hand, when the noise is Poisson, we apply the Anscombe transform (see Section 2) in order to stabilize its variance to 1.

Quantitative and qualitative results

Synthetic image We propose, first, to give some numerical results comparing the three proposed methods. These preliminary tests allow us to calibrate the parameters. We consider symlets wavelets of length 16. For method 1, Sureshrink 24 thresholding is applied followed by the RL algorithm. For method 2, we invert the PSF as described by (9) with ǫ = 0.001, then make an estimation of the noise level (using the robust median estimator) and apply Sureshrink thresholding. Concerning method 3, φ = 0.15| . | and in the anisotropic case, we replace the z wavelet by the Haar wavelet. In all cases, we choose the Anscombe transform as pre-processing stage. As predicted, method 2 does not lead to good numerical results because it has to be better adapted. The anisotropic form of the wavelet transform allows to improve the I-divergence. The first approach also allows to obtain good results.

I-div Initial value

Phantom image

The objective here is to recover the initial shell dimension. In Fig. 7, we represent the original data (left) the restored data using method 1 (center) and the restored data using method 3 (right). It seems that method 3 leads to a reduced thickness. Let's have a look to the intensity profile represented in Fig. 8. The observed data gives approximately a thickness of 1.2µm whereas method 1 leads to 1µm and method 3 to 0.76µm. The two first measurements are too large and the third one is a little bit too large but it better approximates the reality which means that method 3 allows to better recover the initial shell size. Real image Now, we are interested in 3D real image restoration. The obtained results can be judged and evaluated visually.

In Fig. 9, we have represented two slices of the 3D image shown in Fig. 6. The right column corresponds to the restored images; the first impression can be mitigated but a careful look to these images allows to notice some improvements. In the upper left corner, we can see in the stack, the vitelline membrane which is considered as noise and alveolus structures which correspond to the central cells. In the upper right corner we have the restored image: the vitelline membrane has been removed and some central cells are recovered from noise. In the lower right corner, we can distinguish some filaments in the image which correspond to epidermal cells. In the associated restored image, a part of the noise has been removed allowing to recover some central cells and to better see epidermal cells.

CONCLUSION AND PERSPECTIVES

We have presented, in this paper, different image restoration methods for 3D confocal microscopy images based on wavelets. More sophisticated analyses such as the Dual-tree Wavelet Transform [START_REF] Kingsbury | The dual-tree complex wavelet transform : a new efficient tool for image restoration and enhancement[END_REF][START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF] could improve these restoration results at a price of a low redundancy.

Concerning the presented methods, the hyperparameters have been fixed empirically until now. An objective would be to computed them adaptively [START_REF] Jalobeanu | Hyperparameter estimation for satellite image restoration using a MCMC maximum likelihood method[END_REF] . The main idea we want to follow in the future is to develop directly Poisson statistic adapted methods.

Figure 1 .

 1 Figure 1. Scheme of Method 1.

Figure 2 .

 2 Figure 2. Scheme of Method 2.

Figure 3 .

 3 Figure 3. Analysis filter bank for a 3D Discrete Wavelet Transform.

Figure 4 .

 4 Figure 4. 3D synthetic image. Original (left) and degraded (right).

Figure 5 .

 5 Figure 5. Spherical shell of diameter 15 µm.

Figure 6 .

 6 Figure 6. 3D real image: regenerating mammalian skin. c UMR 6543 CNRS / Institute of Signaling, Developmental Biology and Cancer.

1 ,

 1 Pre-processingOn the one hand, when the noise is estimated (for example, after the PSF inversion in method 2), we use the robust median estimator which leads to the estimated variance σ of σ considering the first level detail wavelet coefficients d

Figure 7 .

 7 Figure 7. Restoration of the shell. Raw data (left), restored images using method 1 (center) and method 3 (right).

Figure 8 .

 8 Figure 8. Intensity profile of the 3 xy images respresented in Fig.7. 1 pixel represents 89nm.

Figure 9 .

 9 Figure 9. Restoration of the the regenerating mammalian skin image. c UMR 6543 CNRS / Institute of Signaling, Developmental Biology and Cancer. Original data (left column) and corresponding restored images using method 3 (right column) .

Table 1 .

 1 Restoration results (I-divergence) considering Poisson noise on the synthetic image.

		Method 1 Method 2 Method 3 Method 3 anis
	6.98	1.17	2.11	1.09	0.90

ACKNOWLEDGMENTS

The authors want to thank J.-C. Olivo-Marin and B. Zhang (Pasteur Institute, France), Z. Kam (Weizmann Institute, Israel) and Arie Feuer (Technion, Israel), P. Pankajakshan (INRIA Sophia-Antipolis, France) for several interesting discussions. We would like to thank Pasteur Institute and Nice University (UMR 6543 CNRS / Institute of Signaling, Developmental Biology and Cancer) for providing us the spherical shell images and the mammalian regenerating skin images, respectively. This research was partly supported by the P2R Franco-Israeli Collaborative Program. The first author would like to thank INRIA for funding her Post-doc fellowship.


Student Member, IEEE N Pustelnik
email: pustelni@univ-mlv.fr

Member, IEEE C Chaux
email: chaux@univ-mlv.fr

Senior Member, IEEE J.-C Pesquet
email: pesquet@univ-mlv.fr

Member, IEEE C Comtat

Parallel algorithm and hybrid regularization for dynamic PET reconstruction

To improve the estimation at the voxel level in dynamic Positron Emission Tomography (PET) imaging, we propose to develop a convex optimization approach based on a recently proposed parallel proximal method (PPXA). This class of algorithms was successfully employed for 2D deconvolution in the presence of Poisson noise and it is extended here to (dynamic) space + time PET image reconstruction. Hybrid regularization defined as a sum of a total variation and a sparsity measure is considered in this paper. The total variation is applied to each temporal-frame and a wavelet regularization is considered for the space+time data. Total variation allows us to smooth the wavelet artifacts introduced when the wavelet regularization is used alone. The proposed algorithm was evaluated on simulated dynamic fluorodeoxyglucose (FDG) brain data and compared with a regularized Expectation Maximization (EM) reconstruction. From the reconstructed dynamic images, parametric maps of the cerebral metabolic rate of glucose (CMRglu) were computed. Our approach shows a better reconstruction at the voxel level.

I. INTRODUCTION

The extraction of physiological parameters at the voxel level (parametric imaging) is a challenging problem in dynamic PET imaging (space+time). In [START_REF] Nichols | Spatiotemporal reconstruction of list-mode PET data[END_REF], [START_REF] Kamasac | Direct reconstruction of kinetic parameter images from dynamic PET data[END_REF], [START_REF] Verhaeghe | Dynamic PET reconstruction using wavelet regularization with adapted basis functions[END_REF] the authors demonstrated the advantage of taking into account time-frame characteristics during the reconstruction for improving the image quality and the physiological parameter estimation. More recently, the space+time information was used by considering waveletframe decomposition [START_REF] Verhaeghe | Dynamic PET reconstruction using wavelet regularization with adapted basis functions[END_REF], [START_REF] Zhou | Iterative PET image reconstruction using translation invariant wavelet transform[END_REF], [START_REF] Sureau | Temporal wavelet denoising of PET sinograms and images[END_REF]. To efficiently deal with such representations, the criterion to minimize often involves a ℓ 1norm (regularization term) which leads to a non-differentiable problem. New convex optimization algorithms (within the class of proximal methods) have been developed in [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF] aiming at minimizing criteria which are non-necessarily differentiable. The first application combining wavelet-frame representations and proximal algorithms in dynamic PET image reconstruction was suggested by Verhaeghe et al. [START_REF] Verhaeghe | Dynamic PET reconstruction using wavelet regularization with adapted basis functions[END_REF] by considering Forward-Backward (FB) iterations [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] to minimize a criterion involving a Gaussian data fidelity term and a wavelet spline regularization. In [START_REF] Sureau | Temporal wavelet denoising of PET sinograms and images[END_REF], the authors performed Temporal Activity Curves (TAC) denoising in the sinogram space before reconstruction by using the Douglas-Rachford (DR) algorithm so as to deal with a Kullback-Leibler divergence as a data fidelity term, which is well-adapted to Poisson noise. To perform simultaneously denoising and reconstruction with wavelet regularization, we previously designed in [START_REF] Pustelnik | Adapted convex optimization algorithm for wavelet-based dynamic PET reconstruction[END_REF] a nested iterative algorithm (FB-DR algorithm). This method allowed us to improve the reconstruction quality, but its two main drawbacks were (i) a high computational time due to the required subiteration and (ii) the presence of artifacts introduced by wavelet regularization. To tackle these problems, we consider in this paper a more recent algorithm called Parallel ProXimal Algorithm (PPXA) [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] which is well-adapted to dynamic PET due to (i) its parallel structure and (ii) the possibility of dealing with hybrid regularization such as the combination of a wavelet-frame penalty and a total variation measure [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. This method was presented for image deblurring applications in the presence of Poisson noise [START_REF] Pustelnik | Parallel ProXimal Algorithm for image restoration using hybrid regularization[END_REF] and its application to dynamic PET is now investigated.

This paper is organised as follows. First, we will present the considered degradation model and the criterion we aim at minimizing. Then, convex optimization tools such as the proximity operator and PPXA will be introduced. The main technical difficulty in this work consists of computing the proximity operator of the data fidelity term. A solution to overcome this difficulty will be proposed. Finally, numerical experiments will be provided by considering simulated and real dynamic PET [ 18 F]-FDG exams.

II. PROBLEM

We consider a finite parameterization of the original image denoted by y t = (y i,t ) 1≤i≤N which corresponds to the spatial activity distribution for a time-frame t. These data are observed through the linear system associated with the matrix A = (A j,i ) 1≤j≤M,1≤i≤N where each element A j,i represents the probability for a voxel i to be detected in the tube of response (TOR) j. Moreover, during the acquisition process, data are contaminated by Poisson noise. The effect of the noise is denoted by P. The dynamic PET data z t = (z j,t ) 1≤j≤M are related to the original images by the following degradation model:

(∀t ∈ {1, . . . , T })

z t = P(Ay t ) (1)
where T denotes the number of time-frame. The reconstruction process requires to find the finite parameterization ŷ = (ŷ t ) 1≤t≤T closest to y = (y t ) 1≤t≤T from the measured data z = (z t ) 1≤t≤T .

In order to perform this task, we consider a variational formulation where the criterion is made up of a data fidelity term related to Poisson noise, a hybrid regularization term (total variation + sparsity term in the wavelet-frame domain) and a dynamic range constraint.

The most suited data fidelity term in the presence of Poisson noise is the generalized Kullback-Leibler divergence (related to the Poisson minus log-likelihood) denoted in the following by D KL (z t , •), for z t ∈ R M . Concerning the regularization term, the total variation (denoted by tv) is applied to each time-frame and the wavelet-frame regularization deals simultaneously with the whole space+time data. The total variation allows us to smooth the artifacts introduced by wavelet-frame regularization (ℓ 1 -norm). In the following, F * ∈ R N T ×L denotes a general tight wavelet-frame synthesis operator such that F * • F = νId for some ν ∈]0, +∞[. Thus, the image ȳ can be expressed as y = F * x = (F * t x) 1≤t≤T , where x ∈ R L denotes the wavelet-frame coefficients and F * t ∈ R N ×L . The convex constraint C ⊂ R N T limits the data value range of the solution (for example, it may be a positivity constraint). Below, ι C denotes the indicator function of C such that, for every y ∈ R N T , ι C (y) = 0 if y ∈ C and +∞ otherwise. The resulting minimization problem amounts to finding the wavelet-frame coefficients such that

x ∈ Argmin x∈R L T t=1 D KL (z t , AF * t x) + ϑtv(F * t x) + κ x 1 + ι C (F * x). (2)
ϑ > 0 and γ > 0 denote the regularization parameters. They allow us to adjust the impact of the total variation and sparsity penalties on the solution. The reconstructed image is ŷ = F * x.

Due to the definition of each functional, this criterion is convex but non-differentable. To solve efficiently such a minimization problem, we will consider recent convex optimization tools.

III. CONVEX OPTIMIZATION

One of the most famous convex optimization algorithm is the Projection Onto Convex Sets (POCS) algorithm [START_REF] Bregman | The method of successive projection for a common point of convex sets[END_REF] based on alternated projections. A generalization of projections was proposed by Moreau in [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF], which is known as the proximity operator. This tool is defined for a convex, lower semi-continuous and proper function ϕ : R X → ]-∞, +∞], at a point v ∈ R X and it is such that

prox ϕ : R X → R X : v → arg min u∈R X 1 2 u -v 2 + ϕ(u).
Particularly, if C is a nonempty closed convex set of R X , prox ιC reduces to the projection P C onto C.

In the recent literature, proximal algorithms (that is based on proximity operators) have been proposed to minimize efficiently non-differentiable convex criteria. One of the most famous ones is the FB algorithm [14] also known as thresholded Landweber or iterative soft-thresholding [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] when a ℓ 1 + ℓ 2 criterion is considered. More generally, the FB algorithm allows us to minimize a sum of two functions one of which is Lipschitz differentiable. This framework is however too restrictive to solve the minimization problem (2) due to the non-Lipschitz differentiability of D KL [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF]. Another proximal algorithm known as PPXA [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] allows us to relax the Lipschitz assumption and to solve any minimization problem formulated as

x ∈ Argmin x∈R L Q q=1 f q (x) (3)
where Q ∈ N * and where (f q ) 1≤q≤Q are convex, lower semicontinuous and proper functions from R L to ]-∞, +∞]. The iterations are recalled in Algorithm 1.

Algorithm 1 General form of PPXA 1) Set γ ∈ ]0, +∞[. 2) For every q ∈ {1, . . . , Q}, set (ω q ) 1≤q≤Q ∈]0, 1] Q such that Q q=1 ω q = 1. 3) Set (u q,0 ) 1≤q≤Q ∈ (R L ) Q and x 0 = Q q=1 ω q u q,0 . 4) For n = 0, 1, . . .             4a) For q = 1, . . . , Q ⌊ p q,n = prox γfq u q,n 4b) p n = Q q=1 ω q p q,n 4c) Set λ n ∈ ]0, 2[ 4d) For q = 1, . . . , Q ⌊ u q,n+1 = u q,n + λ n (2 p n -x n -p q,n ) 4e) x n+1 = x n + λ n (p n -x n )
The convergence of the sequence (x n ) n∈N generated by Algorithm 1 to a solution x of (3) is established under the following assumption.

Assumption III.1 (i) lim x →+∞ f 1 (x) + • • • + f Q (x) = +∞. (ii) Q q=1 rint dom f q = ∅. 1 (iii) n∈N λ n (2 -λ n ) = +∞.
The main difficulty of this algorithm lies in Step 4a) which requires to compute the proximity operator of each function f q .

Fortunately, for a large class of functions, this operator takes a closed form [14]. For example, the ℓ 1 -norm proximity operator corresponds to a soft-thresholding. The way to proceed with tv is detailed in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]: it results that tv can be split in a sum of four functions, denoted by (tv i ) 1≤i≤4 when Haar filters are considered. However, a remaining difficulty in using this algorithm for dynamic PET reconstruction is to obtain an explicit form at Step 4a) for D KL (z t , A•). Subsequently, we assume that for every r ∈ {1, . . . , R}, I r is a partition of {1, . . . , M } nonempty sets. Additionally, we can write D KL (z t , •) as a sum of M functions Ψ j,t : R → ]-∞, +∞] operating on each j-th component z j,t of z t . By regrouping the terms in this sum such that j ∈ I r (for every r ∈ {1, . . . , R}, Ψ (r) t (z j,t ) j∈Ir = j∈Ir Ψ j,t (z j,t )), D KL (z t , A.) can be decomposed as a sum of R functions involving reduced-size linear operators A (r) = (A j,i ) j∈Ir,1≤i≤N containing non-overlapping and thus orthogonal rows of A. For every r ∈ {1, . . . , R}, I r denotes a set of TOR involving distinct pixels. The resulting PPXA iterations are detailed in Algorithm 2.

Algorithm 2 PPXA for dynamic PET reconstruction 1) Set γ ∈ ]0, +∞[. 2) Set ω 1 = . . . = ω R = 1 4R , ω R+1 = . . . = ω R+4 = 1 16
and ω R+5 = ω R+6 = 1 4 . 3) Set, for every q ∈ {1, . . . , R + 6}, u q,0 = F z and x 0 = R+6 q=1 ω q u q,0 . 4) For n = 0, 1, . . .

                              4a) For t = 1, . . . , T              For r = 1, . . . , R ⌊ p r,t,n = prox γ ωq Ψ (r) t •A (r) •F * t u r,n p R+1,t,n = prox γ ω R+1 tv1•F * t u R+1,n p R+2,t,n = prox γ ω R+2 tv2•F * t u R+2,n p R+3,t,n = prox γ ω R+3 tv3•F * t u R+3,n p R+4,t,n = prox γ ω R+4 tv4•F * t u R+4,n p R+5,n = prox γ ω R+5 • 1 u R+5,n p R+6,n = prox ιC •F * u R+6,n 4b) p n = R+6 q=1 ω q p q,n 4c) Set λ n ∈ ]0, 2[ 4d) For q = 1, . . . , Q ⌊ u q,n+1 = u q,n + λ n (2 p n -x n -p q,n ) 4e) x n+1 = x n + λ n (p n -x n )

IV. MATERIALS AND METHODS

A. Simulated data

The method was evaluated on simulated 2D+time data, using one slice of the Zubal brain phantom [START_REF] Zubal | Computerized three-dimensional segmented human anatomy[END_REF]. A two compartment model with k 4 = 0 was used to simulate a dynamic [18F]-FDG brain study, with 16 time frames with a duration varying between 30 seconds for the first time-frames to 5 minutes for the last ones. Moreover, small local variations of K i were created in the cortex and putamen (see Figure 2, top-left). 288 (radial) × 144 (angles) sinograms with a radial sampling of 2.247 mm were simulated by analytically projecting the phantom and adding Poisson noise. Attenuation, random and scattered coincidences were not simulated in this example. The number of events varies from 48 for the first time-frame to 26804 for the last time-frame.

The PPXA algorithm is implemented in parallel with the OpenMP library and run over 400 iterations of 10s each. The relaxation parameter is here equal to λ n ≡ 1 and the stepsize is γ = 6.10 5 . The system matrix is split in R = 432 subsets and the function number is Q = R + 6. The chosen separable orthonormal wavelets correspond to symmlet filters of length 6 over 3 resolution levels of spatial decomposition. PPXA 1 refers to our approach without temporal regularization and PPXA 2 considers a two-level decomposition with Daubechies-6 wavelets on the interval for temporal decomposition [START_REF] Cohen | Wavelets on interval and fast wavelet transforms[END_REF]. The latter choice is motived by the small number of time-frames and aims at avoiding boundary effects. The constraint on the data range is C = [0, 10 5 ] N T Bq/cc and the chosen regularization parameters are κ = 0.1 and ϑ = 10 -5 . A limited value of the parameter ϑ was chosen so as to avoid "cartoon effects". The proposed solution without/with temporal regularization in PPXA 1/PPXA 2 is compared with EM-ML and EM-ML with post-smoothing (Sieves). SIEVES 1 denotes a Gaussian filtering optimized to yield the minimum MSE and SIEVES 2 uses a Gaussian kernel leading to a similar noise level than our approach. For EM-ML and SIEVES 1/2, the iteration number is 250. A full width at half maximum (FWHM) of 12 mm x 12 mm (resp. 14 mm x 14 mm) is used for SIEVES 1 (resp. SIEVES 2).

B. Pre-clinical data

We also processed real data resulting from a [ 18 F]-FDG exam on the ecat HR+ scanner (Siemens medical solutions, Knoxville, TN, USA). 284 MBq were injected to a 4.5 kg baboon. The injection was done under camera. A dynamic acquisition was performed with septa extended. The prompt and delayed coincidences were registered separately in 16 frames with a duration varying from 1.5 minutes to 4 minutes. The global exam duration was 54 minutes. Transmission was realized before injection by a rotative source of germanium 57.

We reconstructed prompt data by considering three different approaches: filtered backprojection (FBP), EM-ML reconstruction with Gaussian (12 mm x12 mm FWHM) postfiltering (Sieves) and our approach (PPXA). During the reconstruction process, attenuation effects, normalization effects, and random coincidences have been taken into account. However, scattered coincidences have not been considered (this choice is justified by the septa extended). The system matrix used for SIEVES and PPXA does not model resolution (PSF) effects.

V. RESULTS

A. Simulated data

Figures 1 and 2 display the reconstruction results for different temporal-frames (t = 4: 1.5-2 min and t = 7: 8-13 min). During first frames when the noise level is particularly high, PPXA 2 restores quite well the brain activity even if artifacts can be observed (for example in the thalamus area). Table I compares the normalized MSE for different areas of interest such as the cortex, the thalamus and the striatum. It can be observed that PPXA 2 leads to the smallest normalized MSE for the first frames and equivalent normalized MSE for the last ones. Figures 3,4, and 5 represent Temporal Activity Curves for individual pixels in different area of interest. This allows us to evaluate the impact of the spatio-temporal reconstruction and the good results obtained with PPXA 2. Using Pmod software (PMOD Technologies Ltd, Zürich, CH), parametric images of CMRglu were computed from the EM-ML, Sieves, and PPXA reconstructions. The results are shown in Figure 6. The differences between the CMRglu images is not as much important as the TAC results or the reconstructed images. This is mainly due to the Patlak model which essentially uses the last frames for CMRglu computation, which are less sensitive to noise.

B. Pre-clinical data

Figure 7 shows the reconstructed images obtained with the different reconstruction methods for temporal-frames t = 4 (4-6 min), t = 7 (14-18 min), and t = 14 (42-46 min). These results illustrate the good behaviour and the quantitative correctness of the proposed method on real data.

VI. CONCLUSION

We have presented a convex optimization approach to deal with spatio-temporal PET reconstruction. The minimized criterion is composed with different terms: generalized Kullback-Leibler divergence (minus Poisson log-likelihood) term, ℓ 1norm for wavelet-based regularization, total variation and dynamic range constraint. To minimize efficiently this criterion, we employed PPXA which is a recent algorithm belonging to the class of proximal methods. The convergence of this algorithm is guaranteed and it is implementable in parallel. The effectiveness of the proposed approach was demonstrated on simulated and real pre-clinical data.

Fig. 1 .

 1 Fig. 1. Reconstructed images (in Bq/cc) for t = 4.

Fig. 3 .

 3 Fig. 3. TAC in the thalamus.

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. TAC in the cortex.

TABLE I MSE

 I FOR DIFFERENT AREAS.

Let ϕ : R X → ]-∞, +∞] be a convex, l.s.c. and proper function. dom ϕ = {u ∈ R X | ϕ(u) < +∞} and the interior of C ⊂ R X relative to its affine hull is the set denoted by rint C.