Saima Ben

Hadj -Laure Blanc-Féraud -Gilles Aubert

Saima Ben Hadj

Laure Blanc-Féraud

Gilles Aubert

Space Variant Blind Image Restoration

Keywords: Blind restoration, JMAP, space-variant PSF, SGP JMAP, PSF spatialement variante, restauration aveugle, SGP

In this report, we are interested in blind restoration of optical images that are degraded by a space-variant (SV) blur and corrupted with Poisson noise. For example, blur variation is due to refractive index mismatch in three dimensional fluorescence microscopy or due to atmospheric turbulence in astrophysical images. In our work, the SV Point Spread Function (PSF) is approximated by a convex combination of a set of space-invariant (SI) blurring functions. The problem is thus reduced to the estimation of the set of SI PSFs and the true image. For that, we rely on a Joint Maximum A Posteriori (JMAP) approach where the image and the PSFs are jointly estimated by minimizing a given criterion including l 1 and l 2 norms for regularizing the image and the PSFs. Our contribution is to provide a functional for the SV blind restoration problem allowing to simultaneously estimate the PSFs and the image. We show the existence of a minimizer of such a functional in the continuous setting. We describe an algorithm based on an alternate minimization scheme using a fast scaled gradient projection (SGP) algorithm. The efficiency of the proposed method is shown on simulated and real images.

Introduction

1.1 General framework and state of the art Since the first telescopes invented in the XVth century, optical imaging systems have known a considerable advance providing high resolution and good contrast photography. Nevertheless, till now they suffer from some artifacts mainly due to the inherent limitations of the optical instruments as well as the imaging environment. First, optical images are affected by a blur coming from different sources such as light diffraction, defocus, motion, etc. This resulting blur may vary according to the imaged point in the scene due to different reasons, for example, atmospheric turbulence in astrophysical imaging or moving objects in a real scene imaging. In three-dimensional (3D) fluorescence microscopy, images are affected by a variant blur due to refractive index mismatch between the different mediums composing the system. In addition to this blur, optical images are affected by random errors usually modeled by Gaussian or Poisson noise (or both). Facing such optical limits, it is essential to remove both blur and noise using numerical processing. In a mathematical point of view, the considered degradation can be modeled by the following equation:

g = N (H ( f )) (1)
where f denotes the original image that we are looking for, H is a linear blurring operator, N is the noise operator and g is the observed image. Knowing only the observation g and the noise statistic N, our goal is to recover the sharp image f . When the blurring operator H is additionally known, this inverse problem is known as image restoration. When the blur is space-invariant (SI), the operator H(.) = h * . is a convolution with a blurring function h commonly called Point Spread Function (PSF) or the system's impulse response. In this case, the considered problem is known as deconvolution. When the blurring operator is unknown, the problem is called blind restoration (BR). Blind deconvolution (BD) is thus a particular case of BR where the blur is assumed to be SI. Actually, this problem is more difficult than the non-blind case since the number of unknowns increases without any increase of the number of observed variables. The solution of such a problem is not unique and the problem is thus ill-posed. Literature on BD is wider than literature on BR. A review of the main BD methods can be found for example in [START_REF] Campisi | Blind image deconvolution: theory and applications[END_REF], [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF]) and the references therein. Despite the numerous references on that topic, BD remains a difficult under-determinate inverse problem which is still open in image processing.

RR n°8073

Nevertheless, dealing with significant space variation of the blur such as in 3D fluorescence microscopy (F. [START_REF] Gibson | Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy[END_REF], [START_REF] Pawley | Fundamental limits in confocal microscopy[END_REF]), a rigorous SV model should be used. In this case, the convolution is no longer a valid representation since the system response is different at each point. The degradation operation should be computed at each point separately, which makes computations extremely time consuming. In this context, different methods dealing with the SV BR problem have been developed. Most of them can be seen as an extension of the previous SI BD methods to the SV case. Let us recall some of these works.

In the first BR methods in a SV blur framework, the PSF is assumed to be piecewise-invariant in such a way that the image can be segmented into subimages where the blur is assumed to be invariant, see [START_REF] Lagendijk | Block-adaptive image identification and restoration[END_REF], [START_REF] Yp Guo | Blind restoration of images degraded by spacevariant blurs using iterative algorithms for both blur identification and image restoration[END_REF]). Then, the blur is locally estimated by applying SI blur identification methods. However, this usually leads to domain decomposition artifacts.

Recent methods consider a smoothly varying blur where the SV PSF is modeled by a combination of SI ones. Two different approximation models of the SV PSF have been developed in the literature, one in [START_REF] Preza | Depth-variant maximum-likelihood restoration for three-dimensional fluorescence microscopy[END_REF], [START_REF] Hirsch | Efficient filter flow for spacevariant multiframe blind deconvolution[END_REF], [START_REF] Harmeling | Space-variant single-image blind deconvolution for removing camera shake[END_REF]) and the other in [START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF], [START_REF] Bardsley | Blind iterative restoration of images with spatially-varying blur[END_REF]). They essentially differ in the way of interpolating the SI PSFs. Study and comparison of these two modeling classes are discussed in [START_REF] Denis | Fast model of space-variant blurring and its application to deconvolution in astronomy[END_REF], Ben [START_REF] Hadj | Modeling and removing depth-variant blur from 3D fluorescence microscopy[END_REF]). Nevertheless, none of these aforementioned methods automatically estimate the SI PSF locations or the combination functions used to compute the SV PSF from the SI PSFs. They usually assume that combination functions and SI PSF locations are known and that the accuracy of the results can be enhanced by increasing the SI PSF number. In our knowledge, the only fully automatic work is the one presented in [START_REF] Bar | Restoration of images with piecewise spacevariant blur[END_REF]). However, in that method, it is assumed that there are only two different blurs and that the blur difference is very significant. In our work, both abrupt and smooth blur variations can be estimated. We also assume that the SI PSF number is known. The SI PSF locations and combination functions are thus fixed.

Other works consider a continuously varying PSF as in [START_REF] Moritz Blume | A new and general method for blind shift-variant deconvolution of biomedical images[END_REF]). In other words, the PSF is estimated at every point of the image. This actually leads to a very accurate result at the expense of a huge computational time. The practical application of these methods, especially in 3D imaging, remains thus prohibitive.

In essence, all existing BR methods differ in the way of regularizing the SV PSF and the image in order to compensate for the severe lack of information en-countered in the BR problem. For example in (You andKaveh [1996, 1999]), the authors consider a piecewise smoothness of both the image and the PSF through an anisotropic diffusion based regularization. In [START_REF] An Rajagopalan | MRF model-based identification of shiftvariant point spread function for a class of imaging systems[END_REF]), the SV PSF is parameterized by a single blur parameter, smoothness constraints on the spatial variations of the blur parameter is also incorporated using Markov Random Field (MRF) model. In [START_REF] Yp Guo | Blind restoration of images degraded by spacevariant blurs using iterative algorithms for both blur identification and image restoration[END_REF]), an autoregressive (AR) model and a Tichonov term are used for the image regularization. In [START_REF] Bar | Restoration of images with piecewise spacevariant blur[END_REF]), the Mumford-Shah regularization for the image is used for the SV deblurring problem.

Problem and main contribution

In this report, we propose a solution for the BR problem of optical images that are degraded by a SV blur and corrupted with Poisson noise. In order to avoid the extensive computations due to the SV blur, we approximate the SV PSF by a convex combination of a set of SI PSFs that we denote by h 1 , ..., h M with M ≥ 1 the number of PSFs. That is, we consider an expression of the SV blurred image of the following form:

H ( f ) = H f , h 1 , ..., h M (2)
An explicit expression of H is given in sub-section 2.2. Thanks to this approximation, the blind restoration problem boils down to the estimation of the image f and the set of SI PSFs h 1 , ..., h M . For that, we rely on a Joint Maximum A Posteriori (JMAP) approach where an estimate of f , h 1 , ..., h M is given by minimizing a criterion of the following from:

J f , h 1 , ..., h M = J d f , h 1 , ..., h M + J p f ( f ) + J p h h 1 , ..., h M (3)
J d f , h 1 , ..., h M is a data fidelity term arising from the considered image formation model. J p f ( f ) and J p h h 1 , ..., h M are respectively prior terms on the image and the SI PSF set resulting from a given probability distribution. In our work, we use l 1 norm to built J p f ( f ) and l 2 norm to built J p h h 1 , ..., h M . These norms are applied on the gradient of the variables in order to promote smoothness on the solution: l 1 norm is applied on the gradient of the image in order to obtain an edge preserving smoothing and l 2 norm is applied on the gradient of the PSF in order to capture as much as possible blur in the PSF and not in the recovered image.

