Franc ¸ois Orieux

Jean-Franc ¸ois Giovannelli

Thomas Rodet

Bayesian estimation of regularization and PSF parameters for Wiener-Hunt deconvolution

This paper tackles the problem of image deconvolution with joint estimation of PSF parameters and hyperparameters. Within a Bayesian framework, the solution is inferred via a global a posteriori law for unknown parameters and object. The estimate is chosen as the posterior mean, numerically calculated by means of a Monte-Carlo Markov chain algorithm. The estimates are efficiently computed in the Fourier domain and the effectiveness of the method is shown on simulated examples. Results show precise estimates for PSF parameters and hyperparameters as well as precise image estimates including restoration of high-frequencies and spatial details, within a global and coherent approach.

Introduction

Image deconvolution has been an active research field for several decades and recent contributions can be found in papers such as [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF][START_REF] Molina | Blind deconvolution using a variational approach to parameter, image, and blur estimation[END_REF][START_REF]Blind Image Deconvolution[END_REF]. Examples of application are medical imaging, astronomy, nondestructive testing and more generally imagery problems. In these applications, degradations induced by the observation instrument limit the data resolution while the need of precise interpretation can be of major importance. For example, this is particularly critical for long-wavelength astronomy (see e.g., [START_REF] Rodet | Data inversion for over-resolved spectral imaging in astronomy[END_REF]). In addition, the development of a high quality instrumentation system must rationally be completed by an equivalent level of quality in the development of data processing methods. Moreover, even for poor performance systems, the restoration method can be used to bypass instrument limitations.

When the deconvolution problem is ill-posed a possible solution relies on regularization, i.e., introduction of information in addition to the data and the acquisition model [START_REF] Tikhonov | Solutions of Ill-Posed Problems[END_REF][START_REF] Twomey | On the numerical solution of Fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature[END_REF]. As a consequence of regularization, deconvolution methods are specific to the class of image in accordance with the introduced information. From this standpoint, the present paper is dedicated to relatively smooth images encountered for numerous applications in imagery [START_REF] Rodet | Data inversion for over-resolved spectral imaging in astronomy[END_REF][START_REF] Jalobeanu | Hyperparameter estimation for satellite image restoration by a MCMC maximum likelihood method[END_REF][START_REF] O'sullivan | Roughness penalties on finite domains[END_REF]. The second order consequence of ill-posedness and regularization is the need to balance the compromise between different sources of information.

In the Bayesian approach [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF][START_REF] Demoment | Image reconstruction and restoration: Overview of common estimation structure and problems[END_REF], information about unknowns is introduced by means of probabilistic models. Once these models are designed, the next step is to build the a posteriori law, given the measured data. The solution is then defined as a representative point of this law and the two most classical are (1) the maximizer, and (2) the mean. From a computational standpoint, the first leads to a numerical optimization problem and the latter leads to a numerical integration problem. However, the resulting estimate depends on two sets of variables in addition to the data.

1. Firstly, the estimate naturally depends on the response of the instrument at work, namely the point spread function (PSF). The literature is predominantly devoted to deconvolution in the case of known PSF. On the contrary, the present paper is devoted to the case of unknown or poorly known PSF and there are two main strategies to tackle its estimation from the available data set (without extra measurements).

(i) In most practical cases, the instrument can be modeled using physical operating description. It is thus possible to find the equation for the PSF, at least in a first approximation. This equation is usually driven by a relatively small number of parameters. It is a common case in optical imaging where a Gaussian-shaped PSF is often used [START_REF] Pankajakshani | Blind deconvolution for thin-layered confocal imaging[END_REF]. It is also the case in other fields: interferometry [START_REF] Thibaut | Strict a priori constraints for maximum likelihood blind deconvolution[END_REF], magnetic resonance force microscopy [START_REF] Dobigeon | Hierarchical bayesian sparse image reconstruction with application to MRFM[END_REF], fluorescence microscopy [START_REF] Zhang | Gaussian approximations of fluorescence microscope point-spread function models[END_REF],. . . Nevertheless, in real experiments, the parameter values are unknown or imperfectly known and need to be estimated or adjusted in addition to the image of interest: the question is namely myopic deconvolution. (ii) The second strategy forbears the use of the parametric PSF deduced from the physical analysis and the PSF then naturally appears in a non-parametric form. Practically, the non-parametric PSF is unknown or imperfectly known and needs to be estimated in addition to the image of interest: the question is referred to as blind deconvolution for example in interferometry [START_REF] Mugnier | MISTRAL: a myopic edge-preserving image restoration method, with application to astronomical adaptive-optics-corrected long-exposure images[END_REF][START_REF] Thiébaut | MiRA: an effective imaging algorithm for optical interferometry[END_REF][START_REF] Fusco | Myopic deconvolution method for adaptive optics images of stellar fields[END_REF][START_REF] Conan | Myopic deconvolution of adaptive optics images by use of object and point-spread function power spectra[END_REF].

From an inference point of view, the difficulty of both myopic and blind problems lies in the possible lack of information resulting in ambiguity between image and PSF, even in the noiseless case. In order to resolve the ambiguity, information must be added [START_REF]Blind Image Deconvolution[END_REF][START_REF] Likas | A variational approach for Bayesian blind image deconvolution[END_REF] and it is crucial to make inquiries based on any available source of information. To this end, the knowledge of the parametric PSF represents a precious means to structure the problem and possibly resolve the degeneracies. Moreover, due to instrument design process, a nominal value as well as an uncertainty are usually available for the PSF parameters.

In addition, from a practical and algorithmic standpoint, the myopic case, i.e., the case of parametric PSF, is often more difficult due to the non-linear dependence of the observation model with respect to the PSF parameters. On the contrary, the blind case, i.e., the case of non-parametric PSF, yields a simpler practical and algorithmic problem since the observation model remains linear w.r.t. the unknown elements given the object. Despite the superior technical difficulty, the present paper is devoted to the myopic format since it is expected to be more efficient than the blind format from an information standpoint. Moreover, the blind case has been extensively studied and a large amount of paper is available [START_REF] Bishop | Blind restoration of blurred photographs via AR modelling and MCMC[END_REF][START_REF] Lam | Iterative statistical approach to blind image deconvolution[END_REF][START_REF] Xu | Maximum a posteriori blind image deconvolution with Huber-Markov random-field regularization[END_REF], while the myopic case has been less investigated, though it is of major importance. 2. Secondly, the solution depends on the probability law parameters named hyperparameters (means, variances, parameters of correlation matrix,. . . ). These parameters adjust the shape of the laws and in the same time they tune the compromise between the information provided by the a priori and the information provided by the data. In real experiments, their values are unknown and need to be estimated: the question is namely unsupervised deconvolution.

