© The Authors 2023

1 Introduction

One of the main challenges for ground-based astronomical observations at visible wavelengths is to mitigate the effects of atmospheric turbulence. Without any special treatment, the angular resolution of the biggest telescopes on Earth is naturally limited to an equivalent 10–20 cm telescope in the visible due to atmospheric turbulence, leading to a resolution no smaller than 0.5–1 arcsecond on the astronomical targets (Fried 1966). Development in adaptive optics (AO) has improved the image quality by offering corrections to the perturbations of the atmosphere, partially restoring the angular resolution close to the diffraction limit of the biggest telescopes (e.g. Davies & Kasper 2012). Nevertheless, these AO corrections are only partial (Conan et al. 1994), and the final point spread function (PSF) can have a complex shape.

For instance, if the AO system were to perfectly correct all the aberrations within the range of its deformable mirror, the PSF would be the combination of the airy function near the optical axis and the remains of the extended seeing-limited wings for focal positions above the correction range. In reality, the AO system is not perfect, and suffers from measurement noise, temporal or aliasing errors, among others. These errors impact the PSF shape within the correction range and strongly depend on wavelength, field position, and time. Consequently, the low spatial frequencies (global shape) of the observed scene are conserved but the high spatial frequencies (the smallest details or sharp edges) are strongly attenuated or even lost, resulting in blurring of the images (Fusco et al. 2000). Precise image processing is therefore mandatory to retrieve these fine details in the observed objects, and deconvolution is a process designed to restore the high spatial frequencies to or near to their original level (Starck et al. 2002).

Deconvolution requires knowledge of the PSF, and the quality of the deconvolved image directly depends on the accuracy of its PSF model (Jolissaint et al. 2008; Davies & Kasper 2012). For example, Fétick et al. (2019a) and Marchis et al. (2021) used observations of asteroids from the European Southern Observatory (ESO) Large Program (ID 199.C- 0074, PI P. Vernazza, Vernazza et al. 2021) to show that the classical method of getting reference PSFs by observing calibration stars before or after the science observations leads to unacceptable errors in more than 50% of the cases. Among those cases, the deconvolution products resulted in strong artefacts at the asteroids’ edges, highlighted by the presence of a bright corona. This situation is clearly the consequence of inaccurate reference stars, differing from the actual PSF. Given that the strength and spatial structure of the atmospheric turbulence are constantly evolving, a PSF calibrator acquired with the same instrument setting before or after the science exposure can be completely uncorrelated to the actual PSF.

To overcome this issue, Fétick et al. (2020) introduced an approach whereby (i) the PSF is calculated using a simplified analytical model relying on a handful of parameters and (ii) these PSF parameters are estimated directly from the object thanks to a parametric marginal method. Once the PSF is estimated, it can then be used by standard deconvolution methods like Wiener filtering and Richardson–Lucy deconvolution (Richardson 1972; Lucy 1974), where the PSF is fixed.

In this paper, we aim to develop a more robust method for PSF estimation and deconvolution for astronomical images. In Sect. 2.2, we first investigate the extent to which we need to know the PSF parameters for proper deconvolution under different signal-to-noise ratios (S/Ns). Then, based on these requirements, we explore the performance of the PSF-parameter estimation directly from the image under different configurations in Sect. 3. In particular, we extend the work of Fétick et al. (2020) by considering different object shapes and AO-correction levels. We highlight the limitation of the previous method when dealing with elliptical objects, and we propose a generalisation of the process in Sect. 3.2. Finally, we illustrate the process with simulations and science observations from VLT-SPHERE in Sect. 4.

2 Requirement for deconvolution

2.1 Deconvolution method

During the deconvolution process, we introduced regularisation terms to prevent a drastic increase in noise. l¹ and l² regularisation are common in image processing. The l¹ (total variation) regularisation scales with the gradient of the object, which preserves the edge or the high-frequency contents but is insensitive to low-amplitude noise or features compared to the l² (Wiener like) approach. The l² penalises the criterion by scaling quadratically with the object’s gradient, and tends to over-smooth the edges or sharp features, as opposed to the l¹. We used MISTRAL (Mugnier et al. 2004) because it uses the l²–l¹ norm. This norm is adapted to restore objects with both smooth features and sharp edges. The criterion of the deconvolution is defined as $$J_{mistral}\left(o,h;\delta ,k\right)={||\frac{i-o⊛h}{{\sigma }_n}||}^2+{\delta }^2{\displaystyle {\displaystyle \sum }_{i,j}}\left|\frac{\nabla o_{i,j}}{\delta k}\right|-\ln \left(1+\left|\frac{\nabla o_{i,j}}{\delta k}\right|\right),$$(1)

where δ is the threshold between l²-like and l¹-like behaviour, κ is a scaling factor of the object’s gradient, σ_n is the noise, and ∇ is the isotropic gradient operator defined as, $$\nabla o_{i,j}=\sqrt{{\left(o_{i+1,j}-o_{i,j}\right)}^2-{\left(o_{i,j+1}-o_{i,j}\right)}^2}.$$(2)

The parameter δ acts as a threshold, where the criterion switches between a quadratic and the linear behaviour depending whether the gradient of the object is larger or smaller than δ. Therefore, the combination of l²–l¹ deconvolution is beneficial for high-angular-resolution objects with sharp edges because it cancels the penalties for the large gradients and ensures a good smoothing for the small gradients. We chose (κ = 1, δ = 1.5) based on visual inspections of the deconvolution with the true PSF and applied the same value for our tolerance study.

