On a three-dimensional Compton scattering tomography system with fixed source

J Cebeiro¹, C Tarpau^2,3,4, M A Morvidone¹, D Rubio¹ and M K Nguyen² ¹ Centro de Matemática Aplicada, Universidad Nacional de San Martín (UNSAM), Buenos Aires, Argentina. ² Equipes Traitement de l’Information et Systèmes, CY Cergy Paris University, ENSEA, CNRS UMR 8051, Cergy-Pontoise, France. ³ Laboratoire de Physique Théorique et Modélisation, CY Cergy Paris University, CNRS UMR 8089, Cergy-Pontoise, France. ⁴ Laboratoire de Mathématiques de Versailles, Université de Versailles Saint-Quentin, CNRS UMR 8100, Versailles, France. [email protected]

Abstract

Compton scattering tomography is an emerging scanning technique with attractive applications in several fields such as non-destructive testing and medical imaging. In this paper, we study a modality in three dimensions that employs a fixed source and a single detector moving on a spherical surface. We also study the Radon transform modeling the data that consists of integrals on toric surfaces. Using spherical harmonics we arrive to a generalized Abel’s type equation connecting the coefficients of the expansion of the data with those of the function. We show the uniqueness of its solution and so the invertibility of the toric Radon transform. We illustrate this through numerical reconstructions in three dimensions using a regularized approach. A preliminary version of the algorithm for discrete spherical harmonic expansion is available in the code repository [35].

1 Introduction

Compton scattering imaging takes advantage of scattered radiation in order to explore the interior of an object under study. The basis of Compton scattering imaging is Compton effect where a photon of energy $E_{0}$ interacts with an electron undergoing a change in its direction and losing part of its of energy. The angular deviation $\omega$ experienced by the photon and its remaining energy $E(\omega)$ are related by the well-known Compton formula

E(\omega)=\frac{E_{0}}{1+\frac{E_{0}}{mc^{2}}(1-\cos\omega)},

where $E_{0}$ is the energy of the incident photon and $mc^{2}$ is the energy of the electron at rest (511 keV). Transmission Compton scattering tomography consists in illuminating the object under study with a source of radiation and registering the scattered photons in its neighborhood. In this case the function to recover is a map of the electronic density of the object [1]. The interest in Compton scattering imaging is supported by its numerous potential applications namely non-destructive testing [2], airport security [3] and the study of composite materials, for instance those used in aircraft industry [4]. Recent works in medical imaging suggest potential advantages in diagnosis since images may exhibit better contrast in certain scenarios such as lung tumors [5, 6].

The first formulation of Compton tomography based on a Radon type transform has been introduced by Norton in [1]. The mathematical foundations arising from the seminal works of Allan M. Cormack on circular Radon transforms provide a formal framework to study these modalities [7, 8, 9]. Depending on the setup, forward models are based on integral transforms on different manifolds. In a two dimensional configuration, where collimated sources are used to restrict photons to a given scanning plane, manifolds are circular arcs [1, 2, 10, 11, 12, 13, 14, 15]. If uncollimated detectors are used, manifolds consist in toric sections [3, 16, 17].

In three-dimensional Compton scattering tomography, where uncollimated source and detectors are employed, photons recorded at a given place and energy carry information about the electronic density all over a toric surface. In these setups scanning is performed in a volume rather than in a plane and the stacking of planar slices after reconstruction is no longer necessary. The problem in three dimensions has recently been addressed in [18] where a comprehensive framework was introduced together with four three-dimensional modalities. This framework takes attenuation into account and provides a back-projection algorithm for contour reconstruction. For some of these configurations, there are also results on injectivity, uniqueness and numerical reconstructions of the electronic density [19, 3].

The design of an operating Compton scattering tomography system raises numerous challenges involving size, detectors, sources, shielding and precise mechanisms to move components, sometimes in a synchronous fashion. In this respect, the suppression of the need of synchronized movements, the ability to scan a large object without surrounding it and the reduction of the number of moving components are advantageous design features. For instance, employing one fixed source turns valuable, particularly when using electrically powered sources that are expected to be more difficult to handle because of its weight and wiring.

Here, we study a Compton scattering tomography in three dimensions with one fixed central source and unique detector moving on a sphere. In particular, we show the invertibility of the toric transform modeling data acquisition. Instead of forcing the object under study to be strictly inside the detection sphere [18], we placed it strictly outside. Through this relative positioning of the specimen, the scanning can be carried out using a compact system without the need of enclosing or surrounding the object. This makes the configuration suitable for scanning large bodies or objects attached to major structures. Recently, a two-dimensional reduction of this model has been studied in [17]. This reduction enabled to solve the reconstruction problem through an analytic reconstruction formula. Nevertheless, the setup requires translational relative movement between the object and the source in order to scan planar sections to be stacked in a three-dimensional reconstruction. There are 3D CST modalities for which object reconstruction has been shown feasible by means of the invertibility of its Radon transform. Two configurations, both registering back-scattered photons [19, 3], are suitable for scanning an object without surrounding it. While [19] requires the rotation of the pair source-detector, the design in [3] requires translational movement of the source that can be alternatively accomplished employing an array of multiple sources. Regarding the performance in terms of photon counts, the efficiency of the configuration studied here may be comparable to these setups since all of them employ the same type of photons (back-scattered). There is also a setup in [19] that employs forward-scattered photons to scan a body inside the sphere enclosed by the source-detector pair. In this case, comparative efficiency depends upon the range of energies considered.

In order to model the forward problem, we formulate a Radon transform on tori and study its properties, particularly those connected to the uniqueness of its solution. The paper is organized as follows. In section 2 we explain briefly the setup, the scanning protocol and the toric Radon transform used to model image formation. We also write the transform as a spherical harmonics expansion and obtain an explicit Abel’s type equation relating the coefficients of the expansion of the data with those of the function that represents the electronic density of the sought function. In section 3, we analyze its kernel and show the uniqueness of its solution proving, thus, the invertibility of the toric Radon transform. Numerical reconstructions based on discrete spherical harmonics and Tikhonov regularization are shown in section 4. We discuss the results and the specific advantages and limitations of the modality as well as some new perspectives of our work in section 5. Finally, we close the paper with some conclusions. Details of the spherical harmonics expansion and the discrete algorithm we used to implement it are given in A and B.

2 A setup for three dimensional Compton scattering tomography

We introduce a setup for Compton scattering tomography based on a central source of radiation $S$ located at the origin of coordinates. A single detector $D$ moves at a constant distance $R$ of the origin describing a sphere, as shown in figure 1. The property to recover is the electronic density $f(x,y,z)$ of an object placed strictly outside of this spherical surface. Detector $D$ registers backscattered photons of energy $E(\omega)$ that have interacted with electrons of the object and deviated an angle $\omega$ from its original path. When detector $D$ is at position labeled by angles $(\alpha,\beta)$ and registers a photon of energy $E(\omega)$ , the interaction site is located somewhere on the surface of an apple torus through points $S$ and $D$ , which is characterized by $\alpha,\,\beta$ and $\omega$ . Thus, the flux of photons of energy $E(\omega)$ registered by $D$ at sites $(\alpha,\beta)$ is proportional to the integral of function $f$ on a toric surface $\mathcal{T}^{\omega,\alpha,\beta}$ , with $\alpha\in[0,2\pi)$ , $\beta\in[0,\pi]$ and $\omega\in(\pi/2,\pi)$ , see figures 1 and 2. Notice that the referred torus may be generated by the rotation of a toric section around segment $\overline{SD}$ , see figure 2.

Refer to caption — Figure 1: a: setup of the new CST modality: source ( $S$ ), detector ( $D$ ), detection sphere (mesh). b and c: two views of an apple torus $\mathcal{T}$ (blue) and the detection sphere (yellow).

Before giving the definition of the Radon transform modeling the forward problem, let us introduce some notation. The coordinates of the detector $D$ are given by

D(\alpha,\beta)=R(\begin{array}[]{c}\cos\alpha\sin\beta,\sin\alpha\sin\beta,\cos\beta\end{array})^{T},\quad\alpha\in[0,\,2\pi),\>\beta\in[0,\pi].

We define a parametrization of any torus $\mathcal{T}^{\omega,\alpha,\beta}$ as:

\Phi^{\omega,\alpha,\beta}(\gamma,\psi)=r^{\omega}(\gamma)\Theta^{\alpha,\beta}(\gamma,\psi),

(1)

with $\gamma\in(0,2\omega-\pi)$ , $\psi\in(0,2\pi)$ . Here, the radial part is given by

r^{\omega}(\gamma)=R\frac{\sin(\omega-\gamma)}{\sin\omega}

