Algorithms in Tomography

Frank Natterer

Institut für Numerische und instrumentelle Mathematik Universität Münster Einsteinstrae 62, D-48149 Münster (Germany)

1 Introduction

The basic problem in computerized tomography is the reconstruction of a function from its line or plane integrals. Applications come from diagnostic radiology, astronomy, electron microscopy, seismology, radar, plasma physics, nuclear medicine and many other fields. More recent kinds of tomography replace the straight line model by an inverse problem for a partial differential equation.

The outline of this paper is as follows. In section 2 we survey the mathematical models used in tomography. In section 3 we give a fairly detailed survey on 2D reconstruction algorithm which still are the work horse of tomography. In section 4 we describe recent developments in 3D reconstruction. In section 5 we make a few remarks on the beginning development of algorithms for nonstraight-line tomography.

2 Mathematical Models in Tomography

In the description of the mathematical models we restrict ourselves to those features which are important for the mathematical scientist. For physical and medical aspects see [35], [6], [32].

(a) Transmission CT.

This is the original and simplest case of CT. In transmission tomography one probes an object with non diffracting radiation, e.g. x-rays for the human body. If I_{0} is the intensity of the source, a(x) the linear attenuation coefficient of the object at point x, L the ray along which the radiation propagates, and I the intensity past the object, then

I = I_{0} e^{− ∫_{L} a(x)dx} .  (2.1)

In the simplest case the ray L may be thought of as a straight line. Modeling L as strip or cone, possibility with a weight factor to account for detector inhomogeneites may be more appropriate. In (2.1) we neglegt the dependence of a from the energy (beam hardening effect) and other non linear phenomena (e.g. partial volume effect).

The mathematical problem in transmission tomography is to determine a from measurements of I for a large set of rays L. If L is simply the straight line connecting the source x_{0} with the detector x_{1}, (2.1) gives rise to the integrals

ℓn\frac{I}{I_{0}} = − ∫_{x_{0}}^{x_{1}}a(x)dx  (2.2)

where dx is the restriction to L of the Lebesgue measure in R^{n}. We have to compute a in a domain Ω ⊆ R^{n} from the values of (2.2) where x_{0}, x_{1} run through certain subsets of ∂Ω.

For n = 2, (2.2) is simply a reparametrization of the Radon transform

(Ra)(θ, s) = ∫_{x·θ=s} a(x)dx  (2.3)

where θ ∈ S^{n−1}. Thus our problem is in principle solved by Radon’s inversion formula

a = R^{∗}Kg ,  g = Ra  (2.4)

where

(R^{∗}g)(x) = ∫_{S^{n−1}} g(θ, x · θ)dθ  (2.5)

is the so-called backprojection and

