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

Alix Yan; Laurent M. Mugnier; Jean-François Giovannelli; Romain Fétick; Cyril Petit

doi:10.1117/1.JATIS.9.4.048004

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

Alix Yan, Laurent M. Mugnier, Jean-François Giovannelli, Romain Fétick, Cyril Petit

Author Affiliations +

Journal of Astronomical Telescopes, Instruments, and Systems, Vol. 9, Issue 4, 048004 (November 2023). https://doi.org/10.1117/1.JATIS.9.4.048004

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.

1. 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.¹

The 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,² which leads to a degenerate solution³^,⁴ 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.³ 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.⁵

The 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.⁶^,⁷ To 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.⁸^,⁹ Then, 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

2.1.

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¹⁰

Eq. (1)

i = h * o + n .

This can also be written in the matrix form as $i = H o + 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¹¹ with the Zimpol imaging polarimeter⁵ and ground-based satellite observation at the Côte d’Azur Observatory (OCA) with the ODISSEE AO system.¹²

2.2.

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,¹³ which has been designed specifically for describing an AO-corrected PSF with few physical parameters. Roddier¹⁴ shows that the long-exposure PSF $h$ re-writes as the convolution of three PSFs

Eq. (2)

h = h_{stat} * h_{\det} * h_{atm} .

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¹⁵ 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_{ϕ}$

Eq. (3)

h_{atm} = {FT}^{- 1} (\exp (- σ_{ϕ}^{2}) \times \exp ({FT}^{- 1} (W_{ϕ}))),

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

\sum h_{atm} = 1

. 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}