while the angular part is expressed as $\Theta^{\alpha,\beta}(\gamma,\psi)=u(\alpha)a(\beta)\Theta(\gamma,\psi),$ where

\Theta(\gamma,\psi)=(\begin{array}[]{c}\cos\psi\sin\gamma,\sin\psi\sin\gamma,\cos\gamma\end{array})^{T}

is a point on $S^{2}$ , the unit sphere in $\mathbb{R}^{3}$ , and

u(\alpha)=\left(\begin{array}[]{ccc}\cos\alpha&-\sin\alpha&0\\ \sin\alpha&\,\,\,\cos\alpha&0\\ 0&0&1\end{array}\right)\,\,\textrm{and}

a(\beta)=\left(\begin{array}[]{ccc}\cos\beta&0&\sin\beta\\ 0&1&0\\ -\sin\beta&0&\cos\beta\\ \end{array}\right)

(2)

are a rotation of angle $\alpha$ about the $z$ -axis and a rotation of angle $\beta$ about the $y$ -axis, respectively. Finally, given positive numbers $r_{M},\,r_{m}$ such that $r_{M}>r_{m}>R>0$ , we define the spherical shell $S_{h}(r_{m},r_{M})$ as

S_{h}(r_{m},r_{M})=\{(x,y,z)\in\mathbb{R}^{3}:r_{m}\leq\sqrt{x^{2}+y^{2}+z^{2}}\leq r_{M}\}.

We are now ready to introduce the toric Radon transform associated to this modality.

Definition. Let $f(x,y,z)$ be a compact supported function with support contained in $S_{h}(r_{m},r_{M})$ . We define the toric Radon transform $\mathcal{R}_{\mathcal{T}}f$ of function $f$ as $\mathcal{R}_{\mathcal{T}}f(\alpha,\beta,\omega)=\int_{\mathcal{T^{\omega,\alpha,\beta}}}dS_{\mathcal{T}}\ f(x,y,z)$ . Explicitly

\mathcal{R}_{\mathcal{T}}f(\alpha,\beta,\omega)=\int_{0}^{2\omega-\pi}d\gamma\int_{0}^{2\pi}d\psi\ f\left(\Phi^{\omega,\alpha,\beta}(\gamma,\psi)\right)\,r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega},

(3)

where $\omega\in(\frac{\pi}{2},\pi)$ , $\alpha\in[0,2\pi)$ , $\beta\in[0,\pi]$ .

Equation (3) is the forward operator that models the data recorded, this integral transform is rotational invariant. Its manifold is the part of an apple torus through the origin outside of the sphere of radius $R$ . In [19], the authors studied another rotational invariant toric transform whose manifold was an apple torus that did not pass through the origin. Interesting issues on the toric Radon transform (3) such as spherical harmonics expansion [20], rotational invariance and invertibility [21] can be studied, in what follows we address some of them.

2.1 Spherical harmonics expansion

In this section, we explicit the connection between the components of the spherical harmonics expansion of a function $f$ and those of its toric Radon transform $\mathcal{R}_{\mathcal{T}}f$ . Spherical harmonics of degree $l$ and order $m$ are defined as:

Y_{l}^{m}(\gamma,\psi)=(-1)^{m}\sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}}P^{m}_{l}(\cos\gamma)e^{im\psi},

(4)

where $\gamma\in[0,\pi]$ , $\psi\in[0,2\pi)$ and $P^{m}_{l}(x)$ is the Legendre polynomial of degree $l$ and order $m$ , see [22, 20] for details. The set $\{Y_{l}^{m}\}$ , for $l\in\mathbb{N}$ and $|m|\leq l$ is a complete orthonormal system in $S^{2}.$ Any function $f\in C^{\infty}(\mathbb{R}^{3})$ can be expanded in terms of $Y_{l}^{m}(\gamma,\psi)$ according to

f(r\Theta(\gamma,\psi))=\sum_{l=0}^{\infty}\sum_{|m|\leq l}f_{lm}(r)Y_{l}^{m}(\gamma,\psi),

(5)

where

f_{lm}(r)=\langle f,Y_{l}^{m}\rangle=\int_{0}^{2\pi}\int_{0}^{\pi}f(r\Theta(\gamma,\psi))\overline{Y_{l}^{m}(\gamma,\psi)}\sin\gamma\,d\gamma d\psi,

and $\langle,\rangle$ is the scalar product in $L^{2}(S^{2})$ , the overline denotes complex conjugation, and $r\in\mathbb{R}^{+}_{0}$ is fixed.

Following the ideas in [19], it can be shown that the spherical harmonics expansion of $\mathcal{R}_{\mathcal{T}}f$ can be written as

\mathcal{R}_{\mathcal{T}}f(\alpha,\beta,\omega)=\sum_{l\in\mathbb{N}}\sum_{|m|\leq l}(\mathcal{R}_{\mathcal{T}}f)_{lm}(\omega)Y_{l}^{m}(\alpha,\beta),

(6)

where coefficients in the expansion are given by the following lemma, see A for a proof.

Lemma 2.1.

The Fourier coefficients of data $(\mathcal{R}_{\mathcal{T}}f)_{lm}$ in (6) and those of the object $f_{lm}$ are related by

(\mathcal{R}_{\mathcal{T}}f)_{lm}(\omega)=2\pi\int_{0}^{2\omega-\pi}d\gamma\,r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,f_{lm}\left(r^{\omega}(\gamma)\right)P^{0}_{l}(\cos\gamma),

(7)

where $P^{0}_{l}(.)$ is the zero order associated Legendre polynomial of degree $l$ .

3 Uniqueness of the solution of the toric Radon transform

We address now the question of the uniqueness of solutions of (3) using spherical harmonics expansions. In section 3.1 we show that (7) is a generalized Abel type integral equation with a kernel with zeros on its diagonal. In section 3.2 we analyze the properties of the kernel to prove invertibility.

3.1 An alternative expression for the integral equation

We split up (7) in two parts having integration range $(0,{\omega-\frac{\pi}{2}})$ and $({\omega-\frac{\pi}{2}},2\omega-\pi)$ , respectively and perform the substitution $\gamma=\omega-\sin^{-1}\left(\frac{r\sin\omega}{R}\right)$ in the first integral and $\gamma=\omega+\sin^{-1}\left(\frac{r\sin\omega}{R}\right)-\pi$ in the second one. After some calculations, making use of the identities $\cos(\sin^{-1}x)=\sqrt{1-x^{2}}$ and $P_{l}^{0}(-x)=(-1)^{l}P_{l}^{0}(x)$ (see, for example, [23]), we arrive to

	$\displaystyle(\mathcal{R}_{\mathcal{T}}f)_{lm}(\omega)=\int_{R}^{\frac{R}{\sin\omega}}dr\,f_{lm}\left(r\right)\frac{1}{\sqrt{(\frac{R}{\sin\omega})-r}}\times$
	$\displaystyle\frac{2\pi}{\sin\omega}{\sum_{\sigma=\pm 1}}\frac{r\sin\left(\sin^{-1}\left(\frac{r\sin\omega}{R}\right)-\sigma\omega\right)}{\sqrt{(\frac{R}{\sin\omega})+r}}(\sigma)^{l}P^{0}_{l}\left(\cos\left(\omega-\sigma\sin^{-1}\left(\frac{r\sin\omega}{R}\right)\right)\right).$
			(8)

The support of function $f_{lm}(r)$ enables to replace the lower integration limit $R$ by $r_{m}$ . Making the substitution $p=R/\sin\omega$ ¹¹1Physically $p$ is the diameter of the circles that generate the torus as a surface of revolution. and keeping the notation for readability²²2Strictly speaking, we should have changed the function after this substitution, for instance $(\mathcal{R}_{\mathcal{T}}f)_{lm}(\sin^{-1}R/p)=\widetilde{(\mathcal{R}_{\mathcal{T}}}f)_{lm}(p)$ ., the integral equation reads:

\displaystyle{(\mathcal{R}_{\mathcal{T}}f)_{lm}}(p)=\int_{r_{m}}^{p}dr\,f_{lm}\left(r\right)\frac{1}{\sqrt{p-r}}K_{l}(p,r),

(9)

where the kernel is

\displaystyle K_{l}(p,r)=\frac{2\pi}{R}{\sum_{\sigma=\pm 1}}\sigma^{l}\frac{p\,r\sin\left(\sin^{-1}\left(\frac{r}{p}\right)-\sigma\sin^{-1}\left(\frac{R}{p}\right)\right)}{\sqrt{p+r}}P^{0}_{l}\left(\cos\left(\sin^{-1}\left(\frac{R}{p}\right)-\sigma\sin^{-1}\left(\frac{r}{p}\right)\right)\right).

