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

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

A. González111Institute of Information and Communication Technologies, Electronics and Applied Mathematics (ICTEAM), Université catholique de Louvain (UCL), B-1348 Louvain-la-Neuve, Belgium. adriana.gonzalez@uclouvain.be, laurent.jacques@uclouvain.be,  V. Delouille222Royal Observatory of Belgium, Avenue Circulaire 3, 1180 Bruxelles, Belgium. v.delouille@oma.be  and L. Jacques
July 4, 2019
Abstract

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 (read-out, photon-counting). The problem of estimating both the PSF 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, 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 non-parametric 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 non-parametric methods.

1 Introduction

Deconvolution is an ubiquitous data processing method that arises in a variety of applications, e.g., biomedical imaging [JSM02], astronomy [PCZBB13, She13], remote sensing [DFXL11] and video and photo enhancement [FSH06]. Formally, this problem amounts to recovering a signal from blurred and corrupted observations :

(1)

where stands for the convolution operation between and some other signal , while is a general function accounting for some corrupting noise, e.g., under an additive corruption model with some (Gaussian) noise .

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

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

Since there are infinite combinations of image and filter that are compatible with the distorted observations, the blind deconvolution problem is severely ill-posed. One way to reduce the number of unknowns is to introduce a forward (parametric) model of the PSF. For instance, Oliveira et al. [OFB07] 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. [BMK09] have modeled smoothly varying PSFs by using a simultaneous autoregressive prior, where the only parameter to estimate is the variance of a Gaussian function. However, such methods are limited to specific applications as they can only recover the expected model of the PSF, excluding the estimation of a generic non-parametric PSF. Furthermore, due to the ill-posedness of the blind deconvolution problem, a slight mismatch between the specified model and the true PSF can lead to poorly deconvolved images.

Blind-deconvolution is also very sensitive to the noise present in the observations. While some works [AD88, GSM06] 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 [FWBP95].

Robust blind deconvolution methods, on the other hand, optimize a data fidelity term, which is based on the acquisition model, stabilized by some additional regularization terms. If the observations are corrupted by an additive Gaussian noise, the deconvolution problem can be solved by a regularized least-squares (LS) minimization [AA10]. In the presence of Poisson noise, the problem can be formulated using the Kullback-Leibler (KL) divergence as the data fidelity term [FWBP95, PCZBB13, PCBB13], 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) [DFS09, SFHG12, She13]. 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. [ABRS10] in the context of additive Gaussian noise. This algorithm allows handling LS problems for a large amount of regularization functionals. The advantage of this method over usual alternating minimization approaches, e.g., Fish et al. [FWBP95, AA10], 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. [ABRS10].

In optical telescopes, the PSF originates from various instrumental effects, such as: optical aberrations (spherical, astigmatism, coma), diffraction (produced by, e.g., an entrance filter mesh), scattering (from, e.g., micro-roughness of the mirrors resulting in long-range diffuse illumination), charge spreading, etc. All these effects ‘spread’ light, i.e., incident photons which would otherwise focus to a single point of the focal place, 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.,  [DMW09].

The PSF of solar telescopes is usually modeled using pre-flight instrument specifications [MAL95, GYW13]. 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. [GSM06] fitted the central pixels of the TRACE PSF to a Moffat function, while DeForest et al. [DMW09] modeled the long-range effect of the TRACE PSF as the sum of a measured diffraction pattern with a circularly symmetric scattering profile. The PSFs from the four EUVI instrument channels on STEREO-B were studied by Shearer et al. [SFHG12, She13], where the long-range scattering effect was assumed to follow a parametric piece-wise power-law model. A similar method was used to estimate the PSF of the SWAP instrument on board the PROBA2 satellite [SDS13]. Finally, Poduval et al. [PDSP13] provided a semi-empiric model for the Atmospheric Imaging Assembly (AIA) onboard SDO, similar to the one used by DeForest et al. [DMW09] for TRACE. In all these works, the complete set of parameters that were fitted are in the range of five to eleven 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 [DMW09, SFHG12, PDSP13] have used the information from a celestial transit to help inferring the PSF of a given instrument. A celestial transit, shown in Figure 1, is an astronomical event where the Moon or a planet move 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.