After defining the criterion to be minimized, we show the existence of a solution to such a minimization problem in a given set of functions that we define in the following section. To the best of our knowledge, the present report is the first one that provides such a theoretical proof. In fact, all existing BD methods based on a similar approach either attempt to minimize the criterion assuming the existence of such a minimum, either assume some hypothesis on the objective function to be minimized. For example, in [START_REF] Bolte | Alternating proximal algorithm for blind image recovery[END_REF]) a coercivity condition on a functional similar to (3) is assumed for the minimization problem. Note that in that work of [START_REF] Bolte | Alternating proximal algorithm for blind image recovery[END_REF]), a Gaussian noise model is considered. Furthermore, we derive an algorithm for the numerical resolution of the considered problem, an algorithm which is appropriately designed to overcome the computational time limit. To this end, we use a scaled gradient projection (SGP) algorithm as proposed in [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF], [START_REF] Zanella | Efficient gradient projection methods for edge-preserving removal of poisson noise[END_REF]) in an alternate way on the image f and the SI PSFs h 1 , ..., h M . In fact, SGP algorithm is known to provide fast convergence by approximating Newton method. We show on 2D and 3D simulated and real images that the SV blur model as well as the proposed algorithm are efficient for solving the considered BR problem. In particular, we apply our method on digital camera images where the blur is assumed to be SV due to spherical aberrations coming from the non-ideal shape of lenses (i.e. difference of the lens shape between the center and the edges). We also apply our method on 3D confocal microscopy images where the blur is varying along the z-axis (or the depth of the imaged point). Comparisons with the Maximum Likelihood Expectation-Maximization (MLEM) [START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF]) embedded with the alternate minimization procedure (AM), are also reported to assess the accuracy as well as the computational gain provided by the proposed method.

We now outline the contents of this report. Section 2 is devoted to notations, definitions, preliminary lemmas and problem formulation. In section 3, we prove the existence of a minimizer of a functional that we define in (14). In section 4, we describe the proposed SGP based optimization method embedded with the AM algorithm for solving the considered problem. Experimental tests on simulated and real images are reported and discussed in section 5. The report concludes in section 6 with a summary and some perspectives of our research work.

Problem formulation

Let O and I be two bounded sets of R n (with n = 2, 3 the image dimension), standing for the object space (i.e. true scene) and image space (i.e. observed scene) respectively. Let f ∈ L 1 (O) and k ∈ L 1 (I ) be two functions modeling respectively the original and blurred images. We denote by h ∈ L 1 (I × O) the SV PSF i.e. h(., u ) refers to the system response to a point source placed at u , a location in the object space O. Neglecting the noise effect, the SV blurred image can be modeled by the following equation:

k (u) = O h u, u . f u du , ∀u ∈ I (4)
In spite of the accuracy and the effectiveness of model ( 4), we cannot use it in the estimation process because it is prohibitive in terms of computational time.

Approximation of this model is needed in order to avoid such extensive computations. For that, we first recall some preliminary notations and introduce some necessary lemmas which will be used throughout this report. Then, we present the proposed approximate SV blur model based on a linear interpolation of a set of SI PSFs. Afterwards, relying on a JMAP approach, we focus on inverting that model by defining a criterion to be minimized.

Notations and preliminary lemmas

Let I = [-1, 1] n , n = 2, 3, be the support of the observed 2D or 3D image, and B(O, R) be a closed disc or sphere centered on O, and of radius R ∈ [0, 1] such that B(O, R) I . Consider the following spaces:

• H 1 (I ) = h ∈ L 2 (I ) , ∇h ∈ L 2 (I ) n • H = {h ∈ H 1 (I ) ; h ≥ 0; I h(u)du = 1; supp(h) ⊂ B(O, R)} where supp(h)
refers to the support of function h. When I h(u)du = 1, the PSF h is said to be normalized on I .

• F = { f ∈ BV (R n ); f ≥ 0; I -periodic} where BV (R n
) is the space of functions of bounded variations, see [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF]).

•

F c = { f ∈ F ; I f (u)du = c}, where 0 < c < +∞ is a positive constant. RR n°8073
In the case of SI blur, the degraded image observed on a bounded set I , is assumed to be obtained by convolving the image f with a SI PSF that we denote by h. In a practical point of view, the convolution is computed in the Fourier domain. Consequently, the image f is assumed to be infinite and periodic of period I . That is why we consider I -periodic functions in F . Then, the convolution operation

h * f (u) = R n h(u -u ) f (u )du , ∀u ∈ I = R n f (u -u ) h(u )du , ∀u ∈ I (5)
is defined for f ∈ F and h ∈ H 1 (I ).

We consider in this work a photon conservative optical system. That is, the number of photons entering the optical system is the same as the number of photons leaving the optical system. This property is known as flux conservation. We prove that it is equivalent to PSF normalization, first in the SI case, then in the SV case for the considered approximate SV blur model. In fact, these two latter properties (i.e. flux conservation and PSF normalization) are important for reducing the ambiguity due to the ill-posedness of the problem. For example, one can note that by imposing PSF normalization, one can discard many possible PSF solutions obtained by rescaling the true one.

Lemma 1. Let f ∈ F and h ∈ H 1 (I ) , supp (h) ⊂ B(O, R), we have:

I f * h (u) du = I f (u) du ⇐⇒ B(O, R) h (u) du = 1 (6)
Proof. First, we need to prove the following lemma. Lemma 2. Let f ∈ L 2 (I ), and I -periodic, then we have the following result:

I -{z} f (u)du = I f (u)du, ∀z ∈ R n (7)
Proof. Let prove the lemma 2 in the 1D case. The generalization to the n-dimensional case is straightforward. Consider I = [-1, 1], z ≥ 0 and f is 2-periodic. We have

I -{z} f (u)du = 1-z -1-z f (u)du = -1 -1-z f (u)du + 1 -1 f (u)du - 1 1-z f (u)du (8) RR n°8073
The first term of this latter equation can be expressed as:

-1 -1-z f (u)du = 1 1-z f (y- 2)dy = 1 1-z f (u)du.
Replacing this latter result in (8), we get:

1-z -1-z f (u)du = 1 -1 f (u)du.
For z ≤ 0, the proof is similar. Now, we return to the proof of the equivalence between flux conservation and PSF normalization in the SI case:

I f * h (u) du = I R n h u -u f u du du = I B(O, R) h (z) f (u -z) dz du = B(O, R) h (z) I f (u -z) du dz = B(O, R) h (z) I -{z} f (u) du dz
Using lemma 2, we obtain:

I f * h (u) du = B(O, R) h (z) dz I f (u) du
From the above equation, it is obvious the equivalence between flux conservation and PSF normalization on B(O, R).

Remark. To reduce boundary artifacts in image restoration, usually the image is extended by using a reflected version of it or simply by padding to it zeros (if the image boundaries are zeros) [START_REF] Aghdasi | Reduction of boundary artifacts in image restoration[END_REF]). The new domain Ω on which the image is defined, is thus more extended than I (i.e. I ⊂ Ω). In this case, the convolution computed using Fourier transform is periodic of period Ω. Therefore, the flux conservation no longer holds on the domain I , but on the domain Ω. period considered in the circular convolution.

SV blur model

Approximate SV blur models based on a linear interpolation of a set of SI PSFs were previously studied and compared in (Denis et al. [2011], Ben Hadj and[START_REF] Hadj | Modeling and removing depth-variant blur from 3D fluorescence microscopy[END_REF]). In our work, we consider an efficient approximation similar to that proposed in [START_REF] Preza | Depth-variant maximum-likelihood restoration for three-dimensional fluorescence microscopy[END_REF], [START_REF] Hirsch | Efficient filter flow for spacevariant multiframe blind deconvolution[END_REF], [START_REF] Harmeling | Space-variant single-image blind deconvolution for removing camera shake[END_REF]) given by the following equation:

k (u) = ∑ 1≤i≤M h i * ψ i . f (u) , ∀u ∈ I (9)
where k stands for the blurred image, f ∈ F , h i ∈ H , 1 ≤ i ≤ M is a set of SI PSFs, and

ψ i ∈ H , such that ψ i : R n → [0, 1], 1 ≤ i ≤ M is a set of weighting func- tions verifying ∑ 1≤i≤M ψ i (u ) = 1, ∀u ∈ R n .
Each weighting function is associated with a SI PSF allowing to alleviate transitions between neighboring SI PSFs. According to the variations of the considered weighting functions ψ i , i = 1, ..., M, the transitions between the different SI PSFs is more or less smooth. The expression of the SV PSF associated with model ( 9) is shown to be as follows:

h u, u = ∑ 1≤i≤M ψ i (u )h i (u -u ), ∀u ∈ I , u ∈ R n (10)
An other model, proposed in [START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF], [START_REF] Bardsley | Blind iterative restoration of images with spatially-varying blur[END_REF]) consists in interpolating the SI PSFs in the image domain i.e. ψ i vary with the image coordinates or u in our notations. It is proved in [START_REF] Denis | Fast model of space-variant blurring and its application to deconvolution in astronomy[END_REF], Ben Hadj and Blanc-Féraud [2012]) that such an interpolation is less realistic and less efficient than model (9). In fact, we can easily verify that the SV PSF given by ( 10) is normalized and the flux conservation property is satisfied when the SI PSFs h i , i = 1, ..., M are normalized (see lemma 3 and 4). Whereas, this is not the case for the other model of [START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF], [START_REF] Bardsley | Blind iterative restoration of images with spatially-varying blur[END_REF]). Consequently, in the sequel of this report, we consider the approximate degradation model given by (9). Lemma 3. Let f ∈ L 2 (I ), h i ∈ H , 1 ≤ i ≤ M, and B(u , R) a bounded disc or sphere of radius R ∈ [0, 1], centered on u ∈ R n , we have the following result:

B(u , R) ∑ 1≤i≤M ψ i (u )h i (u -u )du = 1, ∀u ∈ R n (11) Proof. B(u , R) ∑ 1≤i≤M ψ i (u )h i (u -u )du = ∑ 1≤i≤M ψ i (u ) B(u , R) h i (u -u )du = ∑ 1≤i≤M ψ i (u ) B(O, R) h i (t)dt = ∑ 1≤i≤M ψ i (u ) = 1 RR n°8073 Lemma 4. Let f ∈ L 2 (I ), h i ∈ H , 1 ≤ i ≤ M, we have: I ∑ 1≤i≤M h i * ψ i . f (u) du = I f (u) du (12)
Proof. We have:

I ∑ 1≤i≤M h i * ψ i . f (u) du = ∑ 1≤i≤M I h i * ψ i . f (u) du
As h i , 1 ≤ i ≤ M are normalized, we obtain from lemma 1:

∑ 1≤i≤M I ψ i (u) . f (u) du = I ∑ 1≤i≤M ψ i (u) . f (u) du = I f (u) du

JMAP criterion

After defining the SV blur model and checking some necessary properties, we now introduce an objective functional allowing to solve the SV BR problem. In addition to the SV blur, we take into account two kinds of noise: First, the background noise usually appears in astrophysical and microscopy images. For example, in fluorescence microscopy, the background noise is a result of the mounting media either auto-fluorescing or reflecting and scattering the imaged specimen's fluorescence. Similarly, in astronomy, the background noise comes from the diffusion of the incoming light from nearby sources from the sky or from the telescope itself. In fact, even if no visible astronomical objects are present in given part of the sky, there always is some low luminosity present, due to light diffusion from the atmosphere. comes from an error in the calibration of the imaging system so that there is usually an offset in the detector that can be assimilated to a non-negative constant b g > 0. This kind of noise can be assimilated to a non-negative constant b g > 0 that can be easily estimated from a dark area of the observed image which does not contain a signal [START_REF] Pankajakshan | Blind Deconvolution for Confocal Laser Scanning Microscopy[END_REF]). Then, b g is assumed to be known in our model.

Second, a Poisson noise originates from the stochastic nature of the photon detection process at the sensor, usually appears in astronomical and confocal mi-RR n°8073 croscopy images. Consequently, the observed image denoted by g follows a Poisson statistic of the following mean:

g ∼ P ∑ 1≤i≤M h i * ψ i . f (u) + b g (13)
For a fixed set of weighting functions, our goal is thus to jointly estimate both the SI PSF set h i , 1 ≤ i ≤ M and the true image f from the observation g. This problem is ill-posed since the solution ( f , h 1 , ..., h M ) is not unique. In particular, one trivial solution is (g, δ , ..., δ ) i.e. the SI PSFs are Dirac functions and the recovered image is the same as the degraded image. Hence, regularization should be introduced for both image and SI PSFs. In our work, we propose to use a total variation (TV) regularization for the image f [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]) and a quadratic term for each of the SI PSFs h i , 1 ≤ i ≤ M i.e. l 2 norm of the gradient of each of the SI PSFs. In other words, we are interested in minimizing w.r.t ( f , h 1 , ..., h M ) the following JMAP criterion:

J( f , h 1 , ..., h M ) = I ∑ 1≤i≤M h i * ψ i . f (u) + b g du - I g(u).log ∑ 1≤i≤M h i * ψ i . f (u) + b g du +α I |D f | + ∑ 1≤i≤M β i ∇h i 2 2 (14)
The two first terms corresponds to the data fidelity component related to the Poisson statistic [START_REF] Zanella | Efficient gradient projection methods for edge-preserving removal of poisson noise[END_REF], [START_REF] Bertero | A discrepancy principle for poisson data[END_REF]). The third term is the total variation function [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]) which allows to smooth homogeneous areas of the recovered images while preserving sharp edges. The last term, introduced for regularizing the PSFs, allows to promote the largest possible PSF and subsequently prevents a part of the blur to be associated with the recovered image. α and β i , 1 ≤ i ≤ M are regularizing weight parameters that quantify the trade off between these regularizing components and the data fidelity term. As the PSFs h i , 1 ≤ i ≤ M are different, the regularizing parameters β i , 1 ≤ i ≤ M could also be different. Furthermore, as it is stated previously, in addition to these regularizing terms, we take into account some other constraints in order to reduce the degree of freedom of this problem: On the one hand, positivity and normalization constraints are imposed on the SI PSFs. On the other hand, positivity and flux conservation constraints are imposed on the image. According to lemma 1, constraints of flux conservation and PSF normalization are redundant. However, both of them are useful for regularizing the ill-posed problem and need to be included in the estimation process (see section 4). Finally, we propose to solve the SV BR problem by solving the following constrained optimization problem:

f , ĥ1 , ..., ĥ M = Arg Min

F c ×H M J f , h 1 , ..., h M sub. to                f ≥ 0 h i ≥ 0, 1 ≤ i ≤ M h i 1 = 1, 1 ≤ i ≤ M f 1 = c (15)
where c = g 1 = I g(u)du, and F c and H M are as defined in sub-section 2.1. It is worth stressing that minimizing such a criterion is a difficult problem since it is not jointly convex w.r.t f , h 1 , ..., h M . Whereas, it can be shown that it is convex w.r.t. each of the variables separately, strictly convex if and only if the data g is strictly positive. Nevertheless, one of our contributions is to prove the existence of a minimum ( f , h 1 , ..., h M ) min of the considered objective function in the continuous setting. The proof remains also valid in the SI case. In fact, it suffices to take M = 1 and ψ 1 (u) = 1, ∀u ∈ R n .

Existence of a minimizer

Our goal is to prove the existence of ( f , ĥ1 , ..., ĥ M ) ∈ F c × H M satisfying equation (15).

According to [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF]), it suffices to show the existence a minimizing sequence

f n , h 1 n , ..., h M n n of functional J that converges to f 0 , h 1 0 , ..., h M 0 ∈ F c × H M such that: lim n→+∞ J f n , h 1 n , ..., h M n ≥ J f 0 , h 1 0 , ..., h M 0 .
For this end, some necessary lemmas should be proved:

Lemma 5. ∃a > 0 and b ∈ R such that ∀t > 0, we have:

t -g(u) log(t) ≥ at + b, ∀u ∈ I . ( 16
) RR n°8073
Proof. Let us define the following function for fixed a and b in R:

F(t) = (a -1)t + g(u)log(t) + b, ∀t > 0 (17)
The first derivative of this function is given by:

F (t) = a -1 + g(u) t . We know that F (t) ≥ 0 ⇔ g(u) t ≥ 1 -a. Consequently, for a ∈]0, 1[, we have F (t) ≥ 0 ⇔ g(u)
1-a ≥ t. In other words, F is increasing from

F(0) = -∞ to F( g(u) 1-a ) in ]0, g(u) 1-a ] and decreasing from F( g(u) 1-a ) to F(+∞) = -∞ in [ g(u)
1-a , +∞[. Furthermore, we have:

F( g(u) 1 -a ) = b -g(u) (1 + log(1 -a) -log(g(u))) .
Hence, for a ∈]0, 1[, we have:

F( g(u) 1-a ) ≤ 0, ∀u ∈ I . If we choose: b ≤ In f u∈I g(u) (1 + log(1 -a) -log(g(u)))
we have

F (t) ≤ 0, ∀u ∈ I i.e. t -g(u) log(t) ≥ at + b, ∀t > 0, ∀u ∈ I . In particular, for t = ∑ 1≤i≤M h i * ψ i . f (u) + b g , we have: ∑ 1≤i≤M h i * ψ i . f (u) + b g -g(u) log( ∑ 1≤i≤M h i * ψ i . f (u) + b g ) ≥ a ∑ 1≤i≤M h i * ψ i . f (u)) + b g + b
Corollary 1. From the previous lemma, we can deduce that:

J( f , h) ≥ (a b g + b) |I | + a I ∑ 1≤i≤M h i * ψ i . f (u) du +α I |D f | + ∑ 1≤i≤M β i ∇h i 2 2 (18)
where |I | is the measure of I .

RR n°8073

Now, let

f n , h 1 n , ..., h M n ∈ F c ×H M be a minimizing sequence of J f , h 1 , ..., h M i.e. f n , h 1 n , ..., h M n verifies: lim n→+∞ J f n , h 1 n , ..., h M n = In f F c ×H M J f , h 1 , ..., h M (19)
Moreover, assume that this infimum is finite. That is, there exists a constant

K 1 > 0 such that J f n , h 1 n , ..., h M n ≤ K 1 .
Hence, from this assumption and ( 18), we deduce the following lemma. Lemma 6. There exist constants C 1 > 0 and C 2 > 0 such that:

∇h i n 2 2 ≤ C 1 , 1 ≤ i ≤ M (20) I |D f n | ≤ C 1 (21) I h i n * ψ i . f n (u)du ≤ C 2 , 1 ≤ i ≤ M (22)
Corollary 2. For 1 ≤ i ≤ M, there exist h i 0 ∈ H and a sub-sequence of (h i n ) n which we also denote by (h i n ) n that weakly converges to h i 0 in H 1 (I ) and strongly converges to h i 0 in L 2 (I ). Proof. We know that each

h i n , 1 ≤ i ≤ M is bounded in L 1 (I ) since h i n ≥ 0 and I h i n (u)du = 1. From the Poincaré inequality, we have h i n -1 |I | I h i n (u)du L 2 (I ) ≤ C ∇h i n L 2 (I )
with C a positive constant. From the inequality (20), we deduce that

h i n - 1 |I | I h i n (u)du L 2 (I ) ≤ CC 1
Hence, h i n is bounded in L 2 (I ). Thus, there exist h i 0 ∈ H 1 (I ) and a subsequence of h i n that weakly converges to h i 0 in H 1 (I ) and strongly converges to h i 0 in L 2 (I ). Furthermore, we can easily prove that h i 0 ∈ H since we can easily verify by taking the limit that h i 0 ≥ 0,

I h i 0 (u)du = 1 and supp(h i 0 ) ⊂ B(O, R). Lemma 7. There exists a constant C 4 such that | f n | L 1 (I ) ≤ C 4 .
Proof. To show this lemma, we first need to prove the following preliminary lemma.

Lemma 8. Let h ∈ H , and v ∈ L 2 (I ), I -periodic, then we have:

|h * v| L 1 (I ) ≤ |v| L 1 (I ) (23)
Proof. Let v ∈ L 2 (I ), I -periodic, we have:

|h * v| L 1 (I ) = I |h * v(u)| du = I R 2 h(u -u )v(u )du du ≤ I R 2 h(u -u ) v(u ) du du .
Thanks to the I -periodicity of |v(u )| and lemma 1, we obtain |h * v| L 1 (I ) ≤

I |v(u)| du. Now, we focus on proving that the sequence ( f n ) n is bounded. For that, let consider the following two functions:

w n = 1 I I f n (u)du v n = f n -w n
We can easily verify that I v n (u)du = 0 and Dv n = D f n . Consequently, from (21), we obtain

I |Dv n (u)| = I |D f n | ≤ C 1 . Furthermore, from the Poincaré inequality, there exists a constant C > 0 such that v n - 1 |I | I v n (u)du L 1 (I ) ≤ C |Dv n | L 1 (I ) .
Hence, since I v n (u)du = 0, we obtain

|v n | L 1 (I ) ≤ C 3 ( 24
)
with

C 3 = CC 1 . Besides, we know that h i n * ψ i w n = h i n * ψ i f n -h i n * ψ i v n . Hence, we obtain h i n * ψ i w n L 1 (I ) ≤ h i n * ψ i f n L 1 (I ) + h i n * ψ i v n L 1 (I ) .
From inequality (22) and lemma 8, we get:

h i n * ψ i w n L 1 (I ) ≤ C 2 + |v n | L 1 (I )
. From inequality (24), we obtain for all i ∈ {1, ...., M}:

h i n * ψ i w n L 1 (I ) ≤ C 2 +C 3 = C 4 (25) RR n°8073
Besides, thanks to the periodicity of ψ i , we prove that

h i n * ψ i w n L 1 (I ) = I f n (u) du.
Indeed,

h i n * ψ i w n L 1 (I ) = I R 2 h n (u -u ) ψ i w n (u )du du.
Since h i n , ψ i and w n are positive, we obtain:

h i n * ψ i w n L 1 (I ) = I R 2 h n (u -u ) ψ i (u ) w n (u )du du = 1 |I | I f n (u )du I R 2 h n (u -u ) ψ i (u ) du du Since R 2 h i n (u -u )du = 1, ∀u ∈ I
, and ψ i is I -periodic, we obtain using lemma 1 the following result

h i n * ψ i w n L 1 (I ) = 1 |I | I f n (u) du I ψ i (u) du.
Summing over i and taking into account that

∑ 1≤i≤M ψ i (u) = 1, ∀u ∈ I , we obtain, ∑ 1≤i≤M h i n * ψ i w n L 1 (I ) = 1 |I | I f n (u) du I ∑ 1≤i≤M ψ i (u) du = I f n (u) du. Consequently, I f n (u) du ≤ C 4 i.e. | f n | L 1 (I ) ≤ C 4 . We proved that f n is bounded in BV (I ) i.e. | f n | L 1 ≤ C 4 and |D f n | L 1 ≤ C 1 .
Consequently, there exists f 0 such that f n weakly converges to f 0 in BV (I ) and strongly converges to f 0 in L 1 (I ).

Now, let us prove that f 0 ∈ F c . By taking the limit of f n , it is easy to show that I f 0 (u) du = c. It suffices to show that f 0 is I -periodic. For sake of clarity, we present the proof in the 1D case i.e. I = [-1, 1], the considered period is thus 2.

Let f 0 be the I -periodic extension of f 0 and φ ∈ L 1 (I ). We show that for each cell

I k = [-1 + 2k, 1 + 2k], k ∈ Z, f n weakly converges to f 0 in L 1 (I ) lim n→+∞ 1+2k -1+2k f n (u)φ (u)du = lim n→+∞ 1 -1 f n (u + 2k)φ (u + 2k)du = lim n→+∞ 1 -1 f n (u)φ (u + 2k)du = 1 -1 f 0 (u)φ (u + 2k)du = 1+2k -1+2k f 0 (u + 2k)φ (u)du = 1+2k -1+2k f 0 (u)φ (u)du
Thus, f n weakly converges to f 0 in L 1 (I ). In a similar way, we prove the periodicity of D f 0 .

Theorem 1. The problem In f

F c ×H M J f , h 1 , ..., h M has at least one solution.
Proof. Thanks to the previous lemma, we have that any minimizing sequence f n , h 1 n , ..., h M n n converges (up to a subsequence) to f 0 , h 1 0 , ..., h M 0 ∈ F c ×H M . We have to prove that lim inf

n→+∞ J f n , h 1 n , ..., h M n ≥ J f 0 , h 1 0 , ..., h M 0 . From the convexity of functions I |∇h i (u)| 2 du, 1 ≤ i ≤ M and I |D f (u)|, we infer that: lim inf n→+∞ I |∇h n (u)| 2 du ≥ I |∇h 0 (u)| 2 du (26) lim inf n→+∞ I |D f n | ≥ I |D f 0 |du (27)
It remains to show that

lim n→+∞ I G ∑ 1≤i≤M h i n * ψ i . f n (u) + b g du ≥ I G ∑ 1≤i≤M h i 0 * ψ i . f 0 (u) + b g du
with G(t) = tg(u)log(t), ∀t > 0 and u ∈ I .

RR n°8073

We know that G(t) has a lower bound since

G(t) ≥ u -g(u)log(u) ≥ In f u∈I u - g(u)log(u) = m. Assuming that m > -∞, the sequence G ∑ 1≤i≤M h i n * ψ i . f n (u) + b g n
has a finite lower bound. Hence, using Fatou lemma, we obtain:

lim n→+∞ I G ∑ 1≤i≤M h i n * ψ i . f n (u) + b g du ≥ I lim n→+∞ G ∑ 1≤i≤M h i n * ψ i . f n (u) + b g du That is, lim inf n→+∞ I G ∑ 1≤i≤M h i n * ψ i . f n (u) + b g du ≥ I In f G ∑ 1≤i≤M h i 0 * ψ i . f 0 (u) + b g du
Thereby, we proved that

lim inf n→+∞ J f n , h 1 n , ..., h M n ≥ J f 0 , h 1 0 , ..., h M 0
and hence we showed that

f 0 , h 1 0 , ..., h M 0 = Arg Min F c ×H M J f , h 1 , ..., h M .

Numerical minimization

We focus in this section on the numerical approximation of a solution of problem (15). For that, we use an alternate minimization (AM) procedure. Each of the problem variables is thus estimated alternately by setting the others to their RR n°8073 previous estimates. The global minimization problem is thus divided into many elementary minimization problems. To solve each of them, we use a recent scaled gradient projection (SGP) method since it is shown to provide a high convergence speed compared to other gradient based optimization methods [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]). For sake of compactness, we call this method SGP based alternate minimization algorithm and we denote it in the sequel by SGPAM. In this section, we first introduce the discrete version of the considered problem and describe the proposed AM procedure. Then, we introduce the SGP method in its general form. Afterwards, we show how we fit this latter algorithm to image and PSF estimation problems.

Discrete notations and AM scheme

First of all, we need to introduce some discrete notations. We denote vectors by bold symbols and matrices by uppercase letters. The 2D or 3D image can thus be represented by an N-dimensional vector f f f = ( f 1 , ..., f N ) T ∈ R N considering lexicographical order of pixels (or voxels), with N the image size. The SI blur operator can be represented by a matrix H ∈ R N×N which is block circulant with circulant blocks for 2D images and block circulant with circulant block circulant blocks for 3D images. The circular convolution is thus computed by the following matrix-vector multiplication H f f f . Roles can be inverted and the circular convolution can be computed as Fh h h where F ∈ R N×N is a matrix representing the image and h h h ∈ R N is N-dimensional vector representing the SI PSF. The discrete function that we are interested in, is also denoted by J f f f , h h h 1 , ..., h h h M and is expressed as follows:

J f f f , h h h 1 , ..., h h h M = 1 1 1 T ∑ 1≤i≤M H i ψ i f f f + b b b g g g -g g g T log ∑ 1≤i≤M H i ψ i f f f + b b b g g g +α ∇ f f f 1 + ∑ 1≤i≤M β i ∇h h h i 2 2 ( 28
)
where 1 1 1 ∈ R N stands for N-size vector whose components are all equal to 1, b b b g g g ∈ R N is a strictly positive constant vector (all its components are equal) modeling the background noise,

