© A. González et al., Published by EDP Sciences 2016

1. Introduction

Deconvolution is a ubiquitous data processing method that arises in a variety of applications, e.g., biomedical imaging (Jefferies et al. 2002), astronomy (Prato et al. 2012; Shearer 2013), remote sensing (Dong et al. 2011), and video and photo enhancement (Fergus et al. 2006). Formally, this problem amounts to recovering a signal x from blurred and corrupted observations y:$$\mathbf{y}=\mathrm{\Theta }(\mathbf{h}\otimes \mathbf{x}),$$(1)where ⊗ stands for the convolution operation between x and some other signal h, while Θ is a general function accounting for some corrupting noise, e.g., Θ(u) ≔ u + n under an additive corruption model with some (Gaussian) noise n.

For instance, when a scene is captured by an optical instrument, the observation is a blurred or degraded version of the original image, corrupted by noise and by some effects due to the instrument (motion, out-of-focus, light scattering, etc.). In such a case, the imaging system is usually assumed to be well represented by the sensing model (1) where the ideal image x is blurred by a Point Spread Function (PSF) h. Implicitly, the PSF describes the response of the imaging device to a Dirac point source, i.e., the impulse response of the instrument, while Θ can distinguish different sources of noise contaminating the image: read-out, photon counting, multiplicative and compression noise, among others.

If the true PSF can be determined in advance, then it is possible to recover the original, undistorted image by convolving the acquired image with the inverse PSF. This is the so-called “known-PSF deconvolution”. In most practical applications, however, finding the true PSF is impossible and an approximation must be made. As the acquired image is corrupted by various sources of noise, both the PSF and a denoised version of the image should be recovered together. This process is known as “blind deconvolution”.

Since there are infinite combinations of image and filter that are compatible with the distorted observations, the blind deconvolution problem is severely ill-posed. One way to reduce the number of unknowns is to introduce a forward (parametric) model of the PSF. For instance, Oliveira et al. (2007) have modeled the motion blur by a straight line, where the number of unknowns is reduced to two: the line length and angle. Babacan et al. (2009) have modeled smoothly varying PSFs by using a simultaneous autoregressive prior, where the only parameter to estimate is the variance of a Gaussian function. However, such methods are limited to specific applications as they can only recover the expected model of the PSF, excluding the estimation of a generic non-parametric PSF. Furthermore, due to the ill-posedness of the blind deconvolution problem, a slight mismatch between the specified model and the true PSF can lead to poorly deconvolved images.

Blind deconvolution is also very sensitive to the noise present in the observations. While some works (Ayers & Dainty 1988; Gburek et al. 2006) have used fast and efficient techniques such as a simple inverse filtering to recover both the PSF and the undistorted image from the blurred observations, these methods present low tolerance to noise (Fish et al. 1995).

Robust blind deconvolution methods, on the other hand, optimize a data fidelity term, which is based on the acquisition model, stabilized by some additional regularization terms. If the observations are corrupted by an additive Gaussian noise, the deconvolution problem can be solved by a regularized least-squares (LS) minimization (Almeida & Almeida 2010). In the presence of Poisson noise, the problem can be formulated using the Kullback-Leibler (KL) divergence as the data fidelity term (Fish et al. 1995; Prato et al. 2012, 2013), which represents a more complex function to minimize. To avoid dealing with the KL divergence, the Poisson corrupted data can also be handled through a Variance Stabilization Transform (VST) (Dupé et al. 2009; Shearer et al. 2012; Shearer 2013). The VST provides an approximated Gaussian noise distributed data and thus allows working with a regularized LS formulation.

In this work, blind deconvolution is used to recover both the PSF and the undistorted image from blurred observations acquired by solar extreme ultraviolet (EUV) telescopes. We use the proximal alternating minimization method recently proposed by Attouch et al. (2010) in the context of additive Gaussian noise. This algorithm allows handling LS problems for a large class of regularization functionals. The advantage of this method over usual alternating minimization approaches, e.g., Fish et al. (1995); Almeida & Almeida (2010), is that it provides theoretical convergence guarantees and it is general enough to include a wide variety of prior information. To our knowledge, it is the first time that blind deconvolution is solved using the proximal algorithm of Attouch et al. (2010).

In optical telescopes, the PSF originates from various instrumental effects, such as: optical aberrations (spherical, astigmatism, coma), diffraction (produced by, e.g., an entrance filter mesh), scattering (from, e.g., micro-roughness of the mirrors resulting in long-range diffuse illumination), charge spreading, etc. All these effects “spread” light, i.e., incident photons which would otherwise focus to a single point of the focal plane, may get detected at a location that is a bit shifted or even far away. This does affect different types of measurements, such as the photometry of fine features, see, e.g., DeForest et al. (2009).

The PSF of solar telescopes is usually modeled using pre-flight instrument specifications (Martens et al. 1995; Grigis et al. 2013). After building a specific model for a given instrument, few parameters need to be estimated in order to recover the PSF. For instance, Gburek et al. (2006) fitted the central pixels of the TRACE PSF to a Moffat function, while DeForest et al. (2009) modeled the long-range effect of the TRACE PSF as the sum of a measured diffraction pattern with a circularly symmetric scattering profile. The PSFs from the four EUVI instrument channels on STEREO-B were studied by Shearer et al. (2012) and by Shearer (2013), where the long-range scattering effect was assumed to follow a parametric piece-wise power-law model. A similar method was used to estimate the PSF of the SWAP instrument on board the PROBA2 satellite (Seaton et al. 2013). Finally, Poduval et al. (2013) provided a semi-empirical model for the Atmospheric Imaging Assembly (AIA) on board SDO, similar to the one used by DeForest et al. (2009) for TRACE. In all these works, the complete set of parameters that were fitted is in the range of 5–11 values, just a few compared to the resolution of the PSF (millions of pixels).

In addition to building a parametric model of the PSF using the specifications of the instruments, some works (DeForest et al. 2009; Shearer et al. 2012; Poduval et al. 2013) have used the information from a celestial transit to help in inferring the PSF of a given instrument. A celestial transit, shown in Figure 1, is an astronomical event where the Moon or a planet moves across the disk of the sun, hiding a part of it, as seen from the telescope. Since the EUV emission of transiting celestial bodies is zero, any apparent emission is known to be caused by the instrument. Therefore, this type of event provides strong prior information that allows the regularization of the blind deconvolution problem.

Let us note that all these parametric methods present the disadvantage of being specific for a given instrument and depending on a good characterization of the PSF, which is not always possible (Meftah et al. 2014).

1.1. Contribution

Central to our work is the information provided by the transit of a celestial object whose apparent boundary on the recorded image can be predicted with high (subpixel) accuracy. This is typically the case for the observation of Moon or Venus transits for which ephemeris allows us to precisely know their apparent boundaries at the recording time, provided of course that (i) the small variations of the object topography around a perfect disk (e.g., mountains, craters) are smaller than a pixel width and (ii) that such a transiting object has no atmosphere that could blur its apparent boundary (e.g., as for an Earth transit). As will be clear below, these hypotheses have been respected for at least two cases: for the instrument SECCHI/EUVI during the 2007 Lunar transit and for SDO/AIA during the 2012 Venus transit. For these two situations, and for any other transit respecting the conditions above, we know a priori the values for a set of pixels in the image and their exact location. Our work uses such information as constraint in the function to be optimized for blindly deconvolving images. Notice that our approach could also be applied to other situations where pixel values and exact location are known a priori (e.g., dark calibration patterns in Computer Vision applications).

Contrarily to previous works in solar physics, the blind deconvolution method demonstrated in this paper is not based on a specific model of the PSF, but infers it from the observed data. This allows the method to be used for estimating the PSF of any instrument provided it has a prior information on a set of pixel values. Since no model is imposed on the PSF, the large number of unknowns makes the problem computationally intractable. We thus focus on the estimation of the PSF core, which is defined as the central pixels encompassing at least 99% of the total PSF energy. Note that this thresholding follows the radial energy characterization of the SDO/AIA PSF in the paper of Poduval et al. (2013) and a similar study of the available PSFs of SECCHI/EUVI (e.g., Shearer et al. 2012). The PSF core accounts for charge spreading, optical aberrations, and diffraction effects. Note that the same algorithm would be able to handle a parametric model of the PSF, provided a precise instrument characterization is available. For simplicity, we consider an average additive Gaussian noise model that gathers all sources of noise present in the observation. A detailed discussion on how to handle more general noise models is provided in Section 7.

The proposed method is first tested through simulated realistic scenarios. Results show the ability to estimate a given PSF with high quality. We also show the importance of considering multiple transit observations in order to provide a better conditioning to the filter estimation problem.

We validate our method on AIA and EUVI observations using solar transits. The validation of the recovered filters is based on the reduction of the celestial body apparent emissions. The estimated PSFs were also tested for images containing active regions. Results show that, when recovering the PSF core, the proposed non-parametric scheme can achieve similar results to parametric methods. Also, we consider the case where the parametric PSF is incorporated in the imaging model as an additional filter convolving the actual image (Shearer et al. 2012; Shearer 2013). We show that the resulting parametric/non-parametric PSF provides better results than both parametric and non-parametric PSFs estimated independently. Let us note that, since parametric methods are based on the physics of the instrument, they tend to provide a more precise PSF containing both the core and the diffraction peaks. However, when the telescope presents unknown or unexpected behavior, an accurate model of the filter cannot be built. In such cases, we require non-parametric models, as the one presented in this paper, which due to their flexibility can handle better some instrument’s properties.