Figure 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 [MHI14].

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 (sub-pixel) 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 transitting 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 also 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 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. [PDSP13] and a similar study of the available PSFs of SECCHI/EUVI (e.g.,  [SFHG12]). The PSF core accounts for charge spreading, optical aberrations, and diffraction effects. Note that the same algorithm would be able to handle a parametric model of the PSF provided a precise instrument characterization is available. For simplicity, we consider an average additive Gaussian noise model that gathers all sources of noise present in the observation. A detailed discussion on how to handle more general noise models is provided in Section 7.

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

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

1.2 Outline

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

1.3 Notations

We denote by , and the sets of natural, integer and real numbers, respectively. The set of the non-negative real numbers is denoted by . Most of domain dimensions are denoted by capital roman letters, e.g., Vectors and matrices are associated with bold symbols, e.g., or , while lowercase light letters are associated with scalar values. The component of a vector reads either or , while the vector may refer to the element of a vector set. The vector of ones in is denoted by . The set of indices in is . The support of a vector is defined as . The cardinality of a set , measuring the number of elements of the set, is denoted by . The convolution between two vectors for some dimension is denoted equivalently by . We denote the class of proper, convex and lower-semicontinuous functions from a finite dimensional vector space (e.g., ) to  [CoP11]. The (convex) indicator function of the set is equal to if and otherwise. The 2-D discrete delta function, denoted as , is equal to 1 if , and otherwise. The transposition of a matrix reads . For any , the -norm of is . The Frobenius norm of is given by . The Gaussian distribution of mean and variance is denoted by .

2 Problem statement

The telescope imaging process can be mathematically modeled as the instrument’s PSF convolving the true image , i.e., . 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 images containing a celestial body. In order to limit the computation time, the effective Field-of-View (FoV) of each observation is restricted to an image patch centered on the transiting celestial body with a size several times bigger than the body apparent diameter. The effects of this FoV truncation will be discussed in detail in Section 3 and Section 6. In our work, all images are represented as vectors so that linear transformations of images are identifiable with matrices. Each observed patch is a collection of values gathered in a vector , with . 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 are uncompressed and have been corrected for CCD effects (dark current removal, flat fielding, despiking). In the case of SDO/AIA, this corresponds to level 1 data processing.

In 1-D, the discrete circular convolution of a vector with a filter can always be described by the product of with a circulant matrix

(3)

i.e., a matrix where every row is a right cyclic shift of the kernel  [Gra06]. In 2-D and in particular for the model (2), the convolution of by can still be represented by the multiplication of by matrix that is now block-circulant with circulant blocks, i.e., a matrix made of blocks, each block being a circulant matrix representing the action of one row of the 2-D filter in (2).

If we gather the ‘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 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 [Ans48, DFS09, SFHG12] in the acquisition model in (4), both on the observations and on the convolution result (see Section 7 for more details).

Under the AWGN assumption, the energy of the residual noise is known to be bounded using the Chernoff-Hoeffding bound [Hoe63]:

(5)

with high probability for . A first guess of the noise variance 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. 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 (5) considers an average noise model and since (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 .

3 Blind deconvolution problem

In this section, we aim at reconstructing both and from the noisy observations . 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:

(6)

with and .

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

3.1 Prior information on the image

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

(a) Zero disk intensity:

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

(b) Analysis-based sparsity model:

As commonly done with ill-posed inverse problems associated to image restoration tasks [SM94, Don95, SF09], we regularize our method with a 2-D wavelet prior. Compared to a frequency representation, as obtained by the Fourier transform, a 2-D wavelet transform allows to represent images with a multi-resolution scheme [Mal08], 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) [SM94, Don95]. More recently, better inversion results have been obtained using analysis-based sparsity models [CMW12, SJDAJ14], which rather promote sparse projections of the estimated image over a redundant system of wave-functions, i.e., a dictionary. This set can be much larger than an orthonormal basis and hence it can efficiently capture much more different image features.