K = \frac{1}{2}(2π)^{1−n} \{ \begin{matrix}(−1)^{(n−2)/2} H \frac{∂^{n−1}}{∂s^{n−1}} & , & n even \\ (−1)^{(n−1)/2} \frac{∂^{n−1}}{∂s^{n−1}} & , & n odd \end{matrix}  (2.6)

with H the Hilbert transform [39]. In fact the numerical implementation of (2.4) leads to the filtered backprojection algorithm which is the standard algorithm in commercial CT scanners, see section 3.

For n = 3, the relevant integral transform is the x-ray transform

(P a)(θ, x) = ∫_{R^{1}} a(x + sθ)ds

where θ ∈ S^{n−1} and x ∈ θ^{⊥}. P admits a similar inversion formula as R, to wit

a = P^{∗}Kg  g = P f  (2.7)

with K very similar to (2.6) and

(P^{∗}g)(x) = ∫_{S^{n−1}} g(θ, E_{θ}x)dθ

where E_{θ} is the orthogonal projection onto θ^{⊥}. Unfortunately, (2.7) is not as useful as (2.4). The reason is that (2.7) requires g for all θ and y ∈ θ^{⊥}, i.e. (2.2) has to be available for all x_{0}, x_{1} ∈ ∂Ω. This is not practical. Also, it is not necessary for unique reconstruction of a. In fact it can be shown that a can be recovered uniquely from (2.2) with sources x_{0} on a circle surrounding supp(a) and x_{1} ∈ S^{n−1}. Unfortunately the determination of a in such an arrangement is, though uniquely possible, highly unstable. The condition of stability is the following: Each plane meeting supp(a) must contain at least one source [16]. This condition is obviously violated for sources on a circle. Cases in which the condition is satisfied include the helix and a pair of orthogonal circles. A variety of inversion formulae has been derived, see section 4.

If scatter is to be included, a transport model is more appropriate. Let u(x, θ) be the density of the particles at x travelling (with speed 1) in direction θ. Then,

θ · ∇u(x, θ) + a(x)u(x, θ) = ∫_{S^{n−1}} η(x, θ, θ′)u(x, θ′)dθ′
+δ(x − x_{0}) .  (2.8a)

Here, η(x, θ, θ′) is the probability that a particle at x travelling in direction θ is scattered in direction θ′. Again we neglect dependence on energy. δ is the Dirac δ-function modeling a source of unit strength. (2.2a) holds in a domain Ω of R^{n} (n = 2 or 3), and x_{0} ∈ ∂Ω. Since no radiation comes in from outside we have

u(x, θ) = 0 , x ∈ ∂Ω , ν_{x} · θ ≤ 0  (2.8b)

where ν_{x} is the exterior normal on ∂Ω at x ∈ ∂Ω. (4.1) is now replaced by

I(x_{1}, x_{0}, θ) = I_{0}u(x_{1}, θ) , x_{1} ∈ ∂Ω , ν_{x_{1}} · θ ≥ 0 .  (2.8c)

The problem of recovering a from (2.8) in much harder. An explicit formula for a such as (2.4) has not become known and is unlikely to exist. Nevertheless numerical methods have been developed for special choices of η [3], [7]. The situation gets even more difficult if one takes into account that, strictly speaking, η is object dependent and hence not known in advance. (2.8) is a typical example of an inverse problem for a partial differential equation. In an inverse problem one has to determine the differential equation - in our case a, η - from information about the solution - in our case (2.8c).

(b) Emission CT.

In emission tomography one determines the distribution f of radiating sources in the interior of an object by measuring the radiation outside the object in a tomographic fashion. Let again u(x, θ) be the density of particles at x travelling in direction θ with speed 1, and let a be the attenuation distribution of the object. (This is the quantity which is sought for in transmission CT). Then,

θ · ∇u(x, θ) + a(x)u(x, θ) = f(x) .  (2.9a)

This equation holds in the object region Ω ⊆ R^{n} for each θ ∈ S^{n−1}. Again there exists no incoming radiation, i.e.

u(x, θ) = 0 , x ∈ ∂Ω , ν_{x} · θ ≤ 0  (2.9b)

while the outgoing radiation

u(x, θ) = g(x, θ) , x ∈ ∂Ω , ν_{x} · θ ≥ 0  (2.9c)

is measured and hence known. (2.9) again constitutes an inverse problem for a transport equation. For a known, (2.9a-b) is readily solved to yield

u(x, θ) = ∫_{x−∞·θ}^{x} e^{− ∫_{y}^{x} ads}f(y)dy .

Thus (2.9c) leads to the integral equation

g(x, θ) = ∫_{x−∞·θ}^{x} e^{− ∫_{y}^{x} ads}f(y)dy  (2.10)

for f. Apart from the exponential factor, (2.10) is identical - up to notation - to the integral equation in transmission CT. Except for very special cases - e.g. a constant in a known domain [52], [42] no explicit inversion formulas are available. Numerical techniques have been developed but are considered to be slow. Again the situation becomes worse if scatter is taken into account. This can be done by simply adding the scattering integral in (2.8a) to the right hand side of (2.9a).

What we have described to far is called SPECT (= single particle emission CT). In PET (= positron emission tomography) the sources eject the particles pairwise in opposite directions. They are detected in coincidence mode, i.e. only events with two particles arriving at opposite detectors at the same time are counted. (1.10) has to be replaced by

g(x, θ) = ∫_{x−∞·θ}^{x} e^{− ∫_{y}^{x} ads− ∫_{x−∞·θ}^{y} ads}f(y)dy

= e^{− ∫_{x−∞·θ}^{x} ads} ∫_{x−∞·θ}^{x} f(y)dy .  (2.11)

Thus PET is even closer to the case of transmission CT. If a is known, we simply have to invert the x-ray transform. Inversion formulas which can make use of the data collected in PET are available, see section 3.

The main problems in emission CT are unknown attenuation, noise, and scatter. For the attenuation problem, the ideal mathematical solution would be a method for determining f and a in (2.9) simultaneously. Under strong assumption on a (e.g. a constant in a known region [26], a affine distortion of a prototype [37], a close to a known distribution [8]) encouraging results have been obtained. Theoretical results based on the transport formulation have been obtained, even for models including scatter [43]. But a clinically useful way of determining a from the emission data has not yet been found.

Noise and scatter are stochastic phenomena. Thus, besides models using integral equations, stochastic models have been set up for emission tomography [48]. These models are completely discrete. We subdivide the reconstruction region into m pixels or voxels. The number of events in pixel/voxel j is a Poisson random variable φ_{j} whose mathematical expection f_{j} = Eφ_{j} is a measure for the activity in pixel/voxel j. The vector f = (f_{1}, …, f_{n}) is the sought - for quantity. The vector g = (g_{1}, …, g_{n}) of measurements is considered as realization of the random variable γ = (γ_{1}, …, γ_{n}) where γ_{i} is the number of events detected in detector i. The model is determined by the (n, m)-matrix A = (a_{ij}) whose elements are

a_{ij} = P (event in pixel/voxel j detected in detector i)

where P denotes probability. We have E(γ) = Af. f is determined from g by the maximum liklihood method. A numerical method for doing this is the EM (= expection maximation) algorithm. In its basic form it reads

f^{k+1} = f^{k}A^{∗} \frac{g}{Af^{k}} , k = 0, 1, …

where division and multiplication are to be understood component wise. The problem with the EM-algorithm is that it is only semi-convergent, i.e. noise is amplified at high iteration numbers. This is known as the checkerboard effect. Various suggestions have been made to get rid of this effect. The most exciting and interesting ones use “prior information” and attempt to maximize “posterior likelihood”. Thus f is assumed to have a prior probability distribution, called a Gibbs-Markov random field π(f), which gives preference to certain functions f [46], [21], [18]. Most prior π simply add a penalty term to the likelihood function to account for correction between neighboring pixels and do not use biological information. However if π is carefully chosen so that piecewise constant functions f with smooth boundaries forming the region of constancy are preferred then the noise amplification at high iteration numbers can be avoided. The question remains as to whether this conclusion will remain valid for function f which are assigned low probability by π - or, more to the point - whether “real” emission densities f will be well-resolved by this Bayesian method. An ROC study (e.g. double-blind trials where radiologists are to find lesions from images produced by two different algorithms) concluded that maximum likelihood methods were superior to the filtered backprojection algorithm (see section 3) in certain clinical applications. The same type of study is needed to determine whether or not Gibbs priors will improve the maximum likelihoof reconstruction (stopped short of convergence to avoid noise amplification) on real data.

(c) Ultrasound CT

X-rays travel along straight lines. For other sources of radiation, such as ultrasound and microwaves, this is no longer the case. The paths are no longer straight, and their exact shape depends on the internal structure of the object. We can no longer think in terms of simple projections and linear integral equations. More sophisticated non linear models have to be used.

In the following we consider an object Ω ⊆ R^{n} with refractive index n. We assume n = 1 outside the object. The object is probed by a plane wave

e^{−ikt}u_{θ}(x) , u_{θ}(x) = e^{ikx·θ}

with wave number k = \frac{2π}{λ}, λ the wave length, travelling in the direction θ. The resulting wave e^{−ikt}u(x) satisfies the reduced wave equation

Δu + k^{2}(1 + f)u = 0 , f = n^{2} − 1 ,  (2.12a)

plus suitable boundary conditions at infinity. The inverse problem to be solved is now the following. Assume that

g(x, θ) = u(x) , θ ∈ S^{n−1}  (2.12b)

is known outside Ω. Determine f inside Ω!

Uniqueness and stability of the inverse problem (2.12) has recently been settled [36]. However, stability is only logarithmic [2], i.e. a data error of size δ results in a reconstruction error 1/ log(1/δ). Numerical algorithms did not emerge from this work.

Numerical methods for (2.12) are mostly based on linearizations, such as the Born and Rytov approximation [13]. In order to derive the Born approximation, one rewrites (2.12a) as

u(x) = u_{θ}(x) − k^{2} ∫_{Ω} G(x − y)f(y)u(y)dy  (2.13)

where G is an appropriate Green’s function. For n = 3, we have

G(x) = \frac{e^{ik|x|}}{4π|x|} .  (2.14)

The Born approximation is now obtained by assuming u ∼ u_{θ} in the integral in (2.13). With this approximation, (2.12b) reads

g(x, θ) = u_{θ}(x) − k^{2} ∫_{Ω} G(x − y)u_{θ}(y)f(y)dy , x ∉ Ω .  (2.15)

This is a linear integral equation for f, valid for all x outside the object and for all measured directions θ.

Numerical methods based on (2.15) - and a similar equation for the Rytov approximation - have become known as diffraction tomography. Unfortunately, the assumptions underlying the Born- and Rytov approximations are not satisfied in medical imaging. Thus, the reconstructions of f obtained from (2.15) are very poor. However, we may use (2.15) to get some encouraging information about stability. For |x| large, (2.15) assumes the form

g(x, θ) = u_{θ}(x) − \frac{k^{2}}{4π|x|}e^{ik|x|} ∫ e^{ik(θ− \frac{x}{|x|})·y}f(y)dy

= u_{θ}(x) − \frac{k^{2}}{4π|x|}e^{ik|x|}(2π)^{3/2}f̂(k(\frac{x}{|x|} − θ))  (2.16)

with f̂the Fourier transform of f. (2.16) determines f̂within a ball of radius \sqrt{2}k from the data (2.12b) in a completely stable way. We conclude that the stability of the inverse problem (2.12) is much better than logarithmic. If the resolution is restricted to spatial frequencies below \sqrt{2}k - which is perfectly reasonable from a physical point of view - then we can expect (2.12) to be perfectly stable.

So far we considered plane wave irradiation at fixed frequency, and we worked in the frequency domain. Time domain methods are conceivable as well. We start out from the wave equation

\frac{∂^{2}u}{∂t^{2}} = c^{2}Δu  (2.17a)

with the propagation speed c assumed to be 1 outside the object. With x_{0} a source outside the object we consider the initial conditions

u(x, 0) = 0 , \frac{∂u}{∂t}(x, 0) = δ(x − x_{0}) .  (2.17b)

We want to determine c inside the object from knowing

g(x_{0}, x_{1}, t) = u(x_{1}, t) , t > 0  (2.17c)

for many sources x_{0} and receivers x_{1} outside the object. In the one dimensional case, the inverse problem (2.17) can be solved by the famous Gelfand-Levitan method in a stable way. It is not clear how Gelfand-Levitan can be extended to dimensions two and three. The standard methods use sources and receivers on all of the boundary of the object. This is not practical in medical imaging. However, for reduced data sets, comparable to those in 3D X-ray tomography, we do not know how to use Gelfand-Levitan, nor do we know anything about stability.

Of course one can always solve the nonlinear problem (2.17) by a Newton type method. Such methods have been developed [23], [30]. They suffer from excessive computing time and from their apparent inability to handle large wave numbers k.

(d) Optical tomography

Here one uses NIR (= near infra-red) lasers for the illumination of the body. The process is now described by the transport equation

\frac{∂u}{∂t}(x, θ, t) + θ · ∇u(x, θ, t) + a(x)u(x, θ, t)

= b(x) ∫_{S^{n−1}} η(θ · θ′)u(x, θ′, t)dθ′ + f(x, θ, t)  (2.18a)

for the density u(x, θ, t) of the particles at x ∈ Ω flying in direction θ ∈ S^{n−1} at time t. a and b are the sought - for tissue parameters. The scattering kernel η is assumed to be known. The source term f is under the control of the experimenter. Together with the initial and boundary conditions

u(x, θ, 0) = 0 in Ω × S^{n−1} ,

(2.18b)

u(x, θ, t) = 0 on ∂Ω × S^{n−1} × R^{1} , v_{x} · θ ≤ 0

(1.18a) has a unique solution under natural conditions on a, b, η and f. As in (1.1) we pose the inverse problem. Assume that we know the outward radiation

g(x, θ, t) = u(x, θ, t) on ∂Ω × S^{n−1} × R^{1} , ν_{x} · θ ≥ 0 ,  (2.18c)

can we determine one or both the quantities a, b?

There are essentially three methods for illuminating the object, i.e. for choosing the source term f in (2.18a). In the stationary case one puts f = δ(x − x_{0}) where x_{0} ∈ ∂Ω is a source point. u is considered stationary, too. A second possibility is the light flash f = δ(x − x_{0})δ(t). Finally one can also use time harmonic illumination, in which case f = δ(x − x_{0}) e^{iωt}. This case reduces to the stationary case with a replaced by a+iω. In all three cases, the data function g of (2.18c) is measured at x ∈ ∂Ω, possibly averaged over one or both of the variables θ, t.

Light tomography is essentially a scattering phenomenon. This means that the scattering integral in (2.18a) is essential. It can no longer be treated merely as a perturbation as in x-ray CT. Thus the mathematical analysis and the numerical methods are expected to be quite different from what we have seen in other types of tomography.

The mathematical theory of the inverse problem (2.18) is in a deplorable state. There exist some Russian papers on uniqueness [4]. General methods have been developed, too, but apparently they have been applied to 1D problems only [44]. Nothing seems to be known about stability. The numerical methods which have become known are of the Newton type, either applied directly to the transport equation or to the so-called diffusion approximation [5], [31]. The diffusion approximation is an approximation to the transport equation by a parabolic differential equation. Since inverse problems for parabolic equations are severely ill-posed, this approach is questionable. Higher order approximations [29], [20] are hyperbolic, making the inverse problem much more stable.

As an alternative to the transport equation one can also model light tomography by a discrete stochastical model [22]. In the 2D case, break up the object into a rectangular arrangement of pixels labelled by indices i, j with a ≤ i ≤ b and c ≤ j ≤ d. Attach to each pixel the quantities f_{ij}, b_{ij}, r_{ij}, ℓ_{ij} meant to denote the probability of a forward, backward, rightward or leftward transition out of the pixel i, j with respect to the direction used to get into this pixel. For each pair of boundary pixels i, j and i′, j′ let P_{ij,i′} be the probability that a particle that enters the object at pixel i, j will eventually leave the object at pixel i′, j′. The problem is to determine the quantities f_{ij}, b_{ij}, r_{ij}, ℓ_{ij} from the values of P_{ij,i′j′} for all boundary pixels. Preliminary numerical tests show that this is possible, at least in principle. However, the computations are very time consuming. More seriously, they reveal a very high degree of instability.

(e) Electrical Impedance Tomography

Here, the sought-for quantity is the electrical impedance σ of an object Ω. Voltages are applied via electrodes on ∂Ω, and the resulting currents at these electrodes are measured. With u the potential in Ω, we have

div(σ · ∇u) = 0 in Ω

(2.19)

u = g , σ \frac{∂u}{∂ν} = f on ∂Ω .

Knowing many voltage - current pairs g, f on ∂Ω, we have to determine σ from (2.19).

Uniqueness for the inverse problem (2.19) has recently be settled [36]. Unfortunately, the stability properties are very bad [2]. Numerical methods based on Newton’s method, linearization, simple backprojection, layer stripping have been tried. All these methods suffer from the severe ill-posedness of the problem. There seems to be no way to improve stability by purely mathematical means.

(f) Magnetic Resonance Imaging (MRI)

The physical phenomena exploited here is the precession of the spin of a proton in a magnetic field of strength H about the direction of that field. The frequency of this precession is the Larmor frequency γH where γ is the gyromagnetic ratio. By making the magnetic field H space dependent in a controled way the local magnetization M_{0}(x) (together with the relaxation times T_{1}(x), T_{2}(x)) can be imaged. In the following we derive the imaging equations [27].

The magnetization M(x, t) caused by a magnetic field H(x, t) satisfies the Bloch equation

\frac{∂M}{∂t} = γM × H − \frac{1}{T_{2}} (M_{1}e_{1} + M_{2}e_{2}) − \frac{1}{T_{1}} (M_{3} − M_{0})e_{3} .  (2.20)

Here, M_{i} is the i-th component of M, and e_{i} is the i-th unit vector i = 1, 2, 3. The significance of T_{1}, T_{2}, M_{0} become apparent if we solve (2.20) with the static field H = H_{0}e_{3} and with initial values M(x, 0) = M^{0}(x). We obtain with ω_{0} = γH_{0}

M_{1}(x, t) = e^{−t/T_{2}}  (M_{1}^{0} cos(ω_{0}t) + M_{2}^{0} sin(ω_{0}t))
M_{2}(x, t) = e^{−t/T_{2}}  (−M_{1}^{0} sin(ω_{0}t) + M_{2}^{0} cos(ω_{0}t))  (2.21)
M_{3}(x, t) = e^{−t/T_{1}}  M_{3}^{0} + (1 − e^{−t/T_{1}})M_{0}

Thus the magnetization rotates in the x_{1} − x_{2}-plane with Larmor frequency ω_{0} and returns to the equilibrium condition (0, 0, M_{0}) with speed controlled by T_{2} in the x_{1} − x_{2}-plane and by T_{1} in the x_{3}-direction.

In an MRI scanner one generates a field

H(x, t) = (H_{0} + G(t) · x)e_{3} + H_{1}(t)(cos(ω_{0}t)e_{1} + sin(ω_{0}t)e_{2})

where G and H_{1} are under control. In the jargon of MRI, H_{0}e_{3} is the static field, G the gradient, and H_{1} the radio frequency (RF) field. The input G, H_{1} produces in the detecting system of the scanner the output signal

S(t) = −\frac{d}{dt} ∫_{R^{3}} M(x, t)B(x)dx  (2.22)

where B characterizes the detecting system. Depending on the choice of H_{1} various approximation to S can be derived.

(i) H_{1} is constant in the small interval [0, τ] and γ ∫_{0}^{τ} H_{1}dt = \frac{π}{2} (Short \frac{π}{2} pulse).

In that case,

S(t) = ∫_{R^{3}} M_{0}(x)e^{−iγ ∫_{0}^{t} G(t′)dt′·x−t/T_{2}(x)}dx .

Choosing G constant for τ ≤ t ≤ τ + T and zero otherwise we get for T ≪ T_{2}

S(t) = ∫_{R^{3}} M_{0}(x)e^{−iγ(t−τ)G·x}dx

(2.23)

= (2π)^{3/2}M̂_{0}(γ(t − τ)G)

where M̂_{0} is the 3D Fourier transform of M_{0}.

From here we can proceed in several ways. We can either use (2.23) to determine the 3D Fourier transform M̂_{0} of M_{0} and to compute M_{0} by an inverse 3D Fourier transform. This requires M̂_{0} to be known on a Cartesian grid, what can be achieved by a proper choice of the gradients or by interpolation. We can also evoke the central slice theorem to obtain the 3D Radon transform RM_{0} of M_{0} by a series of 1D Fourier transforms. M_{0} is recovered in turn by inverting the 3D Radon transform.

The main numerical problem in MRI is to reconstruct a function from inperfect measurements of its Fourier transform [34]. Not much been achied so far.

(ii) H_{1} is the shaped pulse

H_{1}(t) = ϕ(tγG)e^{iγGx_{3}t}

where ϕ is a smooth positive function supported in [0, τ]. Then, with x′, G′ the first to components of x, G, respectively, we have

S(t) = ∫_{R^{2}} M′_{0}(x′, x_{3})e^{−iγ ∫_{0}^{t} G′(t′)dt′·x′−t/T_{2}(x′,x_{3})}dx′  (2.24)

where

M′_{0}(x′, x_{3}) = ∫ M_{0}(x′, y_{3})Q(x′_{3} − y_{3})dy_{3}

with a function Q essentially supported in a small neighborhood of 0. (2.24) is the 2D analogue of (2.23). So we have again the choice between Fourier imaging (i.e. computing the 2D Fourier transform from (2.24) and doing an inverse 2D Fourier transform) and projection imaging (i.e. doing a series of 1D Fourier transforms in (2.24) and inverting the 2D Radon transform).

Some of the mathematical problems of MRI, e.g. interpolation in Fourier space, are common to other techniques in medical imaging. An interesting mathematical problem occurs if the magnets to not produce suffiently homogeneous fields [33]. It calls for the reconstruction of a function from integrals over slightly curved manifolds. Even though this is a problem of classical integral geometry it has not yet found a satisfactory solution.

(g) Vector Tomography

If the domain under consideration contains a moving fluid, then the Doppler shift can be used to measure the velocity u(x) of motion. Assume the time harmonic signal e^{iω_{0}t} is transmitted along the oriented line L. This signal is reflected by particles travelling with speed ν in the direction of L as e^{i(ω_{0}−kν)t}, where k = 2ω_{0}/c, c the speed of the probing signal, i.e. kν is the Doppler shift. Let S(L, ν) be the Lebesgue measure of these particles on L which move with speed < ν, i.e. u(x) · e_{L} < ν, e_{L} the tangent vector on L. Then the total response is

g(L, t) = ∫_{−∞}^{+∞}e^{i(ω_{0}−kν))t}S(L, ν)dν .