1.2. Outline

In Section 2, we describe the discrete forward model and the noise estimation. In Section 3, we present the formulation of the blind deconvolution problem based on the image and filter prior information. Section 4 describes the alternating minimization method used for estimating the filter and the image. It briefly presents the method, including the initialization, the parameter estimation, the stopping criteria, and the numerical reconstruction. In Section 5, we present a non-blind deconvolution method that is used for the experimental validation. In Section 6, we present the results on both synthetic and experimental data for the SDO/AIA and SECCHI/EUVI telescopes. Finally, Section 7 presents a discussion of the obtained results, providing some perspectives for future research work.

1.3. Notations

We denote by , , and the sets of natural, integer, and real numbers, respectively. The set of the non-negative real numbers is denoted by . Most of the domain dimensions are denoted by capital roman letters, e.g., M, N, … Vectors and matrices are associated with bold symbols, e.g., or , while scalar values are associated with lowercase light letters. The ith component of a vector u reads either u_i or (u)_i, while the vector u_i may refer to the ith element of a vector set. The vector of 1’s in is denoted by 1_D = (1, … ,1)^T. The set of indices in is [D] = {1, …, D}. The support of a vector is defined as supp u = {i ∈ [D] : u_i ≠ 0}. The cardinality of a set , measuring the number of elements of the set, is denoted by . The convolution between two vectors u, for some dimension is denoted equivalently by u ⊗ v = v ⊗ u. We denote by the class of proper, convex and lower-semicontinuous functions from a finite dimensional vector space (e.g., ) to (–∞, +∞] (Combettes & Pesquet 2011). The (convex) indicator function of the set is equal to 0 if and +∞ otherwise. The 2-D discrete delta function, denoted as δ₀ (m, n), is equal to 1 if m = n = 0, and 0 otherwise. The transposition of a matrix Φ reads Φ^T. For any p ≥ 1, the ℓ_p-norm of u is . The Frobenius norm of Φ is given by . The Gaussian distribution of mean 0 and variance σ ² is denoted by .

2. Problem statement

The telescope imaging process can be mathematically modeled as the instrument’s PSF h convolving the true image x, i.e., h ⊗ x. The convolution operation, represented by ⊗, consists of integrating portions of the actual scene weighted by the PSF. In the following, we consider a discrete setting where the convolution integration is represented by a discrete sum.

2.1. Discrete model

Let us consider that, during the transit, the EUV telescope acquires a set of P images containing a celestial body. In order to limit the computation time, the effective Field-of-View (FoV) of each observation is restricted to an image patch centered on the transiting celestial body with a size several times bigger than the body apparent diameter. The effects of this FoV truncation will be discussed in detail in Sections 3 and 6. In our work, all n × n images are represented as N = n ² vectors so that linear transformations of images are identifiable with matrices. Each observed patch is a collection of values gathered in a vector , with j = 1, … , P. The observed values are modeled as a “true” image convolved by the instrument’s PSF (). The PSF is assumed to be spatially and temporally invariant, thus is the same for each observed patch. For simplicity, we consider that all sources of noise present in the observation can be gathered in an average noise model, which is assumed to be additive, white, and Gaussian. The observations are then assumed to be affected by an Additive White Gaussian Noise (AWGN), , with . The acquisition process of the telescope is thus modeled as follows:$${\mathbf{y}}_j=\mathbf{h}\otimes {\mathbf{x}}_j+{\mathbf{n}}_j.$$(2)

This model assumes that the observed images y_j are uncompressed and have been corrected for charge-coupled device (CCD) effects (dark current removal, flat fielding, despiking). In the case of SDO/AIA, this corresponds to level 1 data processing.

In 1-D, the discrete circular convolution of a vector with a filter can always be described by the product of u with a circulant matrix$$\mathbf{\Phi }\left(\mathbf{g}\right)=\left(\begin{array}{llll}{g}_1& g_2& \cdots & g_n\\ g_n& g_1& \cdots & g_{n-1}\\ ⋮& & & ⋮\\ g_2& g_3& \cdots & g_1\end{array}\right)\in {\mathbb{R}}^{n\times n},$$(3)a matrix where every row is a right cyclic shift of the kernel g (Gray 2006). In 2-D and in particular for the model (2), the convolution of x_j by h can still be represented by the multiplication of x_j by matrix that is now block-circulant with circulant blocks, i.e., a matrix made of n × n blocks, each block being a n × n circulant matrix representing the action of one row of the 2-D filter h in (2).

If we gather the P “true” images in a single matrix , the acquisition model can be transformed into the following matrix form:$$\mathbf{Y}=\mathbf{\Phi }\left(\mathbf{h}\right) \mathbf{X}+\mathbf{N},$$(4)with two matrices that similarly gather the P observed images and their noises, respectively.

2.2. Noise estimation

For the sake of simplicity and to keep fast numerical methods, the model (4) implicitly considers an additive noise model. A way to generalize our method to Poisson noise would be to include a VST (Anscombe 1948; Dupé et al. 2009; Shearer et al. 2012) in the acquisition model in (4), both on the observations Y and on the convolution result Φ(h) X (see Sect. 7 for more details).

Under the AWGN assumption, the energy of the residual noise is known to be bounded using the Chernoff-Hoeffding bound (Hoeffding 1963):$${\left|\left|\mathbf{Y} - \mathbf{\Phi }\left(\mathbf{h}\right)\mathbf{X}\right|\right|}_F^2={\left|\left|\mathbf{N}\right|\right|}_F^2<{\epsilon }^2:={\sigma }^2\left(\mathrm{NP}+c\sqrt{\mathrm{NP}}\right),$$(5)with high probability for . A first guess of the noise variance σ ² can be estimated using the Robust Median Estimator () (Donoho & Johnstone 1994), which is based on the assumption that the noise is an additional high-frequency component in the observed signal. More details on this variance estimation are provided in Section 6.2.

Let us note that the variance could also be obtained from the physical noise characteristics using, for instance, dark areas in the images to estimate the dark noise. However, such computations may not consider all the noise sources present in the observations, providing an underestimation of the actual variance.

In this work, we prefer to adopt another strategy. Since Eq. (5) considers an average noise model and since Eq. (4) is an AWGN approximation of the actual data noise corruption, we aim to find an adaptive value of σ that optimizes the AWGN assumption, i.e., which minimizes somehow the proximity of this simple model to the true corruptions. Section 6.2 describes in detail this adaptive strategy and demonstrates its advantage over the fixed choice σ = σ_RME.

3. Blind deconvolution problem

In this section, we aim at reconstructing both X and h from the noisy observations Y. Since the data is corrupted by an AWGN, we can estimate the image and the filter using a least-squares minimization, i.e., by finding the values that minimize the energy of the noise:$${\min }_{\tilde{\mathbf{X}},\tilde{\mathbf{h}}} \frac{1}{2}{\left|\right|\mathbf{Y}-\mathbf{\Phi }\left(\tilde{\mathbf{h}}\right)\tilde{\mathbf{X}}\left|\right|}_F^2,$$(6)with and .

The blind deconvolution problem in (6) is ill-posed and can have infinite possible solutions. By regularizing this problem we aim to reduce the amount of possible solutions to those that are meaningful. In the coming sections, we first describe the prior information and constraints on the image and the filter and then use this information to formulate a regularized blind deconvolution problem.

3.1. Prior information on the image

In the following, we present in detail the available prior information on the image candidate.

(a) Zero disk intensity: Each image of the transit contains the observation of a celestial body, Venus or the Moon. For each observation y_j, we know that the celestial bodies’ EUV emission is zero and therefore can be represented as a black disk of constant radius R and center c_j, both being known from external astronomical observations: (Pesnell et al. 2012) for SDO/AIA and (J.-P. Wuelser private communication) for SECCHI/EUVI. The set of pixels of x_j inside the disk is denoted by Ω_j. In the minimization problem we can set to zero the pixels of the jth image that are inside the set {Ω_j}, i.e., X_ij = (x_j)_i = 0 for all i ∈ Ω_j. This prior information is crucial in the regularization of the blind deconvolution problem as it allows us to remove the issue of interchangeability between the filter and the image. Furthermore, it prevents “oppositely” translated solutions of the image and the filter, a problem occurring because the convolution process is blind to such translations, i.e., , with a translation by .

(b) Analysis-based sparsity model: As commonly done with ill-posed inverse problems associated to image restoration tasks (Starck & Murtagh 1994; Donoho 1995; Starck & Fadili 2009), we regularize our method with a 2-D wavelet prior. Compared to a frequency representation, as obtained by the Fourier transform, a 2-D wavelet transform allows representing images with a multi-resolution scheme (Mallat 2008), i.e., with a linear combination of localized wave atoms, with variable sizes and locations, whose number is much smaller than the initial number of pixels. Most inverse problems in imaging are regularized using a synthesis wavelet prior, i.e., promoting the synthesis of the estimated image with as few “wavelets” as possible (sparsity criterion) (Starck & Murtagh 1994; Donoho 1995). More recently, better inversion results have been obtained using analysis-based sparsity models (Carrillo et al. 2012; Sudhakar et al. 2015), which rather promote sparse projections of the estimated image over a redundant system of wave-functions, i.e., a dictionary. This set can be much larger than an orthonormal basis and hence it can efficiently capture much more different image features.

