Issue 
A&A
Volume 615, July 2018



Article Number  A59  
Number of page(s)  11  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201731765  
Published online  13 July 2018 
Solar hard Xray imaging by means of compressed sensing and finite isotropic wavelet transform
^{1}
Dipartimento di Matematica, Università degli Studi di Genova, Via Dodecaneso 35, 16146 Genova, Italy
email: duvalpoo@dima.unige.it; piana@dima.unige.it
^{2}
CNRSPIN, Via Dodecaneso 33, 16146 Genova, Italy
email: annamaria.massone@cnr.it
Received:
12
August
2017
Accepted:
23
April
2018
Aims. 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 uses the finite isotropic wavelet transform with the Meyer function as the mother wavelet. Furthermore, 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 measurements made with the Spectrometer/Telescope Imaging Xrays (STIX) and experimental observations provided by the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI). The performances of the method are qualitatively validated by comparing some morphological properties of the reconstructed sources with those of the corresponding synthetic configurations. Furthermore, the results concerning experimental data are compared with those obtained by applying other visibilitybased reconstruction methods.
Conclusions. The results show that when the new method is applied to synthetic STIX visibility sets, it provides reconstructions with a spatial accuracy comparable to the accuracy provided by the most popular method in hard Xray solar imaging and with a higher spatial resolution. Furthermore, when it is applied to experimental RHESSI data, the reconstructions are characterized by reliable photometry and by a notable reduction of the ringing effects caused by the instrument point spread function.
Key words: Sun: flares / Sun: Xrays, gamma rays / techniques: image processing
© ESO 2018
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, which are picked up at specific points of the Fourier plane, named (u, ν) 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 2020. 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 also provides a sampling of the Fourier transform of the photon flux in the (u, ν) 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). More recently, novel techniques have been introduced, which rely on a method that has been widely applied in astronomical imaging in the past decade: compressed sensing (Donoho 2006; Candes & Wakin 2008; Bobin et al. 2008). In particular, the present paper describes a possible use of compressed sensing for regularized deconvolution in RHESSI and STIX imaging. In order for compressed sensing to work (Starck & Bobin 2010), we need a sparse representation of the solution of the image reconstruction problem and incoherency between the sampling domain and the sparsitypromoting domain. In our method, we search for sparsity according to a wavelet transformation, which also guarantees a multiscale signal representation. Furthermore, both RHESSI and STIX provide measurements in the Fourier domain, which is incoherent with respect to the wavelet domain.
Most wavelet implementations are associated with multiresolution analysis (MRA; Mallat 1999), mainly because of their computational effectiveness. However, such implementations are far from optimal for applications like filtering and deconvolution, because they are not redundant, meaning that 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 used in radiointerferometry (Li et al. 2011; Garsden et al. 2015), and this for two main reasons: first, it provides redundancy, and second, it is more appropriate for restoring astronomical sources (e.g. stars, galaxies, and 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. For this reason, DuvalPoo et al. (2007) built a 2D isotropic wavelet transform and applied it to two synthetic STIX datasets. Specifically, this 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 such as hard Xray solar imaging. The wavelet system is built in the frequency domain using a 2D isotropic extension of the 1D Meyer mother wavelet (Mallat 1999), but other functions can be used with similar results (Mallat 1999; Portilla & Simoncelli 2000). Following our initial work in DuvalPoo et al. (2007), the present paper provides a systematic formulation of the method comprising both conceptual and technical details, validates the method against five STIX synthetic datasets, computes a quantitative assessment of this analysis, applies the method to several RHESSI experimental observations, and compares the results with the ones provided by other three imaging methods. We also point out that our approach, which we call the finite isotropic wavelet compressed sensing method (VIS_WV), adopts the analysis prior formulation instead of the synthesis prior formulation that is 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 custom Gaussian basis according to the assumption that sources can be represented as linear combinations of a number of Gaussian distributions (Felix et al. 2017). In contrast, VIS_WV proposes a new wavelet implementation that does not rely on a priori information on the specific shapes of the hard Xray sources.
The paper is organized as follows. In Sect. 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 Sect. 5 discusses the results obtained starting from STIX and RHESSI visibilities. Our conclusions are offered in Sect. 6.
2. Hard Xray reconstruction problem
This section provides a quick overview of the data formation model for RHESSI and STIX.
In the energy domain, RHESSI data range from some keV to some MeV with an energy resolution of about 1 keV. The imaging module is composed of 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 of identical grids characterized by a given pitch, different from those 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 (u, ν) plane according to nine circles, as shown in Fig. 1a. It is worth mentioning that the number of samples in each circle is not fixed, but is determined in an optimal way during the computation procedure of the visibilities.
Fig. 1.
Sampling of the (u, ν) plane performed by RHESSI (panel a) and STIX (panel b). For better visualization, panel a shows the sampling provided by seven over nine detectors (those that involve the lowfrequency part of the Fourier domain). 
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. In this way, STIX recording hardware allows sampling the frequency domain in 30 different (u, ν) points placed on spirals, as shown in Fig. 1b. A rigorous description of the data formation process in STIX can be found in Giordano et al. (2015).
As a consequence of these types 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 N × N photon flux image that is to be reconstructed is transformed into the Mdimension vector f ∈ ℝ^{M}, with M = N^{2}, by means of the usual columnwise concatenation procedure. Furthermore, ν ∈ ℂ^{M} is a sparse vector whose nonzero components correspond to the measured visibilities, F ∈ ℂ^{M×M} is the discretized Fourier transform, H is a sparse binary vector that realizes the sampling in the (u, ν) plane according to the way either RHESSI or STIX sample such plane, and · denotes the entrywise product. When we apply the discretized inverse Fourier transform F^{−1} on both sides of Eq. (1), we obtain(2)
where * is the convolution operator. Therefore, the reconstruction of the flux image f from the given visibilities νu can be essentially viewed as a deconvolution problem, where F^{−1}H is the point spread function (PSF) and F^{−1}ν is the dirty map from which the PSF blurring effect must be subtracted.
Fig. 2.
2D isotropic extensions of the Meyer scaling function and mother wavelet with no translation. The functions are represented in the Fourier (top) and space (bottom) domain for three scales. For better visibility, the space domain images are zoomed ×4. 
3. Finite isotropic wavelet transform
This paper introduces two novelties in image reconstruction from visibilities. The first is that a specific wavelet transform is an effective way to implement compressed sensing for hard Xray imaging. Specifically, the reconstruction approach introduced in this paper is a regularization method with an l_{1} penalty term realized by means of a wavelet transform. Therefore we start with describing how the our wavelet transform works.
Let ψ_{a,t} be a family of functions defined by the translation and dilation of a mother wavelet function ψ(x) ∈ L^{2}(ℝ^{2}),(3)
where x = (x_{1}, x_{2}) is a point in ℝ^{2}, and t = (t_{1}, t_{2}) ∈ ℝ^{2} and a ∈ ℝ^{+} are the translation and dilation parameters, respectively. The normalization factor a^{−1/2} ensures that ψ_{a,t} = ψ, where  ·  is the norm in L^{2}(ℝ^{2}). The wavelet transform W(f) of a function f ∈ L^{2}(ℝ^{2}) is defined by Mallat (1999) (4)
where 〈f, ψ_{a,t}〉 is the scalar product in L^{2}(ℝ^{2}), and and are the Fourier transform of f and ψ_{a,t}, respectively. If the admissibility condition is satisfied,(5)
where is the Fourier transform of the mother wavelet ψ, then it is possible to define the inverse wavelet transform as(6)
The second novelty introduced in this work is more technical. Specifically, here we constructed a new specific wavelet transform, named FIWT, using the 2D isotropic extension, in the Fourier domain, of a 1D symmetric function. We started from the 1D Meyer mother function ψ_{M}(x) (Mallat 1999) and constructed the 2D isotropic mother wavelet function in the frequency domain as(7)
Consequently, from Eq. (3) we obtained(8)
Furthermore, in order to span the whole (u, ν) plane, we included in the wavelet framework the scaling function ϕ and constructed it again in the Fourier domain as(9)
where ϕ_{M}(x) is the 1D Meyer scaling function. The FIWT is finally defined from Eq. (4),(10)
where ℱ^{−1} is the inverse Fourier transform. We point out that the wavelet transform defined in Eqs. (7)–(10) is new, although it comes from a very straightforward 2D extension of the 1D formulation.
In the imaging framework the flux distribution f(x) is pixelized at positions(11) (12)
where (c_{1}, c_{2}) is the image center, FOV is the image field of view (FOV), and δ is the pixel dimension in arcsec. With the same reordering as used in Sect. 2, we obtain the vector f = (f_{1}, … , f_{M}) from f((x_{1})_{p}, (x_{2})_{q}) p, q = 0, … , N − 1, where M = N^{2} is the number of pixels of the image to restore and the components correspond to the flux intensity in each pixel. Analogously to what was done in Sect. 2, the vector f is again the unknown of the image reconstruction problem.
In the wavelet transform, we consider j_{0} scales obtained according to the discretization(13)
On the other hand, the discretization of the translation parameter t is made according to(14) (15)
Then, for each scale a_{j} and for each translation ((t_{1})_{n}, (t_{2})_{l}), we construct the Mdimensional vector ψ_{j,n,l} , which corresponds to the pixelization and reordering of the wavelet Eq. (3) for a = a_{j}, t_{1} = (t_{1})_{n}, and t_{2} = (t_{2})_{l}. Analogously, for each translation, we construct the Mdimensional vector ϕ_{n,l}, which is the result of the pixelization and reordering of the scaling function ϕ. This leads to the construction of the set of M · (j_{0} + 1) Mdimensional vectors {u_{i}, i = 1, … , M · (j_{0} + 1)} = { ψ_{j,n,l} j = 0, … , j_{0} − 1; n, l = 0, … , N − 1}⋃} {ϕ_{n,l} n, l = 0, … , N − 1}. By means of a straightforward extension of the proof given for the 1D case in Daubechies (1992) for the Meyer mother function, it can be shown that this set provides a Parseval frame, that is, given f ∈ ℝ^{M}, then(16)
where (·, ·) denotes the canonical inner product in ℝ^{M}. Finally, W is the M · (j_{0} + 1) × M matrix whose rows are made of the vectors u_{i} for i = 1, … , M · (j_{0} + 1). The fact that is a Parseval frame implies that W^{T}W = I and that the forward and inverse discretized FIWT can be written as(18)
respectively, where w is the column vector of dimension M · (j_{0} + 1), whose components are(19)
We conclude by observing that the computational complexity of the FIWT is O(N^{2} log N).
4. Image reconstruction method
The reconstruction of f 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 possible approach in hard Xray imaging is to regularize the inverse problem with the ℓ_{1} norm in some transformation domain. The method we propose in this paper, which we named finite isotropic wavelet transform compressed sensing (VIS_WV), 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, W f_{1}, is designed to penalize an estimate that would not exhibit the sparsity property with respect to FIWT. The regularization parameter, λ > 0, 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 algorithm for the optimization in ℓ_{1} that is widely used in many fields, mainly because of its reliability and rapid convergence. The scheme is illustrated in Table 1 above, where(21)
is the shrinkage operator. Specifically, here FISTA is implemented by imposing the conservation of the flux (line 6 of Table 1, where the subtraction is intended as componentwise) at each iteration. The assessment of the flux is here made by considering the maximum value of the real part of the visibilities provided by detector 9, which is has the smallest radius in the (u, ν) plane: this is the best experimental estimation of the DC component. Furthermore, the FISTA iteration is stopped using a standard optimization rule, which provides k_{opt}.
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 that depends on both the criteria adopted for assessing the quality of the reconstructed image and on the amount of information known about the original image and its noise. Several strategies have been applied in the literature to properly estimate the regularization parameter, of which 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 used most frequently. 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(22)
are known or can be estimated from the dirty map, then the regularization parameter can be chosen as λ = ε/E. In order to satisfy conditions (22), the bounds {ε, E} are estimated by performing the first iteration f_{1} of the FISTA algorithm with λ = 0. The regularization parameter is then set equal to(23)
We compared the performance of the Miller method with that 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; 36–41 keV). More specifically, Figure 3 compares the reconstructions provided by VIS_WV for four values of the regularization parameter: λ = λ_{0} = 0, the value λ = λ_{L} provided by the Lcurve method, the value λ = λ_{M} provided by the Miller method, and the value λ = λ_{1} = 0.1 realizing overregularization. The values of λ_{L} and λ_{M} are very close, and therefore the corresponding reconstructions are very similar, although the Miller method requires less computational effort.
Fig. 3.
July 23, 2002 solar flare image reconstruction by VIS_WV with different choices for λ: panel a: λ = 0 (nonregularized solution), panel b: λ = 3.1 × 10^{−3} as provided by the Lcurve method, panel c: λ = 5.7 × 10^{−3} as provided by the Miller method, and panel d: λ = 0.1 (overregularized solution). The considered time interval is 00:29:10–00:30:19 UT, while the energy channel is 36–41 keV. Detectors from 2 through 9 have been employed. 
Algorithm for the finite isotropic wavelet compressed sensing method (VIS_WV).
5. Experimental results
In this section we experimentally assess the proposed reconstruction algorithm. First, we evaluate our method with STIX synthetic visibilities, where the simulations are performed by means of the STIX Data Processing Software^{1}. Next, VIS_WV is validated against experimental visibilities produced by real flare measurements captured by RHESSI.
5.1. STIX synthetic data
We simulated five flaring events with different configurations. In the images in Fig. 4, panels a, b, and e, the overall photon flux is 10^{4} photons cm^{−2}, but in panel a, the two footpoints have different intensity and different size, in panel b, the brightness and dimension of the two sources are the same, and in panel e, we considered three footpoints where the one to the right of the FOV is more intense. The remaining images (Fig. 4, panels c and d) are characterized by a higher overall intensity (10^{5} photons cm^{−2}) and contain two loops with significantly different curvatures.
Fig. 4.
Reconstructed images for four different STIX simulated flare events. Left column: original source shape. Middle column: reconstructions provided by VIS_CLEAN. Right column: reconstructions provided by VIS_WV. 
Figure 4 compares the images reconstructed by VIS_WV with those provided by VIS_CLEAN (Högbom 1974), and Table 2 contains some physical parameters characterizing the simulated and reconstructed images (in this table the positions and the FWHM are measured in arcsec). These results allow a quantitative comparison between one of the most standard image reconstruction methods for hard Xray astronomy and the new method in terms of spatial accuracy, spatial resolution (assessed through the FWHM) and photometry (assessed through the flux ratio). As far as spatial accuracy is concerned, the behavior of the two methods is quite similar: the reconstruction error associated with VIS_CLEAN ranges from 0.1 to 1.4 arcsec, and the error associated with VIS_WV ranges from 0 to 2 arcsec (although the uncertainties are almost systematically larger in the new algorithm). On the other hand, the resolution power of the new method is undoubtedly better, since it is able to restore the FWHM with a significantly better accuracy than is obtained by VIS_CLEAN (VIS_WV sometimes even overresolves the sources).
5.2. RHESSI experimental measurements
We first considered five flaring events at specific energy channels and time intervals and compared the reconstructions provided by VIS_WV with those obtained by two standard visibilitybased methods, namely VIS_CLEAN (Högbom 1974) and UV_SMOOTH (Massone et al. 2009), and a newly introduced compressed sensing method using a custom basis of predefined Gaussian images (Felix et al. 2017; this method is called VIS_CS in the following). For all experiments we used RHESSI detectors from 2 to 9 and a threescale decomposition for the waveletbased deconvolution methods. In particular, we considered flaring events characterized by rather different morphologies such as 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). For all cases, we used notconjugated, normalized, and edited visibilities with natural weighting, with at least 6 roll bins and 64 roll bins at most.
The results in Fig. 5 show that sparsity promotion in the two compressed sensing methods reduces the artifacts and provides a higher spatial resolution. However, VIS_CS seems to systematically shrink the reconstructed sources according to a behavior rather similar to that of MEM_NJIT (Bong et al. 2006; Massone et al. 2009) and particularly significant in the event on April 15, 2002.
Fig. 5.
Reconstructed images from different flaring events using three reconstruction methods: VIS_CLEAN (left), UV_SMOOTH (center left), VIS_CS (center right), and VIS_WV (right). RHESSI detectors from 2 to 9 have been used to generate the visibilities. 
Figures 6, 7, and 8 focus on datasets acquired by RHESSI during the event on July 23, 2002, and show the potential of the new method well. We are aware that the application of imaging methods to this flare, which has been extensively studied even in a special 2003 issue of “The Astrophysical Journal Letters” may be critical because this event is sufficiently intense to assume that pulse pileup may occur at intermediate energies. However, this issue is not critical for this present analysis, which is devoted to introduce the potentialities and limitations of a newly introduced imaging algorithm. The results shown in these figures agree remarkably closely with those obtained in several papers, in particular, by Emslie et al. (2003) and by Massone et al. (2009). Figure 6 reproduces the same analysis as was performed by Emslie et al. (2003) using VIS_CLEAN and clearly points out how VIS_WV better preserves the source morphologies along the energy increase, reduces the artifacts in between the different sources, and maintains the image reliability particularly at high energies. Figure 7 shows the total flux extracted from the images in Fig. 6 and compares it with the flux estimated from the visibilities collected by detector 9. The photometric behavior of VIS_WV is compared with the behaviors associated with two other visibilitybased methods (VIS_CLEAN and VIS_CS) and to two methods using count modulation profiles as input data, EM (Benvenuto et al. 2013) and PIXON (Metcalf et al. 1996). This comparison suggests that VIS_WV guarantees a rather good photometry at least for this specific analysis. On the other hand, Fig. 8 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 has been argued by Massone et al. (2009), but in agreement with the results in Emslie et al. (2003).
Fig. 6.
Images reconstructed by VIS_WV for the July 23, 2002, flaring event through different energy ranges (00:29:10–00:30:19 UT). Visibilities collected by RHESSI collimators from 2 to 9 have been used for the reconstructions. 
Fig. 7.
Total flux extracted from the reconstructions in Fig. 6 plotted vs. energy and compared to the flux estimated from visibilities collected by detector 9 and to the flux extracted from the reconstructions provided by other methods. Left panel: comparison with visibilitybased methods. Right panel: comparison with methods using count modulation profiles. 
Fig. 8.
Images reconstructed by VIS_WV for the July 23, 2002, flaring event through different time intervals (36–41 keV energy channel). Visibilities collected by RHESSI collimators from 2 to 9 have been used for the reconstructions. 
5.3. Computational performance
The experiments in this paper have been performed by means of a secondgeneration 4core Intel i72600 (3.40 GHz) CPU using IDL 8.4 on Ubuntu 16.04. The computational time complexity of VIS_WV is O(KN^{2} log N), where N is the image size and K is the number of iterations required by FISTA for solving the optimization problem. We did not include the number of employed scales j_{0} for FIWT in the time complexity analysis, since j_{0} is usually a small constant number.
Figure 9 shows the computational performance of VIS_WV under different configurations in the case of the first STIX simulation (Fig. 4a). Figure 9a shows the linear relation between the employed time and the number of iterations K for the reconstruction of an image map of size N = 128 and considering j_{0} = 3. On the other hand, Fig. 9b shows the increase in computational time when VIS_WV recovers images with increasingly greater size, where in this case, the number of iterations and considered scales are fixed at K = 30 and j_{0} = 3, respectively. Since most reconstructions are performed with an image size of 128 pixels and our method usually solves the reconstruction problem with fewer than 100 iterations, we can conclude that the average time required by VIS_WV for reconstructing a single standard image is between 1.5 and 2.5 s.
Fig. 9.
Computational performance of VIS_WV. The computational time is compared versus panel a the number of FISTA iterations and panel b the size of the reconstructed image. The test is made using the synthetic STIX visibilities used in Fig. 4a. 
6. Conclusions
Hard Xray solar space telescopes typically provide their experimental measurements in the form of visibilities, that is, sparse samples of the Fourier transform of the incoming flux. Producing images in this setting therefore requires applying reconstruction methods that realize the inversion of the Fourier transform from limited data. Several methods have been used so far, but just recently, some of them have explored the possibility of exploiting compressed sensing, which is 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 listed below.
It relies on a isotropic wavelet transform because Xray sources are either isotropic or characterized by a slow change of shape.
It avoids the use of a computationally demanding catalogbased 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.
We point out that hard Xray sources are not properly isotropic, and we might therefore consider other types of multiscale decomposition, like shearlets, that account for anisotropies in the image (DuvalPoo et al. 2015). However, this approach would lead to a computationally more demanding algorithm, and since the Xray flaring sources are characterized by rather smooth changes in directionality, we assumed that an isotropic transformation should be effective. The accuracy of the reconstructions confirmed the reliability of this assumption. We are currently investigating the improvement that the use of shearlets would provide to the reconstructions.
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. Furthermore, the method has a notable automation degree, since the only parameter to optimize is the regularization parameter (the number of scales is fixed a priori) and the computational burden required is low. Therefore, this approach is expected to be competitive for an extended use on a large amount of data. The implementation of the algorithm within Solar SoftWare is completed^{2}, and this will allow the systematic use of this approach on RHESSI observations and future application on STIX measurements.
Acknowledgments
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, M. J., & Schmahl, E. 2002, Sol. Phys., 210, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Beck, A., & Teboulle, M. 2009, SIAM J. Imag. Sci., 2, 183 [CrossRef] [Google Scholar]
 Benvenuto, F., Schwartz, R., Piana, M., & Massone, A. M. 2013, A&A, 555, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benz, A., Krucker, S., & Hurford, G., et al. 2012, in Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray, Proc. SPIE, 8443, 84433L [CrossRef] [Google Scholar]
 Bertero, M., & Boccacci, P. 1998, Introduction to Inverse Problems in Imaging (Boca Raton: CRC Press) [Google Scholar]
 Bobin, J., Starck, J. L., & Ottensamer, R. 2008, IEEE J. Sel. Top. Signal Process., 2, 718 [NASA ADS] [CrossRef] [Google Scholar]
 Bong, S.C., Lee, J., Gary, D. E., & Yun, H. S. 2006, ApJ, 636, 1159 [NASA ADS] [CrossRef] [Google Scholar]
 Candes, E. J., & Wakin, M. B. 2008, Signal Process. Mag., 25, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Cornwell, T. J., & Evans, K. 1985, A&A, 143, 77 [Google Scholar]
 Daubechies, I. 1992, Ten Lectures on Wavelets (Philadelphia: SIAM), 61 [Google Scholar]
 Donoho, D. L. 2006, IEEE Trans. Inf. Theory, 52, 1289. [Google Scholar]
 DuvalPoo, M. A., Massone, A. M., & Piana, M. 2007, in 2017 International Conference on SampTA, IEEE, 677 [Google Scholar]
 DuvalPoo, M. A., Odone, F., & De Vito, E. 2015, IEEE Trans. Image Process., 24, 3768 [CrossRef] [Google Scholar]
 Emslie, A. G., Kontar, E. P., Krucker, S., & Lin, R. P. 2003, ApJ, 595, L107 [NASA ADS] [CrossRef] [Google Scholar]
 Engl, H. W., Hanke, M., & Neubauer, A. 1996, Regularization of Inverse Problems (Dordrecht: Springer Science & Business Media), 375 [Google Scholar]
 Felix, S., Bolzern, R., & Battaglia, M. 2017, ApJ, 849, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Garsden, H., Girard, J. M., Starck, J.L., et al. 2015, A&A, 575, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Giordano, S., Pinamonti, N., Piana, M., & Massone, A. M. 2015, SIAM J. Imaging Sci., 8, 1315 [Google Scholar]
 Golub, G. H., Heath, M., & Wahba, G. 1979, Technometrics, 21, 215 [CrossRef] [Google Scholar]
 Hansen, P. C. 1992, SIAM Rev., 34, 561 [CrossRef] [MathSciNet] [Google Scholar]
 Häuser, S., & Steidl, G. 2014, ArXiv eprints [arXiv:1202.1773v2] [Google Scholar]
 Högbom, J. 1974, A&AS, 15, 417 [Google Scholar]
 Hurford, G. J., Schmahl, E. J., Schwartz, R. A., et al. 2002, Sol. Phys., 210, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Li, F., Cornwell, T. J., & de Hoog, F. 2011, A&A, 528, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lin, R., Dennis, B., Hurford, G., et al. 2002, Sol. Phys., 210, 3 [Google Scholar]
 Mallat, S. 1999, A Wavelet Tour of Signal Processing (Waltham: Acaic Press) [Google Scholar]
 Massone, A. M., Emslie, A. G., Hurford, G. J., et al. 2009, ApJ, 703, 2004 [NASA ADS] [CrossRef] [Google Scholar]
 Metcalf, T. R., Hudson, H. S., Kosugi, T., Puetter, R. C., & Pina, R. K. 1996, ApJ, 466, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, K. 1970, SIAM J. Math. Anal., 1, 52 [CrossRef] [MathSciNet] [Google Scholar]
 Morozov, V. A. 1984, Methods for Solving Incorrectly Posed Problems (New York: Springer) [Google Scholar]
 Oster, G., Wasserman, M., & Zwerling, C. 1964, JOSA, 54, 169 [CrossRef] [Google Scholar]
 Portilla, J., & Simoncelli, E. P. 2000, Int. J. Comput. Vis., 40, 49 [CrossRef] [Google Scholar]
 Starck, J.L., & Bobin, J. 2010, Proc. IEEE, 98, 1021 [CrossRef] [Google Scholar]
 Starck, J.L., Fadili, J., & Murtagh, F. 2007, IEEE Trans. Image Process., 16, 297 [Google Scholar]