Thus S can be recovered from g by a Fourier transform. The problem is to recover u from S.

Not much is known about uniqueness. However, the first moment of S,

∫_{−∞}^{+∞}νS(L, ν)dν = ∫_{L} u(x) · e_{L}dx = (Ru)(L)  (2.25)

is similar to the Radon transform. One can show that curl u can be computed from Ru, and an inversion formula similar to the Radon inversion formula exists. Numerical simulations are given in [51].

(h) Tensor Tomography

As an immediate extension of transmission CT to non-isotropic media we consider a matrix valued attenuation a(x) = (a_{ij}(x)), i, j = 1, …, n. We solve the vector differential equation

\frac{du(t)}{dt} = −a(x(t))u(t) , x(t) = (1 − t)x_{0} + tx_{1} , 0 ≤ t ≤ 1  (2.26)

for the vector valued function u(t) = (u_{i}(t))_{i=1,...,n} . Let a be defined in a convex domain Ω, and let x_{0}, x_{1} ∈ ∂Ω. Then, u(x_{1}) = U(x_{0}, x_{1})u(x_{0}) with a nonlinear map U(x_{0}, x_{1}) depending on a. The problem is to recover a in Ω from the knowledge of U(x_{0}, x_{1}) for x_{0}, x_{1} ∈ ∂Ω. For n = 1 we regain (2.1). Applications of (2.25) for n > 1 have become known in photoelasticity [1], but applications to medicine are not totally out of question.

