Image reconstruction from radially incomplete spherical Radon data

TL;DR

This paper presents a new method for uniquely reconstructing images from incomplete spherical Radon data, with practical algorithms demonstrated through numerical examples in medical and radar imaging.

Contribution

It introduces a novel inversion formula for radially incomplete spherical Radon data and provides a robust computational algorithm for practical image reconstruction.

Findings

Unique reconstruction in spherical shells with incomplete data

Reconstruction formulas applicable inside, outside, and across the data sphere

Algorithm demonstrates high accuracy and efficiency in numerical tests

Abstract

We study inversion of the spherical Radon transform with centers on a sphere (the data acquisition set). Such inversions are essential in various image reconstruction problems arising in medical, radar and sonar imaging. In the case of radially incomplete data, we show that the spherical Radon transform can be uniquely inverted recovering the image function in spherical shells. Our result is valid when the support of the image function is inside the data acquisition sphere, outside that sphere, as well as on both sides of the sphere. Furthermore, in addition to the uniqueness result our method of proof provides reconstruction formulas for all those cases. We present a robust computational algorithm based on our inversion formula and demonstrate its accuracy and efficiency on several numerical examples.

Equations

$$x_1=r\cos {\phi }_1$$

$$x_2=r\sin {\phi }_1\cos {\phi }_2$$

$$x_3=r\sin {\phi }_1\sin {\phi }_2\cos {\phi }_3$$

$$⋮$$

$$x_{n-1}=r\sin {\phi }_1\sin {\phi }_2\cdots \sin {\phi }_{n-2}\cos {\phi }_{n-1}$$

$$x_n=r\sin {\phi }_1\sin {\phi }_2\cdots \sin {\phi }_{n-2}\sin {\phi }_{n-1},$$

$$(\cos {\alpha }_1,\sin {\alpha }_1\cos {\alpha }_2,\cdots ,\sin {\alpha }_1\sin {\alpha }_2\cdots \sin {\alpha }_{n-2}\cos {\alpha }_{n-1},\sin {\alpha }_1\sin {\alpha }_2\cdots \sin {\alpha }_{n-2}\cos {\alpha }_{n-1}).$$

$$Rf(\rho ,\alpha )=g(\rho ,\alpha )={\int }_{S(\rho ,\alpha )}f\mathrm{d}\Omega ,$$

$$A(R_1,R_2)=\{x\in {\mathbb{R}}^n:R_1<|x|<R_2\}.$$

$$A(R_1,R_2)=\{(r,\phi ):R_1<r<R_2,0\le {\phi }_i\le \pi \text{mbox}for1\le i\le n-2\text{mbox}and0\le {\phi }_{n-1}\le 2\pi \}.$$

$$f(r,\phi )=\sum _{l=0}^{\infty }{f}_l(r)Y_l(\phi )$$

$$g(\rho ,\alpha )=\sum _{l=0}^{\infty }{g}_l(\rho )Y_l(\alpha ).$$

$$g(\rho ,\alpha )=\underset{S(\rho ,\alpha )}{\int }f\mathrm{d}\Omega =\underset{S(\rho ,\alpha )}{\int }\sum _{l=0}^{\infty }{f}_l(r)Y_l(\phi )\mathrm{d}\Omega .$$

$$g(\rho ,\alpha )=\sum _{l=0}^{\infty }\underset{S(\rho ,\alpha )}{\int }{f}_l(r)Y_l(\phi )\mathrm{d}\Omega .$$

$$d\Omega ={\rho }^{n-1}{\sin }^{n-2}{\tilde{\phi }}_1{\sin }^{n-3}{\tilde{\phi }}_2\cdots {\sin }^2{\tilde{\phi }}_{n-3}\sin {\tilde{\phi }}_{n-2}\mathrm{d}{\tilde{\phi }}_1\cdots \mathrm{d}{\tilde{\phi }}_{n-1}.$$

$$A_i=\vec{P}\cdot {\vec{C}}_i\text{mbox}for1\le i\le n.$$

$$\det \left(\begin{array}{c}{\nabla }_{\phi }{A}_1\\ ⋮\\ {\nabla }_{\phi }{A}_{n-1}\end{array}\right)=\left(\vec{P}\cdot {\vec{C}}_n\right){\sin }^{n-2}{\phi }_1{\sin }^{n-3}{\phi }_2\cdots {\sin }^2{\phi }_{n-3}\sin {\phi }_{n-2}\mathrm{d}{\phi }_1\cdots \mathrm{d}{\phi }_{n-1}.$$

$$\det (AB)=\sum _J\det (A(J))\det (B(J))$$

$$J=j_1,j_2,\cdots ,j_m,1\le j_1<j_2\cdots <j_k\le m$$

