© ESO, 2012

1. Introduction

1.1. Literature overview

The gamma-ray sky has been studied with unprecedented sensitivity and image capability thanks to the Large Area Telescope (LAT), the main instrument of the Fermi Gamma-ray Space Telescope (Atwood et al. 2009), in an energy range between 20 MeV to greater than 300 GeV. The detection of gamma-ray point sources is difficult for two main reasons: the Poisson noise and the instrument’s point spread function (PSF). The Poisson noise is due to the fluctuations in the number of detected photons. Moreover, the effect of Poisson noise is strong because of the weakness of the fluxes of celestial gamma rays, especially outside the Galactic plane and far away from intense sources. The PSF width is strongly energy-dependent, varying from about 3.5 at 100 MeV to less than (68% containment) at 10 GeV. Owing to large-angle multiple scattering in the tracker, the PSF has broad tails, for which the 95%/68% containment ratio may be as large as 3.

An extensive literature exists on Poisson noise removal and the interested reader may refer to Schmitt et al. (2011) and Starck et al. (2010) for a thorough review. Motivated by new X-ray and gamma-ray data challenges, several restoration methods have been released in astrophysics that are based on wavelets (Movit 2009; Starck et al. 2009; Faÿ et al. 2011; Schmitt et al. 2010) or the Bayesian machinery (Conrad et al. 2007; Norris et al. 2010). Wavelets have also been used for source detection in Fermi data (Abdo et al. 2010), and a first Poisson denoising algorithm for spherical data was proposed in Schmitt et al. (2010). Starck et al. (2009) developed a denoising approach that effectively handles multichannel data acquired on a cartesian grid, and where the third dimension can be any physically meaningful index such as wavelength, energy, or time. While traditional techniques integrate over all the third dimension in order to improve the signal-to-noise ratio (S/N) of the sources, their approach allows us to detect the sources while preserving all of their spectral information without sacrificing any the sensitivity. Nonetheless, none of these methods take into account the point spread function (PSF) of the instrument. Deconvolution can be very helpful in many situations such as source identification or flux estimation in the low energy bands where the resolution degrades severely. To the best of our knowledge, no general method for both multichannel denoising and deconvolving on the sphere has been developed in the literature.

1.2. Contributions

We propose a general framework for denoising and deconvolution that complies with all the above requirements, namely deals with

• 1.

  multichannel data;

• 2.

  acquired on the sphere;

• 3.

  and contaminated by Poisson noise.

Our approach builds upon the concept of variance stabilization applied to the spherical wavelet transform coefficients (Starck et al. 2006). This approach gives a multiscale representation of the Poisson data with variance-stabilized coefficients, which can be treated as if they were contaminated by a zero-mean white Gaussian noise. The developed algorithms are validated on simulated Fermi HEALPix multichannel cubes (nside = 256) with energy bands ranging from 50 MeV to 50 GeV.

1.3. Paper organization

The rest of the paper is organized as follows. In Sect. 2, we present a new wavelet transform for multichannel data on the sphere, and we show how the variance stabilization transform can be introduced into the decomposition. Section 3 details the way that this transform can be used for Gaussian and Poisson noise removal. Section 4 describes our deconvolution algorithm on the sphere for both mono-channel and multichannel spherical data. In Sect. 5, we finally draw some conclusions and give possible perspectives of this work.

1.4. Notations

For a real discrete-time filter whose impulse response is h [i] , is its time-reversed version. The discrete circular convolution product of two signals is written  ⋆ , where the term circular represents periodic boundary conditions. The symbol δ [i]  is the Kronecker delta.

For the wavelet representation, the low-pass analysis filter is denoted h and the high-pass is taken as g = δ − h throughout the paper. We denote the up-sampled version of h as h ^↑ j [l]  = h [l]  if l/2^j ∈ Z and 0 otherwise. We define for and h⁽⁰⁾ = δ.

The scaling and wavelet functions used for the analysis (respectively, synthesis) are denoted φ (with $\mathit{\phi }\mathrm{\left(}\frac{\mathit{x}}{\mathrm{2}}\mathrm{\right)}\mathrm{=}{\sum }_{\mathit{k}}\mathit{h}\mathrm{\left[}\mathit{k}\mathrm{\right]}\mathit{\phi }\mathrm{(}\mathit{x}\mathrm{-}\mathit{k}\mathrm{)}\mathit{,x}\mathrm{\in }R \mathrm{and} \mathit{k}\mathrm{\in }Z$) and ψ (with $\mathit{\psi }\mathrm{\left(}\frac{\mathit{x}}{\mathrm{2}}\mathrm{\right)}\mathrm{=}{\sum }_{\mathit{k}}\mathit{g}\mathrm{\left[}\mathit{k}\mathrm{\right]}\mathit{\phi }\mathrm{(}\mathit{x}\mathrm{-}\mathit{k}\mathrm{)}\mathit{,x}\mathrm{\in }R$and   k ∈ Z) (respectively, $\begin{array}{}{\mathrm{􏽥}}^{}\\ \mathit{\phi }\end{array}$ and $\begin{array}{}{\mathrm{􏽥}}_{}\\ \mathit{\psi }\end{array}$).