For both families of parameters (PSF parameters and hyperparameters), two approaches are available. In the first one, the parameter values are empirically tuned or estimated in a preliminary step (with Maximum Likelihood [START_REF] Jalobeanu | Hyperparameter estimation for satellite image restoration by a MCMC maximum likelihood method[END_REF] or calibration [START_REF] Cannon | Blind deconvolution of spatially invariant image blurs with phase[END_REF] for example), then the values are used in a second step devoted to image restoration given the parameters. In the second one, the parameters and the object are jointly estimated [START_REF] Molina | Blind deconvolution using a variational approach to parameter, image, and blur estimation[END_REF][START_REF] Bishop | Blind restoration of blurred photographs via AR modelling and MCMC[END_REF].

For the myopic problem, Jalobeanu et al. [START_REF] Jalobeanu | Estimation of blur and noise parameters in remote sensing[END_REF] address the case of a symmetric Gaussian PSF. The width parameter and the noise variance are estimated in a preliminary step by Maximum-Likelihood. A recent paper [START_REF] Chen | An empirical identification method of Gaussian blur parameter for image deblurring[END_REF] addresses the estimation of a Gaussian blur parameter, as in our experiment, with an empirical method. They found the Gaussian blur parameter by minimizing the absolute derivatives of the restored images Laplacian.

The present paper addresses the myopic and unsupervised deconvolution problem. We propose a new method that jointly estimates the PSF parameters, the hyperparameters, and the image of interest. It is built in a coherent and global framework based on an extended a posteriori law for all the unknown variables. The posterior law is obtained via the Bayes rule, founded on a priori laws: Gaussian for image and noise, uniform for PSF parameters and gamma or Jeffreys for hyperparameters.

Regarding the image prior law, we have paid special attention to the parametrization of the covariance matrix in order to facilitate law manipulations such as integration, conditioning or hyperparameter estimation. The possible degeneracy of the a posteriori law in some limit cases is also studied.

The estimate is chosen as the mean of the posterior law and is computed using Monte-Carlo simulations. To this end, Monte-Carlo Markov chain (MCMC) algorithms [START_REF] Robert | Monte-Carlo Statistical Methods[END_REF] enable to draw samples from the posterior distribution despite its complexity and especially the non-linear dependence w.r.t. the PSF parameters.

The paper is structured in the following manner. Sec. 2 presents the notations and states the problem. The three following sections describe our methodology: firstly the Bayesian probabilistic models are detailed in Sec. 3; then a proper posterior law is established in Sec. 4; an MCMC algorithm to compute the estimate is described in Sec. 5. Numerical results are shown in Sec. 6. Finally, Sec. 7 is devoted to conclusion and perspectives.

Notations and convolution model

Consider N pixels real square images represented in lexicographic order by vector x ∈ R N , with generic elements x n . The forward model writes

y = H w x + ǫ (1)
where y ∈ R N is the vector of data, H w a convolution matrix, x the image of interest and ǫ the modelization errors or the noise. Vector w ∈ R P stands for the PSF parameters, such as width or orientation of a Gaussian PSF.

The matrix H w is block-circulant with circulant-block (BCCB) for computational efficiency of the convolution in the Fourier space. The diagonalization [START_REF] Hunt | A matrix theory proof of the discrete convolution theorem[END_REF] of H w writes Λ H = F H w F † where F is the unitary Fourier matrix and † is the transpose conjugate symbol. The convolution, in the Fourier space, is then

• y = Λ H • x + • ǫ (2)
where

• x = F x,
• y = F y and • ǫ = F ǫ are the 2D discrete Fourier transform (DFT-2D) of image, data and noise, respectively.

Since Λ H is diagonal, the convolution is computed with a term-wise product in the Fourier space. There is a strict equivalence between a description in spatial domain (Eq. ( 1)) and in Fourier domain (Eq. ( 2)). Consequently, for coherent description and computational efficiency, all the developments are equally done in the spatial space or in the Fourier space.

For notational convenience, let us introduce the component at null-frequency

•
x 0 ∈ R and the vector of component at non-null frequencies

•

x * ∈ C N -1 so that the whole set of components writes

• x = [ • x 0 , • x * ].
Let us note 1 the vector of N components equal to 1/N, so that 1 t x is the empirical mean level of the image. The Fourier components are the • 1 n and we have:

• 1 0 = 1 and • 1 n = 0 for n = 0.
Moreover, Λ 1 = F 11 t F † is a diagonal matrix with only one non-null coefficient at null frequency.

Bayesian probabilistic model

This section presents the prior law for each set of parameters. Regarding the image of interest, in order to account for smoothness, the law introduces high-frequency penalization through a differ-ential operator on the pixel. A conjugate law is proposed for the hyperparameters and a uniform law is considered for the PSF parameters.

Moreover, we have paid a special attention to the image prior law parametrization. In the next section we present several parametrization in order to facilitate law manipulations such as integration, conditioning or hyperparameter estimation. Moreover, the correlation matrix of the image law may become singular in some limit cases resulting in a degenerated prior law (when p(x) = 0 for all x ∈ R N ). Based on this parametrization, Sec. 4 studies the degeneracy of the posterior in relation with the parameters of the prior law.

3.A. Image prior law

The probability law for the image is a Gaussian field with a given precision matrix P parametrized by a vector γ. The pdf reads

p(x|γ) = (2π) -N/2 det[P ] 1/2 exp - 1 2 x t P x . (3)
For computational efficiency, the precision matrix is designed (or approximated) in a toroidal manner, and it is diagonal in the Fourier domain Λ P = F P F † . Thus, the law for x also writes

p(x|γ) = (2π) -N/2 det[F ] det[Λ P ] 1/2 det[F † ] exp - 1 2 x t F † Λ P F x (4)
= (2π) -N/2 det[Λ P ] 1/2 exp - 1 2 • x † Λ P • x ( 5
)
and it is sometimes referred to [START_REF] Calder | Introduction to Whittle (1953) 'The analysis of multiple stationary time series[END_REF] as a Whittle approximation (see also [28, p.133]) for the Gaussian law. The filter obtained for fixed hyperparameters is also the Wiener-Hunt filter [START_REF] Hunt | Deconvolution of linear systems by constrained regression and its relationship to the Wiener theory[END_REF], as described in Sec. 5.A. This paper focuses on smooth images, thus on positive correlation between pixels. It is introduced by high-frequencies penalty using any circulant differential operator: p-th differences between pixels, Laplacian, Sobel. . . The differential operator is denoted by D and its diagonalized form by Λ D = F DF † . Then, the precision matrix writes P = γ 1 D t D and its Fourier counterpart writes

Λ P = γ 1 Λ † D Λ D = diag 0, γ 1 | • d 1 | 2 , . . . , γ 1 | • d N -1 | 2 (6)
where γ 1 is a positive scale factor, diag builds a diagonal matrix from elementary components and

