Inverting and calculating the Radon transform using the Fourier transform

You might also want to listen to our podcast on the Fourier transform.

This section is much more mathematically sophisticated than the main article. Most of the mathematics is at university level. If you want to look at examples of how tomography is used in practice, return to the main article. However, if you are feeling brave, then read on as this section contains some really lovely mathematical ideas.

One of the most useful mathematical techniques ever invented is called the Fourier transform. To motivate this, imagine that you are listening to a concert, and that you record the intensity of the sound $u(t)$ of the orchestra as a function of the time $t$. The orchestra is composed of instruments that all make sounds of a frequency $\omega$ with each such sound having an intensity $U(\omega)$. A sound of frequency $\omega$ has the mathematical expression $$e^{i\omega t},$$ where $i$ is the square root of -1. The total sound that we hear is the combination of all of these sounds and it can be expressed as an integral in the form $$u(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} U(\omega) e^{i \omega t} d\omega.$$ This integral, which links the two functions $U(\omega)$ and $u(t)$, is called the inverse Fourier transform after the French mathematician Joseph Fourier. Remarkably, it is easy to find the inverse to this process. This is called the Fourier transform and it is given by $$U(\omega) = \int_{-\infty}^\infty u(t) e^{-i \omega t} dt.$$ The Fourier transform has countless applications, ranging from telecommunications to crystallography, from speech recognition to astronomy, from radar to mobile phones, and from meteorology to archaeology. It even has important applications in cryptography. The whole signal processing industry (indeed much of the digital revolution) owes its existence to the Fourier transform.

The reason is that it links the intensity of a function to the waves that make it up, and as so much of what we do involves sound or light waves, its applications are universal. It is so important that calculating the Fourier transform was one of the first tasks given to the early computers in use in the 1960s. However, these implementations had a lot of difficulties in calculating the Fourier transform and were so slow that they were absorbing a huge amount of computing time. Roughly speaking, if you wanted to find the values of the Fourier transform at $N$ points, you had to do $N^2$ calculations. Unfortunately for accuracy you have to take large values of $N$, which means that $N^2$ is very large indeed; too large to be calculated easily.

However, in 1965 there was a remarkable breakthrough when a technique called the fast Fourier transform, or FFT was invented. This was much faster, taking a time proportional to $Nlog(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 in which it can be used to analyse images. A typical image is represented by pixels, with a pixel at the point $(x,y)$ having intensity $u(x,y)$. Just as music is characterised by waves of different frequencies, so an image can be characterised by repeating patterns - "waves" - in the two directions $x$ and $y$, with repetition frequency $\alpha$ and $\beta$ respectively. 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)dx dy.$$ 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 using the FFT. The Fourier transform of an image has many applications, including 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.$$ In our tomography example, $R(\rho, \theta)$ measures the attenuation experienced by the ray corresponding to distance $\rho$ and angle $\theta$. As we vary $\rho$, we get a picture of how the attenuation changes. The Fourier transform breaks this picture up into its constituent "waves": $r(\omega, \theta)$ measures the contribution of a repetition frequency $\omega$ to this picture. 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}ds d\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) \;\; 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)} dx dy.$$ 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 all this together we get $$r(\omega, \theta) = U(\omega \cos(\theta), \omega \sin(\theta)).$$ This splendid 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. Firstly, it gives us 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 land-mines described later on in the main article, 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 \cos(\theta), \omega \sin(\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 solving 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 a change of variable formula we obtain $$u(x,y) = \frac{1}{4\pi^2}\int \int U(\alpha, \beta) e^{i(x\alpha + y\beta)}d\alpha d\beta = \frac{1}{4\pi^2}\int \int r(\omega, \theta)e^{i(x\omega\cos(\theta)+y\omega\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\rho 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, a rather more stable algorithm that is often used in practical scanning devices.