Non-parametric PSF estimation from celestial transit solar images using blind deconvolution

Context: Characterization of instrumental effects in astronomical imaging is important in order to extract accurate physical information from the observations. Optics are never perfect and the non-ideal path through the telescope is usually represented by the convolution of an ideal image with a Point Spread Function (PSF). Other sources of noise (read-out, Photon) also contaminate the image acquisition process. The problem of estimating both the PSF filter and a denoised image is called blind deconvolution and is ill-posed. Aims: We propose a blind deconvolution scheme that relies on image regularization. Contrarily to most methods presented in the literature, it does not assume a parametric model of the PSF and can thus be applied to any telescope. Methods: Our scheme uses a wavelet analysis image prior model and weak assumptions on the PSF filter's response. We use the observations from a celestial body transit where such object can be assumed to be a black disk. Such constraints limits the interchangeability between the filter and the image in the blind deconvolution problem. Results: Our method is applied on synthetic and experimental data. We compute the PSF for SECCHI/EUVI instrument using the 2007 Lunar transit, and for SDO/AIA with the 2012 Venus transit. Results show that the proposed non-parametric blind deconvolution method is able to estimate the core of the PSF with a similar quality than 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 non-parametric methods.


Introduction
Precise characterization of instrumental effects in solar telescope is important in order to retrieve the full dynamics in observations. In normal-incidence EUV telescopes, the incident flux is converted to digital numbers (DN) through a series of steps, each introducing some noise.
The photon emission in the solar atmosphere itself is random and incoherent in nature. It is usually modeled by a Poisson process so that the number of incident photons falling on a position x between time t and t + T obeys a Poisson law Po(L(x, t, T )) of mean L(x, t, T ) where the radiance L is space and time dependent. For a sufficiently large L (i.e. L > 1000), the Poisson distribution may be approximated by a Gaussian function. As the photon beam impinges the imperfect optical system, its dispersion from an ideal optical path leads to a blurred final image. At the same time, a spectral selection is operated on the signal before it reaches the CCD detector. The latter has an heterogeneous response across its surface which may produce an inhomogeneous flat-field and consequently multiplicative noise [DCH08]. Finally, the camera electronics converts photon counts into DN; this conversion adds Gaussian read-out noise.
The non-ideal optical path through the instrument is usually assumed to be well represented by a convolution between an image obtained from an optically perfect instrument and a point spread function (PSF). The PSF describes the response of the entire optics to a Dirac point source. If the true PSF can be determined, then it is possible to recover an original, undistorted image by convolving the acquired image with an estimate of the inverse PSF, this is the so called 'known-PSF deconvolution'. In practice, 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 filter PSF and a denoised version of the image should be recovered altogether. This process is known as 'blind deconvolution'.
In the case of X-rays and EUV solar telescopes, a PSF approximation can be either found using pre-flight instrumental specifications or inferred from science images. Three main components can be distinguished: (1) the core PSF which accounts for effects such as CCD pixelization, charge spreading, and optical scattering, (2) the diffraction from the entrance filter mesh, and possibly from the focal plane filter mesh, and (3) the long-range diffuse scattering due to micro-roughness of the telescope. Since the diffraction pattern is particularly noticeable during bright events, it can be estimated using solar flare images. The long-range scattered component is more elusive. As it is too diffuse to be clearly observable during flares, its estimation typically requires an occulting body such as the Moon or a planet.
Previous studies aimed at providing a dedicated method to correct one given instrument's imperfections. [ME81] described a scheme for estimating the PSF core from Skylab X-ray data using blind iterative deconvolution, while [MAL95] determined a spatially variable PSF for the Yohkoh/soft X-ray telescope using pre-flight calibration. The TRACE EUV PSF was studied by various groups. [LNT01] studied the diffraction pattern thanks to bright flare kernels. These results were further used by [GSM06], who fitted the core PSF to a parametric model, and uses the algorithm of [AD88] for blind deconvolution. [DMW09] aim at estimating the long-range scattering effect of TRACE PSF using the 2004 Venus transit. The scattering portion of the PSF was modeled as the sum of a measured diffraction pattern, similarly to [GSM06], with a circularly symmetric scattering profile. Wiener filtering was used to regularize the deconvolution problem.
The PSFs from the four EUVI instrument channels on STEREO-B were studied by [SFHG12] using a lunar transit. The goal was also to estimate the long-range scattering effect by assuming a parametric piece-wise power-law model. A maximum likelihood approach was chosen to estimate both the filter PSF and the denoised image. A similar method was used to estimate the PSF of the SWAP instrument onboard PROBA2 [SDS + 13]. Finally, [GYW13] used instrumental specifications to obtain the PSFs for all channels of the Atmospheric Imaging Assembly (AIA) onboard SDO. [PDSP13] provided an independent estimation of the same PSF. From a flare image they recovered the diffraction pattern, whereas for the diffuse component a parametric model was fit using an iterative forward model of the lunar interior in the SDO/AIA images of 2011 March 4 lunar transit.
As we have seen, several methods have been used in the literature for the estimation of the filter PSF. Some works have used linear methods such as Wiener filtering [PDSP13] or iterative inverse filtering [AD88] to recover the PSF from the blurred observations. However, compared to non linear algorithms, linear methods have been proven to provide lower quality reconstructions and to present less noise tolerance [FWBP95,DJ98].
Existing non-linear blind deconvolution methods are based on the alternated optimization between the image and the filter. The algorithm is usually driven by the optimization of a data fidelity term, which is based on the acquisition model, and some additional regularization terms. If the observations are corrupted by an additive Gaussian noise, the deconvolution problem is formulated as a regularized least-squares minimization. A proximal alternated minimization method, recently proposed by [ABRS10], allows to handle such a problem for a large amount of regularization functionals. The advantage of this method over usual alternated minimization approaches (e.g., [AA10]) is that it provides theoretical convergence guarantees. To our knowledge, this method have not be used before to deal with the blind deconvolution problem.
On the other hand, in the presence of Poisson noise, the data fidelity term is given by the Kullback-Leibler (KL) divergence, which represents a more complex function to minimize and algorithms such as the one proposed by [ABRS10] cannot be used. To avoid dealing with the KL divergence function, [SFHG12] have approximated the Poisson corrupted data to have a Gaussian noise distribution using a variance stabilization transform (VST). They use a variable projection method to estimate the filter and the image from the blurred observations. [FWBP95] proposed an iterative framework for minimizing the KL divergence using the wellknown method of Richardson-Lucy (RL). The algorithm provides good reconstruction quality with high noise tolerance. However, the method is computationally inefficient and it does not present theoretical guarantees. [PCBB13] have replaced the RL algorithm by the Scale Gradient Projection (SGP) method in order to increase the efficiency and guarantee the convergence of the algorithm to a critical point. The proposed method is highly dependent on a manually tuned stopping criteria and also on a good PSF initialization.

