© The Authors 2024

1 Introduction

The performance of ground-based instruments is limited by the atmospheric turbulence that corrugates the incident wavefront of the observed target (Roddier 1981). For short exposures, the diffraction-limited point spread function (PSF) of the telescope breaks down to a random speckle field. For long exposures, their average is equivalent to the PSF of a 10 to 20 cm telescope (Fried 1966), strongly degrading the sensitivity and the resolution.

Introduced in the 1990s, adaptive optics (AO) systems are now commonly deployed in instruments of ground-based observatories to mitigate the effects of the turbulence (Tyson 2015).

The coupling of (i) PSF prediction to optimise the instrument design and performance (Fusco et al. 2006; Dohlen et al. 2016), with (ii) PSF modelling (Jolissaint et al. 2006; Fétick et al. 2019b; Berdeu et al. 2023) to improve the data postprocessing via model-fitting or deconvolution (Beltramo-Martin et al. 2020), has resulted in a leap forward in numerous fields such as photometry and astrometry of stellar populations (Turri et al. 2017; Monty et al. 2018), spectroscopy and kinematics of distant galaxies (Förster Schreiber et al. 2018; Bianchin et al. 2022), or the study of extended object surface topology (Rimmele et al. 2021; Vernazza et al. 2021).

In the last decade, the arrival of extreme AO systems pushed the performance in high-contrast and high-resolution imaging even further (Jovanovic et al. 2015). In this context, PSF prediction and modelling reach their limits, unable to faithfully reproduce the complexity of observed AO-corrected PSFs (AO-PSFs). Point spread function estimation or reconstruction techniques must then be considered (Beltramo-Martin et al. 2020). But direct estimation of the PSF parameters from the AO telemetry (Véran et al. 1997; Clénet et al. 2008) or reference PSFs obtained on calibration sources (internal or natural stars) before or after the observation (Mugnier et al. 2004) are not always sufficient, and can for example lead to strong deconvo-lution artefacts (Marchis et al. 2006; Fétick et al. 2019a, 2020; Lau et al. 2023).

Among others, the challenge arises from the evolving nature of AO-PSFs. By essence, the turbulence is random and the AO loops are not perfect, suffering from measurement noise, or temporal and aliasing errors (Rigaut et al. 1998). Finally, flexions or thermal dilations in the instruments induce quasi-static speckles that slowly drift with time, changing the PSF’s small and faint structures (Milli et al. 2016; Vigan et al. 2019). As a consequence, a given observation represents one realisation of an AO-PSF that is barely reproducible, or even impossible to reproduce.

To overcome this random nature of the AO-PSF, the only solution is to extract and reconstruct it directly from the data of interest, a problem known as blind deconvolution (Stockham et al. 1975; Thiébaut & Conan 1995; Soulez et al. 2012; Fétick et al. 2020). The challenge is to overcome the degeneracy between the estimation of the PSF and the object (Little & Rubin 1983; Blanco & Mugnier 2011) and to solve an ill-posed problem with potentially twice more unknowns than measurements if both the images of the PSF and the object are to be retrieved from a single data image. As a consequence, no method is strictly ‘blind’, and there is a need to make some assumptions about the PSF or the deconvolved object (or both) to avoid local minima of the problem and converge on the best ‘PSF and object’ pair.

One solution is to implement marginal approaches (Blanc et al. 2003; Beltramo-Martin et al. 2020). The problem is rewritten to split the PSF and the object contributions, in order to integrate over the space of the object parameters to marginalise them out of the problem. The deconvolution of the object is only performed as a second step. Marginalisation has proved to be efficient when combined with parametric PSF models that strongly limit the number of unknowns to fit (Fétick et al. 2019b). These methods nonetheless imply a need to know some prior on the object, such as the structure of its power spectrum density (Fétick et al. 2020; Yan et al. 2023) or its overall shape (Lau et al. 2023). In addition, such PSF models, relying on a limited number of parameters, cannot fully grasp the complexity of real and potentially broadband AO-PSFs of high-contrast and high-resolution imagers. The challenge will be even greater with the incoming extremely large telescopes (ELTs). Their segmented primary mirror and large spiders holding the secondary mirror will produce intricate and structured PSFs that are very sensitive to AO residuals (Hippler et al. 2019; Neichel et al. 2020; Simioni et al. 2020; Hedglen et al. 2022). More versatile methods than parametric approaches are consequently needed.

Finally, in the context of asteroid study, beyond the deconvo-lution of the images of asteroids’ surfaces (Vernazza et al. 2021), the approximate knowledge of the PSF also limits the study of asteroids’ close vicinity and the detection of faint moons. Similar constraints apply when observing the moon systems around the giant planets of our Solar System (Showalter & Lissauer 2006; Assafin et al. 2008). To get rid of the bright halo induced by the PSF extensions around the main object, techniques were adapted from exoplanet detection algorithms, based on local averaging or median filters (Marchis et al. 2006; Assafin et al. 2008; Pajuelo et al. 2018). But the problem to solve is quite different, the halo being produced not by a coherent source (an unresolved star) but by an extended object (the resolved asteroid or planet). These techniques consequently do not account for the physics at the origin of the halo that is poorly estimated. They thus inherently suffer from self-subtraction problems and can bias the moon photometry or even prevent its detection (Yang et al. 2016).

It is thus essential to properly recover the AO-PSF extensions. But compared to the PSF core, their intensity is several orders of magnitude fainter. This echoes high-contrast imaging challenges where outliers must be carefully handled to avoid any corruption of the model. Some outliers, such as dead pixels, can be identified through proper calibration (Berdeu et al. 2020), handled, and removed from the problem. But others may be random or specific to a given acquisition, such as hot pixels, cosmic rays, or model errors (signals not predicted by the forward model of the problem). This second kind of outlier must be identified on the fly, directly in the data of interest.

To reduce the incidence of such outliers, robust penalisation approaches have been developed (see Zoubir et al. 2018; Flasseur 2019, Chap. 5, for detailed overviews). They consist of replacing the conventional quadratic penalisation of the problem by a robust estimator that is approximately quadratic around zero but that grows sub-quadratically for large deviations to reduce their impact (Hogg 1979; Huber 1996).

One solution is to directly change the cost function of the problem (Yohai 1987; Hubert et al. 2008; Huber 2011). Another solution, so-called iterative reweighted least squares (IRLS, Holland & Welsch 1977; Sigl 2016), is to stay in the least squares framework but to solve a sequence of least squares problems whose weights are iteratively updated with a robust estimator. The mathematical formalism associated with such robust penalisations will be further detailed in Sect. 2.3.

This paper is a continuation of a previous pragmatic method that was initially developed to detect and characterise the third moon of (130) Elektra (Berdeu et al. 2022). In this previous work, the performance of an AO system being better in the infrared, a simple parametric model was sufficient to describe the PSF, to deconvolve the object, and to properly model and subtract its bright halo. In addition, the PSF was directly estimated on a bright moon in the data of interest, acting as a reference point source.

In this work, I intend to provide a more general and blind approach for more complex AO-PSFs with limited priors on the object and the PSF, while removing the need for a point source to be present in the field of view (FoV). There are three main objectives: (i) to reconstruct the AO-PSF in its complexity, reproducing its faint structured extensions, so-called wings, directly from the data; (ii) to deconvolve the object image with this blindly estimated AO-PSF; and (iii) to carefully model the AO residual halo in order to remove it from the data and enhance the signal of potential moons.

The proposed method, which alternates between object and PSF deconvolution, is described in Sect. 2. It is then validated on realistic simulated data in Sect. 3 and applied to real data in Sect. 4. Finally, in Sect. 5, I discuss several ways of further improving the proposed method. For readers interested in knowing more about the algorithms and their implementation, the pseudo-codes of the different steps are detailed in Appendix A.

2 Proposed method

In the isoplanetic domain (the PSF does not depend on the position in the FoV), the forward model of the problem canonically states that the image data, d (discretised onto the sensor pixel grid), is the (discrete) convolution ⋆ of the extended object, o, with the long exposure PSF, p, that combines the telescope and instrument response and the AO residuals (Fétick et al. 2020), $$\mathit{d}=(\mathit{o}\star \mathit{p})+\mathit{n},$$(1)

where n is a nuisance term. This equation is pictured in Fig. 1. As is seen in Fig. 1d, the nuisance term encompasses the classical detector readout noise and photon shot noise, but also outliers such as defective pixels (‘salt & pepper’ noise) and cosmic rays that may hit the sensor during the acquisition (orange arrows). I emphasise here that the way the problem is posed – the convolution of an extended object with a PSF – implies that the unresolved moons (point sources), are also considered to be outliers, and thus belong to the nuisance term (coloured circles). As was stated above, the objective is to split o and p only from the knowledge of d without being corrupted by n.

