Fast Fourier Transforms

In 1965 there was a remarkable breakthrough when a technique called the Fast Fourier Transform or FFT was invented to evaluate the Fourier Transform. This was much faster, taking a time proportional to N \mbox{ log}{(N)}, which is a lot smaller than N^2.  With the FFT available to calculate the Fourier transform quickly, its applications became almost unlimited.

One application, of great importance to us, is the way that it can be used to analyse images. A typical image is represented by pixels so that a pixel at the point (x,y) has intensity u(x,y). The Fourier Transform of such an image is then given by the double-integral.

U(\alpha,\beta) = \int\!\!\!\int e^{-i(\alpha x + \beta y)} u(x,y)\ dxdy.

This transformation also has an inverse which is given by

u(x,y) = \frac{1}{4\pi^2} \int\!\!\!\int e^{i(\alpha x + \beta y)} U(\alpha,\beta)\ d\alpha d\beta.

Both the Fourier Transform of an image and its inverse can be calculated quickly by using the FFT. The Fourier transform of an image also has many applications. Some of the most important of these are the removal of noise from an image, and deblurring a blurred image. However, what is of most importance to us is a link between the Fourier and the Radon transform.

To find this link we need to fix \theta and find the Fourier Transform r(\omega,\theta) of the Radon transform. This is given by

r(\omega,\theta) = \int R(\rho,\theta)e^{-i\omega\rho} d\rho.

If we then substitute in the expression for the Radon Transform we have

r(\omega,\theta) - \int\!\!\!\int u(\rho \cos(\theta) - s \sin(\theta), \rho \sin(\theta) + s \cos(\theta))e^{-i \omega \rho}\ dsd\rho.

Now, we can change variables in this integral by setting

(x,y) = (\rho \cos(\theta) - s \sin(\theta), \rho \sin(\theta) + s \cos(\theta)),
 

so that

\rho = x\cos(\theta) + y\sin(\theta)\mbox{ and }dxdy = dsd\rho

This then gives

r(\omega,\theta) = \int\!\!\!\int u(x,y)e^{-i(\omega \cos(\theta)x\, +\, \omega \sin(\theta)y)}\ dxdy

But, this is a very familiar integral. It is precisely the two-dimensional Fourier Transform of u(x,y) evaluated at a particular point. Putting this all together we get

r(\omega,\theta) = U(\omega \cos(\theta), \omega \sin(\theta))

This formula is called the Fourier Slice Theorem. It tells us that the Fourier Transform of the function u(x,y) is the same as its Fourier Transform evaluated at a particular point.

There are many ways that we can use this formula. First, it provides a quick way to calculate the Radon Transform of the function u. Of course, in medical imaging, the Radon Transform is "calculated" automatically by measuring the attenuation of the X-Rays through the body. However, in other applications such as the detection of landmines described later, it is vital to calculate the Radon Transform as quickly as possible. We can do this by using a combination of the FFT and the Fourier Slice Theorem as follows

  • Calculate the Fourier Transform U(\alpha,\beta) of u(x,y) using the FFT.
  • Set (\alpha,\beta) = (\omega \cos(\theta), \omega \sin(\theta)) and use the FFT to calculate the inverse Fourier Transform of the function r(\omega,\theta) to find R(\rho,\theta).

As each stage of this method uses the FFT it is a lot quicker than calculating R directly by performing a lot of integrals.

The second important application is that the Fourier Slice Theorem gives us a way of inverting the Radon Transform, so that we can find the function u(x,y) if we know R. In particular, substituting the formula for the inverse Fourier Transform into the Fourier Slice Theorem and applying the change of variables, we obtain

u(x,y) = \frac{1}{4\pi^2} \int\!\!\!\int U(\alpha,\beta) e^{i(\alpha x + \beta y)}\ d\alpha d\beta = \frac{1}{4\pi^2} \int\!\!\!\int r(\omega,\theta) e^{i \omega (x \cos(\theta) + y \sin(\theta))}\ \omega\ d\omega d\theta

Or, putting everything together we get the inversion formula

u(x,y) = \frac{1}{4\pi^2} \int\!\!\!\int\!\!\!\int R(\rho,\theta) e^{i\omega(x \cos(\theta) + y \sin(\theta)-\rho)} \omega\ d\omega d\theta


Bingo! Now knowing R we can find u every time. Well not quite. Like all inverse problems in tomography and other applications, the formula only works well if we know R very accurately and there is not too much noise in the results. Furthermore, there is quite a lot of work to do in calculating all of the terms in this integral, and errors can accumulate quite easily. In practice, this formula can be a bit unstable, and it is hard to implement accurately. However, it gives the basis for all other formulae used to find u from R. Indeed, one way to interpret the inversion formula is to note that x \cos(\theta) + y \sin(theta) - \rho = 0 is the formula for the straight line at angle \theta and distance \rho from the origin and the formula is combining the contribution of R along each of these lines. Developing this link leads to the filtered back projection method which is a rather more stable algorithm and is often used in practical scanning devices.