Following this recent trend, we decide to adopt such an analysis prior that, in opposition to the synthesis framework, also allows us to work directly on the image domain without increasing the dimensionality of the optimization problem. As a dictionary, we use the system associated to the Undecimated Discrete Wavelet Transform (UDWT) [SMF10]. 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 [SED04, EFM10], 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 vectors in , the wavelet coefficients are sparse, i.e., the coefficients matrix has few important values and its -norm (computed over all its entries) is expected to be small. However, the wavelet coefficients whose support touches the boundary of the Moon or Venus disk are not sparse. Hence, we do not consider in this prior model those coefficients that are affected by the occulting body. Also, we follow a common practice in the field which removes the (unsparse) scaling coefficients [Mal08] from the -norm computation. The set of detail coefficients that are not affected by the black disk is denoted by , with . See Section 6.1 for further details on how to compute the set . The matrix is the corresponding selection operator that extracts from any coefficients vector only its components indexed by . The global prior information can thus be exploited by promoting a small -norm on the wavelet coefficients belonging to . The rationale of this prior is to enforce noise canceling and avoid artifacts in the reconstructed image.

(c) Image non-negativity:

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

3.2 Constraints on the PSF

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

We first consider that, since the PSF corresponds to an observation of a point, it is non-negative everywhere. Additionally, by assuming that the amount of light entering in the instrument is preserved, the -norm of the filter candidate must be equal to one. This is explained by taking the sum over the pixels of one observed patch in (4), which is equal to in a noiseless process. Since the filter is non-negative, taking its -norm equal to one is equivalent to , which results in , i.e., the light is preserved. Therefore, the filter candidate must belong to the Probability Simplex [PB14], 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 , i.e., the filter has only important values in the central pixels and is negligible beyond the considered support. This allows us to work with a patch (only a part of the observation) and not with the complete FoV of the instrument in order to reduce the computation time.

