The Radon Transform over Cones with Verticeson the Sphere and Orthogonal Axes

The Radon Transform over Cones with Vertices on the Sphere and Orthogonal Axes

Abstract

Recovering a function from its integrals over circular cones recently gained significance because of its relevance to novel medical imaging technologies such emission tomography using Compton cameras. In this paper we investigate the case where the vertices of the cones of integration are restricted to a sphere in -dimensional space and symmetry axes are orthogonal to the sphere. We show invertibility of the considered transform and develop an inversion method based on series expansion and reduction to a system of one-dimensional integral equations of generalized Abel type. Because the arising kernels do not satisfy standard assumptions, we also develop a uniqueness result for generalized Abel equations where the kernel has zeros on the diagonal. Finally, we demonstrate how to numerically implement our inversion method and present numerical results.

Keywords: Computed tomography; Radon transform; SPECT; Compton cameras; conical Radon transform; uniqueness of reconstruction; spherical harmonics decomposition; series expansion; generalized Abel equations; first kind Volterra equations with zeros in diagonal.

AMS Subject Classification: 44A12, 45D05, 92C55.

\setitemize

label= \setenumeratelabel=()

1 Introduction

Many tomographic imaging modalities are based on the inversion of Radon transforms, which map a function onto its integrals over certain surfaces in (see, for example, [23, 31]). The most basic example is the classical Radon transform which maps the function onto its integrals over hyperplanes and which is the mathematical basis of X-ray CT. Another well investigated example is the spherical Radon transform, which maps a function onto its integrals over hyper-spheres and which finds application in the recently developed photoacoustic tomography [13, 14, 24]. In this article we consider the conical Radon transform that maps a function to its integrals over circular half cones. The conical Radon transform recently gained increased interest, mainly due to its relevance for SPECT using Compton cameras (see, for example, [1, 6, 8, 19, 21, 28, 29, 33, 38, 42]).

\psfrag

B \psfragw

1.1 SPECT using Compton cameras

Single-photon emission computed tomography (SPECT) is a well established medical imaging technology for functional imaging. In SPECT, weakly radioactive tracers are given to patient and participate in the physiological processes. The radioactive tracers can be detected through the emission of gamma ray photons which provide information about the interior of patient. In order to obtain location information on the emitted photons, the standard approach in SPECT uses collimators which only record photons that enter the detector vertically. As illustrated in Figure 1(a), such data provide values of line integrals of the tracer distribution.

A major drawback of using collimators is that they remove most photons. Therefore the number of recorded photons is low and the noise level high. Typically, only one out of photons emitted from the patient is actually detected with this standard approach. In order to increase the number of recorded photons, the concept of Compton cameras has been developed in [11, 37, 41]. As illustrated in Figure 1(b), a Compton camera consists of a scatter detector array D1 and an absorption detector array D2. A photon emitted in the direction of the camera undergoes Compton scattering in D1, and is absorbed in D2. The required distinction of individual photons is obtained by coincidence detection. Both detectors are position and energy sensitive, and the measured energies can be used to determine the scattering angle [37]. Using such information, one concludes that the detected photon must have been emitted on the surface of a circular cone, where the vertex is given by the position at D1, the central axis points from the position on D2 to the position on D1, and the opening angle equals the Compton scattering angle. Consequently, for a distribution of tracers, the Compton camera approximately provides integrals of the marker distribution over conical surfaces.

1.2 Inversion of the conical Radon transform

As outlined above, SPECT with Compton cameras yields to the conical Radon transform that maps a function modeling the marker distribution to the surface integrals over right circular half cones

 C(z,β,ψ)={z+rω∣r≥0 and ω∈S2 with ω\vbox\scalebox{.6}{∙}β=cosψ}.

Here is the vertex of the cone, the direction of the central axis, and the half opening angle. Variants of the conical Radon transform in are known as V-line or broken-ray transforms. These transforms appear in emission tomography with one-dimensional Compton cameras [5, 22], or in the recently developed single scattering optical tomography [16]. In this paper, we consider the conical Radon transform in general dimension and further include a radial weight, that can be adjusted to a particular application at hand. In previous work on Compton camera imaging [38], models with and without radial weight have been proposed and used.

The conical Radon transform depends on six parameters , whereas the function only depends on three spatial coordinates. Therefore the problem of reconstructing from its integrals over all circular cones is highly overdetermined. Several authors have studied the problem of inverting the function from integrals over particular subsets of all cones. In SPECT with Compton cameras, the vertex is naturally fixed to the scattering surface D1. In the case where D1 is a plane and the axis is fixed to , Fourier reconstruction formulas have been derived in [8, 32]. Formulas of the filtered backprojection type have been derived in [19, 28]. The case of variable axis and vertices restricted to a surface has been considered in [6, 21, 27, 33, 38, 39, 42]. See also [1, 3, 15, 17] for related results on different conical transforms.

To the best of our knowledge, if the set of vertices is different from a plane and any vertex is associated with a single symmetry axis, no results are known for reconstructing a function from its integrals over such cones. In this paper we develop an inversion approach for the case when is a sphere and the symmetry axes of the cones are orthogonal to the sphere. We derive a reconstruction procedure based on spherical harmonics decomposition and show invertibility of the considered transform. Spherical harmonics decompositions have been previously used for studying other Radon transforms. See, for example, [7, 9, 26, 31] for the classical Radon transform, [35] for a weighted Radon transform over planes, [2] for the circular Radon transform, or [3, 4] for a broken ray transform with vertices in a disc. In these works, the arising generalized Abel equations satisfy all conditions needed in order to apply standard well-posedness results. For the transform we study, one basic assumption of these results is violated, namely the kernels turn out to have zeros on the diagonal (see Theorem 3.2). Nevertheless, we are able to show solution uniqueness; see Theorem 3.5.

1.3 Outline

The paper is organized as follows. In Section 2 we define the conical Radon transform with vertices on the sphere and orthogonal axis, and derive some elementary results for that transform. Our main results are stated in Section 3. By using expansions in spherical harmonics, we are able the reduce the conical Radon transform to a set of explicitly given one-dimensional integral equations of the Abel type (see Theorem 3.2). The invertibility of the one dimensional integral operators will be given in Theorem 3.5. For that purpose, in Appendix A we derive uniqueness results for first kind Volterra equations (Theorem A.2) and generalized Abel equations with kernels having zeros on the diagonal (Theorem 3.4). Theorem 3.5 in particular implies injectivity of the considered transform and additionally yields an efficient inversion method. In Section 4, we develop such a reconstruction procedure based on our theoretical findings and present some numerical results. Finally, in Section 5 we present a short summary and discuss possible lines of future research.

We start this section with defining the conical Radon transform that integrates a function over cones with vertices on the unit sphere and central axis orthogonal to . For and , we denote by

 C(z,ψ)={z+rω∣r≥0 and ω∈Sn−1 % with −ω\vbox\scalebox{.6}{∙}z=cosψ},

the surface of a right circular half cone in with vertex , central axis and half opening angle . We denote by the set of all infinitely times differentiable functions with , where denotes the unit ball in . Further, we denote by the set of all orthogonal matrices and set .

Definition 2.1 (The conical Radon transform Rmf).

Let . We define the conical Radon transform (with vertices on the sphere, orthogonal axis and weighting factor ) of by

 Rmf:Sn−1×(0,π/2)→R:(z,ψ)↦∫C(z,ψ)f(x)∥x−z∥mdS(x). (2.1)

The problem under study is recovering the function from its conical Radon transform . We start by deriving explicit expressions for .

Lemma 2.2.

Let and .

1. If and , then .

2. For every , we have

 (Rmf)(e1,ψ) =∫20rm(rsin(ψ))n−2 (2.2) ×∫Sn−2f(1−rcos(ψ),rsin(ψ)η)dS(η)dr, (Rmf)(e1,ψ) =∫π−ψ0(sin(ψ))n−1(sin(α))m+n−2(sin(α+ψ))m+n (2.3) ×∫Sn−2f(sin(ψ)sin(α+ψ)(cos(α),sin(α)η))dS(η)dα.
Proof.

1 For every and every , we have

 (Rmf)(Qz,ψ) =∫Q(C(z,ψ)) f(x)∥x−Qz∥mdS(x) =∫C(z,ψ) f(Qx)∥Qx−Qz∥mdS(x) =∫C(z,ψ) (f∘Q)(x)∥x−z∥mdS(x) =Rm(f∘Q)(z,ψ).

2 Let be any parametrization of , where is an open subset of . Then

 Ψ:D×(0,∞)→Rn:(r,β)↦(1−rcos(ψ),rsin(ψ)Φ(β))

is a parametrization of . Elementary computation shows that the Gramian determinant of is given by . Consequently,

 (Rmf)(e1,ψ) =∫C(e1,ψ)f(x)∥x−e1∥mdS(x) =∫∞0rm(rsin(ψ))n−2∫Df(1−rcos(ψ),rsin(ψ)Φ(β)) ×√det(Φ′(β)TΦ′(β))dβdr =∫∞0rm(rsin(ψ))n−2∫Sn−2f(1−rcos(ψ),rsin(ψ)η)dS(η)dr,

which is (2.2). Substituting , we have and ; this yields (2.3). ∎

Next we state the continuity of with respect to the -norms for . Similar results could of course be obtained for any .

Lemma 2.3 (Continuity of Rm).

Let , and .

1. If , then .

2. If , then .

3. If , then , for constants and independent of .

Proof.

1 Let and satisfy . By Lemma 2.2, we have

Using the Cauchy-Schwarz inequality, we obtain

 ∥(Rmf)(z,⋅)∥2L2≤|Sn−2|(∫20r2m+n−3dr)(∫π/20(sin(ψ))2(n−2)×∫20∫Sn−2rn−1|(f∘Q)(1−rcos(ψ),rsin(ψ)η)|2dS(η)drdψ).

The first integral equals , and the second can be bounded by . Consequently, . Integration over yields the claimed estimate.

2, 3: Analogous to 1. ∎

3 Analytic inversion of Rm

In this section, first we derive an explicit decomposition of the conical Radon transform in one-dimensional integral operators (see Theorem 3.2). Second, we show the solution uniqueness of the corresponding generalized Abel equations (see Theorem 3.5), which implies the invertibility of . For these results we will use the spherical harmonic decompositions

 f(rθ) =∞∑ℓ=0N(n,ℓ)∑k=1fℓ,k(r)Yℓ,k(θ), (3.1) (Rmf)(z,ψ) =∞∑ℓ=0N(n,ℓ)∑k=1(Rmf)ℓ,k(ψ)Yℓ,k(z). (3.2)

Here , for and , denote spherical harmonics [30, 36] of degree forming a complete orthonormal system in . The set of all with and will be denoted by .

3.1 Integral equations for fℓ,k

Let denote the Gegenbauer polynomials normalized in such a way that . We derive three different relations between and . The first one is as follows.

Lemma 3.1.

Let , and let and for be as in (3.1) and (3.2). Then

 ∀ψ∈(0,π/2):(Rmf)ℓ,k(ψ)=|Sn−2|∫π−ψ0fℓ,k(sin(ψ)sin(α+ψ))×(sin(ψ))n−1(sin(α))m+n−2(sin(α+ψ))m+nC(n−2)/2ℓ(cos(α))dα. (3.3)
Proof.

Fix and let be any rotation with . Using the delta distribution and applying the Funk-Hecke theorem, for any we have

 ∫Sn−2Yℓ,k(Q(cos(α),sin(α)η))dS(η)=∫Sn−1Yℓ,k(Qη)δ(e1\vbox\scalebox{.% 6}{∙}η−cos(α))(1−(e1\vbox\scalebox{.6}{∙}η)2)−(n−3)/2dS(η)=∫Sn−1Yℓ,k(η)δ(z\vbox\scalebox{.6}{∙}η−cos(α))(1−(z\vbox\scalebox{.6}{∙}η)2)−(n−3)/2dS(η)=|Sn−2|Yℓ,k(z)∫1−1δ(t−cos(α))C(n−2)/2ℓ(t)dt=|Sn−2|Yℓ,k(z)C(n−2)/2ℓ(cos(α)). (3.4)

Together with Lemma 2.2, this yields

 Rm[x↦fℓ,k(|x|)Yℓ,k(x/|x|)](z,ψ)=|Sn−2|(∫π−ψ0fℓ,k(sin(ψ)sin(α+ψ))×(sin(ψ))n−1(sin(α))m+n−2(sin(α+ψ))m+nC(n−2)/2ℓ(cos(α))dα)Yℓ,k(z).

The linearity of gives (3.3). ∎

Theorem 3.2 (Generalized Abel equation for fℓ,k).

Let and let and be as (3.1) and (3.2) for . Then, for ,

 (Rmf)ℓ,k(ψ)=|Sn−2|sin(ψ)−m∫1sin(ψ)fℓ,k(ρ)ρKℓ(ψ,ρ)√ρ2−(sin(ψ))2dρ, (3.5)

with the kernel functions

 Kℓ(ψ,ρ)\coloneqqρm+n−2∑σ=±1σℓsin(arcsin(sin(ψ)/ρ)−σψ)m+n−2×C(n−2)/2ℓ(cos(arcsin(sin(ψ)/ρ)−σψ)). (3.6)
Proof.

We split the integral in Lemma 3.1 in one integral over and one over . For we substitute . We have and therefore

In the case , we substitute . Repeating the above computations and using shows

 ∫π−ψπ/2−ψfℓ,k(sin(ψ)sin(α+ψ))(sin(ψ))n−1(sin(α))m+n−2(sin(α+ψ))m+nC(n−2)/2ℓ(cos(α))dα=(−1)ℓ(sin(ψ))−m∫1sinψfℓ,k(ρ)C(n−2)/2ℓ(cos(arcsin(sin(ψ)/ρ)+ψ))×(sin(arcsin(sin(ψ)/ρ)+ψ))m+n−2ρm+n−1dρ√ρ2−(sin(ψ))2.

Together with (3.3), this yields the claim. ∎

The relation between and given in Theorem 3.2 is well suited for the numerical implementation, see Section 4. For showing uniqueness of a solution, the following equivalent form will be more appropriate.

Lemma 3.3.

Let and let and be as (3.1) and (3.2). Further, for every denote

1. ;

2. ;

3. .

Then and are related via:

 ∀t∈[0,1]:^gℓ,k(t)=∫t0^fℓ,k(s)Fℓ(t,s)√t−sds. (3.7)
Proof.

Substituting in (3.5) and using the trigonometric sum and difference identities shows

 1|Sn−2|(Rmf)ℓ,k(arcsin(w))=w−m∫1wfℓ,k(ρ)ρm+n−2∑σ=±1σℓsin(arcsin(w/ρ)−σarcsin(w))m+n−2×C(n−2)/2ℓ(cos(arcsin(w/ρ)−σarcsin(w)))ρdρ√ρ2−w2=w−m∫1wfℓ,k(ρ)ρm+n−2∑σ=±1σℓ(w/ρ√1−w2−σw√1−w2/ρ2)m+n−2×C(n−2)/2ℓ(√1−w2/ρ2√1−w2+σw2/ρ)ρdρ√ρ2−w2=wn−2∫1wfℓ,k(ρ)∑σ=±1σℓ(√1−w2−σ√ρ2−w2)m+n−2×C(n−2)/2ℓ(√ρ2−w2√1−w2+σw2ρ)ρdρ√ρ2−w2

Next we set and make the substitution . Then we have , and , which shows

 (1−t)−(n−2)/2|Sn−2|(Rmf)ℓ,k(arccos(√t))=12∫t0fℓ,k(√1−s)×∑σ=±1σℓ(√t−σ√t−s)m+n−2C(n−2)/2ℓ(√t√t−s+σ(1−t)√1−s)ds√t−s.

This together with 1-3 yields (3.7). ∎

3.2 Solution uniqueness

Any of the integral equations (3.7) is of generalized Abel type. Using the symmetry of the Gegenbauer polynomials, we see that . Since the Gegenbauer polynomials have zeros in , so has the function . Consequently, standard theorems on well-posedness do not apply to (3.7), because such results require a non-vanishing diagonal.

To investigate unique solvability of (3.7) (and, as consequence, of (3.5)), we derive a uniqueness result for generalized Abel equations of the form

 ∀t∈[a,b]:∫taF(t,s)√t−sf(s)ds=g(t), (3.8)

where corresponds to given data and , with , is a continuous kernel.

Theorem 3.4 (Solution uniqueness of Abel equations with kernel having zeros on the diagonal).

Suppose that , where , satisfies the following:

1. .

2. is finite and consists of simple roots.

3. For every , the gradient satisfies

 1+12β1β1+β2>0. (3.9)

Then, for any , equation (3.8) has at most one solution .

Proof.

See Appendix A. ∎

To the best of our knowledge, Theorem 3.4 is new; we are not aware of similar results for generalized Abel equations with zeros in the diagonal of the kernel. We derive this result by exploiting a well-posedness theorem due to Volterra and Pérès for first kind Volterra equations (see Lemma A.1) together with a standard procedure of reducing generalized Abel equations to Volterra integral equations of the first kind. We now apply Theorem 3.4 to the integral equation (3.7):

Theorem 3.5 (Uniqueness of recovering fℓ,k).

Suppose . For any and any , the spherical harmonic coefficient of can be recovered as the unique solution of

 ∀ψ∈(0,π/2):(Rmf)ℓ,k(ψ)=|Sn−2|sin(ψ)−m∫1sin(ψ)fℓ,k(ρ)ρKℓ(ψ,ρ)dρ√ρ2−(sin(ψ))2,

with the kernel functions defined by (3.6).

Proof.

Let vanish outside a ball of Radius . According to Lemma 3.3, it is sufficient to show that (3.7) has a unique solution. To show that this is indeed the case, we apply Theorem 3.4 by verifying that satisfies conditions 1-3.

Ad 1: Using the abbreviations and , the kernel can be written in the form

 ∀(t,s)∈Δ(a,1):Fℓ(t,s)=∑σ=±1σℓ(√t−σ√t−s)qC(√t√t−s+σ(1−t)√1−s).

From this expression it is clear that is smooth on . Further, by using one sees that is an even polynomial in . This shows that is also smooth on the diagonal .

Ad 2: Next, consider the restriction of the kernel to the diagonal. As an orthogonal polynomial, has a finite number of isolated and simple roots. We conclude that the same holds true for .

Ad 3: Let be a zero of and set . Then

 β1+β2=v′(s0)=−sq/20√1−s0C′(√1−s0). (3.10)

Next we compute . We have

 Fℓ(s0,s0−ϵ) =∑σ=±1σℓ(√s0−σ√ϵ)qC(√s0√ϵ+σ(1−s0)√1−s0+ϵ) =∑σ=±1σℓ(sq/20−σqs(q−1)/20√ϵ+q(q−1)2s(q−2)/20ϵ) ×(C′(σ√1−s0)√s0√1−s0√ϵ+(C′′(σ√1−s0)s02(1−s0) −σ√1−s0C′(σ√1−s0))ϵ)+O(ϵ2) =sq/20√1−s0(−(2q+1)C′(√1−s0)+s0C′′(√1−s0)√1−s0)ϵ+O(ϵ2).

Here for the last equality we used the symmetry properties and for the first and second derivatives of the Gegenbauer polynomials. Because is a solution of the differential equation

 (1−x2)C′′(x)−(n−1)xC′(x)+ℓ(ℓ+n−2)C(x)=0

and is a zero of , we have the identity . We conclude that . Together with (3.10) we obtain

 β1=(−2q+n−3)sq/20√1−s0C′(√1−s0). (3.11)

From (3.10) and (3.11) it follows that

 1+β12(β1+β2)=1+2q−n+32=m+n+12>0.

This shows 3. Consequently, Theorem 3.4 implies that is the unique solution of the integral equation (3.3). ∎

Theorem 3.5 immediately implies the following uniqueness result for the conical Radon transform .

Corollary 3.6 (Invertibility of Rm).

Suppose . If are such that , then .

Proof.

Let satisfy for all . According to Theorem 3.5, the integral equation (3.5) has the unique solution , which implies . The linearity of gives the claim. ∎

4 Numerical implementation

Theorems 3.2 and 3.5 are the basis of the following inversion method for the conical Radon transform :

1. Compute the expansion coefficients in (3.2).

2. Recover from by solving (3.5).

3. Compute .

In this section, we show how to implement this reconstruction procedure. We restrict ourselves to two spatial dimensions () and the case ; extensions to general cases are straightforward.

4.1 Basic procedure for numerically inverting the conical Radon transform

In two spatial dimensions, the conical Radon transform with  can be written in the form

 (Rf)(φ,ψ)\coloneqq∑σ=±1∫∞0f((cos(φ),sin(φ))−r(cos(φ−σψ),sin(φ−σψ)))dr. (4.1)

Because consists of integrals of over V-shaped lines, the 2D version is also known as the V-line Radon transform. In the 2D situation, the spherical harmonics expansion equals the common Fourier series expansion, and we obtain the following reconstruction procedure:

Algorithm 1 (Series expansion for inverting the V-line transform).

Goal: Recover from the V-line transform .

1. Compute