Marginalized myopic deconvolution of adaptive optics corrected images using Markov chain Monte Carlo methods

Abstract. Adaptive optics (AO) corrected image restoration is particularly difficult, as it suffers from the lack of knowledge on the point spread function (PSF) in addition to usual difficulties. An efficient approach is to marginalize the object out of the problem and to estimate the PSF and (object and noise) hyperparameters only, before deconvolving the image using these estimates. Recent works have applied this marginal myopic deconvolution method, based on the maximum a posteriori estimator, combined with a parametric model of the PSF, to a series of AO-corrected astronomical and satellite images. However, this method does not enable one to infer global uncertainties on the parameters. We propose a PSF estimation method, which consists in choosing the minimum mean square error estimator and computing the latter as well as the associated uncertainties thanks to a Markov chain Monte Carlo algorithm. We validate our method by means of realistic simulations, in both astronomical and satellite observation contexts. Finally, we present results on experimental images for both applications: an astronomical observation on Very Large Telescope/spectro-polarimetric high-contrast exoplanet research with the Zimpol instrument and a ground-based LEO satellite observation at Côte d’Azur Observatory’s 1.52 m telescope with Office National d'Etudes et de Recherches Aérospatiales’s ODISSEE AO bench.


Introduction
Ground-based high angular resolution imaging in the visible has numerous applications, such as astronomy and satellite observation.The observations are limited by atmospheric turbulence, which can be corrected in real time by adaptive optics (AO).However, the correction is partial and residual blurring remains, impacting high spatial frequencies of the observed object.Therefore, the observation system includes post-processing to restore the high frequencies. 1he residual blurring is described by the system point-spread function (PSF), which is not entirely known, so both the observed object and the PSF are estimated.The historical way to proceed is to estimate them jointly, 2 which leads to a degenerate solution 3,4 in the absence of strong constraints, leading to the sharpest PSF thus the smoothest object.Another way to proceed is to first estimate the PSF by "marginalizing" over the object, i.e., by integrating the joint probability density function over all possible objects with a given prior probability density function and then to deconvolve the image with the estimated PSF. 3 In our case, the PSF has a physical parametric model, and the object is described by a Gaussian prior with a parametric model for its power spectral density (PSD), whose parameters are also estimated along with PSF parameters.The method we have been using so far, AMIRAL (standing for automatic myopic image restoration algorithm), combines PSF and PSD parametrization as well as a marginal maximum a posteriori (MAP) estimator. 5he method we propose in the present paper uses another Bayesian estimator, which minimizes the mean square error on the sought parameters, corresponding to the mean of the marginal posterior distribution.From this method, we also infer uncertainties on PSF and PSD parameters. 6,7To do so, we include prior distributions for PSF and PSD parameters and we compute the marginal posterior distribution by stochastic sampling.We first introduce a Markov chain Monte Carlo (MCMC) algorithm to sample this posterior distribution. 8,9Then, we validate our method on simulated data and we finally apply our method to experimental data, for both astronomical and satellite observation contexts, corresponding to different instruments and turbulence conditions.
2 Imaging Model and MMSE Estimator

Imaging Equation
We consider that the image i results from the 2D discrete convolution of the object o with the PSF h, to which noise n (mostly photon noise and detector readout noise) is added, giving the following imaging model 10 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 4 ; 4 6 0 This can also be written in the matrix form as i ¼ Ho þ n, with H the convolution matrix corresponding to the convolution of the object by h.In this study, we simulate and restore both astronomical and satellite AO-corrected images, implying two different contexts: astronomical images taken on a Very Large Telescope (VLT)/spectro-polarimetric high-contrast exoplanet research (SPHERE)-like instrument 11 with the Zimpol imaging polarimeter 5 and ground-based satellite observation at the Côte d'Azur Observatory (OCA) with the ODISSEE AO system. 12

