Reconstruction algorithm for single photon emission computed tomography and its numerical implementation

A.S Fokas, A Iserles, V Marinakis

Abstract

The modern imaging techniques of positron emission tomography and of single photon emission computed tomography are not only two of the most important tools for studying the functional characteristics of the brain, but they now also play a vital role in several areas of clinical medicine, including neurology, oncology and cardiology. The basic mathematical problems associated with these techniques are the construction of the inverse of the Radon transform and of the inverse of the so-called attenuated Radon transform, respectively. An exact formula for the inverse Radon transform is well known, whereas that for the inverse attenuated Radon transform was obtained only recently by R. Novikov. The latter formula was constructed by using a method introduced earlier by R. Novikov and the first author in connection with a novel derivation of the inverse Radon transform. Here, we first show that the appropriate use of that earlier result yields immediately an analytic formula for the inverse attenuated Radon transform. We then present an algorithm for the numerical implementation of this analytic formula, based on approximating the given data in terms of cubic splines. Several numerical tests are presented which suggest that our algorithm is capable of producing accurate reconstruction for realistic phantoms such as the well-known Shepp–Logan phantom.

Keywords:

1. Introduction

Positron emission tomography (PET) and single photon emission computed tomography (SPECT) are two modern imaging techniques with a wide range of medical applications. Although these techniques were originally developed for the study of the functional characteristics of the brain, they are now used in many diverse areas of clinical medicine. For example, a recent editorial in the New England Journal of Medicine (Koh et al. 2003) emphasized the importance of PET in oncologic imaging. Other medical applications of PET and SPECT are presented in Lauritzen & Olesen (1984), Vorstrup et al. (1986), Mazziotta et al. (1987), Lee et al. (1988), Junck et al. (1989), Reiman et al. (1989), Mazziotta (1992), Tyler & Byme (1992), Minoshima et al. (1995), Jonides et al. (1996), Mark et al. (1996), Andreasen (1997), Hutton (1997), Tjuvajev et al. (1999), Wackers (1999), Yu et al. (2000), Bengel et al. (2001), Doubrovin et al. (2001), Green et al. (2002), Lardinois et al. (2003) and Ost et al. (2003).

The first step in PET is to inject the patient with a dose of a suitable radiopharmaceutical. For example, in brain imaging a typical such radiopharmaceutical is flurodeoxyglucose (FDG), which is a normal molecule of glucose attached artificially to an atom of radioactive fluorine. The cells in the brain, which are more active have a higher metabolism, need more energy, thus will absorb more FDG. The fluorine atom in the FDG molecule suffers a radioactive decay, emitting a positron. When a positron collides with an electron it liberates energy in the form of two beams of gamma rays travelling in opposite direction, which are picked by the PET scanner. In SPECT one uses photon emitting radiopharmaceuticals.

Let x1 and x2 be the Cartesian coordinates in the plane. In both PET and SPECT the radiating sources are inside the body, and the aim is to determine the distribution g(x1,x2) on a given plane, of the relevant radiopharmaceutical from measurements made outside the body of the emitted radiation. By using a sufficiently large number of planes, the three dimensional distribution can be determined. If f(x1,x2) is the X-ray attenuation coefficient of the body, then it is straightforward to show (Natterer 1986) that the intensity I outside the body measured by a detector which picks up only radiation along the straight line L is given byEmbedded Image(1.1)where τ is a parameter along L, and L(x) denotes the section of L between the point (x1,x2) and the detector. The attenuation coefficient f(x1,x2) is precisely the function measured by the usual computed tomography. Thus, the basic mathematical problem in SPECT is to determine the function g(x1,x2) from the knowledge of the ‘transmission’ function f(x1,x2) (determined via computed tomography) and the ‘emission’ function I (known from the measurements).