In a further extension we let a depend on the direction ξ = (x_{1} − x_{0})/|x_{1} − x_{0}|. Such problems occur in the polarization of harmonic electromagnetic and elastic waves in anisotropic media. In linearized form these problems give rise to the transverse x-ray transform

(Ja)(x, θ, ω) = ∫ ω^{T} a(x + tθ)ωdt , ω ⊥ θ  (2.27)

and to the longitudenal x-ray transform

(Ja)(x, θ) = ∫ θ^{T} a(x + tθ)θdt .  (2.28)

(2.27) can easily be reduced to the (n − 1)-dimensional x-ray transform in the plane H_{ω,s} = {y : y · ω = s}. We only have to introduce the function a_{ω}(y) = ω^{T} a(y)ω. Then, (2.27) provides for x ∈ H_{ω,s} all the line integrals of a_{ω} on H_{ω,s}. For (2.28), the situation is not so easy. We decompose a in its solenoidal and potential part, i.e.

a = a_{1} + ∇a_{2} div a_{1} = 0 a_{2} = 0 on ∂Ω .

It can be shown that a_{1} can be recovered uniquely from (2.28), but a_{2} is completely undetermined [47]. This is reminiscent of vector tomography.

3 Basic Algorithms in 2D Tomography

The numerical implementation of Radon’s inversion formula (2.6) is now well understood [39], [25], [28]. We consider only the simplest case of parallel scanning. This means that g(θ, s) = (Ra)(θ, s) is sampled at