$$\left(\begin{array}{c}{\nabla }_{\phi }{A}_1\\ ⋮\\ {\nabla }_{\phi }{A}_{n-1}\end{array}\right)=\left(\begin{array}{c}{\vec{C}}_1\\ ⋮\\ {\vec{C}}_{n-1}\end{array}\right)\left(\begin{array}{c}\frac{\partial {\vec{P}}^t}{\partial {\phi }_1}\cdots \frac{\partial {\vec{P}}^t}{\partial {\phi }_{n-1}}\end{array}\right).$$

$$\left(\begin{array}{c}\frac{\partial {\vec{P}}^t}{\partial {\phi }_1}\cdots \frac{\partial {\vec{P}}^t}{\partial {\phi }_{n-1}}\end{array}\right)=$$

$$\left(\begin{array}{ccccc}-\sin {\phi }_1& 0& \cdots & \cdots & 0\\ \cos {\phi }_1\cos {\phi }_2& -\sin {\phi }_1\sin {\phi }_2& \cdots & \cdots & 0\\ \cos {\phi }_1\sin {\phi }_2\cos {\phi }_3& \sin {\phi }_1\cos {\phi }_2\cos {\phi }_3& -\sin {\phi }_1\sin {\phi }_2\sin {\phi }_3& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & 0\\ \cos {\phi }_1\sin {\phi }_2\cdots \sin {\phi }_{n-2}\cos {\phi }_{n-1}& \cdots & \cdots & \cdots & -\sin {\phi }_1\cdots \sin {\phi }_{n-1}\\ \cos {\phi }_1\sin {\phi }_2\cdots \sin {\phi }_{n-1}& \cdots & \cdots & \cdots & \sin {\phi }_1\cdots \sin {\phi }_{n-2}\cos {\phi }_{n-1}\end{array}\right).$$

$$={\sin }^{n-2}{\phi }_1{\sin }^{n-3}{\phi }_2\cdots \sin {\phi }_{n-2}\det \left(\begin{array}{cccc}{v}_1& v_2& \cdots & v_{n-1}\end{array}\right)$$

$$\det \left(\begin{array}{cccc}{v}_1& \cdots & v_{n-1}& \vec{P}\end{array}\right)=\pm 1,$$

$$\sum _{i=1}^n(-1{)}^{n+i}{\vec{P}}_i\cdot M_{in}=\pm 1,$$

$$\left((-1{)}^{1+n}{M}_{1n},\cdots ,(-1{)}^{2n}{M}_{nn}\right)=\pm \vec{P}.$$

$$\det \left(\begin{array}{c}{\vec{C}}_1\\ ⋮\\ {\vec{C}}_{n-1}\\ {\vec{C}}_n\end{array}\right)=1.$$

$$\left((-1{)}^{1+n}{\tilde{M}}_{n1},\cdots ,(-1{)}^{2n}{\tilde{M}}_{nn}\right)=\pm {\vec{C}}_n.$$

$$\mathrm{d}\Omega =\frac{{\rho }^{n-2}{r}^2}{\left|r-A_1\right|}{\sin }^{n-2}{\phi }_1\cdots \sin {\phi }_{n-2}\mathrm{d}{\phi }_1\cdots \mathrm{d}{\phi }_{n-1},$$

$$d\Omega ={\rho }^{n-1}{\sin }^{n-2}{\tilde{\phi }}_1{\sin }^{n-3}{\tilde{\phi }}_2\cdots {\sin }^2{\tilde{\phi }}_{n-3}\sin {\tilde{\phi }}_{n-2}\mathrm{d}{\tilde{\phi }}_1\cdots \mathrm{d}{\tilde{\phi }}_{n-1}.$$

$$\cos {\tilde{\phi }}_1=\frac{\left(r{\vec{P}}_1-{\vec{C}}_1\right)\cdot {\vec{C}}_1}{|r\vec{P}-{\vec{C}}_1|}=\frac{r{\vec{P}}_1\cdot {\vec{C}}_1-1}{\rho }=\frac{r{A}_1-1}{\rho }.$$

$$\left(\begin{array}{c}\cos {\tilde{\phi }}_1\\ ⋮\\ \cos {\tilde{\phi }}_{n-1}\end{array}\right)=\left(\begin{array}{c}\frac{r{A}_1-1}{\rho }\\ \frac{A_2}{\sqrt{1-A_1^2}}\\ ⋮\\ \frac{A_{n-1}}{\sqrt{1-\left(A_1^2+\cdots +A_{n-2}^2\right)}}\end{array}\right).$$

$$({\phi }_1,\cdots ,{\phi }_{n-1})\to (\cos {\tilde{\phi }}_1,\cdots ,\cos {\tilde{\phi }}_{n-1}).$$

Full text

1. Introduction

The spherical Radon transform (SRT) maps a function of $n$ variables to its integrals over a family of spheres in $\mathbb{R}^{n}$. Such transforms naturally appear in mathematical models of various imaging modalities in medicine [10, 24, 27, 33, 35, 36, 40, 43, 44, 48], geophysical applications [16, 32], radar [15], as well as in some purely mathematical problems of approximation theory [1, 3, 30], PDEs [1, 2, 17, 18, 19, 25, 28, 29] and integral geometry [5, 6, 7, 8, 11, 12, 20, 22, 34, 41, 42].

