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 
is sampled at
The reconstruction region is 
.
We need some tools from sampling theory. For a function f in 
we define the Fourier transform by
We call a function (essentially) 
-band-limited if 
vanishes (approximately) for 
. Shannon's sampling theorem states that an -band-limited function can be recovered exactly from the samples f(hk), 
, if 
. If 
,
both are -band-limited, then
The convolution of two functions is denoted by 
, independently of the dimension.
Rather than using (2.6) we start out from
if 
. Here, 
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 
, the Dirac 
-function. One can show that this is the case iff
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
where 
is a low-pass filter, i.e. 
for 
. For the ideal low-pass 
, 
, we have
with the 
-function 
. 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
The convolution on the right hand side of (3.3) is now approximated by the trapezoidal rule
At this stage we assume a to be essentially -band-limited. Since
where 
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 
, as can be seen from (3.2). For the backprojection we use the trapeziodal rule again, obtaining as approximation to a(x)
It can be shown that 
as a function of 
, 
, is essentially 
-band-limited. Thus, (3.7) is satisfied with good accuracy provided that 
.
The sampling theorem requires a to be evaluated on a grid with stepsize 
. Thus the evaluation of (3.6) needs O 
operations. This number can be reduced to O 
by introducing the functions
evaluating these function for 
and computing 
, 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 
(and no smaller ones) can be accurately reconstructed provided that 
and 
. Thus the number of data needed for reconstructing an essentially -band-limited function in 
is essentially 
.
It can be shown that 
pieces of data suffice. Without loosing resolution we can drop 
for 
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 
[15].
2. The condition has to be satisfied strictly. If it is violated, (3.5) is not valid even approximately, making the reconstruction completely unacceptable.
3. If is not satisfied, then the number p of directions permits artefact free reconstructions only within a circle of radius 
where 
. This means that artefacts occur in a distance of 
from high density objects. These can be avoided
if we sacrifice resolution by choosing 
such that 
in the filter .
4. Only filters with vanishing kernel sum, i.e.
should be used.
5. Choosen p, q in an optimal way leads to 
. 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 
, 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 
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].