c ESO 2017

Astronomy & Astrophysics manuscript no. ms-2017-08-02
August 15, 2017

Solar hard X-ray imaging by means of
Compressed Sensing and

arXiv:1708.03877v1 [astro-ph.SR] 13 Aug 2017

Finite Isotropic Wavelet Transform
M. A. Duval-Poo1 , M. Piana1, 2 , and A. M. Massone2
1

Dipartimento di Matematica, Università degli Studi di Genova, Via Dodecaneso 35, 16146
Genova, Italy
e-mail: duvalpoo@dima.unige.it, e-mail: piana@dima.unige.it

2

CNR-SPIN, Via Dodecaneso 33, 16146 Genova, Italy
e-mail: annamaria.massone@cnr.it

August 15, 2017

ABSTRACT

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 X-ray 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 sparsity-promoting
regularized objective function by means of the Fast Iterative Shrinkage-Thresholding 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 X-rays (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 visibility-based 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.
Key words. Sun: flares – Sun: X-rays, gamma rays – Techniques: image processing

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,
Article number, page 1 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

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, picked up
at specific points of the Fourier plane, named (u, v) 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
2018. 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 provides as well a sampling of the Fourier
transform of the photon flux in the (u, v) 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), 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 X-ray image configuration is made of few, simple shapes (mainly, foot-points and loops)
and therefore it is straightforward to represent it as the superposition of a small number of basis
functions. It is well-known that an advantageous approach to realize compressed sensing is to use
a wavelet transform since wavelets can provide, in addition to compression, a multi-scale signal
representation.
Most wavelet implementations are associated to multi-resolution 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 radio-interferometry (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 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. Therefore in the present paper we built a 2D isotropic wavelet transform that
Article number, page 2 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

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 ill-posedness of the deconvolution problem, in this paper we formulated a sparsity-enhancing regularized version of the FIWT multi-scale sparse decomposition,
which we called the Finite Isotropic waVElet Compressed Sensing method (5-CS) and we used it
to obtain reconstructions of hard X-ray 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 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 dictionary made of few shapes, replicated many times for different scales and positions, and
then to realize sparsity with respect to this dictionary. However, 5-CS 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, catalogue-based 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 X-ray 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 (u, v) 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 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
Article number, page 3 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

(a) RHESSI

(b) STIX

Fig. 1: Sampling of the (u, v) plane performed by RHESSI (left panel) and STIX (right panel). For
better visualization, panel (a) shows the sampling starting from detector 3.

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 (u, v) 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
H · Ff = v,

(1)

where the original N × N photon flux image to reconstruct is lexico-graphically re-ordered to define
the vector f ∈ R M , with M = N 2 . Further, v ∈ C M is a sparse vector whose non-zero components
correspond to the measured visibilities, F ∈ C M×M is the discretized Fourier transform, H is a
sparse binary matrix that realizes the sampling in the (u, v) plane, and · denotes the entry-wise
product. If we apply the discretized inverse Fourier transform F−1 on both sides of (1) we obtain
F−1 H ∗ f = F−1 v,

(2)

where ∗ is the convolution operator. Therefore, the reconstruction of the flux image f from the given
visibilities v can be essentially viewed as a deconvolution problem, where F−1 H is the point spread
function and F−1 v is the dirty map from which the PSF blurring effect must be subtracted.
Article number, page 4 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

3. The Finite Isotropic Wavelet Transform

Let ψa,t be a family of functions defined by the translation and dilation of a mother wavelet function
ψ(x) ∈ L2 (R2 ), i.e.
ψa,t (x) = a−1/2 ψ

x − t

(3)

a

where x = (x1 , x2 ) is a point in R2 , t = (t1 , t2 ) ∈ R2 and a ∈ R+ are the translation and dilation
parameters, respectively. The normalization factor a−1/2 ensures that kψa,t k = kψk, where k · k is the
norm in L2 (R2 ). The wavelet transform W( f ) of a function f ∈ L2 (R2 ) is defined by (Mallat 1999)
W( f )(a, t) = h f, ψa,t i
= h fˆ, ψ̂a,t i
Z
fˆ(u, v)ψ̂a,t(u, v) dudv
=

(4)

R2

where h f, ψa,t i is the scalar product in L2 (R2 ) and fˆ and ψ̂ are the Fourier transform of f and ψ,
respectively. If the admissibility condition is satisfied, i.e.
C=

Z

R2

|ψ̂(u, v)|2
dudv < ∞,
u 2 + v2

(5)

then it is possible to define the inverse wavelet transform as
f (x) =

1
C

Z Z
R2

∞
0

W( f )(a, t)ψa,t(x) dt

da
.
a3

(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 ψ M (x) (Mallat 1999) and constructed
the 2D isotropic mother wavelet function as
√
ψ̂(u, v) = ψ̂ M ( u2 + v2 ) .

(7)

Consequently, from (3) we obtained

 √
ψ̂a,t (u, v) = a1/2 ψ̂ M a u2 + v2 e−2πi(u,v)·t .

(8)

Further, in order to span the whole (u, v) plane, we included in the wavelet framework the scaling
function φ and constructed it again in the Fourier domain as
√
φ̂(u, v) = φ̂ M ( u2 + v2 ) ,

(9)
Article number, page 5 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

where φ M (x) is the 1D Meyer scaling function. The FIWT is finally defined from (4), i.e.
 √


W( f )(a, t) = a1/2 F −1 fˆ(u, v)ψ̂ M a u2 + v2 (t) ,

(10)

where F −1 is the inverse Fourier transform.
In the imaging framework the flux distribution f (x) is pixelized at positions
(x1 ) p = c1 −

FOV
+ pδ p = 0, . . . , N − 1 ,
2

(11)

(x2 )q = c2 −

FOV
+ qδ q = 0, . . . , N − 1 ,
2

(12)

where (c1 , c2 ) is the image center, FOV is the image field-of-view and δ is the pixel dimension
in arcsec. With the same lexico-graphical re-ordering used in section 2 we obtain the vector f =
( f1 , . . . , f M ) from f ((x1 ) p , (x2 )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. Accordingly to
what done is Section 2, the vector f is again the unknown of the image reconstruction problem.
In the wavelet transformation we will consider j0 scales obtained according to the discretization
a j = 2− j ,

j = 0, . . . , j0 − 1 .

(13)

On the other hand, the discretization of the traslation parameter t is made according to
(t1 )n = nδ , n = 0, . . . , N − 1 ,

(14)

(t2 )l = lδ , l = 0, . . . , N − 1 .

(15)

Then, for each scale a j and for each translation ((t1 )n , (t2 )l ), we construct the M−dimensional vector
ψ j,n,l which corresponds to the pixelization and re-ordering of the wavelet (3) for a = a j , t1 =
(t1 )n , and t2 = (t2 )l . Analogously, for each translation we construct the M−dimensional vector φn,l ,
which is the result of the pixelization and re-ordering of the scaling function φ. This leads to the
construction of the set of M · ( j0 + 1) M-dimensional vectors {ui , i = 1, . . . , M · ( j0 + 1)} =
{ψ j,n,l j = 0, . . . , j0 − 1 ; n, l = 0, . . . , N − 1}∪ {φn,l

n, l = 0, . . . , N − 1}. This set provides a

M

Parseval frame, i.e., given f ∈ R , then
kfk2 =

M·(
j0 +1)
X
i=1

|(f, ui)|2 ,

(16)

and
f=

M·(
j0 +1)
X

(f, ui )ui ,

(17)

i=1

Article number, page 6 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

φ̂(u, v)

ψ̂a=1 (u, v)

ψ̂a=2−1 (u, v)

ψ̂a=2−2 (u, v)

Fig. 2: 2D isotropic extensions of the Meyer scaling function and mother wavelet with no translation. The functions are represented in the Fourier domain for 3 scales.
where (·, ·) denotes the canonical inner product in R M . Finally, W is the M · ( j0 + 1) × M matrix

M·( j0 +1)
whose rows are made of the vectors ui for i = 1, . . . , M · ( j0 + 1). The fact that {ui }i=1
is a

Parseval frame implies that WT W = I and that the forward and inverse discretized FIWT can be
written as
w = Wf , f = WT w ,

(18)

respectively, where w is the column vector of dimension M · ( j0 + 1) whose components are
wi = (f, ui ) i = 1, . . . , M · ( j0 + 1) .

(19)

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

4. The image reconstruction method
The reconstruction of f from v 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 valid approach in hard X-ray imaging is to regularize the inverse problem
with the ℓ1 norm in some transformation domain. In fact, the method we propose in this paper,
which we named Finite Isotropic waVElet transform Compressed Sensing (5-CS), addresses the
optimization problem
o
n
min ||H · Ff − v||22 + λ||Wf||1 .

(20)

f

The data term of the objective function to minimize, ||H·Ff−v||22, quantifies the prediction error with
respect to the measurements. The regularization term, ||Wf||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 non-linear optimization. Here we used the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) (Beck & Teboulle
2009), which is an ℓ1 solver widely used in many fields mainly because of its reliability and rapid
Article number, page 7 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

(b) λL = 3.1 × 10−3

(a) λ0 = 0

(c) λ M = 5.7 × 10−3

(d) λ1 = 0.1

Fig. 3: July 23 2002 solar flare image reconstruction by 5-CS with different choices for λ: (a)
λ = 0 (non-regularized solution), (b) λ = 3.1 × 10−3 as provided by the L-curve method, (c)
λ = 5.7 × 10−3 as provided by the Miller method, (d) λ = 0.1 (over-regularized solution). The time
interval is between 00:29:10-00:30:19 UT while the energy channel is 36 − 41 keV. Detectors from
2 through 9 have been employed.

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 L-curve 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
||H · Ff − v||22 ≤ ε,

||Wf||1 ≤ E

(21)

are known, or can be estimated from the dirty map, then the regularization parameter can be chosen
as λ = ε/E. In order to satisfy conditions (21), the bounds {ε, E} are estimated performing the first
iteration f1 of the FISTA algorithm with λ = 0. The regularization parameter is then set equal to
λ=

||H · Ff1 − v||22
.
||Wf1 ||1

(22)

We compared the performance of the Miller method with the one 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 5-CS for four
values of the regularization parameter: λ = λ0 = 0, the value λ = λL provided by the L-curve
method, the value λ = λ M provided by the Miller method, and the value λ = λ1 = 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.
Article number, page 8 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

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 Software1 . Next, 5-CS 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 104 photons/cm2, but in the first image the two foot-points 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 (105 photons/cm2) and contain two loops with significantly different curvatures.
Figure 4 compares the images reconstructed by 5-CS 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 5-CS and CLEAN recover the physical parameters
with comparable accuracy, although 5-CS is able to significantly reduce the impact of spurious
artifacts.

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 5-CS with the ones obtained by two standard visibility-based 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 foot-points
(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 3-scale decomposition for the wavelet-based 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 5-CS 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).
1

https://stix.cs.technik.fhnw.ch/
Article number, page 9 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

(a) Simulated Double Footpoint Flare (SDFF1)

(b) Simulated Double Footpoint Flare (SDFF2)

(c) Simulated Loop Flare (SLF1)

(d) Simulated Loop Flare (SLF2)

Fig. 4: Reconstructed images for four different STIX simulated flare events. Left column: the original source shape. Middle column: reconstructions provided by CLEAN using visibilities. Right
column: reconstructions provided by 5-CS.

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
5-CS is near 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 j0
for FIWT on the time complexity analysis, since j0 is usually a small number.
Article number, page 10 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

Table 1: Physical parameters inferred from the reconstructed maps in (Figure 4) and compared to
the ones in the simulated maps (position and FHWM are measured in arcsec). The reconstructions
have been made 10 times, in correspondence with 10 realizations of the STIX simulated data. The
Table indicates the average values of the parameters and the corresponding standard deviations.
Simulated Map

CLEAN

5-CS

Peak 1

X
FWHM

15.0

Peak 2

Simulated Double Footpoint Flare (SDFF1)
20.0

Y

-20.0

X

-5.0

Y

5.0

FWHM

10.0

Flux Ratio

2.00

20.4 ±0.9

19.2 ±1.5

-18.9 ±0.9

-21.8 ±1.1

-4.0 ±0.7

-5.0 ±1.2

12.8 ±0.8

11.9 ±1.6

9.6 ±3.4
4.8 ±0.6

9.0 ±4.4
3.5 ±2.9

1.41 ±0.10

1.78 ±0.42

-28.9 ±0.6

-29.3 ±0.8

13.2 ±0.9

9.6 ±0.8

Peak 1

X
FWHM

10.0

Peak 2

Simulated Double Footpoint Flare (SDFF2)
-30.0

Y

0.0

X

30.0

Y

0.0

FWHM

10.0

Flux Ratio

1.00

1.0 ±0.7

29.5 ±0.5
-0.1 ±0.6

12.9 ±0.6

Peak

10.0

Y

-15.0

Peak

10.0

Y

0.0

-1.4 ±0.5
9.3 ±0.6

1.02 ±0.09

9.0 ±0.5

9.6 ±0.5

-13.6 ±0.5

-14.6 ±0.7

9.8 ±0.4

8.7 ±0.7

Simulated Loop Flare (SLF2)
X

28.0 ±0.8

0.99 ±0.10

Simulated Loop Flare (SLF1)
X

0.4 ±0.7

-0.1 ±0.3

-0.4 ±0.7

Figure 8 shows the computational performance of 5-CS 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 K for the reconstruction of an image
map of size N = 128 and considering j0 = 3. On the other hand, Figure 8 (b) shows how the
computational time increases when 5-CS recovers images with higher and higher size, where in this
case, the number of iterations and considered scales are fixed at K = 30 and j0 = 3, respectively.
Since most reconstructions are performed with an image size of 128 pixels and our method usually
solves the reconstruction problem with less than 100 iterations, we can conclude that the average
time required by 5-CS for reconstructing a single standard image is between 1.5 and 2.5 seconds.

6. Conclusions
Hard X-ray 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
Article number, page 11 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

(a) February 20, 2002 – 11:06:02-11:06:34 UT – (22-26 keV)

(b) April 15, 2002 – 00:05:00-00:10:00 UT – (12-14 keV)

(c) December 02, 2003 – 22:54:00-22:58:00 UT – (18-20 keV)

(d) August 31, 2004 – 05:33:00-05:38:00 UT – (10-12 keV)

(e) May 13, 2013 – 02:04:16-02:04:48 UT – (6-12 keV)

Fig. 5: Reconstructed images from different flaring events using three reconstruction methods:
CLEAN (left), uv_smooth (center), and 5-CS (Right). RHESSI detectors from 2 to 9 have been
used to generate the visibilities.

Article number, page 12 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing

Fig. 6: Images reconstructed by 5-CS 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: Images reconstructed by 5-CS 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.
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 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 that
– It relies on a continuous isotropic wavelet transform, coherently to the fact that X-ray sources
are either isotropic or characterized by a slow change of shape.
Article number, page 13 of 15

4

20

3

15
Time (seconds)

Time (seconds)

Duval Poo et al: Hard X-ray imaging and compressed sensing

2

10

5

1

0

0

20

40

60

80

100

64

128

Iterations

256
Map size (pixels)

(a)

(b)

Fig. 8: Computational performance of 5-CS. The computational time is compared versus (a) the
number of FISTA iterations and (b) the size of the reconstructed image. The test is made using the
synthetic STIX visibilities used in Figure 4(a).
– It avoids the use of a computationally demanding catalogue-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.
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 B2-C1"

References
Aschwanden, M. J. & Schmahl, E. 2002, Solar Physics, 210, 193
Beck, A. & Teboulle, M. 2009, SIAM journal on imaging sciences, 2, 183
Benz, A., Krucker, S., Hurford, G., et al. 2012, in SPIE Astronomical Telescopes+ Instrumentation, International Society
for Optics and Photonics, 84433L–84433L
Bertero, M. & Boccacci, P. 1998, Introduction to inverse problems in imaging (CRC press)
Bobin, J., Starck, J. L., & Ottensamer, R. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 718
Bong, S. C., Li, L. J., Gary, D. E., & Yun, H. S. 2006, Astrophysical Journal, 636, 1159
Candes, E. J. & Wakin, M. B. 2008, Signal processing magazine, 25, 21
Cornwell, T. J. & Evans, K. 1985, Astronomy and Astrophysics, 143, 77
Donoho, D. L. 2006, IEEE Transactions on Information Theory, 52, 1289
Emslie, A. G., Kontar, E. P., Krucker, S., & Lin, R. P. 2003, The Astrophysical Journal, 595, L107
Engl, H. W., Hanke, M., & Neubauer, A. 1996, Regularization of inverse problems, Vol. 375 (Springer Science & Business
Media)

Article number, page 14 of 15

Duval Poo et al: Hard X-ray imaging and compressed sensing
Garsden, H., Girard, J. M., L, S. J., & S, C. 2015, Astronomy & Astrophysics, 575, A90
Giordano, S., Pinamonti, N., Piana, M., & Massone, A. M. 2015, SIAM Journal on Imaging Sciences, 8, 1315
Golub, G. H., Heath, M., & Wahba, G. 1979, Technometrics, 21, 215
Hansen, P. C. 1992, SIAM review, 34, 561
Häuser, S. & Steidl, G. 2014, ArXiv
Högbom, J. 1974, Astronomy and Astrophysics Supplement Series, 15, 417
Hurford, G. J., Schmahl, E. J., Schwartz, R. A., et al. 2002, Solar Physics, 210, 61
Li, F., Cornwell, T. J., & de Hoog, F. 2011, Astronomy & Astrophysics, 528, A31
Lin, R., Dennis, B., Hurford, G., et al. 2002, Solar Physics, 210, 3
Mallat, S. 1999, A wavelet tour of signal processing (Acaic press)
Massone, A. M., Emslie, A. G., Hurford, G., et al. 2009, The Astrophysical Journal, 703, 2004
Miller, K. 1970, SIAM Journal on Mathematical Analysis, 1, 52
Morozov, V. A. 1984, Methods for solving incorrectly posed problems (Springer New York)
Oster, G., Wasserman, M., & Zwerling, C. 1964, JOSA, 54, 169
Portilla, J. & Simoncelli, E. P. 2000, International journal of computer vision, 40, 49
Starck, J.-L., Fadili, J., & Murtagh, F. 2007, IEEE Transactions on Image Processing, 16, 297

Article number, page 15 of 15

