Circular harmonic algorithms can be derived from the inversion formula of Cormack (1963). For the 2D problem Rf = g it is obtained by Fourier expansions
One can show that
is the Chebysheff polynomial of the first kind of order
.
In principle this defines an inversion formula for R: The Fourier coefficients
of
g = Rf determine the Fourier coefficients
of f via (5.1), hence f
is determined by g.
The formula (5.1) is useless for practical calculations since increases exponentially with
outside [-1,+1]. Cormack (1964) derived also a stable version of his inversion formula. It reads
where is the Chebysheff polynomial of the second kind of order
. This formula does not suffer from exponential increase. It is the starting point of the circular harmonic reconstruction algorithms of Hansen (1981), Hawkins and Barrett (1986), Chapman and Cary (1986).
We take a different route and start out from (2.1) again. We consider only the
case n = 2. Putting ,
, in (2.4) we obtain
The j sum is a convolution. In order to make this convolution cyclic we extend by putting
,
in accordance with the eveness property of the Radon transform. Then,
Defining
we have
This defines the circular harmonic algorithm.
Algorithm 4 (Circular harmonic algorithm for standard parallel geometry.)
g the 2D Radon transform of f.
The second step of the algorithm has to be done with a fast Fourier transform (FFT)
in order to make the algorithm competitive with filtered backprojection. In that case step 2 requires O operations. This is slightly more (by the factor
) than what is needed in the filtered backprojection algorithm. The third step needs
additions.
Algorithm 4 can be used almost without any changes for the interlaced parallel geometry,
i.e. for with
odd missing (p even). One simply puts
for
odd and doubles
in step 2.
Circular harmonic algorithms are also available for standard fan beam data. (2.9) with ,
,
reads
where g is now the fan beam data function from (2.8).
Algorithm 5 (Circular harmonic algorithm for standard fan beam geometry.)
The complexity of Algorithm 5 is again O .
A few remarks are in order.