The PSF presented in Fig. 1c was obtained on a star, as part of the European Southern Observatory (ESO) Large Programme ID 199.C-0074 (PI: Vernazza et al. 2021) and using the imaging mode of the Zurich IMaging POLarimeter instrument (ZIMPOL, Schmid et al. 2018). This figure emphasises all the complexity of an AO-PSF. The AO cut-off frequency (orange annulus) delimits the classical two regimes of AO-PSFs: (i) the outer region (red), dominated by the atmospheric turbulence halo and left uncor-rected by the AO system, and (ii) the inner region dominated by the AO residuals. The roughly symmetric core of the PSF (white and grey) is surrounded by a wind-driven halo (green) elongated in the wind direction (temporal error of the AO system, so called servo-lag) and a speckle field produced by the non-common path aberrations (NCPAs, Vigan et al. 2019) between the science part of the instrument and the AO system (yellow and orange). The diffraction spikes produced by the spiders holding the secondary mirror are also visible (yellow cross along the diagonals). All these features are responsible for the extended and structured bright halo visible in Fig. 1a that hides the neighbouring moons.

I propose to solve Eq. (1) step by step, using minimal assumptions on the object and the PSF, summarised in Table 1. (I) Firstly, as is pictured in Fig. 2, the problem can be approximated, ‘from a distance’, as a sharp-edged flat object, Fig. 2b, convolved with a simple PSF core, Fig. 2c. In the following, such an object will qualify as 'binary' in the sense of morphological operations (Gonzalez et al. 2020). This strongly simplifies the problem in order to get a first estimated separation between the object, o, and the PSF, p. (II) Secondly, as is pictured in Fig. 3, this PSF core was used to deconvolve the main extended object, Fig. 3b. (III) Thirdly, as is pictured in Fig. 4, the deconvolution paradigm is reversed and the object was used to deconvolve the faint PSF extensions, so-called wings, Fig. 4c. Finally, after the removal of the halo model, Fig. 4a, the moon can be seen in the residuals, Fig. 4d. A detection algorithm could then be applied to find the signal of potential faint moons further hiding in the noise. This is nonetheless beyond the scope of this paper. All of these steps are further detailed in the following.

Fig. 1 Simulation of the forward model of the deconvolution. The noisy data, d (Panel a), is the convolution of an extended object, o (Panel b), with the AO-corrected PSF of the instrument, p (Panel c), plus a nuisance term, n (Panel d). This term is composed of the acquisition noises (readout + photon), defective pixels (‘salt & pepper’ pattern), cosmic ray impacts (orange arrows), and signals from potential moons orbiting the main object (coloured circles). Panels a,b: to emphasise both the main body and the surrounding halo, a dual linear scale was used to insert the main body in its halo, as is noted by the ‘/’ in the colour bars. Panels b: photo of 67P/Churyumov–Gerasimenko by Rosetta (European Space Agency). Panels c: PSF obtained with ZIMPOL on HD16469 (P. Vernazza), normalised to peak at one for the display.

2.1 Step 1. Estimation of the PSF core

The main objective of this step is to estimate a first separation between the object, o, and the PSF, p, as is shown in Fig. 2. As a consequence, it must depend on a very limited number of parameters to be robustly constrained by the data.

Concerning the PSF, the method assumes that its core can be approximated by a parametric description. Such a simple model does not include the physics of an AO system and cannot properly model the turbulent halo. It thus dilutes the energy and is consequently not quantitative (Fétick et al. 2020). Parametric model-based profiles can partially solve this issue (Fétick et al. 2019b) but they imply some knowledge of the instrument (such as, for example, its pupil shape, the number of actuators of its AO system, or the wavelength). As a consequence, they are harder to tune. But such precision is not needed for this given step. Indeed, Fétick et al. (2019a) showed that a Moffat profile (Moffat 1969) is a general but sufficient approximation to qualitatively retrieve the morphology of the observed objects, while depending only on a limited number of parameters. A two-dimensional (2D) Moffat pattern is defined as $$m(\mathit{x},\mathit{\alpha },\beta ,\theta )={\left(1+r_1^2/{\alpha }_1^2+r_2^2/{\alpha }_2^2\right)}^{-\beta }$$(2)

where β is the power parameter of the pattern, θ is its orientation, r₁ = x₁ cos θ + x₂ sin θ and r₂ = −x₁ sin θ + x₂ cos θ are the 2D coordinates in its rotated frame, and α = (α₁, α₂) are its elongations along its two axes.

As we shall see later in Sects. 3 and 4, this Moffat approximation worked well with all the tested datasets. Nonetheless, in the hypothetical case that this model is not sufficient for a given dataset, Eq. (2) can be replaced by another more adapted parametric model if needed. Doing so does not change the presented methodology as long as the number of parameters remains limited.

The PSF core is thus simply proportional to this 2D pattern: $${\mathit{p}}^{\text{core }}\left({\mathit{x}}_0,\mathit{\alpha },\beta ,\theta ,\gamma \right)\triangleq \gamma m\left(\mathit{x}-{\mathit{x}}_0,\mathit{\alpha },\beta ,\theta \right)\text{. }$$(3)