Now, we rewrite kernel (3.1) in a suitable form for the evaluation of the existence of a unique solution of (7). Using elementary geometry, we can express the kernel as

	$\displaystyle K_{l}(p,r)={\sum_{\sigma=\pm 1}}\left(Q_{1}-\sigma Q_{2}\sqrt{p-r}\right)\sigma^{l}P^{0}_{l}\left(Q_{3}\sqrt{p-r}+\sigma Q_{4}\right),$
			(11)

where

	$\displaystyle Q_{1}=Q_{1}(p,r)=\frac{2\pi r^{2}\sqrt{p^{2}-R^{2}}}{Rp\sqrt{p+r}}$		(12)
	$\displaystyle Q_{2}=Q_{2}(p,r)=\frac{2\pi r}{p}$		(13)
	$\displaystyle Q_{3}=Q_{3}(p,r)=\frac{\sqrt{p^{2}-R^{2}}\sqrt{p+r}}{p^{2}}$		(14)
	$\displaystyle Q_{4}=Q_{4}(p,r)=\frac{Rr}{p^{2}}.$		(15)

Using the binomial formula to expand each term of $P^{0}_{l}(x)=a_{l}x^{l}+...+a_{1}x+a_{0},$ with $x=Q_{3}\sqrt{p-r}+\sigma Q_{4}$ , we get

	$\displaystyle K_{l}(p,r)={\sum_{n=0}^{l}}a_{n}{\sum_{k=0}^{n}}{{n}\choose{k}}Q_{3}^{n-k}Q_{4}^{k}\left[\left((-1)^{n-k}+1\right)Q_{1}\sqrt{p-r}^{n-k}+\left((-1)^{n-k}-1\right)Q_{2}\sqrt{p-r}^{n-k+1}\right].$
			(16)

We perform some calculations

	$\displaystyle K_{l}(p,r)={\sum_{n=0}^{3}}a_{n}{\sum_{k=0}^{n}}{{n}\choose{k}}Q_{3}^{n-k}Q_{4}^{k}\left[\left((-1)^{n-k}+1\right)Q_{1}\sqrt{p-r}^{n-k}+\left((-1)^{n-k}-1\right)Q_{2}\sqrt{p-r}^{n-k+1}\right]$
	$\displaystyle\phantom{doremi}+\sum_{n=4}^{l}a_{n}\left\{{\sum_{k=0}^{n-4}}{{n}\choose{k}}Q_{3}^{n-k}Q_{4}^{k}\left[\left((-1)^{n-k}+1\right)Q_{1}\sqrt{p-r}^{n-k}+\left((-1)^{n-k}-1\right)Q_{2}\sqrt{p-r}^{n-k+1}\right]\right.$
	$\displaystyle\phantom{doremi}\left.+{\sum_{k=n-3}^{n}}{{n}\choose{k}}Q_{3}^{n-k}Q_{4}^{k}\left[\left((-1)^{n-k}+1\right)Q_{1}\sqrt{p-r}^{n-k}+\left((-1)^{n-k}-1\right)Q_{2}\sqrt{p-r}^{n-k+1}\right]\right\}$
	$\displaystyle\phantom{doremi}=\sum_{n=4}^{l}a_{n}{\sum_{k=0}^{n-4}}{{n}\choose{k}}Q_{3}^{n-k}Q_{4}^{k}\left[\left((-1)^{n-k}+1\right)Q_{1}\sqrt{p-r}^{n-k}+\left((-1)^{n-k}-1\right)Q_{2}\sqrt{p-r}^{n-k+1}\right]$
	$\displaystyle\phantom{doremi}+(p-r)^{2}\left[-2a_{3}{{3}\choose{0}}Q_{3}^{3}\,Q_{2}-2\sum_{n=4}^{l}a_{n}{{n}\choose{n-3}}Q_{3}^{3}\,Q_{4}^{n-3}\,Q_{2}\right]$
	$\displaystyle\phantom{doremi}+2(p-r)\left\{\left[Q_{3}^{2}\,Q_{1}\sum_{n=4}^{l}a_{n}{{n}\choose{n-2}}Q_{4}^{n-2}-Q_{3}\,Q_{2}\sum_{n=4}^{l}a_{n}{{n}\choose{n-1}}Q_{4}^{n-1}\right]\right.$
	$\displaystyle\phantom{doremi}\left.+Q_{3}^{2}\,Q_{1}\left[a_{2}{{2}\choose{0}}+a_{3}{{3}\choose{1}}Q_{4}\right]-Q_{3}\,Q_{2}\left[a_{1}{{1}\choose{0}}+a_{2}{{2}\choose{1}}Q_{4}+a_{3}{{3}\choose{2}}Q_{4}^{2}\right]\right\}$
	$\displaystyle\phantom{doremi}+2Q_{1}\left[\sum_{n=4}^{l}a_{n}{{n}\choose{n}}Q_{4}^{n}+a_{3}{{3}\choose{3}}Q_{4}^{3}+a_{2}{{2}\choose{2}}Q_{4}^{2}+a_{1}{{1}\choose{1}}Q_{4}+a_{0}{{0}\choose{0}}\right],$

and finally we arrive to

	$\displaystyle K_{l}(p,r)=2Q_{1}{\sum_{n=0}^{l}}a_{n}Q_{4}^{n}+2(p-r)\left[\frac{1}{2}Q_{3}^{2}Q_{1}{\sum_{n=2}^{l}}a_{n}n(n-1)Q_{4}^{n-2}-Q_{3}Q_{2}{\sum_{n=1}^{l}}a_{n}n\,Q_{4}^{n-1}\right]$		(17)
	$\displaystyle+{\sum_{n=3}^{l}}a_{n}{\sum_{k=0}^{n-3}}{{n}\choose{k}}Q_{3}^{n-k}Q_{4}^{k}\left[\left((-1)^{n-k}+1\right)Q_{1}\sqrt{p-r}^{n-k}+\left((-1)^{n-k}-1\right)Q_{2}\sqrt{p-r}^{n-k+1}\right].$		(18)

			(19)

Notice that summations on the first line correspond to Legendre polynomials and their first and second order derivatives evaluated in $Q_{4}$ . Also, terms with odd powers of $\sqrt{p-r}$ vanish and only powers of $(p-r)$ greater than or equal to 2 remain. Then, using a compact notation for Legendre polynomials $P(.)=P^{0}_{l}(.)$ and its derivatives $P^{\prime}(.)$ and $P^{\prime\prime}(.)$ , and splitting up summations in order to define a new subindex, we have

	$\displaystyle K_{l}(p,r)=2Q_{1}P(Q_{4})+2(p-r)\left[{\textstyle\frac{1}{2}}Q_{3}^{2}Q_{1}P^{\prime\prime}(Q_{4})-Q_{3}Q_{2}P^{\prime}(Q_{4})\right]$		(20)
	$\displaystyle+2{\sum_{n=3}^{l}}a_{n}\left[Q_{1}{\sum_{i=2}^{\lfloor{\textstyle\frac{n}{2}}\rfloor}}{{n}\choose{n-2i}}Q_{3}^{2i}Q_{4}^{n-2i}(p-r)^{i}-Q_{2}{\sum_{i=2}^{\lfloor{\textstyle\frac{n+1}{2}}\rfloor}}{{n}\choose{n-2i+1}}Q_{3}^{2i-1}Q_{4}^{n-2i+1}(p-r)^{i}\right].$		(21)

Equation (20) is useful for studying the properties of the kernel in order to prove invertibility. Particularly, (20) is more suitable for differentiation on the diagonal $p=r$ according Theorem 3.1 in next section.

3.2 On the invertibility of $\mathcal{R}_{\mathcal{T}}f$ in terms of spherical harmonics

Equation (9) is a generalized Abel type equation. Notice that kernel (20) has zeros on the diagonal $p=r$ , since the Legendre polynomials have zeros in the interval $[0,1]$ . Recently, Schiefeneder and Haltmeier stated conditions for the uniqueness of the solution of this kind of equations with kernels with zeros on its diagonal [24]. For the sake of completeness, we present their results in the following lemma.

Lemma 3.1 ([24], Theorem 3.4).

Consider the generalized Abel type integral equation

\forall t\in[a,b]:g(t)=\int_{a}^{t}ds\,f(s)\frac{1}{\sqrt{t-s}}K(t,s)

(23)