One of the most important questions related to SRT is the possibility of its stable inversion. Since the family of all spheres in $\mathbb{R}^{n}$ has $n+1$ dimensions, the problem of inversion from the set of integrals along all spheres is overdetermined. Hence it is customary to consider the problem of inverting the SRT from the restriction of the full set of integrals to an $n$-dimensional subset. While one can come up with several different choices of such subsets, a common approach (especially in imaging applications) is to restrict the centers of integration spheres to a hypersurface in $\mathbb{R}^{n}$.

For example, a simple model of thermoacoustic tomography (TAT) can be described as follows. A biological object under investigation is irradiated with a short pulse of electromagnetic waves. Certain part of that radiation gets absorbed in the body heating up the tissue leading to its thermoelastic expansion. The latter generates ultrasound waves, which propagate through the body and are registered by transducers placed on its surface. Under a simplifying assumption of constant speed $c$ of ultrasound waves in the tissue, at any moment of time $t$, a single transducer records a superposition of signals generated at locations that are at the fixed distance $ct$ from the transducer. In other words, the transducer measurements can be modeled as integrals of a function along spheres centered at the transducer location and of different radii (depending on time). By moving the transducer around the surface of the object (or equivalently using an array of such transducers) one can essentially measure a 3-dimensional family of spherical integrals of the unknown image function. Hence to recover the image in this simple TAT model, one would need to invert the SRT in the setup described above. Similar mathematical problems arise also in various models of ultrasound reflection tomography, as well as in sonar and radar imaging.

While our work is motivated by its potential applications in imaging problems, we study the spherical Radon transform in $\mathbb{R}^{n}$ for any $n\geq 3$. We discuss the inversion of SRT from integrals of a function $f$ along spheres whose centers lie on the surface of the unit (data acquisition) sphere111Our results carry over with little difficulty when the centers of the SRT data lie on a sphere of radius $R$.. With the additional restriction on the set of radii of integration spheres, we prove the uniqueness as well as derive reconstruction formulas for $f$ from such data. We provide several results that hold for the cases when the support of $f$ is inside, outside, or on both sides of the unit sphere. More precisely, for the case when the support of a function $f$ is inside the unit sphere, our result shows that in order to reconstruct $f$ in the spherical shell $\{r<|x|<1\}$ for any $r<1$, we only need SRT data with centers on the unit sphere and for all radii $\rho$ such that $0<\rho<1-r$. Analogous statements can be made for the case when the support of $f$ is outside or on both sides of the unit sphere. In connection with this, we mention the result [18, Theorem 5], where it was shown that for a bounded open connected set $D$ in $\mathbb{R}^{n}$ for $n$ odd, a function $f$ supported in $\overline{D}$ can be reconstructed from SRT data with centers on $\partial D$ and all radii $\rho$ such that $\rho\in[0,\mbox{diam}(D)/2]$. One of the consequences of our work is a generalization of this result for the case of even dimensions, as well as when the support of the function lies inside, outside or on both sides for the case when $D$ is a sphere. We emphasize here that the uniqueness result [17, Theorem 5] was already generalized for variable sound speeds in [43] in all space dimensions; see Prop. 2 in that paper. If one is interested in uniqueness results alone, unique continuation arguments as in [43] or analytic microlocal analysis methods as in [4, 39] could be used222We thank Plamen Stefanov for bringing this as well as the result stated in the previous sentence to our attention., although to the best of our knowledge, even for the case of spherical acquisition surface and for functions supported outside or on both sides of the sphere, such results have not been published. The advantage of our work, in the specific setting where the acquisition geometry is the unit sphere, is that it provides in addition to uniqueness results, inversion formulas using radially partial data.

The paper is organized as follows. The main results are stated in Section 2 and the proofs are presented in Section 3. In Section 4 we write down the inversion formulas for the special case of $n=3$. In Section 5 we discuss the numerical algorithm based on the product integration method. In Section 6, we provide numerical examples illustrating the accuracy and efficiency of the proposed inversion algorithms.

2. Main results

We consider the usual spherical coordinate system:

[TABLE]

where $0\leq\varphi_{i}\leq\pi$ for $1\leq i\leq n-2$ and $0\leq\varphi_{n-1}\leq 2\pi$. For simplicity, from now on, we will denote $\varphi=(\varphi_{1},\cdots,\varphi_{n-1})$. Let us consider the unit sphere centered at the origin in $\mathbb{R}^{n}$ and fix an arbitrary point $C$ on this sphere. We will denote $C$ in the above spherical coordinates by $\alpha$, where $\alpha=(\alpha_{1},\cdots,\alpha_{n-1})$. The Cartesian coordinates of the point $C$ will then be