θ = θ_{j} = ( \begin{matrix}cos φ_{j} \\ sin φ_{j} \end{matrix} ) ,  φ_{j} = π_{j}/p , j = 0, …, p − 1 ,

(3.1)

s = s_{ℓ} = ℓ\frac{ρ}{q} ,  ℓ = −q, …, q .

The reconstruction region is |x| < ρ.

We need some tools from sampling theory. For a function f in IR^{n} we define the Fourier transform by

f̂(ξ) = (2π)^{−n/2} ∫ e^{−ix·ξ}f(x)dx .

We call a function (essentially) Ω-band-limited if f̂(ξ) vanishes (approximately) for |ξ| > Ω. Shannon’s sampling theorem states that an Ω-band-limited function can be recovered exactly from the samples f(hk), k ∈ Z^{n}, if |h| ≤ π/Ω. If f_{1}, f_{2} both are Ω-band-limited, then

∫_{IR^{n}} f_{1}f_{2}dx = h^{n} ∑_{k∈Z^{n}} f_{1}(hk)f_{2}(hk) .  (3.2)

The convolution of two functions is denoted by f_{1} ⋆ f_{2}, independently of the dimension.

Rather than using (2.6) we start out from

V ⋆ a = R^{∗}(v ⋆ g)  (3.3)

if V = R^{∗}v. Here, R^{∗} is the 2D backprojection from (2.5), and the convolution on the right hand side is with respect to the second argument only. In order to make (3.3) equivalent to (2.6) one has choose v such that V = δ, the Dirac δ-function. One can show that this is the case iff

v̂(σ) = \frac{1}{2}(2π)^{−3/2}|σ| .  (3.4)

Thus v is a distribution rather than a function. For the numerical evaluation of (3.3) we have to regularize v, i.e. we replace it by

v̂_{Ω}(σ) = \frac{1}{2}(2π)^{−3/2}|σ|ϕ̂(σ/Ω)

where ϕ is a low-pass filter, i.e. π(σ) = 0 for |σ| ≥ 1. For the ideal low-pass ϕ(σ) = 1, |σ| ≤ 1, we have

v_{Ω}(s) = \frac{Ω^{2}}{4π^{2}} u(Ωs) , u(s) = sin(s) − \frac{1}{2}(sin(s/2))^{2}

with the sin c-function sin c(x) = sin(x)/x. This filter produces reconstructions with maximal spatial resolution but with unpleasant artefacts. Hence one prefers filter with a smooth transition from non zero to zero values, such as the much used Shepp and Logan filter

ϕ̂(σ) = \{ \begin{matrix}sin c|σ|π/2) & , & 0 ≤ |σ| ≤ 1 , \\ 0 & , & |σ| ≥ 1 . \end{matrix}

The convolution on the right hand side of (3.3) is now approximated by the trapezoidal rule

(v_{Ω} ⋆ g)(θ, s) = h ∑_{ℓ=−q}^{q} v_{Ω}(s − s_{ℓ})g(θ, s_{ℓ}) .  (3.5)

At this stage we assume a to be essentially Ω-band-limited. Since

(Ra)^{∧}(θ, σ) = (2π)^{1/2}â(σθ)  (3.6)

where (Ra)^{∧} means the 1D Fourier transform with respect to the second argument, Ra is essentially Ω-band-limited, too, and (3.5) is correct with good accuracy for \frac{ρ}{q} ≤ \frac{π}{Ω}, as can be seen from (3.2). For the backprojection we use the trapeziodal rule again, obtaining as approximation to a(x)

a(x) = \frac{πh}{p} ∑_{j=0}^{p−1} ∑_{ℓ=−q}^{q} v_{Ω}(x · θ_{j} − s_{ℓ})g(θ_{j}, s_{ℓ}) .  (3.7)