ψ i ∈ R N×N , 1 ≤ i ≤ M is a set of diagonal matrices referring
to weighting coefficients for the SV blur, and H i ∈ R N×N , 1 ≤ i ≤ M is a set of matrices modeling the SI blur operators. Numerical computations of the gradients as well as l 1 and l 2 norms used in the considered regularizing terms are given in appendix A.

The constrained optimization problem that we are interested in, is the following:

f f f , ĥ h h 1 , ..., ĥ

h h M = Arg Min ( f f f , h h h 1 , ..., h h h M ) J f f f , h h h 1 , ..., h h h M sub. to                          f f f ≥ 0 f f f 1 = c h h h i ≥ 0, 1 ≤ i ≤ M h h h i 1 = 1, 1 ≤ i ≤ M supp h h h i ⊂ B, 1 ≤ i ≤ M (29)
where c = g 1 is a positive constant and B is a given index set corresponding to the support of the SI PSFs. As it is shown in lemma 1, constraints h h h i 1 = 1 and f f f 1 = c are redundant. However, it is important to incorporate both of them in the minimization algorithm since proceed by alternate minimization does not necessarily preserve these two properties even if one of them is satisfied (see subsection 4.2). One trivial way for minimizing J is to stack the image and the PSF vectors into the same vector X and then apply any optimization method on this variable. However, this is not the optimal way for dealing with such a multivariate optimization problem. In fact, this should be extremely slow especially when using gradient-descent-based optimization methods since the gradients of the criterion w.r.t. the image and w.r.t. the PSF vectors may have different orders of magnitude. AM scheme is known to be appropriate to such a situation i.e. the objective function to be minimized has two or more unknowns with different orders of magnitude. It consists in splitting the problem into two or more stages and alternating between them. First, keeping the PSF vectors constant, estimate the image. Then, update each of the SI PSFs one after the other by fixing the remaining PSFs and the image to their previous estimates. That is to say, starting from an initial guess of the PSFs ĥ h h 1 (0) , ..., ĥ h h M (0) , the algorithm is given by the two following steps:

• Image estimation step:

f f f (k+1) = Arg Min f J f f f , ĥ h h 1 (k) , ..., ĥ h h M (k) sub. to        f f f ≥ 0, ∑ 1≤i≤N f i = c (30)
• PSF estimation step:

ĥ h h 1 (k+1) = Arg Min h h h 1 J f f f (k+1) , h h h 1 , ĥ h h 2 (k) , ..., ĥ h h M (k) sub. to              h h h 1 ≥ 0 ∑ 1≤i≤N h 1 i = 1 supp h h h 1 ⊂ B (31)
. . . k+1) , ĥ h h 1 (k+1) , ..., ĥ

ĥ h h M (k+1) = Arg Min h h h M J f f f (
h h M-1 (k+1) , h h h M sub. to              h h h M ≥ 0 ∑ 1≤i≤N h M i = 1 supp h h h M ⊂ B (32)
One can remark that at each AM step, the global energy monotonously decreases: 1) , ĥ h h 1 (1) , ..., ĥ h h M (0) ≥ . . . 1) , ĥ h h 1 (1) , ..., ĥ h h M (1) ≥ . . .

J f f f (1) , ĥ h h 1 (0) , ..., ĥ h h M (0) ≥ J f f f (
≥ J f f f (
≥ J f f f (K) , ĥ h h 1 (K) , ..., ĥ h h M (K)
where K > 1 is the number of iterations.

RR n°8073

Several convex optimization methods have been developed in the literature. In particular, gradient-search techniques such as RLTV algorithm [START_REF] Dey | Richardson-lucy algorithm with total variation regularization for 3d confocal microscope deconvolution[END_REF]) are widely used. Here, we propose to use an efficient SGP method [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]) very suitable to our constrained problem and presenting the advantage of fast convergence. In the following sub-sections, we describe this algorithm and show how we use it to solve problems (30), ( 31) and (32).

SGP algorithm

Recently, SGP method is shown to be efficient for solving convex constrained optimization problem of the following form:

Min J J J(x x x) x x x∈Ω (33)
where x x x = (x 1 , ..., x N ) T ∈ R N is an N-dimensional vector, Ω ⊂ R N is a closed convex set, and J J J : Ω → R is a differentiable function. The method was firstly applied for image deblurring by minimizing a non linear convex function arising from the maximum likelihood approach, subject to non-negativity and flux conservation constraints [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]) (no regularization term used in the objective function). Then, it was applied for Poisson noise removal by minimizing an objective function including a total variation term and the Kullback-Leibler divergence [START_REF] Zanella | Efficient gradient projection methods for edge-preserving removal of poisson noise[END_REF]). The main success of this method consists in exploiting effective scaling strategies and step-length updating rules, appropriately designed for improving the convergence rate and thus the computational time also. We limit ourselves to the presentation of the main lines of this method. Details of the algorithm can be found in [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]). The method solves problem (33) by approximating the following fixed point:

x x x * = P Ω,S (x x x * -δ S∇J J J(x x x * )) ( 34
)
where δ is a positive scalar referring to the step-length of the proposed descent method and S is symmetric positive definite N × N matrix which corresponds to the scaling matrix. (δ S) used in (34) allows to approximate the inverse of the Hessian matrix of J J J in order to enforce quasi-Newton properties and thus provide good convergence rate. Appropriate selection of these two latter parameters i.e. the scaling matrix S and the step length δ is thus an important task. In [START_REF] Zanella | Efficient gradient projection methods for edge-preserving removal of poisson noise[END_REF]), a technique for selecting these parameters is derived. In particular, the scaling matrix S is chosen from the KKT condition (see sub-sections 4.3 and RR n°8073

4.4) and the step-length δ is chosen by using an adaptive alternation of Barzilai and Borwein rules [START_REF] Barzilai | Two-point step size gradient methods[END_REF], [START_REF] Frassoldati | New adaptive stepsize selections in gradient methods[END_REF], [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]). Furthermore, the convergence of the proposed method is proved in [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]) when the elements (s i j ) 1≤i, j≤N of S verify:

1 L ≤ s i j ≤ L, L > 1 ( 35
)
P Ω,S : R N → Ω used in ( 34) is the projection operator expressed as follows:

P Ω,S (x x x) = Arg min y y y∈Ω y y yx x x S

where . S indicates the vector norm associated with the symmetric positive definite matrix S: x x x S = √ x x x T Sx x x. To apply the SGP method for image and PSF estimations, it suffices to specify the objective function J J J, its derivative ∇J J J, the scaling matric S, and the step-length δ for the considered problems. This will be the goal of the two following subsections. For further details about the algorithm and its technical implementation, we invite the reader to refer to [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]).

SGP algorithm for image estimation

We are interested in this sub-section in solving problem (30). As the PSF set is supposed to be fixed at this step, the objective function to be minimized is reduced to:

J J J 0 ( f f f ) = 1 1 1 T ∑ 1≤i≤M H i ψ i f f f + b b b g g g -g g g T log ∑ 1≤i≤M H i ψ i f f f + b b b g g g + α ∇ f f f 1 (36)
We apply on this functional the SGP algorithm with positivity and flux conservation constraints on the variable f f f by following the same steps as in [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]. Nevertheless, since we are considering a different objective functional from that used in [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]), some variables need to be expressed, namely the gradient of the functional ∇J J J 0 ( f f f ) and the scaling matrix S used in equation ( 34). The step-length δ depends on S but it can be easily updated using equations ( 29) and ( 30) of [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]). In fact, in the functional that we propose, we consider an additional regularizing term which is not used in the deblurring problem of [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]), only a similar data term including a simple blur operator is considered. In [START_REF] Zanella | Efficient gradient projection methods for edge-preserving removal of poisson noise[END_REF]), a similar regularizing term is used for image denoising but without considering any blurring operator. In the functional that we propose, we consider both of the blurring operator and the regularization term.

We now focus on expressing ∇J J J 0 ( f f f ), then we infer an expression of the scaling matrix S:

∇J J J 0 ( f f f ) = ∑ 1≤i≤M ψ i H i T 1 1 1 -Y -1 g g g + α∇J J J R 0 ( f f f ) ( 37
)
where Y is a diagonal matrix with the following entries:

Y = ∑ 1≤i≤M H i ψ i f f f + b b b g g g (38)
As the PSFs are normalized, then we can easily prove that ∑ 1≤i≤M

ψ i H i T 1 1 1 = 1 1 1.
Expression (37) can thus be simplified as:

∇J J J 0 ( f f f ) = 1 1 1 -∑ 1≤i≤M ψ i H i T Y -1 g g g + α∇J J J R 0 ( f f f ) (39)
The gradient of the regularization term ∇J J J R 0 ( f f f ) can be split into positive and negative parts:

-∇J J J R 0 ( f f f ) = U R 0 ( f f f ) -V R 0 ( f f f ) (40)
where V R 0 ( f f f ) and U R 0 ( f f f ) are non-negative vectors in R N . A similar splitting to that proposed in [START_REF] Zanella | Efficient gradient projection methods for edge-preserving removal of poisson noise[END_REF]), for 2D and 3D images is given in appendix A.