PSF Model
Throughout this work, we will consider having long-exposure PSFs, meaning that the exposure time is greater than the typical variation time of turbulence.For the PSF, we use the PSFAO19 model, 13 which has been designed specifically for describing an AO-corrected PSF with few physical parameters.Roddier 14 shows that the long-exposure PSF h re-writes as the convolution of three PSFs E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 4 ; 2 8 0 The first one is called the "static" PSF h stat and corresponds to the static aberrations, the second one h det is the "detector" PSF describing the integration of the image on the detector's pixels, and the last one is the "atmospheric" PSF h atm corresponding to the impact of atmospheric turbulence.Conan 15 shows that this description is still valid in the case of an AO-corrected PSF.Both static PSF and detector contributions are supposed well known (and static) compared to the highly variable atmospheric PSF.h atm can then be described by the phase PSD W ϕ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 4 ; 1 8 3 with ϕ the turbulence-induced phase, coming from the residual aberrations [not corrected by AO] in the pupil of the telescope, and σ 2 ϕ ¼ FT −1 ðW ϕ ð0ÞÞ, the phase variance, so that The two main parameters of W ϕ , thus of the PSFAO19 model, are the Fried parameter r 0 taken at the imaging wavelength (850 nm), describing the turbulence's strength, and the variance of the residual turbulent phase v ϕ , describing the quality of AO correction.Indeed, the residual phase PSD W ϕ can be separated in two different spatial frequency zones, depending on the AO cutoff frequency f AO E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 2 ; 1 1 7 ; 6 2 4 For the corrected spatial frequencies, a Moffat model is used to describe the core of the PSD.The main parameter is the amplitude A, which is very close to the residual phase variance: v ϕ ≈ A þ CA AO , with A AO the AO-corrected area in the spatial frequency domain (for a circular AO-corrected area, A AO ¼ πf 2 AO ).C is a constant giving the AO-corrected phase PSD background, useful to model the AO-corrected PSD near AO cutoff frequency (where the Moffat function is close to zero).The parameters α (giving the width of the Moffat function) and β (Moffat's power law) do not impact the computation of the residual phase variance, thus they have a less important impact on the PSF.Throughout this work, α, β, and C will be considered as known parameters, and their value will be fixed to the values in Table 1.Finally, N α;β is a normalization factor, which is used to normalize the integral of the Moffat function over the corrected spatial frequencies E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 7 ; 4 5 5 which requires that β > 1.For the high spatial frequencies, the theoretical Kolmogorov model of turbulence is used, the main parameter being the Fried parameter r 0 describing the turbulence's strength, at the imaging wavelength.This model has been validated on several AO systems and on different telescopes. 12,13

Prior Distributions
Noise is taken independently from the object and is approximated as zero-mean, additive, white, and Gaussian, which is a fine description given the flux levels in typical images.Moreover, in this paper, we approximate the noise precision (inverse variance) as homogeneous, and we denote it by γ n .Therefore the noise covariance matrix is R n ¼ 1∕γ n I, with I the identity matrix and the noise PSD is S n ¼ 1∕γ n .An example of simulated astronomical observation is given in Fig. 1, with the true object on the left and the simulated image on the right.The image is simulated using the PSFAO19 model, and with uniform zero-mean additive white Gaussian noise.
As a prior for the object, we consider a Gaussian model described by its mean m o and its PSD S o .Given that we have little information on the the mean object m o , it is taken uniform on all pixels, estimated at the average value of the image considering that P h ¼ 1, supposing flux conservation.For the object PSD, we use the following parametric model: and f ¼ jf j the radial frequency.This circularly symmetric model is a slightly modified writing of Matérn's model. 16,17In this model, γ o sets the global PSD level, p is the PSD decrease rate at high frequencies, and k gives the breakpoint between the two regimes of the model.In previous works, 5 attempts to estimate hyperparameter p jointly with the other parameters have been shown to strongly decrease PSF parameter estimation accuracy.Therefore, we choose to work in a "mostly unsupervised" mode, where p is fixed to a standard value.In the case of astronomical observations of asteroids, a well-fitting empirical value is around p ¼ 3, whereas for satellite observation a standard value for p would be around 2.5 to 2.6.Let P ¼ N 2 the image size in pixels.Given previous assumptions and the matrix form of the imaging model of Eq. ( 6) where we consider i and o as vectors and H as a P × P matrix, the likelihood writes E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 4 ; 4 9 3 pðijo; In Eq. ( 1), we consider o and h as arrays of the same size and approximate them as periodic.Thus, the likelihood in Eq. ( 6) can also be written in the Fourier domain given previous approximations E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 4 ; 4 0 2 pðijo; where : denotes the discrete Fourier transform (DFT) and the product on f is on all pixels in the spatial frequency domain.For the object and the image, the DFT is normalized so as to comply with Parseval's theorem.For the PSF, the DFT is normalized so that hð0Þ equals the sum of the PSF on the numerical array.Moreover, as said previously, this value is set to 1 by convention, to express flux conservation.h is also called the (discretized) optical transfer function (OTF).Given the previous approximations, the object covariance matrix is circulant-block-circulant and we can write the object prior distribution as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 4 ; 2 7 4 Similarly to the likelihood, the object prior in Eq. ( 8) can also be written in the Fourier domain.Indeed, given the structure of the object covariance matrix R o , the latter is diagonalized in the discrete Fourier domain, with S o as written in Eq. ( 5) on its diagonal, so that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 4 ; 2 0 0 Thus, given that the mean object is taken uniform in the spatial domain, it corresponds to a delta function in the Fourier domain.Regarding PSF parameters as well as noise and object PSD parameters, hereafter called parameters, we consider that each parameter γ n , γ o , k, r 0 , and v ϕ can take any value in a given range.Therefore, following the Laplace rule (or principle of insufficient reason), we use uniform priors for each of them. 18The prior interval for each parameter is taken large enough.For γ n , γ o , and k, the prior interval is from 0.1 to 10 times the usual value of the parameter (given empirical knowledge on them).The prior intervals taken for PSF parameters are the following: for r 0 we take [5 cm; 30 cm] and for v ϕ we take ½0.5 rad 2 ; 3.0 rad 2 , which correspond to a large range of values taking into account the global knowledge on the AO system and the turbulence.These values are summarized in Table 2.
In Fig. 2, we provide the chosen hierarchical model, which sums up the variable interdepend-ency.[Ref.19, chap.8] Each upper node (parent) is connected with an edge to a node below (child), and the model says that a child's distribution, given all nodes above, is only a function of its parents.In our model, it means for example that pðijo; Therefore pði; γ o ; kjo; γ n ; r 0 ; v ϕ Þ ¼ pðijo; γ n ; r 0 ; v ϕ Þpðγ o ÞpðkÞ, which means that the image i and object PSD parameters γ o and k are independent conditionally to the object and the other parameters.Additionally, the object, the noise variance, and the PSF parameters are independent conditionally to object PSD parameters meaning pðojγ n ; γ o ; k; r 0 ; v ϕ Þ ¼ pðojγ o ; kÞ.Moreover, as the hierarchical model reads, all parameters (γ o , k, r 0 , v ϕ and γ n ) are modeled as a priori independent.

Marginal Estimator
The joint distribution is, following the conditioning rule, the multiplication of the likelihood by the prior distributions.Given the hierarchical model of Fig. 2 and as explained previously From it, we can derive the expression of any target distribution.As explained above, a way to estimate the object and the parameters is to first estimate the parameters by computing the so called marginalized posterior probability, meaning integrating the posterior density over the object E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 7 ; 1 6 3 In practice, we write the marginal posterior distribution following the Bayes rule, from the marginal likelihood and the priors taken as uniform and independent, as mentioned in Sec.2.3: e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 7 ; 1 0 2  Given that the noise is taken Gaussian, white, homogeneous and a priori independent from the object considered Gaussian, the image being a linear combination of both is also Gaussian.Therefore, the marginal likelihood writes: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 4 ; 7 0 0 with image PSD S i and mean image m i E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 4 ; 6 4 7 [Given Eq. ( 13), maximizing the marginal likelihood, as in the previous method, 3,5 can be interpreted as finding the parameters for the image PSD model S i of Eq. ( 14) that best fit the empirical PSD j ĩðf Þ − mi ðf Þj 2 .]The marginal posterior distribution for the parameters can then easily be written using Eqs.( 12)-( 14) with U x ðxÞ the uniform probability distribution for parameter x, in the range defined for x.In our case, these minimum and maximum values are given in Table 2.

MMSE Estimator and Sampling
The minimum mean square error estimator is known to be the mean of the posterior distribution, whereas the MAP estimator, computed by AMIRAL, is its mode. 20Given the complexity of the posterior, there is no known analytical way to calculate it.A way to compute it is to draw samples under the posterior distribution using a MCMC method for instance and compute the sample mean.The posterior distribution being complex, it is not possible to sample it directly, therefore we use a Metropolis-Hastings algorithm to bypass the problem. 21It consists, for each iteration, in drawing samples under a chosen proposal distribution qðθjθ 0 Þ and accepting the samples (else, duplicating the previous value) with a prescribed probability α.For the k'th iteration, α writes E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 4 ; 2 9 8 α ¼ pðθ prop jiÞ pðθ ðk−1Þ jiÞ with θ prop a sample drawn from the proposal distribution.Several versions are possible: in particular, we can either draw all the parameters simultaneously (standard Metropolis-Hastings) or separately (Metropolis-Hastings-within-Gibbs).Drawing the parameters together can make the acceptance probability fall (except if we use more advanced, e.g., gradient-based algorithms, such as Metropolis-adjusted Langevin algorithm or Hamiltonian Monte Carlo methods), [21][22][23] whereas drawing parameters individually can slow down the algorithm as it changes parameters one by one and requires more likelihood computations.For simplicity, we use here the second version.In a standard Gibbs algorithm, each parameter is drawn under its own conditional posterior distribution, which is proportional to the prior of the considered parameter times the marginal likelihood of Eq. ( 13).The conditional posterior distribution for each parameter writes E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 4 ; 1 2 8 where θ n is the considered parameter and θ m≠n the four other parameters.Note that pðiÞ is not needed due to the fact that it cancels out in the acceptance ratio α computed in Eq. ( 16).
As mentioned above, because drawing the parameters under their conditional posterior distribution is difficult, we use a Metropolis-Hastings-within-Gibbs algorithm instead of the standard Gibbs.Asymptotically, the samples are under the marginal posterior distribution for all parameters, and the sample mean tends towards the mean of the distribution. 21n our case, we use a random walk Metropolis-Hastings algorithm: the proposed sample for each parameter is drawn under a symmetric (Gaussian) distribution qðθ n jθ n 0 Þ around the current value of the parameter.The proposal distribution being symmetric, qðθ in α [in Eq. ( 16)] simplifies.The standard deviation of the Gaussian proposal distribution σ θ n is chosen to be 0.01 times the allowed range of the prior.Precisely, the tuning for each parameter is given in Table 2.This choice is based on the empirical sensitivity of the PSF (or noise PSD, or object PSD) to the parameters.Moreover, this choice only impacts the convergence time, which is an issue we do not tackle in this paper.A typical number of iterations needed to reach convergence (including the boiling time) would be around 30,000 iterations, corresponding to an hour (on an ordinary laptop).The random walk Metropolis-Hastings-within-Gibbs algorithm we use in this paper is provided in Algorithm 1.
3 Results on Simulated Astronomical Data

Simulation Conditions
The obtained results are shown for the simulated image displayed in Fig. 1, using as the true object the synthetic view of asteroid Vesta, built by the OASIS software, 24 on a dark background of size 512 × 512 pixels.True PSF parameters are r 0 ¼ 0.15 m and v ϕ ¼ 1.3 rad 2 at the imaging wavelength λ ¼ 550 nm, which correspond to realistic turbulence and correction conditions.The AO system is a "SPHERE-like" AO system and its parameters are taken identical to those used with the previous method, 5 for comparison purposes.Noise is taken zero-mean, additive, white, and Gaussian with a variance equal to the empirical mean value of the object as a first approximation of the photon noise.The total flux of the object is set to F o ¼ 10 9 ph (photons), typical from VLT/SPHERE/Zimpol asteroid observations (ESO Large Program ID 199.C-0074), therefore γ n ¼ P∕F o ¼ 2.62 × 10 −4 ph −2 .The PSF and PSD parameters are estimated following the proposed method, except the mean object m o , which are estimated to the average value of the image, and the object PSD power, which is fixed to p ¼ 3, which corresponds to a reasonable default value of p for asteroids.The Gibbs sampler is run for 100,000 iterations, which corresponds to a few hours, to verify the convergence.

Results on the Estimated Parameters and Derived Uncertainties
In Fig. 3, we plot the samples chains and the corresponding histograms for γ n , r 0 , and v ϕ .The inspection of Fig. 3 suggests that chains have a short burn-in period, followed by a stationary state.As expected from Markov chains, for each parameter the samples are correlated.Moreover, the samples are concentrated in a small interval relatively to their prior interval.
The sample mean values m, corresponding to our estimates, and standard deviations σ, corresponding to our predicted uncertainties, for each parameter are displayed in Table 3. First, we can note that the error on the parameters is small: the noise precision is very precisely estimated, with an error smaller than 0.2%, and PSF parameters are also well estimated, with a 5% error on r 0 and a 10% error on v ϕ .Additionally, the estimated r 0 and v ϕ are very close to the previous results obtained with AMIRAL: for similar conditions, 5 the estimated PSF parameters were r 0 ¼ 0.142 m and v ϕ ¼ 1.13 rad 2 (compared to r 0 ¼ 0.142 m and v ϕ ¼ 1.17 rad 2 in Table 3).

Results on the OTF, on Object and Image PSD
We also compare the resulting OTF to the true OTF in Fig. 4. The slight underestimation of r 0 leads to the lowering of the global OTF level, the impact of which can mainly be seen at low frequencies.Concerning v ϕ , its mild underestimation leads to a slower decrease of the OTF and impacts the slope of the latter at medium-high frequencies. 5Thus, we notice that the errors on  Concerning the uncertainties derived from our method, we notice in Table 3 that the true value for parameter r 0 is in the range ½m r 0 AE 2σ r 0 , and the true v ϕ is in the interval ½m v ϕ AE 5σ v ϕ , therefore the uncertainties on PSF parameters seem under-estimated.We can also compute uncertainties directly on the sought OTF: for each sample ðr 0 ; v ϕ Þ, we compute the corresponding OTF to compute its sample mean m h and standard deviation σ h, meaning the mean and standard deviation for each frequency of the OTF.As shown in Fig. 4, the true OTF is within the interval ½m h AE 2σ h, for all frequencies.Therefore, even though the uncertainties on PSF parameters are somewhat under-estimated, our method gives a very satisfactory uncertainty estimation on the OTF itself.
In Fig. 5(a), we perform an important sanity check of the method to verify that our model for the image PSD of Eq. ( 14), which combines object PSD, PSF and noise PSD, accurately fits the empirical image PSD averaged azimuthally [cf.Eq. ( 13)].Moreover, given the fact that the true object is not the realization of a Gaussian random field following our PSD model, a way to check γ o and k's estimation accuracy is to look at the fitting of our model to the object empirical PSD, averaged azimuthally.As displayed in Fig. 5  empirical object PSD, the slight overestimation being consistent with the slight underestimation of the OTF.

Results on the Restored Image
After having estimated the PSF and hyperparameters, the object is restored by maximizing the joint distribution, given the PSF and hyperparameters, as in a classical non-myopic deconvolution framework.Given the expression of the joint distribution in Eq. ( 10), maximizing it with respect to the object is equivalent to maximizing the product of the likelihood in Eq. ( 7) and the object prior (corresponding to a quadratic regularization) in Eq. ( 9) Without any specific constraint on the object, the solution of Eq. ( 18) corresponds to the Wiener filtering [Ref. 1, chap.4] with a non-null prior mean m o E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 1 1 4 ; However, it is also possible to minimize the criterion in Eq. ( 18) under some constraints, such as positivity (on the pixels value), which is a natural constraint given the context.It is also possible to use not solely a quadratic regularization but a L1-L2 regularization as an "edge-preserving" regularization.
Figure 6 shows the image in Fig. 1 restored with the estimated OTF, using a quadratic regularization (which hyperparameters are the ones estimated by the method) with positivity constraint.Many details of the Vesta surface can be seen, that were not visible on the data.Particularly, with our method we retrieve sharp edges of the asteroid from which one can estimate the object volume and sphericity, as well as main crater and albedo features.

Posterior Coupling Between Parameters
Sampling the whole posterior distribution, instead of computing a single point of it (for example, the maximum), enables us to study the a posteriori coupling of the parameters.In Figs. 7 and 8, we display the scatter graph of the samples, after boiling time, for two different couples of parameters: (r 0 , v ϕ ) and (r 0 , γ o ).Most couples of parameters have a scatter graph similar to Fig. 7, where the 2D-histogram is rather along the axis suggesting that most parameters are not correlated a posteriori.
The only pair of coupled parameters that does not have an elliptical-like scatter graph, but instead shows a strong a posteriori correlation, is r 0 and γ o .We explain this correlation by the fact that as shown in Ref. 5, r 0 impacts the global level of the OTF whereas γ o gives the global level of the object PSD.Therefore, given the expression of the image PSD in Eq. ( 14), both (r 0 , γ o ) have a similar impact on the global level of the image PSD, which is fitted by our method; that explains their strong correlation.

Test on Several Noise Realizations
To test the robustness of our method to noise, we ran the algorithm for 10 different noise realizations, in the simulation conditions described above.We compute the bias and standard deviation of the estimated parameters on these 10 noise realizations, as well as the maximum error.We also compute the minimum and maximum predicted uncertainty (i.e., the standard deviation of the posterior distribution).These values are summed up in Table 4.
In these 10 cases, we notice very little variation in the estimates: the computed standard deviations (fifth column in Table 4) are small with respect to the true values (second column).Moreover, the estimates are satisfactory: first, the errors on the estimated parameters are quite  column), particularly on the noise precision (error is <1%).Moreover, the predicted uncertainties for γ n are close to the empirical average error made on the noise precision.For the PSF parameters, the error is smaller than 11%.Concerning parameter r 0 , the predicted uncertainty is very satisfactory: the true value is always within the interval ½m r 0 AE 2σ r 0 .For parameter v ϕ , we notice that the error is here dominated by the bias, which is more than 10 times greater than the standard deviation (which is not the case for the other parameters).Our interpretation is that this bias is due to the choice of p, which will be further discussed in 3.7.
Finally, even though the uncertainties are under-estimated for v ϕ with the default p, concerning the OTF itself the uncertainties are always well estimated: for all 10 cases, the true OTF is within the interval ½m h AE 2σ h, as shown in Fig. 9.Moreover, the RMSE on the OTF is smaller than 1.3 times the posterior standard deviation on OTF (averaged on noise realizations), for all frequencies.