It can be shown that (v_{Ω} ⋆ g)(θ, θ · x) as a function of φ, θ = (cos φ, sin φ)^{T}, is essentially Ωρ-band-limited. Thus, (3.7) is satisfied with good accuracy provided that \frac{π}{p} ≤ \frac{π}{Ωρ}. The sampling theorem requires a to be evaluated on a grid with stepsize \frac{π}{Ω}. Thus the evaluation of (3.6) needs O(Ω^{4}) operations. This number can be reduced to O(Ω^{3}) by introducing the functions

b_{j}(s) = ∑_{ℓ=−q}^{q} v_{Ω}(s − s_{ℓ})g(θ_{j}, s_{ℓ}) ,  (3.8)

evaluating these function for s = s_{k} and computing b_{j}(x · θ_{j}), which is needed in (3.7), by interpolation. Practical experience shows that linear interpolation sufficies.

(3.6) is called the filtered backprojection algorithm. The algorithm depends on the parameters p, q, Ω and the filter ϕ. We make some remarks concerning the role of these parameters.

1. Ω controls the spatial resolution of the algorithm. According to the sampling theorem, details of size 2π/Ω (and no smaller ones) can be accurately reconstructed provided that p ≥ Ωρ and q ≥ \frac{1}{π} Ωρ. Thus the number of data needed for reconstructing an essentially Ω-band-limited function in |x| ≤ ρ is essentially 2qp = \frac{2}{π} Ω^{2}ρ^{2}.

It can be shown that \frac{1}{π} Ω^{2}ρ^{2} pieces of data suffice. Without loosing resolution we can drop g(θ_{j}, s_{ℓ}) for j + ℓ = p mod 2. In that case the interpolation step for (3.7) has to be done with care, e.g. with a stepsize much smaller than ρ/q [15].

2. The condition \frac{ρ}{q} ≤ \frac{π}{Ω} has to be satisfied strictly. If it is violated, (3.5) is not valid even approximately, making the reconstruction completely unacceptable.

3. If \frac{π}{p} ≤ \frac{π}{Ωρ} is not satisfied, then the number p of directions permits artefact free reconstructions only within a circle of radius ρ′ < ρ where \frac{π}{p} = \frac{π}{Ωρ′}. This means that artefacts occur in a distance of 2ρ′ from high density objects. These can be avoided if we sacrifice resolution by choosing Ω′ < Ω such that \frac{π}{p} = \frac{π}{Ω′ρ} in the filter ϕ.

4. Only filters with vanishing kernel sum, i.e.

∑_{ℓ} v_{Ω}(s_{ℓ}) = 0

should be used.

5. Choosen p, q in an optimal way leads to p = πq. In practice this relation is not strictly observed. Usually, p is chosen smaller. According to 3. this leads to artefacts. But these artefacts are usually outside the region of interest.

Besides the filtered backprojection algorithm there exists a plethora of other algorithms which, however, are much less used. The direct Fourier algorithm is based on (3.6). Theoretically the complexity of this algorithm is O(Ω^{2} log Ω), but the accurate and efficient implementation is in no way easy. It seems that the gridding method [45] works satisfactorily, but the theory behind this method is not well understood. Direct algebraic algorithms [39] compute a minimal norm solution in L_{2} for finitely many data by FFT techniques. Iterative methods work on discrete versions of linear reconstruction problems [10], [24]. ART (= algebraic reconstruction technique [24], [25]) is based on the Kaczmarz iteration for linear systems and can be viewed as an SOR method [14]. The EM iteration is described in section 2(b).

The filtered backprojection algorithm can also be used for fan beam geometry, in which the X-ray tube sits on a circle surrounding the patient, with a detector array on the opposite side. The sampling requirements are discussed in [38].

4 Formulas for 3D reconstruction

In 3D CT one reconstructes a from the values of

g(θ, v) = ∫ a(v + tθ)dt

where v runs through the source curve V outside supp(a) and θ ∈ S^{2}. Most reconstruction formulas make use of an intermediate function

G(θ, v) = ∫_{S^{2}} g(ω, v)h(ω · θ)dω  (4.1)

where h is homogeneous of degree −2 and a function F on S^{2} × IR^{1} derived from G by

F (θ, v · θ) = G(v, θ) .  (4.2)

The reconstruction formula then reads

a(x) = \frac{1}{16π^{3}} ∫_{S^{2}} (F ⋆ k)(θ, θ · x)dθ  (4.3)

where the convolution with the 1D function k in the second argument. It can be shown that (4.1-3) is in fact a reconstruction formula provided that

ĥ(σ)k̂(σ) = σ^{2} + odd function of σ .  (4.4)

For the proof we start out from

G(θ, v) = (h ⋆ Ra)(θ, v · θ)

with the 3D Radon transform [49], hence F = h ⋆ Ra. Then, (4.3) is compared with (2.4) for n = 3, i.e.

a(x) = −\frac{1}{8π^{2}} ∫_{S^{2}} (Rf)′′(θ, x · θ)dθ .

This coincides with (4.3) if (4.4) holds.

The simplest choice for h, k is ĥ = −(−2π)^{−1/2}iσ, k̂ = (2π)^{1/2}iσ. In this case, both h = −δ′, k = 2πδ′ are local. We obtain Grangeat’s inversion formula [19]

G(θ, v) = ∫_{S^{2}∩θ^{⊥}} \frac{∂}{∂θ} g(v, ω)dω ,  (4.5)

a(x) = \frac{1}{8π^{2}} ∫_{S^{2}} F′(ω, ω · x)dω  (4.6)

where \frac{∂}{∂θ} is the derivative in direction θ with respect to the second argument and F′ is the partial derivative with respect to the second argument. In order to apply (4.5-6), F (ω, s) is needed for each plane ω · x = s meeting supp(a). Since F is obtained from G by means of (4.2) we need for each such plane a source v in that plane. In view of (4.5) this means that g is available for a neighbourhood of the fan in that plane converging to the source v. This is Grangeat’s completeness condition.

The inversion formulas of Tuy [53], B. Smith [50] and Gelfand and Goncharov [17] can be obtained by putting ĥ(σ) = (2π)^{1/2}(|σ| − σ), k̂(σ) = (2π)^{−1/2}σ and ĥ(σ) = (2π)^{1/2}|σ|, k̂(σ) = (2π)^{−1/2}|σ|, respectively [12]. These formulas are not as useful as Grangeat’s formula since h is no longer local.

In practice g(v, ·) is measured on a detector plane D_{v} = p_{v} − v + (p_{v} − v)^{⊥} where p_{v} is the orthogonal projection of v onto D_{v}. Putting g_{v}(y) = g(v, y − v), (4.1) assumes the form