-∇J J J 0 ( f f f ) = ∑ 1≤i≤M ψ i H i T Y -1 g g g + αU R 0 ( f f f ) -1 1 1 + αV R 0 ( f f f ) (41)
Now, we infer an expression of S from the first KKT condition for a minimum f f f * of function (36), i.e. f f f * ∇J J J 0 ( f f f * ) = 0. That is, one can use the following minimizing iteration:

f f f (k+1) = Z (k) -1 f f f (k) ∑ 1≤i≤M ψ i H i T Y (k) -1 g g g + αU R 0 f f f (k) (42) RR n°8073
where Z (k) is a diagonal matrix with entries 1 1 1 + αV R 0 f f f (k) . Remark that the above equation can be seen as a regularized version of the Rishardson-Lucy (RL) or Maximum Likelihood Expectation Maximization (MLEM) algorithm [START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF]) fitted to a SV PSF. Note that this version has an important advantage w.r.t. the widely used RLTV algorithm [START_REF] Dey | Richardson-lucy algorithm with total variation regularization for 3d confocal microscope deconvolution[END_REF]). Thanks to the proposed decomposition of the gradient of the regularization term, the algorithm preserves the positivity of the solution which is not guaranteed by RLTV algorithm. Furthermore, we can prove that (42) can be expressed as:

f f f (k+1) = f f f (k) -Z (k) -1 f f f (k) ∇J J J 0 f f f (k) (43)
This iteration corresponds to the SGP algorithm for δ (k) = 1 and S (k) a diagonal matrix as follows:

S (k) = diag s s s (k) , s s s (k) = (s i (k) ) 1≤i≤N = Z (k) -1 f f f (k) (44)
As the scaling matrix should satisfy condition (35), we choose:

s i (k) = min{L, max{ 1 L , s i (k) }}, i = 1, ..., N, L > 1. ( 45
)
As it is suggested in [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]), L should be set to a high value (e.g. L = 10 10 in our tests). By computing ∇J J J 0 ( f f f ) as in ( 41) and choosing S as in ( 44) and ( 45), one can easily apply SGP algorithm on f f f .

SGP algorithm for PSF estimation

Now, to solve problems ( 31) and ( 32), we follow a similar strategy as previously.

Consider for example the optimization problem w.r.t. the PSF h h h j , j ∈ {1, ..., M}. Denoting by F i , i ∈ {1, ..., M} the matrix obtained from the vector ψ i f f f , we look for minimizing the following function w.r.t. h h h j , j ∈ {1, ..., M}:

J J J j h h h j = 1 1 1 T ∑ 1≤i≤M F i h h h i + b b b g g g -g g g T log ∑ 1≤i≤M F i h h h i + b b b g g g + β j J J J R j h h h j (46) with J J J R j h h h j = ∇h h h j 2 2
. To apply SGP algorithm on the above functional with positivity, PSF normalization, and PSF support constraints, we need to compute RR n°8073 ∇J J J j h h h j and S for this problem. Note that SGP algorithm has never been used for such a problem where the regularizing term is quadratic and the operator applied on the variable to be estimated, is not normalized (F i 1 1 1 = 1 1 1).

Let us first express the gradient of functional J J J R j h h h j w.r.t. h h h j :

∇J J J j h h h j = F j T 1 1 1 -W -1 g g g + β j ∇J J J R j h h h j (47) with W = diag ∑ 1≤i≤M F i h h h i + b b b g g g a diagonal matrix.
In an analogous way as in ( 40), we decompose ∇J J J R j h h h j into positive and negative parts (see appendix A):

-∇J J J R j h h h j = U R j h h h j -V R j h h h j
Now, we can derive an expression of S from the following iteration of the regularized RL algorithm:

h h h j (k+1) = E (k) -1 h h h j (k) F j T W (k) -1 g g g + β j U R j h h h j (k) ( 48
)
where E (k) is a diagonal matrix with entries:

F j 1 1 1 + β j V R j h h h j (k)
and

W (k) = diag F j h h h j (k) + ∑ i∈{1,...,M}-{ j} F i h h h i + b b b g g g is a diagonal matrix.
We can prove that equation ( 48) is equivalent to the following one:

h h h j (k+1) = h h h j (k) -E (k) -1 h h h j (k) ∇J J J j h h h j (49)
Consequently, an expression of S (k) = diag s s s (k) , that can be used for this problem is:

s i (k) = min{L , max{ 1 L , s i (k) }}, i = 1, ..., N, L > 1. ( 50
)
with s i (k) , i = 1, ..., N are entries of the following vector:

s s s (k) = E (k) -1 h h h j (k) and
L is a constant of high value.

RR n°8073

Experimental tests

We test our method on two kinds of data: 2D images of a digital camera and 3D images of confocal microscopy. We first validate our method on simulated data of both systems, then we apply it on real data. According to the considered application, we define appropriate weighting functions dealing with smooth or non smooth blur variations. The choice of these weighting functions was previously discussed in [START_REF] Yuan | 3D fluorescence microscopy imaging accounting for depthvarying point-spread functions predicted by a strata interpolation method and a principal component analysis method[END_REF]), where linear weighting functions as those considered in our experiments are compared with those obtained by a Principle component analysis (PCA) method [START_REF] Arigovindan | A parallel productconvolution approach for representing the depth varying point spread functions in 3D widefield microscopy based on principal component analysis[END_REF]). No very significant difference in the precision provided by both models is noticed in that report. Furthermore, in our experiments, the SGP parameters are set as in [START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]) for both image and PSF estimations: L = L = 10 10 are used in ( 45) and ( 50) to bound the scaling matrices, lower and upper bounds δ min and δ max of the step-length δ are set as follows: δ min = 10 -5 , δ max = 10 5 . The global iterative algorithm is stopped when the global energy remains unchanged (i.e. the normalized energy error between two successive iterations is less than a small threshold) or when a maximum number of iterations of the global algorithm is achieved (e.g. in our tests, the maximum number of iterations used is 1000). The maximum number of iterations of the SGP algorithm is set to 10. All our tests are performed on a machine having a processor frequency of 1.86 GHz. In order to assess the performance of the proposed BR method, we measure two different evaluation criteria:

• the relative reconstruction error defined as:

RRE f f f = f f f -f f f 2 f f f 2 (51)
where f f f refers to the image vector to be reconstructed and f f f is an estimate of this latter.

• structural similarity index (SSIM) allowing to detect changes in structures between two images [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF]), expressed as:

SSIM f f f = 2µ f f f µ f f f + c 1 2σ f f f f f f + c 2 µ f f f 2 µ f f f 2 + c 1 σ f f f 2 σ f f f 2 + c 2 (52) RR n°8073
where µ f f f and µ f f f are respectively averages of f f f and f f f , σ f f f 2 and σ f f f 2 are their variances and σ f f f f f f is the covariance of f f f and f f f . c 1 and c 2 are two constants to stabilize the division that can be set as

c 1 = (k 1 R) 2 , c 2 = (k 2 R) 2
, with R the dynamic range of the pixel-values, and k 1 and k 2 two constants that can be chosen as k 1 = 0.01 , k 2 = 0.03. More the SSIM value is close to 1, more the two image structures are similar.

In order to better evaluate our method (SGPAM), we compare it with an other AM based algorithm where elementary minimizations are performed using a regularized version of RL method (see equations ( 42) and ( 48)). We denote this latter method by RRLAM. In order to show the interest of the proposed SV blur model, we present some BD results using one SI PSF.

Test on 2D digital camera images

The first experiment is related to optical distortions in digital camera which cannot be considered as uniform, simply because an optical device (lenses, mirror, etc) is not flat and does not have the same imaging properties at the center and at the edges due to its spherical shape [START_REF] Born | Principles of optics[END_REF], [START_REF] Sivak | Spherical aberration of the crystalline lens[END_REF]). Optical deformations are not identical over its entire surface. Nevertheless, these distortions could be considered homogeneous on a given crown whose center is the center of optics. That is why we consider circular variations of the PSF in our experiments. In Fig. 1 (a), we show an example of two weighting functions used in our experiments. The dominant noise in digital camera can be considered as Gaussian because of the the high illumination. Nevertheless, our method designed for Poisson noise corrupted images can be applied on these images. In fact, for high intensity level, the Poisson distribution is very similar to the Gaussian distribution. In [START_REF] Stagliano | Analysis of an approximate model for poisson data reconstruction and a related discrepancy principle[END_REF]), the authors prove that the Poisson criterion can be approximated by the Gaussian criterion. Our tests are performed on two different simulated images of a digital camera system, with different blur variations and different noise levels. First, consider the sharp image of 256 × 256 pixels, depicted on Fig. 2 (a), that we denote by the letter "A". We blur this image with a SV PSF constructed from a combination of two Gaussian PSFs displayed on Fig. 2 (e) and (i) in a logarithmic scale. For that, we use two circular varying weighting functions as in Fig. 1 (a). Before adding the Poisson noise to this blurred image, we add to it a background constant of b g = 10 -2 . The final degraded image that we obtained, depicted in RR n°8073 Fig. 2 (b) has a peak signal to noise ratio (PSNR) of 18 dB. Starting from a random initial guess of the two PSFs which is not very far from the true ones (cf. Fig. 2 (f) and (j)), we run our SGPAM algorithm on that image. The regularizing parameters are set in an empirical way by test/error in order to visually obtain a satisfactory reconstructed image, say α = 0.005 and β i = 10 3 , i = 1, 2 for SV and SI restoration tests. Note that α and β i have an antagonist impact on the recovered image: a very high value of β i leads to a wide PSF and a sharp recovered image (ringing artifacts could appear on the recovered image) while a low value of α leads to a smooth image and a narrow PSF. The estimated image and SI PSFs are respectively depicted in Fig. 2 (c), (g) and (k). In order to evaluate the proposed method, we compare in terms of quality of reconstruction and speed of convergence with RRLAM algorithm. Constraints of positivity, flux conservation and PSF normalization are added to that algorithm by projecting each of the estimates on a convenient set. In fact, without these two constraints, the algorithm leads to disappointing results. The obtained image and PSF estimates by this latter algorithm are displayed in Fig. 2 (d), (h), and (l). In order to show the interest and the efficiency of the proposed BR method, we compare our result with that obtained by considering one single SI PSF. In table 2, we summarize the main comparison results. For the different considered tests, we present in that table the RRE and the SSIM measures between the recovered images and the original one. We also give the mean computing time per iteration (t m ), the global computing time (t g ), and the required number of iterations (nb it).