• d n is the n-th DFT-2D coefficient of D.
Under this parametrization of P , the first eigenvalue is equal to zero corresponding to the absence of penalty for the null frequency • x 0 , i.e., no information accounted for about the empirical mean level of the image. As a consequence, the determinant vanishes det[P ] = 0 resulting in a degenerated prior. To manage this difficulty, several approaches have been proposed. Some authors [START_REF] Molina | Blind deconvolution using a variational approach to parameter, image, and blur estimation[END_REF][START_REF] Mardia | Multivariate Analysis[END_REF] still use this prior despite its degeneracy and this approach can be analyzed in two ways.

1. On the one hand, it can be seen as a non-degenerated law for • x * , the set of non-null frequency components only. In this format, the prior does not affect any probability to the null frequency component and the Bayes rule does not apply to this component. Thus, this strategy yields an incomplete posterior law, since the null frequency is not embedded in the methodology. 2. On the other hand, it can be seen as a degenerated prior for the whole set of frequencies. The application of the Bayes rule is then somewhat confusing due to degeneracy. In this format, the posterior law cannot be guaranteed to remain non-degenerated.

Anyway, none of the two standpoints yields a posterior law that is both non-degenerated and addressing the whole set of frequencies.

An alternative parametrization relies on the energy of x. An extra term γ 0 I, tuned by γ 0 > 0, in the precision matrix [START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF], introduces information for all the frequencies including • x 0 . The precision matrix writes

Λ P = γ 0 I + γ 1 Λ † D Λ D = diag γ 0 , γ 0 + γ 1 | • d 1 | 2 , . . . , γ 0 + γ 1 | • d N -1 | 2 (7)
with a determinant

det[Λ P ] = N -1 n=0 γ 0 + γ 1 | • d n | 2 . (8)
The obtained Gaussian prior is not degenerated and undoubtedly leads to a proper posterior. Nevertheless, the determinant Eq. ( 8) is not separable in γ 0 and γ 1 . Consequently, the conditional posterior for these parameters is not a classical law and future development will be more difficult. Moreover, the non-null frequencies

•

x * are controlled by two parameters γ 0 and γ 1 p(

• x|γ 0 , γ 1 ) = p( • x 0 |γ 0 )p( • x * |γ 0 , γ 1 ). (9)
The proposed approach to manage the degeneracy relies on the addition of a term for the null frequency only Λ 1 = diag (1, 0, . . . , 0)

Λ P = γ 0 Λ † 1 Λ 1 + γ 1 Λ † D Λ D . ( 10
) = diag γ 0 , γ 1 | • d 1 | 2 , . . . , γ 1 | • d N -1 | 2 .
The determinant has a separable expression

det[Λ P ] = γ 0 γ N -1 1 N -1 n=1 | • d n | 2 , (11)
i.e., the precision parameters have been factorized. In addition, each parameter controls a different set of frequencies: p(

• x|γ 0 , γ 1 ) = p( • x 0 |γ 0 )p( • x * |γ 1 ) ,
γ 0 drives the empirical mean level of the image

•

x 0 and γ 1 drives the smoothness

•

x * of the image. With the Fourier precision structure of Eq. ( 10), we have the non-degenerated prior law for the image that addresses separately all the frequencies with a factorized partition function w.r.t. (γ 0 , γ 1 )

p(x|γ 0 , γ 1 ) = (2π) -N/2 N -1 n=1 | • d n | γ 1/2 0 γ (N -1)/2 1 exp - γ 0 2 • x 0 2 - γ 1 2 Λ D * • x * 2 . ( 12
)
where Λ D * is obtained from Λ D without the first line and column. The next step is to write the a priori law for the noise in an explicit form and the other parameters, including the law parameters γ and the instrument parameters w.

3.B. Noise and data laws

From a methodological standpoint, any statistic can be included for errors (measurement and model errors). It is possible to account for correlations in the error process or to account for a non-Gaussian law, e.g., Laplacian law, generalized Gaussian law, or other laws based on robust norm,. . . In the present paper, the noise is modeled as zero-mean white Gaussian vector with unknown precision parameter

γ ǫ p(ǫ|γ ǫ ) = (2π) -N/2 γ N/2 ǫ exp - γ ǫ 2 ǫ 2 . (13)
Consequently, the likelihood for the parameters given the observed data writes

p(y|x, γ ǫ , w) = (2π) -N/2 γ N/2 ǫ exp - γ ǫ 2 y -H w x 2 . ( 14
)
It naturally depends on the image x, on the noise parameter γ ǫ and on the PSF parameters w embedded in H w . It clearly involves a least squares discrepancy that can be rewritten in the Fourier domain:

y -H w x 2 = • y -Λ H • x 2 .

3.C. Hyperparameters law

A classical choice for hyperparameter law relies on conjugate prior [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF]: the conditional posterior for the hyperparameters is in the same family as its prior. It results in practical and algorithmic facilities: update of the laws amounts to update of a small number of parameters.

The three parameters γ 0 , γ 1 and γ ǫ are precision parameters of Gaussian laws Eq. ( 12) and ( 14) and a conjugate law for these parameters is the Gamma law (see Appendix B). Given parameters (α i , β i ), for i = 0, 1 or ǫ, the pdf reads