Impact of Hyperparameter p
In Ref. 25, we have tested our method in the exact same conditions, for another value for hyperparameter p, tuned slightly differently, towards the "best" value 5 in the supervised mode p ¼ 2.91.(Both p ¼ 2.9 and p ¼ 2.91 were tested with our method, giving the same results.)The results on estimated parameters and derived uncertainties are summed up in Table 5.With a value of p ¼ 2.9 instead of 3.0, the error on v ϕ becomes smaller and the estimated uncertainties are then satisfactory (the true v ϕ is then in the interval ½m v ϕ AE 2σ v ϕ ).Similarly to the correlation between r 0 and γ o showed in 3.5, we interpret these results as another strong correlation between v ϕ and the fixed hyperparameter p, due to the similar impact they have on the image PSD.Indeed, v ϕ impacts the slope of the OTF in medium-high frequencies 5 whereas p corresponds to the slope of the in medium-high frequencies.Therefore, both v ϕ and p tune the decrease of the image PSD in medium-high frequencies, which can explain their strong posterior correlation.However the differences on the restored image, as displayed in Figs. 10 and 11, are quite small, at most around 10 times smaller than the global image level.with a More Realistic Noise We now simulate the observation of Vesta using the same conditions than in Sec.3.1, except that we now simulate a more realistic noise, using a Poisson distribution to mimic the photon noise with the same total flux as previously, namely 10 9 ph, and an additive stationary Gaussian noise for the read-out noise with a standard deviation of 20 photo-electrons.The results are summed up in Table 6.
Even if the simulated noise does not exactly match the stationary Gaussian noise model, all estimated parameters are still very close to the previous estimations, with a difference between the two estimations smaller than the derived uncertainties σ θ , except for γ n , which does not have a true value because of the inhomogeneous simulated noise.The PSF parameters are still well estimated, with an error around 3% for r 0 and around 8% for v ϕ .Their associated uncertainties are also still satisfactory for r 0 and, similarly, slightly underestimated for v ϕ as discussed previously.Thus deviating from the stationary Gaussian noise model does not impair the results with our method, given our simulation conditions.The small impact of deviating from the stationary Gaussian noise model was also shown by Fétick, 5 using the marginal MAP estimator.
The prior mean object chosen in this work can also be questioned: indeed, taking a uniform prior mean object equal to the mean image makes sense as we have little information on it, but this choice is somewhat arbitrary, and above all, it depends on the data.To check that this choice has little impact on the solution, we perform another reconstruction with a Poisson + Gaussian noise, but changing this time the prior mean object value, which is estimated to zero (and not the image mean value).The results are given in Table 7.
The results for all parameters do not change much, as previously the estimates for each parameter change less than a standard deviation σ, except for γ n .We thus conclude that this uniform prior on the mean object does not impact much the results on the estimated parameters.Simulation Conditions We now show results for a simulated satellite image, using as the true object a synthetic view of the SPOT satellite on a dark background of size 512 × 512 pixels. 26We simulate its observation using the ODISSEE AO system at OCA, 12 and with true PSF parameters r 0 ¼ 0.10 m and v ϕ ¼ 1.85 rad 2 , at the imaging wavelength λ ¼ 850 nm, corresponding to a stronger turbulence, and to a more modest correction than for the astronomical simulation because of a less complex AO system.The noise is taken as zero-mean, additive, white and Gaussian, and its variance is taken equal to the mean value of the object.Here, the mean flux is 10 4 photons per image pixel, corresponding to a somewhat optimistic value.The pixel sampling is close to the Shannon-Nyquist criterion, with slightly more than two pixels per λ∕D.The true object and the simulated data are shown in Fig. 12.
The object PSD power p is fixed to an empirical standard value for satellites p ¼ 2.6, to fit the empirical object PSD.The Gibbs sampler is run for 100,000 iterations.

