Technical Article
Nonparametric PSF estimation from celestial transit solar images using blind deconvolution
^{1}
Institute of Information and Communication Technologies, Electronics and Applied Mathematics (ICTEAM), Université catholique de Louvain (UCL), 1348
LouvainlaNeuve, Belgium
^{2}
Royal Observatory of Belgium, Avenue Circulaire 3, 1180
Bruxelles, Belgium
^{*} Corresponding author: adriana.gonzalez@uclouvain.be
Received:
19
December
2014
Accepted:
12
November
2015
Context: Characterization of instrumental effects in astronomical imaging is important in order to extract accurate physical information from the observations. The measured image in a real optical instrument is usually represented by the convolution of an ideal image with a Point Spread Function (PSF). Additionally, the image acquisition process is also contaminated by other sources of noise (readout, photoncounting). The problem of estimating both the PSF and a denoised image is called blind deconvolution and is illposed.
Aims: We propose a blind deconvolution scheme that relies on image regularization. Contrarily to most methods presented in the literature, our method does not assume a parametric model of the PSF and can thus be applied to any telescope.
Methods: Our scheme uses a wavelet analysis prior model on the image and weak assumptions on the PSF. We use observations from a celestial transit, where the occulting body can be assumed to be a black disk. These constraints allow us to retain meaningful solutions for the filter and the image, eliminating trivial, translated, and interchanged solutions. Under an additive Gaussian noise assumption, they also enforce noise canceling and avoid reconstruction artifacts by promoting the whiteness of the residual between the blurred observations and the cleaned data.
Results: Our method is applied to synthetic and experimental data. The PSF is estimated for the SECCHI/EUVI instrument using the 2007 Lunar transit, and for SDO/AIA using the 2012 Venus transit. Results show that the proposed nonparametric blind deconvolution method is able to estimate the core of the PSF with a similar quality to parametric methods proposed in the literature. We also show that, if these parametric estimations are incorporated in the acquisition model, the resulting PSF outperforms both the parametric and nonparametric methods.
Key words: Point Spread Function / Blind deconvolution / Venus transit / Moon transit / SDO/AIA / SECCHI/EUVI / Proximal algorithms / Sparse regularization
© A. González et al., Published by EDP Sciences 2016
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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:(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, outoffocus, 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: readout, 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 socalled “knownPSF 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 illposed. 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 nonparametric PSF. Furthermore, due to the illposedness 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 leastsquares (LS) minimization (Almeida & Almeida 2010). In the presence of Poisson noise, the problem can be formulated using the KullbackLeibler (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., microroughness of the mirrors resulting in longrange 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 preflight 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 longrange 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 STEREOB were studied by Shearer et al. (2012) and by Shearer (2013), where the longrange scattering effect was assumed to follow a parametric piecewise powerlaw 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 semiempirical 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.
Fig. 1.
Accumulation of observations of the Venus transit captured by SDO/AIA on June 5th–6th, 2012. 
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 nonparametric 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/nonparametric PSF provides better results than both parametric and nonparametric 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 nonparametric 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 nonblind 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 nonnegative 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 lowersemicontinuous 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 2D discrete delta function, denoted as δ_{0} (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 σ ^{2} 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 FieldofView (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 ^{2} 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:(2)
This model assumes that the observed images y_{j} are uncompressed and have been corrected for chargecoupled device (CCD) effects (dark current removal, flat fielding, despiking). In the case of SDO/AIA, this corresponds to level 1 data processing.
In 1D, the discrete circular convolution of a vector with a filter can always be described by the product of u with a circulant matrix(3)a matrix where every row is a right cyclic shift of the kernel g (Gray 2006). In 2D 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 blockcirculant 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 2D 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:(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 ChernoffHoeffding bound (Hoeffding 1963):(5)with high probability for . A first guess of the noise variance σ ^{2} can be estimated using the Robust Median Estimator () (Donoho & Johnstone 1994), which is based on the assumption that the noise is an additional highfrequency 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 leastsquares minimization, i.e., by finding the values that minimize the energy of the noise:(6)with and .
The blind deconvolution problem in (6) is illposed 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) Analysisbased sparsity model: As commonly done with illposed inverse problems associated to image restoration tasks (Starck & Murtagh 1994; Donoho 1995; Starck & Fadili 2009), we regularize our method with a 2D wavelet prior. Compared to a frequency representation, as obtained by the Fourier transform, a 2D wavelet transform allows representing images with a multiresolution 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 analysisbased sparsity models (Carrillo et al. 2012; Sudhakar et al. 2015), which rather promote sparse projections of the estimated image over a redundant system of wavefunctions, 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 Ψ^{T}X has few important values and its ℓ_{1}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 ℓ_{1}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 ℓ_{1}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 nonnegativity: Note that the image corresponds to a measure of a photon emission process, therefore, its pixel values must be nonnegative 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 nonnegative everywhere. Additionally, by assuming that the amount of light entering in the instrument is preserved, the ℓ_{1}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 nonnegative, taking its ℓ_{1}norm equal to one is equivalent to 1^{T}Φ(h) = 1^{T}, which results in 1^{T}y_{j} = 1^{T}x_{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)^{2}, 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 longrange 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 longrange 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)^{2}. We assume that the actual PSF is composed by two different filters: the PSF core with support on Γ, denoted by h_{Γ}, which accounts for shortrange effects, and another one, denoted with support on Γ^{C}, accounting for longrange effects, i.e., . An accurate estimation of the longrange 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 longrange 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 powerlaw 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:(7)with , and cancels the longrange part of the filter. The regularization parameter, denoted by ρ, controls the tradeoff 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(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 nonconvex 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:
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 costtomove 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:(9)where λ_{x} and λ_{h} are the costtomove penalization parameters that control the variations between two consecutive iterations. Section 4.1 describes how these parameters are tuned.
The nonconvexity 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. Costtomove parameters
The presence of costtomove 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 highquality solution has been obtained. Let us define each residual image at iteration k as the difference between the observed image and the blurred estimate:(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 2D autocorrelation of each residual image is approximately the function δ_{0}(m, n). Let us note that, for the computation of the 2D 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 δ_{0} as follows:(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:(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 nonsmooth convex functions, can be solved using the primaldual algorithm of Chambolle & Pock (2011) summarized in Appendix C.1.
Algorithm 1: Proximal alternating minimization algorithm
Initialize: X^{(1)} = X_{0}; h^{(1)} = h_{0}; ; Δ = 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 nonsmooth 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 nonsmooth 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 nonblind 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 leastsquares minimization problem regularized using the prior information available on the image, that is, that the image is nonnegative and has a sparse representation on a suitable wavelet basis . The proposed nonblind deconvolution method reads as follows:(13)where is the selection operator of the set Δ, with , which contains the detail coefficients. The convex set is defined as . The nonblind image estimation is solved using an extension of the primaldual algorithm of Chambolle & Pock (2011) (see Sect. C.1 for more details). Since this algorithm is able to minimize a sum of three nonsmooth 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 i5650 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^{2}). 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.
Fig. 2.
Realistic images: (A) first image x_{1}, (B) second image x_{2}, and (C) third image x_{3} 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 longrange assumption on the PSF.
Fig. 3.
Realistic filters in a (2b + 1) × (2b + 1) window with b = 16. (A) Anisotropic Gaussian Filter. (B) X Filter. 
The measurements are simulated according to (2), where additive white Gaussian noise is added in order to simulate a realistic scenario. The noise variance (σ^{2}) is set such that the blurred signaltonoise ratio BSNR = 10log_{10} (var(Φ(h)X)/σ^{2}) 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^{1} 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.
Reconstruction quality for the different number of simulated transit observations using the synthetic anisotropic Gaussian filter.
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.
Fig. 4.
Results for the anisotropic Gaussian filter and P = 3. (A) Noisy observation of image 1 (y_{1}). (B) Filter reconstruction (h*) using blind deconvolution, RSNR = 14.05 dB. (C) Image Reconstruction () using nonblind deconvolution, ISNR = 1.83 dB. 
Figure 4C presents the reconstructed image using the estimated filter from Figure 4B in the nonblind 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.
Fig. 5.
Results for the X filter and P = 3. (A) Noisy observation of image 1 (y_{1}). (B) Filter reconstruction (h*) using blind deconvolution, RSNR = 7.52 dB. (C) Image reconstruction () using nonblind deconvolution, ISNR = 2.55 dB. 
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 highquality 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 STEREOB 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 nontransit 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 longrange 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 highresolution image. Since the blind deconvolution takes benefit mainly from the black disk in the transit images, we select a patch of size N = 256^{2} 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):(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 nonparametric 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).
Fig. 6.
Logarithm of the nonparametric filters estimated for the SDO/AIA telescope taken inside the set Γ. The filters are calculated for several values of σ providing different residual whiteness (measured by ). (A) σ = σ _{RME}, , (B) σ = 2σ_{RME}, , (C) σ = 3σ_{RME}, . 
The estimated nonparametric 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 nonparametric 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.
Fig. 7.
Logarithm of filters for the SDO/AIA telescope taken inside the set Γ. (A) Parametric PSF estimated by Poduval et al. (2013), i.e., h_{p1}. (B) Diffraction mesh pattern modeled using the parameters estimated by Poduval et al. (2013), i.e., h_{p2}. (C) Combined parametric/nonparametric filters using the parametric PSF shown on the center, i.e., . 
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/nonparametric 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 nonparametric part () is estimated using the proposed blind deconvolution approach. Figure 7C presents the resulting parametric/nonparametric PSF, i.e., .
We can observe that the resulting parametric/nonparametric 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 nonblind 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).
Fig. 8.
Image reconstruction results for the Venus transit observed by SDO/AIA at 00:02 UT on June 6th, 2012. Results are shown for the nonparametric filter, i.e., . (A) Observed image, (B) 2D reconstruction result, and (C, best viewed in color) 1D profile along y = 129. Figures (A) and (B) contain a green line indicating the 1D profiles shown on figure (C). 
Fig. 9.
Image reconstruction results for the Moon transit observed by SDO/AIA at 13:00 UT on March 4th, 2011. Results are shown for the nonparametric filter, i.e., . (A) Observed image, (B) 2D reconstruction result, and (C, best viewed in color) 1D profile along y = 186. Figures (A) and (B) contain a green line indicating the 1D profiles shown on figure (C). 
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 nonparametric PSF of Figure 6B, i.e., ; and (3) the parametric/nonparametric PSF of Figure 7C, i.e., .
Disk intensity ratio for the different filters.
We observe that, for both transits, the parametric/nonparametric model reaches lower disk apparent emissions. We also note that, compared to the nonparametric filter estimation, the parametric model provides slightly better results in terms of reducing the disk apparent emissions. However, the nonparametric 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 nonparametric PSF, i.e., . Figures 8A and 9A present the observed patches; Figures 8B and 9B depict the 2D estimated images; and Figures 8C and 9C present 1D profiles that allow the observation of the disks of Venus and the Moon, respectively. Due to lack of space, only the nonparametric PSF validation is illustrated.
Let us note that, if the longrange 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 nontransit images. For this, the nonblind deconvolution formulation from Section 5, using an adaptive ρ, was applied to the modified observations with μ = 43.3 DN. The results are shown for a nontransit 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 2D estimated image using the nonparametric PSF, i.e., , and Figure 10C shows a 1D profile along y = 187.
Fig. 10.
Reconstruction results for a SDO/AIA image containing an active region. The image has been captured at 10:00 UT on August 8th, 2011. Results are shown for the nonparametric filter, i.e., . (A) Observed image, (B) 2D reconstruction result, and (C, best viewed in color) 1D profile along y = 187. Figures (A) and (B) contain a green line indicating the 1D profiles shown on figure (C). 
We observe that the nonparametric 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 pixelby pixel division), as shown in Figure 11. We notice that the deconvolution resulting from the combined parametric/nonparametric PSF presents a higher correction from the observations than the ones obtained using the other filters. In dark regions, the nonparametric PSF provides results similar to those obtained by the parametric PSF, however, in brighter areas the nonparametric PSF seems to be able to recover more details.
Fig. 11.
Reconstruction results for the image in Figure 10A. The figure (best viewed in color) shows the ratio (computed using a pixel by pixel division) between the deconvolved and the observed images along y = 187, for the different considered filters: , , and . The horizontal axis is in correspondence with the horizontal axis of Figure 10C. 
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^{2}) 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 longrange 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 nonparametric PSF, i.e., ; (4) the parametric/nonparametric PSF obtained by incorporating in the acquisition model the parametric PSF given by euvi_psf.pro, i.e., ; and (5) the parametric/nonparametric PSF obtained by incorporating in the acquisition model the parametric PSF given by euvi_deconvolve.pro, i.e., . The core of the parametric and nonparametric filters is presented in Figure 12 and the core of the resulting parametric/nonparametric filters is depicted in Figure 13.
Fig. 12.
Logarithm of the filters for the SECCHI/EUVI telescope taken inside the set Γ. (A) Parametric PSF given by the euvi_psf.pro procedure of SolarSoft, i.e., . (B) Parametric PSF given by the euvi_deconvolve.pro procedure of SolarSoft, i.e., . (C) Nonparametric PSF, i.e., . 
Fig. 13.
Logarithm of the combined parametric/nonparametric filters for the SECCHI/EUVI telescope taken inside the set Γ. (A) Using the parametric PSF given by the euvi_psf.pro procedure of SolarSoft, i.e., . (B) Using the parametric PSF given by the euvi_deconvolve.pro procedure of SolarSoft, i.e., . 
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 nonparametric 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/nonparametric 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 nonblind 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 nonblind 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.
Fig. 14.
Image reconstruction results for the Moon transit observed by SECCHI/EUVI at 08:02 UT on February 25th, 2007. Results are shown for the nonparametric filter, i.e., . (A) Observed image, (B) 2D reconstruction result, and (C, best viewed in color) 1D profile along y = 257. Figures (A) and (B) contain a green line indicating the 1D profiles shown on figure (C). 
Moon disk intensity ratio for the different filters.
We notice that, as for the Venus transit images, when the longrange 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 longrange effect is considered through the parameter μ, the nonparametric 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/nonparametric PSFs reach lower disk apparent emissions.
In Figure 14 we illustrate the reconstruction results when using the nonparametric PSF, i.e., . Figure 14A presents the observed patch, Figure 14B depicts the 2D estimated image, and Figure 14C depicts a 1D profile along y = 257 (to allow the observation of the lunar disk). Due to lack of space, only the nonparametric PSF validation is illustrated.
We can observe that, in Figure 14, part of the offlimb portion of the image is forced to zero. Since this part of the image is fainter than the lunar disk, the proposed nonblind 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 offlimb region (e.g., by selecting other sparsity priors) is beyond the scope of this work.
The estimated filters were also used to deconvolve nontransit images. For this, the nonblind deconvolution formulation from Section 5, using an adaptive ρ, was applied to the modified observations with μ = 12.5 DN. The results are shown for a nontransit 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 2D estimated image using the nonparametric PSF, i.e., , and Figure 15C shows a 1D profile along y = 141.
Fig. 15.
Reconstruction results for a SECCHI/EUVI image containing an active region. The image has been captured at 04:02 UT on February 25th, 2007. Results are shown for the nonparametric filter, i.e., . (A) Observed image, (B) 2D reconstruction result, and (C) best viewed in color) 1D profile along y = 141. Figures (A) and (B) contain a green line indicating the 1D profiles shown on the right. 
Figure 15 shows that the estimated nonparametric 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 pixelbypixel division). Results are depicted in Figure 16. We observe that the nonparametric 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/nonparametric 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.
Fig. 16.
Reconstruction results for the image in Figure 15A. The figure shows the ratio (computed using a pixel by pixel division) between the deconvolved and the observed images along y = 141. (A, best viewed in color) Ratio for filters , , and . (B, best viewed in color) Ratio for filters and . The horizontal axis is in correspondence with the horizontal axis of Figure 15C. 
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 nonparametric 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 longrange 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 readout 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:(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 ℓ_{2}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 nonlinear 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 loglikelihood 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 longrange 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 ℓ_{1}norm or a weighted ℓ_{∞}ball constraint) with closedform (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 longrange PSF.
Regarding the approximation of the longrange 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 nonparametric one was considered in this work to further stabilize the nonconvexity 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 nonparametric 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 nonparametric 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 nonparametric 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 nonconvex 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 FRSFNRS fund. VD acknowledges support from the Belgian Federal Science Policy Office through the ESAPRODEX program, Grant No. 4000103240. The authors would like to thank Craig DeForest and JeanFrançois Hochedez for inspiring discussions and JeanPierre 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.
The wavelet operator and the other operators used in this work were numerically implemented using the SPARCO framework (Berg et al. 2007).
References
 Almeida, M., and L. Almeida. Blind and semiblind deblurring of natural images. IEEE Trans. Image Process., 19 (1), 36–52, 2010, DOI: 10.1109/TIP.2009.2031231. [CrossRef] (In the text)
 Almeida, M., and M. Figueiredo. Deconvolving images with unknown boundaries using the alternating direction method of multipliers. IEEE Trans. Image Process., 22 (8), 3074–3086, 2013a, DOI: 10.1109/TIP.2013.2258354. [CrossRef] (In the text)
 Almeida, M., and M. Figueiredo. Parameter estimation for blind and nonblind deblurring using residual whiteness measures. IEEE Trans. Image Process., 22 (7), 2751–2763, 2013b, DOI: 10.1109/TIP.2013.2257810. [CrossRef] (In the text)
 Anconelli, B., M. Bertero, P. Boccacci, M. Carbillet, and H. Lanteri. Reduction of boundary effects in multiple image deconvolution with an application to LBT LINCNIRVANA. A&A, 448 (3), 1217–1224, 2006, DOI: 10.1051/00046361:20053848. [NASA ADS] [CrossRef] [EDP Sciences]
 Anscombe, F.J. The transformation of Poisson, binomial and negativebinomial data. Biometrika, 35 (3–4), 246–254, 1948, DOI: 10.1093/biomet/35.34.246. [CrossRef] (In the text)
 Anthoine, S., J.F. Aujol, Y. Boursier, and C. Mélot. Some proximal methods for Poisson intensity CBCT and PET. Inverse Prob. Imaging, 6 (4), 565–598, 2012, DOI: 10.3934/ipi.2012.6.565. [CrossRef] (In the text)
 Attouch, H., J. Bolte, P. Redont, and A. Soubeyran. Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the KurdykaŁojasiewicz inequality. Math. Oper. Res., 35 (2), 438–457, 2010, DOI: 10.1287/moor.1100.0449. [CrossRef] [MathSciNet] (In the text)
 Ayers, G.R., and J.C. Dainty. Iterative blind deconvolution method and its applications. Opt. Lett., 3 (7), 547–549, 1988, DOI: 10.1364/OL.13.000547. [NASA ADS] [CrossRef] [PubMed] (In the text)
 Babacan, S., R. Molina, and A. Katsaggelos. Variational Bayesian blind deconvolution using a total variation prior. IEEE Trans. Image Process., 18 (1), 12–26, 2009, DOI: 10.1109/TIP.2008.2007354. [CrossRef] (In the text)
 Beck, A., and M. Teboulle. Gradientbased algorithms with applications to signal recovery problems. In: D., Palomar, and Y. Eldar, Editors. Convex Optimization in Signal Processing and Communications, Cambridge University Press, 42–88, 2010.
 Berg, E.V., M.P. Friedlander, G. Hennenfent, F. Herrmann, R. Saab, and Ö. Yılmaz. Sparco: A testing framework for sparse reconstruction, Tech. Rep. TR200720, Dept. Computer Science, University of British Columbia, Vancouver, 2007.
 Carrillo, R.E., J.D. McEwen, and Y. Wiaux. Sparsity Averaging Reweighted Analysis (SARA): a novel algorithm for radiointerferometric imaging. Mon. Not. R. Astron. Soc., 426 (2), 1223–1234, 2012, DOI: 10.1111/j.13652966.2012.21605.x. [NASA ADS] [CrossRef] [MathSciNet] (In the text)
 Chambolle, A., and T. Pock. A firstorder primaldual algorithm for convex problems with applications to imaging. J. Math. Imaging Vision, 40 (1), 120–145, 2011, DOI: 10.1007/s1085101002511. [CrossRef] [MathSciNet] (In the text)
 Chang, S., B. Yu, and M. Vetterli. Adaptive wavelet thresholding for image denoising and compression. IEEE Trans. Image Process., 9 (9), 1532–1546, 2000, DOI: 10.1109/83.862633. [NASA ADS] [CrossRef]
 Combettes, P.L., and J.C. Pesquet. Proximal splitting methods in signal processing. In: H.H., Bauschke, R.S. Burachik, P.L. Combettes, V. Elser, D.R. Luke, and H. Wolkowicz, Editors. Fixedpoint algorithms for inverse problems in science and engineering, vol. 49 of Springer optimization and its applications, Springer, New York, 185–212, ISBN: 9781441995681, 2011, DOI: 10.1007/9781441995698_10. [CrossRef] [MathSciNet] (In the text)
 DeForest, C.E., P.C.H. Martens, and M.J. WillsDavey. Solar coronal structure and stray light in TRACE. Astrophys. J., 690, 1264–1271, 2009, DOI: 10.1088/0004637X/690/2/1264. [NASA ADS] [CrossRef] (In the text)
 Dong, W., H. Feng, Z. Xu, and Q. Li. A piecewise local regularized RichardsonLucy algorithm for remote sensing image deconvolution. Opt. Laser Technol., 43 (5), 926–933, 2011, DOI: 10.1016/j.optlastec.2010.12.012. [CrossRef] (In the text)
 Donoho, D. Denoising by softthresholding. IEEE Trans. Inf. Theory, 41 (3), 613–627, 1995, DOI: 10.1109/18.382009. [CrossRef] [MathSciNet] (In the text)
 Donoho, D.L., and I.M. Johnstone. Adapting to unknown smoothness via wavelet shrinkage. J. Amer. Statist. Assoc., 90 (432), 1200–1224, 1995, DOI: 10.1080/01621459.1995.10476626. [CrossRef] [MathSciNet]
 Donoho, D.L., and J.M. Johnstone. Ideal spatial adaptation by wavelet shrinkage, Biometrika, 81 (3), 425–455, 1994, DOI: 10.1093/biomet/81.3.425. [CrossRef] [MathSciNet] (In the text)
 Dupé, F.X., J. Fadili, and J.L. Starck. A proximal iteration for deconvolving Poisson noisy images using sparse representations. IEEE Trans. Image Process., 18 (2), 310–321, 2009, DOI: 10.1109/TIP.2008.2008223. [NASA ADS] [CrossRef] (In the text)
 Elad, M., M. Figueiredo, and Y. Ma. On the role of sparse and redundant representations in image processing. Proc. IEEE, 98 (6), 972–982, 2010, DOI: 10.1109/JPROC.2009.2037655. [CrossRef] (In the text)
 Fergus, R., B. Singh, A. Hertzmann, S.T. Roweis, and W.T. Freeman. Removing camera shake from a single photograph. ACM Trans. Graph., 25 (3), 787–794, 2006, DOI: 10.1145/1141911.1141956. [CrossRef] (In the text)
 Fish, D.A., J.G. Walker, A.M. Brinicombe, and E.R. Pike. Blind deconvolution by means of the RichardsonLucy algorithm. J. Opt. Soc. Am. A, 12 (1), 58–65, 1995, DOI: 10.1364/JOSAA.12.000058. [CrossRef] (In the text)
 Gburek, S., J. Sylwester, and P. Martens. The trace telescope point spread function for the 171 Å filter. Sol. Phys., 239, 531–548, 2006, DOI: 10.1007/s1120700619940. [NASA ADS] [CrossRef] (In the text)
 González, A., L. Jacques, C.D. Vleeschouwer, and P. Antoine. Compressive optical deflectometric tomography: a constrained totalvariation minimization approach. Inverse Prob. Imaging, 8 (2), 421–457, 2014, DOI: 10.3934/ipi.2014.8.421. [CrossRef]
 Gray, R.M. Toeplitz and circulant matrices: a review. Found. Trends Commun. Inf. Theory, 2 (3), 155–239, 2006, DOI: 10.1561/0100000006. [CrossRef] (In the text)
 Grigis, P., S. Yingna, and M. Weber. AIA PSF characterization and image deconvolution. Tech. Rep., AIA team, 2013. (In the text)
 Hoeffding, W. Probability inequalities for sums of bounded random variables. J. Amer. Statist. Assoc., 58 (301), 13–30, 1963, DOI: 10.1080/01621459.1963.10500830. [CrossRef] [MathSciNet] (In the text)
 Howard, R.A., J.D. Moses, A. Vourlidas, J.S. Newmark, D.G. Socker, et al. Sun Earth Connection Coronal and Heliospheric Investigation (SECCHI). Space Sci. Rev., 136, 67–115, 2008, DOI: 10.1007/s1121400893414. [NASA ADS] [CrossRef] (In the text)
 Jefferies, S., K. Schulze, C. Matson, K. Stoltenberg, and E.K. Hege. Blind deconvolution in optical diffusion tomography. Opt. Express, 10 (1), 46–53, 2002, DOI: 10.1364/OE.10.000046. [CrossRef] (In the text)
 Lemen, J.R., A.M. Title, D.J. Akin, P.F. Boerner, and C. Chou. The Atmospheric Imaging Assembly (AIA) on the Solar Dynamics Observatory (SDO). Sol. Phys., 275 (1–2), 17–40, 2012, DOI: 10.1007/s1120701197768. [NASA ADS] [CrossRef] (In the text)
 Mallat, S. Wavelet Tour of Signal Processing, Third Edition: The Sparse Way, 3rd edn., Academic Press, ISBN: 0123743702, 9780123743701, 2008. (In the text)
 Martens, P.C., L.W. Acton, and J.R. Lemen. The point spread function of the soft Xray telescope aboard YOHKOH. Sol. Phys., 157 (1–2), 141–168, 1995, DOI: 10.1007/BF00680614. [CrossRef] (In the text)
 Meftah, M., J.F. Hochedez, A. Irbah, A. Hauchecorne, P. Boumier, et al. Picard SODISM, a space telescope to study the sun from the middle ultraviolet to the near infrared. Sol. Phys., 289 (3), 1043–1076, 2014, DOI: 10.1007/s112070130373x. [NASA ADS] [CrossRef] (In the text)
 Nightingale, R.W. AIA/SDO FITS keywords for scientific usage and data processing at levels 0, 1.0, and 1.5. Tech. Rep., LockheedMartin Solar and Astrophysics Laboratory (LMSAL), 2011. (In the text)
 Oliveira, J.A.P., M.A.T. Figueiredo, and J.M. BioucasDias. Blind estimation of motion blur parameters for image deconvolution. In: J., Mart, J. Bened, A. Mendona, and J. Serrat, Editors. Pattern Recognition and Image Analysis, vol. 4478 of Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, 604–611, ISBN: 9783540728481, 2007, DOI: 10.1007/9783540728498_76. (In the text)
 Parikh, N., and S. Boyd. Proximal algorithms. Found. Trends Optim., 1 (3), 127–239, 2014, DOI: 10.1561/2400000003. [CrossRef] [MathSciNet] (In the text)
 Pence, W., R. Seaman, and R. White. A comparison of lossless image compression methods and the effects of noise. In: Astronomical Data Analysis Software and Systems XVIII, ASP Conference Series, Vol. 411, 25, 2009. (In the text)
 Pesnell, W.D., B.J. Thompson, and P.C. Chamberlin. The Solar Dynamics Observatory (SDO). Sol. Phys., 275 (1–2), 3–15, 2012, DOI: 10.1007/s1120701198413. [NASA ADS] [CrossRef] (In the text)
 Poduval, B., C.E. DeForest, J.T. Schmelz, and S. Pathak. Pointspread functions for the extremeultraviolet channels of SDO/AIA Telescopes. Astrophys. J., 765 (2), 144, 2013, DOI: 10.1088/0004637X/765/2/144. [NASA ADS] [CrossRef] (In the text)
 Prato, M., A.L. Camera, S. Bonettini, and M. Bertero. A convergent blind deconvolution method for postadaptiveoptics astronomical imaging. Inverse Prob., 29 (6), 065017, 2013, , DOI: 02665611296065017. [CrossRef] (In the text)
 Prato, M., R. Cavicchioli, L. Zanni, P. Boccacci, and M. Bertero. Efficient deconvolution methods for astronomical imaging: algorithms and IDLGPU codes. A&A, 539, A133, 2012, DOI: 10.1051/00046361/201118681. [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Puy, G., and P. Vandergheynst. Robust image reconstruction from multiview measurements. SIAM J. Imaging Sci., 7 (1), 128–156, 2014, DOI: 10.1137/120902586. [CrossRef] (In the text)
 Raguet, H., J. Fadili, and G. Peyré. A generalized forwardbackward splitting. SIAM J. Imaging Sci., 6 (3), 1199–1226, 2013, DOI: 10.1137/120872802. [CrossRef] [MathSciNet] (In the text)
 Seaton, D.B., A. Degroof, P. Shearer, D. Berghmans, and B. Nicula. SWAP observations of the longterm, largescale evolution of the extremeultraviolet solar corona. Astrophys. J, 777 (1), 72, 2013, DOI: 10.1088/0004637X/777/1/72. [CrossRef] (In the text)
 Shearer, P., R.A. Frazin, A.O. Hero III, and A.C. Gilbert. The first stray light corrected extreme ultraviolet images of solar coronal holes. Astrophys. J., 749, L8, 2012, DOI: 10.1088/20418205/749/1/L8. [CrossRef] (In the text)
 Shearer, P.R. Separable inverse problems, blind deconvolution, and stray light correction for extreme ultraviolet solar images. Ph.D. thesis, The University of Michigan, 2013. (In the text)
 Starck, J.L., M. Elad, D. Donoho. Redundant multiscale transforms and their application for morphological component separation. Adv. Imaging Electron Phys., 132, 287–348, 2004, DOI: 10.1016/S10765670(04)320069. [CrossRef] (In the text)
 Starck, J.L. and M. Fadili. An overview of inverse problem regularization using sparsity. Image Processing (ICIP), 2009 16th IEEE International Conference on, Cairo, 1453–1456, 2009, DOI: 10.1109/ICIP.2009.5414556. (In the text)
 Starck, J.L., and F. Murtagh. Image restoration with noise suppression using the wavelet transform. A&A, 288 (1), 342–348, 1994. (In the text)
 Starck, J.L., F. Murtagh, and J. Fadili. Sparse image and signal processing: wavelets, curvelets, morphological diversity, Cambridge University Press, New York, NY, USA, ISBN: 0521119138, 9780521119139, 2010. [CrossRef] (In the text)
 Stein, C.M. Estimation of the mean of a multivariate normal distribution. Ann. Stat., 9 (6), 1135–1151, 1981, DOI: 10.1214/aos/1176345632. [CrossRef] [MathSciNet]
 Sudhakar, P., L. Jacques, X. Dubois, P. Antoine, and L. Joannes. Compressive imaging and characterization of sparse light deflection maps. SIAM J. Imaging Sci., 8 (3), 1824–1856, 2015, DOI: 10.1137/140974341. [CrossRef] (In the text)
Appendix A
On the approximation of the longrange 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 longrange 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 longrange effect, i.e., a filter built by taking a standard parametric PSF containing longrange 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.
Fig. A.1.
Validating the approximation of the longrange PSF impact by a constant inside the disk of Venus. (A) Logarithm of the filter used for the validation: parametric filter estimated by Grigis et al. (2013) with the central pixels set to zero inside a window of size (2b + 1) × (2b + 1), with b = 64. (B) Observed image inside a window containing the disk of Venus. (C) Convolved image inside a window containing the disk of Venus. The intensity within the window is almost constant. 
Fig. A.2.
Validating the approximation of the longrange PSF impact by a constant inside the lunar disk. (A) Logarithm of the filter used for the validation: parametric filter given by the euvi_psf.pro procedure of SolarSoft with the central pixels set to zero inside a window of size (2b + 1) × (2b + 1), with b = 64. (B) Observed image inside a window containing the lunar disk. (C) Convolved image inside a window containing the lunar disk. The window intensity is almost constant. 
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 lowfrequency images in Figure A.1C and Figure A.2C validate the use of a constant μ to approximate the longrange effect. Let us note that this is just a first approach to estimate the longrange 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_{0}), the filter (h_{0}), and the regularization parameter (ρ_{0}). 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 δ_{0}. 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 tradeoff between the image regularization and the fidelity to the observations. Several works have studied the choice of this parameter when solving general illposed 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 stateoftheart 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:(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 overregularized 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_{0} = X^{(l)}, h_{0} = h^{(l)} and ρ_{0} = ρ^{(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(C.1)which is uniquely defined for any function , for some (Combettes & Pesquet 2011). Proximal algorithms allow the minimization of nonsmooth functions such as the ℓ_{1}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(C.2)a problem that contains a sum of four functions belonging to , for some dimensions . For this, we use the ChambollePock (CP) primaldual algorithm (Chambolle & Pock 2011), which is commonly used for solving the minimization of two functions in Γ^{0}. 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:(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(C.4)with a blockcirculant matrix of kernel . This problem has the following shape(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:(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:(D.1)
The regularized blind deconvolution problem in (7) transforms as follows:(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.
Cite this article as: González A, Delouille V & Jacques L. Nonparametric PSF estimation from celestial transit solar images using blind deconvolution. J. Space Weather Space Clim., 6, A1, 2016, DOI: 10.1051/swsc/2015040.
All Tables
Reconstruction quality for the different number of simulated transit observations using the synthetic anisotropic Gaussian filter.
All Figures
Fig. 1.
Accumulation of observations of the Venus transit captured by SDO/AIA on June 5th–6th, 2012. 

In the text 
Fig. 2.
Realistic images: (A) first image x_{1}, (B) second image x_{2}, and (C) third image x_{3} from the simulated transit. 

In the text 
Fig. 3.
Realistic filters in a (2b + 1) × (2b + 1) window with b = 16. (A) Anisotropic Gaussian Filter. (B) X Filter. 

In the text 
Fig. 4.
Results for the anisotropic Gaussian filter and P = 3. (A) Noisy observation of image 1 (y_{1}). (B) Filter reconstruction (h*) using blind deconvolution, RSNR = 14.05 dB. (C) Image Reconstruction () using nonblind deconvolution, ISNR = 1.83 dB. 

In the text 
Fig. 5.
Results for the X filter and P = 3. (A) Noisy observation of image 1 (y_{1}). (B) Filter reconstruction (h*) using blind deconvolution, RSNR = 7.52 dB. (C) Image reconstruction () using nonblind deconvolution, ISNR = 2.55 dB. 

In the text 
Fig. 6.
Logarithm of the nonparametric filters estimated for the SDO/AIA telescope taken inside the set Γ. The filters are calculated for several values of σ providing different residual whiteness (measured by ). (A) σ = σ _{RME}, , (B) σ = 2σ_{RME}, , (C) σ = 3σ_{RME}, . 

In the text 
Fig. 7.
Logarithm of filters for the SDO/AIA telescope taken inside the set Γ. (A) Parametric PSF estimated by Poduval et al. (2013), i.e., h_{p1}. (B) Diffraction mesh pattern modeled using the parameters estimated by Poduval et al. (2013), i.e., h_{p2}. (C) Combined parametric/nonparametric filters using the parametric PSF shown on the center, i.e., . 

In the text 
Fig. 8.
Image reconstruction results for the Venus transit observed by SDO/AIA at 00:02 UT on June 6th, 2012. Results are shown for the nonparametric filter, i.e., . (A) Observed image, (B) 2D reconstruction result, and (C, best viewed in color) 1D profile along y = 129. Figures (A) and (B) contain a green line indicating the 1D profiles shown on figure (C). 

In the text 
Fig. 9.
Image reconstruction results for the Moon transit observed by SDO/AIA at 13:00 UT on March 4th, 2011. Results are shown for the nonparametric filter, i.e., . (A) Observed image, (B) 2D reconstruction result, and (C, best viewed in color) 1D profile along y = 186. Figures (A) and (B) contain a green line indicating the 1D profiles shown on figure (C). 

In the text 
Fig. 10.
Reconstruction results for a SDO/AIA image containing an active region. The image has been captured at 10:00 UT on August 8th, 2011. Results are shown for the nonparametric filter, i.e., . (A) Observed image, (B) 2D reconstruction result, and (C, best viewed in color) 1D profile along y = 187. Figures (A) and (B) contain a green line indicating the 1D profiles shown on figure (C). 

In the text 
Fig. 11.
Reconstruction results for the image in Figure 10A. The figure (best viewed in color) shows the ratio (computed using a pixel by pixel division) between the deconvolved and the observed images along y = 187, for the different considered filters: , , and . The horizontal axis is in correspondence with the horizontal axis of Figure 10C. 

In the text 
Fig. 12.
Logarithm of the filters for the SECCHI/EUVI telescope taken inside the set Γ. (A) Parametric PSF given by the euvi_psf.pro procedure of SolarSoft, i.e., . (B) Parametric PSF given by the euvi_deconvolve.pro procedure of SolarSoft, i.e., . (C) Nonparametric PSF, i.e., . 

In the text 
Fig. 13.
Logarithm of the combined parametric/nonparametric filters for the SECCHI/EUVI telescope taken inside the set Γ. (A) Using the parametric PSF given by the euvi_psf.pro procedure of SolarSoft, i.e., . (B) Using the parametric PSF given by the euvi_deconvolve.pro procedure of SolarSoft, i.e., . 

In the text 
Fig. 14.
Image reconstruction results for the Moon transit observed by SECCHI/EUVI at 08:02 UT on February 25th, 2007. Results are shown for the nonparametric filter, i.e., . (A) Observed image, (B) 2D reconstruction result, and (C, best viewed in color) 1D profile along y = 257. Figures (A) and (B) contain a green line indicating the 1D profiles shown on figure (C). 

In the text 
Fig. 15.
Reconstruction results for a SECCHI/EUVI image containing an active region. The image has been captured at 04:02 UT on February 25th, 2007. Results are shown for the nonparametric filter, i.e., . (A) Observed image, (B) 2D reconstruction result, and (C) best viewed in color) 1D profile along y = 141. Figures (A) and (B) contain a green line indicating the 1D profiles shown on the right. 

In the text 
Fig. 16.
Reconstruction results for the image in Figure 15A. The figure shows the ratio (computed using a pixel by pixel division) between the deconvolved and the observed images along y = 141. (A, best viewed in color) Ratio for filters , , and . (B, best viewed in color) Ratio for filters and . The horizontal axis is in correspondence with the horizontal axis of Figure 15C. 

In the text 
Fig. A.1.
Validating the approximation of the longrange PSF impact by a constant inside the disk of Venus. (A) Logarithm of the filter used for the validation: parametric filter estimated by Grigis et al. (2013) with the central pixels set to zero inside a window of size (2b + 1) × (2b + 1), with b = 64. (B) Observed image inside a window containing the disk of Venus. (C) Convolved image inside a window containing the disk of Venus. The intensity within the window is almost constant. 

In the text 
Fig. A.2.
Validating the approximation of the longrange PSF impact by a constant inside the lunar disk. (A) Logarithm of the filter used for the validation: parametric filter given by the euvi_psf.pro procedure of SolarSoft with the central pixels set to zero inside a window of size (2b + 1) × (2b + 1), with b = 64. (B) Observed image inside a window containing the lunar disk. (C) Convolved image inside a window containing the lunar disk. The window intensity is almost constant. 

In the text 