[TABLE]

Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be a continuous function of compact support. Consider a sphere $S(\rho,\alpha)$ of radius $\rho$ centered at $C$. The spherical Radon transform of $f$ along the sphere $S(\rho,\alpha)$ for $\rho>0$ and $\alpha=(\alpha_{1},\cdots,\alpha_{n-2},\alpha_{n-1})\in[0,\pi]\times\cdots\times[0,\pi]\times[0,2\pi]$ is defined as

[TABLE]

where $\mathrm{d}\Omega$ is the usual surface measure on the sphere $S(\rho,\alpha)$.

Finally, let us denote by $A(R_{1},R_{2})$ the spherical shell lying between the spheres of radii $R_{1}$ and $R_{2}$ centered at the origin. Expressing in Cartesian coordinates:

[TABLE]

Note that this spherical shell can also be expressed in spherical coordinates by

[TABLE]

We now state the main results.

Theorem 2.1** (Exterior support).**

Let $f(r,\varphi)$ be a $C^{\infty}$ function supported inside $A(1,\,3)$. If $Rf(\rho,\alpha)$ is known for all $(\rho,\alpha)$ with $0<\rho<R_{1}$ where $0<R_{1}<2$ and $\alpha\in[0,\pi]\times\cdots\times[0,\pi]\times[0,2\pi]$, then $f(r,\varphi)$ can be uniquely recovered in the spherical shell $A(1,\,1+R_{1})$ with an iterative reconstruction procedure.

Theorem 2.2** (Interior support).**

Let $f(r,\varphi)$ be a $C^{\infty}$ function supported inside $A(\varepsilon,1)$. If $Rf(\rho,\alpha)$ is known for all $(\rho,\alpha)$ with $0<\rho<1-\varepsilon$, where $0<\varepsilon<1$ and $\alpha\in[0,\pi]\times\cdots\times[0,\pi]\times[0,2\pi]$, then $f(r,\varphi)$ can be uniquely recovered in the spherical shell $A(\varepsilon,1)$ with an iterative reconstruction procedure.

Theorem 2.3** (Interior and exterior support).**

Let $f(r,\varphi)$ be a $C^{\infty}$ function supported inside the ball $B(0,R_{2})$ centered at the origin and of radius $R_{2}>2$. Define $R_{1}=R_{2}-2$. If $Rf(\rho,\alpha)$ is known for all $\rho$ with $R_{2}-1<\rho<R_{2}+1$ and $\alpha\in[0,\pi]\times\cdots\times[0,\pi]\times[0,2\pi]$, then $f(r,\varphi)$ can be uniquely recovered in the spherical shell $A(R_{1},R_{2})$ with an iterative reconstruction procedure.

3. Proofs

Let $\{Y_{l}\}$ be the full set of spherical harmonics forming an orthonormal basis for $L^{2}$ functions on $\mathbb{S}^{n-1}$. We expand $f$ and $g$ into a series involving $\{Y_{l}\}$. We have

[TABLE]

Due to rotational invariance of the spherical Radon transform, the spherical harmonics expansion of $f$ and $g$ leads to diagonalization of the transform, that is, for each fixed $l\geq 0$ the coefficient $g_{l}(\rho)$ depends only on $f_{l}(r)$. Our main goal in the following calculations is to find that relationship, and express $f_{l}(r)$ through $g_{l}(\rho)$.

Using (2) in (1), the spherical Radon transform is expressed as

[TABLE]

Since $f$ is a $C^{\infty}$ function of compact support, by straightforward modifications of the arguments in [26]333The result in [26] shows that the spherical harmonics series of a sufficiently smooth function $h$ on the unit sphere converges uniformly to $h$. One can adapt the same arguments to show that the spherical harmonics series of a compactly supported smooth function $f$ on $(0,\infty)\times\mathbb{S}^{n-1}$ converges uniformly to $f$ in the radial and angular variables., we have that the spherical harmonics series of $f$ converges uniformly to $f$. Hence we can interchange the sum and the integral, and we have

[TABLE]

We denote by $\vec{C}_{1}$ the vector pointing from the origin to the fixed point $C$ on the unit sphere in $\mathbb{R}^{n}$. Let us fix an orthonormal coordinate system for the plane $\vec{C}_{1}^{\perp}$, which we denote by $\vec{C}_{2},\cdots,\vec{C}_{n}$. Reordering the vectors if necessary, we assume that $(\vec{C}_{1},\cdots,\vec{C}_{n-1},\vec{C}_{n})$ is an oriented orthonormal coordinate system for $\mathbb{R}^{n}$. We can consider spherical coordinates with respect to this coordinate system, and denote them by $(\widetilde{r},\widetilde{\varphi})$ where $\widetilde{\varphi}=(\widetilde{\varphi}_{1},\cdots,\widetilde{\varphi}_{n-1})$.

The surface measure $\mathrm{d}\Omega$ on the sphere $S(\rho,\alpha)$ in this coordinate system is

[TABLE]

Our next goal, which we state as Proposition 3.3 below, is to find the above surface measure $\mathrm{d}\Omega$ in the spherical coordinate system $(r,\varphi)$. To this end, we first prove a lemma.