Our second test is performed on an other image simulated with a more realistic blur model. In fact, in the previous test, the PSFs used are Gaussian while in this test we use real PSFs of a digital camera system. Furthermore, in the previous test, the PSF is assumed to be the same on a given crown. However, an accurate model should take into account the variation of the blur on the same crown [START_REF] Sivak | Spherical aberration of the crystalline lens[END_REF], [START_REF] Cox | Soft contact lens-induced longitudinal spherical aberration and its effect on contrast sensitivity[END_REF]).

Let consider the original image of Fig. 3 (a) of size 500 × 500 pixels. We denote this test image by letter "B" in table 2. We decompose it into three crowns, then each of these crowns unless the central one is sectioned into four pieces. In Fig. 1 (b), we present such a decomposition. As the blur variation between two neighboring regions is not significant, we consider a SI PSF at each region. That is, weighting functions are chosen to be indicator functions (ψ i = 1 on the considered region and ψ i = 0 outside that region). Note that it is possible to chose other weighting functions, for example angular and radial varying functions in order to take into account the transitions between neighboring PSFs. Having nine measured PSFs, the distorted image is obtained by applying each of these PSFs on the corresponding region, adding a background constant of b g = 10 -5 and finally adding the Poisson noise. The resulting image depicted in Fig. 3 (b) has a PSNR of 29 dB. Restoration results using SGPAM and RLAM are respectively depicted in Fig. 3 (c) and (d). The used regularizing parameters are set in an empirical way by test/error: α = 0.05 and β i = 1, i = 1, ..., 9 for both SV and SI tests.

Test on 3D confocal microscopy images

Confocal laser scanning microscopy (CLSM) [START_REF] Inoué | Foundations of confocal scanned imaging in light microscopy[END_REF], [START_REF] Pawley | Fundamental limits in confocal microscopy[END_REF]) allows to observe a 3D biological living specimen at a resolution of about 100 nm. However, 3D images coming from this system are distorted by a blur varying along the z-axis i.e. the depth of the imaged point in the specimen. This varying blur is basically due to light refraction phenomenon when crossing mediums of different refractive indexes in this system. As the system provides threedimensional images, the PSF, the image of a point source is thus 3D. Consequently, we deal with 3D image restoration. The PSF variation mainly manifests by three aspects when increasing the depth: radial and axial spread of the 3D PSF, axial dissymmetry, and axial shift of the main lobe w.r.t the central plane (see Fig. 4 (e) and (i) which show axial slices of two PSFs at 0 µm and 14.5 µm of depths respectively). Knowing all the imaging parameters such as light wavelength, refractive indexes of the different system mediums, the PSF at different depths can be computed using the theoretical model developed in [START_REF] Gibson | Diffraction by a circular aperture as a model for threedimensional optical microscopy[END_REF], F. [START_REF] Gibson | Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy[END_REF]). Nevertheless, some of these parameters are inaccessible, especially the refractive index (RI) which is very sensible to changes of the temperature. That is why a blind or semi-blind restoration method need to be developed for this system.

In what follows, we show test results on both simulated and real CLSM images. The simulated CLSM image is generated using the theoretical PSF model presented in (F. [START_REF] Gibson | Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy[END_REF]).

Test on a simulated CSLM image

Consider a 3D image of 100 × 100 × 100 voxels, of three spherical beads, whose axial slice is depicted in Fig. 4 (a). We denote this image by letter "C". The beads are assumed to be embedded in a medium of RI n s = 1.48. The CLSM imaging system assumed to have a magnification of 100X, a numerical aperture of 1.4 and an oil immersed lens with a RI of n i = 1.5. The cover-slip chosen to have a RI very close to that of the objective immersion medium (n c 1.5), so that aberrations induced by the RI mismatch between these two mediums are negligible. The excitation and emission wavelengths are respectively assumed to be 560 nm and 600 nm. The pinhole of the confocal microscope is considered to be very small so that it can be approximated by a Dirac function in the PSF generation model of (F. [START_REF] Gibson | Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy[END_REF]). Radial and axial sampling steps are respectively set to 50 nm and 145 nm, respecting Nyquist sampling. Knowing all these acquisition parameters, we generate 100 PSFs, each at a different depth, using the theoretical PSF model in (F. [START_REF] Gibson | Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy[END_REF]). Then, we generated the blurred image using a discrete version of model ( 4). Adding to it a background noise of bg = 10 -4 and Poisson noise, we obtain the simulated observation whose axial slice is presented in Fig. 4 (b) with a PSNR of 13 dB.

In order to restore that image, we approximate the depth-variant PSF by a combination of two SI PSFs, taken at the top and the bottom of the sample. The considered weighting functions are varying along the z-axis, and constant along x, y axis, cf. Fig. 1 (c) that shows a plot along z of the used weighting functions. The two SI PSFs to be estimated are initialized using the theoretical model (F. [START_REF] Gibson | Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy[END_REF]), parameters included in this model are set in an approximate way (some random errors are introduced in these parameters). The regularizing parameters are set as follows α = 10 -3 and β i = 0.1, i = 1, 2 for all tests. The estimated PSFs as well as the recovered image using SGPAM and RLAM methods are respectively shown in Fig. 4 (c), (d), (g), (h), (k), and (l).

In order to better evaluate the proposed BR method, we show in Fig. 5 (a) the plots of intensity profiles along the z-axis passing through the centers of the original beads, the observed beads, and the restored beads with SGPAM and RLAM methods. We also show in Fig. 5 (b) and (c) the intensity profiles along the z-axis of the true, initial, and final estimated PSFs.

Test on a real CSLM image

Consider a real image of a bead shell of diameter 5 µm, embedded in a polymer medium. Radial and axial slices of this 3D image of size 140 × 140 × 70 voxels, are depicted in Fig. 6 (a) and (e). We denote this image by letter "D". This image is observed with a Zeiss LSM 780 microscope having a numerical aperture of 1.4, an oil immersed lens and a magnification of 63X. The radial and axial stepsizes are respectively 54 nm and 150 nm. In order to test our algorithm, we also consider a combination of two SI PSFs. To generate the initial PSF estimates, we use the theoretical PSF model (F. [START_REF] Gibson | Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy[END_REF]). For that, we consider the following approximate RI values of the polymer n s = 1.45, of the cover-slip RR n°8073 n c = 1.5, and of the oil immersion medium n i = 1.5. The regularizing parameters used in this test are fixed as follows α = 0.05 and β i = 10, i = 1, 2 for all tests. The restored images using SGPAM and RRLAM methods are respectively displayed in Fig. 6 (b), (f), (c) and (g). We also present in Fig. 6 (d) and (h) a restoration result using a SI PSF.

Discussions

From these experiments, we can say that: First, one can notice from table 2 that SGPAM algorithm provides faster convergence and more accurate results than RRLAM algorithm. Although one iteration of SGPAM is longer than one iteration of RRLAM, SGPAM converges faster since it requires less iterations than RRLAM to reach its convergence. Because of the static step-length value in RRLAM, the algorithm rapidly stagnates to a not necessarily optimal solution. However, in SGPAM method, thanks to the dynamic step-length selection, SGPAM gives more accurate results than RRLAM. To illus- trate this analysis, we present in Fig. 7 the global energy variation with iterations for the test image A. Second, from the presented tests, one can remark that the proposed BR gives more accurate result than BD when the blur variation is significant. This is illustrated in Fig. 5 (a) which shows the intensity profiles passing through the centers of each of the three restored beads using a SV PSF (in magenta) and a SI PSF (in black). In particular, the last bead is not well restored by the BD algorithm since the blur at the bottom of the sample (i.e. for high z) is more significant than at the top of the sample (i.e. for low z). This also can be seen on the PSF images presented in Fig. 4 (e) and (i) as well as their intensity profiles along the z axis in Fig. 5 (b) and (c). One can notice that the second true PSF (Fig. 5 (c)) is more spread than the first PSF (Fig. 5 (b)). Furthermore, using the proposed BR approach, one can efficiently recover the spread of these two different PSFs (see plots in magenta in Fig. 5 (b) and (c)). However, one can remark from the RRE and SSIM values in table 2 that restoration with SI and SV PSFs have practically the same accuracy for the test image B. This is also can be noticed by visually , (g) estimated ĥ1 with SGPAM, (h) estimated ĥ1 with RRLAM, (i) second true PSF h 2 , (j) initial PSF ĥ2 (0) , (k) estimated ĥ2 with SGPAM, (l) estimated ĥ2 with RRLAM. comparing images of Fig. 3 (c) and (d). This is can be explained by the fact that the blur variation across the processed image of Fig. 3 (b) is not very significant. Hence, when the blur variation is not very important, it is not worth applying our BR algorithm since it requires more computational time than BD without any gain in the restoration quality. Figure 5: (a) Intensity profiles along the z-axis passing through the centers of the three beads: the blue plots correspond to the original beads, the red plots correspond to the degraded beads, the magenta plots correspond to the restoration using SGPAM, the green plots correspond to the restoration using RRLAM, and the black plots correspond to the restoration using a SI PSF and SGPAM method.