where $g\in C([a,b])$ and $K\in C(\bigtriangleup(a,b))$ , with $\bigtriangleup(a,b):=\{a\leq s\leq t\leq b\}$ , is a continuous kernel having zeros on the diagonal. Suppose $K:\bigtriangleup(a,b)\rightarrow\mathbb{R}$ , where $a<b$ , satisfies the following assumptions:

1.

$K\in C^{3}(\bigtriangleup(a,b)).$
2.

$N_{K}:=\{s\in[a,b)|K(s,s)=0\}$ is finite and consists of simple roots.
3.

For every $s\in N_{K}$ , the gradient $(\kappa_{1},\kappa_{2})=\nabla K(s,s)$ satisfies

$1+\frac{1}{2}\frac{\kappa_{1}}{\kappa_{1}+\kappa_{2}}>0.$

Then, for any $g\in C([a,b])$ , equation (23) has at most one solution $f\in C([a,b])$ .

We are now ready to introduce our claim.

Theorem 3.1 (Invertibility of equation (9)).

For any $(\mathcal{R}_{\mathcal{T}}f)_{lm}\in C([r_{m},r_{M}])$ , equation (9) has at most one solution $f_{lm}\in C([r_{m},r_{M}])$ .

Proof.

In order to show the uniqueness of the solution of (9), we use lemma 3.1. So, we must check assumptions $(i)$ , $(ii)$ and $(iii)$ . We define the triangle $\bigtriangleup(r_{m},r_{M}):=\{r_{m}\leq r\leq p\leq r_{M}\}$ .

1.

As it is evident from (3.1), the kernel is smooth inside $\bigtriangleup(r_{m},r_{M})$ . Notice in (3.1) that terms containing odd powers of $\sqrt{p-r}$ always vanish. Moreover, the double summation involves only non-negative integer powers of $p-r$ times products of powers of factors $Q_{1}$ , $Q_{2}$ , $Q_{3}$ and $Q_{4}$ . Since all the functions involved are smooth on the diagonal $p=r$ , so is $K_{l}$ .
2.

From (20) we obtain the expression for the kernel on the diagonal

$K_{l}(r,r)={\textstyle\frac{\sqrt{8}\pi}{R}}{\sqrt{r}}\sqrt{r^{2}-R^{2}}\,P\left(\frac{R}{r}\right).$

So, $K_{l}(r,r)$ has the same zeros as $P\left({\textstyle\frac{R}{r}}\right)$ since $0<R<r_{m}\leq r$ . Given that $P$ is an orthogonal polynomial, $P\left({\textstyle\frac{R}{r}}\right)$ has a finite number of simple roots. Then, the values $r_{0}$ such that $P\left({\textstyle\frac{R}{r_{0}}}\right)=0$ are the zeros of the kernel on the diagonal, i.e. $K_{l}(r_{0},r_{0})=0$ .

We obtain an expression for $\nabla K_{l}(r,r)$ and evaluate it at zeros $r_{0}$ . We see in (20) that the derivatives of terms containing factors $(p-r)^{i}$ with $i=2,3,...$ vanish on the diagonal $p=r$ . So, we define a new function

	$\displaystyle\overline{K}_{l}(p,r)=2Q_{1}P(Q_{4})+2(p-r)\left[{\textstyle\frac{1}{2}}Q_{3}^{2}Q_{1}P^{\prime\prime}(Q_{4})-Q_{3}Q_{2}P^{\prime}(Q_{4})\right],$
			(24)

which verifies $\nabla K_{l}(r,r)=\nabla\overline{K}_{l}(r,r)$ and is easier to handle. Then, we differentiate $\overline{K}_{l}$ with respect to $p$

	$\displaystyle\frac{\partial\overline{K}_{l}}{\partial p}(p,r)=2\frac{\partial Q_{1}}{\partial p}P(Q_{4})+2Q_{1}\frac{\partial}{\partial p}P(Q_{4})+2\left[{\textstyle\frac{1}{2}}Q_{3}^{2}Q_{1}P^{\prime\prime}(Q_{4})-Q_{3}Q_{2}P^{\prime}(Q_{4})\right]$		(25)
	$\displaystyle+2(p-r)\frac{\partial}{\partial p}\left[{\textstyle\frac{1}{2}}Q_{3}^{2}Q_{1}P^{\prime\prime}(Q_{4})-Q_{3}Q_{2}P^{\prime}(Q_{4})\right].$
			(26)

Using the identity $P^{\prime\prime}(x)=[2xP^{\prime}(x)-l(l+1)P(x)]/(1-x^{2}),$ see for example [23, 25], and the fact that $Q_{4}(p,r)<1,$ we substitute $P^{\prime\prime}(Q_{4})$ in the third term in (25), and we arrive to

	$\displaystyle\frac{\partial\overline{K}_{l}}{\partial p}(p,r)=\left(2\frac{\partial Q_{1}}{\partial p}-{\textstyle\frac{Q_{3}^{2}Q_{1}l(l+1)}{1-Q_{4}^{2}}}\right)P(Q_{4})+2(p-r)\frac{\partial}{\partial p}\left[{\textstyle\frac{1}{2}}Q_{3}^{2}Q_{1}P^{\prime\prime}(Q_{4})-Q_{3}Q_{2}P^{\prime}(Q_{4})\right]$		(27)
	$\displaystyle+2\left[Q_{1}\frac{\partial Q_{4}}{\partial p}-Q_{3}Q_{2}+{\textstyle\frac{Q_{3}^{2}Q_{1}Q_{4}}{1-Q_{4}^{2}}}\right]P^{\prime}(Q_{4}).$
			(28)

We proceed identically for the other variable $r$

	$\displaystyle\frac{\partial\overline{K}_{l}}{\partial r}(p,r)=2\frac{\partial Q_{1}}{\partial r}P(Q_{4})+2Q_{1}\frac{\partial}{\partial r}P(Q_{4})-2\left[Q_{3}^{2}Q_{1}{\textstyle\frac{1}{2}}P^{\prime\prime}(Q_{4})-Q_{3}Q_{2}P^{\prime}(Q_{4})\right]$		(29)
	$\displaystyle+2(p-r)\frac{\partial}{\partial r}\left[{\textstyle\frac{1}{2}}Q_{3}^{2}Q_{1}P^{\prime\prime}(Q_{4})-Q_{3}Q_{2}P^{\prime}(Q_{4})\right]$
			(30)

and we obtain

	$\displaystyle\frac{\partial\overline{K}_{l}}{\partial r}(p,r)=\left(2\frac{\partial Q_{1}}{\partial r}+{\textstyle\frac{Q_{3}^{2}Q_{1}l(l+1)}{1-Q_{4}^{2}}}\right)P(Q_{4})+2(p-r)\frac{\partial}{\partial r}\left[{\textstyle\frac{1}{2}}Q_{3}^{2}Q_{1}P^{\prime\prime}(Q_{4})-Q_{3}Q_{2}P^{\prime}(Q_{4})\right]$		(31)
	$\displaystyle+2\left[Q_{1}\frac{\partial Q_{4}}{\partial r}+Q_{3}Q_{2}-{\textstyle\frac{Q_{3}^{2}Q_{1}Q_{4}}{1-Q_{4}^{2}}}\right]P^{\prime}(Q_{4}).$
			(32)

The explicit forms (27) and (31) simplify the calculation of $\frac{\partial K_{l}}{\partial r}(r,r)$ and $\frac{\partial K_{l}}{\partial p}(r,r)$ . Evaluating at $p=r=r_{0},$ the first and second terms in (27) and (31) vanish. In fact, the first terms vanish because $P(Q_{4}(r_{0},r_{0}))=P(\frac{Rr_{0}}{r_{0}^{2}})=P(\frac{R}{r_{0}})=0,$ while the second terms vanish because of the factor $(p-r).$ Thus, we have

	$\displaystyle(\kappa_{1},\kappa_{2})=\nabla K_{l}(r_{0},r_{0})=\nabla\overline{K}_{l}(r_{0},r_{0})=$		(33)
	$\displaystyle 2\left(\left[Q_{1}\frac{\partial Q_{4}}{\partial p}-Q_{3}Q_{2}+{\textstyle\frac{Q_{3}^{2}Q_{1}Q_{4}}{1-Q_{4}^{2}}}\right],\left[Q_{1}\frac{\partial Q_{4}}{\partial r}+Q_{3}Q_{2}-{\textstyle\frac{Q_{3}^{2}Q_{1}Q_{4}}{1-Q_{4}^{2}}}\right]\right)P^{\prime}(Q_{4})\biggr{\|}_{p=r=r_{0}}.$
			(34)

We evaluate (33) using the definition of functions $Q_{1},\,Q_{2},\,Q_{3}$ and $Q_{4}$ , then the expressions for the gradient components are