EUV images have a high dynamic range. Hence, due to long-range effects, a given pixel value can be affected by the value of another high intensity pixel located as far as 100 arcsec away for SDO/AIA [PDSP13] and 1500 arcsec away for SECCHI/EUVI [SFHG12]. 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 , with . We assume that the actual PSF is composed by two different filters: the PSF core with support on , denoted by , which accounts for short-range effects, and another one, denoted with support on , accounting for long-range effects, i.e., . An accurate estimation of the long-range PSF is out of the scope of this work. Its effect inside the disk of the celestial body is modeled as a constant , i.e., (see Appendix A for more details). For a given patch, the acquisition model in (2) can then be approximated as: In this work, we aim at estimating only the PSF core, which is called hereafter. The effect of the long-range PSF, denoted by , is computed in a preprocessing stage by the average value inside the center of several observed transit disks (see 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 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 Section 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 , for some upper bound on such as a specific power law decay. However, we will not consider this possibility here as we want to stay as agnostic as possible on the properties of the filter to reconstruct.

3.3 Final Formulation

From Sections 3.1 and 3.2, we can formulate the following regularized blind deconvolution problem:

(7)
s.t.

with , and cancels the long-range part of the filter. The regularization parameter, denoted by , controls the trade-off between the sparsity of the image projection in a wavelet dictionary and the fidelity to the observations. This essential parameter is estimated in Appendix B.2.

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

The problem of estimating the image and the filter in (7) can also be written as

(8)

where is the objective function, with the modified observations. The convex sets and are defined as follows: ; .

4 Proximal Alternating Minimization

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

The problem defined in (8) is non-convex with respect to both and , but it is convex with respect to (resp. ) if (resp. ) 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 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 and are the cost-to-move penalization parameters that control the variations between two consecutive iterations. Section 4.1 describes how these parameters are tuned.

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

4.1 Cost-to-move parameters

The presence of cost-to-move parameters (, ) different than zero ensures the convergence of the algorithm (9) to a critical point of  [ABRS10]. 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 and Vandergheynst [PV14]. The iterations start with high values of and , keeping the image and the filter estimations closer to the initial values. When the number of iterations increases, the filter and the image estimates become more and more accurate and the parameters and 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 and Figueiredo [AF13b], we use the spectral characteristics of the noise to analyze the quality of the estimation. This allows stopping the minimization when a high quality solution has been obtained. Let us define each residual image at iteration 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 has a good quality, otherwise the residual contains structured artifacts. Therefore, the iterations can be stopped when the residual is spectrally white, i.e., when the 2-D autocorrelation of each residual image is approximately the function . Let us note that, for the computation of the 2-D autocorrelation function, each residual vector needs to be previously transformed into a matrix of pixels and normalized to zero mean and unit variance.

Considering a window, we compute the distance between the autocorrelation and the Kronecker function as follows:

(11)

which is higher for whiter residuals. As suggested by Almeida and Figueiredo [AF13b], in our experiments we have used . 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 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 in [AF13b].

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.

1: ; ; ; ; MaxIter =20
2:for  to MaxIter do
3:  1st step, image estimation:
4:     
5:  2nd step, filter estimation:
6:     
7:  Parameters update:
8:     
9:     
10:  Compute residuals and whiteness measure:
11:     
12:     Compute using (12)
13:  Stop when residual is spectrally whiter:
14:     if   then break.      
15:Return and .
Algorithm 1 Proximal Alternating Minimization Algorithm

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

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 [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., . The APG algorithm is briefly described in Appendix C.2.

The algorithms used to solve the minimization problems described above were chosen because they are well suited to solve each problem optimally. Since the cost functions are different in each step, the same algorithm could not be used to optimally find a solution of both problems. The presence of two non-smooth functions in the first step, prevents the use of the APG algorithm. For the second step, the algorithm of Chambolle and Pock [ChP11] 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 and Figueiredo [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 D.

5 Filter validation through non-blind deconvolution

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

To obtain this reconstructed image, we formulate a least squares minimization problem regularized using the prior information available on the image, that is, that the image is non-negative and has a sparse representation on a suitable wavelet basis . The proposed non-blind deconvolution method reads as follows:

(13)
s.t.

where is the selection operator of the set , with , which contains the detail coefficients. The convex set is defined as . The non-blind image estimation is solved using an extension of the primal-dual algorithm of Chambolle and Pock [ChP11] (see Section C.1 for more details). Since this algorithm is able to minimize a sum of three non-smooth convex functions, the constraint can be handled as a convex indicator function. Notice that other algorithms could be used such as the Generalized Forward Backward algorithm [RFP13].

6 Experiments

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

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

6.1 Synthetic Data

Three synthetic images are selected to test the reconstruction (). They are defined on a pixel grid (). 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 pixels each. The center of the cutouts are 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 pixels, centered at the pixel , which represents the celestial body. Figure 2 depicts the images from the simulated transit.

Figure 2: Realistic images: (left) first image , (center) second image and (right) third image 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 and to be defined using disks of radius and , respectively, with one pixel difference with respect to the actual object’s radius, i.e., and .

Two kinds of discrete filters are selected and they have a limited support on a pixel grid (). 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 Figure 3-(left)). The second filter is simulated by a X shape (see Figure 3-(right)), hereafter called the X filter. These functions help to demonstrate the capacity of the proposed method to recover not only diffusion filters such as the anisotropic Gaussian, but also diffraction filters such as the X example. Let us note that, in these synthetic experiments, there is no need to test the long-range assumption on the PSF.

Figure 3: Realistic filters in a window with . (left) Anisotropic Gaussian Filter. (right) 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 () is set such that the blurred signal to noise ratio is equal to 30 dB, which corresponds to a realistic BSNR in the actual observations.

During the experiments, the wavelet transformation used for the sparse image representation is the redundant Daubechies wavelet basis with two vanishing moments and three levels of decomposition [Mal08]. 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 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 indexes of the detail coefficients that are equal to zero.

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

Let us first analyze the behavior of the algorithm when we increase the number of observations 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 3-(left). The results correspond to an average over four trials and they are presented for = 1, 2 and 3 observations. For we use a warm start by initializing the algorithm with the results obtained for .

ISNR [dB] RSNR [dB]
1 1.27 9.63 -0.25
2 1.33 13.07 -0.16
3 1.36 14.05 -0.03
Table 1: Reconstruction quality for 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 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 hours. 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 4-(left) depicts the first noisy observed patch and Figure 4-(center) presents the resulting PSF when reconstructing the anisotropic Gaussian filter for . We can observe how the algorithm is able to estimate the filter with a high quality, even when no hard constraints are introduced on the filter shape. The reconstruction quality can be further improved if stronger conditions are assumed on the filter.

Figure 4: Results for the anisotropic Gaussian filter and . (left) Noisy observation of image 1 (). (center) Filter Reconstruction () using blind deconvolution, RSNR = 14.05 dB. (right) Image Reconstruction () using non-blind deconvolution, ISNR = 1.83 dB.

Figure 4-(right) presents the reconstructed image using the estimated filter from Figure 4-(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 , 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 observations using the X filter of Figure 3-(right). 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.

Figure 5: Results for the X filter and . (left) Noisy observation of image 1 (). (center) Filter Reconstruction () using blind deconvolution, RSNR = 7.52 dB. (right) Image Reconstruction () using non-blind 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 high quality Kronecker delta function with RSNR = 66 dB. This result shows the capability of the algorithm to recover highly localized filters of one pixel radius.

6.2 Experimental Data

The proposed blind deconvolution method was tested on two experimental sets of data. A first data set corresponds to the observations of the Venus transit on June 5th - 6th 2012 by the Atmospheric Imaging Assembly (AIA) [LTA12] 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) [HMV08], a part of the Sun Earth Connection Coronal and Heliospheric Investigation (SECCHI) instrument suite on board the STEREO-B spacecraft. These images are compressed using the RICE algorithm [Nig11], which is lossless [PB14] and hence does not introduce additional errors in the PSF estimation.

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

6.2.1 SDO/AIA - Venus transit

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

The radius of the Venus disk is known to be 49 pixels and hence, the disk represents a small area inside the high resolution image. Since the blind deconvolution takes benefit mainly from the black disk in the transit images, we select a patch of size centered on the Venus disk. As explained in Appendix D, our method considers unknown border pixels based on the filter support (completely defined by 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 , the value estimated by the Robust Median Estimator (RME). Then, the blind deconvolution method is performed for different values of that are multiples of . 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 of the signal features. The variance is computed as follows [DJ94]:

(14)

where is a vector containing the finest scale wavelet coefficients of the observed vector , i.e., , with the selection operator of the set that contains the finest scale wavelet coefficients. This estimation resulted in 2.81 DN.

Figure 6 shows the resulting PSFs for different values of . We can observe that the whitest residual is obtained for 5.62 DN and we thus select this value for our non-parametric PSF estimate depicted in Figure 6-(center). 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 Figure 6-(left)). 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 Figure 6-(right)).

, , ,
Figure 6: Logarithm of the non-parametric 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 ).

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