G(v, θ) = |p_{v} − v| ∫_{D_{v}} g_{v}(y)h(θ · (y − v))dy .

Introducing an orthogonal system W_{v} = (w_{1}, w_{2}) in (p_{v} − v)^{⊥} we obtain in the Grangeat case (4.5)

G(v, θ) = \frac{|p_{v} − v|}{|θ_{v}|^{3}} θ_{v} · (Rg′_{v}) (\frac{θ_{v}}{|θ_{v}|}, \frac{(v − p_{v}) · θ}{|θ_{v}|} )  (4.7)

where θ_{v} = w^{T} θ, R is the 2D Radon transform and g′_{v} the gradient of g_{v}. Thus Grangeat’s formula can be implemented by computing line integrals in the detector plane, followed by a 3D backprojection (4.6). An implementation analogous to the filtered backprojection algorithm of 2D tomography can be found in [12].

In 3D emission CT, the requirements are quite different. In PET one puts the object into a vertical cylinder whose interior surface is covered by detectors. With such an arrangement one measures the X-ray transform for all lines joining two points on the mantle of the cylinder. In principle one could do the reconstruction layer by layer, using only horizontal lines in each layer. However, all the information contained in the oblique rays would be lost.

A formula which at least partially copes with this situation is

f(x) = −\frac{1}{4π^{2}} Δ ∫_{G} ∫_{θ^{⊥}} \frac{P f(θ, x − y)}{|y|L(θ, y)}dydθ  (4.8)

where G ⊆ S^{2} is a spherical zone around the equator and L(θ, y) is the length of the intersection of G and the plane spanned by θ, y [41]. With (4.8) one still has problems near the openings of the cylinder. More satisfactory reconstruction formulas based on the principle of the stationary phase have been given in [11].

5 Algorithm for more general inverse problems

The problem of CT is probably the most simple inverse problem. In general one has to solve the inverse problem for a partial differential equation. The numerical analysis of such inverse problems is still in its infancy. Of course one always can use a Newton method for solving the inverse problem simply as a nonlinear problem [31]. A method for solving the inverse problem (2.12) of ultrasound CT by an ART type method (compare section 3) is as follows [40]. For each direction θ_{j}, j = 1, …, p, define the nonlinear (Radon-type) operator R_{j} = L_{2}(Ω) → L_{2}(∂Ω)

R_{j}(f) = u|_{∂Ω_{j}}

where u is the solution to (2.12a) for θ = θ_{j} with boundary values \frac{∂}{∂ν} u(x) = \frac{∂}{∂ν} g(x, θ_{j}). Then one has to find an approximate solution f to

R_{j}(f) = g_{j} , g_{j}(x) = g(x, θ_{j}) on ∂Ω , j = 1, …, p .

This can be done by a Kaczmarz-type iteration. Starting out from an initial approximation f_{0} we put

f_{j} = f_{j−1} − ωR′_{j}(f_{j−1})^{∗}C_{j}(R_{j}(f_{j−1}) − g_{j}) , j = 1, …, p  (5.1)

with some positive definite operator C_{j} and a relaxation parameter ω. After p steps one repeats the whole process. It turns out that for practical purposes, one can take C_{j} to be the identity. In this form each step of (5.1) requires the solution of (2.12a) with boundary values as specified above (this yields R_{j}(f_{j−1})), and the solution of an adjoint boundary value problem of the same structure for the application of R′_{j}(f_{j−1})^{∗}. In the engineering literature this is known as the adjoint field method [9]. Nothing is known about convergence, and the performance - even though it is much faster than Newton-type methods - leaves much to be desired. Presently, 2D problems with a resolution 128 × 128 can be solved in a few minutes on a workstation, but 3D problems in a clinical environment are totally out of question. The same applies to equations such as (2.8) and to the equations of optical tomography (2.18).

Bibliography

1. Aben, H.K.: Integrated Photoelasticity, Valgus, Talin 1975 (Russian).

2. Allesandrini, G.: Stable Determination of Conductivity by Boundary Mesurements, Appl. Analysis, 27, 153-172, (1988).

3. Anikonov, D.S. - Prokhorov, I.V. - Kovtanyuk, E.E.: Investigation of Scattering and Absorbing Media by the Methods of X-ray Tomography, J. Inv. Ill-Posed Problems, 1, 259-281; (1993).

4. Anikonov, D.S.: Uniqueness of Simultaneous Determination of two Coefficients of the Transport Equation, Soviet Math. Dokl. 30, 149-151, (1984).

5. Arridge, S.R. - Van der Zee, P. - Cope, M. - Delpy, D.T.: Reconstruction Methods for Infrared Absorption Imaging, Poc. SPIE 1431, 204-215, (1991).

6. Barrett, H.H. - Swindell, S.: Radiological Imaging, Vol. I, II, Academic Press 1981.

7. Bondarenko, A.- Antyufeev, V.: X-Ray Tomography in Scattering Media, Institute of Mathematics, Novosibirsk, Russia (1990).

8. Bronnikov, A.V.: Degration Transform in Tomography, Pattern Recognition Letters 15, 527-592, (1994).

9. Bui, H.D.: Inverse Problems in the Mechanics of Materials. CRC Press 1994.

10. Censor, Y.: Finite Series-Expansion Reconstruction Methods, Proc. IEEE 71, 409-419, (1983).

11. Defrise, M. et.al.:Performance Study of 3D Reconstruction Algorithms for Positron Emission Tomography, International Meeting on Fully Three-Dimensional Image Reconstsruction in Radiology and Nuclear Medicine, June 23-25, 1993, Snowbird, Utah, USA.

12. Defrise, M. - Clack, R.: A Cone-Beam Reconstruction Algorithm Using Shift-Variant Filtering and Cone-Beam Backprojection, IEEE Transactions on Medical Imaging 13, 186-195 (1995).

13. Devaney, A.J.: A Filtered Backpropagation algorithm for Diffraction Tomography, Ultrasonic Imaging, 4, 336-350, (1982).

14. Elving, T.: Block-iterative methods for consistent and inconsistent linear equations, Numer. Math. 35, 1-12 (1980).

15. Faridani, A.: Reconstructing From Efficiently Sampled Data in Parallel-Beam Computed Tomography, in: Inverse Problems and Imaging, G.F. Roach, ed., Pitman Res. Notes Math. Ser. 245, 68-102, (1991).

16. Finch, D.V.: Cone Beam Reconstruction with Sources on a Curve. SIAM J. Appl. Math., 45, 665-673,(1985).

