Solar hard Xray imaging by means of Compressed Sensing and
Finite Isotropic Wavelet Transform
Key Words.:
Sun: flares – Sun: Xrays, gamma rays – Techniques: image processingAbstract
Context:
Aims:This paper shows that compressed sensing realized by means of regularized deconvolution and the Finite Isotropic Wavelet Transform is effective and reliable in hard Xray solar imaging.
Methods:The method utilizes the Finite Isotropic Wavelet Transform with Meyer function as the mother wavelet. Further, compressed sensing is realized by optimizing a sparsitypromoting regularized objective function by means of the Fast Iterative ShrinkageThresholding Algorithm. Eventually, the regularization parameter is selected by means of the Miller criterion.
Results:The method is applied against both synthetic data mimicking the Spectrometer/Telescope Imaging Xrays (STIX) measurements and experimental observations provided by the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI). The performances of the method are compared with the results provided by standard visibilitybased reconstruction methods.
Conclusions:The results show that the application of the sparsity constraint and the use of a continuous, isotropic framework for the wavelet transform provide a notable spatial accuracy and significantly reduce the ringing effects due to the instrument point spread functions.
1 Introduction
Imaging spectroscopy is a powerful tool for exploring the physics underlying particle acceleration and transport in solar flares. In order to realize imaging spectroscopy in the hard Xray range, in 2002 NASA launched the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI) (Lin et al. 2002), whose data have resulted in hard Xray images of unprecedented angular and energy resolution. In a nutshell, RHESSI rotating collimators modulate the Xray flux coming from the Sun, providing as a result sparse samples of its Fourier transform, named visibilities, picked up at specific points of the Fourier plane, named plane in this context.
The Spectrometer/Telescope for Imaging Xrays (STIX) (Benz et al. 2012) is one of the ten instruments in the payload of Solar Orbiter, which will be launched by ESA close to the Sun in 2018. Analogously to RHESSI the main goal of STIX is to measure hard Xray photons emitted during solar flares in order to determine the intensity, spectrum, timing and location of accelerated electrons near the Sun. The imaging system that characterizes this device relies on the Moiré pattern concept (Oster et al. 1964) and, similarly to RHESSI, it provides as well a sampling of the Fourier transform of the photon flux in the plane (Giordano et al. 2015).
For both RHESSI and STIX, image reconstruction is needed to determine the actual spatial photon flux distribution from the few Fourier components acquired by the hard Xray collimators and several methods have been realized to this goal (Högbom 1974; Cornwell & Evans 1985; Aschwanden & Schmahl 2002; Bong et al. 2006; Massone et al. 2009), but none of them exploits a methodology that has been widely applied in astronomical imaging in the last decade, i.e. compressed sensing (Donoho 2006; Candes & Wakin 2008; Bobin et al. 2008). The present paper describes a possible use of compressed sensing for regularized deconvolution in RHESSI and STIX imaging. In order to work, compressed sensing requires data incoherency and a sparse representation of the solution of the image reconstruction problem. Both RHESSI and STIX sample the Fourier domain in a way characterized by a notable level of incoherency; on the other hand, a typical hard Xray image configuration is made of few, simple shapes (mainly, footpoints and loops) and therefore it is straightforward to represent it as the superposition of a small number of basis functions. It is wellknown that an advantageous approach to realize compressed sensing is to use a wavelet transform since wavelets can provide, in addition to compression, a multiscale signal representation.
Most wavelet implementations are associated to multiresolution analysis (MRA) (Mallat 1999), mainly because of their computational effectiveness. However such implementations are far from optimal for applications like filtering and deconvolution, owing to the fact that they are not redundant, i.e., the dimension of the image decreases at coarser scales (Starck et al. 2007). As an alternative to MRA, the Isotropic Undecimated Wavelet Transform (IUWT) is rather often utilized in radiointerferometry (Li et al. 2011; Garsden et al. 2015) and this for two main reasons: first, it provides redundancy and, second, it is better appropriate for restoring astronomical sources (e.g.: stars, galaxies, flares), which are mostly isotropic or quasiisotropic.
IUWT relies on a discrete wavelet transform path, which is not fully appropriate when the input data are provided by a rather sparse sampling of the frequency domain, as in visibilitybased hard Xray imaging. Therefore in the present paper we built a 2D isotropic wavelet transform that follows the continuous wavelet transform path. Specifically, our Finite Isotropic Wavelet Transform (FIWT) is inspired by the shearlet transform implementation proposed by Häuser & Steidl (2014) and is a redundant transform, which can be effectively applied in deconvolution problems like the RHESSI and STIX ones. This wavelet system is built in the frequency domain by using a 2D isotropic extension of the 1D Meyer mother wavelet (Mallat 1999) but other functions can be used with comparable results (Mallat 1999; Portilla & Simoncelli 2000). In order to reduce the numerical instability due to the illposedness of the deconvolution problem, in this paper we formulated a sparsityenhancing regularized version of the FIWT multiscale sparse decomposition, which we called the Finite Isotropic waVElet Compressed Sensing method (5CS) and we used it to obtain reconstructions of hard Xray images from synthetic STIX and experimental RHESSI data. We also point out that our approach adopts the analysis prior formulation instead of the synthesis prior formulation followed in the case of radiointerferometry visibilities (Li et al. 2011).
Finally, it is worth noting that an alternative way to realize compressed sensing in imaging is to use a dictionary made of few shapes, replicated many times for different scales and positions, and then to realize sparsity with respect to this dictionary. However, 5CS realizes compressed sensing utilizing a continuous wavelet transformation and therefore it does not need any catalogue of basis images to work. This has two advantages: first, the construction of a catalogue requires to know in advance all possible source shapes and, second and more importantly, cataloguebased compressed sensing is computationally more demanding.
The paper is organized as follows. In Section 2 we formally set up the imaging problem for RHESSI and STIX. Section 3 defines the FIWT and its implementation. Section 4 describes the image reconstruction method based on FIWT and compressed sensing and Section 5 discusses the results obtained starting from STIX and RHESSI visibilities. Our conclusions are offered in Section 6.
2 The hard Xray reconstruction problem
This section provides a quick overview of the model of data formation for RHESSI and STIX.
In the energy domain, RHESSI data range from some keV to some MeV with energy resolution of around 1 keV. The imaging module is composed by nine pairs of Rotating Modulation Collimators (RMCs), each one formed by a pair of equally spaced fine grids, placed in front of the detecting device. Each pair is composed by identical grids characterized by a given pitch, different from the ones characterizing all the other pairs of grids. The rotational motion of the spacecraft around its own axis causes a periodic modulation of the incident flux. As a result (Hurford et al. 2002) the instrument samples the plane according to nine circles, as shown in Figure 1 (a). It is worth mentioning that the number of samples in each circle is not fixed but determined in an optimal way during the computation procedure of the visibilities.
STIX is formed by 30 detectors recording Xray photons in the range 4 – 150 keV. On each detector, the incident flux is modulated by means of a subcollimator formed by two distant grids with slightly different pitches and slightly different orientations. The effect of this grid configuration is to create the superposition of two spatial modulations, named Moiré pattern. The recording process on the detector associated with each Moiré pattern provides a specific STIX visibility. Therefore STIX recording hardware allows sampling the frequency domain in 30 different points placed on spirals as shown in Figure 1 (b). A rigorous description of the data formation process in STIX can be found in Giordano et al. (2015).
As a consequence of these kinds of hardware design, the mathematical model for data formation in RHESSI and STIX can be formulated in a matrix form as
(1) 
where the original photon flux image to reconstruct is lexicographically reordered to define the vector , with . Further, is a sparse vector whose nonzero components correspond to the measured visibilities, is the discretized Fourier transform, is a sparse binary matrix that realizes the sampling in the plane, and denotes the entrywise product. If we apply the discretized inverse Fourier transform on both sides of (1) we obtain
(2) 
where is the convolution operator. Therefore, the reconstruction of the flux image from the given visibilities can be essentially viewed as a deconvolution problem, where is the point spread function and is the dirty map from which the PSF blurring effect must be subtracted.
3 The Finite Isotropic Wavelet Transform
Let be a family of functions defined by the translation and dilation of a mother wavelet function , i.e.
(3) 
where is a point in , and are the translation and dilation parameters, respectively. The normalization factor ensures that , where is the norm in . The wavelet transform of a function is defined by (Mallat 1999)
(4) 
where is the scalar product in and and are the Fourier transform of and , respectively. If the admissibility condition is satisfied, i.e.
(5) 
then it is possible to define the inverse wavelet transform as
(6) 
We constructed a new specific wavelet transform, named the Finite Isotropic Wavelet Transform (FIWT), using the 2D isotropic extension, in the Fourier domain, of a 1D symmetric function. Specifically, we started from the 1D Meyer mother function (Mallat 1999) and constructed the 2D isotropic mother wavelet function as
(7) 
Consequently, from (3) we obtained
(8) 
Further, in order to span the whole plane, we included in the wavelet framework the scaling function and constructed it again in the Fourier domain as
(9) 
where is the 1D Meyer scaling function. The FIWT is finally defined from (3), i.e.
(10) 
where is the inverse Fourier transform.
In the imaging framework the flux distribution is pixelized at positions
(11) 
(12) 
where is the image center, is the image fieldofview and is the pixel dimension in arcsec. With the same lexicographical reordering used in section 2 we obtain the vector from , where is the number of pixels of the image to restore and the components correspond to the flux intensity in each pixel. Accordingly to what done is Section 2, the vector is again the unknown of the image reconstruction problem.
In the wavelet transformation we will consider scales obtained according to the discretization
(13) 
On the other hand, the discretization of the traslation parameter is made according to
(14) 
(15) 
Then, for each scale and for each translation , we construct the dimensional vector which corresponds to the pixelization and reordering of the wavelet (3) for , , and . Analogously, for each translation we construct the dimensional vector , which is the result of the pixelization and reordering of the scaling function . This leads to the construction of the set of dimensional vectors . This set provides a Parseval frame, i.e., given , then
(16) 
and
(17) 
where denotes the canonical inner product in . Finally, is the matrix whose rows are made of the vectors for . The fact that is a Parseval frame implies that and that the forward and inverse discretized FIWT can be written as
(18) 
respectively, where is the column vector of dimension whose components are
(19) 
We conclude by observing that the computational complexity of the FIWT is around .
4 The image reconstruction method
The reconstruction of from is an illposed problem (Bertero & Boccacci 1998) and therefore some prior knowledge about the image is required to regularize the reconstruction problem (Engl et al. 1996). A valid approach in hard Xray imaging is to regularize the inverse problem with the norm in some transformation domain. In fact, the method we propose in this paper, which we named Finite Isotropic waVElet transform Compressed Sensing (5CS), addresses the optimization problem
(20) 
The data term of the objective function to minimize, , quantifies the prediction error with respect to the measurements. The regularization term, , is designed to penalize an estimate that would not exhibit the sparsity property with respect to FIWT. The regularization parameter, , provides a tradeoff between fidelity to the measurements and sparsity.
Problem (20) can be numerically solved by means of any algorithm for nonlinear optimization. Here we used the Fast Iterative ShrinkageThresholding Algorithm (FISTA) (Beck & Teboulle 2009), which is an solver widely used in many fields mainly because of its reliability and rapid convergence. Specifically, here FISTA is implemented by imposing the conservation of the flux at each iteration and by utilizing a standard rule for optimally stopping the iterations.
The only open parameter for method (20) is the regularization parameter , which plays a crucial role in the reconstruction process. Finding an appropriate value for is often not a trivial problem, depending on both the criteria adopted for assessing the quality of the reconstructed image and the amount of information known about the original image and its noise. There exist several strategies in the literature to properly estimate the regularization parameter, where the discrepancy principle (Morozov 1984), the Miller method (Miller 1970), the generalized crossvalidation (GCV) method (Golub et al. 1979), and the Lcurve method (Miller 1970; Hansen 1992) are the most used ones. In this work, we chose the Miller method (Miller 1970) because of its simplicity and because it is not computationally demanding. According to this selection criterion, if the following bounds
(21) 
are known, or can be estimated from the dirty map, then the regularization parameter can be chosen as . In order to satisfy conditions (21), the bounds are estimated performing the first iteration of the FISTA algorithm with . The regularization parameter is then set equal to
(22) 
We compared the performance of the Miller method with the one of the Lcurve criterion in the case of the reconstruction of an image of the July 23 2002 event (00:29:10  00:30:19 UT; keV). More specifically, Figure 3 compares the reconstructions provided by 5CS for four values of the regularization parameter: , the value provided by the Lcurve method, the value provided by the Miller method, and the value realizing overregularization. The values of and are very close and therefore the corresponding reconstructions are very similar, although the Miller method requires less computational effort.
5 Experimental results
In this Section we performed an experimental assessment of the proposed reconstruction algorithm. First, we evaluated our method with STIX synthetic visibilities where the simulations are performed by means of the STIX Data Processing Software^{1}^{1}1https://stix.cs.technik.fhnw.ch/. Next, 5CS was validated against experimental visibilities produced by real flare measurements captured by RHESSI.
5.1 The case of Stix synthetic data