2.2 Deconvolution sensitivity

The sensitivity of the PSF parameters in deconvolution is often overlooked because there is no generic metric to describe the quality of the end product. Without access to the ground truth, it is difficult to derive a quantitative statement for the sensitivity of the results to the level of our knowledge of the PSF. While having a precise statement to evaluate the quality of the deconvolved object is almost impossible, here we aim to investigate the sensitivity of the process with respect to the input PSF parameters. This is achieved based on observation simulation, and the use of an analytical AO-PSF model. The advantage of the parametric PSF lies in the flexibility of its parameters, which can be adjusted to better estimate and quantify the impact of the PSF parameters on the final deconvolution accuracy.

Fétick et al. (2019a) provides the first quantitative estimate of the impact of our level of knowledge of the PSF and the physical meaning of the parameters for deconvolution of ground-based AO observation of an asteroid surface. The PSF model in Fétick et al. (2019a) was relatively simple, as they created synthetic PSFs using a Moffat profile. In particular, their PSF model does not include any physical parameters or the extended seeing halo. Nevertheless, these authors show that for PSFs with equivalent, large full width at half maximum (FWHM, or too low a Strehl ratio SR), the deconvolution induces the corona artefacts. In other words, overestimating the FWHM or underestimating the SR eventually leads to deconvolution artefacts. Conversely, for small FWHMs, or large SRs, the resulting images appear underdeconvolved, with the extreme situation being one in which a dirac-PSF is used and there is no deconvolution at all. Therefore, inaccurate or insufficient knowledge of the PSF impacts the deconvolution accuracy with different effects on the final image depending on the under- or overestimation of the PSF performance.

In this work, we pursue the analysis of Fétick et al. (2019a), with a more representative and accurate AO PSF model called PSFAO19 (Fétick et al. 2019b). We reiterate the analytical expression of the PSF, h, defined as $$W_{\varphi }\left(f\right)={\left[{\sigma }^2{M}_A\left(f_x,f_y\right)+C\right]}_{f<f_{\text{AO}}}+{\left[0.023r_0^{-5/3}{f}^{-11/3}\right]}_{f>f_{\text{AO}}},$$(3) $$h\left(f\right)={ℱ}^{-1}\left\{{\tilde{h}}_T{e}^{{\text{B}}_{\varphi }\left(0\right)}{e}^{ℱ-1\left\{W_{\varphi }\left(f\right)\right\}}\right\},$$(4)

where M_A is a Moffat-like function with fixed parameters, B_ϕ is the phase auto-correlation function, and W_ϕ(f) is the phase power spectrum density (PSD), which is characterised by physical parameters: the fried parameter (r₀) and phase residual variance (σ²). Beyond the AO cut-off frequency f_AO, the PSD is simply defined by a power law (Kolmogorov spectrum), the amplitude of which is set by the fried parameter (r₀). The AO cut-off frequency can also be related to the number of actuators used by the AO deformable mirror (DM). Indeed, the AO DM corrects over a limited frequency range. We therefore introduce a cutoff frequency of the DM such that $$f_{\text{AO}}=\frac{1}{2d_{\text{DM}}},$$(5)

where d_DM is the inter-actuator pitch projected into the pupil, often called the DM pitch, and is defined as $$d_{\text{DM}}=\frac{D}{\text{Nact}},$$(6)

where D is the telescope primary mirror diameter and Nact is the linear number of actuators projected onto the telescope primary mirror. Figure 1 shows different PSF profiles, as expected for an AO system installed on an 8m telescope working at visible wavelengths, when Nact changes from 20 to 40 actuators across the pupil (top figure), and when r₀ or σ² change for a given Nact.

We then took this PSF model with the parameters from Table 1 to simulate science observations. We simulated observations of one of the moons of Jupiter, Ganymede, with a reference PSF generated by PSFAO19 model and a high-resolution object taken from NASA (2021), and eventually binned down to fit the angular size expected from the ground, and the equivalent angular resolution for an 8m telescope working at visible wavelengths. Finally, we performed l²–l¹ deconvolution (Fusco et al. 2003; Mugnier et al. 2004) with 225 different PSFs – varying r₀ and σ² with respect to the true PSF parameters – in order to carry out the deconvolution tolerance study.

We use the rms error (RMSE) to evaluate the performance of the deconvolution with respect to PSF parameters. The RMSE is

defined as $${\text{RMSE}}_i=\sqrt{\frac{1}{n}{\displaystyle \sum {\left({\overline{o}}_i-{\overline{o}}_{\text{ref}}\right)}^2}},$$(7)

where n is the number of pixels, ō_i is the deconvolved object with 225 PSFs, and ō_ref is the deconvolved object with respect to the reference PSF.

Table 1 sets the reference PSF parameters for the simulations, with r₀ = 13.5 cm and σ² = 1.3 rad². We choose these PSF parameters as they are representative of the typical turbulence statistics and performance of the SPHERE extreme AO system (SAXO). Figure 2 shows the RMSE maps when the PSF parameters (r₀,σ²) used for deconvolution are changing. Each RMSE on the map is the sum over the pixel of Eq. (7). As expected, the RMSE between the real object and the deconvolved image is minimal when the exact same PSF is used for deconvolution. When the PSF parameters are changing, the deconvolution is degraded and the RMSE increases. We considered two cases of S/N, with a high S/N (=100) case (Fig. 2 left) and a low S/N (=10) case (Fig. 2 right), where the S/N is defined per pixel over the object. From Fig. 2, it is interesting to note that (1) the slope of the RMSE as a function of r₀ is steeper than σ² regardless of whether there is over- or underestimation of the parameters, and (2) the RMSE is asymmetric with respect to the PSF parameters. Overestimating r₀ and underestimating σ² will always result in a better deconvolution than the opposite scenario. Our results agree with the findings of Fétick et al. (2020), as increasing r₀ or reducing σ² corresponds to a better-quality PSF. This leads to underdeconvolution, whereas overdeconvolution results in strong artefacts, which increases the RMSE dramatically.