Following this recent trend, we decide to adopt such an analysis prior that, in opposition to the synthesis framework, also allows us to work directly on the image domain without increasing the dimensionality of the optimization problem. As a dictionary, we use the system associated to the Undecimated Discrete Wavelet Transform (UDWT; Starck et al. 2010). The UDWT can also be seen as the union of all translations of an orthonormal DWT. Conversely to the DWT, this makes the UDWT translation invariant, i.e., it enables an efficient characterization of all image features whatever their location (Starck et al. 2004; Elad et al. 2010), hence providing a better image reconstruction than the traditional DWT.

In this context, we assume that, if we represent the image candidate on a redundant wavelet dictionary made of W vectors in , the wavelet coefficients are sparse, i.e., the coefficient matrix Ψ^TX has few important values and its ℓ₁-norm (computed over all its entries) is expected to be small. However, the wavelet coefficients whose support touches the boundary of the Moon or Venus disk are not sparse. Hence, we do not consider in this prior model those coefficients that are affected by the occulting body. Also, we follow a common practice in the field which removes the (unsparse) scaling coefficients (Mallat 2008) from the ℓ₁-norm computation. The set of detail coefficients that are not affected by the black disk is denoted by Θ, with . See Section 6.1 for further details on how to compute the set Θ. The matrix is the corresponding selection operator that extracts from any coefficient vector only its components indexed by Θ. The global prior information can thus be exploited by promoting a small ℓ₁-norm on the wavelet coefficients belonging to Θ. The rationale of this prior is to enforce noise canceling and avoid artifacts in the reconstructed image.

(c) Image non-negativity: Note that the image corresponds to a measure of a photon emission process, therefore, its pixel values must be non-negative everywhere.

3.2. Constraints on the PSF

Since the filter is not known a priori and it changes from one instrument to another, we are interested in adding soft constraints that are common to most solar EUV telescopes. We are not aiming at reconstructing all of the PSF components but rather at estimating its core, which is responsible for most of the diffused light.

We first consider that, since the PSF corresponds to an observation of a point, it is non-negative everywhere. Additionally, by assuming that the amount of light entering in the instrument is preserved, the ℓ₁-norm of the filter candidate must be equal to one. This is explained by taking the sum over the pixels of one observed patch in (4), which is equal to in a noiseless process. Since the filter is non-negative, taking its ℓ₁-norm equal to one is equivalent to 1^TΦ(h) = 1^T, which results in 1^Ty_j = 1^Tx_j, i.e., the light is preserved. Therefore, the filter candidate must belong to the Probability Simplex (Parikh & Boyd 2014), defined as .

One important assumption in our proposed method is that the filter candidate is of size , for any , and is centered in the spatial origin of the discrete grid. This means that the filter candidate has a limited support inside a set Γ of size (2b + 1)², i.e., the filter has only important values in the central pixels and is negligible beyond the considered support. This allows us to work with a patch (only a part of the observation) and not with the complete FoV of the instrument in order to reduce the computation time.

EUV images have a high dynamic range. Hence, due to long-range effects, a given pixel value can be affected by the value of another high intensity pixel located as far as 100 arcsec away for SDO/AIA (Poduval et al. 2013) and 1500 arcsec away for SECCHI/EUVI (Shearer et al. 2012). In such cases, despite a PSF that can rapidly decay far from its center, the high intensity pixel can induce a long-range effect that is not taken into account by a filter with truncated support. Let us define the complementary set of Γ as Γ^C = [N]\Γ, with |Γ^C| = N – (2b + 1)². We assume that the actual PSF is composed by two different filters: the PSF core with support on Γ, denoted by h_Γ, which accounts for short-range effects, and another one, denoted with support on Γ^C, accounting for long-range effects, i.e., . An accurate estimation of the long-range PSF is out of the scope of this work. Its effect inside the disk of the celestial body is modeled as a constant μ, i.e., (see Appendix A for more details). For a given patch, the acquisition model in (2) can then be approximated as: In this work, we aim at estimating only the PSF core, which is called h hereafter. The effect of the long-range PSF, denoted by μ, is computed in a preprocessing stage by the average value inside the center of several observed transit disks (see Sect. 6.2). This simple estimation is motivated by the presence of a systematic intensity background at the center of each observed disk transit. This one is quite independent of the Sun activity around such disk and its mean intensity (few tens of DNs) is also far above the estimated noise level (few DNs). In a future work we plan to replace this rough evaluation by a joint estimation of the value of μ with the filter and the image (see Sect. 7 for some discussions).

Notice that our method could allow the addition of other convex constraints, e.g., if we know that the filter is sparse or that 0 ≤ h_i ≤ g_i, for some upper bound g_i on h_i such as a specific power-law decay. However, we will not consider this possibility here as we want to stay as agnostic as possible on the properties of the filter to reconstruct.

3.3. Final formulation

From Sections 3.1 and 3.2, we can formulate the following regularized blind deconvolution problem:$$\begin{array}{cc}{\min }_{\tilde{\mathbf{X}},\tilde{\mathbf{h}}}& {\rho \left|\right| {\mathbf{S}}_{\mathrm{\Theta }}{\mathbf{\Psi }}^{\mathrm{T}} \tilde{\mathbf{X}}\left|\right|}_1+\frac{1}{2}\left|\right| \mathbf{Y}- \mathbf{\Phi }\left(\tilde{\mathbf{h}}\right) \tilde{\mathbf{X}}- \mu {\mathbf{1}}_N{\mathbf{1}}_P^{\mathrm{T}}{\left|\right|}_F^2\\ \mathrm{s}.\mathrm{t}. & ({\tilde{\mathbf{x}}}_j{)}_i=0 \mathrm{if} i\in {\mathrm{\Omega }}_j; ({\tilde{\mathbf{x}}}_j{)}_i\ge 0 \mathrm{otherwise}\\ & \tilde{\mathbf{h}}\in \mathcal{PS}; \mathrm{supp} \tilde{\mathbf{h}}=\mathrm{\Gamma } \end{array}$$(7)with , and cancels the long-range part of the filter. The regularization parameter, denoted by ρ, controls the trade-off between the sparsity of the image projection in a wavelet dictionary and the fidelity to the observations. This essential parameter is estimated in Appendix B.2.

It is important to note that the prior information related to the darkness of the occulting body ( if ), which is crucial in the regularization of the blind deconvolution problem, is taken into account in the first constraint in (7).

The problem of estimating the image and the filter in (7) can also be written as$$\left\{{\mathbf{X}}^{\mathrm{*}},{\mathbf{h}}^{\mathrm{*}}\right\}={\mathrm{argmin}}_{\tilde{\mathbf{X}},\tilde{\mathbf{h}}} \mathcal{L}\left(\tilde{\mathbf{X}},\tilde{\mathbf{h}}\right)\hspace{1em}\mathrm{s}.\mathrm{t}.\hspace{1em}\tilde{\mathbf{X}}\in {\mathcal{P}}_0,\hspace{1em} \tilde{\mathbf{h}}\in \mathcal{D},$$(8)where is the objective function, with the modified observations. The convex sets and are defined as follows: ; .

4. Proximal alternating minimization

In this section, we describe the alternating minimization method used for solving (8) and hence estimating the image and the PSF from the noisy observations.

The problem defined in (8) is non-convex with respect to both X and h, but it is convex with respect to X (resp. h) if h (resp. X) is known. Such a problem is usually solved iteratively, by estimating the image and the filter alternately (Almeida & Almeida 2010). The general behavior of this type of algorithm is as follows:$$\left({\mathbf{X}}^{\left(k\right)},{\mathbf{h}}^{\left(k\right)}\right)\to \left({\mathbf{X}}^{(k+1)},{\mathbf{h}}^{\left(k\right)}\right)\to \left({\mathbf{X}}^{(k+1)},{\mathbf{h}}^{(k+1)}\right).$$

It has been shown in Almeida & Almeida (2010) that this algorithm is able to provide fairly good results. However, there are no theoretical guarantees for its convergence. Recently, Attouch et al. (2010) proposed a proximal alternating minimization algorithm where cost-to-move functions are added to the common alternating algorithm. These functions are quadratic costs that penalize variations between two consecutive iterations. They guarantee the algorithm convergence provided some mild conditions are met on the regularity of the objective function and the constraints in (8). The algorithm iterates as follows:$$\begin{array}{c}{\mathbf{X}}^{(k+1)}={\mathrm{argmin}}_{\tilde{\mathbf{X}}} \mathcal{L}\left(\tilde{\mathbf{X}}, {\mathbf{h}}^{\left(k\right)}\right) +\frac{{\lambda }_x^{\left(k\right)}}{2}{\left|\right|\tilde{\mathbf{X}}-{\mathbf{X}}^{\left(k\right)}\left|\right|}_F^2,\\ {\mathbf{h}}^{(k+1)}={\mathrm{argmin}}_{\tilde{\mathbf{h}}} \mathcal{L}\left({\mathbf{X}}^{(k+1)},\tilde{\mathbf{h}}\right)+\frac{{\lambda }_h^{\left(k\right)}}{2}{\left|\right| \tilde{\mathbf{h}}-{\mathbf{h}}^{\left(k\right)}\left|\right|}_2^2,\end{array}$$(9)where λ_x and λ_h are the cost-to-move penalization parameters that control the variations between two consecutive iterations. Section 4.1 describes how these parameters are tuned.

The non-convexity of (8) prevents attaining a global minimizer. However, it can be shown that, since the objective function and the constraints in (8) meet the conditions stated in Attouch et al. (2010), the algorithm (9) converges to a critical point of the problem. Later in this section we describe how to stop the iterations so that the vicinity of a critical point of (8) is reached.