p(γ i ) = 1 β α i i Γ(α i ) γ α i -1 i exp (-γ i /β i ) , ∀γ i ∈ [0, +∞[ (15)
In addition to computational efficiency, the law allows for non-informative priors. With specific parameter values, one obtains two improper non-informative prior : the Jeffreys' law p(γ) = 1/γ and the uniform law p(γ) = U [0,+∞[ (γ) with (α i , β i ) set to (0, +∞) and (1, +∞), respectively. Jeffreys' law is a classical law for the precisions and is considered as non-informative [START_REF] Kass | The selection of prior distributions by formal rules[END_REF]. This law is also invariant to power transformations: the law of γ n [START_REF] Kass | The selection of prior distributions by formal rules[END_REF][START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF] is also a Jeffreys' law. For these reasons development is done using the Jeffreys' law.

3.D. PSF parameters law

Regarding the PSF parameters w, we consider that the instrument design process or a physical study provides a nominal value w with uncertainty δ, that is to say w ∈ [wδ , w + δ]. The "Principle of Insufficient Reason" [START_REF] Kass | The selection of prior distributions by formal rules[END_REF] leads to a uniform prior on this interval

p(w) = U w,δ (w) (16)
where U w,δ is a uniform pdf on [wδ , w + δ]. Nevertheless, within the proposed framework, the choice is not limited and other laws, such as Gaussian, are possible. Anyway other choices do not allow easier computation because of the non-linear dependency of the observation model w.r.t. PSF parameters.

Proper posterior law

At this point, the prior law of each parameter is available: the PSF parameters, the hyperparameters and the image. Thus, the joint law for all the parameters is built by multiplying the likelihood Eq. ( 14) and the a priori laws Eq. ( 12), ( 15) and ( 16)

p( • x, γ ǫ , γ 0 , γ 1 , w, • y) = p( • y| • x, γ ǫ , w)p( • x|γ 0 , γ 1 )p(γ ǫ )p(γ 0 )p(γ 1 )p(w) (17)
and explicitly

p( • x, γ ǫ , γ 0 , γ 1 , w, • y) = (2π) -N N -1 n=1 | • d n | β αǫ ǫ Γ(α ǫ ) β α 0 0 Γ(α 0 ) β α 1 1 Γ(α 1 ) γ αǫ+N/2-1 ǫ γ α 0 -1/2 0 γ α 1 +(N -1)/2-1 1 exp - γ ǫ β ǫ - γ 0 β 0 - γ 1 β 1 U w,δ (w) exp - γ ǫ 2 • y -Λ H • x 2 - γ 0 2 • x 0 2 - γ 1 2 Λ D • x 2 . ( 18
)
According to the Bayes rule, the a posteriori law reads

p( • x, γ ǫ , γ 0 , γ 1 , w| • y) = p( • x, γ ǫ , γ 0 , γ 1 , w, • y) p( • y) (19)
where p(

• y) is a normalization constant p( • y) = p( • y, • x, γ, w) d • x dγ dw. (20)
As described before, setting γ 0 = 0 leads to degenerated prior and joint laws. However, when the observation system preserves the null frequency γ 0 can be considered as a nuisance parameter. In addition, only prior information on the smoothness is available.

In Bayesian framework, a solution to eliminate the nuisance parameters is to integrate them out in the a posteriori law. According to our parametrization Sec. 3.A, the integration of γ 0 is the integration of a Gamma law. Application of Appendix B.B on γ 0 in the a posteriori law Eq. ( 19) provides

p( • x, γ ǫ , γ 1 , w| • y) = p( • x 0 )p( • y, • x * , γ ǫ , γ 1 , w| • x 0 ) p( • x 0 )p( • y, • x * , γ ǫ , γ 1 , w| • x 0 ) dγ ǫ dγ 1 dw d • x * d • x 0 (21)
with p(

• x 0 ) = p( • x 0 |γ 0 )p(γ 0 ) dγ 0 =   1 + β 0 • x 2 0 2   -α 0 -1/2 . ( 22
)
Now the parameter is integrated, the parameters α 0 and β 0 are set to remove the null frequency penalization. Since we have α 0 > 0 and β 0 > 0 we get

(1 + β 0 • x 2 0 /2) -α 0 -1/2 ≤ 1 and the joint law is majored   1 + β 0 • x 2 0 2   -α 0 -1/2 p( • y, • x * , γ ǫ , γ 1 , w| • x 0 ) ≤ p( • y, • x * , γ ǫ , γ 1 , w| • x 0 ). (23)
Consequently, by the dominated convergence theorem [START_REF] Lang | Real and functional analysis[END_REF], the limit of the law with α 0 → 1 and β 0 → 0 can be placed under the integral sign at the denominator. Then the null-frequency penalization p( • x 0 ) from the numerator and denominator are removed. It is equivalent with the integration of the γ 0 parameter under a Dirac (see appendix B). The equation is simplified and the integration with respect to

• x 0 in the denominator Eq. (20) R p( • y| • x, γ ǫ , w)p( • x * |γ 1 )p(γ 1 , γ ǫ , w) d • x 0 ∝ R p( • y 0 | • x 0 , γ ǫ , w) d • x 0 (24) ∝ R exp - γ ǫ 2 • y 0 - • h 0 • x 0 2 d • x 0 (25)
converges if and only if

• h 0 = 0: the null frequency is observed. If this condition is met, Eq. ( 21)

with β 0 = 0 and α 0 = 1 is a proper posterior law for the image, the precision parameters and the PSF parameters. In other words, if the average is observed, the degeneracy of the a priori law is not transmitted to the a posteriori law.

Then, the obtained a posteriori law writes

p( • x, γ ǫ , γ 1 , w| • y) = p( • x, γ ǫ , γ 1 , w, • y) p( • y) ∝ γ αǫ+N/2-1 ǫ γ α 1 +(N -1)/2-1 1 U w,δ (w) exp - γ ǫ 2 • y -Λ H • x 2 - γ 1 2 Λ D * • x * 2 exp - γ ǫ β ǫ - γ 1 β 1 . (26)
Finally, inference is done on this law Eq. ( 26). If the null frequency is not observed, or information must be added, the previous Eq. ( 19) can be used.

Posterior mean estimator and law exploration

This section presents the algorithm to explore the posterior law Eq. ( 19) or [START_REF] Hunt | A matrix theory proof of the discrete convolution theorem[END_REF] and to compute an estimate of the parameters. For this purpose, Monte Carlo Markov chain is used to provide samples. Firstly, the obtained samples are used to compute different moments of the law. Afterwards, they are also used to approximate marginal laws as histograms. These two representations are helpful to analyse the a posteriori law, the structure of the available information and the uncertainty. They are used in Sec. 6.C.2 to illustrate the mark of the ambiguity in the myopic problem.

Here, the samples of the a posteriori law are obtained by a Gibbs sampler [START_REF] Robert | Monte-Carlo Statistical Methods[END_REF][START_REF] Brémaud | Markov Chains. Gibbs fields, Monte Carlo Simulation, and Queues[END_REF][START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF]: it consists in iteratively sampling the conditional posterior law for a set of parameters given the other parameters (obtained at previous iteration). Typically, the sampled laws are the law of • x, γ i and w. After a burn-in time, the complete set of samples are under the joint a posteriori law. The three next sections present each sampling step.

5.A. Sampling the image

The conditional posterior law of the image is a Gaussian law

• x (k+1) ∼ p • x| • y, γ (k) ǫ , γ (k) 0 , γ (k) 1 , w (k) (27)
∼ N µ (k+1) , Σ (k+1) . ( 28
)
The covariance matrix is diagonal and writes

Σ (k+1) = γ (k) ǫ |Λ (k) H | 2 + γ (k) 0 |Λ 1 | 2 + γ (k) 1 |Λ D | 2 -1 (29)
and the mean

µ (k+1) = γ (k) ǫ Σ (k+1) Λ † H (k) • y. ( 30
)
where † is the transpose conjugate symbol. The vector µ (k+1) is the regularized least square solution at the current iteration (or the Wiener-Hunt filter). Clearly, if the null-frequency is not observed

• h 0 = 0 and if γ 0 = 0, the covariance matrix Σ is not invertible and the estimate is not defined as described Sec. 4.

Finally, since the matrix is diagonal, the sample

• x (k+1
) is obtained by a term-wise product of F ǫ (where ǫ is white Gaussian) with the standard deviation matrix Σ (k+1) 1/2 followed by the addition of the mean µ (k+1) also computed with term-wise products Eq. [START_REF] Mardia | Multivariate Analysis[END_REF]. Consequently, the sampling of the image is effective even with high-dimensional object.

5.B. Sampling precision parameters

The conditional posterior laws of the precisions are Gamma corresponding to their prior law with parameters updated by the likelihood

γ (k+1) i ∼ p γ i | • y, • x (k+1) , w (k) (31)
∼ G γ i |α (k+1) i , β (k+1) i
.

For γ ǫ , γ 0 and γ 1 the parameters law are, respectively,

α (k+1) ǫ = α ǫ + N/2 and β (k+1) ǫ = β -1 ǫ + 1 2 • y -Λ (k) H • x (k+1) 2 -1 , (33)
α (k+1) 0 = α 0 + 1/2 and β (k+1) 0 = β -1 0 + 1 2 • x (k+1) 0 2 -1 , (34) α (k+1) 1
= α 1 + (N -1)/2 and β

(k+1) 1 = β -1 1 + 1 2 Λ D • x (k+1) 2 -1 . (35)
In the case of Jeffreys' prior, the parameters are

α (k+1) ǫ = N/2 and β (k+1) ǫ = 2/ • y -Λ (k) H • x (k+1) 2 , (36)
α (k+1) 0 = 1/2 and β (k+1) 0 = 2/ • x (k+1) 0 2 , (37) α (k+1) 1
= (N -1)/2 and β

(k+1) 1 = 2/ Λ D • x (k+1) 2 . ( 38
)
Remark 1 -If the a posteriori law Eq. ( 26) without γ 0 is considered, there is no need to sample this parameter (Eq. ( 34) and [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF] are not useful) and γ (k) 0 = 0 in Eq. ( 29).

5.C. Sample PSF parameters

The conditional law for PSF parameters writes

w (k+1) ∼ p w| • y, • x (k+1) , γ (k+1) ǫ ( 39
) ∝ exp - γ (k+1) ǫ 2 • y -Λ H,w • x (k+1) 2 (40)
where parameters w are embedded in the PSF Λ H . This law is not standard and intricate: no algorithm exists for direct sampling and we use the Metropolis-Hastings (M.-H.) method to bypass this difficulty. In M.-H. algorithm, a sample w p is proposed and accepted with a certain probability. This probability depends on the ratio between the likelihood of the proposed value and the likelihood of the current value w (k) . In practice, in the independent form described in appendix C, with prior law as proposition law, it is divided in several steps.

PROPOSITION: Sample a proposition

w p ∼ p(w) = U [a b] (w). (41)
2. PROBABILITY OF ACCEPTATION: Calculate the criterion

J w (k) , w p = γ (k+1) ǫ 2 • y -Λ H,w (k) • x (k+1) 2 - • y -Λ H,wp • x (k+1) 2 . ( 42
)
3. UPDATE: Sample t ∼ U [0 1] and takes

w (k+1) =    w p if log t < J w (k) otherwise. ( 43
)

5.D. Empirical mean

The sampling of

•

x, γ and w are repeated iteratively until the law has been sufficiently explored.

These samples

• x (k) , γ (k) , w (k)
follow the global a posteriori law of Eq. [START_REF] Bishop | Blind restoration of blurred photographs via AR modelling and MCMC[END_REF]. By the large numbers law, the estimate, defined as the posterior mean, is approximated by

x = F † E[ • x] ≈ F † 1 K K-1 k=0 • x (k) . (44)
As described by Eq. ( 44), to obtain an estimate of the image in the spatial space, all the computation are achieved recursively in the Fourier space with a single IFFT at the end. An implementation example in pseudo code is described Fig. 9.

Deconvolution results

This section presents numerical results obtained by the proposed method. In order to completely evaluate the method, true value of all parameters x, w, γ ǫ but also γ 1 , γ 0 is needed. In order to achieve this, an entirely simulated case is studied: image and noise are simulated under their respective prior laws Eq. ( 12) and ( 13) with given values of γ 0 , γ 1 and γ ǫ . Thanks to this protocol, all experimental conditions are controlled and the estimation method is entirely evaluated. The method has also been applied in different conditions (lower signal to noise ratio, broader PSF, different and realistic (non-simulated) images, . . . ) and showed similar behaviour. However, in the case of realistic images, since the true value of the hyperparameters γ 0 and γ 1 is unknown, the evaluation cannot be complete.

6.A. Practical experimental conditions

Concretely, a 128 × 128 image is generated in the Fourier space as the product of a complex white Gaussian noise and the a priori standard deviation matrix Σ = (γ 0 Λ

† 1 Λ 1 + γ 1 Λ † D Λ D ) -1/2
, given by Eq. [START_REF] Pankajakshani | Blind deconvolution for thin-layered confocal imaging[END_REF]. The chosen matrix Λ D results from the FFT-2D of the Laplacian operator [0 1 0; 1 -4 1; 0 1 0] /8 and the parameter values are γ 0 = 1 and γ 1 = 2.

These parameters provide the image shown in Fig. 1(a) : it is an image with smooth features similar to a cloud. Pixels have numerical values between -100 and 150, and the profile line 68 shows fluctuations around a value of -40.

The a priori law for the hyperparameters are set to the non-informative Jeffreys' law by fixing the (α i , β i ) to (0, +∞), as explained in Sec. 3.C. In addition, the PSF is obtained in the Fourier space by discretization of a normalized Gaussian shape

• h(ν α , ν β ) = exp -2π 2 ν 2 α (w α cos 2 ϕ + w β sin 2 ϕ) + ν 2 β (w α sin 2 ϕ + w β cos 2 ϕ) + 2ν α ν β sin ϕ cos ϕ (w α -w β ) (45)
with frequencies (ν α , ν β ) ∈ [-0.5; 0.5] 2 . This low-pass filter, illustrated in Fig. 2, is controlled by three parameters:

• two width parameters w α and w β set to 20 and 7, respectively. Their a priori laws are uniform: p(w α ) = U The empirical mean level of the image is correctly observed (the null frequency coefficient of H w is • h 0 = 1) so the parameter γ 0 is considered as a nuisance parameter. Consequently it is integrated out under a Dirac (see Sec. 4). This is equivalent to fix its value to 0 in the algorithm Fig. 9, line 4.

Finally, the method is evaluated on two different situations.

1. The unsupervised and non-myopic case: the parameters w are known. Consequently, there is no Metropolis-Hastings step (Sec. 5.C): lines 9 to 16 are ignored in the algorithm of Fig. 9 and w is set to its true value. To obtain sufficient law exploration, the algorithm is run until the difference between two successive empirical means is less than 10 -3 . In this case, 921 samples are necessary and they are computed in approximately 12 seconds on a processor at 2.66 GHz with Matlab, 2. The unsupervised and myopic case: all the parameters are estimated. To obtain sufficient law exploration, the algorithm is run until the difference between two successive empirical means is less than 5 × 10 -5 . In this case, 18 715 samples are needed and they are computed in approximately 7 minutes.

Remark 2 -The algorithm has also been run for up to 1 000 000 samples, in both cases, without perceptible qualitative changes.

6.B. Estimation results

6.B.1. Images

The two results for the image are given Figs. 1(c) and 1(d) for the non-myopic and the myopic cases, respectively. The effect of deconvolution is notable on the image, as well as on the shown profile. The object is correctly positioned, the orders of magnitude are respected and the mean level is correctly reconstructed. The image is restored, more details are visible and the profiles are closer matching to the true image than data. More precisely, the pixels 20-25 of the 68-th line in Fig. 1 show the restoration of the original dynamic whereas it is not visible in the data. Between pixels 70 and 110, fluctuations not visible in data are also correctly restored.

In order to visualize and study the spectral contents of the images, circular average of empirical power spectral density is considered and called "spectrum" hereafter. The subjacent spectral variable is a radial frequency f such as f 2 = ν 2 α +ν 2 β . The spectrum of the true object, data and restored object are shown Figs. 3(a) and 3(b) in non-myopic and myopic cases, respectively. It is clear that the spectrum of the true image is correctly retrieved, in both cases, up to the radial frequency f ≈ 0.075. Above this frequency, noise is clearly dominant and information about the image is almost lost. In other words, the method produces correct spectral equalization in the properly observed frequency band. The result is expected from a Wiener-Hunt method but the achievement is the joint estimation of hyperparameter and instrument parameters in addition to the correct spectral equalization.

Concerning a comparison between non-myopic and myopic cases, there is no visual differences. The spectrum Figs. 3(a) and 3(b) in non-myopic and myopic cases respectively are visually indistinguishable. This is also the case when comparing Figs. 1(c) and 1(d) and especially 68-th line. From a more precise quantitative evaluation, a slight difference is observed and detailed below.

In order to quantify performances, a normalized euclidean distance

e = x -x * / x * (46)
between an image x and the true image x * is considered. It is computed between true image and estimate images as well as between true image and data. Results are reported in Tab. 1 and confirm that the deconvolution is effective with an error of approximately 6 % in myopic case compared to 11 % with data. Both non-myopic and myopic deconvolution reduce error by a factor 1.7 with respect to the observed data.

Regarding a comparison between non-myopic and myopic case, the errors are almost the same, with a slightly lower value for the non-myopic case, as expected. This difference is coherent with the intuition: more information are injected in the non-myopic case through the true PSF parameters values.

6.B.2. Hyperparameters and instrument parameters

Concerning the other parameters, their estimates are close to the true values and are reported in Tab. 2. The γ ǫ estimate is very close to the true value with γ ǫ = 0.49 instead of 0.5 in the two cases. The error for the PSF parameters are 0.35%, 2.7% and 1.9% for w α , w β and ϕ, respectively. The value of γ 1 is underestimated in the two cases with approximately 1.7 instead of 2. All the true values fall in the μ ± 3σ interval.

In order to deepen the numerical study, the paper evaluates the capability of the method to accurately select the best values for hyperparameters and instrument parameters. To this end, we compute the estimation error Eq. ( 46) for a set of "exhaustive" values of the parameters [γ ǫ , γ 1 , w α , w β , ϕ]. The protocol is the following: 1) choose a new value for a parameter (γ ǫ for example) and fix the other parameters to the value provided by our algorithm, 2) compute the Wiener-Hunt solution (Sec. 5.A) and 3) compute the error index.