In order to quantify the effects described above, Fig. 2 (bottom right) shows the normalised RMSE values for different PSFs selected over the 225 sample, for which we also computed their SR. The conclusion remains as before; if the PSF is barely known, it is always more desirable to overestimate the PSF quality before using it for deconvolution. Figure 2 (bottom row) also shows that using the SR as a unique metric for quantifying the PSF quality is insufficient for deconvolution, as one can have many PSFs with the same SR but very different RMSE outcome. In other words, the PSF accuracy cannot be constrained with the SR-estimation accuracy alone, as the morphology of the PSF core (in our case constrained by σ²) also impacts the deconvolution. On the other hand, using the proposed analytical PSF model with constraints on σ² and r₀ seems appropriate for obtaining constraints on the required PSF accuracy for deconvolution. Finally, it is important to highlight that the overall tendency of the deconvolution sensitivity is independent of the S/N; therefore this latter conclusion is to be applied during the PSF estimation procedure, as presented below.

The final step in assessing deconvolution quality is a visual inspection, as normally there is no ground truth to compare with science observations. We inspect the impacts of over- and underestimating the PSF parameters on the deconvolution in Fig. 3 and their respective RMSE are labelled by diamonds in Fig. 2. Again with overestimated PSF parameters, we obtain better visual results in the deconvolved objects. The RMSE map (Fig. 2) and its visualisation (3) demonstrate that, in general, in order to achieve a good deconvolution, the estimated PSF parameters may fall in the ranges −15% < r₀ < 30% and −30% < σ² < 15%. Though these results do not define or set the generic requirement for deconvolution sensitivity, they provide quantitative guidelines as to the accuracy required on PSF parameters in order to achieve good-quality deconvolution.

Fig. 1 PSF profiles from PSFAO19 for different parameter choices. Top: PSFs for three different values of Nact. Bottom: PSF profiles with 40 by 40 actuators when the parameters r0 and σ2 are changed.

Fig. 2 Deconvolution sensitivity with two different S/Ns. Top left: Deconvolution RMSE with high S/N = 100. Top right: deconvolution RMSE with low S/N = 10. The colour intensity represents the strength of the RMSE. Redder colour in the RMSE map corresponds to lower RMSE values. The square represents the true PSF parameter and the diamond symbols represent the visual results displayed in Fig. 3. Bottom left: Strehl ratio error with the PSF parameters in the deconvolution RMSE map. Bottom right: Strehl ratio error versus normalised RMSE in the deconvolution RMSE map with a high S/N.

3 Direct PSF estimation from science observation

One of the traditional methods for blind deconvolution is joint-estimation (Markham & Conchello 1999; Mugnier et al. 2004), whereby simultaneous estimations of the PSF and the object of the image are performed. The difficulty with this method comes from the fact that we have twice as many unknowns (object and PSF) as data points, which is a very poor statistical contrast (Blanco & Mugnier 2011), and this can end up in coupling between the PSF and the object. Another approach, called marginal estimator Blanco & Mugnier (2011), is to separate the object o from the problem by integrating over the probability of the object distribution. This estimator is consistent and non-biased (Lehmann & Casella 2006) when there are more data or the S/N increases. Fétick et al. (2020) reduces the statistical contrast by PSF parmeterisation, $$P\left(\gamma |i\right)={\displaystyle \int P\left(i|o,\gamma \right)P\left(\gamma \right)P\left(o\right)\text{d}o,}$$(8)

where γ are the PSF parameters. To compute the marginalisation integral, we assume a Gaussian prior probability for the object and assume that the total noise is white (homogeneous overall), stationary, and follows Gaussian statistics.

AO images are dominated by photon noise and read-out noise. Mugnier et al. (2004) showed that the Poisson distribution may be approximated by Gaussian noise when the flux is typically > 10 photons/pixel. The read-out noise is Gaussian, and dominates the image when the S/N is low, and so in the end, one can state that the noise is Gaussian over a wide range of S/N. The stationary distribution is a strong assumption, as when photon noise dominates (for high-S/N cases), the spatial distribution of the noise will follow the object intensity. However, these assumptions are mandatory in order to write a simple analytical expression in the Fourier domain for minimisation, and we demonstrate in Sect. 3.1 that they do not introduce bias in the results. Along with the Gaussian statistics of the object, the probability function is $$P\left(i|\gamma \right)\propto \frac{1}{\sqrt{\det R_i}}{\exp }^{-\frac{1}{2}{\left(i-i_m\right)}^t{R}_i^{-1}\left(i-i_m\right)},$$(9)

where R_i is the image covariance matrix and i_m is the convolution of estimated PSF and the mean object. With the assumption of the Gaussian stationary noise, the image covariance matrix is $$R_i=H_{\gamma }{R}_o{H}_{\gamma }^t+\langle {\sigma }_n^2\rangle I_d,$$(10)