We simulated four flaring events with different configurations. In the first two images (Figure 4, (a) and (b)) the overall photon flux is photons/cm, but in the first image the two footpoints have different intensity and different size, while in the second one the brightness and dimension of the two sources are the same. The other two images (Figure 4, (c) and (d)) are characterized by a higher overall intensity ( photons/cm) and contain two loops with significantly different curvatures.
Figure 4 compares the images reconstructed by 5CS with the ones provided by CLEAN (Högbom 1974) and Table 1 contains some physical parameters characterizing the simulated and reconstructed images (in this Table the positions and the full width at half maximum (FWHM) are measured in arcsec). These results imply that 5CS and CLEAN recover the physical parameters with comparable accuracy, although 5CS is able to significantly reduce the impact of spurious artifacts.
Simulated Map  CLEAN  5CS  
Simulated Double Footpoint Flare (SDFF1)  
X  20.0  20.4 0.9  19.2 1.5  
Y  20.0  18.9 0.9  21.8 1.1  
{sideways}Peak 1  FWHM  15.0  9.6 3.4  9.0 4.4 
X  5.0  4.0 0.7  5.0 1.2  
Y  5.0  4.8 0.6  3.5 2.9  
{sideways}Peak 2  FWHM  10.0  12.8 0.8  11.9 1.6 
Flux Ratio  2.00  1.41 0.10  1.78 0.42  
Simulated Double Footpoint Flare (SDFF2)  
X  30.0  28.9 0.6  29.3 0.8  
Y  0.0  1.0 0.7  0.4 0.7  
{sideways}Peak 1  FWHM  10.0  13.2 0.9  9.6 0.8 
X  30.0  29.5 0.5  28.0 0.8  
Y  0.0  0.1 0.6  1.4 0.5  
{sideways}Peak 2  FWHM  10.0  12.9 0.6  9.3 0.6 
Flux Ratio  1.00  0.99 0.10  1.02 0.09  
Simulated Loop Flare (SLF1)  
X  10.0  9.0 0.5  9.6 0.5  
{sideways}Peak  Y  15.0  13.6 0.5  14.6 0.7 
Simulated Loop Flare (SLF2)  
X  10.0  9.8 0.4  8.7 0.7  
{sideways}Peak  Y  0.0  0.1 0.3  0.4 0.7 
5.2 The case of Rhessi experimental measurements
We first considered five flaring events at specific energy channels and time intervals and compared the reconstructions provided by 5CS with the ones obtained by two standard visibilitybased methods, namely CLEAN (Högbom 1974) and uv smooth (Massone et al. 2009). In particular we have considered flaring events characterized by rather different morphologies like double footpoints (20 February 2002), loops (15 April 2002; 13 May 2013), extended sources (31 August 2004), and extended plus compact sources (2 December 2003). The results in Figure 5 show that sparsity promotion and the use of a continuous wavelet formulation reduce the artifacts and provide a higher spatial resolution. For all experiments we used RHESSI detectors from 2 to 9 and a 3scale decomposition for the waveletbased deconvolution methods. Figures 6 and 7 focus on datasets acquired by RHESSI during the July 23 2002 event. In particular, Figure 6 reproduces the same analysis performed by Emslie et al. (2003) using CLEAN and clearly points out how 5CS better preserves the sources’ morphology along the energy increase, reduces the artifacts in between the different sources and maintains the image reliability particularly at high energies. On the other hand, Figure 7 shows the time evolution of the flaring emission at a fixed energy channel and probably seems to reject the presence of emission along a curved locus joining the northern and souther sources, in contrast to what argued by Massone et al. (2009), but accordingly to the results in (Emslie et al. 2003).