Results are reported in Fig. 4. In each case, smooth variation of error is observed when varying hyperparameters and instrument parameters and an unique optimum is visible. By this way, one can find the value of the parameters that provide the best Wiener-Hunt solution when the true image x ⋆ is known. It is reported on Tab. 1 and shows almost imperceptible improvement: optimization of the parameters (based on the true image x ⋆ ) allow negligible improvement (smaller than 0.02 % as reported in Tab. 1).

So, the main conclusion is that, the unsupervised and myopic proposed approach is a relevant tool in order to tune parameters: it works (without the knowledge of the true image), as well as an optimal approach (based on the knowledge of the true image).

6.C. A posteriori law characteristics

This section describes the a posteriori law using histograms, means and variances of the parameters. The sample histograms, Figs. 5 and6, provide an approximation of the marginal posterior law for each parameter. Tabs. 1 and 2 report the variance for the image and law parameters respectively and thus allow to quantify the uncertainty.

6.C.1. Hyperparameter characteristics

The histograms for γ ǫ and γ 1 , Fig. 5, are concentrated around a mean value in both non-myopic and myopic cases. The variance for γ ǫ is lower than the one for γ 1 and it can be explained as follows.

The observed data are directly impacted by noise (present at the system output) whereas they are indirectly impacted by the object (present at the system input). The convolution system damages the object and not the noise: as a consequence, the parameter γ ǫ (that drives noise law) is more reliably estimated than γ 1 (that drives object law).