where $<{\sigma }_n^2>$ is the average noise variance of the pixels, I_d is the identity matrix, and H_γ is the convolution operator. The marginal criterion J can be written in the Fourier domain by taking the logarithm of P(γ|i) (Blanco & Mugnier 2011), $$\begin{array}{c}{J}_{\text{marg}}\left(\gamma ;S_{\text{obj}},<{\sigma }_n^2>\right)=\frac{1}{2}{{\displaystyle \sum }}^{\text{}}\ln \left(S_{\text{obj}}|{\tilde{h}}_{\gamma }{|}^2+<{\sigma }_n^2>\right)\\ +\frac{1}{2}{{\displaystyle \sum }}^{\text{}}\frac{{\left|\tilde{i}-\tilde{h}{\tilde{o}}_m\right|}^2}{\left(S_{\text{obj}}|\tilde{h}{|}^2+<{\sigma }_n^2>\right)}-\ln \left(\gamma \right),\end{array}$$(11)

where S_obj is the spatial power spectrum density (PSD) of the object and ln γ represents the prior information on the PSF parameters (if there is any). This method requires prior information on the object PSD, which can be provided via an axis-symmetric model (Conan et al. 1998), $$S_{\text{obj}}=\frac{k}{1+{\left(f/{\rho }_0\right)}^p},$$(12)

where k is the object PSD value at f = 0 (which is close to the square of the flux of the image), ρ₀ is a parameter related to the size of the object, and p is the index of the power law. This object PSD model is called ‘PSD Conan’ in the remainder of the paper. The validity of this approach was demonstrated by Fétick et al. (2020), who showed that a fully unsupervised algorithm (estimation of both the PSF and the PSD parameters of the object) could provide reasonable-quality PSFs for deconvolution of large objects (e.g. asteroid Vesta). The authors also show that the PSF estimation can be improved with mostly unsupervised approaches, where p is fixed, and all the remaining parameters are ‘free’ during the minimisation. We estimate the parameters by minimising Eq. (11) with the variable metric method (VMLM-B) descent algorithm developed by Thiebaut (2002). The PSF model and object PSD model with their respective analytical gradients are provided to the algorithm for the minimisation process.

Fig. 3 Simulated images produced from deconvolutions of the objects shown schematically in Fig. 2. The simulated image before deconvolution is shown in the top right corner of each panel of nine images. Left: deconvolved objects marked by the diamond symbols in the top left of Fig. 2 with high S/N (=100), right: deconvolved objects marked by the diamond symbols in the top right of Fig. 2 with low S/N (=10).

Fig. 4 Simulated observations of our study. Left: Simulated observation of Ganymede. Right: Simulated observation of the asteroid Vesta. These observations are shown with arrays of equal size, namely 300×300 pixels.

3.1 Limitation of the current method on PSF retrieval

Based on the case study of Vesta, we began including cases with less information for PSF estimation (i.e. smaller objects) and with higher-order AO corrections. Figure 4 shows two different objects of interest, Ganymede when it is almost 100% illuminated (case 1), and Vesta (case 2) with a larger angular size than Ganymede on the image. We use these two examples to explore the effect of different angular size on the performance of the method proposed by Fétick et al. (2020). We generated these cases with both Poisson (photon) noise and Gaussian (read-out) noise as expected in real data; this is different from the Gaussian noise assumption for the criterion Eq. (11). By doing so, we test whether the Gaussian noise assumption is valid within this marginal framework to estimate the PSF parameters. The telescope and the PSF configuration of the simulations are provided in Table 1. Before directly applying the object PSD model with the marginal estimator, we tested the robustness of the estimator by introducing the true object and noise into the model in Fig. 5. The supervised criterion map (or supervised map) means that all parameters are fixed, and we compute the value of the criteria for different PSF parameter values, and different pairs of r₀ and σ². Figure 5 shows that the shape of the marginal criterion is not convex, but it seems to exhibit only one minimum. A gradientdescent algorithm is therefore suitable for solving this criterion as long as it is fine-tuned to converge precisely in the sharp valley. In the case of Fig. 5, we substitute the object PSD model with the modulus of the Fourier transform of the object squared, S_obj = ‖F(obj)‖². In this ideal case, the minimum of the criteria map indeed corresponds to the actual PSF parameters. In other words, if we had perfect knowledge of the object PSD characteristics, we would be able to perfectly estimate the PSF from the image. Therefore, it is crucial to provide the estimator with an accurate object description to retrieve the PSF parameters for deconvolution. This example also shows that the assumptions made in terms of the noise statistics (Gaussian) and spatial distribution (stationary noise) are not biasing the criterion, and validates these required simplifications.

The next step is to reproduce the same analysis in Fétick et al. (2020), computing the supervised maps for all objects in Fig. 4 but using the axis-symmetric object PSD model from Eq. (12). The object PSD parameters being fitted empirically with the data.

The supervised maps computed for Ganymede (case 1) and Vesta (case 2) are shown in Fig. 6. Compared with Fig. 5, the estimated PSF parameters (minimum of the supervised map) can significantly differ from the actual PSF parameters. In the case of Vesta (first row of Fig. 6), we see that when the number of actuators increases, the estimated parameters slowly diverges away from the true solution. In this case, the situation is not dramatic, as the worst r₀ estimated would still fall within the required sensitivity for deconvolution as defined in Sect. 2.2. We then look at case 1 (Ganymede), which represents a more challenging configuration, with less illuminated pixels (and therefore less information) across the object than for Vesta. Results are shown in the second row of Fig. 6. We clearly see that the PSFparameter estimation fails, as the minimum of the criteria map notably deviates from the actual PSF parameters. At the lowest number of actuators (20 × 20), σ² is overestimated by 13%, which would still be acceptable for deconvolution purposes. However, when the AO-correction level is improved (Nact with a higher actuator density), the PSF-parameter estimation strongly diverges from the simulated input PSF-parameters. For both cases, it is important to remember that if one uses the real object PSD representation (and not the axis-symmetric approximation), then the PSF-parameter estimation is accurate, as demonstrated in Fig. 5. This is clear evidence that the source of error lies in the object PSD description, and that the axis-symmetric model assumption used in Fétick et al. (2020) does not provide a proper object representation for these cases.