4.1. Cost-to-move parameters

The presence of cost-to-move parameters (λ_x, λ_h) different than zero ensures the convergence of the algorithm (9) to a critical point of (X, h) (Attouch et al. 2010). However, their values are not crucial in the reconstruction results. In the following, we describe the tuning criteria used in this work which is based on the works of Puy & Vandergheynst (2014). The iterations start with high values of λ_x and λ_h, keeping the image and the filter estimations closer to the initial values. When the number of iterations k increases, the filter and the image estimates become more and more accurate and the parameters λ_x and λ_h can be progressively decreased.

4.2. Stopping criteria

It is important to find an automatic criterion for stopping the iterations in (9) when the solution reaches the vicinity of a critical point of (8). A usual stopping criterion is based on the quality of the reconstruction with respect to the Ground Truth image and filter, which are in general unavailable. As proposed by Almeida & Figueiredo (2013b), we use the spectral characteristics of the noise to analyze the quality of the estimation. This allows stopping the minimization when a high-quality solution has been obtained. Let us define each residual image at iteration k as the difference between the observed image and the blurred estimate:$${\mathbf{R}}^{\left(k\right)}=\mathbf{Z}-\mathbf{\Phi }\left({\mathbf{h}}^{(k+1)}\right) {\mathbf{X}}^{(k+1)}=\left({\mathbf{r}}_1^{\left(k\right)},\dots ,{\mathbf{r}}_P^{\left(k\right)}\right).$$(10)

Since the observation model is degraded by an additive white noise, we know that the residual image is spectrally white if the estimation at iteration k has a good quality, otherwise the residual contains structured artifacts. Therefore, the iterations can be stopped when the residual is spectrally white, i.e., when the 2-D autocorrelation of each residual image is approximately the function δ₀(m, n). Let us note that, for the computation of the 2-D autocorrelation function, each residual vector needs to be previously transformed into a matrix of n × n pixels and normalized to zero mean and unit variance.

Considering a (2L + 1) × (2L + 1) window, we compute the distance between the autocorrelation and the Kronecker function δ₀ as follows:$$\mathcal{M}\left({\mathbf{r}}_j^{\left(k\right)}\right)=-{\sum }_{(m,n)=(-L,-L)}^{(L,L)}{\left({\delta }_0(m,n)-C_{{\mathbf{r}}_j{\mathbf{r}}_j}(m,n)\right)}^2=-{\sum }_{\begin{array}{c}\left(m,n\right)=\left(-L,-L\right)\\ (m,n)\ne \left(\mathrm{0,0}\right)\end{array}}^{(L,L)}{\left(C_{{\mathbf{r}}_j{\mathbf{r}}_j}(m,n)\right)}^2,$$(11)which is higher for whiter residuals. As suggested by Almeida & Figueiredo (2013b), in our experiments we have used L = 4. We then average over the measures obtained for each residual image:$$\mathcal{M}\left({\mathbf{R}}^{\left(k\right)}\right)=\frac{1}{P}\sum _{j=1}^P\mathrm{}\mathcal{M}\left({\mathbf{r}}_j^{\left(k\right)}\right).$$(12)

In practice, we observe that has large negative values at the beginning of the iterations when the image and the filter are not properly estimated. Then, the value of increases with k until it reaches a maximum close to zero. Finally, it starts to decrease again after the algorithm has converged to the best estimations for a fixed value of ρ and starts overfitting the noise. Therefore, the iterations can be stopped when the maximum of the whiteness measure () is reached. A similar behavior has also been observed by Almeida & Figueiredo (2013b).

4.3. Numerical reconstruction

The proximal alternating minimization algorithm is summarized in Algorithm 1. This algorithm requires to set the initial values of the image, the filter, and the regularization parameter. The initialization of the image and the filter is discussed in Appendix B.1, and the regularization parameter is tuned iteratively as explained in Appendix B.2.

The first step in Algorithm 1 estimates the image by minimizing a sum of three convex functions with the image constrained to the set . This constraint can be handled by adding the convex indicator function on the set, i.e., . The resulting optimization, containing a sum of two smooth and two non-smooth convex functions, can be solved using the primal-dual algorithm of Chambolle & Pock (2011) summarized in Appendix C.1.

Algorithm 1: Proximal alternating minimization algorithm

Initialize: X⁽¹⁾ = X₀; h⁽¹⁾ = h₀; ; Δ = 0.75; MaxIter = 20

 1: for k = 1 to MaxIter do

     1st step, image estimation :

 2:   

     2nd step, filter estimation :

 3:  

     Parameters update :

 4:  

 5:  

     Compute residuals and whiteness measure :

 6:  

 7:   Compute using (12)

     Stop when residual is spectrally whiter :

 8:   break.

 9: Return and .

The second step in Algorithm 1 estimates the filter by minimizing a sum of two smooth convex functions with the filter constrained to the set . The optimization problem can be solved using the Accelerated Proximal Gradient (APG) method (Parikh & Boyd 2014). This algorithm is able to minimize a sum of a non-smooth and a smooth convex function, which allows us to handle the constraint as a convex indicator function on the set, i.e., . The APG algorithm is briefly described in Appendix C.2.

The algorithms used to solve the minimization problems described above were chosen because they are well suited to solve each problem optimally. Since the cost functions are different in each step, the same algorithm could not be used to optimally find a solution of both problems. The presence of two non-smooth functions in the first step prevents the use of the APG algorithm. For the second step, the algorithm of Chambolle & Pock (2011) could be used to solve the optimization problem, however, this algorithm presented a slower convergence than APG, which is optimal when the optimized cost contains a differentiable function.

Let us note that, in the numerical reconstruction, the convolution operator is implemented using the Fast Fourier Transform (FFT). This operation assumes that the boundaries of the image are periodic, a condition that is far from reality. In actual imaging, no condition can be assumed on the boundary pixels, i.e., they cannot be assumed to be zero, neither periodic, nor reflexive. Not taking care of the boundaries conditions would introduce ringing artifacts. Based on the work of Almeida & Figueiredo (2013a), in the numerical experiments we consider unknown boundary conditions, i.e., the pixels belonging to the image border are not observed. To take this into account, the problem formulation needs to be appropriately modified as discussed in Appendix D.

5. Filter validation through non-blind deconvolution

This section is dedicated to the validation of the filter estimated using Algorithm 1. Since the true emissions of the Moon and Venus are zero, the deconvolution with a correct filter should remove all the celestial bodies’ apparent emissions. Therefore, the validation process consists of taking an image from a transit observation that was not used for filter estimation (i.e., ), deconvolving this image using the estimated filter h, and verifying that, in the reconstructed image , the pixels inside the celestial body disk are zero.

To obtain this reconstructed image, we formulate a least-squares minimization problem regularized using the prior information available on the image, that is, that the image is non-negative and has a sparse representation on a suitable wavelet basis . The proposed non-blind deconvolution method reads as follows:$$\begin{array}{c}{\mathbf{x}}_j^{\mathrm{*}}={\mathrm{argmin}}_{{\tilde{\mathbf{x}}}_j} \rho {\left|\right| {\mathbf{S}}_{\Delta } {\mathbf{\Psi }}^{\mathrm{T}} {\tilde{\mathbf{x}}}_j\left|\right|}_1+\frac{1}{2}{\left|\right| {\mathbf{z}}_j-\mathbf{h}\otimes {\tilde{\mathbf{x}}}_j\left|\right|}_2^2\\ \mathrm{s}.\mathrm{t}. {\tilde{\mathbf{x}}}_j\in \mathcal{P}\end{array}$$(13)where is the selection operator of the set Δ, with , which contains the detail coefficients. The convex set is defined as . The non-blind image estimation is solved using an extension of the primal-dual algorithm of Chambolle & Pock (2011) (see Sect. C.1 for more details). Since this algorithm is able to minimize a sum of three non-smooth convex functions, the constraint can be handled as a convex indicator function. Notice that other algorithms could be used such as the Generalized Forward Backward algorithm (Raguet et al. 2013).

6. Experiments

In this section, we first present some synthetic results that allow us to validate the effectiveness of the proposed blind deconvolution method to recover the correct filter. Later we present experimental results on images taken by the SDO/AIA and SECCHI/EUVI telescopes.

All algorithms were implemented in MATLAB and executed on a 3.2 GHz Intel i5-650 CPU with 3.7 GiB of RAM, running on a 64 bit Ubuntu 14.04 LTS operating system.

6.1. Synthetic data

Three synthetic images are selected to test the reconstruction (P = 3). They are defined on a 256 × 256 pixel grid (N = 256²). To simulate actual images from the celestial transit, we take an image of the Sun (the image observed by SDO/AIA at 00:00 UT on June 6th, 2012) and select three different cutouts of N pixels each. The center of the cutouts is arbitrarily selected such that the corresponding regions are sufficiently distant from each other and correspond to a sample of a full Venus trajectory. In each image, we add a black disk of radius R = 48 pixels, centered at the pixel [129, 129], which represents the celestial body. Figure 2 depicts the images from the simulated transit.

During the solar transit, the position and size of the celestial body are known at all moments. However, when this event is imaged, the coordinates of the center represented in the discrete image grid may contain an error of maximum one pixel. This uncertainty is simulated in our synthetic experiments by considering the sets Ω_j and Θ to be defined using disks of radii r_Ω and r_Θ, respectively, with one pixel difference with respect to the actual object’s radius, i.e., r_Ω = R − 1 and r_Θ = R + 1.