A second observation is the smaller variance for γ 1 in the non-myopic case Fig. 5(c) than in the myopic case Fig. 5(d). It is the consequence of the addition of information in the non-myopic case w.r.t. the myopic one, through the value of the PSF parameters. In the myopic case, the estimates are founded on the knowledge of an interval for the values of the instrument parameters, whereas in the non-myopic case, the estimates are founded on the true values for the instrument parameters.

6.C.2. PSF parameter characteristics

Fig. 6 gives histograms for the three PSF parameters and their appearances are quite different from the one for hyperparameters. The histograms for w α and w β , Figs. 6(a) and 6(b) are not as concentrated as the one of Fig. 5 for hyperparameters. Their variances are quite large with regards to the interval of the prior law. On the contrary, the histogram for the parameter ϕ, Fig. 6(c), has the smallest variance. It is analyzed as a consequence of a larger sensitivity of the data w.r.t. the parameter ϕ than w.r.t. the parameters w α and w β . In an equivalent manner, the observed data are more informative about the parameter ϕ than about the parameters w α and w β .

6.C.3. Mark of the myopic ambiguity

Finally, a correlation between parameters (γ 1 , w α ) and (γ 1 , w β ) is visible on their joint histograms Fig. 7. It can be interpreted as a consequence of the ambiguity in the primitive myopic deconvolution problem, in the following manner: the parameters γ 1 and w both participate in the interpretation of the spectral content of data, γ 1 as a scale factor and w as a shape factor. An increase of w α or w β results in a decrease of the cutoff frequency of the observation system. In order to explain the spectral content of a given data set, the spectrum of the original image must contain more high frequencies, i.e., a smaller γ 1 . This is also observed on the histogram illustrated Fig. 7(a).