In PET the situation is simpler. Indeed, since the sources eject particles pairwise in opposite directions and the radiation in opposite directions is measured simultaneously, equation (1.1) is replaced byEmbedded Image(1.2)where L+, L are the two half-lines of L with endpoint x. Since, L++L=L, equation (1.2) becomesEmbedded ImageWe recall that the line integral of the function f(x1,x2) along L is precisely what is known from the measurements in the usual computed tomography. Thus, since both I and the integral of f(x1,x2) are known (from the measurements of PET and of computed tomography, respectively), the basic mathematical problem of PET is to determine g(x1,x2) from the knowledge of its line integrals. This mathematical problem is identical with the basic mathematical problem of computed tomography.

1.1 Notation

  1. A point of a line L making an angle θ with the x1-axis is specified by the three real numbers (τ,ρ,θ), where τ is a parameter along L, −∞<τ<∞, ρ is the distance from the origin to the line, −∞<ρ<∞, and 0≤θ≤2π (figure 1).

    Figure 1

    Local coordinates for the mathematical formulation of PET and SPECT.

  2. The above parameterization implies that, for a fixed θ, the Cartesian coordinates (x1,x2) can be expressed in terms of the local coordinates (τ, ρ) by the equations (see §2)Embedded Image(1.3)A function f(x1,x2) rewritten in local coordinates will be denoted by F(τ,ρ,θ),Embedded ImageThus, F(τ,ρ,θ) and G(τ,ρ,θ) will denote the X-ray attenuation coefficient f(x1,x2) and the distribution of the radiopharmaceutical g(x1,x2), rewritten in local coordinates.

  3. The line integral of a function f is called its Radon transform and will be denoted by Embedded Image. In order to compute Embedded Image, we first write f in local coordinates and then integrate with respect to τ,Embedded Image(1.4)The line integral of the function g with respect to the weight f appearing in equation (1.1) is called the attenuated Radon transform of g (with the attenuation specified by f) and will be denoted by Embedded Image. In order to compute Embedded Image, we write both g and f in local coordinates and then evaluate the following integral:Embedded Image(1.5)

1.2 Mathematical methods

The basic mathematical problem of both computed tomography and PET is to reconstruct a function f from the knowledge of its Radon transform Embedded Image, i.e. to solve equation (1.4) for f(x1,x2) in terms of Embedded Image. The relevant formula is called the inverse Radon transform and is given byEmbedded Image(1.6)where −∞<xj<∞, j=1, 2 and ∮ denotes principal value integral.

A novel approach for deriving equation (1.6) was introduced in Fokas & Novikov (1991), and is based on the analysis of the equationEmbedded Image(1.7)where λ is a complex parameter different than zero. The application of this approach to a slight generalization of equation (1.7) can be used to reconstruct a function g from the knowledge of its attenuated Radon transform Embedded Image, i.e. this approach can be used to solve equation (1.5) for g(x1,x2) in terms of Embedded Image and f(x1,x2). The relevant formula, called the inverse attenuated Radon transform, was obtained by Novikov (2002) by analysing, instead of equation (1.7), the equationEmbedded Image(1.8)

1.3 Organization of the paper

In §2 we first review the analysis of equation (1.7), and then show that if one uses the basic result obtained in this analysis, it is possible to construct immediately the inverse attenuated Radon transform. This provides a dramatic simplification of the derivation of Novikov (2002). In §3 we present a new numerical reconstruction algorithm for SPECT. This algorithm is based on approximating the given data in terms of cubic splines (a similar reconstruction algorithm for PET is given in appendix A). We recall that both the exact inverse Radon transform as well as the exact inverse attenuated Radon transform involve the Hilbert transform of the data functions. For example, the inverse Radon transform involves the functionEmbedded Image(1.9)Existing numerical approaches use the convolution property of the Fourier transform to compute the Hilbert transform and employ appropriate filters to eliminate high frequencies. It appears that our approach has the advantage of simplifying considerably the mathematical formulas associated with these techniques. Furthermore, accurate reconstruction is achieved, for noiseless data, with the additional use of an averaging or of a median filter. Several numerical tests are presented in §4. One of these tests involves the Shepp–Logan phantom (Shepp & Logan 1974) (see figure 3c).