Let us consider an arbitrary point $P$ in the coordinate system $(r,\varphi)$ and denote it as $r\vec{P}$ with $\vec{P}\in\mathbb{S}^{n-1}$. We recall the orthonormal coordinate system $(\vec{C}_{1},\cdots,\vec{C}_{n-1},\vec{C}_{n})$ with respect to the arbitrary fixed point $C\in\mathbb{S}^{n-1}$ introduced above and define

[TABLE]

Lemma 3.1**.**

We have the following formula:

[TABLE]

The proof of the above lemma relies on the following result due to Cauchy and Binet.

Theorem 3.2** (Cauchy-Binet).**

Let $A$ be an $m\times n$ and $B$ be an $n\times m$ matrix. Then

[TABLE]

with

[TABLE]

and $A(J)$ denotes the matrix formed from $A$ with the columns $J$ with the order preserved and $B(J)$ denotes the matrix formed from $B$ with the rows $J$ with the order preserved.

Proof of Lemma 3.1.

Since $A_{i}=\vec{P}\cdot\vec{C}_{i}$, we have that

[TABLE]

We are interested in calculating the determinant of $(n-1)\times(n-1)$ matrix that is written as a product of $(n-1)\times n$ matrix with an $n\times(n-1)$ matrix (the first matrix comprising of $\vec{C}_{i}$ and the second one involving the derivatives with respect to $\varphi$ of $\vec{P}$).

We have

[TABLE]

[TABLE]

The determinant of this matrix is

[TABLE]

with the vectors $v_{i}$ for $1\leq i\leq n-1$ being an orthonormal collection of $n-1$ vectors perpendicular to the vector $\vec{P}$. Note that each of these vectors is perpendicular to $\vec{P}$ because each $v_{i}$ is obtained by differentiating $\vec{P}$ with respect to $\varphi_{i}$.

Now we have

[TABLE]

since the matrix belongs to $O(n)$. We can write the above determinant as

[TABLE]

where $\vec{P}_{i}$ denotes the $i^{\mathrm{th}}$ component of $i$ and $M_{in}$ denotes the corresponding minor.

Since $\begin{pmatrix}v_{1}&\cdots&v_{n-1}&\vec{P}\end{pmatrix}\in O(n)$, this implies that

[TABLE]

Since $\vec{C}_{i}$ for $1\leq i\leq n-1$ are orthonormal and oriented, we have that

[TABLE]

The same argument as above shows that the vector with the minors $\widetilde{M}_{ni}$ coming from this matrix satisfies

[TABLE]

Now using Cauchy-Binet theorem, (7) is proved. ∎

Now we find the surface measure (5) with respect to the coordinate system $(r,\varphi)$.

Proposition 3.3**.**

The surface measure $\mathrm{d}\Omega$ on the sphere $S(\rho,\alpha)$ with respect to the spherical coordinate system $(r,\varphi)$ is given by

[TABLE]

where $A_{1}$ is defined in (6).

Proof.

We have

[TABLE]

We express $\cos\widetilde{\varphi}_{i}$ for $1\leq i\leq n-1$ in terms of the coordinates $(\varphi_{1},\cdots,\varphi_{n-1})$.

We have:

[TABLE]

Furthermore it is easy to see that

[TABLE]

Let us compute the determinant of the Jacobian of the transformation

[TABLE]

Since

[TABLE]

differentiating this equation, we get

[TABLE]

The Jacobian matrix of (8) is

[TABLE]

Here $\nabla_{\varphi}$ denotes the $(n-1)$-vector $(\frac{\partial}{\partial\varphi_{1}},\cdots,\frac{\partial}{\partial\varphi_{n-1}})$.

The determinant of the matrix above is the same as the determinant of the matrix

[TABLE]

Recall that we are interested in expressing

[TABLE]

in terms of $\mathrm{d}{\varphi}_{1}\cdots\mathrm{d}{\varphi}_{n-1}$. Using Lemma 3.1, we have,

[TABLE]

Note that

[TABLE]

Therefore we have

[TABLE]

Hence

[TABLE]

Now we have

[TABLE]

Multiplying and dividing the left hand side of (3), by $\sqrt{1-(A_{1}^{2}+\cdots+A_{n-4}^{2})}$ and then by $(1-(A_{1}^{2}+\cdots+A_{n-5}^{2}))$ and continuing this way, we get

[TABLE]

Since we are interested in the absolute value of the determinant of the Jacobian of the transformation in (8), we finally have

[TABLE]

This completes the proof. ∎

3.1. Exterior problem

In this section, we prove Theorem 2.2.

We have

[TABLE]

and

[TABLE]

Using Proposition 3.3, we can write the above surface measure $\mathrm{d}\Omega=\frac{\rho^{n-2}r^{2}}{\left|r-A_{1}\right|}\;\mathrm{d}\Omega(\varphi)$, where $\mathrm{d}\Omega(\varphi)=\sin^{n-2}\varphi_{1}\cdots\sin\varphi_{n-2}\mathrm{d}\varphi_{1}\cdots\mathrm{d}\varphi_{n-1}.$ Then