We interpret our finding that the order of correction impacts the PSF-parameter estimation as follows: the higher the correction, the higher the spatial frequencies corrected, and therefore the higher the energy within the AO-corrected area (less energy is ‘lost’ in the halo). With the higher correction, the mismatch between the object PSD model and the object PSD is amplified, as more energy is translated into the object instead of the residual wavefront. A straight-line power law may no longer be sufficient to describe the features in the spatial frequencies domain. Similarly, when the angular size of the object is too small, only a few pixels are available from the image to estimate the PSF, and the process becomes very sensitive to the assumption made on the object PSD.

In conclusion, even though Fétick et al. (2020) further reduced the dimensionality of the problem from 2N² + to N² + the number of PSF parameters + the number of object PSD model parameters by introducing the analytical PSF model into the estimator, our study shows that with the current setup of the marginal estimator and the object PSD model (axis-symmetric), the PSF parameters cannot be estimated when (1) the object (or its angular size) is too small and (2) the order of AO correction is too high. Understanding these limitations, one may improve the object PSD model to ensure a better description for the marginal estimator.

Fig. 5 Supervised criterion map of Ganymede with 40×40 actuators, which is the worst case from Fig. 6. We compute the map with ‖F(obj)‖2 and noise from the simulation.

3.2 Improved PSD model for the object

In Sect. 3.1, we show that the description of the object PSD prevents the marginal estimator from providing a more realistic PSF evaluation, as the estimated minimum diverges from the simulated parameters. The estimator converges towards the true PSF parameters when the ground truth of the object and the noise are provided, as opposed to the case of Ganymede in Fig. 6 with 40 × 40 actuators. A superior object PSD model should increase the adaptability of the estimator, for instance, allowing higherorder AO correction and fewer pixels across the object.

Figure 7 shows the PSD circular average for Ganymede. From Fig. 7, one can already see that the axis-symmetric object PSD model fails to describe the ring-like characteristics, as well as some of low-spatial-frequency behaviour, which might explain the breakdown of the marginal estimator in Fig. 6. Therefore, we propose the introduction of an elliptical function to describe the object PSD, defined as $$\begin{array}{c}{\rho }^2=\left(a^2{\sin }^2\theta +{\left(r\times a\right)}^2{\cos }^2\theta \right)f_x^2\\ +2\left({\left(ra\right)}^2-a^2\right)\sin \theta \cos \theta f_x{f}_y\\ +\left(a^2{\cos }^2\theta +{\left(r\times a\right)}^2{\sin }^2\theta \right)f_y^2,\end{array}$$(13)

where r is ratio between semi-major and semi-minor axes, a is the semi-major axis, and r × a is the semi-minor axis. We choose an elliptical model as it is more general than a simple disk, and will better represent observations of extended sources in the Solar System. The final elliptical object PSD model can eventually be described by $$S_{\text{obj}}\left(\rho \right)=k\times {\left|\frac{J_1\left(2\pi \rho \right)}{2\pi \rho }\right|}^2,$$(14)

where k is the object PSD value at f = 0, ρ is the spatial frequency, J₁ is the Bessel function of the first kind, and S_obj (a, r, θ,µ) is the object PSD model in the frequency domain. We adopt the same convention for k as Blanco & Mugnier (2011), and define $\mu =\langle {\sigma }_n^2\rangle /k$. In practice, we find it necessary to add a Gaussian filter to the object PSD, and so the elliptical object PSD model is $$S_{\text{obj}}\left(\rho \right)=k\times {\left|\frac{J_1\left(2\pi \rho \right)}{2\pi \rho }\right|}^2⊛𝒢,$$(15)

where 𝒢 is a Gaussian function. The need for the Gaussian filter comes from the sharpness of the airy-like pattern. Sharp rings lead to some pixel values close to zero, driving the marginal criterion in an aberrant fashion. We find that a fixed Gaussian filter is adequate for all objects, and so this term remains constant. The Gaussian filter 𝒢 we use is $$𝒢=\frac{1}{2\pi {\sigma }^2}{e}^{-\left(x^2+y^2\right)/\left(2{\sigma }^2\right)},$$(16)

where x and y are the coordinate of the array, and σ is the standard deviation of the distribution, the value of which is fixed to 2 for the filter.

With this particular formalism, and unlike the model in Conan et al. (1998), the physical information for the object is compressed into the variables, with a (semi-major axis) and r (ratio between semi-major and semi-minor axes) being directly related to the size of the object. Here, θ relates to the orientation of the object PSD and µ is associated with the noise level of the image. Considering that this model is analytical, minimising the model parameters with their analytical gradients is beneficial to our multi-variable PSF estimation. The object PSD parameters are a secondary product of the marginal criterion because the goal is to retrieve the PSF parameters to build a PSF for deconvolution.

Figure 7 shows how this elliptical model improves the representation of the PSD features. The object PSD parameters are derived by fitting ‖F(obj)‖² with the elliptical model from Eq. (15). The next step is to validate the performance of this PSD model for PSF-parameter estimation with simulations.