\kappa_{1}=-4\pi P^{\prime}\left({\textstyle\frac{R}{r_{0}}}\right){\textstyle\frac{\sqrt{2r_{0}}\sqrt{r_{0}^{2}-R^{2}}}{r_{0}^{2}}}

(35)

and

\kappa_{2}=2\pi P^{\prime}\left({\textstyle\frac{R}{r_{0}}}\right){\textstyle\frac{\sqrt{2r_{0}}\sqrt{r_{0}^{2}-R^{2}}}{r_{0}^{2}}}.

(36)

Finally, we get the ratio

1+\frac{1}{2}\frac{\kappa_{1}}{\kappa_{1}+\kappa_{2}}=2>0,

as we wanted to show.

Thus, the radial components $f_{lm}(r)$ of function $f$ can be recovered uniquely from the coefficients $(\mathcal{R}_{\mathcal{T}}f)_{lm}(p)$ .

∎

We can express data in terms of variable $p$ instead of the scattering angle $\omega$

\mathcal{R}_{\mathcal{T}}f(p,\alpha,\beta)=\sum_{l\in\mathbb{N}}\sum_{|m|\leq l}(\mathcal{R}_{\mathcal{T}}f)_{lm}(p)Y_{l}^{m}(\beta,\alpha).

From the geometrical point of view, $p$ is the diameter of the circles generating the torus, see section 4.1 for the integral definition of $\mathcal{R}_{\mathcal{T}}f(p,\alpha,\beta)$ . The following corollary summarizes our result on uniqueness.

Corollary 3.1 (Invertibility of the $\mathcal{R}_{\mathcal{T}}f$ ).

If $f_{1}$ and $f_{2}$ are compact supported functions in $C^{\infty}(S_{h}(r_{m},r_{M}))$ and $\mathcal{R}_{\mathcal{T}}f_{1}=\mathcal{R}_{\mathcal{T}}f_{2}$ , then $f_{1}=f_{2}$ .

Proof.

Let $f$ satisfy $(\mathcal{R}_{\mathcal{T}}f)_{lm}=0$ for all $l,m$ in the spherical harmonics expansion. According to Theorem 3.1, there is a unique solution $f_{lm}=0$ , which implies $f=0$ . The linearity of $\mathcal{R}_{\mathcal{T}}f$ yields the claim. ∎

4 Numerical Simulations

We use a discrete spherical harmonics expansion of the functions representing the data and the object in order to carry out numerical reconstructions. We employ a product integration approach to model the discrete problem in the spherical harmonic domain. Then, Tikhonov regularization is used to solve a set of normal equations.

4.1 Alternative forward model

An alternative definition of our toric Radon transform useful in numerical simulation is

\mathcal{R}_{\mathcal{T}}f(p,\alpha,\beta)=p^{2}\int_{0}^{\pi}d\gamma\int_{0}^{2\pi}d\psi\cos\left(\gamma-\cos^{-1}\frac{R}{p}\right)\sin\gamma f\left(\Phi^{p,\alpha,\beta}(\gamma,\psi)\right),

(37)

where $p\in(R,+\infty)$ is the diameter of the circles making the torus and

\Phi^{p,\alpha,\beta}(\gamma,\psi)=r\Theta^{\alpha,\beta}(\gamma,\psi)|_{r=p\cos(\gamma-\cos^{-1}\frac{R}{p})}

(38)

is the parametrization of the toric surface labeled by variables $(p,\alpha,\beta)$ . The unit vector $\Theta^{\alpha,\beta}(\gamma,\psi)$ was introduced in the definition (1). The scattering angle $\omega$ and the diameter $p$ are related through $p=R/\sin\omega$ .

4.2 The algebraic problem

In this section we use a discrete spherical harmonics expansion of order $N$ to write the problem as an algebraic product suitable for Tikhonov regularization. This is achieved using numerical algorithms for the Discrete-Inverse Spherical Harmonics transform (DSHT-IDSHT), see an outline in B, with data as well as with the sought function according to equations (39) and (40). The pair DSHT-IDSHT allows us to write the problem in the domain of the spherical harmonics according to

\mathbf{g}_{nk}^{j}\mathop{\rightleftarrows}^{\mathrm{DSHT}}_{\mathrm{IDSHT}}\mathbf{g}_{lm}^{j}

(39)

\mathbf{f}_{nk}^{i}\mathop{\rightleftarrows}^{\mathrm{DSHT}}_{\mathrm{IDSHT}}\mathbf{f}_{lm}^{i}

(40)

where the discrete functions are $\mathbf{g}_{nk}^{j}=\mathcal{R}_{\mathcal{T}}f(p_{j},\alpha_{n},\beta_{k})$ , $\mathbf{g}_{lm}^{j}=\left(\mathcal{R}_{\mathcal{T}}f\right)_{lm}(p_{j})$ with $j=0,...,N_{p}-1$ , $\mathbf{f}_{nk}^{i}=f(r_{i}\cos\psi_{n}\sin\gamma_{k},r_{i}\sin\psi_{n}\sin\gamma_{k},r_{i}\cos\gamma_{k})$ and $\mathbf{f}_{lm}^{i}=f_{lm}(r_{i})$ with $i=0,...,N_{r}-1$ . Variables $p$ and $r$ are in the range $(R,r_{M}^{*}]$ where $r_{M}^{*}$ is such that $r_{M}^{*}\geq r_{M}$ to cover the radial support of function $f$ . In this domain, the components of the vector representing the unknown function $\mathbf{f}_{lm}$ are related to components of the vector for the known data $\mathbf{g}_{lm}$ through the equation $\mathbf{g}_{lm}=A_{l}\mathbf{f}_{lm}$ that is an algebraic relative of equation (9). We are aimed to solve the equation for each combination $l,m$ .

4.3 Matrix generation

Matrix $A_{l}\in\mathbb{R}^{N_{p}\times N_{r}}$ is the key for solving the numerical inverse problem. Given that we use an expansion of order $N$ and the kernel in (9) is $l$ -dependent and $m$ -independent, there are only $N+1$ different matrices. Adopting the convention $Np=Nr=M$ and splitting up the integration range in equation (9), we rewrite

\left(\mathcal{R}_{\mathcal{T}}f\right)_{lm}(p_{j})=\sum_{q=1}^{j}\int_{r_{q-1}}^{r_{q}}dr\,f_{lm}(r)\frac{r}{\sqrt{p_{j}^{2}-r^{2}}}\tilde{K_{l}}(p_{j},r)

(41)

where $r_{q}=R+q(r_{M}^{*}-R)/M$ and $\tilde{K_{l}}(p_{j},r)=\sqrt{p_{j}+r}{K_{l}}(p_{j},r)/r$ . We use product integration [26] but instead of using the mid-point rule we approximate $\tilde{K_{l}}(p_{j},r)$ by its average $\tilde{K_{l}}_{q}(p_{j})$ in ten equidistant points in the interval $[r_{q-1},r_{q}]$ . Thus, the discrete form for the equation is

\mathbf{g}_{lm}^{j}=\left(\mathcal{R}_{\mathcal{T}}f\right)_{lm}(p_{j})\simeq\sum_{q=1}^{j}w_{j,q}\tilde{K_{l}}_{q}(p_{j})f_{lm}(r_{q}),

(42)

where the weighting factor has been calculated analytically according to

w_{j,q}:=\int_{r_{q-1}}^{r_{q}}dr\frac{r}{\sqrt{p_{j}^{2}-r^{2}}}

(43)

with $j,q=1,...,M$ and $w_{j,q}=0$ if $j<q$ . Finally, the entries of the lower-triangular matrix in equation $\mathbf{g}_{lm}=A_{l}\mathbf{f}_{lm}$ are

A_{l}=\left(w_{j,q}\tilde{K_{l}}_{q}(p_{j})\right)_{j,q=1,....M}\in\mathbb{R}^{M\times M}.

(44)

4.4 Overview of the reconstruction algorithm

The aim now is recovering vector $\mathbf{f}_{lm}$ from $\mathbf{g}_{lm}$ using the equation $\mathbf{g}_{lm}=A_{l}\mathbf{f}_{lm}$ . When $A_{l}$ is non-singular and well-conditioned the problem can be easily solved by forward substitution. Because the kernel has zeros on its diagonal, matrix $A_{l}$ may have diagonal entries being zero or close to zero. Thus, solving the system may be ill-conditioned and regularization methods must be applied. The algorithm is summarized in Table 1 and is aimed to solve the matrix problem $\mathbf{g}_{lm}=A_{l}\mathbf{f}_{lm}$ for $l=0,...,N$ and $|m|\leq l$ . Tikhonov regularization requires to solve the normal equations