17. Gelfand, I.M. - Goncharov, A.B.: Recovery of a Compactly Supported Function Starting from Its Integrals over Lines Intersecting a Given Set of Points in Space, Doklady 290 (1986), English Translation in Soviet Math.. Doklady 34, 373-376 (1987).

18. Geman, S. - McClure, D.: Statistical Methods for Tomographic Image Reconstruction, ISI Tokio session, Bull. Int. Statist. Inst., LII(4), 5-21, 1987.

19. Grangeat, P.: Mathematical Framework of Cone Beam 3D Reconstruction via the First Derivative of Radon Transform, in: Herman et al. (eds.): Mathematical methods in tomography, Springer 1991.

20. Gratton, E. et al.: A novel approach to laser tomography, Bioimaging, 1, 40-46 (1993).

21. Green, P.J.: Bayesian Reconstructions from Emission Tomography Data Using a Modified EM Algorithm, IEEE Transactions on Medical Imaging, 9(1), 84-93, März 1990.

22. Grünbaum, F.A. - Kohn, P.D. - Latham, G.A. - Singer, J.R. - Zubelli, J.P.: Diffuse Tomography, Proc. SPIE, 1431, 232-238 (1991).

23. Gutman, S. - Klibanov, M.V.: Regularized Quasi-Newton Method for Inverse Scattering Problems, Mathl. Comput. Modelling 18, No. 1, 5-31, Pergamon Press Ltd. (1993).

24. Herman, G.T. - Lent, A.: Iterative Reconstruction Algorithms, Comput. Biol. Med. 6, 273-294, (1976).

25. Herman, G.: Image Reconstruction From Projections. The Fundamentals of Computerized Tomography. Academic Press 1980.

26. Hertle, A.: The Identification Problem for the Constantly Attenuation Radon Transform, Math. Z. 197, 13-9, (1988).

27. Hinshaw, W.S. - Lent, A.H.: An Introduction to NMR Imaging: From the Bloch Equation to the Imaging Equation, Proc. IEEE 71, 338-350 (1983).

28. Kak, A.C. - Slaney, M.: Principle of Computerized Tomography Imaging. IEEE Press 1987.

29. Kaltenbach, J.-M. - Kaschke, M.: Frequency- and time-domain modelling of light transport in random media, Technical Report, Carl Zeiss, PF 1980, Oberkochen, Germany, 1992.

30. Kleinman, R.E. - van den Berg, P.M.: A Modified Gradient Method for Two-Dimensional Problems in Tomography, J. Comp. Appl. Math., 42, 17-35, (1992).

31. Klibanov, M.V. - Gutman, S. - Barbour, R. - Chang, J. - Malinsky, J. - Alfano, R.R.: Consideration of Solutions to the Inverse Scattering Problem for Biomedical Applications, Proc. SPIE 1887, (1993).

32. Krestel, E. (ed.): Imaging Systems for Medical Diagonostics, Siemens Aktiengesellschaft, 1990.

33. Lai, C-M.: Reconstructing NMR Images from Projections Under Inhomogeneous Magnetic Field and Non-Linear Field Gradient, Phys. Med. Biol. 8, 925-938, (1983).

34. Liang, Z.-P. - Boada, F.E. - Constable, R.T. - Haacke, E.M. - Lauterbur, P.C. - Smith, M.R.: Constrained Reconstruction Methods in MR Imaging, Rev Magn. Reson. Med. 4, 67-185, (1992).

35. Louis, A.K.: Medical Imaging: State of the Art and Future Development, Inverse Problems 8, 709-738 (1992).

36. Nachman, A.I.: Global Uniqueness for a Two-Dimensional Inverse Boundary Value Problem, Department of Mathematics, Preprint Series, Number 19, University of Rochester (1993).

37. Natterer, F.: Determination of Tissue Attenuation in Emission Tomography of Optically Dense Media, Inverse Problems 9, 731-736 (1993).

38. Natterer, F.: Sampling in Fan Beam Tomography, SIAM J. Appl. Mathematics 53, 358-380 (1993).

39. Natterer, F.: The Mathematics of Computerized Tomography. Wiley-Teubner 1986.

40. Natterer, F. - Wübbeling, F.: A propagation-backpropapagation method for ultrasound tomography, Inverse Problems 11, 1225-1232 (1995).

41. Orlov, S.S.: Theory of Three Dimensional Reconstruction. II. The Recovery Operator, Sov. Phys. Crystallogr. 20, 429-433 (1976).

42. Palamodov, V.: An Inversion Method for Attenuated X-Ray Transform in Space, submitted to SIAM J. Appl. Math..

43. Romanov, V.G.: Conditional Stability Estimates for the Problem of Recovering of Absorption Coefficients and Right Hand Side of Transport Equations (Russian), to appear in Siberia Math. J.

44. Sanchez, R. - McCormick, N.J.: General Solutions to Inverse Transport Problems, J. Math. Phys. 22, 847-855, (1981).

45. Schomberg, H. - Timmer, J.: The Gridding Method for Image Reconstruction by Fourier Transformation, IEEE Transactions on Medical Imaging 14, 596-607 (1995).

46. Setzepfandt, B.: ESNM: Ein rauschunterdrückendes EM - Verfahren für die Emissionstomographie. Thesis, Fachbereich Mathematik der Universität Münster, Germany 1992.

47. Sharafutdinov, V.A.: Integral Geometry of Tensor Fields. Nauka, Novosibirsk 1993 (Russian).

48. Shepp, L.A. - Vardi, Y.: Maximum Liklihood Reconstruction for Emission Tomography, IEEE Trans. Med. Imag. 1, 113-121, (1982).

49. Smith, K.T. - Solmon, D.C. - Wagner, S.L.: Practical and Mathematical Aspects of the Problem of Reconstructing Objects From Radiographs, Bull AMS 83, 1227-1270 (1977).

50. Smith, B.D.: Image Reconstruction from Cone-Beam Projections: Necessary and Sufficient Conditions and Reconstruction Methods, IEEE Transactions on Medical Imaging 4, 14-25 (1985).

51. Sparr, G. - Stråklén, K. - Lindström, K. - Persson, W.: Doppler tomography for vector fields, Inverse Problems, 11, 1051-1061 (1995). 1, 113-121, (1982).

52. Tretiak, O.J. - Metz, C.: The Exponential Radon Transform, SIAM J. Appl. Math., 39, 341-354, (1980).

53. Tuy, H.K.: An Inversion Formula for Cone-Beam reconstruction, SIAM J. Appl. Math. 43, 546-552 (1983).