Inversion of the attenuated Radon transform

F. Natterer

Institut für Numerische und instrumentelle Mathematik Westf. Wilhelms-Universität Münster Einsteinstrasse 62, D-48149 Münster, Germany e-mail: nattere@math.uni-muenster.de

Abstract We derive an exact inversion formula for the attenuated Radon transform. The formula is closely related to Novikov’s inversion formula [5], but our derivation is completely different.

1 Introduction

Let a be a sufficiently smooth real valued function in IR^{2} which is decaying sufficiently fast at infinity. We put for x ∈ IR^{2}, θ ∈ S^{1}

(Da)(x, θ) = ∫_{0}^{∞}a(x + tθ)dt .  (1.1)

The attenuated Radon transform R_{a} is defined by

(R_{a}f)(θ, s) = ∫_{x·θ=s} e^{−(Da)(x,θ^{⊥})}f(x)dx .  (1.2)

Here, dx stands for the restriction of the Lebesgue measure in IR^{2} to x · θ = s and

θ^{⊥} = ( \begin{matrix}cos(φ + π/2) \\ sin(φ + π/2) \end{matrix} ) = ( \begin{matrix}− sin φ \\ cos φ \end{matrix} ) for θ = ( \begin{matrix}cos φ \\ sin φ \end{matrix} ) .  (1.3)

For a = 0, R_{0} is the usual Radon transform R. R_{a} is the relevant integral transform in single photon emission tomography (SPECT); see [2]. In emission tomography, f is the activity distribution of a radiopharmaceutical inside the body and a is the attenuation of the tissue. The problem is to recover f from g = R_{a}f, the attenuation map a assumed to be known. For a = 0, the problem is solved by the Radon inversion formula. For a constant and known inside a convex set the problem can be reduced to the exponential Radon transform for which various inversion formulas are available [7], [6].

The purpose of this paper is an explicit inversion formula for R_{a}. Our formula is different but similar to the one given in [5]. Our derivation is based on tools developed in the study of the range of R_{a} in [4]. Our formula has the same structure as Radon’s inversion formula for R and leads immediately to a reconstruction algorithm of the filtered backprojection type for SPECT.

2 Derivation of the formula

The function

h = \frac{1}{2}(I + iH)Ra  (2.1)

where R is the Radon transform and H is the Hilbert transform, i.e.

(Hg)(s) = \frac{1}{π} ∫_{R^{1}} \frac{g(t)}{s − t}dt ,

acting on the second variable of Ra, plays a central role in our and in Novikov’s inversion formula. It also is an essential ingredient of the consistency conditions in the range of R_{a}.

Lemma 2.1 Let θ = ( \begin{matrix}cos φ \\ sin φ \end{matrix} ) and u(x, θ) = h(θ, x·θ)−(Da)(x, θ^{⊥}). Then,

with certain functions u_{ℓ}(x),

u(x, θ) = ∑_{ℓ>0 odd} u_{ℓ}(x)e^{iℓφ} .

For the proof see the proof of Theorem I.6.2 in [3].

Lemma 2.2 Let θ = ( \begin{matrix}cos φ \\ sin φ \end{matrix} ) and \frac{x^{⊥}}{|x|} = ( \begin{matrix}cos ψ \\ sin ψ \end{matrix} ). Then,