Results on the Estimated Parameters and Derived Uncertainties
In Fig. 13, we plot the sample chains and the corresponding histograms for γ n , r 0 , and v ϕ .
Similarly to the previous simulations, the sample mean values m and standard deviations σ of the posterior distribution for each parameter are displayed in Table 8.The noise precision γ n as well as PSF parameters r 0 and v ϕ are relatively well estimated, with an error of 2%, 14%, and 17%, respectively.These results are very close to those obtained with AMIRAL: for similar conditions, the estimated PSF parameters are r 0 ¼ 0.112 m and v ϕ ¼ 2.16 rad 2 .We notice that the error on PSF parameters is greater for the satellite observation than for the astronomical observation.Our interpretation of these results, checked by complementary simulations, is that it is due to the spectrum of the satellite object which is less isotropic than Vesta, and therefore does not fit our isotropic PSD model as well.
Concerning uncertainties, similarly to previous asteroid case, the posterior standard deviation for r 0 seems to be a good uncertainty prediction for this parameter as the true value is within the interval ½m r 0 AE 2σ r 0 .On the contrary, σ v ϕ seems small, giving an under-estimated uncertainty.The reasons for this under-estimation are being investigated, though it should be linked to the more difficult observation conditions simulated here.Additionally, the previous discussion about the correlation between parameters p and v ϕ 25 is still valid and further work on tuning hyper-parameter p for satellite observation should be done.
We also compare the resulting estimated OTF to the true OTF in Fig. 14.Here again, as discussed previously, we notice that the errors on both parameters partially compensate, as a result the normalized RMSE for the OTF is quite low (around 8%).Concerning the uncertainties, we notice again that even though the uncertainties on PSF parameters are under-estimated, the uncertainties on the OTF itself are quite satisfactory as the true OTF is within the interval ½m˜h AE 3σ h.estimations result in a good image PSD fitting, as shown in Fig. 15.Moreover, as displayed in Fig. 15, the object PSD model visually fits well the empirical object PSD.Finally, Fig. 16 shows results from the restoration of the image in Fig. 12 (right) using the estimated OTF.We notice that details of the satellite surface are restored.