Figure 3

Test phantoms for the PET algorithm.

Numerical algorithms based on the filtered back projection are discussed in Kunyansky (2001), Natterer (2001), Guillement et al. (2002) and Guillement & Novikov (2004), while algorithms based on iterative techniques can be found in Hebert et al. (1988), Liang & Hart (1988) and Nuyts & Fessler (2003). Spline based algorithms for computed tomography are presented in La Rivière & Pan (1998).

2. Mathematical methods

We first review the basic result of Fokas & Novikov (1991). It will be shown later that using this result it is possible to derive both the inverse Radon as well as the inverse attenuated Radon transforms in a straightforward manner.

Define the complex variable z byEmbedded Image(2.1)where x1, x2 are the real Cartesian coordinates −∞<xj<∞, j=1, 2, and λ is a complex variable, λ≠0. Assume that the function f(x1,x2) has sufficient decay as |x1|+|x2|→∞. Let μ(x1, x2, λ) satisfy the equationEmbedded Image(2.2)as well as the boundary condition μ=O(1/z) as |x1|+|x2|→∞. Let λ+(θ) and λ(θ) denote the limits of λ as it approaches the unit circle (figure 2) from inside and outside the unit disc, respectively, i.e.Embedded ImageThenEmbedded Image(2.3)where Embedded Image denotes the Radon transform of f, F denotes f in the local coordinates (see §1.1), P± denote the usual projection operators in the variable ρ, i.e.Embedded Image(2.4)anddenotes the principal value integral.

Before deriving this result, we first note that equation (2.1) is a direct consequence of equation (1.7). Indeed, equation (1.7) motivates the introduction of the variable z defined by equation (2.1). Taking the complex conjugate of equation (2.1) we findEmbedded Image(2.5)where Embedded Image and Embedded Image denote the complex conjugates of z and λ, respectively. Equations (2.1) and (2.5) define a change of variables from (x1,x2) to Embedded Image. Using this change of variables to compute Embedded Image and Embedded Image in terms of ∂z and Embedded Image, equation (1.7) becomes (2.2).