2. Multiscale representation for multichannel spherical data with poisson noise

2.1. Fast undecimated 2D-1D wavelet decomposition/reconstruction on the sphere

Our goal is to analyze multichannel data acquired on a sphere with a non-isotropic two-and-one-dimensional (2D-1D) wavelet, where the two first dimensions are spatial (longitude and latitude) and the third dimension is either the time or the energy. Since the dimensions do not have the same physical meaning, it appears natural that the wavelet scale along the third dimension (energy or time) should not be connected to the spatial scale. Hence, we define the wavelet function as $\mathit{\psi }\mathrm{\left(}{\mathit{k}}_{\mathit{\theta }}\mathit{,}{\mathit{k}}_{\mathit{\varphi }}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right)}\mathrm{=}{\mathit{\psi }}^{\mathrm{\left(}\mathit{\theta \phi }\mathrm{\right)}}\mathrm{\left(}{\mathit{k}}_{\mathit{\theta }}\mathit{,}{\mathit{k}}_{\mathit{\varphi }}\mathrm{\right)}{\mathit{\psi }}^{\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{\left(}{\mathit{k}}_{\mathit{t}}\mathrm{\right)}\mathit{,}$(1)where ψ^(θφ) is the spherical two dimensional (2D) spatial wavelet and ψ^(t) the one-dimensional (1D) wavelet along the third dimension. Similarly to Starck et al. (2009), we consider only isotropic and dyadic spatial scales. We build the discrete 2D-1D wavelet decomposition by first taking a spherical 2D undecimated wavelet transform for each k_t, followed by a 1D wavelet transform for each spatial wavelet coefficient along the third dimension.

Hence, for a given multichannel data set on the sphere Y [k_θ,k_ϕ,k_t]  and after applying first the 2D spherical undecimated wavelet transform, we have the reconstruction formula $\begin{array}{ccc}\mathit{Y}\mathrm{\left[}{\mathit{k}}_{\mathit{\theta }}\mathit{,}{\mathit{k}}_{\mathit{\varphi }}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}& & \mathrm{=}{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}}\mathrm{\left[}{\mathit{k}}_{\mathit{\theta }}\mathit{,}{\mathit{k}}_{\mathit{\varphi }}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\\ & & \mathrm{+}\sum _{{\mathit{j}}_{\mathrm{1}}\mathrm{=}\mathrm{1}}^{\mathit{J}\mathrm{1}}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}}\mathrm{\left[}{\mathit{k}}_{\mathit{\theta }}\mathit{,}{\mathit{k}}_{\mathit{\varphi }}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathit{,} \mathrm{\forall }{\mathit{k}}_{\mathit{t}}\mathit{,}\end{array}$(2)where J₁ is the number of spatial scales, a_J1 is the (spatial) approximation subband, and $\mathrm{\{}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}}{\mathrm{\}}}_{{\mathit{j}}_{\mathrm{1}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{1}}}$ are the (spatial) detail subbands. To simplify the notations in the sequel, we replace the two spatial indices by a single index k_r, which corresponds to the pixel index. Equation (2)now reads $\mathit{Y}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{=}{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{+}\sum _{{\mathit{j}}_{\mathrm{1}}\mathrm{=}\mathrm{1}}^{\mathit{J}\mathrm{1}}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathit{,} \mathrm{\forall }{\mathit{k}}_{\mathit{t}}\mathit{.}$(3)For each spatial location k_r and each 2D wavelet scale j₁, we then apply a 1D wavelet transform along t on the spatial wavelet coefficients w_j1 [k_r,·]  such that ${\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{=}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{+}\sum _{{\mathit{j}}_{\mathrm{2}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{2}}}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathit{,} \mathrm{\forall }\mathrm{\left(}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right)}\mathit{,}$(4)where J₂ is the number of scales along t. The approximation spatial subband a_J1 is processed in a similar way, hence yielding ${\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{=}{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{+}\sum _{{\mathit{j}}_{\mathrm{2}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{2}}}{\mathit{w}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathit{,} \mathrm{\forall }\mathrm{\left(}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right)}\mathit{.}$(5)Inserting Eqs. (4)and (5)into (3), we obtain the 2D-1D spherical undecimated wavelet representation of Y: $\begin{array}{ccc}\mathit{Y}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}& & \mathrm{=}{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{+}\sum _{{\mathit{j}}_{\mathrm{1}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{1}}}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\\ & & \mathrm{+}\sum _{{\mathit{j}}_{\mathrm{2}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{2}}}{\mathit{w}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{+}\sum _{{\mathit{j}}_{\mathrm{1}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{1}}}\sum _{{\mathit{j}}_{\mathrm{2}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{2}}}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathit{.}\end{array}$(6)In this expression, four kinds of coefficients can be distinguished:

• Detail-detail coefficients ( and ): $\begin{array}{ccc}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}& \mathrm{=}& \mathrm{(}\mathit{\delta }\mathrm{-}\mathit{h̅}\mathrm{1}\mathrm{D}\mathrm{)}\mathit{\star }\mathrm{(}{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{(}{\mathit{j}}_{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{\star }{\mathit{a}}_{{\mathit{j}}_{\mathrm{1}}\mathrm{-}\mathrm{1}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\\ & & \mathrm{-}{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{(}{\mathit{j}}_{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{\star }{\mathit{a}}_{{\mathit{j}}_{\mathrm{1}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathrm{)}\mathit{.}\end{array}$(7)

• Approximation-detail coefficients (j₁ = J₁ and ): ${\mathit{w}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathrm{=}{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{(}{\mathit{j}}_{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{\star }{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathrm{-}{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{\left(}{\mathit{j}}_{\mathrm{2}}\mathrm{\right)}}\mathit{\star }{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathit{.}$(8)

• Detail-approximation coefficients ( and j₂ = J₂): ${\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathrm{=}{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{\left(}{\mathit{J}}_{\mathrm{2}}\mathrm{\right)}}\mathit{\star }{\mathit{a}}_{{\mathit{j}}_{\mathrm{1}}\mathrm{-}\mathrm{1}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathrm{-}{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{\left(}{\mathit{J}}_{\mathrm{2}}\mathrm{\right)}}\mathit{\star }{\mathit{a}}_{{\mathit{j}}_{\mathrm{1}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathit{.}$(9)

• Approximation-approximation coefficients (j₁ = J₁ and j₂ = J₂): ${\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathrm{=}{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{\left(}{\mathit{J}}_{\mathrm{2}}\mathrm{\right)}}\mathit{\star }{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathit{.}$(10)

2.2. Multi-scale variance stabilizing transform on the sphere (MS-VSTS)

Schmitt et al. (2010) proposed a multiscale variance stabilizing transform adapted for Poisson spherical data. This transform was dubbed the multi-scale variance stabilizing transform on the Sphere (MS-VSTS). The MS-VSTS is a multi-scale decomposition method designed for Poisson noise that is based on a variance stabilizing transform(VST). Poisson noise is indeed signal-dependent, which complicates its removal. The aim of a VST is to get rid of this signal-dependence by transforming a Poisson distribution into a Gaussian distribution of known variance. In a nutshell, the MS-VSTS consists of plugging a VST into a multi-scale transform – the isotropic undecimated wavelet transform on the sphere (IUWTS) – in order to realize (approximately) Gaussian zero-mean multiscale coefficients with constant variance. The noise on coefficients can then be easily removed using Gaussian denoising methods such as wavelet shrinkage, as we see in the next section.

The MS-VSTS scheme is defined recursively by inserting a (nonlinear) square-root VST into the IUWTS steps, that is $\begin{array}{ccc}\mathrm{IUWTS}& & \{\begin{array}{c}\\ & & \\ & & \end{array}\mathrm{=}\mathrm{\Rightarrow }\\ & & \begin{array}{c}\\ \\ \end{array}\{\begin{array}{c}\\ & & \\ & & \end{array}\end{array}$(11)where T_j is the VST operator on scale j${\mathit{T}}_{\mathit{j}}\mathrm{\left(}{\mathit{a}}_{\mathit{j}}\mathrm{\right)}\mathrm{=}{\mathit{b}}^{\mathrm{\left(}\mathit{j}\mathrm{\right)}}\mathrm{sign}\mathrm{(}{\mathit{a}}_{\mathit{j}}\mathrm{+}{\mathit{c}}^{\mathrm{\left(}\mathit{j}\mathrm{\right)}}\mathrm{)}\sqrt{\mathrm{|}{\mathit{a}}_{\mathit{j}}\mathrm{+}{\mathit{c}}^{\mathrm{\left(}\mathit{j}\mathrm{\right)}}\mathrm{|}}\mathit{,}$(12)with the VST constants b^(j) and c^(j) that depend solely on the filter h and the scale level j. Zhang et al. (2008) showed that the MS-VSTS detail coefficients d_j on locally homogeneous parts of the underlying intensity signal follow asymptotically a zero-mean normal distribution with an intensity-independent variance that relies only on the filter h and the current scale j. Consequently, for a given h, both the stabilized variances and the constants b^(j) and c^(j) can be pre-computed once and for all (Schmitt et al. 2010).

2.3. Multichannel MS-VSTS

We now extend the MS-VSTS machinery to the multichannel case. This amounts to plugging the VST into the spherical 2D-1D undecimated wavelet transform introduced in Sect. 2.1. This again gives rise to four types of coefficients that take the following forms:

• Detail-detail coefficients ( and ): $\begin{array}{ccc}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}& & \mathrm{=}\mathrm{(}\mathit{\delta }\mathrm{-}\mathit{h̅}\mathrm{1}\mathrm{D}\mathrm{)}\mathit{\star }({\mathit{T}}_{{\mathit{j}}_{\mathrm{1}}\mathrm{-}\mathrm{1}\mathit{,}{\mathit{j}}_{\mathrm{2}}\mathrm{-}\mathrm{1}}({\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{(}{\mathit{j}}_{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}}\\ & & \mathit{\star }{\mathit{a}}_{{\mathit{j}}_{\mathrm{1}}\mathrm{-}\mathrm{1}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]})\mathrm{-}{\mathit{T}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}\mathrm{-}\mathrm{1}}({\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{(}{\mathit{j}}_{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{\star }{\mathit{a}}_{{\mathit{j}}_{\mathrm{1}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}))\mathit{.}\end{array}$(13)

• Approximation-detail coefficients (j₁ = J₁ and ): $\begin{array}{ccc}{\mathit{w}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}& & \mathrm{=}{\mathit{T}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}\mathrm{-}\mathrm{1}}\mathrm{(}{{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{(}{\mathit{j}}_{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{\star }{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}}^{\mathrm{)}}\\ & & \mathrm{-}{\mathit{T}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{(}{{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{\left(}{\mathit{j}}_{\mathrm{2}}\mathrm{\right)}}\mathit{\star }{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}}^{\mathrm{)}}\mathit{.}\end{array}$(14)

• Detail-approximation coefficients ( and j₂ = J₂): $\begin{array}{ccc}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}& & \mathrm{=}{\mathit{T}}_{{\mathit{j}}_{\mathrm{1}}\mathrm{-}\mathrm{1}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{(}{{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{\left(}{\mathit{J}}_{\mathrm{2}}\mathrm{\right)}}\mathit{\star }{\mathit{a}}_{{\mathit{j}}_{\mathrm{1}}\mathrm{-}\mathrm{1}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}}^{\mathrm{)}}\\ & & \mathrm{-}{\mathit{T}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{(}{{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{\left(}{\mathit{J}}_{\mathrm{2}}\mathrm{\right)}}\mathit{\star }{\mathit{a}}_{{\mathit{j}}_{\mathrm{1}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}}^{\mathrm{)}}\mathit{.}\end{array}$(15)

• Approximation-approximation coefficients (j₁ = J₁ and j₂ = J₂): ${\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathrm{=}{\mathit{h}}_{\mathrm{1}\mathrm{D}}^{\mathrm{\left(}{\mathit{J}}_{\mathrm{2}}\mathrm{\right)}}\mathit{\star }{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}}\mathrm{[}{\mathit{k}}_{\mathit{r}}\mathit{,}\mathrm{·}\mathrm{]}\mathit{.}$(16)

In summary, all 2D-1D wavelet coefficients are now stabilized, and the noise in all these wavelet coefficients is a zero-mean Gaussian with known variance that depends solely on h on the resolution levels (j₁,j₂). As before, these variances can be easily tabulated.

3. Application to multichannel denoising

We define X to be the noiseless data and Y their observed noisy version. In the case of the additive zero-mean white Gaussian noise, we have , and for the Poisson noise . The main objective behind denoising is to estimate X from Y.

3.1. Warm-up: Gaussian noise

We start with the simple and instructive case where the noise in Y is additive white Gaussian. As the spherical 2D-1D undecimated wavelet transform described in Sect. 2.1 is linear, the noise remains Gaussian in the transform domain. Therefore, the thresholding strategies that have been developed for wavelet Gaussian denoising can be applied to the spherical 2D-1D wavelet transform. Denoting the thresholding operator as TH, the denoised 2D-1D estimate of X obtained by thresholding the wavelet coefficients in Eq. (6)reads $\begin{array}{ccc}\begin{array}{}\mathrm{􏽥}\\ \mathit{X}\end{array}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}& \mathrm{=}& {\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{+}\sum _{{\mathit{j}}_{\mathrm{1}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{1}}}\mathrm{TH}\mathrm{\left(}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{\right[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\left]}\mathrm{\right)}\\ & & \mathrm{+}\sum _{{\mathit{j}}_{\mathrm{2}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{2}}}\mathrm{TH}\mathrm{\left(}{\mathit{w}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{\right[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\left]}\mathrm{\right)}\\ & & \mathrm{+}\sum _{{\mathit{j}}_{\mathrm{1}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{1}}}\sum _{{\mathit{j}}_{\mathrm{2}}\mathrm{=}\mathrm{1}}^{{\mathit{J}}_{\mathrm{2}}}\mathrm{TH}\mathrm{\left(}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{\right[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\left]}\mathrm{\right)}\mathit{.}\end{array}$(17)A typical choice of TH is the hard thresholding operator parametrized by the scalar threshold , i.e. $\mathrm{TH}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathrm{\left|}\mathit{x}\mathrm{\right|}\mathit{<}\mathit{\tau ,}\\ & \mathrm{otherwise}\mathit{.}\end{array}$(18)The threshold τ is typically chosen to be between three and five times the noise standard deviation.

3.2. Poisson noise

The treatment of Poisson noise is far more complicated than its Gaussian counterpart, the main reason being that the noise variance is equal to its mean. This is where the MS-VSTS comes into play. The role of the MS-VSTS is indeed to get rid of this dependence of the variance on the mean by ensuring that the transformed coefficients are Gaussian with constant variance (without loss of generality, this variance can be assumed to be equal to 1). In other words, after the MS-VSTS, we are brought to a Gaussian denoising problem where standard thresholding approaches apply.

Nevertheless, denoising is not straightforward because there is no explicit reconstruction formula available because of the form of the non-linear stabilization equations above. Formally, the stabilizing operators T_j1,j2 and the convolution operators along the spatial and the third dimensions do not commute, even though the filter bank satisfies the exact reconstruction formula. To circumvent this difficulty, we propose to solve this reconstruction problem by advocating an iterative reconstruction scheme.

3.3. Iterative reconstruction

We define W to be the transform operator associated with the 2D-1D IUWTS described in Sect. 2.1, and R to be its inverse transform. We define ℳ to be the multiresolution support, which is determined by the set of significant coefficients detected among WY at each scale (j₁,j₂) and location (k_r,k_t), i.e. $\mathrm{ℳ}\mathrm{=}\mathrm{\left\{}\mathrm{\right(}{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\left)}|\mathrm{\right(}\mathbf{W}\mathit{Y}{\mathrm{)}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{issignificant}\mathrm{\}}\mathit{.}$(19)We define M to be the orthogonal projector onto ℳ, i.e. ∀d$\mathrm{\left(}\mathbf{M}\mathit{d}{\mathrm{)}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{\right[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{]}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathrm{\left(}{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right)}\mathrm{\in }\mathrm{ℳ}\mathit{,}\\ & \mathrm{otherwise}\mathrm{.}\end{array}$(20)Our goal is to seek a solution $\begin{array}{}{\mathrm{􏽥}}_{}\\ \mathit{X}\end{array}$ that preserves the significant structures of the original data by reproducing exactly the same wavelet coefficients as those of the input data Y, but only on scales and at positions where significant coefficients have been detected. Furthermore, as Poisson intensity functions are positive by nature, a positivity constraint is imposed on the solution. It is clear that there are many solutions satisfying the positivity and multiresolution support consistency requirements, e.g. Y itself. Thus, our reconstruction problem based solely on these constraints is an ill-posed inverse problem that must be regularized. Typically, the solution in which we are interested must be sparse by involving the lowest budget of wavelet coefficients. Therefore, our reconstruction is formulated as a constrained sparsity-promoting minimization problem over the transform coefficients d$\begin{array}{ccc}& & \underset{\mathit{d}}{\min }\mathrm{\parallel }\mathit{d}{\mathrm{\parallel }}_{\mathrm{1}}\mathrm{subjectto}\\ & & \{\begin{array}{c}\\ \\ \end{array}\end{array}$(21)and the intensity estimate $\begin{array}{}{\mathrm{􏽥}}_{}\\ \mathit{X}\end{array}$ is reconstructed as $\begin{array}{}{}_{\mathrm{􏽥}}\\ \mathit{X}\end{array}\mathrm{=}\mathbf{R}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{d}\end{array}$, where $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{d}\end{array}$ is a global minimizer of Eq. (21). We recall that  ∥·∥ ₁ =  ∑ _j1,j2,kr,kt | d_j1,j2 [k_r,k_t]  |  is the ℓ₁-norm playing the role of regularization and is well-known to promote sparsity (Donoho 2004). This problem can be solved efficiently using the hybrid steepest descent algorithm (Yamada 2001; Zhang et al. 2008), and requires about ten iterations in practice. Transposed into our context, its main steps can be summarized as follows:

where P₊ is the orthogonal projector onto the positive orthant, and ST_βn is the entry-wise soft-thresholding operator with threshold β_n, i.e. for x ∈ R, ST_βn(x) = max(0,1 − β_n/|x|)x.

The final multichannel MS-VSTS Poisson noise removal algorithm is summarized in the following steps:

Fig. 1 Result of the multichannel Poisson denoising algorithm on simulated Fermi data over the energy band 220–360 MeV. Top: simulated intensity skymap. Middle: simulated noisy skymap. Bottom: denoised skymap. Maps are on a logarithmic scale.

Fig. 2 Result of the multichannel Poisson denoising algorithm on simulated Fermi data over the energy band 589–965 MeV. Top: simulated intensity skymap. Middle: simulated noisy skymap. Bottom: denoised skymap. Maps are on a logarithmic scale.

3.4. Experiments

The multichannel MS-VSTS algorithm has been applied to a simulated Fermi data set, with 14 energy bands between 50 MeV and 1.58 GeV. Figures 1 and 2 depict the denoising results for two energy bands. The algorithm is able to recover most of the sources, even the faint ones, on each energy band. Even more importantly, the 2D-1D MS-VSTS denoising algorithm allows us to recover the spectral information for each spatial position, as can be seen from Fig. 3.

4. Deconvolution of spherical data with Poisson noise

We now introduce a wavelet deconvolution approach for monochannel and multichannel data on the sphere with Poisson noise. The main idea underlying the method is to apply the MS-VSTS method described above. We first introduce the deconvolution problem and then describe how the MS-VSTS can be used to solve the deconvolution problem.

4.1. Problem statement

Many problems in signal and image processing can be cast as an inversion of a linear system $\mathit{Y}\mathrm{=}\mathbf{H}\mathit{X}\mathrm{+}\mathit{\epsilon ,}$(22)where is the data to recover, is the degraded noisy observation, ε is an additive noise, and is a bounded linear operator that is typically ill-behaved because it models an acquisition process that encounters loss of information. When H is the identity, it is just a denoising problem that can be treated with the previously described methods. Inverting Eq. (22)is usually an ill-posed problem, which means that there is no unique and stable solution.

Our objective is to remove the effect of the instrument’s PSF. In our case, H is the convolution by a blurring kernel (i.e. PSF) operator that causes Y to lack the high frequency content of X. Furthermore, since the noise is Poisson, ε has a variance profile HX. The problem at hand is then a deconvolution problem in the presence of Poisson noise.

We therefore need to both regularize the problem and handle the Poisson statistics of the noise. To regularize this inversion problem and reduce the space of candidate solutions, one has to add some prior knowledge of the typical structure of the original data X. This prior information accounts for the smoothness of the solution and can range from the uniform smoothness assumption to more complex knowledge of the geometrical structures of X.

In our LAT realistic simulations, the PSF width depends strongly on the energy, from at 50 MeV to better than at 10 GeV and above. Figure 4 shows the normalized profiles of the PSF for different energy bands.

4.2. Monochannel deconvolution

We first consider the single-channel case. In the literature, several algorithms have been proposed to perform image deconvolution on a cartesian grid. The Richardson-Lucy algorithm is certainly the most famous in astrophysics. In this paper, we propose a regularized Richardson-Lucy algorithm to deconvolve data on the sphere data.

Fig. 3 Spectrum of a single gamma-ray point source recovered using the multichannel MS-VSTS denoising algorithm. Top: single gamma-ray source from simulated Fermi data integrated along the energy axis (left: simulated source; middle: noisy source; right: denoised source). Bottom: spectrum of the center of the point source: intensity as a function of the energy band with 14 energy bands between 50 MeV and 50 GeV (black: true simulated spectrum; cyan: restored spectrum from our denoising algorithm.

The Richardson-Lucy algorithm originates from a fixed-point equation obtained by maximizing the Poisson likelihood with respect to X while preserving positivity. This entails a multiplicative update rule, starting at n = 0 and X⁽⁰⁾ = 1 and iterating ${\mathit{X}}^{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{=}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\otimes }\mathrm{(}{{\mathbf{H}}^{\mathit{T}}\mathrm{\left(}\mathit{Y}\mathbf{H}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\right)}}^{\mathrm{)}}\mathit{,}$(23)where  ⊗  (respectively  $\begin{array}{}{}_{\mathrm{}}\\ \mathit{}\end{array}\mathrm{\otimes }\begin{array}{}\mathrm{}\\ \mathit{}\end{array}$) represents the element-wise multiplication (respectively division) between two vectors, and H^T is the transpose of H whose action on an image consists in convolving it with the time-reversed version of the PSF associated with H. However, it is well-known that owing to the lack of regularization, the Richardson-Lucy algorithm tends to amplify the noise after a few iterations.

We define R⁽ⁿ⁾ as the residual at iteration n${\mathit{R}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{=}\mathit{Y}\mathrm{-}\mathbf{H}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathit{,}$(24)where R⁽ⁿ⁾ can be written as the sum of its IUWTS detail subband and the last approximation subband a_J, that is, ${\mathit{R}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathrm{\right]}\mathrm{=}{\mathit{a}}_{\mathit{J}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathrm{\right]}\mathrm{+}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{J}}{\mathit{d}}_{\mathit{j}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathrm{\right]}\mathit{,} \mathrm{\forall }{\mathit{k}}_{\mathit{r}}\mathit{.}$(25)The wavelet transform provides a means of extracting only the significant structures from the residual at each iteration. With increasing number of iterations, a large part of the residual becomes statistically insignificant. The regularized significant residual is then, for a location k_r$\mathit{R̅}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathrm{\right]}\mathrm{=}{\mathit{a}}_{\mathit{J}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathrm{\right]}\mathrm{+}\sum _{\mathrm{\left(}\mathit{j,}{\mathit{k}}_{\mathit{r}}\mathrm{\right)}\mathrm{\in }\mathrm{ℳ}}{\mathit{d}}_{\mathit{j}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathrm{\right]}\mathit{,}$(26)where ℳ is the multiresolution support defined in a similar way to Eq. (19). The regularized Richardson-Lucy scheme then becomes ${\mathit{X}}^{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{=}{\mathbf{P}}_{\mathrm{+}}\mathrm{(}{{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\otimes }\mathrm{(}{{\mathbf{H}}^{\mathit{T}}\mathrm{(}{\mathrm{(}\mathbf{H}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{+}\mathit{R̅}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{)}\mathbf{H}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}}^{\mathrm{)}}}^{\mathrm{)}}}^{\mathrm{)}}\mathit{.}$(27)This algorithm is similar to Murtagh et al. (1995), except that the à trous wavelet transform is replaced by the undecimated wavelet transform. In the next subsection, we describe how the same algorithm can be extended to the multi-channel case.

Fig. 4 Normalized profile of the PSF for different energy bands as a function of the angle in degree. Black: 50−82 MeV. Cyan: 220−360 MeV. Orange: 960−1.6 GeV. Blue: 4.2−6.9 GeV. Green: 19−31 GeV.

Fig. 5 Result of the multichannel deconvolution algorithm for different energy bands. Top: simulated (blurred) intensity skymap. Middle: blurred and noisy skymap. Bottom: deconvolved skymap. Energy band: 82–134 MeV. Maps are on a logarithmic scale.

Fig. 6 Result of the multichannel deconvolution algorithm for different energy bands. Top: simulated (blurred) intensity skymap. Middle: blurred and noisy skymap. Bottom: deconvolved skymap. Energy band: 220–360 MeV. Maps are on a logarithmic scale.

Fig. 7 Result of the multichannel deconvolution algorithm for different energy bands. Top: simulated (blurred) intensity skymap. Middle: blurred and noisy skymap. Bottom: deconvolved skymap. Energy band: 360–589 MeV. Maps are on a logarithmic scale.

Fig. 8 Result of the multichannel deconvolution algorithm for different energy bands. Top: simulated (blurred) intensity skymap. Middle: blurred and noisy skymap. Bottom: deconvolved skymap. Energy band: 589–965 MeV. Maps are on a logarithmic scale.

4.3. Multichannel deconvolution

As the PSF is channel-dependent, the convolution observation model is now $\mathit{Y}\mathrm{[}\mathrm{·}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{]}\mathrm{=}{\mathbf{H}}_{{\mathit{k}}_{\mathit{t}}}\mathit{X}\mathrm{[}\mathrm{·}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{]}\mathrm{+}\mathit{\epsilon }\mathrm{[}\mathrm{·}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{]}$in each channel k_t, where H_kt is the (spatial) convolution operator in channel k_t with known PSF.

The same recipe as in the monochannel case applies with the notable difference that the spherical 2D-1D MS-VSTS is used instead of its monochannel counterpart. The multichannel multiresolution support ℳ is obtained after thresholding these coefficients.

We now define H to be the multichannel convolution¹ operator, which acts on a 2D-1D multichannel spherical data set X by applying H_t to each channel X [·,k_t]  independently². The regularized multichannel Richardson-Lucy scheme is then ${\mathit{X}}^{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{=}{\mathbf{P}}_{\mathrm{+}}\mathrm{(}{{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\otimes }\mathrm{(}{{\mathbf{H}}^{\mathit{T}}\mathrm{(}{\mathrm{(}\mathbf{H}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{+}\mathit{R̅}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{)}\mathbf{H}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}}^{\mathrm{)}}}^{\mathrm{)}}}^{\mathrm{)}}\mathit{,}$(28)where is the regularized (significant) residual $\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\begin{array}{}\end{array}\mathit{R̅}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{=}{\mathit{a}}_{{\mathit{J}}_{\mathrm{1}}\mathit{,}{\mathit{J}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathrm{+}\sum _{\mathrm{\left(}{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right)}\mathrm{\in }\mathrm{ℳ}}{\mathit{w}}_{{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}\mathrm{\left[}{\mathit{k}}_{\mathit{r}}\mathit{,}{\mathit{k}}_{\mathit{t}}\mathrm{\right]}\mathit{.}\end{array}\begin{array}{}\end{array}\end{array}$(29)

4.4. Experiments

The algorithm was applied to the seven energy bands (50 MeV–1.58 GeV) of our simulated Fermi data set. Figures 5 to 8 display the deconvolution results for four energy bands. Figure 9 shows the performance of the multichannel MS-VSTS deconvolution algorithm for a single point source. The deconvolution not only effectively removes the blur and recovers sharply localized point sources, but also allows us to restore all of the spectral information. To get a better visual impression of the performance of the deconvolution algorithm, Fig. 10 depicts the result of the algorithm on a single HEALPix face covering the Galactic plane. We find that the deconvolution is remarkably effective: our MS-VSTS multichannel deconvolution algorithm manages to remove a large part of the blur introduced by the PSF.

5. Conclusion

This paper extends the MS-VSTS framework to deal with monochannel deconvolution, multichannel denoising and multichannel deconvolution. Unlike the monochannel MS-VSTS, the multichannel MS-VSTS fully exploits the information in the 2D-1D data set and allows us to recover the spectral information on the sources. As the PSF strongly depends on the energy, it is very important to have a multichannel method for deconvolution. Multichannel deconvolution using MS-VSTS removes a large part of the PSF blur and significantly improves the sharpness of the spatial localization of point sources.

The software related to this paper, MRS/Poisson, and its full documentation will be included in the next version of ISAP (Interactive Sparse astronomical data Analysis Packages) via the ISAP web site³.

Fig. 9 Profile of a single gamma-ray point source recovered using the multichannel MS-VSTS deconvolution algorithm. Top: single gamma-ray point source on simulated (blurred) Fermi data (energy band: 360–589 MeV) (left: simulated blurred source; middle: blurred noisy source; right: deconvolved source). Bottom: profile of the point source (cyan: simulated spectrum; black: restored spectrum from the deconvolved source.

Fig. 10 View on a single HEALPix face. Result of an application of the deconvolution algorithm to the Galactic plane. Top left: simulated Fermi Poisson intensity. Top right: simulated Fermi noisy data. Bottom: Fermi data deconvolved with multichannel MS-VSTS. Energy band: 360–589 MeV. Pictures are on a logarithmic scale.

Acknowledgments

Some of the results in this paper have been derived using the Healpix (Górski et al. 2005) and the MR/S software (Starck et al. 2006). This work was supported by the French National Agency for Research (ANR-08-EMER-009-01) and the European Research Council grant SparseAstro (ERC-228261).

References

All Figures

	Fig. 1 Result of the multichannel Poisson denoising algorithm on simulated Fermi data over the energy band 220–360 MeV. Top: simulated intensity skymap. Middle: simulated noisy skymap. Bottom: denoised skymap. Maps are on a logarithmic scale.
In the text

	Fig. 2 Result of the multichannel Poisson denoising algorithm on simulated Fermi data over the energy band 589–965 MeV. Top: simulated intensity skymap. Middle: simulated noisy skymap. Bottom: denoised skymap. Maps are on a logarithmic scale.
In the text

Fig. 3 Spectrum of a single gamma-ray point source recovered using the multichannel MS-VSTS denoising algorithm. Top: single gamma-ray source from simulated Fermi data integrated along the energy axis (left: simulated source; middle: noisy source; right: denoised source). Bottom: spectrum of the center of the point source: intensity as a function of the energy band with 14 energy bands between 50 MeV and 50 GeV (black: true simulated spectrum; cyan: restored spectrum from our denoising algorithm.

In the text

	Fig. 4 Normalized profile of the PSF for different energy bands as a function of the angle in degree. Black: 50−82 MeV. Cyan: 220−360 MeV. Orange: 960−1.6 GeV. Blue: 4.2−6.9 GeV. Green: 19−31 GeV.
In the text

	Fig. 5 Result of the multichannel deconvolution algorithm for different energy bands. Top: simulated (blurred) intensity skymap. Middle: blurred and noisy skymap. Bottom: deconvolved skymap. Energy band: 82–134 MeV. Maps are on a logarithmic scale.
In the text

	Fig. 6 Result of the multichannel deconvolution algorithm for different energy bands. Top: simulated (blurred) intensity skymap. Middle: blurred and noisy skymap. Bottom: deconvolved skymap. Energy band: 220–360 MeV. Maps are on a logarithmic scale.
In the text

	Fig. 7 Result of the multichannel deconvolution algorithm for different energy bands. Top: simulated (blurred) intensity skymap. Middle: blurred and noisy skymap. Bottom: deconvolved skymap. Energy band: 360–589 MeV. Maps are on a logarithmic scale.
In the text

	Fig. 8 Result of the multichannel deconvolution algorithm for different energy bands. Top: simulated (blurred) intensity skymap. Middle: blurred and noisy skymap. Bottom: deconvolved skymap. Energy band: 589–965 MeV. Maps are on a logarithmic scale.
In the text

Fig. 9 Profile of a single gamma-ray point source recovered using the multichannel MS-VSTS deconvolution algorithm. Top: single gamma-ray point source on simulated (blurred) Fermi data (energy band: 360–589 MeV) (left: simulated blurred source; middle: blurred noisy source; right: deconvolved source). Bottom: profile of the point source (cyan: simulated spectrum; black: restored spectrum from the deconvolved source.

In the text

	Fig. 10 View on a single HEALPix face. Result of an application of the deconvolution algorithm to the Galactic plane. Top left: simulated Fermi Poisson intensity. Top right: simulated Fermi noisy data. Bottom: Fermi data deconvolved with multichannel MS-VSTS. Energy band: 360–589 MeV. Pictures are on a logarithmic scale.
In the text