RR n°8073 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)
RR n°8073 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)

Conclusion

In this report, we presented a BR method for images degraded with a SV blur and Poisson noise. In this method, the SV PSF is approximated by a convex combination of a set of SI PSFs. In a JMAP framework, we defined a convenient criterion to solve the simultaneous estimation problem of the SI PSF set and the RR n°8073 We proved the existence of a solution of the considered problem, a minimizer of the given criterion. We then proposed an SGPAM algorithm for the numerical resolution of this problem. That is, the multivariate minimization problem is solved by alternating between mono-variate minimizations using an SGP algorithm, an efficient algorithm that provides fast convergence. Comparison of the proposed method with the RRLAM algorithm (a regularized version of RL algorithm embedded with the AM procedure) on simulated and real images of digital camera and confocal microscopy shows that our algorithm is faster and more precise. Furthermore, comparison of our algorithm with a BD algorithm where one SI PSF is used reveals the potential interest of the SV blur model. Nevertheless, some limitations of this work should be reported. In fact, in our method, the number of the considered SI PSFs as well as the set of weighting functions are fixed by the user. Further investigations about the choice of these parameters and their impact on the accuracy of the solution, in a completely automatic BR framework, would be an interesting future work. Moreover, regularizing weight parameters included in our method are set in by try/error. It could be interesting to automatically estimate them during the estimation procedure.

A Gradient decomposition

We give here a decomposition of the gradient of the used regularizing terms for both 2D and 3D images. In this report, the considered regularizing terms are of the following form:

• for 2D images X ∈ R n × R n : J R (X) = ∑ 1≤i, j≤n φ D 2 i, j (53)
where D 2 i, j = (∇X) i, j 2 and the discrete gradient (∇X) i, j is computed by adding one column and one row to X as follows:

(∇X) i, j = X i+1, j -X i, j X i, j+1 -X i, j

with X n+1, j = X 1, j and X i,n+1 = X i,1 ;

• for 3D images

X ∈ R n × R n × R n : J R (X) = ∑ 1≤i, j≤n φ D 2 i, j,k (55)
RR n°8073

where D 2 i, j,k = (∇X) i, j,k 2 and the discrete gradient (∇X) i, j,k is as follows:

(∇X) i, j,k =   X i+1, j,k -X i, j,k X i, j+1,k -X i, j,k X i, j,k+1 -X i, j,k   (56) with X n+1, j,k = X 1, j,k , X i,n+1,k = X i,1,k , and X i, j,n+1 = X i, j,n+1

The function φ (t) is chosen as φ (t) = 2 √ t + ε 2 , ∀t ≥ 0, ε > 0 for the total variation term used in (36) and as φ (t) = t, ∀t ≥ 0 for the quadratic regularizing term used in (46). The parameter ε is a very small scalar used to avoid singularity points in the derivative of the total variation term (e.g. in our simulation we consider ε = 10 -8 as it is suggested in [START_REF] Vogel | Computational methods for inverse problems[END_REF], [START_REF] Zanella | Efficient gradient projection methods for edge-preserving removal of poisson noise[END_REF])). Now, let us first consider 2D images. The gradient of function ( 53) w.r.t. X i, j is expressed as: ∇ i, j J R (X) = φ D 2 i, j 2X i, j -X i+1, j -X i, j+1 +φ D 2 i-1, j X i, j -X i-1, j +φ D 2 i, j-1 X i, j -X i, j-1

with φ (t) the derivative of φ (t). A possible decomposition of the gradient ∇ i, j J R (x) into two positive terms is as follows:

U R (X) = φ D 2 i, j X i+1, j + X i, j+1 + φ D 2 i-1, j X i-1, j +φ D 2 i, j-1 X i, j-1

(58) V R (X) = 2φ D 2 i, j + φ D 2 i-1, j + φ D 2 i, j-1 X i, j (59)
Now, we give similar expressions for 3D images:

U R (X) = φ D 2 i, j,k X i+1, j,k + X i, j+1,k + X i, j,k+1 +φ D 2 i-1, j,k X i-1, j,k + φ D 2 i, j-1,k X i, j-1,k + φ D 2 i, j,k-1 X i, j,k-1

(60) V R (X) = X i, j,k [2φ D 2 i, j,k + φ D 2 i-1, j,k + φ D 2 i, j-1,k + φ D 2 i, j,k-1 ] (61)
RR n°8073

  Figure 1: (a) Circular variations of two considered weighting functions for the test image A: ψ 1 is increasing from the center to the edges of the image while ψ 2 is decreasing from the center to the edges of the image, the presented image intensities vary from 0 in blue to 1 in red, (b) considered decomposition for the test image B. The blur is assumed to be invariant on each region. Nine weighting functions are used for this test, each of them is chosen to be an indicator function (ψ i = 1 on the considered region and ψ i = 0 outside that region), (c) and (d) linear variation of the two considered weighting functions along the z-axis for the test images C and D: the presented plots show only a z-line of the 3D weighting functions.

Figure 2 :

 2 Figure 2: BR results on the test image A: (a) Original image, (b) observation, (c) restoration using SGPAM, (d) restoration using RRLAM, (e) first true PSF h 1 , (f) initial PSF ĥ1 (0), (g) estimated ĥ1 with SGPAM, (h) estimated ĥ1 with RRLAM, (i) second true PSF h 2 , (j) initial PSF ĥ2 (0) , (k) estimated ĥ2 with SGPAM, (l) estimated ĥ2 with RRLAM.

Figure 3 :

 3 Figure 3: BR results on the test image B: (a) original image, (b) degraded image, (c) restoration using SGPAM, (d) restoration using RRLAM

Figure 4 :

 4 Figure 4: BR results on the test image C: a simulated CLSM image. Only (Y, Z) sections of the 3D volumes are presented: (a) Original image, (b) observation, (c) restoration using SGPAM, (d) restoration using RRLAM, (e) first true PSF h 1 , (f) initial PSF ĥ1 (0), (g) estimated ĥ1 with SGPAM, (h) estimated ĥ1 with RRLAM,

  Figure5: (a) Intensity profiles along the z-axis passing through the centers of the three beads: the blue plots correspond to the original beads, the red plots correspond to the degraded beads, the magenta plots correspond to the restoration using SGPAM, the green plots correspond to the restoration using RRLAM, and the black plots correspond to the restoration using a SI PSF and SGPAM method. (b) and (c) intensity profiles along the z-axis of the first and second SI PSFs taken respectively at the top and the bottom of the considered sample. True PSFs are presented in blue, initial PSFs are presented in red, estimated PSFs with SGPAM are presented in magenta, and estimated PSFs with RRLAM are presented in green.

Figure 6 :

 6 Figure 6: BR results on the test image D: a bead shell image acquired by Gilbert Engler and Elie Maalouf at INRA with a Zeiss LSM 780 microscope. The first raw shows (X, Y ) slices and the second raw shows (Y, Z) slices. The first column ((a) and (e)) corresponds to the observation, the second column ((b) and (f)) corresponds to the restoration using SGPAM, the third column ((c) and (g)) corresponds to the restoration using RRLAM method, and the fourth column ((d) and (h)) corresponds to the restoration using a SI PSF and SGPAM method.

Figure 7 :

 7 Figure 7: Energy variation with iterations for the test image A.

Table 1 :

 1 SGPAM method vs. RRLAM method. Tests are presented for 2D of digital camera.

			Image A			Image B
		RRLAM	SGPAM	RRLAM	SGPAM
		SI	SV	SI	SV	SI	SV	SI	SV
	RRE (%) 8.71	6.18 7.62	6.31 1.09	1.10 1.02	1.05
	SSIM	0.933 0.955 0.943 0.955 0.974 0.974 0.975 0.975
	t m (min)	0.010	0.01 0.011	0.02 0.08	0.71 0.22	1.29
	t g (min)	10.59 18.01 6.31 18.97 8.98 72.09 5.07 41.48
	nb it	1000	1000	563	796	101	101	23	33

Table 2 :

 2 SGPAM method vs. RRLAM method. Tests are presented for 3D images of confocal microscopy.

RR n°8073

http://www.inria.fr ISSN 0249-6399

Acknowledgment

We acknowledge Gilbert Engler and Elie Maalouf for collaborating with us within the ANR Diamond project and particularly for acquiring real CLSM images at INRA. The first author acknowledge the ANR Diamond project for the financial support. We also acknowledge ATE for important interactions within the FUI Gyrovision project and for providing images for digital camera tests. We also acknowledge Roberto Cavicchioli for providing us SGP code and Praveen Pankajakshan for confocal PSF code.