We now derive equation (2.3). The derivation is based on the following two steps, which have been used extensively in the field of nonlinear integrable partial differential equations (see, for example, Fokas & Gel'fand 1994).

In the first step (sometimes called the direct problem), we consider equation (2.2) as an equation which defines μ in terms of f, and we construct an integral representation of μ in terms of f, for all complex values of λ. This representation isEmbedded Image(2.6)Indeed, suppose that the function μ(zR, zI) satisfies the equationEmbedded Imageas well as the boundary condition μ=O(1/z) as z→∞. Then, Pompieu's formula (see, for example, Ablowitz & Fokas 1997) impliesEmbedded Image(2.7)This equation is the direct consequence of the fact that since μ in equation (2.7) depends on z only through (z′−z)−1, then its dependence on Embedded Image can be computed explicitly using the formulaEmbedded Image(2.8)In our caseEmbedded Imagethus equation (2.7) becomes (2.6).

In the second step (sometimes called the inverse problem), we analyse the analyticity properties of μ with respect to λ, and we find an alternative representation for μ. This representation involves certain integrals of f called spectral functions. For our problem, this representation is equation (2.3). Indeed, since μ is an analytic function of λ for |λ|≠1 and since μ=O(1/λ) as λ→∞, we can reconstruct the function μ if we know its ‘jump’ across the unit circle:Embedded Image(2.9)whereEmbedded ImageThus, we need to compute the limits of μ as λ tends to λ±. As ϵ→0,Embedded ImageSubstituting this expression in the definition of z (equation (2.1)) and simplifying, we findEmbedded Image(2.10)The right-hand side of this equation can be rewritten in terms of the local coordinates ρ, ρ′, τ, τ′: Let k and k denote two unit vectors along the line L and perpendicular to this line, respectively. ThenEmbedded ImageorEmbedded ImageHence, x1 and x2 are given by equation (1.3). Inverting these equations we findEmbedded Image(2.11)Thus, equation (2.10) becomesEmbedded ImageSubstituting this expression in equation (2.6) and using the fact that the relevant sign equals 1, we findEmbedded Image(2.12)Using the change of variables (x1,x2)↔(τ, ρ) defined by equations (1.3) and (2.11), and noting that the relevant Jacobian is 1, i.e.Embedded Imagewe find that the right-hand side of equation (2.12) equalsEmbedded Image(2.13)In order to simplify this expression we split the integral over dτ′ in the formEmbedded Imageand note that in the first integral τ′−τ<0, while in the second integral τ′−τ>0. Thus, using the second set of equation (2.4) the expression in (2.13) becomesEmbedded ImageFinally, adding and subtracting the integral Embedded Image we findEmbedded ImageThe first two terms in the right-hand side of this equation equal Embedded Image, hence we find (2.3)+. The derivation of equation (2.3) is similar. ▪

Using equation (2.3) it is now straightforward to derive both the inverse Radon and the inverse attenuated Radon transforms. In this respect we note that the result of proposition 2.1 can be rewritten in the formEmbedded Image(2.14)whereEmbedded Image(2.15)

2.1 The inverse radon transform

Equation (2.3) yieldsEmbedded Image(2.16)Equation (2.9) impliesEmbedded ImageSubstituting this expression in equation (1.7) we findEmbedded Image(2.17)Replacing in this equation J by the right-hand side of equation (2.16) we find equation (1.6).

2.2 The attenuated radon transform

Equation (1.8) can be rewritten in the formEmbedded Imagewhere ν is defined by equation (2.15). Hence,Embedded ImageorEmbedded ImageReplacing in this equation Embedded Image by the right-hand side of equation (2.14) we findEmbedded ImageFor the computation of the right-hand side of this equation we use again equation (2.14), where f is replaced by g times the two exponentials appearing in the above relation. Hence,Embedded Image(2.18)Note that the term Embedded Image is independent of τ′; thus this term comes out of the integral Embedded Image, and furthermore the same term appears in the left-hand side of equation (2.18). Hence, when computing the jump μ(x1, x2, λ+)−μ(x1, x2, λ), the second term in the right-hand side of equation (2.18) cancels and we find that the relevant jump in now given byEmbedded Image(2.19)where τ and ρ are expressed in terms of x1 and x2 by equation (2.11).

Equation (2.9) is still valid, furthermore equation (2.17) is valid if f is replaced by g. Hence, replacing f by g in equation (2.17) we findEmbedded Image(2.20)where J is defined by equation (2.19). This formula is equivalent to Novikov's formula.

In summary, let Embedded Image be defined by equation (1.5), let F(τ,ρ,θ) denote the function f(x1,x2) written in local coordinates (see §1.1) and let Embedded Image denote the Radon transform of f(x1,x2) (see equation (1.4)). Then g(x1,x2) is given by equation (2.20) where the function J is explicitly given in terms of Embedded Image and Embedded Image by equation (2.19).

3. Reconstruction algorithm for SPECT

We denote the first exponential term of the right-hand side of equation (2.19) by I(τ,ρ,θ), i.e.Embedded Image(3.1)Note that, we assume compact support; thus the integration domain is finite, i.e. Embedded Image and F(τ,ρ,θ)=0 for |ρ|≥1, or for Embedded Image.

Definition (2.4) becomesEmbedded ImageMoreover,Embedded ImageEmbedded ImageWe introduce the following notation:Embedded Image(3.2)Embedded Image(3.3)Embedded Image(3.4)Using this notation and setting R(τ,ρ,θ)=−J(τ,ρ,θ), after some calculations, equation (2.19) becomesEmbedded Image(3.5)We now setEmbedded Imagethus equation (3.5) becomesEmbedded ImageWe denote the right-hand side of this equation by −ir(τ,ρ,θ). Taking the real part of g(x1,x2) in equation (2.20), we obtainEmbedded Image(3.6)where τ and ρ are given by equation (2.11) andEmbedded Image(3.7)For the numerical calculation of the Hilbert transform we suppose that Embedded Image is given, for every θ, at n equally spaced points ρi∈[−1,1], i.e. we suppose that Embedded Image are known. Moreover, in each interval [ρi,ρi+1] we approximate Embedded Image using the relationEmbedded Image(3.8)whereEmbedded Imageand Embedded Image denotes the second derivative of Embedded Image with respect to ρ, at ρ=ρi. In other words, we approximate Embedded Image by a cubic spline (in ρ) with equally-spaced nodes.

We then writeEmbedded Image(3.9)If ρ=ρi or ρ=ρi+1 the integral in the right-hand side of equation (3.9) can be written asEmbedded ImageThus, after some calculations, we obtainEmbedded Image(3.10)If ρρi and ρρi+1 the integral in the right-hand side of equation (3.9) can be written asEmbedded Imageand after some calculation we obtainEmbedded Image(3.11)where Fi is the right-hand side of equation (3.10).

Taking the real part of equation (1.6) it follows that f(x1,x2) is given byEmbedded Image(3.12)where h(ρ,θ) is defined by equation (1.9). We assume that f(x1,x2) has compact support, namely f(x1,x2)=0, for Embedded Image. Thus, in order to calculate numerically I(τ,ρ,θ) for any x1, x2, θ, we use relation (3.12) and (2.11)2 to obtainEmbedded Imageand consequentlyEmbedded Image(3.13)We can now calculate F(τ,ρ,θ) following the procedure outlined in appendix A. We then calculate I(τ,ρ,θ) using relation (3.1) if τ≥0, alternatively, the relationEmbedded Image(3.14)if τ<0. For the numerical calculation of the integrals appearing in equations (3.1) and (3.14) we use the Gauss–Legendre quadrature with two functional evaluations at every step, i.e.Embedded Imagewhere the abscissas τ1, τ2 and the weights w1, w2 are given byEmbedded ImageWe also note that we have tried subdivision of the interval (α, β) into several intervals and the improvement is very minor. Therefore, we use just one interval, i.e. two function evaluations per quadrature, since the major increase in running time of the program implicit in using panel quadrature is not justified by the modest improvement in accuracy.

For the numerical calculation of the integrals in equations (3.6) and (3.13) we use the formulaEmbedded ImageFor the numerical calculation of the partial derivatives Embedded Image and Embedded Image in equation (3.6) we use the forward difference schemeEmbedded Imagefor the first half of the interval [−1,1], and the backward difference schemeEmbedded Imagefor the second half.

Thus, for the numerical calculation of g(x1,x2) from the data Embedded Image and Embedded Image we apply the following procedure. First, we calculate the second derivatives Embedded Image, using subroutine spline from Press et al. (1992), setting Embedded Image (i.e. we use the natural cubic spline interpolation). Consequently, we calculate h(ρ,θ) using equations (3.9) and (3.10) for all given ρ and θ. We note that if |ρi|=1 then, since we have assumed compact support, Embedded Image, thus the first term in equation (3.9) is absent. We then calculate fcpe(ρ,θ) and fspe(ρ,θ) using equation (3.2), as well as fc(ρ,θ) and fs(ρ,θ) using equation (3.4) (at this stage we use the second data function Embedded Image). Finally we calculate, again using spline, the second derivatives for the natural cubic spline interpolation of the functions fc(ρ, θ) and fs(ρ,θ).

Having calculated all the necessary second derivatives we now proceed as follows. First, we calculate Embedded Image for any x1, x2 (and θ) using equations (2.11) and (3.8). For this purpose we have used subroutine splint from Press et al. (1992). Consequently, we calculate h(ρ,θ) using equation (3.11). Then, we calculate fcme(ρ,θ) and fsme(ρ,θ) using equation (3.3), fc(ρ,θ) and fs(ρ,θ) using splint and finally, hc(ρ,θ) and hs(ρ,θ) using relations similar to equation (3.11). These last six functions are used in equation (3.7). We then calculate I(τ,ρ,θ) as described earlier. Finally, we calculate r(τ,ρ,θ) using equation (3.7) and, consequently, g(x1,x2) using equation (3.6).

4. Numerical tests

The θ points are equally spaced in [0,2π], while the ρ points are equally spaced in [−1,1]. The density plots presented below were drawn by using Mathematica (Wolfram 1999). The dark colour represents zero (or negative) values while the white colour represents the maximum value of the original (or reconstructed) function.

First, we tested the PET algorithm for the three different phantoms shown in figure 3. Figure 3a,b were taken from Kunyansky (2001) and Guillement & Novikov (2004), respectively. These figures depict the X-ray attenuation coefficient for a function f(x1,x2) modelling a section of a human thorax. The small circles represent bones and the larger ellipses the lungs. Figure 3c is the well-known Shepp–Logan phantom, which provides a model of a head section. All these phantoms consist of different ellipses with various densities.

Using the Radon transform (1.4), we computed the data function Embedded Image for 200 points for θ and 100 points for ρ. This computation was carried out by using Mathematica. We then used these data in the numerical algorithm to reevaluate f(x1,x2). Furthermore, in order to remove the effect of the Gibbs–Wilbraham phenomenon, we applied an averaging filter as follows. We first found the maximum value (max) of f(x1,x2) in the reconstructed image. We then set to zero those values of f(x1,x2), which were less than 1/20 max. Finally, we applied the averaging filter with averaging parameter a=0.005. This filtering procedure was applied five times, with the additional elimination of those values of f(x1,x2), which were less than 1/20 max at the end of the procedure. In figures 4 and 5 we present the results before and after the filtering procedure, respectively. The reconstruction took place in a 500×500 grid.

Figure 4

The reconstruction of the phantoms of figure 3 before the filtering procedure.

Figure 5

The reconstruction of the phantoms of figure 3 after the filtering procedure.

We then tested the SPECT algorithm for the three different phantoms shown in figure 6. Figure 6a,b were taken from Kunyansky (2001). In these cases the function f(x1,x2) is given by figure 3a. Figure 6c was taken from Guillement & Novikov (2004). The white ring represents the distribution of the radiopharmaceutical at the myocardium. In this case the function f(x1,x2) is given by figure 3b.

Figure 6

Test phantoms for the SPECT algorithm. In (a) and (b) the function f(x1,x2) is given by figure 3a, while in (c) the function f(x1,x2) is given by figure 3b.

By using the Radon transform (1.4) and the attenuated Radon transform (1.5), we computed the data functions Embedded Image and Embedded Image for 200 values of θ and 100 points of ρ (again using Mathematica). We consequently used these data in our program to reconstruct g(x1,x2). In order to remove the effect of the Gibbs–Wilbraham phenomenon, a median filter was used, with the additional elimination of those values of g(x1,x2), which were less than 1/20 max before and after the application of the filter. The results are shown in figures 7 and 8, before and after the filtering procedure, respectively. The reconstruction took place in a 140×140 grid.

Figure 7

The reconstruction of the phantoms of figure 6 before the filtering procedure.

Figure 8

The reconstruction of the phantoms of figure 6 after the filtering procedure.

For the above phantoms it seems that even a rough estimation of F(τ,ρ,θ) is sufficient for an accurate reconstruction. This means that, in order to compute numerically F(τ,ρ,θ) using equation (3.13), it is sufficient to use 10 equally spaced points for t, rather than 200. This reduces considerably the reconstruction time.

Acknowledgments

V.M. was supported by a Marie Curie Individual Fellowship of the European Community under contract number HPMF-CT-2002-01597. A.S.F. acknowledges support from EPSRC. We are grateful to Prof. B. Hutton for useful suggestions.

Footnotes

    • Received May 12, 2005.
    • Accepted June 30, 2005.

References

View Abstract