6.D. MCMC algorithm characteristics

Globally, the chains of Figs. 5 and6, have a Markov feature (correlated) and explore the parameter space. They have a burn-in period followed by a stationary state. This characteristic has always been observed regardless the initialization. For fixed experimental conditions, the stationary state of multiple runs was always around the same value. Considering different initializations, the only visible change is on the length of the burn-in period.

More precisely, the chain of γ ǫ is concentrated in a small interval, the burn-in period is very short (less than 10 samples) and its evolution seems independent of the other parameters. The chain of γ 1 has a larger exploration, the burn-in period is longer (approximately 200 samples) and the histogram is larger. This is in accordance with the analysis of Section 6.C.1.

About the PSF parameters, the behaviour is different for (w α , w β ) and ϕ. The chain of the two width parameters has a very good exploration with quasi-instantaneous burn-in period. Conversely, the chain of ϕ is more concentrated and its burn-in period is approximately 4 000 samples. This is also in accordance with previous analysis (Section 6.C.2).

Acceptation rates in the Metropolis-Hastings algorithm are reported in Tab. 3: they are quite small, especially for the rotation parameter. This is due to the structure of the implemented algorithm: an independant Metropolis-Hastings algorithm with the prior law as a proposition law. The main advantage of this choice is its simplicity but as a counterpart, a high rejection rate is observed due to a large a priori interval for the angle parameter. A future work will be devoted to the design of more accurate proposition law.

6.E. Robustness of prior image model

Fig. 8 illustrates the proposed method on a more realistic image with heterogeneous spatial structures. The original is the Lena image and the data has been obtained with the same Gaussian PSF and also corruption by white Gaussian noise. The Fig. 8(b) shows that the restored image is closer to the true one than the data. Smaller structures are visible and edges are sharper, for example around pixel 200. The estimated parameters are γ ǫ = 1.98 while the true value is γ ⋆ ǫ = 2. Concerning the PSF parameters, the results are w α = 19.3, w β = 7.5 and ϕ = 1.15 while the true values are respectively w ⋆ α = 20, w ⋆ β = 7 and ϕ ⋆ = 1.05 as in the previous section. Here again, the estimated PSF parameters are close to the true values giving a first assessment of the capability of the method in a more realistic context.

Conclusion and perspectives

This paper presents a new global and coherent method for myopic and unsupervised deconvolution of relatively smooth images. It is built within a Bayesian framework and a proper extended a posteriori law for the PSF parameters, the hyperparameters and the image. The estimate, defined as the posterior mean, is computed by means of an MCMC algorithm in less than a few minutes.

Numerical assessment testifies that the parameters of the PSF and the parameters of the prior laws are precisely estimated. In addition, results also demonstrate that the myopic and unsupervised deconvolved image is closer to the true image than the data and show true restored high-frequencies as well as spatial details.

The paper focuses on linear invariant model often encountered in astronomy, medical imaging, nondestructive testing and especially in optical problems. Non-invariant linear models can also be considered in order to address other applications such as spectrometry [START_REF] Rodet | Data inversion for over-resolved spectral imaging in astronomy[END_REF] or fluorescence microscopy [START_REF] Zhang | Gaussian approximations of fluorescence microscope point-spread function models[END_REF]. The loss of invariance property precludes entirely Fourier-based computations but the methodology remains valid and practicable. In particular, it is possible to draw samples of the image by means of an optimization algorithm [START_REF] Orieux | Stochastic sampling of large dimension nonstationnary gaussian field for image restoration[END_REF].

Gaussian law, related to L 2 penalization, is known for possible excessive sharp edges penalization in the restored object. The use of convex L 2 -L 1 penalization [START_REF] Künsch | Robust priors for smoothing and image restoration[END_REF][START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF][START_REF] Giovannelli | Unsupervised Bayesian convex deconvolution based on a field with an explicit partition function[END_REF] or non convex L 2 -L 0 penalization [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] can overcome this limitation. In these cases a difficulty occurs in the development of myopic and unsupervised deconvolution: the partition function of the prior law for the image is in intricate or even unknown dependency w.r.t. the parameters [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF][START_REF] Jalobeanu | Hyperparameter estimation for satellite image restoration by a MCMC maximum likelihood method[END_REF][START_REF] Descombes | Estimation of Markov random field prior parameters using Markov chain Monte Carlo maximum likelihood[END_REF]. However a recent paper [START_REF] Giovannelli | Unsupervised Bayesian convex deconvolution based on a field with an explicit partition function[END_REF] overcome the difficulty resulting in an efficient unsupervised deconvolution and we plan to extend this work for the myopic case.

Regarding noise, Gaussian likelihood limits robustness to outliers or aberrant data and it is possible to appeal to robust law such as Huber penalization in order to bypass the limitation. Nevertheless, the partition function for the noise law is again difficult or impossible to manage and it is possible to resort to the idea proposed in [START_REF] Giovannelli | Unsupervised Bayesian convex deconvolution based on a field with an explicit partition function[END_REF] to overcome the difficulty.

Finally, estimation of parameters of correlation matrix (cutoff frequency, attenuation coefficients,. . . ) is possible within the same methodological framework. This could be achieved for the correlation matrix of the object or the noise. As for the PSF parameters, the approach could rely on an extended a posteriori law, including the new parameters and a Metropolis-Hastings sampler.

Acknowledgment

The authors would like to thank Professor Alain Abergel in IAS laboratory at Université Paris-Sud 11, France, for fruitful discussions and constructive suggestions. The authors are also grateful to Cornelia Vacar, IMS laboratory, for carefully reading the paper.

A. Law in Fourier space

For a Gaussian vector x ∼ N (µ, Σ), the law for • x = F x (the FFT of x) is also Gaussian whose first two moments are the following:

• The mean is

• µ = E[ • x] = F E[ • x] = F µ. (47)
• The covariance matrix is

• Σ = E[( • x - • µ)( • x - • µ) † ] = F ΣF † . (48)
Moreover, if the covariance matrix Σ is circulant it writes

• Σ = F ΣF † = Λ Σ . (49)
i.e., the covariance matrix

• Σ is diagonal.

B. The Gamma probability density

B.A. Definition

The Gamma pdf for γ > 0, with given parameter α > 0 and β > 0, is written

G(γ|α, β) = 1 β α Γ(α) γ α-1 exp (-γ/β) . (50)
Tab. 4 gives three limit cases for (α, β). The following properties hold:

• The mean is E G [γ] = αβ • The variance is V G [γ] = αβ 2
• The maximiser is β(α -1) if and only if α > 1