∫_{0}^{2π} \frac{θ}{x · θ}e^{iℓφ}dφ = \{ \begin{matrix}0 & , & ℓ  odd , \\ 2πx/|x|^{2} & , & ℓ = 0  , \\ −2πie^{iℓψ}x^{⊥}/|x|^{2} & , & ℓ > 0 even . \end{matrix}

The integral is to be understood as Principle Value.

Proof: We have for each integer ℓ

∫_{0}^{2π} \frac{θ}{x · θ}e^{iℓφ}dφ = −\frac{1}{|x|} ∫_{0}^{2π}( \begin{matrix}cos φ \\ sin φ \end{matrix} ) \frac{e^{iℓφ}dφ}{sin(φ − ψ)}

= −\frac{1}{|x|}e^{iℓψ} ∫_{0}^{2π}( \begin{matrix}cos(φ + ψ) \\ sin(φ + ψ) \end{matrix} ) \frac{e^{iℓφ}dφ}{sin φ}

= −\frac{1}{|x|}e^{iℓψ} ∫_{0}^{2π}( \begin{matrix}cos φ cos ψ − sin φ sin ψ \\ sin φ cos ψ + cos φ sin ψ \end{matrix} ) \frac{e^{iℓφ}dφ}{sin φ} .

For ℓ = 0 we obtain immediately

∫_{0}^{2π} \frac{θ}{x · θ}dφ = −\frac{2π}{|x|} ( \begin{matrix}− sin ψ \\ cos ψ \end{matrix} ) = \frac{2π}{|x|} \frac{x}{|x|} .

For ℓ ≠ 0 we have

∫_{0}^{2π} \frac{θ}{x · θ}e^{iℓφ}dφ = −\frac{1}{|x|}e^{iℓψ} ( \begin{matrix}cos ψ \\ sin ψ \end{matrix} ) ∫_{0}^{2π}cos φ\frac{e^{iℓφ}}{sin φ}dφ .  (2.2)

With the help of

∫_{0}^{2π} \frac{sin(ℓφ)}{sin φ}dφ = 2π , ℓ > 0 odd

(see [1], formula 3.612) we recognize the integrals in (2.2) as 2πi for ℓ > 0 even. Hence, for ℓ > 0 even,

∫_{0}^{2π} \frac{θ}{x · θ}e^{iℓφ}dφ = −\frac{2πi}{|x|}e^{iℓψ} \frac{x^{⊥}}{|x|} .

□

Lemma 2.3 With δ the Dirac function we have

div \frac{x}{|x|^{2}} = 2πδ(x) .

Proof: We have to show that for f ∈ C_{0}^{∞}(IR^{2})

∫_{IR^{2}} \frac{x}{|x|^{2}} · ∇f(x)dx = −2πf(0) .  (2.3)

For each smooth vector field g in |x| ≥ ε we have

∫_{|x|≥ε} g · ∇fdx = − ∫_{|x|≥ε} f div gdx + ∫_{|x|=ε} fg · νds

where ν is the exterior normal on |x| = ε of |x| ≥ ε, i.e. ν = −x/|x|, and ds the surface measure on |x| = ε. Putting g = x/|x|^{2} we have div g = 0 for x ≠ 0, hence

∫_{|x|≥ε} \frac{x}{|x|^{2}} · ∇f(x)dx = − ∫_{|x|=ε} f \frac{x}{|x|^{2}} · \frac{x}{|x|}ds

= −\frac{1}{ε} ∫_{|x|=ε} f(x)ds .

(2.3) follows by letting ε → 0.

□

Now we can prove our inversion formula.

Theorem 2.1 Let g = R_{a}f and h the function (2.1). Then,

f(x) = \frac{1}{4π} Re div ∫_{S^{1}} θe^{(Da)(x,θ^{⊥})}(e^{−h}He^{h}g)(θ, x · θ)dθ .

Proof: It sufficies to prove the theorem for f(x) = δ(x−y). With this choice of f,

(He^{h}g)(θ, s) = \frac{1}{π} ∫_{IR^{1}} \frac{1}{s − t}e^{h(θ,t)} ∫_{x·θ=t} e^{−(Da)(x,θ^{⊥})}δ(x − y)dxdt

= \frac{1}{π} \frac{1}{s − θ · y} e^{h(θ,θ·y)−(Da)(y,θ^{⊥})}
= \frac{1}{π} \frac{1}{s − θ · y} e^{u(y,θ)}

with u as in Lemma 2.1. Using this it remains to show that

δ(x − y) = Re\frac{1}{4π^{2}} div ∫_{S^{1}} \frac{θ}{(x − y) · θ}e^{u(y,θ)−u(x,θ)}dθ .  (2.4)

Putting θ = ( \begin{matrix}cos φ \\ sin φ \end{matrix} ) we have from Lemma 2.1

u(y, θ) − u(x, θ) = ∑_{ℓ>0 odd}(u_{ℓ}(y) − u_{ℓ}(x))e^{iℓφ} .

It follows that with certain functions u_{ℓ}(x, y)

cosh(u(y, θ) − u(x, θ)) = 1 + ∑_{ℓ>0 even} u_{ℓ}(x, y)e^{iℓφ} ,
sinh(u(y, θ) − u(x, θ)) = ∑_{ℓ>0 odd} u_{ℓ}(x, y)e^{iℓφ} .

From Lemma 2.2 we obtain with \frac{(x−y)^{⊥}}{|x−y|} = ( \begin{matrix}cos ψ \\ sin ψ \end{matrix} ) = ω

∫_{0}^{2π} \frac{θ}{(x − y) · θ} cosh (u(y, θ) − u(x, θ)) dφ

= 2π \frac{x − y}{|x − y|^{2}} − 2πi ∑_{ℓ>0even} u_{ℓ}(x, y)e^{iℓψ} \frac{(x − y)^{⊥}}{|x − y|^{2}}

= 2π \frac{x − y}{|x − y|^{2}} − 2πi (cosh(u(y, ω) − u(x, ω)) − 1) \frac{(x − y)^{⊥}}{|x − y|^{2}} ,

∫_{0}^{2π} \frac{θ}{(x − y) · θ} sinh (u(y, θ) − u(x, θ)) dφ = 0 .

Since x · ω = y · ω we have

u(y, ω) − u(x, ω) = −(Da)(y, ω^{⊥}) + (Da)(x, ω^{⊥})
= − ∫_{y}^{x} ads

with integration along the straight line joining y and x. This is real. Hence,

Re ∫_{0}^{2π} \frac{θ}{(x − y) · θ}e^{u(y,θ)−u(x,θ)}dφ = 2π \frac{x − y}{|x − y|^{2}} .

Now (2.4) follows from Lemma 2.3.

□

We remark that in our notation Novikov’s formula reads

f(x) = −\frac{1}{4π} Re div ∫_{S^{1}} θe^{−(Da)(x,θ^{⊥})} (e^{h}He\overline{h}ǧ) (θ, x · θ)dθ

with ǧ(θ, s) = (R_{a}f)(−θ, −s). Using Lemma 2.1 it is easy to show that Novikov’s formula implies ours. Note also that for vanishing attenuation, i.e. a = 0, both formulas reduce to

f(x) = \frac{1}{4π} div ∫_{S^{1}} θ(Hg)(θ, x · θ)dθ

= \frac{1}{4π} ∫_{S^{1}} (Hg′)(θ, x · θ)dθ .

This is Radon’s inversion formula.

3 Implementation

The implementation of the inversion formula is an immediate extension of the familiar filtered backprojection algorithm of computerized tomography [2], [4]: First we have to compute the function

g_{a} = Re e^{−h}He^{h}g , h = \frac{1}{2}(I + iH)Ra .

Putting h = h_{1} + ih_{2}, h_{1} = \frac{1}{2}Ra, h_{2} = \frac{1}{2}HRa we have

g_{a} = e^{−h_{1}} (cos h_{2} He^{h_{1}} cos h_{2} + sin h_{2} He^{h_{1}} sin h_{2}) g .  (3.1)

The Hilbert transform H can be computed by convolution: With

ĝ(σ) = (2π)^{−1/2} ∫_{IR^{1}} e^{−isσ}g(s)ds

the Fourier transform of g we have

(Hg)^{∧}(σ) = \frac{sgn(σ)}{i}ĝ(σ) .

Now let ϕ be a low-pass filter, i.e. ϕ(σ) ∼ 1 for |σ| small and ϕ(σ) = 0 for |σ| > 1, and let Ω > 0 be a bandwidth. We take

(H_{Ω}g)^{∧}(σ) = ϕ(σ/Ω)(Hg)^{∧}(σ)

as an approximation to H. It follows that

H_{Ω}g = v_{Ω} ∗ g ,

(3.2)

v_{Ω}(s) = \frac{1}{2π} ∫_{−Ω}^{Ω} ϕ(σ/Ω)\frac{sgn(σ)}{i}e^{isσ}dσ .

For ϕ the ideal low-pass; i.e. ϕ(σ) = 1 for |σ| ≤ 1 and ϕ(σ) = 0 for |σ| > 1 we obtain

v_{Ω}(s) = \frac{Ω}{π}u(Ωs) , u(s) = (1 − cos s)/s .  (3.3)

Thus (3.1) can be computed for each direction θ by pre-multiplying the emission data g, one-dimensional filtering with v_{Ω} and post-multiplying. v_{Ω} corresponds to the ramp filter in the usual filtered backprojection algorithm.

Once g_{a} is computed we have to carry out the weighted backprojection

f(x) = \frac{1}{4π} div ∫_{S^{1}} θe^{(Da)(x,θ^{⊥})}g_{a}(θ, x · θ)dθ  (3.4)

for each reconstruction point x, followed by applying the div operator. The weighted backprojection can be done exactly as in the usual filtered backprojection algorithm by using the trapezoidal rule, with linear interpolation in the second argument of g_{a}. The div operator can easily be applied by taking finite differences.

Alternatively one can carry out the div operator under the integral sign, leading to the inversion formula

f(x) = \frac{1}{4π} ∫_{S^{1}} e^{(Da)(x,θ^{⊥})}g′_{a}(θ, x · θ)dθ

(3.5)

+\frac{1}{4π} ∫_{S^{1}} \frac{∂}{∂θ}(Da)(x, θ^{⊥})e^{(Da)(x,θ^{⊥})}g_{a}(θ, x · θ)dθ

with \frac{∂}{∂θ} denoting the directional derivative in direction θ with respect to the first argument of Da, and the prime denoting the derivative with respect to the second argument of g_{a}. Both derivatives can easily be computed by finite differences.

We implemented both versions of the algorithm. The difference in accuracy and efficiency is negligible. The complexity of the algorithm is the same as for the usual filtered backprojection algorithm: For the reconstruction on a q × q grid from p projections one needs O(pq^{2}) operations.

For the choice of p, q, Ω we found that they are subject to the same restrictions as in the Radon case [3], which we have to adjust to 360^{∘} scanning: For artefact-free reconstruction of a function of essential bandwidth Ω we have to have p ≥ 2Ω, Ω ≤ πq/ρ where |x| < ρ is the reconstruction region.

References

[1] Gradshteyn, I.S. and Ryzhik, I.M.: Table of Integrals, Series, and Products. Academic Press 1965.

[2] Kak, A.C. and Slaney, M.: Principles of Computerized Tomography. IEEE Press 1987.

[3] Natterer, F.: The Mathematics of Computerized Tomography. Wiley-Teubner 1986.

[4] Natterer, F.: The ill-posedness of Radon’s integral equation, Symposion on Ill-Posed Problems: Theory and Practice, Delaware, October 2-6, 1979.

[5] Novikov, R.G.: An inversion formula for the attenuated X-ray transform, Preprint, Département de Mathématique, Nantes, April 2000.

[6] Palamodov, V.P.: An inversion formula for an attenuated X-ray transform, Inverse Problems 12, 717-730 (1996).

[7] Tretiak, O. and Metz, C.: The exponential Radon transform, SIAM J. Appl. Math. 39, 341-354 (1980).