Fig. 6 Supervised criterion map using the axis-symmetric and elliptical object PSD model. Top: Vesta with the Conan model (S/N = 100) [20 by 20, 30 by 30, 40 by 40], Middle: Ganymede with the Conan model (S/N = 100) [20 by 20, 30 by 30, 40 by 40], Bottom: Ganymede with the elliptical model (S/N = 100) [20 by 20, 30 by 30, 40 by 40].

Fig. 7 Circular average of the normalised intensity versus spatial frequency of ‖F(obj)‖2 (red solid line), the elliptical object model by Eq. (15) (blue dotted-dash line), and the Conan model (green dashed line).

4 Applications: Simulation and on-sky data

4.1 Simulation: Fully supervised case

We first test the marginal estimator with the elliptical PSD model by reproducing the analysis of Fig. 6 for the most problematic case of Ganymede (small object) with the highest AO correction (bottom-right plot of Fig. 6). Compared with the bottom-right plot of Fig. 6, introducing the elliptical PSD model into the marginal estimator leads to a significant improvement, as shown in Fig. 6. The minimum of the criteria is now within less than 10% of the actual PSF parameters. Again, this result confirms that the assumptions made in the noise model do not introduce bias into the estimation process.

For the case of Ganymede, the relevant parameters are found to be as follows: $$\{\begin{array}{l}{r}_0=14.2\text{\hspace{0.17em}cm}\hfill \\ {\sigma }^2=1.46{\text{red}}^2\hfill \\ \mu =9.037\times {10}^{-9}\hfill \\ a=34.0\text{\hspace{0.17em}pix}\hfill \\ r=0.867\hfill \\ \theta =-1.48\text{\hspace{0.17em}rad.}\hfill \end{array}$$(17)

When the object parameters (µ, a, r, θ) are well calibrated, our elliptical model provides a stronger constraint towards the estimation of (r₀, σ²), and therefore allows a degenerate regime to be avoided, as opposed to the axis-symmetric object PSD model shown in the bottom-right plot of Fig. 6.

4.2 Simulation: Fully unsupervised case

Without knowing the identity of a target object, one can perform the marginal estimator along with object parameters. This is what we refer to as ‘fully unsupervised’, when the optimisation process is carried out on all free parameters, estimating both the object PSD parameters and the PSF at the same time. Fétick et al. (2020) showed that such an approach does not provide satisfactory results, mostly because the criterion is flat and because of the couplings between parameters. We reproduced the experiment performed by these latter authors using the elliptical PSD model, but the same limitations are faced, especially with a strong coupling between r₀ and µ.

In the Fourier domain, the PSD of the image is $S{}_{\text{obj}}\times {\left|\tilde{h}\right|}^2+n$, and both r₀ and µ determine the intensity of the image PSD. The estimator therefore has difficulty in separating the contribution of the PSF from the object. This means a stronger constraint is required either on the PSF parameter (r₀) or on the object parameter (µ). In particular, we note that, if we can provide initial guesses that are sufficiently close to the true value (−5% to +5%), then the minimisation can converge towards a good solution. The fully unsupervised case therefore does not appear to be a robust PSF-estimation method.

4.3 Simulation: mostly unsupervised case

For actual observations, we do have access to some information for the object from the observation itself. For instance, our object PSD model allows us to obtain relevant guesses for S_obj via (a, r, θ, µ) directly from the Fourier domain of the image, which was not possible from the axis-symmetric object PSD model. Indeed, these parameters are linked to the physical parameters, which directly come from the size and orientation of the object in the image. Furthermore, µ can be estimated directly from the circular average of the image PSD by calculating the ratio between the maximum and the minimum of this circular average. Therefore, the PSF parameter estimation can be carried out with µ fixed to the guess obtained from the science image, and allowing (a, r, θ) free to vary, starting from educated guesses. Applying such a method, Fig. 8 shows the resulting criteria maps for the same case as shown in Fig. 6, with the reference r₀ = 13.0 cm and σ² = 1.3 rad². The relevant parameters are found to be, $$\{\begin{array}{l}{r}_0=14.9\text{\hspace{0.17em}cm}\hfill \\ {\sigma }^2=1.47{\text{red}}^2\hfill \\ \mu =9.03\times {10}^{-9}\hfill \\ a=34\text{\hspace{0.17em}pix}\hfill \\ r=0.866\hfill \\ \theta =-1.48\text{\hspace{0.17em}rad.}\hfill \end{array}$$(18)

In conclusion here, fixing µ significantly improves the PSF estimation because it removes its coupling with r₀. With the initial guesses of the object parameters obtained from the observation, we are able to retrieve a solution close to the fully supervised case (known object information). This approach provides a robust way of estimating the PSF parameters from the images, even in specific configurations of small and elongated objects. The final step is to validate the approach with real data, as shown in the following section.

Fig. 8 Mostly unsupervised criterion map of Ganymede 40 × 40 using the elliptical PSD model in Eq. (15).

4.4 On-sky: Fully or mostly unsupervised

Here, we make use of actual observations to test and validate the newly proposed algorithm. To test whether or not our marginal estimator with the elliptical object PSD model provides robust estimates of the PSF parameters, we select a test case, namely the observation of the asteroid Kleopatra from 14 July 2017, because its elliptical shape provides an appropriate challenge to our estimator. The asteroid Kleopatra was observed as part of the ESO Large Program (ID 199.C- 0074, PI P. Vernazza – Vernazza et al. 2021), targeting Main Belt asteroids to study their albedo, shape, excavation volume, and surface. This ESO Larger Program used the Zurich IMaging POLarimeter (ZIMPOL) instrument with the N_R filter of the central wavelength λ = 645.9 nm. ZIMPOL is mounted on Spectro Polarimetric High-contrast Ex-oplanet REsearcher (SPHERE) of the Very Large Telescope (VLT) observatory. It is equipped with SAXO, a high-order AO system. The images were previously deconvolved by Marchis et al. (2021). In this case, we have no prior information on the object, and we apply the strategy defined in Sect. 4.3: fix µ and extract first guesses for (a, r, θ) directly from the image. We then perform PSF estimation and use this to deconvolve the image.

For this observation, the PSF parameters found are $$\{\begin{array}{l}{r}_0=12.1\text{\hspace{0.17em}cm}\hfill \\ {\sigma }^2=4.17{\text{red}}^2\hfill \\ \mu =2.42\times {10}^{-9}\hfill \\ a=29\text{\hspace{0.17em}pix}\hfill \\ r=0.55\hfill \\ \theta =-0.92\text{\hspace{0.17em}rad.}\hfill \end{array}$$(19)

We compare the result of our algorithm to other deconvolution strategies in Fig. 9, that is, with: (a) the observation, (b) a deconvolution performed using the PSF estimated with the previous model Conan et al. (1998), (c) a deconvolution using the estimation with the elliptical model, (d) a deconvolution using the calibration star observed after the observation, and (e) the deconvolution by Marchis et al. (2021), who used a Moffat PSF profile. Regarding the deconvolution, we use a l¹ – l² deconvolution scheme Fusco et al. (2003) while keeping all hyperparameters consistent for all the different deconvolved objects. Figure 10 shows the intensity along the longest axis of the image, the deconvolution of the axis-symmetric model, the deconvolution of the elliptical model, and the deconvolution of the previous work. The estimation made by the elliptical model successfully retrieves the morphology of the asteroid while avoiding overdeconvolution. Finally, Fig. 11 shows the PSF and the OTF profile of the estimated PSF with the elliptical model (red solid line), a reference PSF (black dotted line) observed after the asteroid observation, and the estimated PSF with the axis-symmetric model. We can retrieve the overall PSF shape with the elliptical model when comparing with the reference PSF, assuming it provides a good proxy for the PSF.

Fig. 9 Observation and the deconvolutions with different PSFs. (a) Observation, (b) Deconvolution performed using the PSF estimated using a previous model Conan et al. (1998). (c) Deconvolution using the estimation with the elliptical model, (d) Deconvolution using the calibration star observed after the observation, (e) Deconvolution by Marchis et al. (2021) using a Moffat PSF profile. These images are plotted with the same intensity scale.

Fig. 10 Intensity of the image and the deconvolutions. Slice of the asteroid Kleopatra for the image (black line), previous work related to Marchis et al. (2021) (green dotted line), a deconvolution using the marginal estimator with the axis-symmetric model Conan et al. (1998) (blue dotted line), and a deconvolution using the marginal estimator with elliptical model (red solid line).

Fig. 11 PSFs and OTF profiles used for deconvolution. Top: cut of the PSFs and the reference PSF. The PSF estimated with by Conan et al. (1998) is shown with a blue dashed line, the PSF estimated with the elliptical model is shown with a red solid line, and the reference PSF obtained after the observation is shown as a black dotted line. Bottom: cut of the OTFs and the reference OTF. The OTF estimated by Conan et al. (1998) is shown as a blue dashed line, the OTF estimated with the elliptical model is shown as a red solid line, and the reference OTF obtained after the observation is shown as a black dotted line.

5 Conclusion and perspectives

This paper addresses the deployment of advanced data-processing techniques in the framework of AO-assisted observations of extended targets and Solar System objects in particular. AO correction is not perfect, especially at visible wavelengths, and image post-processing is mandatory to retrieve the finer details in the observed objects. Deconvolution is the image postprocessing step to reduce the residual blurring effect on the image, but it requires an accurate description of the actual AO-PSF. In the first part of this paper, we make use of an analytical AO-PSF model to investigate the sensitivity of the deconvolution process to the accuracy of the available knowledge of the PSF. We defined boundaries on the PSF-parameter values in order to achieve a deconvolution without inducing strong artefacts. Our results show that if the estimated PSF parameters fall in the ranges of −15% < r₀ < 30% and −30% < σ² < 15%, then we can achieve deconvolutions without introducing many artefacts and distortions. The boundaries are not symmetrically distributed, and we find that it is always better to overestimate the PSF quality in order to reduce the error in the deconvolution, which agrees with the findings of Fétick et al. (2020). We also conclude that the SR alone is not a good metric for the PSF description for deconvolution.

We then explored algorithms allowing to extract the PSF directly from the science images. Such methods, known as marginal blind-deconvolution, have previously been shown to be efficient, especially for Solar System-like observations (Fétick et al. 2020). In Sect. 3.1, we illustrate the limitations of such a method for specific cases, and in particular for objects of small angular size. We then improved the object description by introducing a new elliptical object PSD model in Sect. 3.2 and we present the results of our tests of this elliptical model in Sect. 4, which use the problematic cases from Sect. 3.1. After testing with simulations, we validated our model with on-sky data from VLT-ZIMPOL. With the elliptical object PSD model, we successfully retrieve the PSF from the image and improve the shape and flux retrieval on the observations.

Other than improving the PSF estimation with the challenging cases for the previous model, more importantly, our model allows us to obtain relevant guesses directly from the Fourier domain of the image. The object parameters required are linked to physical parameters, that can be estimated from the size and orientation of the object in the image. As such, the process is more robust as it starts from educated guesses.

In future work, a first path to improve the PSF estimation would be to use external data available from AO telemetry in order to introduce priors on the PSF parameters. AO telemetry data provide insights into observing conditions and AO performance, from which we can obtain estimations of r₀ and/or σ². These priors may be used to improve the minimisation process even further, and eventually the robustness of the PSF estimation. Another path would be to extend the method to 3D data cubes, for instance from Integral Field Spectrograph observations. Assuming that the structure of an object does not change significantly with wavelength, a single set of object PSDs could be used for all wavelengths, which would further improve the estimation algorithm. Furthermore, one could also take advantage of the known PSF scaling with wavelengths, reducing the number of unknowns with respect to the data even further, which would constrain the output PSF. Finally, further work is needed to adapt our method to different classes of extended objects; for example, galaxies with different topologies. Extending the blind deconvolution method to handle a more diverse class of objects will either need additional inputs from the astronomical model, or machine learning methods could perhaps be used to build more complex models. Access to the Python package of this work is possible upon request.

Acknowledgements

This work benefited from the support of the French National Research Agency (ANR) with WOLF (ANR-18-CE31-0018), APPLY (ANR-19-CE31-0011) and LabEx FOCUS (ANR-11-LABX-0013); the Programme Investissement Avenir F-CELT (ANR-21-ESRE-0008), the Action Spécifique Haute Résolution Angulaire (ASHRA) of CNRS/INSU co-funded by CNES, the ECOS-CONYCIT France-Chile cooperation (C20E02), the ORP H2020 Framework Programme of the European Commission’s (Grant number 101004719) and STIC AmSud (21-STIC-09).

References

All Tables

All Figures

	Fig. 1 PSF profiles from PSFAO19 for different parameter choices. Top: PSFs for three different values of Nact. Bottom: PSF profiles with 40 by 40 actuators when the parameters r0 and σ2 are changed.
In the text

Fig. 2 Deconvolution sensitivity with two different S/Ns. Top left: Deconvolution RMSE with high S/N = 100. Top right: deconvolution RMSE with low S/N = 10. The colour intensity represents the strength of the RMSE. Redder colour in the RMSE map corresponds to lower RMSE values. The square represents the true PSF parameter and the diamond symbols represent the visual results displayed in Fig. 3. Bottom left: Strehl ratio error with the PSF parameters in the deconvolution RMSE map. Bottom right: Strehl ratio error versus normalised RMSE in the deconvolution RMSE map with a high S/N.

In the text

Fig. 3 Simulated images produced from deconvolutions of the objects shown schematically in Fig. 2. The simulated image before deconvolution is shown in the top right corner of each panel of nine images. Left: deconvolved objects marked by the diamond symbols in the top left of Fig. 2 with high S/N (=100), right: deconvolved objects marked by the diamond symbols in the top right of Fig. 2 with low S/N (=10).

In the text

	Fig. 4 Simulated observations of our study. Left: Simulated observation of Ganymede. Right: Simulated observation of the asteroid Vesta. These observations are shown with arrays of equal size, namely 300×300 pixels.
In the text

	Fig. 5 Supervised criterion map of Ganymede with 40×40 actuators, which is the worst case from Fig. 6. We compute the map with ‖F(obj)‖2 and noise from the simulation.
In the text

	Fig. 6 Supervised criterion map using the axis-symmetric and elliptical object PSD model. Top: Vesta with the Conan model (S/N = 100) [20 by 20, 30 by 30, 40 by 40], Middle: Ganymede with the Conan model (S/N = 100) [20 by 20, 30 by 30, 40 by 40], Bottom: Ganymede with the elliptical model (S/N = 100) [20 by 20, 30 by 30, 40 by 40].
In the text

	Fig. 7 Circular average of the normalised intensity versus spatial frequency of ‖F(obj)‖2 (red solid line), the elliptical object model by Eq. (15) (blue dotted-dash line), and the Conan model (green dashed line).
In the text

	Fig. 8 Mostly unsupervised criterion map of Ganymede 40 × 40 using the elliptical PSD model in Eq. (15).
In the text

Fig. 9 Observation and the deconvolutions with different PSFs. (a) Observation, (b) Deconvolution performed using the PSF estimated using a previous model Conan et al. (1998). (c) Deconvolution using the estimation with the elliptical model, (d) Deconvolution using the calibration star observed after the observation, (e) Deconvolution by Marchis et al. (2021) using a Moffat PSF profile. These images are plotted with the same intensity scale.

In the text

	Fig. 10 Intensity of the image and the deconvolutions. Slice of the asteroid Kleopatra for the image (black line), previous work related to Marchis et al. (2021) (green dotted line), a deconvolution using the marginal estimator with the axis-symmetric model Conan et al. (1998) (blue dotted line), and a deconvolution using the marginal estimator with elliptical model (red solid line).
In the text

Fig. 11 PSFs and OTF profiles used for deconvolution. Top: cut of the PSFs and the reference PSF. The PSF estimated with by Conan et al. (1998) is shown with a blue dashed line, the PSF estimated with the elliptical model is shown with a red solid line, and the reference PSF obtained after the observation is shown as a black dotted line. Bottom: cut of the OTFs and the reference OTF. The OTF estimated by Conan et al. (1998) is shown as a blue dashed line, the OTF estimated with the elliptical model is shown as a red solid line, and the reference OTF obtained after the observation is shown as a black dotted line.

In the text