Results on Experimental Astronomical
testing our method on both astronomical and satellite simulated data, therefore for different turbulence conditions and AO systems, we apply it to experimental images.Here we process an experimental image of Vesta 27 taken by SPHERE/Zimpol in the same mostly unsupervised mode as previously where p ¼ 3, and run the Gibbs sampler for 100,000 iterations.Data and restored object are shown in Fig. 17.We recognize the same surface features as from the synthetic view in Fig. 1.In this experimental case, the bright edge corona starts to appear (on the left side), and the image is slightly granular.This may be due to a slight over-deconvolution i.e., to a slight under-estimation of the OTF, as to the quadratic regularization.
Results obtained for the PSF parameters (mean ± standard deviation) are the following: r 0 ¼ 0.26 AE 0.04 m and v ϕ ¼ 2.62 AE 0.06 rad 2 .These values are close to the values obtained with AMIRAL 5 (r 0 ¼ 0.32 m, v ϕ ¼ 2.78 rad 2 ) for the same conditions, the newly estimated r 0 being more likely than the one estimated by AMIRAL according to the known statistics on r 0 . 13e also look at the image PSD model and the empirical image PSD for Vesta in Fig. 18.They fit well, especially at low and medium frequencies, where signal dominates noise.For high frequencies, where the noise is dominant, we see that the noise floor is not flat (whereas we model the noise as white), and believe it may be due to the data reduction by SPHERE/ Zimpol's pipeline.