To solve the blind deconvolution problem, the algorithm requires an average of four iterations on the value of and a total time of 8 hours 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 [PB14]. The computation time can be further reduced if, for instance, the algorithms are implemented in C instead of MATLAB.

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

Figure 7: Logarithm of filters for the SDO/AIA telescope taken inside the set . (left) Parametric PSF estimated by Poduval et al. [PDSP13], i.e., . (center) Diffraction mesh pattern modeled using the parameters estimated by Poduval et al. [PDSP13], i.e., . (right) Combined parametric/non-parametric filters using the parametric PSF shown on the center, i.e., .

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

The estimated filters from Figure 6-(center), Figure 7-(left) and Figure 7-(right) were validated using the non-blind deconvolution formulation from Section 5. Notice that, similarly to what has been done for the blind deconvolution, the value of has been selected adaptively by optimizing the residual whiteness defined in (11). We also observed empirically that the value of must be kept smaller than the one obtained in Algorithm 2. The validation was done on transit images using the modified observations with DN. A first validation was performed on the Venus transit image taken by SDO/AIA at 00:02 UT on June 6th 2012 (see Figure 8-(left)), a patch that has not been used before for estimating . 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 Figure 9-(left)).

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. [PDSP13], i.e., ; (2) the non-parametric PSF of Figure 6-(center), i.e., ; and (3) the parametric/non-parametric PSF of Figure 7-(right), i.e., .