5.3 Computational performance
The experiments in this paper have been performed by means of a second generation 4 core Intel i72600 (3.40 GHz) CPU using IDL 8.4 on Ubuntu 16.04. The computational time complexity of 5CS is near where is the image size and is the number of iterations required by FISTA for solving the optimization problem. We did not include the number of employed scales for FIWT on the time complexity analysis, since is usually a small number.
Figure 8 shows the computational performance of 5CS under different configurations in the case of the first STIX simulation (Figure 4(a)). In Figure 8 (a) we can observe the linear relation between the employed time and the number of iterations for the reconstruction of an image map of size and considering . On the other hand, Figure 8 (b) shows how the computational time increases when 5CS recovers images with higher and higher size, where in this case, the number of iterations and considered scales are fixed at and , respectively. Since most reconstructions are performed with an image size of pixels and our method usually solves the reconstruction problem with less than 100 iterations, we can conclude that the average time required by 5CS for reconstructing a single standard image is between 1.5 and 2.5 seconds.
6 Conclusions
Hard Xray solar space telescopes typically provide their experimental measurements in the form of visibilities, i.e., sparse samples of the Fourier transform of the incoming flux. Therefore producing images in this setting requires the application of reconstruction methods that realize the inversion of the Fourier transform from limited data. Several methods have been utilized so far but none of them explicitly exploited compressed sensing, i.e. the use of a sparsityenhancing penalty term in regularization. Here we introduced a waveletbased deconvolution method promoting sparsity for hard Xray image reconstruction from visibilities. The main aspects of this method are that