[TABLE]

The integrand in (3.1) is to be interpreted as [math] outside a suitable range of $\varphi$. Now since $r=A_{1}+\sqrt{A_{1}^{2}+\rho^{2}-1}$, we have

[TABLE]

Now we apply Funk-Hecke theorem.

Theorem 3.4** (Funk-Hecke).**

If

[TABLE]

then

[TABLE]

where $\lvert\mathbb{S}^{n-2}\rvert$ denotes the surface measure of the unit sphere in $\mathbb{R}^{n-1}$ and $C_{l}^{\frac{n}{2}-1}$ are the Gegenbauer polynomials.

Using this theorem, we have,

[TABLE]

Making the change of variables $r=x+\sqrt{x^{2}+\rho^{2}-1}$, we have

[TABLE]

This can be written in the form

[TABLE]

where

[TABLE]

This is a Volterra integral equation of the first kind (see [46]). The kernel $K_{l}(\rho,r)$ is continuous together with its first derivatives and $K_{l}(\rho,\rho)\neq 0$ on the interval $(0,R_{1})$, where $0<R_{1}<2$. Equations of this type have a unique solution, which can be obtained through modification to a Volterra equation of the second kind, and then using a resolvent kernel given by Picard’s process of successive approximations (see [38, 45]). This completes the proof of Theorem 2.1.

3.2. Interior problem

Next we prove Theorem 2.2.

Our starting point is:

[TABLE]

We split the integral

[TABLE]

where $\mathbb{S}^{n-1}_{\pm}$ corresponds to those points on the unit sphere such that the line passing through it and the origin intersects a point on the sphere $S(\rho,\alpha)$ corresponding to $r=A_{1}\pm\sqrt{A_{1}^{2}+\rho^{2}-1}$. Let us denote the right hand side of the above equation as $I_{1}+I_{2}$. We have

[TABLE]

Applying Funk-Hecke theorem, this integral is

[TABLE]

Similarly

[TABLE]

Making the change of variables $r=x+\sqrt{x^{2}+\rho^{2}-1}$ in $I_{1}$ and $r=x-\sqrt{x^{2}+\rho^{2}-1}$ and summing up the two integrals, we get,

[TABLE]

This is of the form

[TABLE]

where

[TABLE]

Note that $K_{n}(\rho,\rho)$ does not vanishes in the interval $(0,1-\varepsilon)$ and its derivatives exist and are continuous. The rest of the proof follows exactly as in Theorem 2.1.

3.3. Interior/exterior problem

Finally we prove Theorem 2.3.

Since the argument is exactly as in Theorems 2.1 and 2.2, we will only give the final integral identity. Assume that the function $f$ is supported inside the ball $B(0,R_{2})$ centered at the origin and of radius $R_{2}$, where $R_{2}>2$ and $R_{1}=R_{2}-2$. Suppose the spherical Radon transform data is known along all spheres of radius $\rho$ centered on the unit sphere with $R_{2}-1<\rho<R_{2}+1$, then we have the following Volterra-type integral equation:

[TABLE]

Making a change of variable $\hat{\rho}=R_{2}+1-\rho$ we get

[TABLE]

where

[TABLE]

The rest of the proof follows exactly as before.

4. Three dimensional case

In the numerical simulations below, we specialize to the case of $3$-dimensions. Therefore, in this section, we give the formulas derived earlier for the case of $n=3$.

In this section, for the sake of convenience, we rename the vector $\alpha$ as $(\alpha,\beta)$ and the vector $\varphi$ as $(\varphi,\theta)$. Thus in this section and the next, the point $C$ will be denoted by $(\alpha,\beta)$, more precisely, the Euclidean coordinates of the point $C$ on the unit sphere will be denoted by $(\cos\alpha,\sin\alpha\cos\beta,\sin\alpha\sin\beta)$. A point $P$ on the sphere $S(\rho,\alpha,\beta)$ will be denoted by $(r\cos\varphi,r\sin\varphi\cos\theta,r\sin\varphi\sin\theta)$.

Here the spherical harmonics for $f$ and $Rf=g$ are expanded as

[TABLE]

[TABLE]

In the case of $3$-dimensions, we have that $C_{l}^{(\frac{1}{2})}(x)=P_{l}(x)$, where $P_{l}(x)$ are the Legendre polynomials and $C_{l}^{(\frac{1}{2})}(1)=1$. Therefore the relation between the spherical harmonics coefficients in the three cases are as follows:

1. (Exterior case)

[TABLE] 2. (Interior case)

[TABLE] 3. (Interior/exterior case)

[TABLE]

Note, that in all three cases the kernel $K_{l}(\rho,r)$ is bounded and has a continuous first derivative on the support of (corresponding) $F_{l}^{m}$, and $K_{l}(\rho,\rho)\neq 0$. Hence, these Volterra equations of the first kind can be transformed into equations of the second kind and solved using a resolvent kernel given by Picard’s process of successive approximations (see [38, 45]).

5. Numerical Algorithm

5.1. Generating the Radon data

We consider a generic sphere of integration $S(\rho,\alpha,\beta)$ to be centered at $C=(a_{1},b_{1},c_{1})$ and radius $\rho$ where the center $(a_{1},b_{1},c_{1})$ lies on the sphere of radius $R$. For the interior and exterior cases, we choose $R=1$ and thus use the formulas (4) and (4) derived in the previous sections. For the combined interior and exterior case, we use $R=1.49$ and note that (4) can be easily generalized for acquisition spheres of radius $R$. We consider test phantoms $f$ to be disjoint unions of characteristic functions of balls. To find the spherical Radon transform of $f$, we need to find the surface area of intersection of $S(\rho,\alpha,\beta)$ with $f$. This is equivalent to summing up the surface area of intersection of $S(\rho,\alpha,\beta)$ with characteristic function of each ball. Thus, in the forthcoming calculations, we consider a ball $B$ centered at $(a_{2},b_{2},c_{2})$ and radius $a$.

The sphere $S(\rho,\alpha,\beta)$ and the ball $B$ intersect only when the following conditions do not occur:

[TABLE]

and

[TABLE]

To compute the surface area of intersection of $S(\rho,\alpha,\beta)$ with $B$ in these cases, we first determine the center of the circle of intersection of $S$ and $B$, denoted by $(x_{c},y_{c},z_{c})$. The equation of the plane $\mathcal{P}$ passing through the intersection $S$ and $B$ is given as follows

[TABLE]

The equation of the straight line passing through $(x_{c},y_{c},z_{c})$ and perpendicular to the plane $\mathcal{P}$ is given as follows

[TABLE]

We then can compute

[TABLE]

where $Z=(a_{1}-a_{2})^{2}+(b_{1}-b_{2})^{2}+(c_{1}-c_{2})^{2}.$ Let $d$ be the distance between $(x_{c},y_{c},z_{c})$ from the center of $S$. Then, by an elementary calculation, the surface area of intersection of $S(\rho,\alpha,\beta)$ with $B$, denoted by $\mathcal{S}$, is given as

[TABLE]

Thus

[TABLE]

5.2. Evaluating the spherical harmonics coefficients of the Radon data

After obtaining the Radon data $g(\rho,\alpha,\beta)$, we need to determine $g_{l}^{m}(\rho),~{}l=0,\ldots,\infty,$ $~{}m=-l,\ldots,l$ given by

[TABLE]

where $\bar{Y}^{m}_{l}(\alpha,\beta)=(-1)^{-m}Y^{-m}_{l}(\alpha,\beta).$ This is done numerically by following the method employed in [13]. Given $g(\rho,\alpha,\beta)$ at $\alpha_{j}=\pi j/2N,\beta_{k}=\pi k/N,j,k=0,\ldots,2N-1$, we compute

[TABLE]

where

[TABLE]

5.3. Inversion of the integral equations (4), (4) and (4)

To solve the integral equations numerically, we use the product trapezoidal method as found in [37, 40, 47]. In this method, the integral equations are discretized and the integrands are approximated by product trapezoidal rule. This in turn leads to matrix-vector equations and thus, we obtain discrete solutions of the discretized integral equations, provided the matrices are invertible. In the following section, we provide the corresponding matrix-vector equations for solving (4), (4) and (4) corresponding to the exterior, interior and the combined interior/exterior problems respectively and prove their invertibility.

5.3.1. Exterior case

We rewrite (4) as

[TABLE]

where

[TABLE]

We discretize $\rho\in(0,1)$ into $M+1$ equidistant points of interval length $h$ as $\rho_{i},~{}i=0,\ldots,M$. The corresponding matrix-vector equation is given as follows

[TABLE]

where

[TABLE]

and $A_{E}=(a_{ik})$ where

[TABLE]

Note that $A_{E}$ is a lower triangular matrix, and since

[TABLE]

we have that $A_{E}$ is invertible.

5.3.2. Interior Case

We again discretize $\rho\in(0,1)$ into $M+1$ equidistant points as $\rho_{i},~{}i=0,\ldots,M$ and obtain the following matrix-vector equation

[TABLE]

where

[TABLE]

and $A_{I}=(a_{ik})$ where

[TABLE]

Therefore, if $\rho_{0}>0$

[TABLE]

Thus $A_{I}$ is invertible.

5.3.3. Interior/Exterior Case

In a similar way as in the previous two cases, we discretize $\rho\in(0,2R)$ into $M+1$ equidistant points as $\rho_{i},~{}i=0,\ldots,M$ to obtain the following matrix-vector equation

[TABLE]

where

[TABLE]

and $A_{IE}=(a_{ik})$ where

[TABLE]

Therefore, if $\rho_{0}>0$

[TABLE]

Thus $A_{IE}$ is invertible.

The following theorem states the error estimate for the numerical solution of the integral equations (4), (4) and (4) which follows from [31, Thm. 7.2].

Theorem 5.1** (Error Estimates).**

Let $f_{l}^{m,\mathrm{exact}}$ be the $C^{3}$ solution of (4) [(4), (4) in $[0,R]$ and $f_{l}^{m}$ be the solution to 18 [or (20) and (22) resp.]. Then

[TABLE]

where $\|\cdot\rVert_{2}$ represents the discrete version of the continuous $L^{2}$ norm in $[0,R]$ (see for e.g., [14, Ch. 4]).

To solve the matrix equations (18), (20) and (22), we need to invert the matrices $A_{E},A_{I},A_{IE}$. It turns out that the condition numbers of these matrices are greater than $10^{4}$ for almost all values of $l,m$. It is well known that numerically inverting a matrix with condition number$10^{r}$ leads to a loss of $r$ digits of accuracy [23]. Thus for inversion, we use the technique of Truncated Singular Value Decomposition (TSVD), originally proposed in [21]. See also [9, 40].

6. Numerical Results

In this section we show the results of the numerical computations performed for the inversion of spherical transforms described in Section 2 with functions supported in interior, exterior and both interior and exterior of the acquisition sphere. We discretize $\rho\in[\epsilon,R-\epsilon]$, with $\epsilon=0.001$, into 50 equally spaced grid points, $\alpha,\theta\in[0,\pi]$ and $\beta,\phi\in[0,2\pi]$ into 100 equally spaced grid points for all our computations. As mentioned before, for the interior and the exterior cases, the value of $R=1$ whereas for the combined interior and exterior case, the value of $R=1.49$.

6.1. Functions supported inside the acquisition sphere

Figures 2a and 2b show the horizontal and vertical cross sections of a phantom represented by a ball centered at $(0.5,0,0)$ and radius 0.3. Figure 2c and 2d shows the horizontal and the vertical cross sections of the reconstructed phantom. We note the good recovery in this case.

To demonstrate the robustness of our algorithm, we also tested it on the spherical Radon data with 5% multiplicative Gaussian noise. The results are shown in Figures 3a and 3b. We again note the good recovery in presence of noisy data.

We also applied our reconstruction algorithm to a phantom whose support is inside the acquisition sphere and contains the origin. Notice, that our result about uniqueness of the inversion (Thm. 2.2) does not cover this case, since here the kernels $K_{l}(\rho,r)$ of the integral equations appearing in the proof vanish, when $\rho=r=1$ (see eq. (4)). Hence, one may not expect stable recovery in the numerical method. And indeed, the reconstructed image in Figure 4b shows instability near the origin.

6.2. Functions supported outside the acquisition sphere

Figures 5a and 5b show the horizontal and vertical cross sections of a phantom represented by two balls centered at $(-1.5,0,0)$ and $(1.5,0,0)$ with radius 0.2 and 0.3 respectively. Figures 5c and 5d shows the horizontal and the vertical cross sections of the reconstructed phantom. Microlocal analysis arguments show that the entire spherical shell of the balls cannot be constructed stably with the given spherical Radon transform data. We see the presence of an increased number of artifacts in contrast to the interior case. The reconstructions are consistent with this analysis.

6.3. Functions supported on both sides of the acquisition sphere

Figure 6a shows the horizontal cross section of a phantom represented by two balls centered at $(0.5,0,0)$ and $(-2.0,0,0)$ with radius 0.2 and 0.3 respectively. The balls lie on either side of the acquisition sphere, which is centered at the origin with radius 1.49. Figure 6b shows the horizontal cross section of the reconstructed phantoms. Again by microlocal analysis arguments, the ball outside the acquisition sphere cannot be constructed stably whereas the ball inside the acquisition sphere can be constructed stably. This is depicted in the reconstructions.

7. Conclusion

We studied the problem of inverting the spherical Radon transform in spherical geometry of data acquisition with incomplete radial data. Such problems arise in image reconstruction procedures in photo- and thermo-acoustic tomography, ultrasound reflection tomography, as well as in radar and sonar imaging. We considered three distinct scenarios of the location of the support of the image function: strictly inside the acquisition sphere (interior problem), strictly outside (exterior problem), and both inside and outside (interior/exterior problem). For all three cases we provided a constructive proof of the uniqueness of inversion of SRT from incomplete radial data and obtained an iterative procedure to recover the image function. We presented a robust computational algorithm based on our inversion procedure and demonstrated its accuracy and efficiency on several numerical examples.

Acknowledgments

Ambartsoumian was supported in part by US NSF Grants DMS 1109417, DMS 1616564 and Simons Foundation Grant 360357.

Gouia-Zarrad was supported in part by the American University of Sharjah (AUS) research grant FRG3.

Krishnan was supported in part by NSF grants DMS 1109417 and DMS 1616564. He and Roy benefited from support of the Airbus Group Corporate Foundation Chair “Mathematics of Complex Systems” established at TIFR Centre for Applicable Mathematics and TIFR International Centre for Theoretical Sciences, Bangalore, India.