Filter Venus Moon
Table 2: Disk intensity ratio for the different filters.

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

Figures 8 and 9 illustrate the image reconstruction results when using the non-parametric PSF, i.e., . The figure on the left presents the observed patch, the one on the center depicts the 2-D estimated image and the one on the right presents a 1-D profile that allows the observation of the disk of Venus and the Moon, respectively. Due to lack of space, only the non-parametric PSF validation is illustrated.

Figure 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 non-parametric filter, i.e., . (left) Observed image, (center) 2-D reconstruction result, and (right, best viewed in color) 1-D profile along . The figures on the left and on the center contain a green line indicating the 1-D profiles shown on the right.
Figure 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 non-parametric filter, i.e., . (left) Observed image, (center) 2-D reconstruction result, and (right, best viewed in color) 1-D profile along . The figures on the left and on the center contain a green line indicating the 1-D profiles shown on the right.

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

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

Figure 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 non-parametric filter, i.e., . (left) Observed image, (center) 2-D reconstruction result, and (right, best viewed in color) 1-D profile along . The figures on the left and on the center contain a green line indicating the 1-D profiles shown on the right.

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

Figure 11: Reconstruction results for the image in Figure 10-(left). The figure (best viewed in color) shows the ratio (computed using a pixel by pixel division) between the deconvolved and the observed images along , for the different considered filters: , and . The horizontal axis is in correspondence with the horizontal axis of Figure 10-(right).

6.2.2 SECCHI/EUVI - Moon transit

We consider three images from the transit () 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 window () centered around the Moon disk. The filter is assumed to have a limited support inside a pixel grid (). This allows one to obtain the core of the PSF of around pixels, as observed by Shearer et al. [SFHG12], and encompass 99% of the energy of previously estimated PSFs.

The long-range PSF is modeled by a constant 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 = 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. [SFHG12], i.e., , and given by the euvi_deconvolve.pro procedure of SolarSoft; (3) the non-parametric PSF, i.e., ; (4) the parametric/non-parametric PSF obtained by incorporating in the acquisition model the parametric PSF given by euvi_psf.pro, i.e., ; and (5) the parametric/non-parametric PSF obtained by incorporating in the acquisition model the parametric PSF given by euvi_deconvolve.pro, i.e., . The core of the parametric and non-parametric filters are presented in Figure 12 and the core of the resulting parametric/non-parametric filters are depicted in Figure 13.

Figure 12: Logarithm of the filters for the SECCHI/EUVI telescope taken inside the set . (left) Parametric PSF given by the euvi_psf.pro procedure of SolarSoft, i.e., . (center) Parametric PSF given by the euvi_deconvolve.pro procedure of SolarSoft, i.e., . (right) Non-parametric PSF, i.e., .
Figure 13: Logarithm of the combined parametric/non-parametric filters for the SECCHI/EUVI telescope taken inside the set . (left) Using the parametric PSF given by the euvi_psf.pro procedure of SolarSoft, i.e., . (right) 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 degrees. The estimated non-parametric PSF clearly favors one of the two diagonals. Nevertheless, it presents an orientation at  degrees. The source of this difference in orientation is unknown and it should be further investigated. This, however, is beyond the scope of the present paper. Both parametric/non-parametric PSFs present a similar behavior, favoring both diagonals but with a slight dominance of that at  degrees.