on an Experimental Satellite Image
Finally, we test our method on an experimental image of Envisat, taken at the OCA with Office National d'Etudes et de Recherches Aérospatiales (ONERA)'s ODISSEE AO system. 12We fix paramater p to a reasonable value for satellites (p ¼ 2.5) and run the Gibbs sampler for 500,000 iterations.
Our results on the estimated PSF parameters are the following: r 0 ¼ 0.08 m and v ϕ ¼ 0.89 rad 2 , which can be compared to the results obtained with AMIRAL: r 0 ¼ 0.06 m and v ϕ ¼ 1.13 rad 2 , for the same conditions.Concerning the restored image, as shown in Fig. 19, we retrieve some elements of the satellite, and we checked on a computer aided design model of Envisat that the bright spots we obtain on the restored image indeed correspond to instruments and antennas on its surface.
Finally, concerning the image PSD model and the empirical image PSD for Envisat, as shown in Fig. 20, we have a globally good image PSD fitting.The oscillations of the empirical image PSD are likely to come from oscillations of the OTF, which are consistent with the exposure time (≈1 s), which is short with respect to turbulence residuals averaging, and constitutes a deviation to the infinite exposure assumption of our AO-corrected PSF model.These oscillations might also come partly from the spectrum of the object itself, which deserves further studies by means of simulations.

