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.IntroductionGround-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.1 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,2 which leads to a degenerate solution3,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.5 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.6,7 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.8,9 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 Estimator2.1.Imaging EquationWe consider that the image results from the 2D discrete convolution of the object with the PSF , to which noise (mostly photon noise and detector readout noise) is added, giving the following imaging model10 This can also be written in the matrix form as , with the convolution matrix corresponding to the convolution of the object by . 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 instrument11 with the Zimpol imaging polarimeter5 and ground-based satellite observation at the Côte d’Azur Observatory (OCA) with the ODISSEE AO system.12 2.2.PSF ModelThroughout 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. Roddier14 shows that the long-exposure PSF re-writes as the convolution of three PSFs The first one is called the “static” PSF and corresponds to the static aberrations, the second one is the “detector” PSF describing the integration of the image on the detector’s pixels, and the last one is the “atmospheric” PSF corresponding to the impact of atmospheric turbulence. Conan15 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. can then be described by the phase PSD with the turbulence-induced phase, coming from the residual aberrations [not corrected by AO] in the pupil of the telescope, and , the phase variance, so that . The two main parameters of , thus of the PSFAO19 model, are the Fried parameter taken at the imaging wavelength (850 nm), describing the turbulence’s strength, and the variance of the residual turbulent phase , describing the quality of AO correction. Indeed, the residual phase PSD can be separated in two different spatial frequency zones, depending on the AO cutoff frequencyFor the corrected spatial frequencies, a Moffat model is used to describe the core of the PSD. The main parameter is the amplitude , which is very close to the residual phase variance: , with the AO-corrected area in the spatial frequency domain (for a circular AO-corrected area, ). 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 will be considered as known parameters, and their value will be fixed to the values in Table 1. Finally, is a normalization factor, which is used to normalize the integral of the Moffat function over the corrected spatial frequencies which requires that . For the high spatial frequencies, the theoretical Kolmogorov model of turbulence is used, the main parameter being the Fried parameter describing the turbulence’s strength, at the imaging wavelength. This model has been validated on several AO systems and on different telescopes.12,13Table 1Moffat fixed parameters.
2.3.Prior DistributionsNoise 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 . Therefore the noise covariance matrix is , with the identity matrix and the noise PSD is . 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 and its PSD . Given that we have little information on the the mean object , it is taken uniform on all pixels, estimated at the average value of the image considering that , supposing flux conservation. For the object PSD, we use the following parametric model: and the radial frequency. This circularly symmetric model is a slightly modified writing of Matérn’s model.16,17 In this model, sets the global PSD level, is the PSD decrease rate at high frequencies, and gives the breakpoint between the two regimes of the model. In previous works,5 attempts to estimate hyperparameter 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 is fixed to a standard value. In the case of astronomical observations of asteroids, a well-fitting empirical value is around , whereas for satellite observation a standard value for would be around 2.5 to 2.6.Let the image size in pixels. Given previous assumptions and the matrix form of the imaging model of Eq. (6) where we consider and as vectors and as a matrix, the likelihood writes In Eq. (1), we consider and 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 where denotes the discrete Fourier transform (DFT) and the product on 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 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. 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: 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 , the latter is diagonalized in the discrete Fourier domain, with as written in Eq. (5) on its diagonal, so that 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 , , , , and 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.18 The prior interval for each parameter is taken large enough. For , , and , 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 we take [5 cm; 30 cm] and for we take , 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 2Prior intervals and tuning of the Gaussian standard deviation for γn, r0, vϕ, γo, and k.
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 . Therefore , which means that the image and object PSD parameters and 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 . Moreover, as the hierarchical model reads, all parameters (, , , and ) are modeled as a priori independent. 2.4.Marginal EstimatorThe 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 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: 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: with image PSD and mean image[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 of Eq. (14) that best fit the empirical PSD .] The marginal posterior distribution for the parameters can then easily be written using Eqs. (12)–(14) with the uniform probability distribution for parameter , in the range defined for . In our case, these minimum and maximum values are given in Table 2.2.5.MMSE Estimator and SamplingThe 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.20 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.21 It consists, for each iteration, in drawing samples under a chosen proposal distribution and accepting the samples (else, duplicating the previous value) with a prescribed probability . For the ’th iteration, writes with 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–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 where is the considered parameter and the four other parameters. Note that 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.21 In our case, we use a random walk Metropolis–Hastings algorithm: the proposed sample for each parameter is drawn under a symmetric (Gaussian) distribution around the current value of the parameter. The proposal distribution being symmetric, thus in [in Eq. (16)] simplifies. The standard deviation of the Gaussian proposal distribution 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 1Metropolis–Hastings-within-Gibbs algorithm.
3.Results on Simulated Astronomical Data3.1.Simulation ConditionsThe 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 . True PSF parameters are and at the imaging wavelength , 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 (photons), typical from VLT/SPHERE/Zimpol asteroid observations (ESO Large Program ID 199.C-0074), therefore . The PSF and PSD parameters are estimated following the proposed method, except the mean object , which are estimated to the average value of the image, and the object PSD power, which is fixed to , which corresponds to a reasonable default value of 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 UncertaintiesIn Fig. 3, we plot the samples chains and the corresponding histograms for , , and . 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 , 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 and a 10% error on . Additionally, the estimated and are very close to the previous results obtained with AMIRAL: for similar conditions,5 the estimated PSF parameters were and (compared to and in Table 3). Table 3Mean 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.
3.3.Results on the OTF, on Object and Image PSDWe also compare the resulting OTF to the true OTF in Fig. 4. The slight underestimation of leads to the lowering of the global OTF level, the impact of which can mainly be seen at low frequencies. Concerning , its mild underestimation leads to a slower decrease of the OTF and impacts the slope of the latter at medium-high frequencies.5 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 with true OTF and estimated OTF , is quite small (around 7%). Concerning the uncertainties derived from our method, we notice in Table 3 that the true value for parameter is in the range , and the true is in the interval , therefore the uncertainties on PSF parameters seem under-estimated. We can also compute uncertainties directly on the sought OTF: for each sample , we compute the corresponding OTF to compute its sample mean and standard deviation , meaning the mean and standard deviation for each frequency of the OTF. As shown in Fig. 4, the true OTF is within the interval , 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 and ’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. 3.4.Results on the Restored ImageAfter 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 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. 3.5.Posterior Coupling Between ParametersSampling 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: (, ) and (, ). 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. The only pair of coupled parameters that does not have an elliptical-like scatter graph, but instead shows a strong a posteriori correlation, is and . We explain this correlation by the fact that as shown in Ref. 5, impacts the global level of the OTF whereas gives the global level of the object PSD. Therefore, given the expression of the image PSD in Eq. (14), both (, ) 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 RealizationsTo 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 4Summary 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.
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 ). Moreover, the predicted uncertainties for are close to the empirical average error made on the noise precision. For the PSF parameters, the error is smaller than 11%. Concerning parameter , the predicted uncertainty is very satisfactory: the true value is always within the interval . For parameter , 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 , which will be further discussed in 3.7. Finally, even though the uncertainties are under-estimated for with the default , concerning the OTF itself the uncertainties are always well estimated: for all 10 cases, the true OTF is within the interval , 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. 3.7.Impact of HyperparameterIn Ref. 25, we have tested our method in the exact same conditions, for another value for hyperparameter , tuned slightly differently, towards the “best” value5 in the supervised mode . (Both and 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 instead of 3.0, the error on becomes smaller and the estimated uncertainties are then satisfactory (the true is then in the interval ). Similarly to the correlation between and showed in 3.5, we interpret these results as another strong correlation between and the fixed hyperparameter , due to the similar impact they have on the image PSD. Indeed, impacts the slope of the OTF in medium-high frequencies5 whereas corresponds to the slope of the object PSD in medium-high frequencies. Therefore, both and tune the decrease of the image PSD in medium-high frequencies, which can explain their strong posterior correlation. Table 5Mean 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.
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. 3.8.Results with a More Realistic NoiseWe 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 , 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 6Mean 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.
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 , 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 and around 8% for . Their associated uncertainties are also still satisfactory for and, similarly, slightly underestimated for 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. Table 7Mean 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.
The results for all parameters do not change much, as previously the estimates for each parameter change less than a standard deviation , except for . 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 Image4.1.Simulation ConditionsWe 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 .26 We simulate its observation using the ODISSEE AO system at OCA,12 and with true PSF parameters and , at the imaging wavelength , 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 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 . The true object and the simulated data are shown in Fig. 12. The object PSD power is fixed to an empirical standard value for satellites , to fit the empirical object PSD. The Gibbs sampler is run for 100,000 iterations. 4.2.Results on the Estimated Parameters and Derived UncertaintiesIn Fig. 13, we plot the sample chains and the corresponding histograms for , , and . Similarly to the previous simulations, the sample mean values and standard deviations of the posterior distribution for each parameter are displayed in Table 8. The noise precision as well as PSF parameters and 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 and . 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 8Mean value, associated standard deviation and true value, for γn, r0, vϕ, γo, and k, for simulated satellite image (stationary Gaussian noise).
Concerning uncertainties, similarly to previous asteroid case, the posterior standard deviation for seems to be a good uncertainty prediction for this parameter as the true value is within the interval . On the contrary, 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 and 25 is still valid and further work on tuning hyper-parameter 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 . 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. 5.Results on Experimental Astronomical DataAfter 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 Vesta27 taken by SPHERE/Zimpol in the same mostly unsupervised mode as previously where , 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: and . These values are close to the values obtained with AMIRAL5 (, ) for the same conditions, the newly estimated being more likely than the one estimated by AMIRAL according to the known statistics on .13 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. 6.Results on an Experimental Satellite ImageFinally, 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.12 We fix paramater to a reasonable value for satellites () and run the Gibbs sampler for 500,000 iterations. Our results on the estimated PSF parameters are the following: and , which can be compared to the results obtained with AMIRAL: and , 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 (), 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. 7.ConclusionWe 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 , 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 , as we verified it by changing its value.25 Moreover, jointly estimating with the other parameters is difficult as mentioned in earlier studies.5 In the near future, we plan to tackle the joint estimation of . Beyond the choice of , 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.21–23 Code and Data AvailabilityThe 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. AcknowledgmentsThis 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. ReferencesJ. Idier, Bayesian Approach to Inverse Problems, ISTE Ltd. and John Wiley & Sons Inc(
(2008). Google Scholar
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
C. M. Bishop, Pattern Recognition and Machine Learning (Information Science and Statistics), 1st ed.Springer(
(2007). Google Scholar
H. Pishro-Nik, Introduction to Probability, Statistics, and Random Processes, Kappa Research, LLC(
(2014). Google Scholar
C. P. Robert and G. Casella, Monte Carlo Statistical Methods, 2 Springer(
(2004). Google Scholar
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
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
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
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
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
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
ESO, “Observational raw data query form,”
http://archive.eso.org/eso/eso_archive_main.html
().
Google Scholar
LAM, “Asteroids as tracers of solar system formation: probing the interior of primordial main belt asteroids,”
https://observations.lam.fr/astero/
().
Google Scholar
R. Fétick,
“Modelization of the adaptive optics PSF in PYthon (MAOPPY),”
https://gitlab.lam.fr/lam-grd-public/maoppy
().
Google Scholar
BiographyAlix 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. |