Two kinds of discrete filters are selected and they have a limited support on a 33 × 33 pixel grid (b = 16). The first filter is simulated by an anisotropic Gaussian function with standard deviation of two pixels horizontally and four pixels vertically, rotated by 45° (see Fig. 3A). The second filter is simulated by a X shape (see Fig. 3B), hereafter called the X filter. These functions help to demonstrate the capacity of the proposed method to recover not only diffusion filters such as the anisotropic Gaussian, but also diffraction filters such as the X example. Let us note that, in these synthetic experiments, there is no need to test the long-range assumption on the PSF.

The measurements are simulated according to (2), where additive white Gaussian noise is added in order to simulate a realistic scenario. The noise variance (σ²) is set such that the blurred signal-to-noise ratio BSNR = 10log₁₀ (var(Φ(h)X)/σ²) is equal to 30 dB, which corresponds to a realistic BSNR in the actual observations.

During the experiments, the wavelet transformation used for the sparse image representation is the redundant Daubechies wavelet basis¹ with two vanishing moments and three levels of decomposition (Mallat 2008). Other mother wavelets and more vanishing moments can be considered, however, due to the dictionary redundancy, the choice does not have a significant impact on the reconstruction results. Higher levels can also be considered but the computational time is notably higher and the reconstruction quality does not significantly improve.

As discussed in Section 3.1, the wavelet coefficients whose support touches the boundary of the occulting body are not sparse as they spread over all the scales of the wavelet transform. Thus, we need to define the set Θ containing the coefficients that are not affected by the disk. To determine this set, we first generate a set of images with constant background that contain disks of radius r_Θ centered as the celestial body. All the elements inside the disks are set as random values with a standard uniform distribution (). We then compute the wavelet coefficients of this image on the basis Ψ. Finally, the set Θ is composed by the indices of the detail coefficients that are equal to zero.

The reconstruction quality of X* with respect to the true image X ^GT is measured using the increase in SNR (ISNR) defined as . The reconstruction quality of h* with respect to the true filter h ^GT is measured using the reconstruction SNR (RSNR), where .

Let us first analyze the behavior of the algorithm when we increase the number of observations P used for the reconstruction. Table 1 presents the reconstruction quality of the results using the blind deconvolution problem for the anisotropic Gaussian filter of Figure 3A. The results correspond to an average over four trials and they are presented for P = 1, 2, and 3 observations. For P = 2, 3 we use a warm start by initializing the algorithm with the results obtained for P − 1.

We notice that the reconstruction quality of the filter significantly improves as the number of observations P increases. The whiteness measure presents the same behavior, which indicates that the residual image is whiter when more observations are considered. This is explained by an increase of the available information in the filter reconstruction problem for the same number of unknowns.

Regarding the numerical complexity, the algorithm requires an average of five iterations on the value of ρ and a total time of 1.5 h. When the number of iterations on ρ increases, we observed that the estimated residual energy is closer to the actual noise energy.

Let us now show some reconstruction results of the different filters and how they reproduce the true zero pixel values inside the object. Figure 4A depicts the first noisy observed patch and Figure 4B presents the resulting PSF when reconstructing the anisotropic Gaussian filter for P = 3. We can observe how the algorithm is able to estimate the filter with a high quality, even when no hard constraints are introduced on the filter shape. The reconstruction quality can be further improved if stronger conditions are assumed on the filter.

Figure 4C presents the reconstructed image using the estimated filter from Figure 4B in the non-blind deconvolution of Section 5. Note that the majority of the pixels inside the estimated disk are zero, except for some numerical errors. To quantify this validation, we use as measurement the disk intensity, i.e., the sum of the pixel values inside the disk. We compare the ratio between the disk intensity for the deconvolved image, denoted by , and the disk intensity for the observed patch, denoted by . We obtained a disk intensity ratio of . This shows the effectiveness of the estimated filter to recover the true zero emissions inside the disk.

Figure 5 depicts the results for P = 3 observations using the X filter of Figure 3B. Let us note that the reconstruction quality is lower than in the case of the Gaussian filter because the X filter is more complex and hence it is harder to estimate without extra information on its shape. Nevertheless, since it causes less diffusion on the observation, the image reconstruction quality is better. When comparing the values inside the disk between the observations and the estimated images, we observed that the disk intensity ratio is , which validates the zero emissions inside the disk.

Finally, we also considered a case (not displayed here) where the observation is only corrupted by noise. As expected, the filter estimated by the algorithm is a high-quality Kronecker delta function with RSNR = 66 dB. This result shows the capability of the algorithm to recover highly localized filters of one pixel radius.

6.2. Experimental data

The proposed blind deconvolution method was tested on two experimental sets of data. A first data set corresponds to the observations of the Venus transit on June 5th–6th, 2012 by the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) on board the Solar Dynamics Observatory (SDO). The second data set corresponds to the observations of the Moon transit on February 25th, 2007 by the Extreme Ultraviolet Imager (EUVI; Howard et al. 2008), a part of the Sun Earth Connection Coronal and Heliospheric Investigation (SECCHI) instrument suite on board the STEREO-B spacecraft. These images are compressed using the RICE algorithm (Nightingale 2011), which is lossless (Pence et al. 2009) and hence does not introduce additional errors in the PSF estimation.

In the following, we present for each telescope the results obtained when reconstructing the filter using the blind deconvolution approach. Then, we validate the obtained filters with transit images that were not used for the PSF estimation. Finally, we demonstrate how the obtained filters work when deconvolving non-transit images. All results are compared with previously estimated PSFs.

6.2.1. SDO/AIA – Venus transit

We consider three 4096 × 4096 level 1 images from the transit (P = 3) recorded by the 19.3 nm channel of AIA. The filter is assumed to have a limited support inside a 129 × 129 pixel grid (b = 64), which allows encompassing 99% of the energy of previously estimated PSFs. The presence of a long-range PSF was verified in the observations by analyzing the pixels inside the disk of Venus. To estimate this effect, we considered a total of 10 patches and computed the mean intensity value on a disk of radius 10 pixels inside the disk of Venus. The obtained value μ = 43.3 DN was removed from the observations to estimate the filter using the blind deconvolution scheme.

The radius of the Venus disk is known to be 49 pixels and hence, the disk represents a small area inside the high-resolution image. Since the blind deconvolution takes benefit mainly from the black disk in the transit images, we select a patch of size N = 256² centered on the Venus disk. As explained in Appendix D, our method considers unknown border pixels based on the filter support (completely defined by b after removing the value of μ from the observed patch). Therefore, cropping the observed image does not have any effects on the reconstruction results. Furthermore, considering smaller observations keeps the optimization problem computationally tractable.

As explained in Section 2.2, we use an adaptive noise estimation strategy to optimize the AWGN assumption in (4). For this, we start from σ = σ_RME, the value estimated by the Robust Median Estimator (RME). Then, the blind deconvolution method is performed for different values of σ that are multiples of σ_RME. Finally, we take the value of σ that optimizes the whiteness of the residual as defined in (12), i.e., the residual in (10) only contains white noise without any remaining signal features. The variance is computed as follows (Donoho & Johnstone 1994):$${\sigma }_{\mathrm{RME}}^2=\frac{\mathrm{med}| \mathbf{\alpha } |}{0.6745},$$(14)where α is a vector containing the finest scale wavelet coefficients of the observed vector y_j, i.e., , with the selection operator of the set that contains the finest scale wavelet coefficients. This estimation resulted in σ_RME = 2.81 DN.

Figure 6 shows the resulting PSFs for different values of σ. We can observe that the whitest residual is obtained for σ = 2σ_RME = 5.62 DN and we thus select this value for our non-parametric PSF estimate depicted in Figure 6B. Let us note that, if the noise is underestimated, we reconstruct part of the noise in the PSF and image, and the resulting PSF is too noisy (see Fig. 6A). Oppositely, if the noise is overestimated, the algorithm provides the trivial solution, i.e., the PSF tends to a discrete delta and the image to the observations (see Fig. 6C).

The estimated non-parametric PSF in Figure 6B is compared with the parametric PSF estimated by Poduval et al. (2013), depicted in Figure 7A. This PSF, denoted h_p1, was obtained by fitting a parametric model based on the optical characterization of the telescope. We observe that the obtained non-parametric filter is highly localized in the center of the discrete grid and presents a limited support. Let us note that the observation noise level prevents the estimation of the diffraction peaks present in h_p1. These can only be obtained by constraining the shape of the filter in the reconstruction process, as implicitly done by parametric deconvolution methods.

To solve the blind deconvolution problem, the algorithm requires an average of four iterations on the value of ρ and a total time of 8 h on our computer. In order to reduce the computation time, some functions such as the proximity operators may be computed in parallel, as many of the convex functions on which they are defined are separable (Parikh & Boyd 2014). The computation time can be further reduced if, for instance, the algorithms are implemented in C language instead of MATLAB.

In order to obtain a more accurate PSF estimation containing the diffraction peaks, let us incorporate a parametric PSF inside the acquisition model in (2) by considering a combined parametric/non-parametric PSF defined as (Shearer et al. 2012; Shearer 2013). The parametric part is obtained by considering only the mesh diffraction components in the PSF estimated by Poduval et al. (2013). This modified PSF, denoted h_p2, is depicted in Figure 7B. The non-parametric part () is estimated using the proposed blind deconvolution approach. Figure 7C presents the resulting parametric/non-parametric PSF, i.e., .

We can observe that the resulting parametric/non-parametric PSF is a corrected version of the parametric model from Figure 7B. Its central shape is also more diffused similarly to what is observed in the parametric PSF estimated by Poduval et al. (2013).

