We already made use of the relation
for the Radon transform in in section 2.1, and of the corresponding formula for the nD ray transform
in section 3.4. In this section we use these formulas to derive reconstruction algorithms. To fix ideas we consider the case of the 2D Radon transform, sampled as in the standard parallel geometry. This means that g = Rf is given for ,
,
and
,
, as in
section 2.1. f is assumed to vanish outside
.
The idea of Fourier reconstruction is very simple: Do a 1D Fourier transform on g with respect to the second variable for each . According to (6.1) this yields
in all of
. Do a 2D inverse Fourier
transform to obtain f. Even though this seems fairly obvious, the
numerical implementation in a discrete setting is quite intricate. In fact good
Fourier algorithms have been found only quite recently.
To begin with we describe the simplest possible implementation. We warn the reader that this algorithm is quite useless since it is not sufficiently accurate.
Algorithm 6 (Standard Fourier reconstruction.)
g the 2D Radon transform of f.
The algorithm is designed to reconstruct a function f with support in which is essentially
-band-limited. (2.6), (2.7)
have to be satisfied. We stress again that the algorithm as it stands is not to be recommended because of poor accuracy. Better versions are described below.
A few comments are in order. In step 1 we compute an approximation to
Under the assumption (2.6) this approximation is reliable since Fourier transforms are evaluated exactly by the trapezoidal rule if the Nyquist condition, in our case (2.6), is satisfied. According to (6.1), step 1 provides us with the values
In step 2 we compute on the cartesian grid
by nearest neighbour interpolation:
Since f vanishes outside ,
has bandwidth
. Thus
sampling of
on a 2D grid with stepsize
is adequate
by the sampling theorem.
Step 3 is the trapezoidal rule for the 2D inverse Fourier transform,
properly discretized and complying with the Nyquist condition. Hence is an approximation to
.
Step 1 and step 3 of Algorithm 6 are justified by the sampling theorem. Thus the failure of the algorithm must be caused by the interpolation in step 2. This is in fact the case. Of course we can replace the interpolation by a more accurate one, such as linear interpolation. However this doesn't help much.
Inspite of its poor accuracy, Fourier reconstruction is attractive because of its favourable complexity which is due to the fast Fourier transform (FFT).
We have used FFT in the circular harmonic algorithm already, but FFT is so essential for Fourier reconstruction that we say a few words about FFT; for a thourough treatment we refer to Nussbaumer (1982). The discrete Fourier transform of length q is defined to be
An algorithm which computes from
in less then O
, typically O
, operations is a called an FFT. In the circular harmonic algorithm we have used FFT just for the evaluation of (6.3). In Fourier reconstruction
we employ FFT for the evaluation of Fourier integrals
for the functions f in with support in
. Sampling theory tells us
that
has to be discretized with stepsize
(
has bandwidth
). With
the stepsize for f the trapezoidal rule provides the approximation
The range of k is in agreement with the sampling theorem: The stepsize corresponds to the bandwidth
, hence
in (6.4) suffices. Of course (6.4) is a discrete Fourier transform
of length 2q. Sometimes one wants to adjust the stepsizes of f and
differently. Then one has to evaluate
with an arbitrary real parameter c. This can be done by the chirp-z-algorithm, see
Nussbaumer (1982), again with typically O operations.
Assuming that we have a fast Fourier transform (FFT) algorithm which does a discrete Fourier transform of length q with O operations, the complexity of Algorithm 6 is as follows. Step 1 does
p Fourier transforms of length 2q, requiring O
operations. Assuming that the interpolation in step 2 can be done in O(1) operations per
point we get O
as the complexity of step 2. The 2D Fourier transform in
step 3 can be done with O
operations. Hence the complexity of
Algorithm 6 is O
. This is much better than the filtered backprojection algorithm (Algorithm 1), which needs O
operations for a
reconstruction on a
grid.
Presently there exist two Fourier methods with satisfactory accuracy and favourable complexity.
1. The linogram algorithm (Edholm and Herman (1987)).
Here, interpolation in Fourier space is avoided altogether by a clever choice of the
directions . The linogram algorithm works on the data
where j, . Doing a 1D Fourier transform on
results in
where arccot (j/q). For
we get from
(6.1)
Note that this can be done efficiently by the chirp-z-algorithm. The key observation is that the points
form a grid lying on vertical lines with distance , being evenly spaced within
each vertical line (though with different stepsizes in different lines). On such a grid
we can do a 1D FFT in the horizontal direction in the usual way. In the horizontal
direction the stepsize is not what we need for a direct application of the FFT, but
the chirp-z-algorithm is still applicable. This takes care of the 2D inverse Fourier transform in
. For
we proceed analogously
with the data
, evaluating
for
. We remark that the linogram data in Edholm and Herman (1987) is a little different from
ours, namely
,
, respectively.
The use of these detector positions makes the linogram algorithm even simpler in that the
factor
in the right hand side of (6.6) dissapears.
2. The gridding algorithm (O'Sullivan (1985), Kaveh (1987), Schomberg and Timmer (1995).
This algorithm works on the standard parallel data which we used in Algorithm 6. It does the interpolation in step 2 of Algorithm 6 in the following way. Let w be a smooth function in with w = 1 on
which is decaying exponentially at infinity. Put
. Then,
by (6.1). Using a quadrature rule with nodes ,
and weights
we obtain the approximation
to . The method relies on the following assumptions:
If these conditions are met then of (6.3) is a good approximation to
which can be evaluated essentially in O(1) operations for each k. This takes care of step 2. Of course when using (6.7) we have to divide f by w after step 3 to make up for the previous multiplication with w.
It is needless to say that our derivation of the gridding algorithm is purely heuristic. It seems that presently there exists no convincing theoretical analysis of the gridding algorithm.