B.B. Marginalisation

First consider a N dimensional zero-mean Gaussian vector with a given precision matrix γΓ with γ > 0. The pdf reads

p(x|γ) = (2π) -N/2 γ N/2 det[Γ] 1/2 exp -γx t Γx /2 . (51)
So consider the conjugate pdf for γ as a Gamma law with parameter (α, β) (see previous Annex). The joint law for (x, γ) is the product of the pdf given by Eq. ( 50) and Eq. (51): p(x, γ) = p(x|γ)p(γ). The marginalization of the joint law is known [START_REF] Box | Bayesian inference in statistical analysis[END_REF]:

p(x) = R + p(x|γ)p(γ) dγ = β N/2 det[Γ] 1/2 Γ (α + N/2) (2π) N/2 Γ(α) 1 + βx t Γx 2 -α-N/2 (52)
which is a N dimensional t-Student law of 2α degrees of freedom with a βΓ precision matrix. Finally, the conditional law reads:

p(γ|x) = (2π) -N/2 det[Γ] 1/2 β α Γ(α) γ α+N/2-1 exp -γ x t Γx /2 + 1/β . (53)
Thanks to conjugacy, it is also a Gamma pdf with parameters ᾱ , β given by ᾱ = α + N/2 and β-1 = β -1 + 2/(x t Γx).

C. The Metropolis-Hastings algorithm

The Metropolis-Hastings algorithm provides samples of a target law f (w) that cannot be directly sampled but can be evaluated, at least up to a multiplicative constant. Using the so called "instrument law" q w p |w (t) , samples of the target law are obtained by the following iterations.

1. Sample a proposition w p ∼ q w p |w (t) . 2. Compute the probability

ρ = min    f (w p ) f (w (t) ) q w (t) |w p q (w p |w (t) ) , 1    . (54)
3. Take

w (t+1) =    w p with ρ probability w (t) with 1 -ρ probability . (55)
At convergence, the samples follow the target law f (w) [START_REF] Robert | Monte-Carlo Statistical Methods[END_REF][START_REF] Brémaud | Markov Chains. Gibbs fields, Monte Carlo Simulation, and Queues[END_REF]. When q w p |w (t) = q(w p ) the algorithm is named independent Metropolis-Hastings. In addition, if the instrument law is uniform, the acceptance probability gets simpler in 46). The symbol '×' is the minimum and the symbol '.' is the estimated value by our approach. The y-axis of γ ǫ and γ 1 are in logarithmic scale. . Pseudo-code algorithm. gamrnd, rand and randn draw samples of gamma variable, uniform variable, and zero-mean unit-variance white complex Gaussian vector respectively.

ρ = min f (w p ) f (w (t) ) , 1 . (56)

[19 21 ]

 21 (w α ) and p(w β ) = U[6 8] (w α ) corresponding to an uncertainty of about 5% and 15% around the nominal value (see Sec 3.D).• a rotation parameter ϕ set to π/3. The a priori law is also uniform p(ϕ) = U [π/4 π/2] (ϕ) corresponding to 50% uncertainty.Then, the convolution is computed in the Fourier space and the data are obtained by adding a white Gaussian noise with precision γ ǫ = 0.5. Data are shown Fig.1(b): they are naturally smoother than the true image and the small fluctuations are less visible and corrupted by the noise.

Fig. 1 .

 1 Fig. 1. The figure 1(a) represents a 128 × 128 sample of the a priori law for the object with γ 0 = 1 and γ 1 = 2. Fig. 1(b) is the data computed with the PSF shown in Fig. 2. Figs. 1(c) and 1(d) are the estimates with non-myopic and the myopic estimate, respectively. Profiles correspond to the 68-th line.

Fig. 2 .

 2 Fig. 2. PSF with w α = 20, w β = 7 and ϕ = π/3. The x-axis and y-axis are reduced frequency.

Fig. 3 .

 3 Fig.3. Circular average of the empirical power spectral density of the image, the convolued image, the data (convolued image corrupted by noise) and the estimates, in radial frequency with y-axis in logarithmic scale. The x-axis is the radial frequency.

Fig. 4 .

 4 Fig. 4. Computation of the best parameters in the sense e Eq. (46). The symbol '×' is the minimum and the symbol '.' is the estimated value by our approach. The y-axis of γ ǫ and γ 1 are in logarithmic scale.

  γ 1 for myopic case

Fig. 5 .

 5 Fig. 5. Histograms and chains for the non-myopic case in Figs. 5(a)-5(c) and the myopic case in Figs. 5(b)-5(d) for γ ǫ and γ 1 , respectively. The symbol × localizes the initial value and the dashed line corresponds to the true value. The x-axis are iteration's index for the chains and parameter value for the histograms.

Fig. 6 .

 6 Fig.6. Histogram and chain for the PSF parameters w α in Fig.6(a), w β in Fig.6(b) and ϕ in Fig.6(c). The symbol × localizes the initial value and the dashed line corresponds to the true value. The x-axis for the histograms and the y-axis of the chain are limits of a priori law.

Fig. 7 .

 7 Fig. 7. Joint histograms for the couple (γ 1 , w α ) and (γ 1 , w β ) in Figs. 7(a) and 7(b) respectively. The x-axis and y-axis are the parameter value.

Fig. 8 .

 8 Fig. 8. Observed image Fig. 8(a) and restored image Fig. 8(b). Profiles correspond to the 68-th line. The solid line is the true profile. Dashed line correspond to data in Fig. 8(a) and estimated profiles in Fig. 8(b).

  , γ (0) , w (0) , k = 0 gamrnd(α 0 , β 0 ) % Sample of w 9:w p ← rand * (ab) + a if log(rand) < min{J, 0} then 12: w (k) ← w p 13:Λ H ← Λ H,wp 14:

19 :

 19 Fig.9. Pseudo-code algorithm. gamrnd, rand and randn draw samples of gamma variable, uniform variable, and zero-mean unit-variance white complex Gaussian vector respectively.

Table 1 .

 1 Error e (Eq. (46)) and averaged standard deviation σ of the posterior image law. The "Best" error has been obtained with the knowledge of the true image.

		Data	Non-myopic Myopic	Best
	Error (e)	11.092 %	6.241 %	6.253 % 6.235 %
	σ of x law	-	3.16	3.25	-

Table 2 .

 2 Quantitative evaluation: true and estimated values of hyperparameters and PSF parameters.

	γ ǫ ± σ	γ 1 ± σ	w α ± σ	w β ± σ	ϕ ±	σ

Table 3 .

 3 Acceptation rate.

	Parameter	w α	w β	ϕ
	Acceptation rate 14.50 % 9.44 % 2.14 %

Table 4 .

 4 Specific laws obtained as limit of the Gamma pdf.