Contribution
In this work, we propose a non-parametric blind deconvolution scheme to estimate the core PSF of an optical telescope that has captured images from a celestial transit. Since the estimation is not performed in the low intensity regions in the Sun, the acquisition is assumed to follow an additive Gaussian model and the problem is solved through the proximal alternated minimization method of [ABRS10]. The considered Gaussian noise model can be easily generalized to Poisson type of noise with the use of a VST.
The proposed method is first tested through simulated realistic scenarios. Results show the ability of estimating 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, respectively, a Venus and a Moon transit. The validation of the recovered filters is based on the reduction of the celestial body apparent emissions. Results show that, when recovering the core PSF, the proposed non-parametric scheme can achieve similar results than parametric methods. Also, if the parametric PSF is incorporated in the imaging model, the quality of the resulting PSF outperforms both the parametric and non-parametric methods considered 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 present unknown or unexpected behavior, an accurate model of the filter cannot be built. In such cases, we require non-parametric models, as the one presented in this paper, which due to their flexibility can handle better some instrument's imperfections.

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 alternated minimization method used for estimating the filter and the image. It briefly presents the method, including the initialization, the parameter estimation, the stopping criteria and the numerical reconstruction. In Section 5 we present a non-blind deconvolution method that is used for the experimental validation. Finally, in Section 6 we present the results on both synthetic and experimental data for the Moon and Venus transits.