The filters from Figure 12 and Figure 13 were validated using the non-blind deconvolution described in Section 5. Similarly to what has been done for SDO/AIA, the value of has been selected adaptively by optimizing the value of the residual whiteness. Let us note that, as observed empirically, the value of must be kept smaller than the one obtained in Algorithm 2. The filter validation was performed on the Moon transit image taken by SECCHI/EUVI at 08:02 UT on February 25th 2007 (see Figure 14-(left)). Let us note that this patch has not been used before for estimating and that the non-blind deconvolution is applied to the modified observation using 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.

Filter
5.72
5.80
2.02
1.81
1.29
Table 3: Moon disk intensity ratio for the different filters.

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

In Figure 14 we illustrate the reconstruction results when using the non-parametric PSF, i.e., . The figure on the left presents the observed patch, the one on the center depicts the 2-D estimated image and the one on the right depicts a 1-D profile of along (to allow the observation of the lunar disk). Due to lack of space, only the non-parametric PSF validation is illustrated.

Figure 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 non-parametric filter, i.e., . (left) Observed image, (center) 2-D reconstruction result, and (right, best viewed in color) 1-D profile along . The figures on the left and on the center contain a green line indicating the 1-D profiles shown on the right.

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

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

Figure 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 non-parametric filter, i.e., . (left) Observed image, (center) 2-D reconstruction result, and (right, best viewed in color) 1-D profile along . The figures on the left and on the center contain a green line indicating the 1-D profiles shown on the right.

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

Figure 16: Reconstruction results for the image in Figure 15-(left). The figure shows the ratio (computed using a pixel by pixel division) between the deconvolved and the observed images along . (top, best viewed in color) Ratio for filters , , and . (bottom, best viewed in color) Ratio for filters and . The horizontal axis is in correspondence with the horizontal axis of Figure 15-(right).

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 Figure 1). The proposed non-parametric approach is comparable to previous parametric ones with respect to the quality of the PSF core. It presents, however, some limitations in terms of the considered noise model and the estimation of the PSF long-range effects. Fortunately, the optimization techniques we use in this work are flexible enough to open perspectives of improvement in future works.

7.1 Noise model

The incidence of photon flux on a EUV telescope is converted to digital numbers (DN) through a series of steps, each potentially introducing some noise. The beam of photons impinges the optical system where the PSF acts as a blurring operator. Simultaneously, a spectral selection is performed on the signal before it reaches the CCD detector. The latter has an heterogeneous response across its surface. Finally, the camera electronics convert electrons into DN, adding the read-out noise. The pixels in the resulting image can be modeled as the realization from a random variable , whose noise part can be decomposed into additive, Poissonian, and multiplicative degradations. The expectation () and the variance () of then verify:

(15)

where is the standard deviation of the additive component and are parameters.

Shearer et al. [SFHG12] handled Poisson corrupted data through a Variance Stabilization Transform (VST) which provides an approximated Gaussian noise distributed data. In order to generalize our AWGN model to models such as (15), a VST could be included in the -fidelity term in (7) and (13), on both the observations and the convolution result [SFHG12]. 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 [DFS09, AABM12]. Notice, however, that such a stabilization cannot be applied only on the observations as a mere preprocessing of the data, while using afterwards other deconvolution methods assuming additive Gaussian noise corruption. Indeed, the VST being a non-linear process of the form for some  [Ans48], it breaks the convolutive model (2) on which those deconvolution methods rely. In other words, the VST of the convolution of by in (2) is not equivalent to (and hardly approximable by) the convolution of the variance stabilized vectors, which breaks the applicability of those techniques. As an alternative to using a VST, the true Poisson distribution could also be used to design a specific fidelity term, which is based on the negative log-likelihood of the posterior distribution induced by the observation model, i.e., the KL divergence of this distribution [FWBP95, AABM12, PCBB13]. However, it is unclear how to adapt the method proposed by Attouch et al. [ABRS10] to the resulting framework.

7.2 Estimation of long-range effects

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