The estimated filters from Figures 6B, 7A, and 7C were validated using the non-blind deconvolution formulation from Section 5. Notice that, similarly to what has been done for the blind deconvolution, the value of ρ has been selected adaptively by optimizing the residual whiteness defined in (11). We also observed empirically that the value of ρ must be kept smaller than the one obtained in Algorithm 2. The validation was done on transit images using the modified observations with μ = 43.3 DN. A first validation was performed on the Venus transit image taken by SDO/AIA at 00:02 UT on June 6th, 2012 (see Fig. 8A), a patch that has not been used before for estimating h*. A second validation of the filters was done on the Moon transit image taken by SDO/AIA at 13:00 UT on March 4th, 2011 (see Fig. 9A).

We quantify the apparent disk emissions by summing the pixel values inside the disk. Table 2 displays the disk intensity ratio for the images deconvolved using the following filters: (1) the parametric PSF estimated by Poduval et al. (2013), i.e., h_p1; (2) the non-parametric PSF of Figure 6B, i.e., ; and (3) the parametric/non-parametric PSF of Figure 7C, i.e., .

We observe that, for both transits, the parametric/non-parametric model reaches lower disk apparent emissions. We also note that, compared to the non-parametric filter estimation, the parametric model provides slightly better results in terms of reducing the disk apparent emissions. However, the non-parametric scheme presents the advantage of being generally applicable for any optical instrument without the need of an exact model of the imaging process.

Figures 8 and 9 illustrate the image reconstruction results when using the non-parametric PSF, i.e., . Figures 8A and 9A present the observed patches; Figures 8B and 9B depict the 2-D estimated images; and Figures 8C and 9C present 1-D profiles that allow the observation of the disks of Venus and the Moon, respectively. Due to lack of space, only the non-parametric PSF validation is illustrated.

Let us note that, if the long-range effect μ is not removed from the observations, the parametric PSFs taken with a larger support of 2048 × 2048 pixels are not able to eliminate the offset and, hence, are not able to remove the apparent emissions inside the disk of Venus.

Finally, the estimated filters were also used to deconvolve non-transit images. For this, the non-blind deconvolution formulation from Section 5, using an adaptive ρ, was applied to the modified observations with μ = 43.3 DN. The results are shown for a non-transit image taken by SDO/AIA at 10:00 UT on August 8th, 2011. We selected a portion of the original image of 512 × 512 pixels around the active region (see Fig. 10A). Figure 10B depicts the 2-D estimated image using the non-parametric PSF, i.e., , and Figure 10C shows a 1-D profile along y = 187.

We observe that the non-parametric PSF enhances the image, providing brighter active regions and coronal loops, and darker regions of lower intensity than in the original image. To compare those results with those of the other filters mentioned above, we take the ratio between the deconvolved images and the observation (computing a pixel-by- pixel division), as shown in Figure 11. We notice that the deconvolution resulting from the combined parametric/non-parametric PSF presents a higher correction from the observations than the ones obtained using the other filters. In dark regions, the non-parametric PSF provides results similar to those obtained by the parametric PSF, however, in brighter areas the non-parametric PSF seems to be able to recover more details.

6.2.2. SECCHI/EUVI – Moon transit

We consider three 2048 × 2048 images from the transit (P = 3) recorded by the 17.1 nm channel of EUVI. The images are calibrated with the secchi_prep.pro procedure available within the IDL SolarSoft library. Following the practice described for the SDO/AIA Venus transit, each image was cropped using a 512 × 512 window (N = 512²) centered around the Moon disk. The filter is assumed to have a limited support inside a 129 × 129 pixel grid (b = 64). This allows one to obtain the core of the PSF of around 100 × 100 pixels, as observed by Shearer et al. (2012), and encompass 99% of the energy of previously estimated PSFs.

The long-range PSF is modeled by a constant μ = 12.5 DN, estimated by computing the mean intensity value (over five patches) on a disk of radius 10 pixels inside the Moon disk. The estimated noise variance using the RME provided . Following the same procedure as described before for SDO/AIA, we obtained that the value of σ that allows obtaining the whitest residual is σ = 2 σ_RME = 4.04 DN.

Hereafter, we provide a comparison between the following filters: (1) the parametric PSF given by the euvi_psf.pro procedure of SolarSoft, i.e., , the standard PSF used for SECCHI/EUVI image analysis; (2) the parametric PSF estimated by Shearer et al. (2012), i.e., , and given by the euvi_deconvolve.pro procedure of SolarSoft; (3) the non-parametric PSF, i.e., ; (4) the parametric/non-parametric PSF obtained by incorporating in the acquisition model the parametric PSF given by euvi_psf.pro, i.e., ; and (5) the parametric/non-parametric PSF obtained by incorporating in the acquisition model the parametric PSF given by euvi_deconvolve.pro, i.e., . The core of the parametric and non-parametric filters is presented in Figure 12 and the core of the resulting parametric/non-parametric filters is depicted in Figure 13.

We can observe that while the two parametric PSFs favor both diagonals, the one given by the euvi_deconvolve.pro presents a slight dominance of that at 45°. The estimated non-parametric PSF clearly favors one of the two diagonals. Nevertheless, it presents an orientation at −45°. The source of this difference in orientation is unknown and it should be further investigated. This, however, is beyond the scope of the present paper. Both parametric/non-parametric PSFs present a similar behavior, favoring both diagonals but with a slight dominance of that at −45°.

The filters from Figures 12 and 13 were validated using the non-blind deconvolution described in Section 5. Similarly to what has been done for SDO/AIA, the value of ρ has been selected adaptively by optimizing the value of the residual whiteness. Let us note that, as observed empirically, the value of ρ must be kept smaller than the one obtained in Algorithm 2. The filter validation was performed on the Moon transit image taken by SECCHI/EUVI at 08:02 UT on February 25th, 2007 (see Fig. 14A). Let us note that this patch has not been used before for estimating h* and that the non-blind deconvolution is applied to the modified observation using μ = 12.5 DN. We quantify the apparent Moon emissions by summing the pixel values inside the disk. Table 3 displays the disk intensity ratio for the image deconvolved using the different filters presented above.

We notice that, as for the Venus transit images, when the long-range effect μ is not removed from the observations, the parametric PSFs considered in a larger support of 1024 × 1024 pixels are not able to remove the offset and to recover zero emissions inside the lunar disk. If the long-range effect is considered through the parameter μ, the non-parametric PSF presents a better behavior and is able to reduce more of the disk emissions as compared to the parametric PSFs (see Table 3). We also observe that the parametric/non-parametric PSFs reach lower disk apparent emissions.

In Figure 14 we illustrate the reconstruction results when using the non-parametric PSF, i.e., . Figure 14A presents the observed patch, Figure 14B depicts the 2-D estimated image, and Figure 14C depicts a 1-D profile along y = 257 (to allow the observation of the lunar disk). Due to lack of space, only the non-parametric PSF validation is illustrated.

We can observe that, in Figure 14, part of the off-limb portion of the image is forced to zero. Since this part of the image is fainter than the lunar disk, the proposed non-blind deconvolution method in Section 5 is not able to preserve such low intensities. However, finding a deconvolution method that would preserve these low intensities in the off-limb region (e.g., by selecting other sparsity priors) is beyond the scope of this work.

The estimated filters were also used to deconvolve non-transit images. For this, the non-blind deconvolution formulation from Section 5, using an adaptive ρ, was applied to the modified observations with μ = 12.5 DN. The results are shown for a non-transit image taken by SECCHI/EUVI at 04:02 UT on February 25th, 2007. We selected a portion of the original image of 256 × 256 pixels around the active region (see Fig. 15A). Figure 15B depicts the 2-D estimated image using the non-parametric PSF, i.e., , and Figure 15C shows a 1-D profile along y = 141.

Figure 15 shows that the estimated non-parametric PSF is able to enhance the image and provide more details. Similarly to SDO/AIA, we compare these results with the ones from the other PSFs by taking the ratio between the deconvolved images and the observation (computing a pixel-by-pixel division). Results are depicted in Figure 16. We observe that the non-parametric PSF provides similar results to those obtained by the parametric PSF given by the euvi_psf.pro procedure. These two PSFs are able to provide higher details than the parametric PSF given by the euvi_deconvolve.pro procedure. The deconvolution resulting from the combined parametric/non-parametric PSFs present a higher correction from the observations than the ones obtained using the other three filters. This corresponds to what was observed for the Moon transit images in Table 3.

7. Discussion and perspectives

We proposed a blind deconvolution method that allows one to recover the PSF of a given astronomical instrument provided the latter captures a celestial body transit (as observed in Fig. 1). The proposed non-parametric approach is comparable to previous parametric ones with respect to the quality of the PSF core. It presents, however, some limitations in terms of the considered noise model and the estimation of the PSF long-range effects. Fortunately, the optimization techniques we use in this work are flexible enough to open perspectives of improvement in future works.

7.1. Noise model

The incidence of photon flux on a EUV telescope is converted to digital numbers (DN) through a series of steps, each potentially introducing some noise. The beam of photons impinges the optical system where the PSF acts as a blurring operator. Simultaneously, a spectral selection is performed on the signal before it reaches the CCD detector. The latter has an heterogeneous response across its surface. Finally, the camera electronics convert electrons into DN, adding the read-out noise. The pixels in the resulting image can be modeled as the realization from a random variable Y, whose noise part can be decomposed into additive, Poissonian, and multiplicative degradations. The expectation () and the variance () of Y then verify:$$\mathbb{E}\left[Y\right]=x\hspace{1em}\mathrm{and}\hspace{1em}\mathbb{V}\left[Y\right]={\sigma }^2+\beta x+\alpha x^2,$$(15)where σ is the standard deviation of the additive component and α, β are parameters.

Shearer et al. (2012) handled Poisson corrupted data through a Variance Stabilization Transform (VST) which provides an approximated Gaussian noise distributed data. In order to generalize our AWGN model to models such as (15), a VST could be included in the ℓ₂-fidelity term in (7) and (13), on both the observations and the convolution result (Shearer et al. 2012). This can be done provided that we know the conversion between DN and the photon counts (e.g., from instrument specifications). The algorithm described in Section 4 is adaptable to this stabilized fidelity term through specific proximal operators or gradient descent (Dupé et al. 2009; Anthoine et al. 2012). Notice, however, that such a stabilization cannot be applied only on the observations as a mere preprocessing of the data, while using afterwards other deconvolution methods assuming additive Gaussian noise corruption. Indeed, the VST being a non-linear process of the form for some (Anscombe 1948), it breaks the convolutive model (2) on which those deconvolution methods rely. In other words, the VST of the convolution of h by x_j in (2) is not equivalent to (and hardly approximable by) the convolution of the variance stabilized vectors, which breaks the applicability of those techniques. As an alternative to using a VST, the true Poisson distribution could also be used to design a specific fidelity term, which is based on the negative log-likelihood of the posterior distribution induced by the observation model, i.e., the KL divergence of this distribution (Fish et al. 1995; Anthoine et al. 2012; Prato et al. 2013). However, it is unclear how to adapt the method proposed by Attouch et al. (2010) to the resulting framework.

7.2. Estimation of long-range effects

Our paper aimed at preventing strong assumptions on the shape of the central part of the PSF. We provide a general deconvolution method for different types of telescopes, the generality of our method being of particular interest for telescopes exhibiting optical aberrations. This strategy differs from the one adopted by most parametric PSF estimations. Nonetheless, additional convolutive filter regularizations could be included in our scheme by, for instance, promoting the sparsity (in synthesis or in analysis) of this filter in an appropriate basis (e.g., wavelet basis), or by enforcing a certain decaying law of its amplitude in function of the radial distance. Both kind of priors can be expressed as convex costs (e.g., with a ℓ₁-norm or a weighted ℓ_∞-ball constraint) with closed-form (or simple) proximal operators. These adaptations can thus be integrated in Algorithm 1 with additional efforts for limiting the computational time of the more complex deconvolution procedure. Moreover, such additional priors can help in enlarging the support Γ where the filter is truly estimated, hence reaching an estimation of the long-range PSF.

Regarding the approximation of the long-range PSF impact by a constant, we notice that, instead of estimating the constant μ in a preprocessing step (see Sect. 3.2), in future work we could let μ be a free parameter of the deconvolution problem (7). The value of μ can then be optimized jointly with the image and the filter as this does not break the convexity of the subproblems detailed in Section 4.

Finally, inspired by the works of Shearer et al. (2012) and Shearer (2013), a convolutive combination of a known parametric filter with an unknown non-parametric one was considered in this work to further stabilize the non-convexity of the blind deconvolution problem. However, this construction is limited since, from the convolution theorem, no correction of the parametric PSF can be made in the part of its spectrum where it vanishes. Future works should therefore consider additive correction of the parametric PSF as suggested by Poduval et al. (2013), a modification that Algorithm 1 could also include.

8. Conclusion

We have demonstrated how a non-parametric blind deconvolution technique is able to estimate the core of the PSF of an optical instrument with high quality. The quality of the estimated PSF core is comparable with the one provided by parametric models based on the optical characterization of the imaging instrument. We also demonstrate that, if the parametric PSF is incorporated in the acquisition model, the blind deconvolution approach is able to provide a “corrected” PSF such that most of the apparent emissions inside the disk of a celestial body during a solar transit can be removed. Let us note that non-parametric techniques cannot outperform the accuracy of parametric methods, however in situations where the telescope imaging model cannot be obtained due to some instrument’s properties (e.g., PICARD/SODISM; Meftah et al. 2014), the non-parametric method presents a great advantage. Moreover, the proposed method is not specific to a given instrument but can be applied to any optical instrument provided that we have strong knowledge of some image pixel values and their exact location, such as the information available during the transit of the Moon or a planet. We have also shown that the use of the Proximal Alternating Minimization technique proposed by Attouch et al. (2010) allows to efficiently solve the non-convex problem of blind deconvolution with theoretical convergence guarantees. Furthermore, we show the importance of considering multiple observations for the same filter in order to provide a better conditioning to the filter estimation problem.

Acknowledgments

Part of this work is funded by the AOC SPW Project (SKYWIN – 6894). LJ is supported by the Belgian FRS-FNRS fund. VD acknowledges support from the Belgian Federal Science Policy Office through the ESA-PRODEX program, Grant No. 4000103240. The authors would like to thank Craig DeForest and Jean-François Hochedez for inspiring discussions and Jean-Pierre Wuelser for having kindly provided us Moon trajectory data in SECCHI/EUVI images. Thanks are extended to the referees Richard Frazin and Marco Prato for their comments, which substantially improved the quality of the paper. Also, we would like to thank Richard Frazin who provided the MATLAB code that computes only the mesh diffraction component for the PSF of SDO/AIA. The editor thanks Richard Frazin and Marco Prato for their assistance in evaluating this paper.

---

References

Appendix A

On the approximation of the long-range PSF impact by a constant

In this work, in order to compensate the restriction of the estimated filter to its core, we assume that, for moderate solar intensities, the long-range effect of the PSF can be approximated by a constant inside a limited patch. This assumption is compatible with the observations made on the convolution of a solar image with a parametric filter that only contains the long-range effect, i.e., a filter built by taking a standard parametric PSF containing long-range patterns and setting to zero the central pixels inside a window of size (2b + 1) × (2b + 1) with b = 64.

For SDO/AIA, we use the parametric filter estimated by Grigis et al. (2013). The validation was performed using a 4096 × 4096 level 1 image from the Venus transit recorded by the 19.3 nm channel of AIA at 00:00:01 UT on June 6th, 2012. These results are presented in Figure A.1.

For SECCHI/EUVI, we use the parametric filter given by the euvi_psf.pro procedure of SolarSoft. The validation was performed using a 2048 × 2048 image from the Moon transit recorded by the 17.1 nm channel of EUVI at 14:00:00 UT on February 25th, 2007. The convolution results are presented in Figure A.2.

The resulting low-frequency images in Figure A.1-C and Figure A.2-C validate the use of a constant μ to approximate the long-range effect. Let us note that this is just a first approach to estimate the long-range effect and that more sophisticated methods can be used as explained in Section 7.

Appendix B

Algorithm initialization

Algorithm 1 requires an initial value of the image (X₀), the filter (h₀), and the regularization parameter (ρ₀). In this section, we describe in detail how we proceed with this initialization.

B.1. Image and filter

For the first value of ρ, the alternating algorithm (Algorithm 1) is initialized with the trivial solution, where the image is given by the observations Z and the filter is given by the delta function δ₀. For the subsequent iterations on ρ, the image and the filter are initialized using the value of the previous iteration.

B.2. Regularization parameter

The regularization parameter has an important role in solving the first step in Algorithm 1, since it determines the trade-off between the image regularization and the fidelity to the observations. Several works have studied the choice of this parameter when solving general ill-posed inverse problems. In basis pursuit denoising, Donoho & Johnstone (1994) proposed a universal value for ρ given by . Since this value increases with the amount of samples (N), it tends to provide overly smoothed images. Another state-of-the-art choice, proposed by Donoho & Johnstone (1995), is based on the minimization of the Stein’s Unbiased Risk Estimate (SURE; Stein 1981). This method provides accurate denoising results, however it requires complete knowledge of the degradation model which makes it unsuitable for our blind deconvolution problem. In this work, we use the Bayesian approach proposed by Chang et al. (2000), to estimate the parameter ρ using the noise characteristics and the probability distribution function of the wavelet coefficients:$$\rho =\frac{\sqrt{2} {\sigma }^2}{\tau },$$(B.1)where τ is the standard deviation of the wavelet coefficients of the signal to reconstruct. When we assume that they follow a Laplacian probability distribution with zero mean, the value of τ can be estimated as: , where is the matrix containing the image wavelet coefficients inside the set Θ.

Since the Ground Truth signal X is usually not available, the value of τ in (B.1) is determined using the modified observations Z, i.e., , with . Since this value is higher than the actual ρ, it will lead to over-regularized images. Therefore, we propose to refine the value of ρ by an iterative update based on the discrepancy principle between the actual noise energy in (5) and the energy of the residual image (10). By performing a maximum of five updates of the parameter ρ, we assure an appropriate image regularization. The iterative estimation of ρ is summarized in Algorithm 2.

Algorithm 2: Iterative estimation of ρ

Initialize:

 1: for l = 1 to MaxIter do

    Images and Filter estimation:

 2:  Estimate X^(l+1) and h^(l+1) using Algorithm 1 with X₀ = X^(l), h₀ = h^(l) and ρ₀ = ρ^(l)

    Compute residual and whiteness measure:

 3:  

 4:  Compute using (12)

     Parameter update :

 5:  

 6:  

     Stop when residual is spectrally whiter :

 7:  if then break.

 8: Return and