Its parameters mainly control how the object is blurred by the PSF. Most of the signal is then given by the object’s edges rather than its surface texture. If the shape of the object is overall smooth with sharp edges, applying a threshold, $\overline{d}$, on the data, d, defined as $$\left[f^{\mathrm{thr}}(\mathit{d},\overline{d})\right](\mathit{x})=\{\begin{array}{l}1\text{ if }d(\mathit{x})\ge \overline{d}\hfill \\ 0\text{ otherwise }\hfill \end{array},$$(4)

is a sufficient approximation of its support to extract the information needed to estimate the parameters of Eq. (3). Indeed, in a parametric approach with a limited number of arguments, namely the parameters of the core and the threshold on the data, a crude description of the object is sufficient because there are far more measurements than unknowns to recover.

The values of the different parameters of Eqs. (3) and (4) are obtained by minimising the following cost function, $$\begin{array}{l}{\mathcal{C}}_{\mathit{d},\mathit{w}}^{\text{core }}\left(\overline{d},\gamma ,{\mathit{x}}_0,\mathit{\alpha },\beta ,\theta \right)\\ \triangleq {\mathcal{D}}^{\text{wls}}\left(\mathit{d},{\mathit{p}}^{\text{core}}\left({\mathit{x}}_0,\mathit{\alpha },\beta ,\theta ,\gamma \right)\star f^{\text{thr}}(\mathit{d},\overline{d}),\mathit{w}\right),\end{array}$$(5)

which is a pure data fidelity term between the data, d, and the model, ${\mathit{p}}^{\text{core }}\left(x_0,\mathit{\alpha },\beta ,\theta ,\gamma \right)\star f^{\text{thr}}(\mathit{d},\overline{d})$, weighted by w. 𝒟^wls is defined as the weighted least square difference (wls) $$\begin{array}{l}{\mathcal{D}}^{\text{wls}}\left({\mathit{\phi }}_1,{\mathit{\phi }}_2,\mathit{w}\right)\\ \triangleq \frac{1}{2}{\displaystyle \sum _{\mathit{x},x^{\prime }}w}\left(\mathit{x},{\mathit{x}}^{\prime }\right)\left({\phi }_1(\mathit{x})-{\phi }_2(\mathit{x})\right)\left({\phi }_1\left({\mathit{x}}^{\prime }\right)-{\phi }_2\left({\mathit{x}}^{\prime }\right)\right),\end{array}$$(6)

where w is the inverse of the data covariance matrix. Its role is to weigh the measurements in the data fidelity term in terms of confidence.

The readout noise is assumed to be of a Gaussian statistic of variance, υ^ron. On the other hand, for fluxes higher than a few photons, the Poisson statistic of photon shot noise can be approximated by a Gaussian statistic whose variance is proportional to the incident flux (and thus approximately to the data after proper calibration and pre-reduction) by a factor, η, that accounts for the conversion between photons and analogue to digital unit (ADU) and for the quantum efficiency of the pixels (Berdeu et al. 2020). In the general case, these noise terms are independent from one pixel of the sensor to another. As a consequence, the inverse of the data covariance is diagonal, $$w\left(\mathit{x},{\mathit{x}}^{\prime }\right)=\{\begin{array}{l}w(\mathit{x})\text{ if }\mathit{x}={\mathit{x}}^{\prime }\hfill \\ 0\text{ otherwise }\hfill \end{array},$$(7)

and is given by (Mugnier et al. 2004; Fétick et al. 2019b) $$w(\mathit{x})=1/\left(\eta (\mathit{x})d(\mathit{x})+{\upsilon }^{\text{ron}}(\mathit{x})\right).$$(8)

The sum in Eq. (6) is thus a sum on x only $${\mathcal{D}}^{\text{wls}}\left({\mathit{\phi }}_1,{\mathit{\phi }}_2,\mathit{w}\right)\triangleq \frac{1}{2}{\displaystyle \sum _{\mathit{x}}w}(\mathit{x}){\left({\phi }_1(\mathit{x})-{\phi }_2(\mathit{x})\right)}^2$$(9)

I emphasise that η and υ^ron are technically vectors as their values may depend on the position in the FoV, for example because of inhomogeneities in the pixel response or detector gains. Thus, contrary to other methods that approximate the noise level as being homogeneous across the pixels, to allow the marginalisa-tion of the object and the PSF in the Fourier domain (Lau et al. 2023; Yan et al. 2023), the proposed method includes a physical model of the noise statistics that depends on the pixels and that is a function of the data intensities (measured, see Eq. (7), or expected by the model, see Eq. (18) below). This subtlety is critical: for the other methods, the unknown is the deconvolved image of the bright asteroid for which a homogeneous photon noise is a good enough approximation. But the objectives of my method are to additionally retrieve the PSF wings and remove the diffuse halo, while discarding faint outliers, which implies having a fine model of the noise according the intensity level.

As was discussed above, most of the signal is given by the blurring of the object’s edges, and is condensed in the bright core of the halo. To avoid any corruption by the faint extensions of the halo that are not modelled under the parametric PSF core assumption, a threshold is applied to the data to reweight, w, as is detailed in the Algorithm 1 of Appendix A.1. Only the pixels not shaded in red in Fig. 2a are considered in the fit.

Finally, I mention that this first step can be replaced by more classical approaches, for example using reference PSFs obtained on a star or myopic approaches (Mugnier et al. 2004). Indeed, its main objective is to get a good enough PSF core to initialise the object deconvolution of Sect. 2.2. The global PSF will be further refined later, in Sect. 2.3.

Fig. 2 Step 1 – estimation of the PSF core. The data (Panel a, red-shaded pixels excluded from the fit) is approximated by the convolution of a binary object (Panel b, contour highlighted in orange) with a parametric PSF core (Panel c). Panel d: residuals. Colour bars: see Fig. 1.

Fig. 3 Step 2 – object deconvolution. The estimated PSF core (Panel c=Fig. 2c) was used to deconvolve the object (Panel b). Panel a: model of the convolution. Panel d: residuals. Colour bars: see Fig. 1.

2.2 Step 2. Deconvolution of the main object

As is shown in Fig. 3, the objective of this step is to deconvolve the image of the object, o, assuming that the PSF, p, is known, for example from the fit of its core done in Sect. 2.1. It is a classical deconvolution problem solved by minimising the following cost function (see also, Berdeu et al. 2022): $${\mathcal{C}}_{\mathit{d},p,w}^{\text{obj}}(\mathit{o})\triangleq {\mathcal{D}}^{\text{wls}}(\mathit{d},\mathit{p}\star \mathit{o},\mathit{w})+{\mu }^{\text{obj}}{ℛ}^{\text{obj}}(\mathit{o}),$$(10)

under the physical constraint that o ≥ 0. This hard constraint also limits oscillating artefacts close to sharp edges, commonly seen in deconvolution problems (Marchis et al. 2006; Fétick et al. 2020). µ^obj is a hyperparameter to balance the regularisation, ℛ^obj(o), compared to the data fidelity term. ℛ^obj(o) was chosen to favour smooth objects with sharp edges, by encouraging the sparsity of spatial gradients (Rudin et al. 1992; Charbonnier et al. 1997) $${ℛ}^{\text{obj}}(\mathit{\phi })\triangleq {\displaystyle \sum _x\left[\sqrt{\left({\left[{\nabla }_1\phi (x)\right]}^2+{\left[{\nabla }_2\phi (x)\right]}^2+{\left[{ϵ}^{\text{obj}}\right]}^2\right)}-{ϵ}^{\text{obj}}\right]},$$(11)

where ∇₁ and ∇₂ correspond to finite difference operators along the first and second spatial dimensions of the image. є^obj > 0 is a threshold controlling the transition between a ℓ¹ -norm (edge-preserving) and a ℓ²-norm (smoothness). є^obj also ensures that Eq. (11) is differentiable at zero. Such an edge-preserving regularisation or equivalents are common in asteroid deconvolu-tion (Mugnier et al. 2004; Fétick et al. 2020; Berdeu et al. 2022; Lau et al. 2023; Yan et al. 2023). In practice, these hyperparam-eters, µ^obj and є^obj, were manually tuned. µ^obj mainly depends on the object’s characteristics. Since all asteroids share similar properties, it is not expected that this parameter changes a lot for resolved asteroids. є^obj varies with the awaited contrast on the asteroid surface: (i) for a well-resolved asteroid, values of a few percent of the object’s dynamics provide good results, (ii) for poorly resolved targets where surface details are hardly expected, values of around a thousandth or less better favour sharp-edged binary reconstructions.

At this stage, there is still no constraint on potential moons. They are thus in the deconvolved model and disappear from the residuals, as is shown in Fig. 3d. Once the object is deconvolved, it can be segmented with a threshold to remove the moons or other artefacts from the image model, as is described in Appendix A.2, Algorithm 2.

2.3 Step 3. Deconvolution of the extended PSF wings

Now that a de-blurred image of the main object, o, has been correctly estimated, it is possible to revert the problem and use this image as the convolution kernel to deconvolve the AO-PSF wings, as is shown in Fig. 4. In this approach, the moons, which are a priori unknown, are not explainable by a model of an extended object convolved with a PSF, and should thus appear in the residuals; so should defective or hot pixels or random cosmic rays impacting the sensor, which cannot be calibrated in advance. These outliers should not corrupt the reconstruction of the PSF wings. As was introduced in Sect. 1, the mean square error of Eq. (9) is consequently not adapted. I propose to tackle this problem with robust penalisation (Zoubir et al. 2018; Flasseur 2019, Chap. 5). This consists of replacing the quadratic norm of the residuals in Eq. (9) with a function, ρ, that limits the influence of large (and aberrant) values in the minimisation: $${\mathcal{D}}^{\rho }\left({\mathit{\phi }}_1,{\mathit{\phi }}_2,\mathit{w}\right)\triangleq {\displaystyle \sum _x\rho }\left(\sqrt{w(x)}\left({\phi }_1(x)-{\phi }_2(x)\right)\right).$$(12)

Of the different robust estimators, ρ, available in the literature (Holland & Welsch 1977), I chose the Cauchy function: $$\rho (r)\triangleq \frac{{\gamma }^2}{2}\ln \left(1+r^2/{\gamma }^2\right).$$(13)

If the confidence model, w in Eq. (12), is properly scaled and calibrated so that the argument $$r(x)=\sqrt{w(x)}\left({\phi }_1(x)-{\phi }_2(x)\right)$$(14)

follows a normalised Gaussian noise statistics (centred on zero and of a standard deviation equal to one) in the absence of outliers, setting γ = 2.385 ensures that minimising Eqs. (9) or (12) gives similar results. In doing so, the Cauchy function has the advantages of having no tuning parameter and of being differentiable.

In more detail, $$w^{\rho }(r)=\frac{1}{r}\frac{\partial \rho (r)}{\partial r}={\left(1+r^2/{\gamma }^2\right)}^{-1},$$(15)

the so-called robust weight in the following, is the correction factor to apply to the weights in the least square error of Eq. (9) to make it equivalent to the robust penalisation of Eq. (12). This is typically used in the context of IRLS, as was introduced in Sect. 1: $$w^{\text{irls}}(x)=w^{\rho }(r(x))w(\mathit{x})$$(16)

Figure 5a gives an example of a robust weight map, obtained by applying w^ρ of Eq. (15) to a weighted halo residual map when fitting the full PSF, p. The outliers can directly be identified in black, with equivalent weights tending towards zero. It is then possible to discard them on the fly, using a conservative threshold, ${\overline{\omega }}^{\rho }$ (typically 50 %, as is discussed in Appendix A.2), and as is shown in Fig. 5b. This is done by updating Eq. (8), with the additional remark that this equation is in fact only an approximation. Indeed, the noisy data, d, is not the real flux. A better estimate of the weighting confidence term, w, can be obtained by replacing d with the forward model $${\mathit{d}}^{\text{mod }}=\mathit{p}\star \mathit{o}$$(17)

in Eq. (8). In total, as is detailed in Algorithm 2 of Appendix A.2, the updated weight map now reads as $$w(\mathit{x})=\{\begin{array}{l}0\text{ if }{w}^{\rho }(\mathit{x})\le {\overline{w}}^{\rho }\hfill \\ 1/\left(\eta d^{\text{mod }}(\mathit{x})+{\upsilon }^{\text{ron }}\right)\text{ otherwise }\hfill \end{array}.$$(18)

It is iteratively updated in an approach similar to IRLS, Eq. (16), but with (i) a binary update when it comes to the robust estimator, and (ii) a refinement of the expected noise level with the intensity predicted by the forward model. Doing so prevents any further corruption of the deconvolved PSF wings by strong outliers.

The cost function that is iteratively minimised to reconstruct the PSF, p, is $${\mathcal{C}}_{\mathit{d},\mathit{o},\mathit{w}}^{\text{pff}}(\mathit{p})\triangleq {\mathcal{D}}^{\text{wls}}(\mathit{d},\mathit{p}\star \mathit{o},\mathit{w})+{\mu }^{\text{psf}}{ℛ}^{\text{psf}}(\mathit{p}),$$(19)

under the physical constraint that p ≥ 0. µ^psf is a hyperparameter to balance ℛ^psf(p). I define this regularisation term as $${ℛ}^{\text{psf}}(\mathit{\phi })\triangleq {\displaystyle \sum _x{\left[{\nabla }_1[\ln (\phi (x))]\right]}^2}+{\left[{\nabla }_2[\ln (\phi (x))]\right]}^2,$$(20)

which consists of the classical ℓ²-norm on the gradient but applied to the logarithm of the PSF. The former ensures a smooth reconstruction, while the latter acts as a ‘dynamic whitening’ term. Indeed, the challenge in regularising the PSF wings is that p spans multiple orders of magnitude between the PSF core vicinity and the external turbulence halo, as is seen in Fig. 1c. ℛ^psf(p) must act through this high dynamic range. Similar to w in Eqs. (9) and (12), which scales the residuals according to the admissible awaited noise level (noise whitening), the regularisation must scale the PSF wings. Since $${\nabla }_i[\ln (\phi (x))]=\left[{\nabla }_i\phi (x)\right]/\phi (x),$$(21)

I let the opportunity for the gradients in Eq. (20) to evolve proportionally to the local intensity of the total PSF, despite the large dynamics of the unknowns, both refining the PSF core and recovering its wings. As for µ^obj in Sect. 2.2, the hyperparameter µ^psf is a function of the noise level of the data, d, and was manually tuned in this study.

Technically, one could stop here if the main objective were to retrieve the faint moon hidden in the bright halo. Nonetheless, now that a better model of the AO-PSF is known, it is possible to loop back on the object deconvolution step and iterate the steps of Sects. 2.2 and 2.3 to further refine the reconstruction while robustly updating the noise model, w, via Eq. (18) in Eqs. (10) and (19). This implementation is further detailed in Algorithm 2 of Appendix A.2. This proved to be efficient for complex PSF cores, badly described by a 2D Moffat pattern, for example in the presence of motion blur or for bad seeing conditions. On top of a better object deconvolution, this also improves the residual quality close to the object edges, potentially unveiling moons very close to the main body.

Fig. 4 Step 3 – PSF wing deconvolution. The segmented object (Panel b) was used to deconvolve the PSF wings (Panel c). Panel a: model of the convolution. Panel d: residuals. Colour bars: see Fig. 1.

Fig. 5 Automatic removal of outliers with respect to the model. Panel a: equivalent weights of the robust penalisation, $w^{\rho }(\sqrt{\mathit{w}}(\mathit{d}-{\mathit{d}}^{\text{mod}}))$. Panel b: discarding of pixels (black) below a threshold, ${\overline{\omega }}^{\rho }$.

3 Validation on simulations

The method was tested on simulated data representative of real conditions. The simulated object, shown in Fig. 1b, is an image of 67P/Churyumov-Gerasimenko obtained by Rosetta (European Space Agency), downscaled to the resolution and size (between 100 and 200 milliarcseconds) that are typical of the targets observed within the ESO Large Programme described in Sect. 4.1. As was mentioned in Sect. 2, the PSF used to blur the object image is a real AO-PSF obtained with ZIMPOL on a star. Three moons with different contrasts and distances were added to the main body (coloured circles in Fig. 1d). The photon shot noise and the readout noise parameters, η and υ^ron, are assumed to be identical for all the pixels in the FoV: $$\forall x\in \mathrm{FOV},\eta (\mathit{x})=\eta \text{ and }{\upsilon }^{\text{ron}}(\mathit{x})={\upsilon }^{\text{ron}}\text{. }$$(22)

On top of these noises, random pixels are modelled as hot or dead (‘salt & pepper’ noise) and five cosmic rays are randomly simulated (orange arrows).

The results of the different steps of the method are shown in Figs. 2, 3, and 4. Figure 6 focuses on the main body decon-volution, providing zooms on the deconvolved images, o^mod, the model residuals, d − d^mod, and the object residuals, õ − κo^mod. õ stands for the shifted true object, o. Indeed, due to the PSF core deformation, the PSF in Fig. 1c may not be properly centred. By construction (see Appendix A), my method reconstructs a centred PSF. As a consequence, the deconvolved object may be translated compared to the truth. This shift is estimated by fitting a perfect 2d Moffat model in the true PSF. õ is the cubic interpolation of o translated by the position of this fitted 2D Moffat pattern. Figure 6i shows that this shift is vertical, producing the bright arcs at the top and the bottom of the object residuals. To emphasise the surface texture retrieval rather than a bias in the global intensity, the deconvolved o^mod is further normalised with $$\kappa =\underset{\tilde{\kappa }}{\mathrm{argmin}}{\displaystyle \sum _{\mathit{x}}\left|o^{\mathrm{mod}}(\mathit{x})-\tilde{\kappa }\tilde{o}(\mathit{x})\right|},$$(23)

using the absolute value in the sum to be robust to outlier regions, badly deconvolved or smoothed by the deconvolution such as edges or very bright or dark areas.

Figure 2 shows that the oversimplified model of the first step efficiently fits the bright inner part of the halo that is correctly removed from the residuals of Figs. 2d and 6f. By design, the outer part of the structured halo, mainly driven by the AO cutoff frequency, is still present in the residuals, as it cannot be explained by this simple model. As for the object, the threshold applied to the data produces a binary object (see Figs. 2b and 6b, whose edges are consistent with the blurred object, as is shown by the orange contour in Fig. 2a). This is further confirmed by the slices of Fig. 6m, the edges of the binary approximation (dark grey) being within a few pixels of the true slice of the object (light grey). As was expected, the overall object intensity is underestimated by more than 10 % because the energy diluted by the PSF wings is not integrated and of course none of the surface structure is retrieved, Fig. 6j. As for the PSF core, the radially averaged profile of the blindly estimated Moffat (blue curve in Fig. 6n) shows a very good agreement with the theoretical radial profile of the 2D Moffat (red curve) fitted directly in the true PSF.

Figure 3 shows that using the fitted PSF core is indeed sufficient to retrieve the morphology of the object, as is discussed in Sect. 2.1. Structures at the surface of the deconvolved object, Figs. 3b and 6c, are consistent with the simulated truth, Figs. 1b and 6a, as is shown by the good residuals of Fig. 6k. Nonetheless, the fact that the Moffat model cannot properly grasp the deformations of the PSF core produces two kinds of artefacts: (i) a bright edge on the left side of the object in Figs. 6c,k that suggests that this edge is ‘over-deconvolved’, and (ii) non-homogeneous structures in the residuals at the edge’s location, shown in Fig. 6g. This can also be seen in the blue slice in Fig. 6m, which shows unnatural deconvolution peaks at the edges of the reconstructed object. And as was already discussed, the fact that the PSF core does not account for the PSF wings induces an underestimation of the object’s intensity, here by more than 6 %. I also emphasise that, by construction, the halo is found in this deconvolved image of Fig. 3b, as well as the moons and other noise artefacts. They consequently appear in the modelled data, Fig. 3a, and the moons are not visible in the residuals, Fig. 3d.

All of these problems are solved after a few iterations of the alternate algorithm, once the PSF wings are properly deconvolved, as is shown by Fig. 4. Concerning the object, the zooms of Figs. 6d,l and the green slice of Fig. 6m show that, despite being slightly smoothed by the regularisation of Eq. (10), the details are quantitatively retrieved with the correct intensity within 100.2−97.3 ≤ 3 %. The sharp edges are also well defined at the correct location (see the green and light grey curves). Compared with Figs. 6c,k, the ‘over-deconvolution’ problem on the left side of the object is greatly improved. Nonetheless, the combination of (i) the threshold to segment the main object, (ii) the regularisation in Eq. (10), and (iii) the approximate knowledge of the PSF core (or equivalently the cut-off of the optical transfer function of the system) within a blind approach, limits the precision of the object’s edge recovery, as is visible in the data residuals of Fig. 6h. In addition, the spatial regularisa-tion erases the finest details on the surface of the object, as is seen in the model residuals of Fig. 6l. These are known limits of the edge-preserving regularisation (Fétick et al. 2020) but are a necessary trade-off with the PSF reconstruction to control the alternate blind deconvolution convergence. Otherwise, the residuals reach the noise level, below the percent of the data dynamic (see colour bars of Fig. 4d for Fig. 6h). There is no noticeable feature at the location of the object’s surface that would imply an over-regularisation in the deconvolution, indicating that the details smoothed by the regularisation are indeed lost in the noise level.

Concerning the PSF, all the features mentioned in Sect. 2 are quantitatively retrieved in Fig. 4c. The wind-driven halo is well defined. The brightest structures on the AO cut-off frequency ring and the diffraction spikes are clearly visible and deconvolved. Only the resolution on the speckle pattern is smoothed by the regularisation but its global structure has the correct dynamics. In Fig. 6n, the radial profile of the deconvolved PSF perfectly matches the radial profile obtained on the true PSF. As is explained in Sect. 2.3, the moons and other outliers are robustly penalised in the equivalent weights in Fig. 5a and the brightest ones are successfully rejected from the fit, Fig. 5b. Thus, they do not impact the model of the halo, Fig. 4a, and are consequently nicely visible in the residuals, Fig. 4d. Compared with the injected signals in Fig. 1d, the moons are slightly dimmed and erased from the data, showing that the method suffers from slight self-subtraction. This is nonetheless a cost worth paying to efficiently retrieve the halo and constitutes a trade-off in the regularisation of Eq. (19) between the reconstruction of the PSF structures and the moon rejection from their deconvolution.

Fig. 6 Comparison of the deconvolution with the simulated truth. Panel a: simulated true object, o. Panel b: estimated binary object, $\overline{\mathit{d}}$, applying a threshold to Panel e. Panel c: the object’s first deconvolu-tion after the estimation of the PSF core parameters, pcore. Panel d: the object’s final deconvolution after the alternate estimation of the PSF wing, p. Panel e: simulated blurred and noisy data, d. Panels f-h: model residuals, d − dmod, for the binary object (Panel f), after the first deconvolution (Panel g) and after the final deconvolution (Panel h). Panels i–l: thrice the absolute value of the object residuals with the shifted object, 3|õ − κomod|, for the true object (Panel i), for the binary object (Panel j), after the first deconvolution (Panel k), and after the final deconvolution (Panel l). Intensity scale of Panels a–e, i–l: see Fig. 1b. Intensity scale of Panels f–h: see Fig. 2d. Panel m: intensity slices along the dashed coloured line of Panels a–e. Panel n: radial profiles of the true PSF (black), its closest Moffat fit (red), the PSF core estimate (blue), and the PSF wing deconvolution (green) on a logarithmic scale.

4 Examples in real data

4.1 The multiple asteroids (216) Kleopatra and (130) Elektra

The ESO Large Programme ‘Asteroids as tracers of Solar System formation: probing the interior of primordial main belt asteroids’ (ID 199.C-0074, PI: Vernazza et al. 2021, R filter, λ = 645.9 nm, Δλ = 56.7 nm) targeted several asteroids of the Solar System’s main belt in order to characterise the albedo, surface, shape, and mass of the ones accompanied by moons. The data were obtained with ZIMPOL, mounted on the Spec-tro Polarimetric High-contrast Exoplanet REsearcher (SPHERE, Beuzit et al. 2019) of the Very Large Telescope (VLT) observatory, equipped with the SPHERE Adaptive eXtreme Optics system (SAxO, Fusco et al. 2016). Mainly designed for polarimetry acquisitions, ZIMPOL has an imaging mode. Working in visible wavelengths, it offers an unprecedented resolution of small Solar System bodies at the cost of a more complex AO-PSF. I applied my method to two targets of this Large Programme: (216) Kleopatra and (130) Elektra.

Among the 100 km class asteroids, (216) Kleopatra is quite original. It has a dumbbell shape and is orbited by two known moons, which constrain its gravitational field, and thus its mass and density (Ostro et al. 2000; Descamps et al. 2011; Shepard et al. 2018; Marchis et al. 2021). Its elongated shape and the availability of state-of-the-art deconvolutions and 3D reconstructions make it a good case with which to test the deconvolution performance of the proposed method.

(130) Elektra is also an asteroid of the main belt, surrounded by three moons (Fuksa et al. 2023). Its faintest and closest companion was missed until recently due to the lack of proper data reduction and data analysis tools (Yang et al. 2016; Berdeu et al. 2022). It was detected in archival data obtained in 2014 with the Integral Field Spectrograph of SPHERE (IFS, Claudi et al. 2008). The complexity of ZIMPOL AO-PSFs and the variability of the brightness and distance of the different moons relative to the main body represent a challenging case with which to test the proposed method. (130) Elektra was observed at different epochs in 2018 and 2019 by Vernazza et al. (2021), each epoch gathering multiple frames. In order to assess the consistency of the method results, both in terms of deconvolution and moon enhancement, two different frames of four epochs were selected. To show the versatility of the method and the large acceptability of its assumptions, I also selected additional datasets from other instruments: the IFS (ID 60.A-9362(A), PI: Yang et al. 2014, YJH filter, λ from 0.95 to 1.65 µm), the InfraRed Dual-band Imager and Spectrograph (IRDIS, Dohlen et al. 2008, ID 296.C-5038(A), PI: Yang et al. 2017, DK12 filter, λ = 2.11 µm, Δλ = 102 nm) of SPHERE, and the Near Infrared Camera 2 (NIRC2, ID U58N2, PI: de Pater et al. 2005, Kp filter, λ = 1.248 µm, Δλ = 163 nm) of the Keck II telescope.

4.2 Main body deconvolution and surface detail recovery

Figures 7 and 8 gather the results obtained on (130) Elektra and (216) Kleopatra for different epochs with different seeing conditions. The noise model parameters, (η, υ^ron), of Eq. (8), critical to scale the robust penalisation, are empirically adjusted in the data, as is described in Appendix B.

Figure 7 focuses on the final deconvolutions of (130) Elektra, obtained after the alternate algorithm. The reconstructed edges are sharp without any ‘corona artefact' (see below discussion on Fig. 8) or motion blur residual, even for 2019-07-30 (2), which is a challenging case, as is highlighted by the reference decon-volution or the important deformation of the PSF core in Fig. 9. Nonetheless, for the 2019-08-05 frames, the overshoot on the left side of the asteroid also indicates a slight over-deconvolution, as was already seen in the simulations, Sect. 3. Features at the asteroid’s surface are resolved, as is seen in the zooms of Fig. 7. The evolution of the projected shadows on the surface is visible between the frames of the best epochs, 2019-08-05 and 2019-08-06, and the corresponding slices of Fig. 7. This supports the facts that the method is consistent and that these are not deconvolution artefacts. For 2019-08-06, the slices show that the contrast in the structures decreases between the two frames, as the land forms face the Sun more. Except for the main shadow of 2019-08-05, none of these details can be seen in the reference deconvolution obtained by Vernazza et al. (2021).

Figure 8 details the different steps of the blind deconvolution algorithm on (216) Kleopatra. For each epoch, the first columns show a zoom of the data, the fitted binary model, and the deconvolved images based on the knowledge of the PSF core or the PSF wings (and the corresponding residuals) as well as intensity slices along the coloured lines. The deconvolutions obtained by Vernazza et al. (2021) are also given to serve as references, Panels b. Panels c show that the binary threshold of the data provides a good approximation of the object support, both in terms of size and shape. Then, it clearly appears that the mere knowledge of the PSF core is sufficient to retrieve the object body, but the deconvolved images of Panels f suffer from what is sometimes called ‘corona artefact’, a stepwise artefact surrounding the object edges that originates from a coupling between the edge-preserving regularisation and an incorrect description of the PSF core (Marchis et al. 2006; Fétick et al. 2019a, 2020; Lau et al. 2023). These images are similar to the reference reconstructions. For some epochs, such as 2017-07-27, 2017-08-10, and 2017-0822, the artefacts surrounding the object strongly suggest hiccups of the AO loop or even motion blur. This could be explained by a shift in the telescope pointing direction in between the four Detector Integration Times (DITs), which are averaged in one ZIMPOL-reduced frame by the data reduction pipeline (Schmid et al. 2018). This is particularly visible in the blue slices of 2017-07-27 and 2017-08-10 that cut along this motion blur.

Deconvolving the full PSF solves these issues. In the final step of the proposed method, the edges of (216) Kleopatra are de-blurred and sharp, and the duplicated images are collapsed into a single image in the cases of strong motion blur. This is also emphasised by the green slices in Fig. 8. These deconvolved images obtained via a blind deconvolution algorithm with limited priors on the object and the PSF compare well with state-of-the-art marginal deconvolutions and are consistent with the previously derived 3D model of the asteroid (Shepard et al. 2018; Marchis et al. 2021; Vernazza et al. 2021; Lau et al. 2023). The orange slices of Fig. 8 show that the reference method underestimates the object’s intensity and can only retrieve smooth and large-scale structures. Due to these limitations, only a difference in albedo between the two lobes of the asteroid was studied by Marchis et al. (2021). The small-scale structures in the green slices of Fig. 8 suggest that finer details from (216) Kleopatra’s surface can be retrieved with my method.

Interestingly, the deconvolution residuals in Fig. 8 show a line-by-line horizontal stripe pattern. This can be directly linked to the ZIMPOL imaging mode where only one out of two lines of the sensor is exposed. To produce ‘high-resolution’ images, the reduction pipeline interpolates the missing lines (Schmid et al. 2018). This implies that the proposed deconvolution algorithm reaches the level of the artefacts introduced by the data reduction. To limit the impact of this ‘zebra’ in the deconvolved image, I had to increase the regularisation parameter, µ^obj, compared to Sect. 3, limiting the performance of the method of edge and detail recovery. This explains why spatial features can still be seen in the residuals. These reduction features that correlate neighbouring pixels also mean that the assumption of independent noise previously mentioned in Sect. 2.1 for Eq. (7) is not strictly correct and that the diagonal weighting term, w, in Eq. (19) should actually be the inverse of a covariance matrix, as in Eq. (6). This effect was nonetheless neglected in this study.

Fig. 7 Zoom on the deconvolved image of (130) Elektra in Fig. 9 obtained with ZIMPOL. For information, the deconvolutions (green) are compared with the data (red) and the reference (orange) obtained Vernazza et al. (2021). To emphasise the gain in resolution, all the figures are normalised to the same linear intensity scale. For each epoch, two different frames (plain (1) and dashed (2) curves) are given.

4.3 Point spread function reconstruction and halo removal

A quick look at the halo brightness in the data column of Fig. 8 shows that the turbulence conditions change between the different epochs. This directly impacts the deconvolved AO-PSF wings and their contrast. For the two epochs 2017-07-14 and 2017-08-10, the PSF is even mainly dominated by the wind-driven halo. Nonetheless, the blind PSF deconvolution successfully retrieves the AO cut-off limit at the correct location as well as its brightest speckles, despite the absence of an instrumental prior in the PSF model. For good seeing conditions, in spite of the elongated shape of (216) Kleopatra, the reconstructed speckles appear nicely resolved and split. In addition, the orientation and elongation of the wind-driven halo do not appear to correlate with (216) Kleopatra’s orientation (see especially 201707-14 and 2017-08-22), as one could have feared from a potential cross-talk during the blind deconvolution. This overall consistency favours the robustness of the method in connection with the shape of the object and with the turbulence conditions.

Looking at Figs. 9 and C.1, the conclusions are similar with (130) Elektra, with the additional supporting fact that for each epoch the two reconstructed PSFs are extremely consistent in terms of turbulence features (contrast and wind-driven halo). Furthermore, the rotation of the spider diffraction spikes is clearly visible (see 2019-07-30 and 2019-08-06 in particular).

After halo removal, the two moons of (216) Kleopatra are clearly visible at all epochs in the halo residuals, as is shown in the last column of Fig. 8 (coloured circles). The transition between a readout noise-limited regime to a photon shot noise-limited regime is also nicely visible: the noise level increases in the object’s vicinity where the halo is very bright in the data, especially for bad seeing conditions (2017-07-14, 2017-08-10, and 2017-08-22). Looking at the halo residuals in Fig. 9, this transition is also visible, especially for the bad seeing cases of 2019-07-30 and 2019-08-04 with an increased noise level in the object’s vicinity. The three moons of (130) Elektra are visible in all the presented ZIMPOL frames (coloured circles), even when the deconvolution is not perfect (see 2019-07-30). Indeed, by construction, even if the object and the PSF are not perfectly split, moons (point sources) remain outliers of the model and should appear in the residuals. This is an encouraging result that shows the robustness of the approach. The counterparts of the moons in the robust weight maps are also nicely visible. The orange arrows point towards discarded outliers, such as defective pixels or random cosmic rays as well as camera defects impacting some columns of the detector.

The last lines of Fig. 9 present results on data obtained with other instruments and reference PSFs measured on a star. They show that the method is robust for a large variety of AO-PSFs.

Looking at the IFS line, my method is able to reproduce the results of my previous pragmatic approach (Berdeu et al. 2022), extracting the three moons from the diffuse halo. The poor AO performance of this dataset explains why the diffraction rings of the reconstructed PSF differ from the IFS reference.

Nonetheless, the diffraction spikes are visible and the reconstructed multi-lobed core nicely matches the PSF produced by the brightest moon (red circles): it presents a strong secondary lobe, at the top left of the main core. A similar complex PSF core is retrieved with ZIMPOL for the 2019-07-30 (2) frame (secondary lobe on the right). This shows that the method is, to some extent, robust for multi-lobed PSF cores diverging from a simple Moffat pattern, slightly relaxing this potential limitation of the method. The robust weights obtained with the IFS appear quite low. This comes from the wrong noise model fitted in the reduced data via the method presented in Appendix B. Indeed, the reduction performed by PIC of IFS data (Berdeu et al. 2020) is based on a spatial regularisation that biases and correlates the noise in the reconstruction. Such an incorrect noise model can be handled by relaxing the constraint on ${\overline{\omega }}^{\rho }$ to lower values for automatic outlier rejection.

The IRDIS data are very noisy, with a dense pollution with defective pixels (common with infrared detectors), while (130) Elektra only spans a few pixels. Despite this unfavourable situation, the first Airy rings of the PSF are reconstructed with the correct dynamics and with negligible signal in the residuals. The first ring is fragmented due to poor seeing conditions. Unfortunately, only the brightest moon is visible because of the high noise level and high density of unusable pixels.

Finally, with NIRC2, (130) Elektra is barely resolved, with an angular size similar to that of the first diffraction ring of the instrument PSF. Despite this challenging situation, the characteristic hexagonal shape of this ring, due to the fragmentation of the Keck II mirror with hexagonal segments, is nicely recovered. So are the secondary diffraction peaks (two in green at the top right of the core for the first order and a few in yellow for the second order) that produce replicated images of (130) Elektra in the saturated data. Concerning the moons, only the brightest is visible in the residuals, highlighting the need for more frames and epochs and dedicated detection algorithms to increase the signal-to-noise ratio.

Fig. 8 Examples of deconvolution on (216) Kleopatra at different epochs. First column – zoom on the main body. Panels a,b: data (top, red slice) and reference deconvolution (bottom, orange slice) obtained by Vernazza et al. (2021). Panels c,d: estimated binary object after the estimation of the PSF core parameters (top, grey slice) and residuals (bottom). Panels e,f: deconvolution with the PSF core (top, blue slice) and residuals (bottom). Panels g,h: deconvolution after the alternate estimation of the PSF wing (top, green slice) and residuals (bottom). To highlight the faint intensities, all the object representations share the same square root stretch intensity scale normalised to the data’s maximal value. The residuals are displayed with the colour scale of Fig. 2b. Second column - slices along the coloured lines of the first column. Third column - data displayed with the dual linear scale of Fig. 1a. Fourth column – PSF displayed with the colour scale of Fig. 1c (peak normalised to one). Fifth column – residuals of the halo model displayed with the colour scale of Fig. 1d. For information, the deconvolved main body image was inserted.

Fig. 9 Examples of robust halo removal on (130) Elektra with different instruments (see also Fig. C.1). First (resp. second) column - blurred data, d (resp. deconvolved image, o, and halo model, dmod = p ★ o), normalised and displayed with a dual linear scale. Third column - deconvolved PSF, p (peak normalised to one). Fourth column - residuals of the halo model, d − dmod. Fifth column - map of the equivalent robust weights, $w^{\rho }(\sqrt{\mathit{w}}(\mathit{d}-{\mathit{d}}^{\text{mod}}))$. Orange arrows: camera defects, dead pixel clusters, or cosmic rays. When visible, the three moons of (130) Elektra are circled. The red arrows point towards the outer moon when it is outside the FoV. Last line – examples of PSF on a star for each instrument.

5 Conclusions and perspectives

In this paper, I have presented a novel blind deconvolution technique that can retrieve both the object and the faint wings of AO-PSFs. The only main priors on the object are that it be contrasted with a smooth shape and sharp edges. Contrary to most of the existing blind deconvolution algorithms, it is not based on a parametric model of the PSF, allowing one to reconstruct a large variety of shapes and structures, as featured by PSFs after an AO correction. The proposed inverse problem approach is an alternate deconvolution with an edge-preserving regulari-sation for the object and a regularisation on the smoothness of the gradients for the PSF with intensity whitening. Tested on simulated and real data, my method outperforms state-of-the-art techniques (Fétick et al. 2019b; Vernazza et al. 2021; Lau et al. 2023; Yan et al. 2023), by both efficiently reconstructing details at the surface of the objects and faithfully retrieving the features of the AO-PSFs. It produces a physical and realistic model of the bright halo surrounding the object, which can then be carefully removed. This enhances faint moons in its close vicinity with a better contrast and at closer distances than previous halo-subtraction methods (Assafin et al. 2008; Yang et al. 2016; Pajuelo et al. 2018; Berdeu et al. 2022).

Implemented in a general framework and with a limited number of assumptions, my method can be applied to different instruments mounted on different telescopes equipped with different AO systems. It paves the way for the study of resolved asteroids or moon systems around the giant planets of our Solar System via archival data of different instruments for which previous methods could not reveal the companions. It also prepares the arrival of the future generation of ELTs, whose complex and sensitive PSFs and AO systems will require robust and versatile processing techniques. Nonetheless, as is discussed below, there are several ways to further improve this method.

Solving the dual deconvolution problem implies carefully tuning the method. This is mainly done through the hyperparam-eters of the regularisations on the object and the PSF, namely µ^obj, є^obj, µ^psf, and ${\overline{\omega }}^{\rho }$. In particular, one could decide to first keep a strong regularisation on the object to prioritise the PSF reconstruction, and to relax it the second time, once the PSF has been correctly estimated, to better retrieve the object’s details. Providing an automatic tuning of such parameters to get a fully unsupervised algorithm is an active field of research (Thé et al. 2022) and implementing such methods was beyond the scope of this paper. In this work, these parameters were manually tuned and fixed once and for all for each instrument. The fact that they worked for all the tested datasets is encouraging and shows their good admissible range in terms of noise level, object shape, and seeing conditions. This thus limits the need for further manual tuning by the user.

It was emphasised in Sect. 4 that lots of datasets suffer from poor AO-PSF quality or even motion blur. Robustness is obtained via a higher regularisation level, but this can degrade the image quality or produce over-deconvolution. To increase the resolution and the deconvolution quality, and since the targets are very bright, one should work on the individual ZIMPOL DITs rather than the stacked image. The results also show that the method reaches the limits of the ZIMPOL reduction pipeline, whose artefacts dominate the residuals in some datasets. Implementing such changes implies modifying the reduction pipeline and was beyond the scope of this paper.

As is shown in Sect. 4, most of the outliers, such as defective pixels or cosmic rays, are successfully and robustly identified and discarded on the fly by the proposed algorithm. Nonetheless, some of them are missed, and so are the faintest moons. As is discussed in Sect. 3, the proposed method consequently suffers from slight self-subtraction, which could dim the signal of potential moons. Even if it is less pronounced than with other methods (Yang et al. 2016; Pajuelo et al. 2018), some precautions must be taken when fitting the photometry of a moon that can be underestimated, or when announcing the absence of moons below a given noise-limited contrast. It is possible to assess these biases by injecting false moons as is done in the standard exoplanet detection algorithm to quantitatively estimate the contrast performance (Flasseur et al. 2018). This was beyond the scope of the paper, but should be kept in mind for further studies.

Finally, the natural next step of this method is to use to the reconstructed PSF within matching filter algorithms to detect even fainter moons in the residuals. Combined with the refined noise model and with the robust weights to limit false positive detections by discarding outliers, this would push the detection limits even further. Echoing the challenges of exoplanet detection, it comes nonetheless with the additional difficulty of the rapid motion of the moons around the main object (Showalter et al. 2019; Berdeu et al. 2022) and is an open research field.

Acknowledgements

The author would like to warmly thank the Referee Laurent Mugnier for his pertinent feedback and numerous remarks and suggestions that greatly contributed to improving the quality and clarity of the article. The author would like to thank Éric Thiébaut and Ferréol Soulez for the fruitful discussions on robust detection algorithms and adequate prior for high dynamic unknowns as well as Benoît Carry for his careful proofreading. All the ZIMPOL reductions are based on public data acquired for the ESO Large Programme ID 199.C-0074 (Vernazza et al. 2021, ‘Asteroids as tracers of Solar System formation: probing the interior of primordial main belt asteroids’) and available at https://observations.lam.fr/astero/. The IFS and IRDIS reductions are based on public data provided by the ESO Science Archive Facility. The author would like to thank Antoine Kaszczyc and Maud Langlois who provided the PSFs and reduced the IRDIS raw data. The NIRC2 reductions are based data provided by the Keck Observatory Archive (KOA), which is operated by the W.M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration. The author would like to thank Kate Minker who reduced the raw data and provided the PSF. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 101004719.

Appendix A Pseudo-codes of the method

This appendix gathers the pseudo-codes of the different steps of the proposed method, described in Sect. 2. In practice, all the code is implemented in MATLAB within the open-source GlobalBioIm framework (Soubies et al. 2019). I also detail in this appendix how the different cost functions are minimised.

Appendix A.1 Estimation of the point spread function core

The different parameters of the first step of Sect. 2.1 were alternately fitted, as was described in Algorithm 1. Indeed, applying a threshold to the data is not a differentiable operation, complicating the minimisation with iterative gradient descent methods. A first set of parameters of the PSF core was guessed, Line 3, on a coarse grid of threshold values, Line 2. These values were then refined by alternating five fits on (i) the threshold on the data, along with the amplitude parameter, Line 5, and (ii) the PSF core parameters, Line 6. The optimisation problems were solved with the simplex search method of Lagarias et al. (1998). For each problem, 200 iterations were performed. The parameters of the Moffat profiles were initialised with α = 3 and β = 1.6, which are standard values (see Fétick et al. 2019a). To reweight, w, Line 1, a general threshold of 2.5 % was applied to the data. From the perspective of the dynamic of the PSF, this value roughly corresponds to the transition between its wind-driven core and its wings (see Fig. 1c). Only the meaningful pixels in the brightest part of the halo, not shaded in red, Fig. 2a, are thus kept in the data fidelity term.

Appendix A.2 Alternate object and point spread function deconvolution

Algorithm 2 details the implementation of the dual deconvolu-tion algorithm of the object and the PSF wings of Sects. 2.2 and 2.3. The role of Line 8 is to define the main object support. I empirically chose a threshold of ${\overline{d}}^{\text{sup }}$ between 15 % and 25 % of the maximal value of the median-filtered deconvolved object, ƒ^med(o) (kernel of 5 × 5 pixels, as in Appendix 3). Such a threshold is low enough to account for partially illuminated pixels on the object’s edge, while being conservative enough to correctly constrain and tackle the ‘corona artefacts’ mentioned in Sect. 4. In practice, this threshold was applied only in the vicinity of the main body. To further avoid crenellated edges on the segmented image, this support was dilated by one pixel. The rest of the FoV was automatically excluded to avoid bright moons or artefacts, such as cosmic rays in the deconvolved image, exceeding this threshold.

Conversely, for Line 6 and Line 13, the ${\overline{\omega }}^{\rho }$ threshold was applied only in the halo to exclude the moons or artefacts, such as cosmic rays, but not in the main body vicinity, to keep this area (see Fig. 5b). Indeed, this area is critical to fit the PSF core features but hard to fit due to the artificial threshold on the deconvolved image and the edge-preserving regularisation, as is seen in the robust weights of Fig. 9. I typically used the conservative value of 50 %, as it is important to remove any moon suspicion while the deconvolution of the smoothly varying PSF is robust for a high rate of missing data. In the presence of strong outliers (dead pixels) on the object support or its vicinity, as in Fig. 5b, it is still possible to apply a less conservative threshold to this region, typically below 10 %.

The two minimisation problems of Line 8 and Line 14 were solved with 1000 iterations of the variable metric with limited memory-bounded algorithm (VMLM-B, Thiébaut 2002), a limited-memory quasi-Newton method with Broyden-Fletcher-Goldfarb-Shanno updates (BFGS, Nocedal 1980) that handles bound constraints. In the case of strong motion blur, typically n^alt = 30 alternate iterations are looped.

In Line 6 and Line 13, the weights were robustly updated only after a few alternate loops, n^wgt = 5. I would like to remark here that if a proper calibration of the reduced data has been made, the robust weights, Line 2, can be initialised with a map of already known aberrant pixels. This was done for the IRDIS dataset.

Appendix B Empirical noise model

The attentive reader may notice that the proposed method implies knowledge of the noise model parameters, (η, υ^ron), in Eq. (8). If they are known in the simulated case of Sect. 3, they are a priori unknown for the real data of Sect. 4. In principle, the noise model can be obtained through proper calibration and analysis of both calibration and science data (see for example, Berdeu et al. 2020; Denneulin et al. 2021). Nonetheless, most of the data used in this paper are ZIMPOL data from the ESO Large Programme described in Sect. 4.1 that were already reduced (Vernazza et al. 2021). It was beyond the scope of this paper to fully recharacterise the ZIMPOL sensor noise and rerun a reduction procedure. As a consequence, I chose to implement a pragmatic approach that could empirically fit these parameters directly in the data of interest and that could be easily applied to reduced data from other instruments. The principle of the method is pictured in Fig. B.1.

Fig. B.1 Empirical fit of the noise model, ηd(x) + υron, directly from the data. Panel a: Visualisation of the arcs centred on the main object on which the model is fitted. Panel b: Approximated map of the noise after removing a median filter on the data. Panel c: Estimated noise model (black curve) compared with the true model (dashed curve). In blue (resp. red): points that are kept (resp. discarded) by the robust fit. The green box is a zoom on the relevant points.

The noise model of Eq. (8) is basically an afflne law that links the data intensity and its associated variance. With the assumption that the noise parameters are constant over the FoV, Eq. (22), the two noise parameters, (η, υ^ron), can be fitted on the full data image using the diversity of information, as is implemented in Algorithm 3.

First, as is shown in Fig. B.1b, the noise map was approximated by removing from the data, a smoothed image obtained after applying a median filter of 5 × 5 pixels, Line 1. The pixels were then split among n^arc arcs (with a typical width of 5 pixels and length of 20 pixels for the ZIMPOL data) centred on the main body, shown in Fig. B.1a. The intensities (resp. variances) of the model were empirically estimated by averaging the data (resp. by computing the variance of the noise map) on each arc, Line 3 (resp. Line 4). To be robust with regard to the presence of outliers among this point cloud, Fig. B.1c, the affine law was robustly fitted via an approach based on the median absolute deviation of the model from the empirical points, Line 7, a robust estimator of the standard deviation (Huber 2011).

This method was tested on the simulated data of Sect. 3, for which the true noise parameters are known. The results of the fit, given in Fig. B.1c, show that the noise model was successfully estimated.

I point out here that these noise parameters could be iteratively refined during the iterations of the alternate algorithm once the model of the convolution between the object and the PSF, namely the theoretical intensity map, and the model of the halo residuals, namely the theoretical noise map, are better constrained.¹ This was nonetheless beyond the scope of this paper and is kept for a future work. Indeed, the homogeneity of the weight maps in Figs. 9 and C.1 shows that the proposed method is a good enough approximation, correctly whitening both the readout noise and the photon shot noise.

Appendix C Additional frames of (130) Elektra

Figure C.1 comes in addition to Fig. 9 and presents the second frame of the couples introduced in Fig. 7. This figure supports the consistency of my method that produces similar results within each epoch, in terms of PSF structures and dynamics as well as moon enhancement.

Fig. C.1 Additional frames of (130) Elektra with ZIMPOL (see caption of Fig. 9).

References

All Tables

All Figures

Fig. 1 Simulation of the forward model of the deconvolution. The noisy data, d (Panel a), is the convolution of an extended object, o (Panel b), with the AO-corrected PSF of the instrument, p (Panel c), plus a nuisance term, n (Panel d). This term is composed of the acquisition noises (readout + photon), defective pixels (‘salt & pepper’ pattern), cosmic ray impacts (orange arrows), and signals from potential moons orbiting the main object (coloured circles). Panels a,b: to emphasise both the main body and the surrounding halo, a dual linear scale was used to insert the main body in its halo, as is noted by the ‘/’ in the colour bars. Panels b: photo of 67P/Churyumov–Gerasimenko by Rosetta (European Space Agency). Panels c: PSF obtained with ZIMPOL on HD16469 (P. Vernazza), normalised to peak at one for the display.

In the text

	Fig. 2 Step 1 – estimation of the PSF core. The data (Panel a, red-shaded pixels excluded from the fit) is approximated by the convolution of a binary object (Panel b, contour highlighted in orange) with a parametric PSF core (Panel c). Panel d: residuals. Colour bars: see Fig. 1.
In the text

	Fig. 3 Step 2 – object deconvolution. The estimated PSF core (Panel c=Fig. 2c) was used to deconvolve the object (Panel b). Panel a: model of the convolution. Panel d: residuals. Colour bars: see Fig. 1.
In the text

	Fig. 4 Step 3 – PSF wing deconvolution. The segmented object (Panel b) was used to deconvolve the PSF wings (Panel c). Panel a: model of the convolution. Panel d: residuals. Colour bars: see Fig. 1.
In the text

	Fig. 5 Automatic removal of outliers with respect to the model. Panel a: equivalent weights of the robust penalisation, $w^{\rho }(\sqrt{\mathit{w}}(\mathit{d}-{\mathit{d}}^{\text{mod}}))$. Panel b: discarding of pixels (black) below a threshold, ${\overline{\omega }}^{\rho }$.
In the text

Fig. 6 Comparison of the deconvolution with the simulated truth. Panel a: simulated true object, o. Panel b: estimated binary object, $\overline{\mathit{d}}$, applying a threshold to Panel e. Panel c: the object’s first deconvolu-tion after the estimation of the PSF core parameters, pcore. Panel d: the object’s final deconvolution after the alternate estimation of the PSF wing, p. Panel e: simulated blurred and noisy data, d. Panels f-h: model residuals, d − dmod, for the binary object (Panel f), after the first deconvolution (Panel g) and after the final deconvolution (Panel h). Panels i–l: thrice the absolute value of the object residuals with the shifted object, 3|õ − κomod|, for the true object (Panel i), for the binary object (Panel j), after the first deconvolution (Panel k), and after the final deconvolution (Panel l). Intensity scale of Panels a–e, i–l: see Fig. 1b. Intensity scale of Panels f–h: see Fig. 2d. Panel m: intensity slices along the dashed coloured line of Panels a–e. Panel n: radial profiles of the true PSF (black), its closest Moffat fit (red), the PSF core estimate (blue), and the PSF wing deconvolution (green) on a logarithmic scale.

In the text

Fig. 7 Zoom on the deconvolved image of (130) Elektra in Fig. 9 obtained with ZIMPOL. For information, the deconvolutions (green) are compared with the data (red) and the reference (orange) obtained Vernazza et al. (2021). To emphasise the gain in resolution, all the figures are normalised to the same linear intensity scale. For each epoch, two different frames (plain (1) and dashed (2) curves) are given.

In the text

Fig. 8 Examples of deconvolution on (216) Kleopatra at different epochs. First column – zoom on the main body. Panels a,b: data (top, red slice) and reference deconvolution (bottom, orange slice) obtained by Vernazza et al. (2021). Panels c,d: estimated binary object after the estimation of the PSF core parameters (top, grey slice) and residuals (bottom). Panels e,f: deconvolution with the PSF core (top, blue slice) and residuals (bottom). Panels g,h: deconvolution after the alternate estimation of the PSF wing (top, green slice) and residuals (bottom). To highlight the faint intensities, all the object representations share the same square root stretch intensity scale normalised to the data’s maximal value. The residuals are displayed with the colour scale of Fig. 2b. Second column - slices along the coloured lines of the first column. Third column - data displayed with the dual linear scale of Fig. 1a. Fourth column – PSF displayed with the colour scale of Fig. 1c (peak normalised to one). Fifth column – residuals of the halo model displayed with the colour scale of Fig. 1d. For information, the deconvolved main body image was inserted.

In the text

Fig. 9 Examples of robust halo removal on (130) Elektra with different instruments (see also Fig. C.1). First (resp. second) column - blurred data, d (resp. deconvolved image, o, and halo model, dmod = p ★ o), normalised and displayed with a dual linear scale. Third column - deconvolved PSF, p (peak normalised to one). Fourth column - residuals of the halo model, d − dmod. Fifth column - map of the equivalent robust weights, $w^{\rho }(\sqrt{\mathit{w}}(\mathit{d}-{\mathit{d}}^{\text{mod}}))$. Orange arrows: camera defects, dead pixel clusters, or cosmic rays. When visible, the three moons of (130) Elektra are circled. The red arrows point towards the outer moon when it is outside the FoV. Last line – examples of PSF on a star for each instrument.

In the text

Fig. B.1 Empirical fit of the noise model, ηd(x) + υron, directly from the data. Panel a: Visualisation of the arcs centred on the main object on which the model is fitted. Panel b: Approximated map of the noise after removing a median filter on the data. Panel c: Estimated noise model (black curve) compared with the true model (dashed curve). In blue (resp. red): points that are kept (resp. discarded) by the robust fit. The green box is a zoom on the relevant points.

In the text

	Fig. C.1 Additional frames of (130) Elektra with ZIMPOL (see caption of Fig. 9).
In the text