Regarding the approximation of the long-range PSF impact by a constant, we notice that, instead of estimating the constant in a pre-processing step (see Section 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 sub-problems detailed in Section 4.

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

8 Conclusion

We have demonstrated how a non-parametric blind deconvolution technique is able to estimate the core of the PSF of an optical instrument with high quality. The quality of the estimated PSF core is comparable with the one provided by parametric models based on the optical characterization of the imaging instrument. We also demonstrate that, if the parametric PSF is incorporated in the acquisition model, the blind deconvolution approach is able to provide a ‘corrected’ PSF such that most of the apparent emissions inside the disk of a celestial body during a solar transit can be removed. Let us note that non-parametric techniques cannot outperform the accuracy of parametric methods, however in situations where the telescope imaging model cannot be obtained due to some instrument’s properties (e.g., PICARD/SODISM [MHI14]), 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 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. [ABRS10], allows to efficiently solve the non-convex problem of blind deconvolution with theoretical convergence guarantees. Furthermore, we show the importance of considering multiple observations for the same filter in order to provide a better conditioning to the filter estimation problem.

Acknowledgements

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

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

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

For SDO/AIA, we use the parametric filter estimated by [GYW13]. The validation was performed using a 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 17.

Figure 17: Validating the approximation of the long-range PSF impact by a constant inside the disk of Venus. (left) Logarithm of the filter used for the validation: parametric filter estimated by [GYW13] with the central pixels set to zero inside a window of size , with . (center) Observed image inside a window containing the disk of Venus. (right) Convolved image inside a window containing the disk of Venus. The intensity within the window 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 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 18.

Figure 18: Validating the approximation of the long-range PSF impact by a constant inside the lunar disk. (left) 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 , with . (center) Observed image inside a window containing the lunar disk. (right) Convolved image inside a window containing the lunar disk. The window intensity is almost constant.

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

Appendix B Algorithm initialization

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

b.1 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 modified observations and the filter is given by the delta function . For the subsequent iterations on , the image and the filter are initialized using the value of the previous iteration.

b.2 Regularization parameter

The regularization parameter has an important role in solving the first step in Algorithm 1, since it determines the trade-off between the image regularization and the fidelity to the observations. Several works have studied the choice of this parameter when solving general ill-posed inverse problems. In basis pursuit denoising, [DJ94] proposed a universal value for given by . Since this value increases with the amount of samples (), 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:

(16)

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 is usually not available, the value of in (16) is determined using the modified observations , i.e., , with . Since this value is higher than the actual , it will lead to over-regularized images. Therefore, we propose to refine the value of by an iterative update based on the discrepancy principle between the actual noise energy in (5) and the energy of the residual image (10). By performing a maximum of five updates of the parameter , we assure an appropriate image regularization. The iterative estimation of is summarized in Algorithm 2.

1: ; ; ; ; MaxIter = 5
2:for  to MaxIter do
3:  Images and Filter estimation:
4:     Estimate and using Algorithm 1 with , and
5:  Compute residual and whiteness measure:
6:     
7:     Compute using (12)
8:  Parameter update:
9:     
10:     
11:  Stop when residual is spectrally whiter:
12:     if   then break.      
13:Return and
Algorithm 2 Iterative estimation of

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 [ChP11, CoP11, PB14] for a comprehensive understanding.

Proximal algorithms rely on a key element in convex signal analysis, the proximal operator

(17)

which is uniquely defined for any function , for some  [CoP11]. Proximal algorithms allow the minimization of non-smooth functions such as the -norm and the indicator functions present in the image and filter subproblems.

c.1 First subproblem: Image estimation

In the first subproblem, we are interested in finding the image candidate that minimizes

(18)

a problem that contains a sum of four functions belonging to , for some dimensions . For this, we use the Chambolle-Pock (CP) primal-dual algorithm [ChP11], which is commonly used for solving the minimization of two functions in . In order to take into account the sum of more than two functions, the algorithm is extended as in [GJVA14]. If we set for , for , for and for , the CP iterations are: