An approach to the spherical mean Radon transform with detectors on a line

TL;DR

This paper introduces a new consistency-based method for inverting the 2D spherical Radon transform with detectors on a line, providing a practical iterative algorithm for functions supported on one side of the line, relevant to tomography.

Contribution

It proposes a novel consistency method and iterative formula for approximate inversion of the 2D spherical Radon transform with line detectors, addressing the lack of exact formulas.

Findings

Proved local description of the reconstruction process.

Developed a practical iterative algorithm for functions on one side of a line.

Applicable to thermo- and photo-acoustic tomography problems.

Abstract

The article suggests a new approach what is called a consistency method for the inversion of the spherical Radon transform in 2D with detectors on a line. It is known that there is not an exact inversion formula in 2D. By means of the method was proved that the reconstruction has a local description and found a new iteration formula which give an practical algorithm to recover an unknown function supported completely on one side of a line $L$ from its spherical means over circles centered on the line $L$. Such an inversion is required in problems of thermo- and photo-acoustic tomography.

Equations

$$Mf(p,r)=\frac{1}{2\pi }{\int }_{{\mathbf{S}}^1}f(p+r\omega )d\phi ,\text{for}(p,r)\in L\times [0,\infty ).$$

$$(p,r,\phi )\text{and}(x,y,\varphi ),$$

$$G(x,y,\varphi )=G(x,y)$$

$$G(p,r,\phi )\to G(x,y,\varphi )=G(x,y)=f(x,y)$$

$$Z_{n,i}(t)=\sum _{j=1}^{n+i}{z}_j(2n,2i)t^{2j-1},t\in [0,1]$$

$$z_j(2n,2i)=\sum _{k=i,j-i\le k}^n{A}_j(2k,2i)\text{for}1\le j\le n+i.$$

$$f(x,y)=\underset{n\to \infty }{\lim }\left(2(n+1)Mf((x,0),y)+\sum _{i=0}^n{\int }_0^y{y}^{2i-1}{Z}_{n,i}(\frac{u}{y})(Mf((x,0),u){)}_x^{(2i)}du\right),$$

$$f(p,r,\phi ),\phi \in [-\pi ,\pi ]$$

$$f(p,r,\phi )=\sum _{k=0}^{\infty }\left(a_k(p,r)\cos k\phi +b_k(p,r)\sin k\phi \right).$$

$$f(p,r,\phi )=Mf(p,r)+\sum _{k=1}^{\infty }\left(a_k(p,r)\cos k\phi +b_k(p,r)\sin k\phi \right).$$

$$\{\begin{array}{l}x=p+r\cos \varphi \\ y=r\sin \varphi \\ \varphi =\phi +\frac{\pi }{2}\end{array}$$

$${\phi }_{\varphi }^{\prime }=1,r_{\varphi }^{\prime }=\frac{r\sin \phi }{\cos \phi },p_{\varphi }^{\prime }=\frac{r}{\cos \phi }.$$

$$f_{r}r\sin \phi +f_{p}r+f_{\phi }(\cos \phi -r)=0.$$

$$(f(x,y){)}_{\varphi }^{\prime }=(f(p,r,\phi ){)}_{\varphi }^{\prime }=0.$$

$${\int }_{-\pi }^{\pi }\left(f_{r}r\sin \phi +f_{p}r+f_{\phi }(\cos \phi -r)\right)\cos k\phi d\phi =0.$$

$${\int }_{-\pi }^{\pi }\left(f_{r}r\sin \phi +f_{p}r+f_{\phi }(\cos \phi -r)\right)\sin k\phi d\phi =0.$$

$$\begin{gathered} {\int }_{-\pi }^{\pi }([Mf(p,r)+\sum _{i=1}^{\infty }(a_i(p,r)\cos i\phi +b_i(p,r)\sin i\phi ){]}_r^{\prime }r\sin \phi + \\ r[Mf(p,r)+\sum _{i=1}^{\infty }(a_i(p,r)\cos i\phi +b_i(p,r)\sin i\phi ){]}_p^{\prime }+ \\ [Mf(p,r)+\sum _{i=1}^{\infty }(a_i(p,r)\cos i\phi +b_i(p,r)\sin i\phi ){]}_{\phi }^{\prime }(\cos \phi -r))\cos k\phi d\phi =0. \end{gathered}$$

$$r(b_1(p,r){)}_r^{\prime }+b_1(p,r)=-2r(Mf(p,r){)}_p^{\prime };$$

$$r(b_2(p,r){)}_r^{\prime }+2b_2(p,r)+2r(a_1(p,r){)}_p^{\prime }=0;$$

$$r(b_k(p,r){)}_r^{\prime }+k{b}_k(p,r)-r(b_{k-2}(p,r){)}_r^{\prime }+(k-2)b_{k-2}(p,r)+2r(a_{k-1}(p,r){)}_p^{\prime }=0.$$

$$r(a_2(p,r){)}_r^{\prime }+2a_2(p,r)-2r(b_1(p,r){)}_p^{\prime }=2r(Mf(p,r){)}_r^{\prime };$$

$$r(a_k(p,r){)}_r^{\prime }+k{a}_k(p,r)-r(a_{k-2}(p,r){)}_r^{\prime }+(k-2)a_{k-2}(p,r)-2r(b_{k-2}(p,r){)}_p^{\prime }=0.$$

$$\{\begin{array}{l}r(b_1(p,r){)}_r^{\prime }+b_1(p,r)=-2r(Mf(p,r){)}_p^{\prime }\\ r(b_2(p,r){)}_r^{\prime }+2b_2(p,r)+2r(a_1(p,r){)}_p^{\prime }=0\\ r(a_2(p,r){)}_r^{\prime }+2a_2(p,r)-2r(b_1(p,r){)}_p^{\prime }=2r(Mf(p,r){)}_r^{\prime },\end{array}$$

$$\{\begin{array}{l}r(b_k(p,r){)}_r^{\prime }+k{b}_k(p,r)-r(b_{k-2}(p,r){)}_r^{\prime }+(k-2)b_{k-2}(p,r)+2r(a_{k-1}(p,r){)}_p^{\prime }=0\\ r(a_k(p,r){)}_r^{\prime }+k{a}_k(p,r)-r(a_{k-2}(p,r){)}_r^{\prime }+(k-2)a_{k-2}(p,r)-2r(b_{k-2}(p,r){)}_p^{\prime }=0,\end{array}$$

$$f(x,y)=f(x,y,\varphi )=Mf(p,r)+\sum _{k=1}^{\infty }\left(a_k(p,r)\cos k\phi +b_k(p,r)\sin k\phi \right).$$

$$f(x,y)=Mf(p,r)+\sum _{k=1}^{\infty }{a}_k(p,r).$$

$$\sum _{k=1}^{\infty }{a}_k(p,r).$$

$$0=f(p,r,-\pi )=f(p,r,\pi )=Mf(p,r)+\sum _{k=1}^{\infty }\left(a_k(p,r)\cos k\pi +b_k(p,r)\sin k\pi \right),$$

$$Mf(p,r)+\sum _{k=1}^{\infty }(-1{)}^k{a}_k(p,r)=0.$$

$$f(x,y)=2Mf(p,r)+\sum _{k=1}^{\infty }{a}_{2k}(p,r),$$

Full text

affiliation: 1. Russian-Armenian University

2. Institute of Mathematics NAS RA.

e-mail: rafikaramyan@yahoo.com

Keywords: Tomography, thermoacoustic tomography, spherical Radon transform, inverse problem.

1 Introduction and formulation of the problem

Medical tomography has had a huge impact on medical diagnostics. The classical Radon transform maps a function to its integrals over straight lines and serves as the basis of x-ray Computer Tomography. Recently researchers have been developing novel methods that combine different physical types of signals. The most successful example of such a combination is the thermoacoustic tomography (TAT). Thermoacoustic theory has been discussed in many literature reviews such as [13]-[24]. Briefly TAT procedure is: a short-duration electromagnetic (EM) pulse is sent through a biological object with the aim of triggering a thermoacoustic response in the tissue. The amount of energy absorbed at a location $X$ strongly depends on the local biological properties of the cells. Thus, if the energy absorption distribution function $f$ were known, it would provide a great diagnostic tool. The acoustic wave which is the result of the thermoelastic expansion can be measured by transducers placed outside the object (assuming the sound speed c constant). Thus, one effectively measures the integrals of $f$ over all spheres centered at the transducers’ locations. To recover $f$ one needs to invert the so-called spherical Radon transform of $f$ that integrates a function over all such spheres.

We denote by $\mathbf{R}^{n}$ ($n\geq 2$) the Euclidean $n$ - dimensional space. Let ${\mathbf{S}^{n-1}}$ be the $n-1$ dimensional unit sphere in $\mathbf{R}^{n}$ with the center at the origin $O\in\mathbf{R}^{n}$, $\sigma_{n-1}$ its total surface measure. By $S(p,r)$ we denote the sphere of radius $r>0$ centered at $p\in\mathbf{R}^{n}$.

The above motivated the study of the following mathematical problem. For a continuous, real valued function $f$ supported in a compact region $G$, we are interested in recovering $f$ from the mean value $Mf(p,r)$ of $f$ over spheres $S(p,r)$ centered on $L$; that is, given $Mf(p,r)$ for all $p\in L$ and $r>0$, we wish to recover $f$.

In order to implement the TAT reconstruction the following problems arise. For which sets $L$ the data collected by transducers placed along $L$ is sufficient for unique reconstruction of $f$ (set $L$ is called a set of injectivity if the transform (1.1) is injective) and what are inversion formulas.

Agranovsky and Quinto in [1], [2] have proved several significant uniqueness results for the spherical Radon transform. In [1] they gave a complete characterization of sets of uniqueness (sets of centers) for the circular Radon transform on compactly supported functions in the plane. In [5], was provided a complete range description in case of circular Radon transform in 2D.

Obviously any line L (or a hyperplane in higher dimensions) is a non-uniqueness set, since any function $f$ odd with respect to $L$ will clearly produce no signal: $Mf(p,r)=0$. On the other hand (see [26], [27]), if $f$ is supported completely on one side of the the line $L$ (the standard situation in TAT), it is uniquely recoverable from its spherical means centered on $L$, and thus from the observed data.

Exact inversion formulas for the spherical Radon transform are currently known for boundaries of special domains, including spheres, cylinders and hyperplanes ([8], [4], [12], [16], [17], [14], [15], [26]).

In this paper for a continuous, real valued function $f$ defined in $\mathbf{R}^{2}$ and supported completely on one side of a line $L$, we are interested in recovering $f$ from the mean value of $f$ over circles centered on $L$. The article suggests a new approach what is called a consistency method for the inversion of the spherical Radon transform in 2D with detectors on a line. By means of the method a new iteration formula was found which give an practical algorithm to recover an unknown function supported in a compact region from its spherical means over circles centered on a line outstand the region. Also was proved that reconstruction has a local description (see Theorem 1 below).

The consistency method, suggested by the author of the paper, first was applied in [9] (see also [10], [11]) to inverse generalize Radon transform on the sphere.

Note, one can apply the consistency method to inverse the spherical Radon transform for dimensions $n\geq 3$. Also, one can consider to apply the consistency method to inverse the spherical Radon transform for different geometries of transducers.

Now we consider the circular Radon transform on the plane. For a continuous function $f$ supported in the compact region $G\in\mathbf{R}^{2}$ we have (see (1))

[TABLE]

Here $d\varphi$ is the circular Lebesgue measure on ${\mathbf{S}^{1}}$, $\omega=(\sin\varphi,\cos\varphi)$. The value $Mf(p,r)$ is the average of $f$ over the circle $S(p,r)$ with center $p\in L$ and radius $r>0$.

We consider the restriction of $f$ onto the circle ${S}(p,r)$ for $(p,r)\in L\times[0,\infty)$.

A pair say $(S(p,r),q)$ where $q\in{S}(p,r)$ we call a circular flag (in integral geometry there is a concept of a flag an ordered pair of orthogonal unit vectors [6],[7]). There are two equivalent representations of a circular flag $(S(p,r),q)$ where $q\in{S}(p,r)$, dual each other:

[TABLE]

where $p\in L$ and $\varphi\in[-\pi,\pi]$ is the angular coordinate of $q\in{S}(p,r)$ measured from the direction perpendicular to $L$, while $x,y$ are the Euclidean coordinates of $q$ and $\phi$ is the direction (the angular coordinate) of $\overrightarrow{qp}$.

Thus one can represent the restriction of $f$ onto the circle ${S}(p,r)$ by $f(p,r,\varphi)$, where $\varphi\in[-\pi,\pi]$ is the angular coordinate of $q\in{S}(p,r)$.

The idea of the method is the following: for every $(p,r)$ the equation (1.1) reduces to an integral equation on the circle $S(p,r)$. The general solution of the reduced integral equation we write in terms of Fourier series expansion with unknown coefficients. Let $G(p,r,\cdot)$ be a solution of the reduced integral equation for $(p,r)$.

***Definition *1. If $G$ written in dual coordinates satisfies

[TABLE]

(no dependence on the variable $\phi$), then $G$ is called a consistent solutions.

There is a principle: each consistent solutions $G$ produces via the map

[TABLE]

the solution of (1.1), and vice versa. Conversely, restrictions of the solution of (1.1) onto the circles $S(p,r)\,$, ($(p,r)\in L\times[0,\infty)$) is a consistent solutions.

Hence the problem of finding the solution of (1.1) reduces to finding the consistent solutions of the reduced equations (1.1).

In this paper was proved the following theorem. Let $f$ be a continuous, real valued function supported in the compact region $G$ located on one side of the the line $L$. On the plane consider usual cartesian system of coordinate choosing $L$ as the $x$-axis. $Mf(p,r)$ is the average of $f$ over a circle with center $p\in L$ and radius $r>0$.

Theorem 1

Let $f$ be a continuous, real valued function supported in the compact region $G$ located on one side of the line $L$. For $(x,y)\in G$ the value $f(x,y)$ depends on values Mf on a neighborhood of $p=(x,0)\in L$ and $0\leq r\leq y$.

Now we describe the inversion formula. We define a sequence of standard polynomials $Z_{n,i}$ defined on the interval $[0,1]$, where $n,i$ are integers and $0\leq i\leq n$:

[TABLE]

with coefficients

[TABLE]

In §5 was found recurrent relations by means of which one can find the coefficients $A_{j}(2k,2i)$ for integers $k\geq 1$, $0\leq i\leq k$ and $1\leq j\leq k+i$. We call $Z_{n,i}$ standard polynomials because their construction does not depend on $f$.

Theorem 2

Let $f$ be an infinitely differentiable real valued function $f$ supported in the compact region $G$ located on one side of the line $L$. For $(x,y)\in G$ we have

[TABLE]

here $(Mf((x,0),u))^{(2i)}_{x}$ is the derivative of order $2i$ with respect the variable $x$

($(Mf((x,0),u))^{(0)}_{x}=Mf((x,0),u)$).

Theorem 2 suggests a practical algorithm to reconstruct $f$.

2 General solution of the reduced equation

(1.1)

For a fix $(p,r)\in L\times[0,\infty)$ the restriction of $f$ onto the circle $S(p,r)$ we write in the form

[TABLE]

where $\varphi$ is the angular coordinate of $q\in S(p,r)$ measured from the direction perpendicular to $L$. On the plane we consider usual cartesian system of coordinate choosing $L$ as the $x$-axis and below the point $(p,0)\in L$ we will identify with $p$.

It is known that periodic, continuous, with piecewise-continuous first-derivative function $f$ can be written as its Fourier series expansion. For any $(p,r)$ the Fourier series expansion of the restriction $f(p,r,\varphi)$ is

[TABLE]

Taking into account (1.1) we have

[TABLE]

Now we are going to write $f(p,r,\varphi)$ in dual coordinates.

The transform $(p,r,\varphi)\longrightarrow(x,y,\phi)$ (see (1.2)) can be represented by the following system

[TABLE]

We denote the (partial) derivative of a function $f$ with respect to a variable say $v$ by $f^{\prime}_{v}$. From (2.4) we get the following expressions for the derivatives

[TABLE]

3 The consistency condition

Now we consider the coefficients $a_{k}(p,r)$,  $b_{k}(p,r)$ ($k=1,2,...$) in (2.3) as functions of $(p,r)\in(-\infty,\infty))\times[0,\infty)$ and try to find them from the consistency condition. We write $f(p,r,\varphi)$ in dual coordinates and require that the right side should not depend on $\phi$ for every $(x,y)\in\mathbf{R}^{2}$. It follows from (2.4) the following theorem.

Theorem 3

Let $f$ be an infinitely differentiable real valued function $f$ supported in the compact region $G$ located on one side of the line $L$. We have

[TABLE]

Proof of the Theorem 3 follows from (2.5) and the condition that

[TABLE]

All Fourier coefficients of the function on the left side of (3.1) equals [math] as the function identity equals [math]. We have

[TABLE]

and

[TABLE]

Substituting (2.3) into (3.3) we get

[TABLE]

From (3.5) we obtain the following differential equation for $b_{k}(p,r)$.

For $k=1$

[TABLE]

for $k=2$

[TABLE]

for $k>2$

[TABLE]

By analogous way substituting (2.3) into (3.4) for $a_{k}(p,r)$ we get. For $k=2$

[TABLE]

for $k>2$

[TABLE]

Thus we obtain the following system of differential equations for unknown coefficients $a_{k}(p,r)$,  $b_{k}(p,r)$ ($k=1,2,...$)

[TABLE]

for $k=1$, $k=2$ and

[TABLE]

for $k>2$.

4 The consistent solution

Returning to the formula (2.3) for every $(x,y)$ we have

[TABLE]

$\phi=\pi/2$ corresponds to $\varphi=0$ the case that $p\in L$ is the projection of $(x,y)$ onto $L$, hence we have $p=x$. Thus for $p=x$ and $r=y$ from (4.1) we have

[TABLE]

Now the problem is to calculate

[TABLE]

Taking into account that $f$ supported in the compact region $G$ located on one side of the line $L$ for any $(p,r)$ we have

[TABLE]

hence

[TABLE]

From (4.2) and (4.4) we have

[TABLE]

where $p=x$ and $r=y$.

Now the problem is to calculate

[TABLE]

Taking the sums of the second equations of (3.12) for even $k$ we get the following recurrent equations for $a_{2k}(p,r)$.

[TABLE]

and we see that to calculate $a_{2k}(p,r)$ we need to known coefficients $b_{k}(p,r)$ for odd $k$.

Taking the sums of the first equations of (3.12) for odd $k$ we get the following recurrent equations for $b_{2k-1}(p,r)$.

[TABLE]

Using equations (4.7) and (4.8) one can calculate step by step the unknown coefficients $a_{2k}(p,r)$ and the unknown coefficients $b_{2k-1}(p,r)$ (at first we find $b_{1}$, next $a_{2}$, next $b_{3}$ next $a_{4}$ and so on).

5 Solution of the differential equations

Now we are going to find a boundary conditions for the differential equations. We have for $k\geq 1$

[TABLE]

and

[TABLE]

Taking into account that $f$ supported in the compact region $G$ located on one side of the line $L$ we get the following boundary conditions

[TABLE]

Thus we get two systems of differential equations (4.7) and (4.8) with boundary conditions (5.3). The unique solution of (4.7) for $k\geq 1$ is

[TABLE]

[TABLE]

The unique solution for $b_{1}(p,r)$ is

[TABLE]

The unique solution of (4.8) for $k\geq 2$ is

[TABLE]

[TABLE]

Note that it follows from (4.5), (5.4) and (5.6) that for $(x,y)\in G$ the value $f(x,y)$ depends on values Mf on a neighborhood of $p=(x,0)\in L$ and $0\leq r\leq y$. Theorem 1 is proved.

Lemma 1

There are sequences of polynomials defined on $[0,1]$

1) $A_{2k,2i}$ of degree $2k-1+2i$ for integers $k\geq 1$ and $0\leq i\leq k$

[TABLE]

2) $B_{2k-1,2i-1}$ of degree $2k+2i-3$ for integers $k\geq 1$ and $1\leq i\leq k$

[TABLE]

such that

[TABLE]

and

[TABLE]

here and below $(Mf(p,u))^{(j)}$ is the derivative of order $j$ with respect the variable $p$ ($(Mf(p,u))^{(0)}=Mf(p,u)$).

Proof 1

Mathematical induction can be used to prove Lemma 1. It follows from (5.5) that for $k=1$ (5.9) and (5.10) are true. Indeed

[TABLE]

from (5.4) using (5.11) we have

[TABLE]

[TABLE]

Suppose (5.9) and (5.10) are true for some $n=k$. Prove that (5.9) and (5.10) are true for $n=k+1$. From (5.6) we have

[TABLE]

[TABLE]

Substituting the expressions for $b_{2j-1}(p,u)$ and $a_{2j}(p,u)$ from (5.9) and (5.10) into (5.13) we obtain

[TABLE]

Changing the order of summation in (5.14) and the order of integration we get

[TABLE]

After grouping of summands finally we obtain

[TABLE]

where for $i=1$

[TABLE]

with

[TABLE]

for $1<i\leq k$

[TABLE]

with

[TABLE]

for $i=k+1$

[TABLE]

with

[TABLE]

Now lets prove Lemma 1 for $a_{2(k+1)}$. From (5.4) we have

[TABLE]

Substituting the expressions for $b_{2j-1}(p,u)$ and $a_{2j}(p,u)$ from (5.9) and (5.10) into (5.23) we obtain

[TABLE]

Changing the order of summation in (5.24) and the order of integration we get

[TABLE]

After grouping of summands finally we obtain

[TABLE]

where for $i=0$

[TABLE]

with

[TABLE]

for $1\leq i\leq k$

[TABLE]

with

[TABLE]

for $i=k+1$

[TABLE]

with

[TABLE]

Also note that we obtain recurrent relations (5.18), (5.20), (5.22), (5.28), (5.30), (5.32) between coefficients $A_{j}(2k,2i)$ and $B_{j}(2k-1,2i-1)$ for $k>1$ and $i\leq k$.

The first few coefficients are:

[TABLE]

Using recurrent relations one can calculate all coefficients of $A(2n,2i)$ and

$B(2n-1,2i-1)$ for $i\leq n$ by means of coefficients $A(2k,2i)$ and $B(2k-1,2i-1)$ for $1\leq k<n$.

6 Partial sums of the series

We are going to consider the partial sums of the series (see (4.6))

[TABLE]

Taking into account (5.9) we have

[TABLE]

Changing the order of summation in (6.2) we obtain

[TABLE]

Note that here we assume $A_{0,0}\equiv 0$

We denote by $Z_{n,i},$ $0\leq i\leq n$ the following polynomial of degree $2n+2i-1$ defined on the interval $[0,1]$

[TABLE]

Substituting (5.7) into (6.4) we obtain

[TABLE]

Substituting (6.3) and (6.5) into (4.5) we obtain

[TABLE]

where $p=x$, $\,r=y$ and

[TABLE]

are polynomials with coefficients

[TABLE]

Note that one can find the coefficients $A_{j}(2k,2i)$ from recurrent relations (5.18), (5.20), (5.22), (5.28), (5.30), (5.32). Theorem 2 is proved.

7 Implementation of the reconstruction

formula

The problem of reconstructing a function from spherical means is important for many imaging and remote sensing applications (see, for example, [13], [20], [3]). These applications require inversion of the spherical Radon transform. However, those formulas require continuous data, whereas in practical applications only a discrete data set is available. In some tomographic applications iterative reconstruction algorithms are more common. In spite of absence of exact FBP formulas in 2D, approximate ones that preserve all the singularities of the image can be easily written and then improved by successive iterative corrections. However, due to the presence of the derivative, the inversion formulas are sensitive to error in the data $Mf$ (see [18]).

In the present paper, we have established a new iterative reconstruction algorithm to recover a function $f$ supported in a compact region from its spherical means $Mf$ (see Theorem 2) which is different from the existing ones in [14], [15], [25]. Our reconstruction formula can be numerically implemented due to a local description. Thus, when evaluating $f(x,y)$ for $(x,y)\in B(O,R)$ we just compute the integral for the set of frequencies uniformly distributed over the interval $[0,y]$. However, due to the presence of the derivatives of higher order, the inversion formulas are sensitive to error in the data $Mf$. In some cases for the derivatives one can use their analytic expressions. To estimate the iteration speed we use the following known result from the theory of Fourier series expansion. Let f be $2\pi$-periodic, continuous, with piecewise-continuous first-derivative function. Then the Fourier series of $f$ converges uniformly

[TABLE]

where $S_{n}(f,x)$ is the partial sum of the Fourier series of $f$ and $c$ does not depend on $n$.

Note that, we get not only uniform convergence, but also a rate of convergence. Thus by finding polynomials $Z_{n,i}$ (see (1.6)) for large $n$ one can recover a function $f$ by approximation as close as we want.