W_{ϕ} (f) = {\begin{cases} A N_{α, β} (1 + f^{2} / α^{2})^{- β} + C & if f \leq f_{AO} \\ 0.023 r_{0}^{- 5 / 3} f^{- 11 / 3} & else . \end{cases}

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_{ϕ} \approx A + C A_{AO}$ , with $A_{AO}$ the AO-corrected area in the spatial frequency domain (for a circular AO-corrected area, $A_{AO} = π f_{AO}^{2}$ ). $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

Eq. (4)

N_{α, β} ≜ \frac{β - 1}{π α^{2}} {[1 - {(1 + \frac{f_{A O}^{2}}{α^{2}})}^{1 - β}]}^{- 1},

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.¹²^,¹³

Table 1

Moffat fixed parameters.

Parameter	Value
Moffat width $α$ ( $m^{- 1}$ )	0.05
Moffat power law $β$	1.5
AO area constant $C$ ( ${rad}^{2} m^{2}$ )	$10^{- 10}$

2.3.

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.

Fig. 1

(a) Synthetic view of Vesta (true object $o$ ), of size $512 \times 512$ . (b) Simulated image $i$ , with true parameters $r_{0} = 0.15 m$ , $v_{ϕ} = 1.3 {rad}^{2}$ , and $γ_{n} = 2.62 \times 10^{- 4} {ph}^{- 2}$ .

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 $\sum h = 1$ , supposing flux conservation. For the object PSD, we use the following parametric model:

Eq. (5)

S_{o} (f) = \frac{1}{γ_{o}} {\bar{S}}_{o} (f), with {\bar{S}}_{o} (f) = 1 / (k + f^{p}),

and

f = | f |

the radial frequency. This circularly symmetric model is a slightly modified writing of Matérn’s model.¹⁶^,¹⁷ In 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,⁵ 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 \times P$ matrix, the likelihood writes

Eq. (6)

p (i | o, γ_{n}, r_{0}, v_{ϕ}) = \det {(2 π R_{n} (γ_{n}))}^{- 1 / 2} \times \exp (- \frac{1}{2} {(i - H (r_{0}, v_{ϕ}) o)}^{t} R_{n}^{- 1} (γ_{n}) (i - H (r_{0}, v_{ϕ}) o)) .

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

Eq. (7)

p (i | o, γ_{n}, r_{0}, v_{ϕ}) = {(\frac{γ_{n}}{2 π})}^{P / 2} \prod_{f} \exp (- \frac{γ_{n}}{2} | \tilde{i} (f) - \tilde{h} (f) \tilde{o} (f) |^{2}),

where

\tilde{.}

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

\tilde{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.

\tilde{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:

Eq. (8)

p (o | γ_{o}, k) = \det {(2 π R_{o} (γ_{o}, k))}^{- 1 / 2} \exp (- \frac{1}{2} (o - m_{o})^{t} R_{o}^{- 1} (γ_{o}, k) (o - m_{o})) .

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

Eq. (9)

p (o | γ_{o}, k) = {(\frac{γ_{o}}{2 π})}^{P / 2} \prod_{f} ({\bar{S}}_{o} (f)^{- 1 / 2} \exp [- \frac{γ_{o}}{2} {\bar{S}}_{o} (f)^{- 1} {| \tilde{o} (f) - {\tilde{m}}_{o} (f) |}^{2}]) .

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.¹⁸ The 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.

Table 2

Prior intervals and tuning of the Gaussian standard deviation for γn, r0, vϕ, γo, and k.

Parameter	Prior min - max	Step tuning
$γ_{n}$ ( ${ph}^{- 2}$ )	$2.62 \times 10^{- 5}$ to $2.62 \times 10^{- 3}$	$2.62 \times 10^{- 6}$
$r_{0}$ (m)	0.05 to 0.5	0.001
$v_{ϕ} ({rad}^{2})$	0.5 to 3.0	0.01
$γ_{o}$ ( ${ph}^{- 2}$ )	$2.62 \times 10^{- 15}$ to $2.62 \times 10^{- 11}$	$2.62 \times 10^{- 14}$
$K$	0.01 to 10	0.1

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 (i | o, γ_{n}, γ_{o}, k, r_{0}, v_{ϕ}) = p (i | o, γ_{n}, r_{0}, v_{ϕ})$ . Therefore $p (i, γ_{o}, k | o, γ_{n}, r_{0}, v_{ϕ}) = p (i | o, γ_{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 (o | γ_{n}, γ_{o}, k, r_{0}, v_{ϕ}) = p (o | γ_{o}, k)$ . Moreover, as the hierarchical model reads, all parameters ( $γ_{o}$ , $k$ , $r_{0}$ , $v_{ϕ}$ and $γ_{n}$ ) are modeled as a priori independent.

Fig. 2

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

2.4.

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

Eq. (10)

p (i, o, γ_{n}, γ_{o}, k, r_{0}, v_{ϕ}) = p (i | o, γ_{n}, r_{0}, v_{ϕ}) p (o | γ_{o}, k) p (γ_{n}) p (γ_{o}) p (k) p (r_{0}) p (v_{ϕ}) .

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

Eq. (11)

p (γ_{n}, γ_{o}, k, r_{0}, v_{ϕ} | i) = \int p (o, γ_{n}, γ_{o}, k, r_{0}, v_{ϕ} | i) d o = \frac{1}{p (i)} \int p (i, o, γ_{n}, γ_{o}, k, r_{0}, v_{ϕ}) d o .

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:

Eq. (12)

p (γ_{n}, γ_{o}, k, r_{0}, v_{ϕ} | i) = \frac{p (γ_{n}) p (γ_{o}) p (k) p (r_{0}) p (v_{ϕ})}{p (i)} p (i | γ_{n}, γ_{o}, k, r_{0}, v_{ϕ}) .

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:

Eq. (13)

p (i | γ_{n}, γ_{o}, k, r_{0}, v_{ϕ}) = {(2 π)}^{- P / 2} \prod_{f} (S_{i} (f)^{- 1 / 2} \exp [- \frac{1}{2} {| \tilde{i} (f) - {\tilde{m}}_{i} (f) |}^{2} / S_{i} (f)]),

with image PSD

S_{i}

and mean image

m_{i}

Eq. (14)

S_{i} (f) = S_{o} (f) {| \tilde{h} (f) |}^{2} + S_{n} {\tilde{m}}_{i} (f) = \tilde{h} (f) {\tilde{m}}_{o} (f) .

[Given Eq. (13), maximizing the marginal likelihood, as in the previous method,³^,⁵ can be interpreted as finding the parameters for the image PSD model $S_{i}$ of Eq. (14) that best fit the empirical PSD $| \tilde{i} (f) - {\tilde{m}}_{i} (f) |^{2}$ .] The marginal posterior distribution for the parameters can then easily be written using Eqs. (12)–(14)

Eq. (15)

p (γ_{n}, γ_{o}, k, r_{0}, v_{ϕ} | i) = \frac{p (γ_{n}) p (γ_{o}) p (k) p (r_{0}) p (v_{ϕ})}{p (i)} p (i | γ_{n}, γ_{o}, k, r_{0}, v_{ϕ}) = \frac{1}{p (i)} U_{γ_{n}} (γ_{n}) U_{γ_{o}} (γ_{o}) U_{k} (k) U_{r_{0}} (r_{0}) U_{v_{ϕ}} (v_{ϕ}) \times (2 π)^{- N / 2} \prod_{f} {(\frac{1}{γ_{o}} \frac{1}{k + f^{p}} | \tilde{h} (f; r_{0}, v_{ϕ}) |^{2} + \frac{1}{γ_{n}})}^{- 1 / 2} \times \exp [- \frac{1}{2} \frac{{| \tilde{i} (f) - \tilde{h} (f; r_{0}, v_{ϕ}) {\tilde{m}}_{o} (f) |}^{2}}{γ_{o}^{- 1} {(k + f^{p})}^{- 1} {| \tilde{h} (f; r_{0}, v_{ϕ}) |}^{2} + γ_{n}^{- 1}}],

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.

2.5.

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.²⁰ Given 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.²¹ It consists, for each iteration, in drawing samples under a chosen proposal distribution $q (θ | θ^{'})$ and accepting the samples (else, duplicating the previous value) with a prescribed probability $α$ . For the $k$ ’th iteration, $α$ writes

Eq. (16)

α = \frac{p (θ^{prop} | i)}{p (θ^{(k - 1)} | i)} \frac{q (θ^{(k - 1)} | θ^{prop})}{q (θ^{prop} | θ^{(k - 1)})},

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),²¹^–²³ 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

Eq. (17)

p (θ_{n} | i, θ_{m \neq n}) = \frac{p (θ_{n}) p (i | θ)}{p (i)},

where

θ_{n}

is the considered parameter and

θ_{m \neq 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.²¹

In 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} | θ_{n^{'}})$ around the current value of the parameter. The proposal distribution being symmetric, $q (θ_{n}^{(k - 1)} | θ_{n}^{prop}) = q (θ_{n}^{prop} | θ_{n}^{(k - 1)})$ thus $\frac{q (θ_{n}^{(k - 1)} | θ_{n}^{prop})}{q (θ_{n}^{prop} | θ_{n}^{(k - 1)})}$ 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.

Algorithm 1

Metropolis–Hastings-within-Gibbs algorithm.

Define initial

θ^{(0)}

for each iteration $k$ do

for each parameter $θ_{n}$ do

Propose

θ_{n}^{prop} \sim N (θ_{n}^{(k - 1)}, σ_{θ_{n}})

Acceptance rate

α \leftarrow \min (1, \frac{p (θ_{n}^{prop}) p (i | θ_{m < n}^{(k)}, θ_{n}^{prop}, θ_{m > n}^{(k - 1)})}{p (θ_{n}^{(k - 1)}) p (i | θ_{m < n}^{(k)}, θ_{m \geq n}^{(k - 1)})})

▹ Eqs. (16) and (17)

Random acceptance

u \sim U ([0; 1])

if $u < α$ then

Accept the proposal

θ_{n}^{(k)} \leftarrow θ_{n}^{prop}

else

Duplicate previous sample

θ_{n}^{(k)} \leftarrow θ_{n}^{(k - 1)}

end if

end for

3. Results on Simulated Astronomical Data

3.1.

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,²⁴ on a dark background of size $512 \times 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,⁵ 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 \times 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.

3.2.

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.

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.

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,⁵ 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).

Table 3

Mean value, associated standard deviation and true value, for γn, r0, vϕ, γo, and k for simulated astronomical image (stationary Gaussian noise), with p=3 and mo=mi.

Parameter	m±σ	True
$γ_{n}$ ( ${ph}^{- 2}$ )	$2.620 \times 10^{- 4} \pm 0.008 \times 10^{- 4}$	$2.621 \times 10^{- 4}$
$r_{0}$ (m)	$0.142 \pm 0.007$	0.15
$v_{ϕ} ({rad}^{2})$	$1.17 \pm 0.03$	1.30
$γ_{o}$ ( ${ph}^{- 2}$ )	$2.37 \times 10^{- 13} \pm 5.41 \times 10^{- 14}$	—
$K$	$0.768 \pm 0.594$	—

3.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.⁵ Thus, we notice that the errors on both parameters partially compensate. As a result, the normalized root mean squared error (RMSE) for the OTF, computed as $\sqrt{\frac{\sum | m_{\tilde{h}} - \tilde{h} |^{2}}{\sum | m_{\tilde{h}} |^{2}}}$ with true OTF $\tilde{h}$ and estimated OTF $m_{\tilde{h}}$ , is quite small (around 7%).

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).

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}} \pm 2 σ_{r_{0}}]$ , and the true $v_{ϕ}$ is in the interval $[m_{v_{ϕ}} \pm 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_{\tilde{h}}$ and standard deviation $σ_{\tilde{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_{\tilde{h}} \pm 2 σ_{\tilde{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(b), the object PSD model visually fits correctly the empirical object PSD, the slight overestimation being consistent with the slight underestimation of the OTF.

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.

3.4.

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)

Eq. (18)

\hat{o} = \underset{o}{\arg \max} p (i, o, γ_{n}, r_{0}, v_{ϕ}, γ_{o}, k) = \underset{o}{\arg \max} p (i | o, γ_{n}, r_{0}, v_{ϕ}) p (o | γ_{o}, k) = \underset{o}{\arg \max} - \ln (p (i | o, γ_{n}, r_{0}, v_{ϕ}) p (o | γ_{o}, k)) = \underset{o}{\arg \max} \sum_{f} (\frac{1}{2} {| \tilde{i} (f) - \tilde{h} (f) \tilde{o} (f) |}^{2} / S_{n} + \frac{1}{2} {| \tilde{o} (f) - {\tilde{m}}_{o} (f) |}^{2} / S_{o} (f)) .

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}$

Eq. (19)

\hat{o} (f) = \frac{{\tilde{h}}^{*} (f) \tilde{i} (f) + \frac{S_{n}}{S_{o} (f)} {\tilde{m}}_{o} (f)}{{| \tilde{h} (f) |}^{2} + \frac{S_{n}}{S_{o} (f)}} .

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.

Fig. 6

(a), (b) True object and image for simulated asteroid observation, $256 \times 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.

3.5.

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 Gaussian and along the axis suggesting that most parameters are not correlated a posteriori.

Fig. 7

Marginal posterior scatter graph of the samples for ( $r_{0}$ , $v_{ϕ}$ ) after boiling time.

Fig. 8

Marginal posterior scatter graph of the samples for ( $r_{0}$ , $γ_{o}$ ) after boiling time.

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.

3.6.

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.

Table 4

Summary of results on 10 noise realizations: true value, maximum error and bias (if available), standard deviation of estimates for γn, r0, vϕ, γo, and k, and minimum/maximum predicted uncertainty.

Parameter	True	Max. error	Empirical bias	Empirical std. dev.	Predicted uncertainty
$γ_{n}$ ( ${ph}^{- 2}$ )	$2.621 \times 10^{- 4}$	$2.3 \times 10^{- 6}$	$9.5 \times 10^{- 8}$	$9.0 \times 10^{- 7}$	$\in [7.8 \times 10^{- 7}, 8.0 \times 10^{- 7}]$
$r_{0}$ (m)	0.150	0.012	0.006	0.005	$\in [0.006, 0.009]$
$v_{ϕ} ({rad}^{2})$	1.30	0.14	0.12	0.01	$\in [0.02, 0.03]$
$γ_{o}$ ( ${ph}^{- 2}$ )	—	—	—	$4.1 \times 10^{- 14}$	$\in [4.5 \times 10^{- 14}, 8.9 \times 10^{- 14}]$
$K$	—	—	—	0.1	$\in [0.4, 0.7]$

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 small (third 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}} \pm 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_{\tilde{h}} \pm 2 σ_{\tilde{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.

Fig. 9

Results on OTF uncertainties for 10 realizations of a noisy simulated astronomical image: true OTF (in black) and predicted range $[m_{\tilde{h}} \pm 2 σ_{\tilde{h}}]$ (+ and − for upper and lower uncertainty bounds, each color corresponds to a noise realization).

3.7.

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⁵ 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_{ϕ}} \pm 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⁵ whereas $p$ corresponds to the slope of the object PSD 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.

Table 5

Mean value, associated standard deviation and true value, for γn, r0, vϕ, γo, and k for simulated astronomical image (stationary Gaussian noise), with p=2.9 and mo=mi.

Parameter	m±σ	True
$γ_{n}$ ( ${ph}^{- 2}$ )	$2.620 \times 10^{- 4} \pm 0.008 \times 10^{- 4}$	$2.621 \times 10^{- 4}$
$r_{0}$ (m)	0.141 ± 0.006	0.15
$v_{ϕ} ({rad}^{2})$	1.33 ± 0.02	1.30
$γ_{o}$ ( ${ph}^{- 2}$ )	$2.65 \times 10^{- 13} \pm 5.39 \times 10^{- 14}$	—
$K$	$0.619 \pm 0.469$	—

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.

Fig. 10

(a) Restored image, fixing $p = 3$ . (b) Restored image, fixing $p = 2.9$ . (c) Ten times the absolute difference between the two first images.

Fig. 11

Horizontal sectional plot of the restored images (at $N / 2$ ), fixing $p = 3$ and $p = 2.9$ , and their difference.

3.8.

Results 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.

Table 6

Mean value, associated standard deviation and true value, for γn, r0, vϕ, γo, and k for simulated astronomical image with a more realistic noise (Poisson + Gaussian noise), with p=3 and mo=mi.

Parameter	m±σ	True
$γ_{n}$ ( ${ph}^{- 2}$ )	$2.569 \times 10^{- 4} \pm 0.008 \times 10^{- 4}$	—
$r_{0}$ (m)	0.145 ± 0.007	0.15
$v_{ϕ} ({rad}^{2})$	1.20 ± 0.03	1.30
$γ_{o}$ ( ${ph}^{- 2}$ )	$2.49 \times 10^{- 13} \pm 5.57 \times 10^{- 14}$	—
$K$	0.717 ± 0.524	—

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,⁵ 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.

Table 7

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

Parameter	m±σ	True
$γ_{n}$ ( ${ph}^{- 2}$ )	$2.352 \times 10^{- 4} \pm 0.007 \times 10^{- 4}$	—
$r_{0}$ (m)	0.148 ± 0.008	0.15
$v_{ϕ} ({rad}^{2})$	1.18 ± 0.03	1.30
$γ_{o}$ ( ${ph}^{- 2}$ )	$2.86 \times 10^{- 13} \pm 7.32 \times 10^{- 14}$	—
$K$	0.522 ± 0.408	—

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.

4. Results on Simulated Satellite Image

4.1.

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 \times 512 pixels$ .²⁶ We simulate its observation using the ODISSEE AO system at OCA,¹² 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.

4.2.

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_{ϕ}$ .

Fig. 12

Left: Synthetic view of SPOT (true object), of size $512 \times 512$ . Right: simulated image, with true parameters $r_{0} = 0.10 m$ , $v_{ϕ} = 1.85 {rad}^{2}$ , and $γ_{n} = 1.00 10^{- 4} {ph}^{- 2}$ .

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.

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.

Table 8

Mean value, associated standard deviation and true value, for γn, r0, vϕ, γo, and k, for simulated satellite image (stationary Gaussian noise).

Parameter	m±σ	True
$γ_{n}$ ( ${ph}^{- 2}$ )	$1.02 \times 10^{- 4} \pm 3.63 \times 10^{- 7}$	$1.00 \times 10^{- 4}$
$r_{0}$ (m)	0.114 ± 0.008	0.10
$v_{ϕ} ({rad}^{2})$	2.16 ± 0.01	1.85
$γ_{o}$ ( ${ph}^{- 2}$ )	$1.90 \times 10^{- 14} \pm 2.37 \times 10^{- 15}$	—
$K$	3.14 ± 1.45	—

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}} \pm 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_{ϕ}$ ²⁵ 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_{\tilde{h}} \pm 3 σ_{\tilde{h}}]$ .

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).

Additionally, the 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.

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. 16

(a), (b) True object and image for simulated satellite observation, $256 \times 256$ cropped from $512 \times 512$ (Fig. 12). (c) Restored object using a L2-norm regularization, with positivity constraint, also cropped.

5. Results on Experimental Astronomical Data

After 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²⁷ 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.

Fig. 17

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

Results obtained for the PSF parameters (mean ± standard deviation) are the following: $r_{0} = 0.26 \pm 0.04 m$ and $v_{ϕ} = 2.62 \pm 0.06 {rad}^{2}$ . These values are close to the values obtained with AMIRAL⁵ ( $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}$ .¹³ We 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.

Fig. 18

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

6. Results 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.¹² We 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.

Fig. 19

(a) Envisat observed by ODISSEE at the OCA.¹² (b) Restored object using a L2-norm regularization, with positivity constraint.

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 ( $\approx 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.

Fig. 20

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

7. Conclusion

We have presented a marginal 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.²⁵ Moreover, jointly estimating $p$ with the other parameters is difficult as mentioned in earlier studies.⁵ In the near future, we plan to tackle the joint estimation of $p$ . Beyond the choice of $p$ , for objects that are far from isotropic, which is the case for some satellites, it would be worth considering an anisotropic prior model. To enable the estimation of such a richer prior model and improve the PSF estimation quality, we plan to add constraints on the object (namely, support and/or positivity contraints). Indeed, such constraints should help separate the contributions of the object and of the PSF to the image. This would then change the prior model on the object which can not be described by a simple Gaussian distribution anymore.

Finally, we currently sample each parameter individually, using a Metropolis–Hastings-within-Gibbs algorithm and the convergence speed issue has not been investigated. To accelerate the convergence, a possible development is to use a Metropolis–Hastings algorithm to sample all parameters jointly, thus without using the Gibbs sampler, and to use gradient-based methods, such as a Metropolis-adjusted Langevin algorithm.²¹^–²³

Code and Data Availability

The raw SPHERE astronomical data presented in this article are publicly available in Ref. 28, and the deconvolved images as well as the estimated PSF can be found in Ref. 29. The satellite data that support the findings of this article are not publicly available due to privacy concerns. They can be requested from the author at cyril.petit@onera.fr. Supporting code that can be used to generate PSFs according to the PSFAO19 model is publicly available in Ref. 30.

Acknowledgments

This study has been partly funded by the French Aerospace Lab (ONERA) in the framework of the SUSA Project. The satellite images were obtained using ODISSEE on the MéO telescope, thus the authors would like to thank OCA’s GEOAZUR team and the MéO operators, for their support and access to the telescope. Additionally, the authors would like to thank the reviewers of this paper for their careful review and insightful comments, which brought significant improvements to the paper. The authors declare no conflict of interest.

References

1.

J. Idier, Bayesian Approach to Inverse Problems, ISTE Ltd. and John Wiley & Sons Inc( (2008). Google Scholar

2.

L. M. Mugnier, T. Fusco and J.-M. Conan, “MISTRAL: a myopic edge-preserving image restoration method, with application to astronomical adaptive-optics-corrected long-exposure images,” JOSA, 21 1841 –1854 https://doi.org/10.1364/JOSAA.21.001841 JSDKD3 (2004). Google Scholar

3.

L. Blanco and L. M. Mugnier, “Marginal blind deconvolution of adaptive optics retinal images,” Opt. Express, 19 23227 –23239 https://doi.org/10.1364/OE.19.023227 OPEXFF 1094-4087 (2011). Google Scholar

4.

A. Levin et al., “Understanding and evaluating blind deconvolution algorithms,” in IEEE Conf. Comput. Vis. and Pattern Recognit., 1964 –1971 (2009). https://doi.org/10.1109/CVPR.2009.5206815 Google Scholar

5.

R. J.-L. Fétick et al., “Blind deconvolution in astronomy with adaptive optics: the parametric marginal approach,” MNRAS, 496 4209 –4220 https://doi.org/10.1093/mnras/staa1813 (2020). Google Scholar

6.

F. Orieux et al., “Estimating hyperparameters and instrument parameters in regularized inversion. Illustration for Herschel/SPIRE map making,” Astron. Astrophys., 549 A83 https://doi.org/10.1051/0004-6361/201219950 AAEJAF 0004-6361 (2013). Google Scholar

7.

E. Villeneuve and H. Carfantan, “Nonlinear deconvolution of hyperspectral data with MCMC for studying the kinematics of galaxies,” IEEE Trans. Image Process., 23 (10), 4322 –4335 https://doi.org/10.1109/TIP.2014.2343461 IIPRE4 1057-7149 (2014). Google Scholar

8.

F. Orieux, J.-F. Giovannelli and T. Rodet, “Bayesian estimation of regularization and point spread function parameters for Wiener-Hunt deconvolution,” JOSA, 27 1593 –1607 https://doi.org/10.1364/JOSAA.27.001593 JSDKD3 (2010). Google Scholar

9.

A. Yan et al., “Extending AMIRAL’s blind deconvolution of adaptive optics corrected images with Markov chain Monte Carlo methods,” Proc. SPIE, 12185 121853V https://doi.org/10.1117/12.2627414 PSISDG 0277-786X (2022). Google Scholar

10.

G. Demoment, “Image reconstruction and restoration: overview of common estimation structures and problems,” IEEE Trans. Acoust. Speech Signal Process., 37 (12), 2024 –2036 https://doi.org/10.1109/29.45551 IETABA 0096-3518 (1989). Google Scholar

11.

J.-L. Beuzit et al., “SPHERE: the exoplanet imager for the very large telescope,” Astron. Astrophys., 631 A155 https://doi.org/10.1051/0004-6361/201935251 AAEJAF 0004-6361 (2019). Google Scholar

12.

C. Petit et al., “LEO satellite imaging with adaptive optics and marginalized blind deconvolution,” in 21st AMOS Adv. Maui Opt. and Space Surveillance Technol. Conf., (2020). Google Scholar

13.

R. J.-L. Fétick et al., “Physics-based model of the adaptive-optics-corrected point spread function,” Astron. Astrophys., 628 A99 https://doi.org/10.1051/0004-6361/201935830 AAEJAF 0004-6361 (2019). Google Scholar

14.

F. Roddier, “The effects of atmospheric turbulence in optical astronomy,” Progr. Opt., 19 281 –376 https://doi.org/10.1016/S0079-6638(08)70204-X POPTAN 0079-6638 (1981). Google Scholar

15.

J.-M. Conan, “Etude de la correction partielle en optique adaptative,” Thèse de doctorat dirigée par Léna, Pierre Terre, océan, (1994). Google Scholar

16.

J.-M. Conan et al., “Myopic deconvolution of adaptive optics images by use of object and point-spread function power spectra,” Appl. Opt., 37 4614 –4622 https://doi.org/10.1364/AO.37.004614 APOPAI 0003-6935 (1998). Google Scholar

17.

S. Ramani et al., “Nonideal sampling and regularization theory,” IEEE Trans. Signal Process., 56 (3), 1055 –1070 https://doi.org/10.1109/TSP.2007.908997 ITPRED 1053-587X (2008). Google Scholar

18.

R. E. Kass and L. Wasserman, “The selection of prior distributions by formal rules,” J. Am. Stat. Assoc., 91 1343 –1370 https://doi.org/10.1080/01621459.1996.10477003 (1996). Google Scholar

19.

C. M. Bishop, Pattern Recognition and Machine Learning (Information Science and Statistics), 1st ed.Springer( (2007). Google Scholar

20.

H. Pishro-Nik, Introduction to Probability, Statistics, and Random Processes, Kappa Research, LLC( (2014). Google Scholar

21.

C. P. Robert and G. Casella, Monte Carlo Statistical Methods, 2 Springer( (2004). Google Scholar

22.

C. Vacar, J.-F. Giovannelli and Y. Berthoumieu, “Langevin and Hessian with Fisher approximation stochastic sampling for parameter estimation of structured covariance,” in IEEE Int. Conf. Acoust., Speech and Signal Process. (ICASSP), 3964 –3967 (2011). https://doi.org/10.1109/ICASSP.2011.5947220 Google Scholar

23.

C. Vacar, J.-F. Giovannelli and Y. Berthoumieu, “Bayesian texture classification from indirect observations using fast sampling,” IEEE Trans. Signal Process., 64 (1), 146 –159 https://doi.org/10.1109/TSP.2015.2480040 ITPRED 1053-587X (2016). Google Scholar

24.

L. Jorda et al., “OASIS: a simulator to prepare and interpret remote imaging of solar system bodies,” Proc. SPIE, 7533 753311 https://doi.org/10.1117/12.838893 PSISDG 0277-786X (2010). Google Scholar

25.

A. Yan et al., “Restauration d’images astronomiques corrigées par optique adaptative: méthode marginale étendue par algorithme MCMC,” in Proc. GRETSI 2022 Conf., (2022). Google Scholar

26.

L. M. Mugnier et al., “Myopic deconvolution from wave-front sensing,” J. Opt. Soc. Am. A, 18 862 –872 https://doi.org/10.1364/JOSAA.18.000862 JOAOD6 0740-3232 (2001). Google Scholar

27.

R. J.-L. Fétick et al., “Closing the gap between Earth-based and interplanetary mission observations: vesta seen by VLT/SPHERE,” Astron. Astrophys., 623 A6 https://doi.org/10.1051/0004-6361/201834749 AAEJAF 0004-6361 (2019). Google Scholar

28.

ESO, “Observational raw data query form,” http://archive.eso.org/eso/eso_archive_main.html (). Google Scholar

29.

LAM, “Asteroids as tracers of solar system formation: probing the interior of primordial main belt asteroids,” https://observations.lam.fr/astero/ (). Google Scholar

30.

R. Fétick, “Modelization of the adaptive optics PSF in PYthon (MAOPPY),” https://gitlab.lam.fr/lam-grd-public/maoppy (). Google Scholar

Biography

Alix Yan graduated from Ecole Nationale Supérieure d’Électronique, Informatique, Télécommunications, Mathématique et Mécanique de Bordeaux, Talence, France, in 2020. Currently, she is a third-year PhD student at ONERA, Châtillon, France. Her current research interests include high angular resolution imaging and image processing, more specifically Bayesian methods for image restoration. She is a member of SPIE.

Laurent M. Mugnier graduated from Ecole Polytechnique in 1988. He received his PhD in 1992 from Telecom Paris for his work on the reconstruction of incoherent-light holograms. He is currently a research director at ONERA, in charge of data processing activities of the High Angular Resolution Team in the Optics Department. His applications of interest include astronomical and satellite observation from the ground with adaptive optics, space optics, exoplanet detection, Earth observation, retinal imaging, and free-space optical communication.

Jean-François Giovannelli is currently a professor at the University of Bordeaux and a researcher at the Laboratoire de l’Integration du Matériau au Système, Groupe Signal-Image. His research focuses on inverse problems, mainly unsupervised and myopic questions. From a methodological standpoint, the developed regularization methods are both deterministic (penalty, constraints,…) and Bayesian. Regarding the numerical algorithms, the work relies on optimization and stochastic sampling. His application fields essentially concern astronomical, medical, proteomics, radars, and geophysical imaging.

Romain Fétick graduated from the SUPAERO Engineering School (Toulouse) in 2016. He then moved to Laboratoire d'Astrophysique de Marseille and ONERA for his PhD in post-processing of images taken with adaptive optics (AO). Since then, he works at ONERA on AO design and performance analysis, especially for ELT/HARMONI instrument. He contributed to mount and test the PAPYRUS AO instrument installed at the OHP observatory. He is interested in vulgarization, open access, and ethics thinking for science.

Cyril Petit is a senior research scientist at ONERA (France) with 20 years of experience in adaptive optics (AO). He graduated from ENSTA, then obtained his PhD in 2006 at ONERA, where he 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.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Alix Yan, Laurent M. Mugnier, Jean-François Giovannelli, Romain Fétick, and Cyril Petit "Marginalized myopic deconvolution of adaptive optics corrected images using Markov chain Monte Carlo methods," Journal of Astronomical Telescopes, Instruments, and Systems 9(4), 048004 (22 November 2023). https://doi.org/10.1117/1.JATIS.9.4.048004

Received: 22 March 2023; Accepted: 23 October 2023; Published: 22 November 2023

Access the abstract

JOURNAL ARTICLE
22 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

KEYWORDS

Point spread functions

Adaptive optics

Optical transfer functions

Satellites

Image deconvolution

Astronomy

Earth observing sensors

1.

Introduction

2.

Imaging Model and MMSE Estimator

2.1.

Imaging Equation

Eq. (1)

2.2.

PSF Model

Eq. (2)

Eq. (3)

Eq. (4)

Table 1

2.3.

Prior Distributions

Fig. 1

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Table 2

Fig. 2

2.4.

Marginal Estimator

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

2.5.

MMSE Estimator and Sampling

Eq. (16)

Eq. (17)

Algorithm 1

3.

Results on Simulated Astronomical Data

3.1.

Simulation Conditions

3.2.

Results on the Estimated Parameters and Derived Uncertainties

Fig. 3

Table 3

3.3.

Results on the OTF, on Object and Image PSD

Fig. 4

Fig. 5

3.4.

Results on the Restored Image

Eq. (18)

Eq. (19)

Fig. 6

3.5.

Posterior Coupling Between Parameters

Fig. 7

Fig. 8

3.6.

Test on Several Noise Realizations

Table 4

Fig. 9

3.7.

Impact of Hyperparameter p

Table 5

Fig. 10

Fig. 11

3.8.

Results with a More Realistic Noise

Table 6

Table 7

4.

Results on Simulated Satellite Image

4.1.

Simulation Conditions

4.2.

Results on the Estimated Parameters and Derived Uncertainties

Fig. 12

Fig. 13

Table 8

Fig. 14

Impact of Hyperparameter $p$