Conclusion
We have presented myopic deconvolution method extending previous works, using a MCMC algorithm, more precisely a random walk Metropolis-Hastings-within-Gibbs algorithm.In addition to PSF and hyperparameter estimation combined with image restoration, we now have access to the whole posterior distribution.This enables us to compute the optimal estimator minimizing the mean square error.Additionally, with the posterior distribution, we can compute uncertainties based on the posterior standard deviation.This method has been validated on simulated images, giving accurate estimations of noise and object hyperparameters, as well as satisfactory OTF estimations.Two different contexts were simulated: on the one hand, the observation of asteroid Vesta on the VLT, and on the other hand, the observation of the SPOT satellite using ONERA's AO bench on the 1.52 m-telescope at OCA.The satisfactory results obtained in both conditions suggest the broad applicability of the method.Additionally, for the simulated asteroid images, we have computed our estimations for several noise realizations, to check the robustness of our method to noise, both for estimated parameter values and predicted uncertainties.Finally, our method has also been applied to experimental images, in both contexts.
In this work, hyperparameter p, which codes for the decrease of the object PSD, has been fixed to a reasonable value according to the class of the object (either asteroids, or satellites).The PSF estimation quality is sensitive to the choice of p, as we verified it by changing its value. 25oreover, jointly estimating p with the other parameters is difficult as mentioned in earlier studies. 5In the near future, we plan to tackle the joint estimation of p. Beyond the choice  has worked since 2005.He is involved in various research activities in AO applied to space observation, astronomy, free space optical telecommunication, and biomedical imaging.His research interests include design, numerical simulation, experimental implementation, and testing of AO systems, with a focus on control.He is now in charge of the development of the ONERA optical ground station for GEO-feeder links, FEELINGS.