(A_{l}^{T}A_{l}+\lambda I)\mathbf{f}_{lm}=A_{l}^{T}\mathbf{g}_{lm},

(45)

where $I$ is the identity matrix and $\lambda$ is a regularization parameter.

Table 1: Sketch of the reconstruction algorithm.

Matrix $A_{l}$ is precalculated according to (44).
1: Chose a suitable value for $\lambda$ according to signal conditions
2: Perform DSHT (39) to data $\mathbf{g}_{nk}^{j}$ to obtain $\mathbf{g}_{lm}^{j}$
3: For each pair $l,m$ : solve the normal equations in (45)
and obtain an approximation of $\mathbf{f}_{lm}^{i}$
4: Perform IDSHT (40) to obtain the reconstruction $\mathbf{f}_{nk}^{i}$
5: Interpolate to obtain the function in discrete Cartesian coordinates $\tilde{\mathbf{f}}$

For non-vanishing kernel diagonals, the product integration method is convergent [27]. To the best of our knowledge, there are no reported results on the convergence of product integration with vanishing kernel-diagonals. As suggested in [26], numerical evidence indicate that, for a suitable selection of the regularization parameter, a convergence analysis may be possible.

4.5 Results

Data was simulated using equation (37), a system where detector moves on a sphere of radius $R=1/8$ was considered. The object was a $64\times 64\times 64$ volume with two balls with different intensities, and different contrasts, gray on black, white on gray, etc. One object has also a defect (crack) in some planes, see figure 3. The function was supported in a cube of side $L=1$ in the first octant with coordinate $(x_{\min},y_{\min},z_{\min})=(L/64,L/64,R)$ . Discretization parameters of data are: $N_{\alpha}=513$ , $N_{\beta}=256$ and $N_{p}=512$ . The maximal diameter of the circles generating the torus is $r^{*}_{M}=2r_{M}$ where $[R,\,r^{*}_{M}]$ is the radial support of the cube containing the phantom. Numerical integration is performed in variables $\gamma$ and $\psi$ with $\Delta\psi=2\pi/N_{\psi}$ and the variable step $\Delta\gamma=\cos^{-1}\left(p/R\right)/N_{\gamma}$ , with $N_{\gamma}=256$ and $N_{\psi}=256$ . Figure 4 shows simulated data for different values of angle $\alpha$ . The trapezoidal rule was used to perform numerical integration. According to the discretization chosen for data, the order of the spherical harmonics expansion was $N=256$ ( $N_{\alpha}=2N+1$ ). In order to get more realistic simulations, data $\mathbf{g}=\mathbf{g}_{nk}^{j}$ were also corrupted with additive Gaussian noise with zero mean. Noise variance was manually adjusted to get several signal-to-noise ratios: SNR $=10$ , $20$ and $30$ dB. These SNR correspond to $\epsilon=29$ , $10$ and $3\%$ noise levels where $\epsilon=100||\mathbf{\tilde{g}}-\mathbf{g}||_{2}/||\mathbf{g}||_{2}$ , and $\mathbf{\tilde{g}}$ is the corrupted data. In order to asses the quality of reconstruction we used the following measures of error: the Normalized Mean Square Error ( $\%$ )

\mathrm{NMSE}=\frac{100}{N^{3}}\frac{||\mathbf{f}-\tilde{\mathbf{f}}||_{2}^{2}}{\max_{i}\left\{\mathbf{f}_{i}^{2}\right\}}

(46)

and the Normalized Mean Absolute Error ( $\%$ )

\mathrm{NMAE}=\frac{100}{N^{3}}\frac{||\mathbf{f}-\tilde{\mathbf{f}}||_{1}}{\max_{i}\left\{\mathbf{f}_{i}\right\}},

(47)

Figure 3: Original 3D phantom used for simulations. Notice the crack in some planes.

where $\mathbf{f}$ is the original image and $\tilde{\mathbf{f}}$ is the reconstruction. We used the algorithm in Table 1 to carry out reconstructions for two values of the regularization parameter, $\lambda=0.01$ (noiseless data) and $\lambda=0.05$ (noisy). These values were chosen heuristically and are the same for all equations in a set (45), i.e. they do not depend on $l,m$ . Figure 5 shows reconstructions for noiseless data and figures 6 to 8 show reconstructions from corrupted data. No post-processing was applied to images. The original function as well as reconstructions are shown at planes $z=4$ , 14, 15, 22, 26, 31, 38 and 58 from top to bottom and from left to right and error metrics are displayed in figure captions. Reconstructions exhibit acceptable quality: structures inside the object are distinguished, shapes are kept and error metrics seem to be reasonable. Moreover, the crack is visible in all the images where its is expected. There are, however, background artifacts and blurring, particularly in upper planes. As expected, stronger regularization is required for noisy data and error metrics get worse with the level of noise. Although Tikhonov regularization performs well, reconstructions can be improved using more advanced regularization techniques. There are algorithms enforcing total variation minimization in the Cartesian domain [28, 29] that can be applied after step 5 in the algorithm in Table 1.

From an algorithmic point of view, the reconstruction approach allows to split up the problem in several equations reducing the size of the matrices involved with respect to a standard algebraic treatment. In this framework, computation time may be saved through parallelization since the resulting algebraic equations are independent.

Figure 5: Reconstructions from noiseless data. Errors are NMSE=0.32

\%

, NMAE=3.81

\%

. Regularization parameter

\lambda=0.01

Figure 6: Reconstructions from noisy data with SNR=30 dB (relative level of noise 3

\%

). Errors are NMSE=0.34

\%

, NMAE=3.77

\%

. Regularization parameter

\lambda=0.05

Figure 7: Reconstructions from noisy data with SNR=20 dB (relative level of noise 10

\%

). Errors are NMSE=0.40

\%

, NMAE=4.35

\%

. Regularization parameter

\lambda=0.05

Figure 8: Reconstructions from noisy data with SNR=10 dB (relative level of noise 29

\%

). Errors are NMSE=0.92

\%

, NMAE=7.38

\%

. Regularization parameter

\lambda=0.05

5 Discussion

What follow are some final comments on our work.

•

The toric transform we study exhibits rotational invariance and we employed spherical harmonic expansion to address the problem of invertibility. In order to study solution uniqueness, we used a theory for Abel and first kind Volterra equations with kernels having zeros on the diagonal [30]. The proof relied on the analysis of the kernel gradient on its diagonal. In other studies on 3D CST [19], spherical harmonic expansion was used together with an approach to solve weakly singular Volterra equations of the first kind based on the theory in [31].
•

Numerical simulations confirm an additional advantage of this modality: the scanning is feasible when the system is smaller than the object under study. In the proposed simulation we use an object whose side was eight times larger than the radius of the detection sphere.
•

Since the model neglects attenuation, extra analysis should be carried out when addressing the realistic scenario. When attenuation is significantly low, few photon counts are expected leading to a sparse signal. In such case, the algorithm may be adapted to deal with sparse data. The influence of attenuation is modeled through weighted Radon transforms on tori leading to weighted kernels. Under some assumptions, this mathematical framework may be applied with weighted kernels both in the theoretical treatment as well as in reconstructions. Other works address the reconstruction problem when attenuation is present [18, 32].

6 Conclusion

We studied a Compton scattering tomography in three dimensions suitable for scanning large objects. This system has interesting advantages such as fixed and unique source and compactness. Measured data was modeled by a toric Radon transform. We demonstrated the invertibility of the toric transform by proving the uniqueness of the solutions of the one dimensional Abel’s type equation resulting from its spherical harmonics expansion. In addition, we developed a reconstruction method based on the discrete spherical harmonics expansion and Tikhonov regularization. Numerical simulations confirm the feasibility of the approach. Future challenges emerge such as incorporating attenuation to the model and dealing with sparse data. Some studies on that direction are on the way.

Acknowledgements

