## 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.

## 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 *x*_{1} and *x*_{2} 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*(*x*_{1},*x*_{2}) 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*(*x*_{1},*x*_{2}) 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 by(1.1)where *τ* is a parameter along *L*, and *L*(*x*) denotes the section of *L* between the point (*x*_{1},*x*_{2}) and the detector. The attenuation coefficient *f*(*x*_{1},*x*_{2}) is precisely the function measured by the usual computed tomography. Thus, the basic mathematical problem in SPECT is to determine the function *g*(*x*_{1},*x*_{2}) from the knowledge of the ‘transmission’ function *f*(*x*_{1},*x*_{2}) (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 by(1.2)where *L*_{+}, *L*_{−} are the two half-lines of *L* with endpoint *x*. Since, *L*_{+}+*L*_{−}=*L*, equation (1.2) becomesWe recall that the line integral of the function *f*(*x*_{1},*x*_{2}) along *L* is precisely what is known from the measurements in the usual computed tomography. Thus, since both *I* and the integral of *f*(*x*_{1},*x*_{2}) are known (from the measurements of PET and of computed tomography, respectively), the basic mathematical problem of PET is to determine *g*(*x*_{1},*x*_{2}) from the knowledge of its line integrals. This mathematical problem is identical with the basic mathematical problem of computed tomography.

### 1.1 Notation

A point of a line

*L*making an angle*θ*with the*x*_{1}-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).The above parameterization implies that, for a fixed

*θ*, the Cartesian coordinates (*x*_{1},*x*_{2}) can be expressed in terms of the*local coordinates*(*τ*,*ρ*) by the equations (see §2)(1.3)A function*f*(*x*_{1},*x*_{2}) rewritten in local coordinates will be denoted by*F*(*τ*,*ρ*,*θ*),Thus,*F*(*τ*,*ρ*,*θ*) and*G*(*τ*,*ρ*,*θ*) will denote the X-ray attenuation coefficient*f*(*x*_{1},*x*_{2}) and the distribution of the radiopharmaceutical*g*(*x*_{1},*x*_{2}), rewritten in local coordinates.The line integral of a function