Appendix C

Numerical reconstruction algorithms

For the sake of completeness, in this section, we briefly describe the proximal algorithms that we use to solve the two subproblems of Algorithm 1. We refer the reader to Chambolle & Pock (2011); Combettes & Pesquet (2011); Parikh & Boyd (2014), for a comprehensive understanding.

Proximal algorithms rely on a key element in convex signal analysis, the proximal operator$${\mathrm{prox}}_{\nu \phi }\mathbf{z}:={\mathrm{argmin}}_{\mathbf{u}\in {\mathbb{R}}^D} \phi \left(\mathbf{u}\right)+\frac{1}{2\nu }\left|\right|\mathbf{u}-\mathbf{z}{\left|\right|}^2,$$(C.1)which is uniquely defined for any function , for some (Combettes & Pesquet 2011). Proximal algorithms allow the minimization of non-smooth functions such as the ℓ₁-norm and the indicator functions present in the image and filter subproblems.

C.1 First subproblem: Image estimation

In the first subproblem, we are interested in finding the image candidate that minimizes$${\mathbf{X}}^{\left(k+1\right)}={\mathrm{argmin}}_{\tilde{\mathbf{X}}\in {\mathbb{R}}^{N\times P}} \rho \hspace{0.5em}{\left|\right|{\mathbf{S}}_{\mathrm{\Theta }} {\mathbf{\Psi }}^{\mathrm{T}} \tilde{\mathbf{X}}\left|\right|}_1+\frac{1}{2}{\left|\right| \mathbf{Z}-\mathbf{\Phi }\left({\mathbf{h}}^{\left(k\right)}\right)\tilde{\mathbf{X}}\left|\right|}_F^2+\frac{{\lambda }_x^{\left(k\right)}}{2}{\left|\right| \tilde{\mathbf{X}}-{\mathbf{X}}^{\left(k\right)}\left|\right|}_F^2+i_{{\mathcal{P}}_0}\left(\tilde{\mathbf{X}}\right),$$(C.2)a problem that contains a sum of four functions belonging to , for some dimensions . For this, we use the Chambolle-Pock (CP) primal-dual algorithm (Chambolle & Pock 2011), which is commonly used for solving the minimization of two functions in Γ⁰. In order to take into account the sum of more than two functions, the algorithm is extended as in González et al. (2014). If we set for , for , for , and for , the CP iterations are:$$\begin{array}{l}{\mathbf{V}}_1^{(t+1)} ={\mathrm{prox}}_{\nu F_1^{\mathrm{*}}}\left({\mathbf{V}}_1^{\left(t\right)}+\nu {\mathbf{S}}_{\mathrm{\Theta }}{\mathbf{\Psi }}^{\mathrm{T}}{\overline{\mathbf{U}}}^{\left(t\right)}\right),\\ {\mathbf{V}}_2^{(t+1)} ={\mathrm{prox}}_{\nu F_2^{\mathrm{*}}}\left({\mathbf{V}}_2^{\left(t\right)}+\nu \mathbf{\Phi }\left({\mathbf{h}}^{\left(k\right)}\right){\overline{\mathbf{U}}}^{\left(t\right)}\right), \\ \begin{array}{l}{\mathbf{V}}_3^{(t+1)} ={\mathrm{prox}}_{\nu F_3^{\mathrm{*}}}\left({\mathbf{V}}_3^{\left(t\right)}+\nu {\overline{\mathbf{U}}}^{\left(t\right)}\right),\\ {\mathbf{U}}^{(t+1)} ={\mathrm{prox}}_{\frac{\alpha }{3}H}\left\{{\mathbf{U}}^{\left(t\right)}-\frac{\beta }{3}\left[\mathbf{\Psi }{\mathbf{S}}_{\mathrm{\Theta }}^{\mathrm{T}}{\mathbf{V}}_1^{(t+1)}+{\left(\mathbf{\Phi }\left({\mathbf{h}}^{\left(k\right)}\right)\right)}^{\mathrm{T}}{\mathbf{V}}_2^{(t+1)}+{\mathbf{V}}_3^{(t+1)}\right]\right\},\\ {\overline{\mathbf{U}}}^{(t+1)} =2 {\mathbf{U}}^{(t+1)}-{\mathbf{U}}^{\left(t\right)},\end{array}\end{array}$$(C.3)with U ^(t) tending to a minimizer X ^(k+1) of (C.2) for t → +∞. The functions are the convex conjugates of functions F_i, with i = {1, 2, 3}, and their proximal operators are computed via the proximal operator of F_i in (C.1) by means of the conjugation property (Combettes & Pesquet 2011): . The step sizes ν and β are adaptively selected using the procedure described in González et al. (2014).

C.2 Second subproblem: filter estimation

In the second subproblem, we are interested in finding the filter candidate that minimizes$${\mathbf{h}}^{(k+1)}={\mathrm{argmin}}_{\tilde{\mathbf{h}}\in {\mathbb{R}}^N} \frac{1}{2}\sum _{j=1}^P\mathrm{}{\left|\right| {\mathbf{z}}_j-\mathbf{\Phi }\left({\mathbf{x}}_j^{(k+1)}\right) \tilde{\mathbf{h}}\left|\right|}_2^2+\frac{{\lambda }_h^{\left(k\right)}}{2}{\left|\right| \tilde{\mathbf{h}}-{\mathbf{h}}^{\left(k\right)}\left|\right|}_2^2+i_{\mathcal{D}}\left(\tilde{\mathbf{h}}\right),$$(C.4)with a block-circulant matrix of kernel . This problem has the following shape$${\mathrm{argmin}}_{\mathbf{u}\in {\mathbb{R}}^N}F\left(\mathbf{u}\right)+G\left(\mathbf{u}\right),$$(C.5)with and . Since is convex and differentiable with gradient , and to belongs to , the problem can be solved using the Accelerated Proximal Gradient (APG) method (Parikh & Boyd 2014). The APG iterations are:$$\begin{array}{l}{\mathbf{v}}^{(t+1)} ={\mathbf{u}}^{\left(t\right)}+{\beta }^{\left(t\right)}\left({\mathbf{u}}^{\left(t\right)}-{\mathbf{u}}^{(t-1)}\right), \\ {\mathbf{u}}^{(t+1)} ={\mathrm{prox}}_{\nu G}\left({\mathbf{v}}^{(t+1)}-\nu \nabla F\left({\mathbf{v}}^{(t+1)}\right)\right),\\ {\beta }^{\left(t\right)}=\frac{t}{t+3},\end{array}$$(C.6)with u ^(t) tending to a minimizer h ^(k+1) of (C.4) for t → +∞. The step size ν is adaptively selected using the line search of Beck & Teboulle (2010).

Appendix D

Convolutions with unknown boundary conditions

As explained in Section 2, our blind deconvolution approach is realized in a restricted FoV of size n × n. For all our computations, as those realized in Algorithm 1, fast convolution methods exploiting the FFT must be performed. Therefore, it is important to properly handle the frontiers of our images and avoid both implicit FFT frontier periodizations and the influence of unobserved image values integrated by the filter extension. Inspired by the work of of Almeida & Figueiredo (2013a), we implement their unknown boundary conditions, where instead of expanding the observation using zero padding as in Anconelli et al. (2006) and Prato et al. (2012), the boundaries are considered unknown in the convolution process and then, by properly selecting the pixels inside the known boundaries, the observed image is obtained. In a nutshell, the unknown boundary conditions method proceeds by considering every convolution in a bigger space obtained by expanding the original one by the radius of the filter in all directions, i.e., by a border of width b. Later, the optimization methods can freely set the values in this border. We only impose that the correct convolution values are those selected in the former domain.

Mathematically, we consider a larger image and a larger filter , with M = N + 2b. The observation is obtained by selecting the pixels inside the boundaries of the convolved image using a selection operator . The problem in (2) can be rewritten as:$$\mathbf{Y}={\mathbf{S}}_N \mathbf{\Phi }\left(\overline{\mathbf{h}}\right) \overline{\mathbf{X}}+\mathbf{N}.$$(D.1)

The regularized blind deconvolution problem in (7) transforms as follows:$$\begin{array}{cc}\begin{array}{c}{\min }_{\tilde{\mathbf{X}},\tilde{\mathbf{h}}} \\ \mathrm{s.t}. \\ \end{array}& \begin{array}{c}\rho {\left|\right|{\mathbf{S}}_{\mathrm{\Theta }} {\mathbf{\Psi }}^{\mathrm{T}} \tilde{\mathbf{X}}\left|\right|}_1+\frac{1}{2}{\left|\right| \mathbf{Z}- {\mathbf{S}}_N \mathbf{\Phi }\left(\tilde{\mathbf{h}}\right) \tilde{\mathbf{X}}\left|\right|}_F^2\\ ({\tilde{\mathbf{x}}}_j{)}_i=0 \mathrm{if} i\in {\mathrm{\Omega }}_j; ({\tilde{\mathbf{x}}}_j{)}_i\ge 0 \mathrm{otherwise}\\ \tilde{\mathbf{h}}\in \mathcal{P}\mathcal{S}; \mathrm{supp} \tilde{\mathbf{h}}\in \mathrm{\Gamma }\end{array}\end{array}$$(D.2)with and . The actual image X^* and filter h^* are obtained by applying the operator S_N to both the estimated image and filter , respectively, i.e., , .

The unknown boundaries are also taken into consideration in the set Θ by excluding from this set the wavelet coefficients that are affected by the image boundaries.

All Tables

All Figures