J. Cebeiro is supported a CONICET postdoctoral grant ( $\#$ 171800). D. Rubio, J. Cebeiro and M. A. Morvidone are partially supported by SOARD-AFOSR (grant number FA9550-18-1-0523). C. Tarpau research work is supported by grants from Région Ile-de-France (in Mathematics and Innovation) 2018-2021 and LabEx MME-DII (Modèles Mathématiques et Economiques de la Dynamique, de l’Incertitude et des Interactions) (No. ANR-11-LBX-0023-01).

Appendix A Spherical Harmonics Expansion

Here, we present the proof lemma 2.1 explaining how to get form (3) to (7). First, we define some useful rotation properties.

A.1 Rotation properties

Given $h^{-1}=u(\alpha)a(\beta)$ as introduced in (2) we define the rotation operator $\Lambda_{h}$ by its action on a function $f$ as:

\left(\Lambda_{h}f\right)\left(r\Theta(\gamma,\psi)\right)=f\left(rh^{-1}\Theta(\gamma,\psi)\right).

(48)

Applying this linear operator $\Lambda_{h}$ to a function expanded as in (5) we have

(\Lambda_{h}f)(r\Theta(\gamma,\psi))=\sum_{l\in\mathbb{N}}\sum_{|m|\leq l}f_{lm}(r)(\Lambda_{h}Y_{l}^{m})(\gamma,\psi).

(49)

Properties $(i)$ - $(iii)$ will be useful in what follows:

1.

Any rotated spherical harmonic of degree $l$ is a linear combination of $Y^{n}_{l}$ with $|n|\leq l:$

$(\Lambda_{h}Y_{l}^{m})(\gamma,\psi)=\sum_{|n|\leq l}Y_{l}^{n}(\gamma,\psi)D_{n,m}^{(l)}(h),$ (50)
2.

Matrices $D^{(l)}(h)$ verify:

$D_{0,m}^{(l)}(h)=\overline{D_{m,0}^{(l)}}(h^{-1}).$ (51)
3.

From the definition of the matrix entries $D_{n,m}^{(l)}(h)$ [22], it follows that:

$Y_{l}^{m}(\alpha,\beta)=\sqrt{\frac{2l+1}{4\pi}}\overline{D^{(l)}_{m,0}}(h^{-1}).$ (52)

See [22] and [20] for details.

A.2 Proof of lemma 2.1

Proof.

From (3) and (48) we have:

\mathcal{R}_{\mathcal{T}}f(\alpha,\beta,\omega)=\int_{0}^{2\omega-\pi}d\gamma\,\int_{0}^{2\pi}d\psi\,r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,\left(\Lambda_{h}f\right)\left(r^{\omega}(\gamma)\Theta(\gamma,\psi\right)).

(53)

Using (49) and (50), we can write:

	$\displaystyle\mathcal{R}_{\mathcal{T}}f(\alpha,\beta,\omega)={\sum_{l\in\mathbb{N}}\sum_{\|m\|\leq l}}\int_{0}^{2\omega-\pi}d\gamma\,\int_{0}^{2\pi}d\psi\,r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,f_{lm}\left(r^{\omega}(\gamma)\right)(\Lambda_{h}Y_{l}^{m})\left(\gamma,\psi\right)$
	$\displaystyle={\sum_{l\in\mathbb{N}}\sum_{\|m\|\leq l}\sum_{\|n\|\leq l}}\int_{0}^{2\omega-\pi}d\gamma\,\int_{0}^{2\pi}d\psi\,r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,f_{lm}\left(r^{\omega}(\gamma)\right)Y_{l}^{n}\left(\gamma,\psi\right)D^{(l)}_{n,m}(h).$

Now, using (4) and performing the integration on variable $\psi$ , we have that the only non-vanishing term is the $n=0$ term, so:

	$\displaystyle\mathcal{R}_{\mathcal{T}}f(\alpha,\beta,\omega)={\sum_{l\in\mathbb{N}}\sum_{\|m\|\leq l}\sum_{\|n\|\leq l}}$
	$\displaystyle\int_{0}^{2\omega-\pi}d\gamma\,\int_{0}^{2\pi}d\psi\,r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,f_{lm}\left(r^{\omega}(\gamma)\right)(-1)^{n}\sqrt{\frac{(2l+1)(l-n)!}{4\pi(l+n)!}}P^{n}_{l}(\cos\gamma)e^{in\psi}D^{(l)}_{n,m}(h)$
	$\displaystyle={\sum_{l\in\mathbb{N}}\sum_{\|m\|\leq l}}2\pi\int_{0}^{2\omega-\pi}d\gamma r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,f_{lm}\left(r^{\omega}(\gamma)\right)P^{0}_{l}(\cos\gamma)\sqrt{\frac{(2l+1)}{4\pi}}D_{0,m}^{(l)}(h)$
	$\displaystyle={\sum_{l\in\mathbb{N}}\sum_{\|m\|\leq l}}\sqrt{\frac{(2l+1)}{4\pi}}D_{0,m}^{(l)}(h)2\pi\int_{0}^{2\omega-\pi}d\gamma r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,f_{lm}\left(r^{\omega}(\gamma)\right)P^{0}_{l}(\cos\gamma).$		(55)

Then, using properties (51) and (52) we arrive to

\mathcal{R}_{\mathcal{T}}f(\alpha,\beta,\omega)={\sum_{l\in\mathbb{N}}\sum_{|m|\leq l}}Y_{l}^{m}(\alpha,\beta)\,\,\,\,2\pi\int_{0}^{2\omega-\pi}d\gamma r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,f_{lm}\left(r^{\omega}(\gamma)\right)P^{0}_{l}(\cos\gamma).

(56)

Finally, the factor accompanying the $l,m$ -spherical harmonic is the $l,m$ -coefficient of the Fourier expansion of $(\mathcal{R}_{\mathcal{T}}f)$ in terms of the scattering angle $\omega$ and we get (7). This completes the proof.

∎

Appendix B An Algorithm for Discrete Spherical Harmonics Expansion

Here, we describe succinctly our implementation of the forward-inverse Discrete Spherical Harmonics Transform (DSHT-IDSHT). The technique is based on the algorithms for the discrete Fourier and Legendre transforms as described in [33, 34, 20]. Let’s be the discrete function $\textbf{F}_{kn}^{j}=F(r_{j},\theta_{k},\varphi_{n})$ sampled with the following parameters $\Delta r=(r_{\max}-r_{\min})/(N_{r}-1)$ and $\Delta\varphi=2\pi/(2N+1)$ with their corresponding indices $j=0,...,N_{r}-1$ and $n=-N,...,N$ . Discrete values of variables for radial and azimuthal variables are $r_{j}=j\Delta r+r_{\min}$ and $\varphi_{n}=2\pi n/(2N+1)$ . Latitude variable $\theta$ can be arbitrarily sampled in the interval $[0,\pi]$ and is labeled $\theta_{k}$ , $k=1,...,N_{\theta}$ , we used uniform sampling. The relationship between the discrete function $\textbf{F}_{kn}^{j}$ and its representation in the domain of discrete spherical harmonics $\mathbf{F}_{lm}^{j}$ is summarized in the diagram:

\mathbf{F}_{nk}^{j}\mathrel{\mathop{\rightleftarrows}^{\mathrm{DSHT}}_{\mathrm{IDSHT}}}\mathbf{F}_{lm}^{j}\,\,\,\,\,\,:\,\,\,\,\,\,\left\{\textbf{F}_{kn}^{j}\right\}\mathrel{\mathop{\rightleftarrows}^{{\textstyle\frac{1}{N}}\cdot\mathrm{DFT}}_{N\cdot\mathrm{IDFT}}}\left\{{\textbf{ F}}_{km}^{j}\right\}\mathrel{\mathop{\rightleftarrows}^{\mathrm{IDLT}}_{\mathrm{DLT}}}\left\{\textbf{F}_{lm}^{j}\right\},

while the first block is given by the well known discrete Fourier pairs (DFT-IDFT), the second one is given by the associated Discrete Legendre Transform pairs (DLT-IDLT)

\textbf{F}_{km}^{j}=\mathrm{DLT}^{m}\left\{\textbf{F}_{lm}^{j}\right\}=\sum_{l=|m|}^{N}\textbf{F}_{lm}^{j}q_{l}^{m}P_{l}^{m}(t_{k}),

(57)

and

\textbf{F}_{lm}^{j}=\mathrm{IDLT}^{m}\left\{\textbf{F}_{km}^{j}\right\}=\sum_{k=1}^{N_{\theta}}\textbf{F}_{km}^{j}q_{l}^{m}P_{l}^{m}(t_{k})w_{k},

(58)

where $m=-N,...,N$ , $l=|m|,...,N$ , $t_{k}=\cos\theta_{k}$ , $w_{k}$ are the Gaussian quadrature coefficients and

q_{l}^{m}=(-1)^{m}\sqrt{{\textstyle\frac{2l+1}{4\pi}}{\textstyle\frac{(l-m)!}{(l+m)!}}}.

(59)

An alternative implementation can be found in [20]. A preliminary version of the algorithm can be found in the code repository [35].

References

[1] Norton S J 1994 Compton scattering tomography J. Appl. Phys 76 2007–2015
[2] Cebeiro J, Nguyen M K, Morvidone M A and Noumowé A 2017 New ’improved’ Compton scatter tomography modality for investigative imaging of one-sided large objects Inverse Problems in Science and Engineering 25 1676–1696
[3] Webber J and Miller E L 2020 Compton scattering tomography in translational geometries Inverse Problems 36 025007
[4] Tarpau C, Cebeiro J and Nguyen M K 2019 A new bi-imaging NDT system for simultaneous recovery of attenuation and electronic density maps Eleventh Int. Conf. Non Destruct. Test. Aerosp
[5] Jones K C, Redler G, Templeton A, Bernard D, Turian J V and Chu J C 2018 Characterization of Compton-scatter imaging with an analytical simulation method Physics in Medicine & Biology 63 025016
[6] Redler G, Jones K C, Templeton A, Bernard D, Turian J and Chu J C 2018 Compton scatter imaging: A promising modality for image guidance in lung stereotactic body radiation therapy Medical physics 45 1233–1240
[7] Cormack A M 1981 The Radon transform on a family of curves in the plane Proceedings of the American Mathematical Society 83 325 – 330
[8] Cormack A M 1982 The Radon transform on a family of curves in the plane Proceedings of the American Mathematical Society 86 293– 298
[9] Cormack A M 1984 Radon’s problem - old and new SIAM - AMS Proceedings 14
[10] Truong T T and Nguyen M K 2011 Radon transforms on generalized Cormack’s curves and a new Compton scatter tomography modality Inverse Problems 27
[11] Truong T T and Nguyen M K 2011 On new V-line transforms in ${R}^{2}$ and their inversion Journal of Physics A: Mathematical and Theoretical 44 13
[12] Truong T T and Nguyen M K 2012 Recent Developments on Compton Scatter Tomography: Theory and Numerical Simulations, Chapter 6, pp 102-128, in Numerical Simulation - From Theory to Industry (Edited by Mykhaylo Andriychuk, INTECH. ISBN: 978-953-51-0749-1, DOI: 10.5772/50012)
[13] Webber J 2016 X-ray Compton scattering tomography Inverse Problems in Science and Engineering 24 1323–1346
[14] Norton S J 2019 Compton-scattering tomography with one source and one detector: a simple derivation of the filtered-backprojection solution Inverse Problems in Science and Engineering 0 1–12
[15] Tarpau C, Cebeiro J, Morvidone M A and Nguyen M K 2020 A new concept of compton scattering tomography and the development of the corresponding circular radon transform IEEE Transactions on Radiation and Plasma Medical Sciences 4 433–440
[16] Webber J W and Quinto E T 2020 Microlocal analysis of a Compton tomography problem SIAM Journal on Imaging Sciences 13 746–774
[17] Tarpau C, Cebeiro J, Nguyen M K, Rollet G and Morvidone M A Accepted may 2020 Analytic inversion of a Radon transform on double circular arcs with applications in Compton scattering tomography IEEE Transactions on Computational Imaging
[18] Rigaud G and Hahn B 2018 3D Compton scattering imaging and contour reconstruction for a class of Radon transforms Inverse Problems 34
[19] Webber J W and Lionheart W R B 2018 Three dimensional Compton scattering tomography Inverse Problems 34
[20] Driscoll J and Healy D 1994 Computing Fourier transforms and convolutions on the 2-sphere Advances in Applied Mathematics 15 202–250
[21] Quinto E T 1983 The invertibility of rotation invariant Radon transforms Journal of Mathematical Analysis and Applications 91 510 – 522 ISSN 0022-247X
[22] Biedenharn L C and Louck J D 1981 Angular Momentum in Quantum Physics (Reading, MA, USA: Addison Wesley Publishing Company)
[23] Laden H N 1938 An historical and critical development of the theory of Legendre polynomials before 1900 Ph.D. thesis University of Maryland
[24] Schiefeneder D and Haltmeier M 2017 The Radon transform over cones with vertices on the sphere and orthogonal axes SIAM Journal on Applied Mathematics 77 1335–1351
[25] Jackson J D 1999 Classical Electrodynamics, Third Edition (Hoboken, NJ, USA: John Wiley and Sons, Inc.)
[26] Haltmeier M, Moon S and Schiefeneder D 2017 Inversion of the attenuated V-line transform with vertices on the circle IEEE Transactions on Computational Imaging 3 853–863
[27] Weiss R and Anderssen R A 1971 A product integration method for a class of singular first kind volterra equations Numer. Math. 18 442–456
[28] Condat L 2013 A direct algorithm for 1-d total variation denoising IEEE Signal Processing Letters 20 1054–1057
[29] Condat L 2014 A generic proximal algorithm for convex optimization—application to total variation minimization IEEE Signal Processing Letters 21 985–989
[30] Moon S and Haltmeier M 2017 Analytic inversion of a conical Radon transform arising in application of Compton cameras on the cylinder SIAM Journal on Imaging Sciences 10 535–557
[31] Bownds J M 1976 On solving weakly singular Volterra equations of the first kind with Galerkin approximations Mathematics of Computation 30 747–757
[32] Rigaud G and Hahn B N 2020 Reconstruction algorithm for 3d compton scattering imaging with incomplete data Inverse Problems in Science and Engineering 0 1–23
[33] Basko R, Zeng G L and Gullberg G T 1998 Application of spherical harmonics to image reconstruction for the Compton camera Physics in Medicine and Biology 43 887–94
[34] Taguchi K, Zeng G L and Gullberg G T 2001 Cone-beam image reconstruction using spherical harmonics Physics in Medicine and Biology 46 127–138
[35] Cebeiro J, Tarpau C, Morvidone M A, Rubio D and Nguyen M K 2020 Algorithm for Discrete Spherical Harmonic Expansion of 3D functions URL https://github.com/Javiercbr/Discrete-Spherical-Harmonics

	$\displaystyle\mathcal{R}_{\mathcal{T}}f(\alpha,\beta,\omega)={\sum_{l\in\mathbb{N}}\sum_{\|m\|\leq l}\sum_{\|n\|\leq l}}$
	$\displaystyle\int_{0}^{2\omega-\pi}d\gamma\,\int_{0}^{2\pi}d\psi\,r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,f_{lm}\left(r^{\omega}(\gamma)\right)(-1)^{n}\sqrt{\frac{(2l+1)(l-n)!}{4\pi(l+n)!}}P^{n}_{l}(\cos\gamma)e^{in\psi}D^{(l)}_{n,m}(h)$
	$\displaystyle={\sum_{l\in\mathbb{N}}\sum_{\|m\|\leq l}}2\pi\int_{0}^{2\omega-\pi}d\gamma r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,f_{lm}\left(r^{\omega}(\gamma)\right)P^{0}_{l}(\cos\gamma)\sqrt{\frac{(2l+1)}{4\pi}}D_{0,m}^{(l)}(h)$
	$\displaystyle={\sum_{l\in\mathbb{N}}\sum_{\|m\|\leq l}}\sqrt{\frac{(2l+1)}{4\pi}}D_{0,m}^{(l)}(h)2\pi\int_{0}^{2\omega-\pi}d\gamma r^{\omega}(\gamma)\frac{\sin\gamma}{\sin\omega}\,f_{lm}\left(r^{\omega}(\gamma)\right)P^{0}_{l}(\cos\gamma).$		(55)

On a three-dimensional Compton scattering tomography system with fixed source

Abstract

1 Introduction

2 A setup for three dimensional Compton scattering tomography

2.1 Spherical harmonics expansion

Lemma 2.1.

3 Uniqueness of the solution of the toric Radon transform

3.1 An alternative expression for the integral equation

3.2 On the invertibility of ℛ𝒯​f\mathcal{R}_{\mathcal{T}}f in terms of spherical harmonics

Lemma 3.1 ([24], Theorem 3.4).

Theorem 3.1 (Invertibility of equation (9)).

Proof.

Corollary 3.1 (Invertibility of the ℛ𝒯​f\mathcal{R}_{\mathcal{T}}f).

Proof.

4 Numerical Simulations

4.1 Alternative forward model

4.2 The algebraic problem

4.3 Matrix generation

4.4 Overview of the reconstruction algorithm

4.5 Results

5 Discussion

6 Conclusion

Acknowledgements

Appendix A Spherical Harmonics Expansion

A.1 Rotation properties

A.2 Proof of lemma 2.1

Proof.

Appendix B An Algorithm for Discrete Spherical Harmonics Expansion

References

References

3.2 On the invertibility of $\mathcal{R}_{\mathcal{T}}f$ in terms of spherical harmonics

Corollary 3.1 (Invertibility of the $\mathcal{R}_{\mathcal{T}}f$ ).