next up previous contents
Next: Fourier reconstruction Up: Numerical Methods in Tomography Previous: MART (multiplicative algebraic reconstruction

Circular harmonic algorithms

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

eqnarray2041

One can show that

  equation2045

tex2html_wrap_inline4246 is the Chebysheff polynomial of the first kind of order tex2html_wrap_inline3288 . In principle this defines an inversion formula for R: The Fourier coefficients tex2html_wrap_inline4252 of g = Rf determine the Fourier coefficients tex2html_wrap_inline4070 of f via (5.1), hence f is determined by g.

The formula (5.1) is useless for practical calculations since tex2html_wrap_inline4246 increases exponentially with tex2html_wrap_inline3288 outside [-1,+1]. Cormack (1964) derived also a stable version of his inversion formula. It reads

eqnarray2056

where tex2html_wrap_inline4270 is the Chebysheff polynomial of the second kind of order tex2html_wrap_inline3288 . 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 tex2html_wrap_inline4276 , tex2html_wrap_inline3532 , in (2.4) we obtain

displaymath4280

The j sum is a convolution. In order to make this convolution cyclic we extend tex2html_wrap_inline3396 by putting tex2html_wrap_inline4286 , in accordance with the eveness property of the Radon transform. Then,

displaymath4288

Defining

eqnarray2088

we have

displaymath4290

This defines the circular harmonic algorithm.

Algorithm 4 (Circular harmonic algorithm for standard parallel geometry.)

Data:
The values tex2html_wrap_inline3332 , tex2html_wrap_inline3236 , tex2html_wrap_inline3240 ,

g the 2D Radon transform of f.

Step 1:
Precompute the number tex2html_wrap_inline4302 and extend tex2html_wrap_inline4304 to all tex2html_wrap_inline4306 by tex2html_wrap_inline4308 .
Step 2:
For tex2html_wrap_inline4310 , tex2html_wrap_inline3240 carry out the discrete cyclic convolutions

displaymath4314

Step 3:
For tex2html_wrap_inline4310 and tex2html_wrap_inline4318 compute

displaymath4320

Result:
tex2html_wrap_inline4322 is an approximation to tex2html_wrap_inline4324 .

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 tex2html_wrap_inline4326 operations. This is slightly more (by the factor tex2html_wrap_inline4328 ) than what is needed in the filtered backprojection algorithm. The third step needs tex2html_wrap_inline4330 additions.

Algorithm 4 can be used almost without any changes for the interlaced parallel geometry, i.e. for tex2html_wrap_inline4304 with tex2html_wrap_inline3398 odd missing (p even). One simply puts tex2html_wrap_inline4338 for tex2html_wrap_inline3402 odd and doubles tex2html_wrap_inline4342 in step 2.

Circular harmonic algorithms are also available for standard fan beam data. (2.9) with tex2html_wrap_inline4344 , tex2html_wrap_inline4346 , tex2html_wrap_inline4348 reads

displaymath4350

where g is now the fan beam data function from (2.8).

Algorithm 5 (Circular harmonic algorithm for standard fan beam geometry.)

Data:
The values tex2html_wrap_inline3596 , tex2html_wrap_inline3236 , tex2html_wrap_inline3240 with g from (2.8).
Step 1:
Precompute the numbers tex2html_wrap_inline4362 .
Step 2:
For tex2html_wrap_inline4310 , tex2html_wrap_inline3240 carry out the discrete cyclic convolutions

displaymath4368

Step 3:
For tex2html_wrap_inline4310 and tex2html_wrap_inline4372 compute the sums

displaymath4374

Result:
tex2html_wrap_inline4322 is an approximation to tex2html_wrap_inline4324 .

The complexity of Algorithm 5 is again O tex2html_wrap_inline4326 .

A few remarks are in order.

1.
Circular harmonic algorithms compute the reconstruction on a grid in polar coordinates. Interpolation to a cartesian grid (for instance for the purpose of display) is not critical and can be done by linear interpolation.
2.
The resolution of the circular harmonic algorithms is the same as for the corresponding filtered backprojection algorithms in section 2.1.
3.
Even though circular harmonic algorithms are asymptotically a little slower (by a factor of tex2html_wrap_inline4328 ) than filtered backprojection they usually run faster due to their simplicity. This is true in particular for fan beam data because in that case the backprojection is quite time consuming.
4.
Circular harmonic algorithms tend to be more accurate than filtered backprojection because no additional approximations, such as interpolations or homogeneity approximations (in the fan beam case) are used.
5.
The disadvantage of circular harmonic algorithms lies in the fact that they start with angular convolutions. This is considered as unpractical in radiological applications.


next up previous contents
Next: Fourier reconstruction Up: Numerical Methods in Tomography Previous: MART (multiplicative algebraic reconstruction

Frank Wuebbeling
Thu Sep 10 10:51:17 MET DST 1998