All Tables
All Figures
Fig. 1.
Sampling of the (u, ν) plane performed by RHESSI (panel a) and STIX (panel b). For better visualization, panel a shows the sampling provided by seven over nine detectors (those that involve the lowfrequency part of the Fourier domain). 

In the text 
Fig. 2.
2D isotropic extensions of the Meyer scaling function and mother wavelet with no translation. The functions are represented in the Fourier (top) and space (bottom) domain for three scales. For better visibility, the space domain images are zoomed ×4. 

In the text 
Fig. 3.
July 23, 2002 solar flare image reconstruction by VIS_WV with different choices for λ: panel a: λ = 0 (nonregularized solution), panel b: λ = 3.1 × 10^{−3} as provided by the Lcurve method, panel c: λ = 5.7 × 10^{−3} as provided by the Miller method, and panel d: λ = 0.1 (overregularized solution). The considered time interval is 00:29:10–00:30:19 UT, while the energy channel is 36–41 keV. Detectors from 2 through 9 have been employed. 

In the text 
Fig. 4.
Reconstructed images for four different STIX simulated flare events. Left column: original source shape. Middle column: reconstructions provided by VIS_CLEAN. Right column: reconstructions provided by VIS_WV. 

In the text 
Fig. 5.
Reconstructed images from different flaring events using three reconstruction methods: VIS_CLEAN (left), UV_SMOOTH (center left), VIS_CS (center right), and VIS_WV (right). RHESSI detectors from 2 to 9 have been used to generate the visibilities. 

In the text 
Fig. 6.
Images reconstructed by VIS_WV for the July 23, 2002, flaring event through different energy ranges (00:29:10–00:30:19 UT). Visibilities collected by RHESSI collimators from 2 to 9 have been used for the reconstructions. 

In the text 
Fig. 7.
Total flux extracted from the reconstructions in Fig. 6 plotted vs. energy and compared to the flux estimated from visibilities collected by detector 9 and to the flux extracted from the reconstructions provided by other methods. Left panel: comparison with visibilitybased methods. Right panel: comparison with methods using count modulation profiles. 

In the text 
Fig. 8.
Images reconstructed by VIS_WV for the July 23, 2002, flaring event through different time intervals (36–41 keV energy channel). Visibilities collected by RHESSI collimators from 2 to 9 have been used for the reconstructions. 

In the text 
Fig. 9.
Computational performance of VIS_WV. The computational time is compared versus panel a the number of FISTA iterations and panel b the size of the reconstructed image. The test is made using the synthetic STIX visibilities used in Fig. 4a. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.