© 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 X-ray range, in 2002, NASA launched the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI; Lin et al. 2002), whose data have resulted in hard X-ray images of unprecedented angular and energy resolution. In a nutshell, RHESSI rotating collimators modulate the X-ray 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 X-rays (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 X-ray 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 X-ray 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 sparsity-promoting domain. In our method, we search for sparsity according to a wavelet transformation, which also guarantees a multi-scale 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 radio-interferometry (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 quasi-isotropic.

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 visibility-based hard X-ray imaging. For this reason, Duval-Poo 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 X-ray 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 Duval-Poo 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 radio-interferometry 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 X-ray 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 X-ray 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 low-frequency part of the Fourier domain).

STIX is formed by 30 detectors recording X-ray photons in the range 4–150 keV. On each detector, the incident flux is modulated by means of a sub-collimator 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$$\mathit{H}\cdot \mathbf{F}\mathit{f}=\mathit{\upsilon },$$(1)

where the original N × N photon flux image that is to be reconstructed is transformed into the M-dimension vector f ∈ ℝ^M, with M = N², by means of the usual column-wise concatenation procedure. Furthermore, ν ∈ ℂ^M is a sparse vector whose non-zero 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 entry-wise product. When we apply the discretized inverse Fourier transform F⁻¹ on both sides of Eq. (1), we obtain$${\mathbf{F}}^{-1}\mathit{H}*\mathit{f}={\mathbf{F}}^{-1}\mathit{\upsilon },$$(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⁻¹H is the point spread function (PSF) and F⁻¹ν 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 X-ray imaging. Specifically, the reconstruction approach introduced in this paper is a regularization method with an l₁ 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²(ℝ²),$${\psi }_{a,\mathbf{t}}\left(\mathbf{X}\right)=a^{-1/2}\psi \left(\frac{\mathbf{x}-\mathbf{t}}{a}\right),$$(3)

where x = (x₁, x₂) is a point in ℝ², and t = (t₁, t₂) ∈ ℝ² 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²(ℝ²). The wavelet transform W(f) of a function f ∈ L²(ℝ²) is defined by Mallat (1999) $$\begin{array}{c}\mathcal{W}(f)(a,t)=\langle f,{\psi }_{a,t}\rangle \\ =\langle \widehat{f},{\widehat{\psi }}_{a,t}\rangle \\ ={\displaystyle {\int }_{{ℝ}^2}\widehat{f}(u,v)}\overline{{\widehat{\psi }}_{a,t}(u,v)}\hspace{0.17em}\text{d}u\text{d}v\end{array}$$(4)

where 〈f, ψ_a,t〉 is the scalar product in L²(ℝ²), and $\hat{f}$ and ${\hat{\psi }}_{a,\mathbf{t}}$ are the Fourier transform of f and ψ_a,t, respectively. If the admissibility condition is satisfied,$$C={\int }_{{ℝ}^2}\frac{{\left|\hat{\psi }(u,\upsilon )\right|}^2}{u^2+{\upsilon }^2}\mathrm{d}u\mathrm{d}\upsilon <\infty ,$$(5)

where $\hat{\psi }$ is the Fourier transform of the mother wavelet ψ, then it is possible to define the inverse wavelet transform as$$f\left(\mathbf{x}\right)=\frac{1}{C}{\int }_{{ℝ}^2}{\int }_0^{\infty }\mathcal{W}\left(f\right)(a,\mathbf{t}){\psi }_{a,\mathbf{t}}\left(\mathbf{x}\right)\mathrm{d}\mathbf{t}\frac{\mathrm{d}a}{a^3}\mathrm{.}$$(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$$\hat{\psi }(u,\upsilon )={\hat{\psi }}_M\left(\sqrt{u^2+{\upsilon }^2}\right)\mathrm{.}$$(7)

Consequently, from Eq. (3) we obtained$${\hat{\psi }}_{a,\mathbf{t}}(u,\upsilon )=a^{1/2}{\hat{\psi }}_M\left(a\sqrt{u^2+{\upsilon }^2}\right){\mathrm{e}}^{-2\pi \mathrm{i}(u,\upsilon )\cdot \mathbf{t}}\mathrm{.}$$(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$$\hat{\varphi }(u,\upsilon )={\hat{\varphi }}_M\left(\sqrt{u^2+{\upsilon }^2}\right),$$(9)

where ϕ_M(x) is the 1D Meyer scaling function. The FIWT is finally defined from Eq. (4),$$\mathcal{W}\left(f\right)(a,\mathbf{t})=a^{1/2}{F}^{-1}\left(\hat{\mathcal{F}}(u,\upsilon ){\hat{\psi }}_M\left(a\sqrt{u^2+{\upsilon }^2}\right)\right)\left(\mathbf{t}\right),$$(10)

where ℱ⁻¹ 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$$\begin{array}{cc}{\left(x_1\right)}_p=c_1-\frac{FOV}{2}+p\delta & p=0,\dots ,N-1,\end{array}$$(11) $$\begin{array}{cc}{\left(x_2\right)}_q=c_2-\frac{FOV}{2}+q\delta & q=0,\dots ,N-1,\end{array}$$(12)

where (c₁, c₂) 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₁, … , f_M) from f((x₁)_p, (x₂)_q) p, q = 0, … , N − 1, where M = N² 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₀ scales obtained according to the discretization$$\begin{array}{cc}{a}_j=2^{-j},& j=0,\dots ,j_0-1.\end{array}$$(13)

On the other hand, the discretization of the translation parameter t is made according to$$\begin{array}{cc}{\begin{array}{c}\left(t_1\right)\end{array}}_n=n\delta ,& n=0,\dots ,N-1,\end{array}$$(14) $$\begin{array}{cc}{\begin{array}{c}\left(t_2\right)\end{array}}_l=l\delta ,& l=0,\dots ,N-1.\end{array}$$(15)

Then, for each scale a_j and for each translation ((t₁)_n, (t₂)_l), we construct the M-dimensional vector ψ_j,n,l , which corresponds to the pixelization and reordering of the wavelet Eq. (3) for a = a_j, t₁ = (t₁)_n, and t₂ = (t₂)_l. Analogously, for each translation, we construct the M-dimensional 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₀ + 1) M-dimensional vectors {u_i, i = 1, … , M · (j₀ + 1)} = { ψ_j,n,l j = 0, … , j₀ − 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$${\Vert \mathit{f}\Vert }^2=\underset{i=1}{\overset{M\cdot (j_0+1)}{\mathrm{\Sigma }}}{\left|(\mathit{f},{\mathit{u}}_i)\right|}^2,$$(16)

and$$\mathit{f}=\underset{i=1}{\overset{M\cdot (j_0+1)}{\mathrm{\Sigma }}}(\mathit{f},{\mathit{u}}_i){\mathit{u}}_i,$$(17)

where (·, ·) denotes the canonical inner product in ℝ^M. Finally, W is the M · (j₀ + 1) × M matrix whose rows are made of the vectors u_i for i = 1, … , M · (j₀ + 1). The fact that ${\left\{{\mathit{u}}_i\right\}}_{i=1}^{M\cdot (j_0+1)}$ is a Parseval frame implies that W^TW = I and that the forward and inverse discretized FIWT can be written as$$\begin{array}{cc}\mathit{w}=\mathbf{W}\mathit{f},& \mathit{f}={\mathbf{W}}^{\mathrm{T}}\mathit{w},\end{array}$$(18)

respectively, where w is the column vector of dimension M · (j₀ + 1), whose components are$$\begin{array}{cc}{\mathit{w}}_i=(\mathit{f},{\mathit{u}}_i)& i=1,\dots ,M\cdot (j_0+1)\mathrm{.}\end{array}$$(19)

We conclude by observing that the computational complexity of the FIWT is O(N² log N).

4. Image reconstruction method

The reconstruction of f from ν is an ill-posed 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 X-ray imaging is to regularize the inverse problem with the ℓ₁ 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$$\underset{f}{\min }\{{\Vert \mathit{H}\cdot \mathbf{F}\mathit{f}-\mathit{\upsilon }\Vert }_2^2+\lambda {\Vert \mathbf{W}\mathit{f}\Vert }_1\}\mathrm{.}$$(20)

The data term of the objective function to minimize, ${\Vert \mathit{H}\cdot \mathbf{F}\mathit{f}-\mathit{\upsilon }\Vert }_2^2$, quantifies the prediction error with respect to the measurements. The regularization term, ||W f||₁, 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 shrinkage-thresholding algorithm (FISTA; Beck & Teboulle 2009), which is an algorithm for the optimization in ℓ₁ that is widely used in many fields, mainly because of its reliability and rapid convergence. The scheme is illustrated in Table 1 above, where$$\begin{array}{cc}{\begin{array}{c}\left({\mathcal{T}}_{\lambda /2}\left(\mathbf{x}\right)\right)\end{array}}_i={(\left|x_i\right|-\lambda /2)}_{+}\mathrm{sgn}\left(x_i\right)& i=1,\dots ,M(j_0+1)\end{array}$$(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 component-wise) 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 cross-validation (GCV) method (Golub et al. 1979), and the L-curve 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$$\begin{array}{cc}{\Vert \mathit{H}\cdot \mathbf{F}\mathit{f}-\mathit{\upsilon }\Vert }_2^2\le \epsilon ,& {\Vert \mathbf{W}\mathit{f}\Vert }_1\le E\end{array}$$(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₁ of the FISTA algorithm with λ = 0. The regularization parameter is then set equal to$$\begin{array}{c}\lambda =\frac{{\Vert \mathit{H}\cdot \mathbf{F}\mathit{f}-\mathit{\upsilon }\Vert }_2^2}{{\Vert \mathbf{W}{\mathit{f}}_1\Vert }_1}\mathrm{.}\end{array}$$(23)

We compared the performance of the Miller method with that of the L-curve 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, the value λ = λ_L provided by the L-curve method, the value λ = λ_M provided by the Miller method, and the value λ = λ₁ = 0.1 realizing over-regularization. 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 (non-regularized solution), panel b: λ = 3.1 × 10−3 as provided by the L-curve method, panel c: λ = 5.7 × 10−3 as provided by the Miller method, and panel d: λ = 0.1 (over-regularized 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.

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¹. 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⁴ photons cm⁻², 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⁵ photons cm⁻²) 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 X-ray 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 over-resolves 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 visibility-based 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 three-scale decomposition for the wavelet-based 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 not-conjugated, 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 visibility-based 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 visibility-based 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 second-generation 4-core Intel i7-2600 (3.40 GHz) CPU using IDL 8.4 on Ubuntu 16.04. The computational time complexity of VIS_WV is O(KN² 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₀ for FIWT in the time complexity analysis, since j₀ 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₀ = 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₀ = 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 X-ray 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 sparsity-enhancing penalty term in regularization. Here we introduced a wavelet-based deconvolution method promoting sparsity for hard X-ray image reconstruction from visibilities. The main aspects of this method are listed below.

• It relies on a isotropic wavelet transform because X-ray sources are either isotropic or characterized by a slow change of shape.

• It avoids the use of a computationally demanding catalog-based 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 X-ray sources are not properly isotropic, and we might therefore consider other types of multi-scale decomposition, like shearlets, that account for anisotropies in the image (Duval-Poo et al. 2015). However, this approach would lead to a computationally more demanding algorithm, and since the X-ray 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², 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 B2-C1”.

References

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 low-frequency 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 (non-regularized solution), panel b: λ = 3.1 × 10−3 as provided by the L-curve method, panel c: λ = 5.7 × 10−3 as provided by the Miller method, and panel d: λ = 0.1 (over-regularized 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 visibility-based 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