It relies on a continuous isotropic wavelet transform, coherently to the fact that Xray sources are either isotropic or characterized by a slow change of shape.

It avoids the use of a computationally demanding cataloguebased compressed sensing.

It realizes regularization by means of optimization of a minimum problem where the penalty term promotes sparsity.

It realizes numerical optimization by means of a fast iterative algorithm.
The applications against both synthetic STIX and experimental RHESSI visibilities show the reliability of the method in terms of both the spatial resolution achieved and the reduction of spurious artifacts. Finally, the computational burden required by the method is low and competitive with respect to possible big data applications. The implementation of the algorithm within Solar SoftWare, which is under construction, will allow the systematic use of this approach against RHESSI observations and future application against STIX measurements.
Acknowledgements.
We would like to acknowledge Federico Benvenuto for useful discussions. This research work has been supported by the ASI/INAF grant ”Solar Orbiter ILWS  Supporto scientifico per la realizzazione degli strumenti METIS e SWA/DPU nelle fasi B2C1”References
 Aschwanden & Schmahl (2002) Aschwanden, M. J. & Schmahl, E. 2002, Solar Physics, 210, 193
 Beck & Teboulle (2009) Beck, A. & Teboulle, M. 2009, SIAM journal on imaging sciences, 2, 183
 Benz et al. (2012) Benz, A., Krucker, S., Hurford, G., et al. 2012, in SPIE Astronomical Telescopes+ Instrumentation, International Society for Optics and Photonics, 84433L–84433L
 Bertero & Boccacci (1998) Bertero, M. & Boccacci, P. 1998, Introduction to inverse problems in imaging (CRC press)
 Bobin et al. (2008) Bobin, J., Starck, J. L., & Ottensamer, R. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 718
 Bong et al. (2006) Bong, S. C., Li, L. J., Gary, D. E., & Yun, H. S. 2006, Astrophysical Journal, 636, 1159
 Candes & Wakin (2008) Candes, E. J. & Wakin, M. B. 2008, Signal processing magazine, 25, 21
 Cornwell & Evans (1985) Cornwell, T. J. & Evans, K. 1985, Astronomy and Astrophysics, 143, 77
 Donoho (2006) Donoho, D. L. 2006, IEEE Transactions on Information Theory, 52, 1289
 Emslie et al. (2003) Emslie, A. G., Kontar, E. P., Krucker, S., & Lin, R. P. 2003, The Astrophysical Journal, 595, L107
 Engl et al. (1996) Engl, H. W., Hanke, M., & Neubauer, A. 1996, Regularization of inverse problems, Vol. 375 (Springer Science & Business Media)
 Garsden et al. (2015) Garsden, H., Girard, J. M., L, S. J., & S, C. 2015, Astronomy & Astrophysics, 575, A90
 Giordano et al. (2015) Giordano, S., Pinamonti, N., Piana, M., & Massone, A. M. 2015, SIAM Journal on Imaging Sciences, 8, 1315
 Golub et al. (1979) Golub, G. H., Heath, M., & Wahba, G. 1979, Technometrics, 21, 215
 Hansen (1992) Hansen, P. C. 1992, SIAM review, 34, 561
 Häuser & Steidl (2014) Häuser, S. & Steidl, G. 2014, ArXiv
 Högbom (1974) Högbom, J. 1974, Astronomy and Astrophysics Supplement Series, 15, 417
 Hurford et al. (2002) Hurford, G. J., Schmahl, E. J., Schwartz, R. A., et al. 2002, Solar Physics, 210, 61
 Li et al. (2011) Li, F., Cornwell, T. J., & de Hoog, F. 2011, Astronomy & Astrophysics, 528, A31
 Lin et al. (2002) Lin, R., Dennis, B., Hurford, G., et al. 2002, Solar Physics, 210, 3
 Mallat (1999) Mallat, S. 1999, A wavelet tour of signal processing (Acaic press)
 Massone et al. (2009) Massone, A. M., Emslie, A. G., Hurford, G., et al. 2009, The Astrophysical Journal, 703, 2004
 Miller (1970) Miller, K. 1970, SIAM Journal on Mathematical Analysis, 1, 52
 Morozov (1984) Morozov, V. A. 1984, Methods for solving incorrectly posed problems (Springer New York)
 Oster et al. (1964) Oster, G., Wasserman, M., & Zwerling, C. 1964, JOSA, 54, 169
 Portilla & Simoncelli (2000) Portilla, J. & Simoncelli, E. P. 2000, International journal of computer vision, 40, 49
 Starck et al. (2007) Starck, J.L., Fadili, J., & Murtagh, F. 2007, IEEE Transactions on Image Processing, 16, 297