Notation
Most of domain dimensions are denoted by capital roman letters, e.g., M, N, . . . Vectors and matrices are associated to bold symbols, e.g., Φ ∈ R M ×N or u ∈ R M , while lowercase light letters are associated to scalar values. The i th component of a vector u reads either u i or (u) i , while the vector u i may refer to the i th element of a vector set. The vector of ones in R D is denoted by 1 D = (1, · · · , 1) T . The set of indices in R D is [D] = {1, · · · , D}. The cardinality of a set C, measuring the number of elements of the set, is denoted by C. The convolution between two vectors u, v ∈ R D for some dimension D ∈ N is denoted equivalently by u ⊗ v = v ⊗ u. The (convex) indicator function ı C (x) of the set C is equal to 0 if x ∈ C and +∞ otherwise. The 2-D discrete delta function, denoted as δ 0 = δ 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 u p The Frobenius norm of Φ is given by Φ 2 F = i j |φ ij | 2 . The support of a vector u ∈ R D is defined as supp u = {i ∈ N D : u i = 0}.

Problem statement
The telescope imaging process can be mathematically modeled as the instrument's PSF h convolving the actual pure scene x, i.e., h⊗x. The convolution operation, represented by ⊗, consists on 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.

Discrete model
Let us consider that, during the transit, the telescope acquires a set of P observations containing a celestial body. In order to limit the computation time, the effective Field-of-View 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. In our work, all n × n images are represented as N = n 2 vectors so that linear transformation of images are identifiable with matrices. Each observed patch is a collection of values gathered in a vector y j ∈ R N , with j = 1, · · · , P . The observed values are modeled as an image x j ∈ R N convolved by the instrument's PSF (h ∈ R N ). The PSF is assumed to be spatially and temporally invariant, thus is the same for each observed patch. The observations are assumed to be affected by an Additive White and Gaussian Noise (AWGN), n j ∈ R N , with (n j ) i ∼ iid N (0, σ 2 obs ). The acquisition process of the telescope is thus modeled as follows: Let Φ(h) ∈ R N ×N be a circulant matrix of kernel h. If we gather the P pure images in a single matrix X = (x 1 , · · · , x P ) ∈ R N ×P , the acquisition model can be transformed into the following matrix form: with Y, N ∈ R N ×P two matrices that similarly gather the observed images and their noises, respectively, of the P images.

Noise estimation
To keep the model (2) mathematically tractable, we consider additive noise. The Gaussian assumption is justified because we aim at recovering the core of the PSF without focusing on the low-intensity parts of the image. A generalization of our method to include Poisson noise is relatively straightforward thanks to the Variance Stabilization Transform (VST) method [Ans48,DFS09,SFHG12]. Under the AGWN assumption, the energy of the residual noise is known to be bounded using the Chernoff-Hoeffding bound [Hoe63]: with high probability for c = O(1). The noise variance (σ 2 ) can be estimated using the Robust Median Estimator [DJ94], which is based on the assumption that the noise is an additional high frequency component in the observed signal.

Blind deconvolution problem
In this section, we aim at reconstructing both X and h from the noisy observations Y. Since the data is corrupted by an AWGN, we can estimate the image and the filter using a least squares minimization, i.e., by finding the values that minimize the energy of the noise: withX ∈ R N ×T andh ∈ R N . The blind deconvolution problem in (4) is ill-conditioned and without regularization it may only provide trivial solutions (X * , h * ) for the image and the filter, e.g., X * = Y and h * = δ 0 . In the coming sections, we first describe the prior information and constraints on the image and the filter and then we use this information to formulate a regularized blind deconvolution problem.

Prior information on the image
In the following, we present in detail the available prior information on the image candidate: 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. 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 j th 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 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., h ⊗ Sparse image prior: as commonly done in image denoising [SM94,Don95], we regularize our method by a wavelet analysis prior model based on the structure of the image. We assume that, if we represent the image candidate on a suitable wavelet basis Ψ ∈ R N ×W made of W vectors in R N , the wavelet coefficients are sparse, i.e., the coefficients 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 [Mal08] from the 1 -norm computation. The set of detail coefficients that are not affected by the black disk is denoted by Θ, with |Θ| = Q. The matrix S Θ ∈ R Q×W is the corresponding selection operator that extracts from any coefficients vector u ∈ R W 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. Note that we use an Undecimated Discrete Wavelet Transform (UDWT), a redundant dictionary. Since the UDWT captures efficiently all image features whatever their location [SED04,EFM10], it allows a better image reconstruction than the traditional DWT.
Image non-negativity: we should note that the image corresponds to a measure of a photon emission process. Therefore, its pixel values must be non-negative everywhere.

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's components but rather at estimating its core, which is the responsible of most of the diffused light.
We can first consider that, since the PSF corresponds to an observation of a point, it is non-negative everywhere. Additionally, by assuming the amount of light entering in the instrument is preserved, the One important assumption in our proposed method is that the filter candidate is centered in the spatial origin of the discrete grid and 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 Field-of-View of the instrument in order to reduce the computation time.
In high dynamics solar images, a given pixel value can be affected by the value of another high intensity pixel located as far as 100 arcsec away [PDSP13]. In such cases, despite a PSF that can rapidly decay far from its center, the high intensity pixel can induce a long-range effect that is not taken into account by a filter with truncated support. Let us define the complementary set of Γ as Γ C = [N ]\Γ, with |Γ C | = N − (2b + 1) 2 . We assume that the actual PSF is composed by two different filters: one denoted by h Γ with support on Γ, which accounts for short-range effects, and another one, denoted h Γ C with support on Γ C , accounting for long-range effects, i.e., h = h Γ + h Γ C . An accurate estimation of the long-range PSF is out of the scope of this work. Its effect inside the disk of the celestial body is modeled as a constant µ, i.e., (h Γ C ⊗ x j ) i ≈ µ; ∀i ∈ Ω j . For a given patch, the acquisition model in (1) can then be approximated as: In this work, we aim at estimating only the short-range PSF, which is called h hereafter. The effect of the long-range PSF, denoted by µ, is computed in a preprocessing stage by the averaged value of the center of several observed transit disks (See Section 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 level of noise (few DNs). In a future work we envisage to replace this rough evaluation by a joint estimation of the value of µ with the filter and the image.
Notice that our method could allow the addition of other convex constraints, e.g., if we know that 0 ≤ h i ≤ g i , for some upper bound g i on h i such as a specific power law decay. However, we will not consider this possibility here as we want to stay as agnostic as possible on the properties of the filter to reconstruct.

Final Formulation
From Sections 3.1 and 3.2, we can formulate the following regularized blind deconvolution problem: withX ∈ R N ×P ,h ∈ R N and µ1 N 1 T P cancels the long-range part of the filter. The regularization parameter, denoted by ρ, controls the trade-off between the sparsity of the image wavelet coefficients and the fidelity to the observations. This important parameter is estimated in Appendix A.1. The problem of estimating the image and the filter in (5) can also be written as where the modified observations. The convex sets P 0 and D are defined as follows:

Proximal Alternated Minimization
In this section, we describe the alternated minimization method used for solving (6) and hence estimating the image and the PSF from the noisy observations. The problem defined in (6) is non convex with respect to both X and h, but it is convex with respect to X (resp. h) if h (resp. X) is known. Such a problem is usually solved iteratively, by estimating the image and the filter alternately [AA10]. The general behavior of this type of algorithm is as follows: It has been shown in [AA10] that this algorithm is able to provide fairly good results. However, there are not theoretical guarantees for its convergence. Recently, [ABRS10] proposed a proximal alternating minimization algorithm where cost-to-move functions are added to the common alternated algorithm. These functions are quadratic costs that penalize too large variations between two consecutive iterations. They guarantee the algorithm convergence provided some mild conditions are met on the regularity of the objective function L (X, h) and the constraints in (6). The algorithm iterates as follows: where λ x and λ h are the cost-to-move penalization parameters, which control the variations between two consecutive iterations. Section 4.1 describes how these parameters are tuned.
The non-convexity of (6) prevents us to reach its global minimizer. However, it can be shown that, since the objective function L (X, h) and the constraints in (6) meet the conditions stated in [ABRS10], the algorithm (7) converges to a critical point of the problem. Later in this section we describe how to stop the iterations such that a critical point of (6) is reached.

Cost-to-move parameters
The presence of cost-to-move parameters (λ x , λ h ) different than zero ensures the convergence of the algorithm (7) to a critical point of (X, h) [ABRS10]. However, their values are not critical in the reconstruction results. In the following, we describe the tuning criteria used in this work which is based on the works of [PV14]. The iterations start with high values of λ x and λ h , which allows us to keep 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.

Stopping criteria
It is very important to find an automatic criteria that allows to stop the iterations in (7) once a critical point of (6) is reached. An usual stopping criteria 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 [AF13b], we use the spectral characteristics of the noise to analyze the quality of the estimation, which allows stopping the minimization when a high quality solution has been obtained. Let us define each residual image at iteration k as the difference between the observed image and the blurred estimate: Since the observation model is degraded by an additive white noise, we know that the residual image is spectrally white if the estimation at iteration k has a good quality, otherwise the residual contains structured artifacts. Therefore, the iterations can be stopped when the residual is spectrally white, i.e., when the 2-D autocorrelation C rj rj (m, n) of each residual image r j needs to be previously transformed to a matrix of n × n pixels.
which is higher for whiter residuals. We then average over the measures obtained for each residual image: The value of M(R (k) ) is small at the beginning of the iterations when the image and the filter are not properly estimated. Then, the value of M(R (k) ) 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 ρ. The iterations are stopped when the maximum of the whiteness measure (M) is reached.

Numerical Reconstruction
The proximal alternated minimization algorithm is summarized in Algorithm 1. The initialization of the image and filter is discussed in Appendix A. The regularization parameter is tuned iteratively. In a first iteration the value of ρ is set based on the ratio between the variance of the observation noise and the 1 -norm of the wavelet coefficients of the observations. Then, the value of ρ is updated using the ratio between the energy of the noise in (3) and the energy of the residual image in (8). For further details see Appendix A.
The first step in Algorithm 1 allows us to estimate the image by minimizing a sum of three convex functions with the image constrained to the set P 0 . This constraint can be handled by adding the convex indicator function on the set, i.e., ı P0 (X). The resulting optimization, containing a sum of two smooth and two non-smooth convex functions, can be solved using the primal-dual algorithm of [CP11] defined in a Product Space [GJVA14].
The second step in Algorithm 1 allows us to estimate the filter by minimizing a sum of two smooth convex functions with the filter constrained to the set D. The optimization problem can be solved using the Accelerated Proximal Gradient Method [PB14]. This algorithm is able to minimize a sum of a non-smooth and a smooth convex function, which allows us to handle the constraint as a convex indicator function on the set, i.e., ı D (h).
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 [AF13a], 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 B.

Filter validation through non-blind deconvolution
This section is dedicated to the validation of the filter estimated using Algorithm 1. Since the true emissions of the Moon and Venus are zero, the deconvolution with a correct filter should remove all the celestial bodies' apparent emissions. Therefore, the validation process consists on using the estimated filter to deconvolve the observations and verifying that, in the reconstructed image, the pixels inside the object's disk are zero. Let us note that this validation is done using a patch that has not been considered previously during the filter estimation, i.e., on y j / ∈ {y 1 , · · · , y P }. To reconstruct the image we formulate a least squares minimization problem regularized using the prior information available on the image. In this process, we can only assume that the image is non-negative and that it has a sparse representation on a suitable wavelet basis Ψ ∈ R W ×N . No knowledge on the celestial body can be added so that this information can be used for the validation. The proposed non-blind deconvolution method reads as follows: where S ∆ ∈ R R×W is the selection operator of the set ∆, with |∆| = R, which contains the detail coefficients. The convex set P is defined as P = u ∈ R N : u i ≥ 0 . The non-blind image estimation is solved using the primal-dual algorithm of [CP11] defined in a Product Space [GJVA14]. Since this algorithm is able to minimize a sum of three non-smooth convex functions, the constraint can be handled as a convex indicator function. Notice that other algorithms could be used such as the Generalized Forward Backward algorithm [RFP13].

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 of the Venus and Moon transits taken, respectively, by the telescopes SDO/AIA and SECCHI/EUVI.

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 and we shift it three times by 128 pixels horizontally and 5 pixels vertically to have in total three images. In each image, we add a black disk of radius R = 48 pixels, centered at the pixel [129 129], which represents the celestial body. Fig. 1 depicts the images from the simulated transit. During the solar transit, the position and size of the celestial body are known at all moments. However, when this event is imaged, the coordinates of the center represented in the discrete image grid may contain an error of maximum one pixel. This uncertainty is simulated in our synthetic experiments by considering the sets Ω j to be defined using a disk of radius r Ω which is one pixel smaller than the actual object's radius, i.e., 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 2 pixels horizontally and 4 pixels vertically, rotated by 45 degrees (see Fig. 2-(left)). The second filter is simulated by an X shape (see Fi. 2-(right)), hereafter called the X filter. Both functions are similar to actual filter responses that we could expect to have in optical telescopes. Let us note that, in these synthetic experiments, there is no need to test the short/long-range assumption on the PSF. The measurements are simulated according to (1), where additive white Gaussian noise N ∼ iid N (0, σ 2 ) is added in order to simulate a realistic scenario. The noise variance (σ 2 ) is set such that the blurred signal to noise ratio BSNR = 10 log 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 with three levels of decomposition [Mal08]. Higher levels can also be considered, however the computational time is notably higher and the reconstruction quality does not significantly improve.
The reconstruction quality of X * with respect to the true image X GT is measured using the Increase in SNR (ISNR) defined as ISNR = 20 log 10 Y − X GT F / X * − X GT F . The reconstruction quality of h * with respect to the true filter h GT is measured using the Reconstruction SNR (RSNR), where RSNR = 20 log 10 h GT 2 / h GT − h * 2 . Let us first analyze the behavior of the algorithm when we increase the amount 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 Fig. 2-(left). The results correspond to an average over four trials and they are presented for P = 1, 2 and 3 observations. For P = 2, 3 we use a warm start by initializing the algorithm with the results obtained for P − 1.
We notice that the reconstruction quality of the filter significantly improves when the amount of observations P increases. The whiteness measure M(R) presents the same behavior, which indicates that the residual image is whiter when more observations are considered. This occurs because when more observations are considered for the reconstruction, the problem of reconstructing the filter becomes better conditioned, i.e., the amount of information available increases without adding more unknowns. However, the number of unknowns increases for the image estimations, thus there is less improvement in the image reconstruction quality. About the numerical complexity, the algorithm requires an average of four iterations on the value of ρ and a total time of 3 hours on our MATLAB implementation running on a 2.6 GHz Intel computer. When the number of iterations on ρ increases, we observed 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 validate the true zero pixel values inside the object. Fig. 3 depicts the first patch (left) and the resulting PSF (center) 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. 3-(center) in the non-blind deconvolution of Section 5. Note that the majority of the pixels inside the estimated disk are zero, except for some numerical errors. To quantify this validation, we use as measurement the disk intensity, i.e., the sum of the pixel values inside the disk. We compare the ratio between the disk intensity for the deconvolved image, denoted by S X , and the disk intensity for the observed patch, denoted by S Y . We obtained a disk intensity ratio of S X /S Y = 8.16 × 10 −2 . This shows the effectiveness of the estimated filter to recover the true zero emissions inside the disk. Fig. 4 depicts the results for P = 3 observations using the X filter of Fig. 2-Right. Let us note the reconstruction quality is lower than in the case of the Gaussian filter because the X filter presents more complexity and hence it is harder to estimate without extra information on its shape. However, since it causes less diffusion on the observation, the image reconstruction quality is better. When comparing the values inside the disk for both the observations and the estimated image we observed that the disk intensity ratio is S X /S Y = 2.03 × 10 −2 , which validates the zero emissions inside the disk Finally, we also considered a case (not displayed here) where the observation is only corrupted by noise. As expected, the filter estimated by the algorithm is a high quality Kronecker delta function with RSNR = 66 dB. This result shows the capability of the algorithm to recover highly localized filters of one pixel radius.

Experimental Data
The proposed 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) [LTA + 12] onboard 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) [HMV + 08], a part of the Sun Earth Connection Coronal and Heliospheric Investigation (SECCHI) instrument suite onboard the STERE0-B spacecraft. In the following, we present for each telescope the results obtained when reconstructing the filter using the blind deconvolution approach and the validation of the obtained filter. Results are compared with previously estimated PSF of both telescopes.

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 to observe the core of the PSF. The radius of the Venus disk is known to be 49 pixels and hence, the disk represents a small area inside the high resolution image. Since the blind deconvolution takes benefit mainly from the black disk in the transit images, we select a patch of size N = 256 2 centered on the Venus disk. As explained in Appendix B, our method considers unknown border pixels based on the filter support. Therefore, cropping the observed image does not have any effects on the reconstruction results. Furthermore, considering smaller observations keeps the optimization problem computationally tractable.
The presence of a long-range PSF was validated in the observations by analyzing the pixels inside the disk of Venus. To estimate this effect, we summed a total of 10 patches and we estimated 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 observed patch to estimate the filter using the blind deconvolution scheme. The noise variance was computed using a Median Robust Estimator, which resulted in a BSNR equal to 29 dB.
[GYW13] estimated a high quality parametric PSF for AIA using the optical characterization of the telescope. This parametric PSF, denoted by h p , is depicted in Fig. 5-(left) and is used for comparison purposes. Fig. 5-(center) presents the core of the non-parametric PSF (h * np ) estimated using the proposed blind deconvolution approach. We observe that the obtained filter is highly localized in the center of the discrete grid and it presents a limited support. Let us note that the noise level prevents the estimation of the diffraction peaks present in h p . These can only be obtained by constraining the shape of the filter in the reconstruction process, as implicitly done by parametric deconvolution methods.
In order to obtain a more precise PSF estimation, let us incorporate the parametric PSF inside the acquisition model in (1) by considering a combined parametric/non-parametric PSF defined as h p−np = h np ⊗ h p . The non-parametric PSF (h np ) was estimated using the proposed blind deconvolution approach and the resulting PSF is presented in Fig. 5-(right).
To solve the blind deconvolution problem, the algorithm requires an average of 2 iterations on the value of ρ and a total time of 6 hours on our computer. The whiteness measure M(R) was observed to decrease to a value of −0.07 at the final iteration. Moreover, we observed that when the number of iterations on ρ increases, the estimated residual energy is closer to the actual noise energy.
The filters from Fig. 5 were validated using the non-blind deconvolution formulation from Section 5 to estimate the image. Let us note that the patch used for the filter validation has not been used for estimating h * . Fig. 6 depicts the 2-D estimated images and a 1-D profile of those images along y = 129 (to allow the observation of the Venus disk) using: (1) the parametric PSF, (2) the non-parametric PSF and (3) the resulting parametric/non-parametric PSF. The apparent Venus emissions can be quantified by integrating the pixel values inside the disk. For this, we computed the disk intensity ratio for the three deconvolved images. Results are presented in Table 2.
0.61 × 10 −3 Non-parametric filter (h * np ) 3.1 × 10 −3 Parametric/non-parametric filter (h * p−np ) 0 We observe that the parametric/non-parametric model reduces completely the disk apparent emissions and it provides an enhanced image reconstruction. We also note that, compared to the non-parametric filter estimation, the parametric model provides slightly better results in terms of reducing the disk apparent emissions. However, the non-parametric model presents the advantage of being generally applicable for any optical instrument without the need of an exact model of the imaging process.
Let us note that, if the long-range effect (µ) is not removed from the observations, the parametric PSF taken with a larger support of 1024 × 1024 pixels, is not able to eliminate the offset and hence, is not able to remove the apparent emissions inside the Venus disk.

SECCHI/EUVI -Moon transit
We consider images calibrated with the secchi-prep procedure available within the IDL SolarSoft library. Three 2048 × 2048 images from the transit (P = 3) recorded by the 17.1 nm channel of EUVI were used. Following the practice described for the SDO/AIA 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).
The long-range PSF is modeled by a constant µ = 12.5, computed as the mean value inside the Moon disk using the sum of 5 patches. The estimated noise variance provided a BSNR equal to 35 dB. Similarly to the SDO/AIA transit, we provide a comparison between three different PSFs: (1) a parametric PSF given by the euvi-psf procedure of SolarSoft; (2) a non-parametric PSF estimated using the proposed blind deconvolution approach; and (3) a resulting parametric/non-parametric PSF. The different filters are presented in Fig. 7. Due to lack of space, the estimated images using the non-blind deconvolution are not presented. However, we provide in Table 3 the quantitative validation results using the disk intensity ratio estimated on the deconvolution of an observation that was not considered for the filter estimation. We noticed that, as for the Venus transit images, when the long-range effect is not removed from the observations, the parametric PSF (considered in a larger support of 1024 × 1024 pixels) is not able to remove the offset and to recover zero emissions inside the Venus disk. If the long-range effect is considered through the parameter µ, the parametric PSF presents a better behavior and achieves a slightly higher Filter S X /S Y Parametric PSF (h p ) 0.91 × 10 −2 Non-parametric filter (h * np ) 1.57 × 10 −2 Parametric/non-parametric filter (h * p−np ) 0 reduction of the disk emissions compared to the non-parametric filter (See Table 3). We also observed that the parametric/non-parametric model reduces completely the disk apparent emissions. When comparing the results for the Moon transit with the ones obtained for the Venus transit, we can validate the generality of the proposed blind deconvolution approach to deal with different type of imaging instruments.

Conclusion
We have demonstrated how a non-parametric blind deconvolution technique is able to estimate the core of the PSF of an optical instrument with high quality. The quality of the estimated PSF core is comparable with the one provided by parametric models based on the optical characterization of the imaging instrument. We also demonstrate that, if the parametric PSF is incorporated in acquisition model, the blind deconvolution approach is able to provide a "corrected" PSF such that all the apparent emissions inside the disk of a celestial body during a solar transit can be removed. Let us note that non-parametric techniques cannot outperform the accuracy of parametric methods, however in situations where the telescope imaging model cannot be obtained due to some instrument's imperfections, the non-parametric method presents a great advantage. Moreover, the proposed method is not specific to a given instrument but can be applied to any optical instrument provided that we have strong knowledge on some image pixels values, such as the information available during the transit of the Moon or a planet. We also have shown that the use of the Proximal Alternated Minimization technique proposed by [ABRS10], 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. In this work, only the short-range PSF was estimated optimally. In the future, we would like to address the problem of estimating the long-range PSF jointly with the short-range PSF and the image.

5:
ε (l) = R (l) if M (l+1) < M (l) then break. 8: Return X * = X (l+1) and h * = h (l+1) the choice of this parameter when solving general ill-posed inverse problems. In basis pursuit denoising, [DJ94] proposed a universal value for ρ given by σ 2log(N ). Since this value increases with the amount of samples (N ), it tends to provide overly smoothed images. Another state-of-the-art choice, proposed by [DJ95], is based on the minimization of Stein's Unbiased Risk Estimate (SURE) [Ste81]. 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 [CYV00], to estimate the parameter ρ using the noise characteristics and the probability distribution function of the wavelet coefficients: 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: √ 2 P Q P j=1 (A) j 1 , where A = S Θ Ψ T X 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 (12) is determined using the observations Y, i.e., τ = τ Y = √ 2 P Q P j=1 (A ) j 1 , with A = S Θ Ψ T Y. Since this value is higher than the actual ρ, it will lead to over-regularized images. Therefore, we propose to refine the value of ρ by an iterative update based on the discrepancy principle between the actual noise energy in (3) and the energy of the residual image (8). 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.

A.2 Image and filter
For the first value of ρ, the alternated algorithm (Algorithm 1) is initialized with the trivial solution, where the image is given by the observations Y 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 Convolutions with unknown boundary conditions
As explained in Section 2, our blind deconvolution approach is realized in a restricted Field-of-View of size n × n. For all our computations (as those realized in Algorithm 1) fast convolution methods exploiting the Fast Fourier Transform must be performed. Therefore, it is really 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 [AF13a], we implement their unknown boundary conditions. In a nutshell, this 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 imagex j ∈ R M and a larger filterh ∈ R M , with M = N + 2b. The observation is obtained by selecting only the pixels inside the boundaries of the convolved image using a selection operator S N ∈ R N ×M . The problem in (1) can be rewritten as: The regularized blind deconvolution problem in (5) transforms as follows: s.t. (x j ) i = 0 if i ∈ Ω j ; (x j ) i ≥ 0 otherwisẽ h ∈ PS; supph ∈ Γ withX ∈ R M ×P andh ∈ R M . The actual image X * and filter h * are obtained by applying the operator S N to both the estimated imageX * and filterh * , respectively, i.e., X * = S NX * , h * = S Nh * . 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.