*f*is called its*Radon transform*and will be denoted by . In order to compute , we first write*f*in local coordinates and then integrate with respect to*τ*,(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 . In order to compute , we write both*g*and*f*in local coordinates and then evaluate the following integral:(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 , i.e. to solve equation (1.4) for *f*(*x*_{1},*x*_{2}) in terms of . The relevant formula is called the *inverse Radon transform* and is given by(1.6)where −∞<*x*_{j}<∞, *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 equation(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 , i.e. this approach can be used to solve equation (1.5) for *g*(*x*_{1},*x*_{2}) in terms of and *f*(*x*_{1},*x*_{2}). The relevant formula, called the *inverse attenuated Radon transform*, was obtained by Novikov (2002) by analysing, instead of equation (1.7), the equation(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 function(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 3*c*).

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 by*(2.1)*where x*_{1}*, x*_{2} *are the real Cartesian coordinates* −∞<*x*_{j}<∞*, j*=1*,* 2*, and λ is a complex variable,* *λ*≠0. *Assume that the function f*(*x*_{1},*x*_{2}) *has sufficient decay as* |*x*_{1}|+|*x*_{2}|→∞. *Let μ*(*x*_{1}, *x*_{2}, *λ*) *satisfy the equation*(2.2)*as well as the boundary condition* *μ*=*O*(1/*z*) *as* |*x*_{1}|+|*x*_{2}|→∞. *Let λ*^{+}(*θ*) *and λ*^{−}(*θ*) *denote the limits of λ as it approaches the unit circle (**figure 2**) from inside and outside the unit disc,* *respectively,* *i.e.**Then*(2.3)*where* *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.*(2.4)*and* ∮ *denotes 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 find(2.5)where and denote the complex conjugates of *z* and *λ*, respectively. Equations (2.1) and (2.5) define a change of variables from (*x*_{1},*x*_{2}) to . Using this change of variables to compute and in terms of ∂_{z} and , 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 is(2.6)Indeed, suppose that the function *μ*(*z*_{R}, *z*_{I}) satisfies the equationas well as the boundary condition *μ*=*O*(1/*z*) as *z*→∞. Then, Pompieu's formula (see, for example, Ablowitz & Fokas 1997) implies(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 can be computed explicitly using the formula(2.8)In our casethus 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:(2.9)whereThus, we need to compute the limits of *μ* as *λ* tends to *λ*^{±}. As *ϵ*→0,Substituting this expression in the definition of *z* (equation (2.1)) and simplifying, we find(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. ThenorHence,

*x*

_{1}and

*x*

_{2}are given by equation (1.3). Inverting these equations we find(2.11)Thus, equation (2.10) becomesSubstituting this expression in equation (2.6) and using the fact that the relevant sign equals 1, we find(2.12)Using the change of variables (

*x*

_{1},

*x*

_{2})↔(

*τ*,

*ρ*) defined by equations (1.3) and (2.11), and noting that the relevant Jacobian is 1, i.e.we find that the right-hand side of equation (2.12) equals(2.13)In order to simplify this expression we split the integral over d

*τ*′ in the formand 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) becomesFinally, adding and subtracting the integral we findThe first two terms in the right-hand side of this equation equal , 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 form(2.14)where(2.15)

### 2.1 The inverse radon transform

Equation (2.3) yields(2.16)Equation (2.9) impliesSubstituting this expression in equation (1.7) we find(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 formwhere *ν* is defined by equation (2.15). Hence,orReplacing in this equation by the right-hand side of equation (2.14) we findFor 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,(2.18)Note that the term is independent of *τ*′; thus this term comes out of the integral , and furthermore the same term appears in the left-hand side of equation (2.18). Hence, when computing the jump *μ*(*x*_{1}, *x*_{2}, *λ*^{+})−*μ*(*x*_{1}, *x*_{2}, *λ*^{−}), the second term in the right-hand side of equation (2.18) cancels and we find that the relevant jump in now given by(2.19)where *τ* and *ρ* are expressed in terms of *x*_{1} and *x*_{2} 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 find(2.20)where *J* is defined by equation (2.19). This formula is equivalent to Novikov's formula.

In summary, let be defined by equation (1.5), let *F*(*τ*,*ρ*,*θ*) denote the function *f*(*x*_{1},*x*_{2}) written in local coordinates (see §1.1) and let denote the Radon transform of *f*(*x*_{1},*x*_{2}) (see equation (1.4)). Then *g*(*x*_{1},*x*_{2}) is given by equation (2.20) where the function *J* is explicitly given in terms of and 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.(3.1)Note that, we assume compact support; thus the integration domain is finite, i.e. and *F*(*τ*,*ρ*,*θ*)=0 for |*ρ*|≥1, or for .

Definition (2.4) becomesMoreover,We introduce the following notation:(3.2)(3.3)(3.4)Using this notation and setting *R*(*τ*,*ρ*,*θ*)=−*J*(*τ*,*ρ*,*θ*), after some calculations, equation (2.19) becomes(3.5)We now setthus equation (3.5) becomesWe denote the right-hand side of this equation by −i*r*(*τ*,*ρ*,*θ*). Taking the real part of *g*(*x*_{1},*x*_{2}) in equation (2.20), we obtain(3.6)where *τ* and *ρ* are given by equation (2.11) and(3.7)For the numerical calculation of the Hilbert transform we suppose that is given, for every *θ*, at *n* equally spaced points *ρ*_{i}∈[−1,1], i.e. we suppose that are known. Moreover, in each interval [*ρ*_{i},*ρ*_{i+1}] we approximate using the relation(3.8)whereand denotes the second derivative of with respect to *ρ*, at *ρ*=*ρ*_{i}. In other words, we approximate by a cubic spline (in *ρ*) with equally-spaced nodes.

We then write(3.9)If *ρ*=*ρ*_{i} or *ρ*=*ρ*_{i+1} the integral in the right-hand side of equation (3.9) can be written asThus, after some calculations, we obtain(3.10)If *ρ*≠*ρ*_{i} and *ρ*≠*ρ*_{i+1} the integral in the right-hand side of equation (3.9) can be written asand after some calculation we obtain(3.11)where *F*_{i} is the right-hand side of equation (3.10).

Taking the real part of equation (1.6) it follows that *f*(*x*_{1},*x*_{2}) is given by(3.12)where *h*(*ρ*,*θ*) is defined by equation (1.9). We assume that *f*(*x*_{1},*x*_{2}) has compact support, namely *f*(*x*_{1},*x*_{2})=0, for . Thus, in order to calculate numerically *I*(*τ*,*ρ*,*θ*) for any *x*_{1}, *x*_{2}, *θ*, we use relation (3.12) and (2.11)_{2} to obtainand consequently(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 relation(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.where the abscissas *τ*_{1}, *τ*_{2} and the weights *w*_{1}, *w*_{2} are given byWe 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 formulaFor the numerical calculation of the partial derivatives and in equation (3.6) we use the forward difference schemefor the first half of the interval [−1,1], and the backward difference schemefor the second half.

Thus, for the numerical calculation of *g*(*x*_{1},*x*_{2}) from the data and we apply the following procedure. First, we calculate the second derivatives , using subroutine spline from Press *et al*. (1992), setting (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, , thus the first term in equation (3.9) is absent. We then calculate *f*^{cpe}(*ρ*,*θ*) and *f*^{spe}(*ρ*,*θ*) using equation (3.2), as well as *f*^{c}(*ρ*,*θ*) and *f*^{s}(*ρ*,*θ*) using equation (3.4) (at this stage we use the second data function ). Finally we calculate, again using spline, the second derivatives for the natural cubic spline interpolation of the functions *f*^{c}(*ρ*, *θ*) and *f*^{s}(*ρ*,*θ*).

Having calculated all the necessary second derivatives we now proceed as follows. First, we calculate for any *x*_{1}, *x*_{2} (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 *f*^{cme}(*ρ*,*θ*) and *f*^{sme}(*ρ*,*θ*) using equation (3.3), *f*^{c}(*ρ*,*θ*) and *f*^{s}(*ρ*,*θ*) using splint and finally, *h*^{c}(*ρ*,*θ*) and *h*^{s}(*ρ*,*θ*) 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*(*x*_{1},*x*_{2}) 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 3*a*,*b* were taken from Kunyansky (2001) and Guillement & Novikov (2004), respectively. These figures depict the X-ray attenuation coefficient for a function *f*(*x*_{1},*x*_{2}) modelling a section of a human thorax. The small circles represent bones and the larger ellipses the lungs. Figure 3*c* 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 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*(*x*_{1},*x*_{2}). 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*(*x*_{1},*x*_{2}) in the reconstructed image. We then set to zero those values of *f*(*x*_{1},*x*_{2}), 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*(*x*_{1},*x*_{2}), 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.

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

By using the Radon transform (1.4) and the attenuated Radon transform (1.5), we computed the data functions and for 200 values of *θ* and 100 points of *ρ* (again using Mathematica). We consequently used these data in our program to reconstruct *g*(*x*_{1},*x*_{2}). 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*(*x*_{1},*x*_{2}), 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.

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.

- © 2005 The Royal Society