Fig. 2
Fig.2Hierarchical model summing up the inter-dependency between the object, the image, and all parameters.

Algorithm 1
Metropolis-Hastings-within-Gibbs algorithm.Define initial θ ð0Þ for each iteration k do for each parameter θ n do Propose θ prop n

Fig. 3
Fig. 3 From top to bottom: γ n , r 0 , v ϕ , γ o , k .(a) Chain of samples for simulated astronomical image.(b) Corresponding histogram.True values in dashed line.

Fig. 4
Fig. 4 True (in green) and estimated (in blue) OTF for simulated astronomical image, including computed uncertainties (in blue, + and − for upper and lower uncertainty bounds).

Fig. 5
Fig. 5 PSDs for simulated astronomical image.Model in solid line and empirical PSD averaged azimuthally in dashed line.(a) Image PSD and (b) object PSD.

Fig. 6
Fig. 6 (a), (b) True object and image for simulated asteroid observation, 256 × 256 cropped from Fig. 1.(c) Restored object from the estimated PSF and PSD parameters using a L2-norm regularization, with positivity constraint, also cropped.

Fig. 7
Fig. 7 Marginal posterior scatter graph of the samples for (r 0 , v ϕ ) after boiling time.

Fig. 8
Fig. 8 Marginal posterior scatter graph of the samples for (r 0 , γ o ) after boiling time.

Table 4 7 r 1 ×Fig. 9
Fig. 9 Results on OTF uncertainties for 10 realizations of a noisy simulated astronomical image: true OTF (in black) and predicted range ½mh AE 2σh (+ and − for upper and lower uncertainty bounds, each color corresponds to a noise realization).

Fig. 13 4 r
Fig. 13 From top to bottom: γ n , r 0 , v ϕ , γ o , and k .(a) Chains of samples for simulated satellite image.(b) Corresponding histogram.True values in dashed line.

Fig. 14
Fig. 14 True (in green) and estimated (in blue) OTF for simulated satellite image, including computed uncertainties (in blue, + and − for upper and lower uncertainty bounds).

Fig. 15
Fig. 15 PSDs for simulated satellite image.Model in solid line, and empirical PSD averaged azimuthally in dashed line.(a) image PSD and (b) object PSD.

Fig. 17
Fig. 17(a) Vesta observed by SPHERE/Zimpol on the European VLT in Chile. 27(b) Restored object with the estimated PSF using a L2-norm regularization, with positivity constraint.
Fig. 17(a) Vesta observed by SPHERE/Zimpol on the European VLT in Chile. 27(b) Restored object with the estimated PSF using a L2-norm regularization, with positivity constraint.

Fig. 18
Fig. 18 PSD model and Vesta empirical image PSD averaged azimuthally, in dashed line.

Fig. 20
Fig. 20 PSD model and Envisat empirical image PSD averaged azimuthally, in dashed line.

Table 1
Moffat fixed parameters.

Table 3
Mean value, associated standard deviation and true value, for γ n , r 0 , v ϕ , γ o , and k for simulated astronomical image (stationary Gaussian noise), with p ¼ 3 and m o ¼ m i .

Table 7
Mean value, associated standard deviation and true value, for γ n , r 0 , v ϕ , γ o , and k for simulated astronomical image with a more realistic noise (Poisson + Gaussian noise), with p ¼ 3 and m o ¼ 0.