CT Scans and the Radon Transform

Last December, Dublin’s Tallaght Hosptal acquired a new CT scanner, a Toshiba Aquilon Prime, the first of its type in the country. The state-of-the-art scanner is housed in a room with a ‘sky ceiling’ that allows patients to enjoy an attractive outdoor image during the scanning process.

This equipment, which cost €600,000 will undoubtedly result in timely treatment of patients and the saving of lives. The process of generating images from CT scans is described in the latest That’s Maths column (TM016) in the Irish Times.

Toshiba Aquilion Prime CT Scanner installed at Tallaght Hospital in 2012.

CT Scanning

Modern medicine depends on non-invasive imaging techniques that enable us to see inside the body without the risks of exploratory surgery. One of the most valuable methods of “probing our innards” is computer tomography or CT scanning.

To take a CT scan, the patient lies on a table that slides through a hole in a large annular apparatus. An x-ray source rotates around the annulus, and x-rays passing through the patient are detected on the opposite side. This gives a series of readings from many angles and these can be analysed to form a cross-sectional image. Many cross-sections can be combined to form a 3D image.

CT stands for “computed tomography”. The word tomography comes from the Greek tomos, or slice, and a CT scan is made by combining x-ray images of cross-sections or slices through the body. From these, a 3-D representation of internal organs can be built up.

In a CT scan, multiple x-ray images are taken from different directions. The x-ray data are then fed into a tomographic reconstruction program to be processed by a computer. The deduction of the tissue structure from the x-rays is done using a technique first obtained in 1917 by Johann Radon.

Radon, an Austrian mathematician, was studying the mathematical properties of the operation that we now call the Radon transform. He was motivated by purely theoretical interest, and could not have anticipated the great utility of his work in the practical context of CT. Reconstruction techniques have grown in sophistication, but are still founded on Radon’s work.

X-ray Absorption

The total attenuation, or dampening, of an x-ray as it passes through the body tissue is expressed as a “line integral”, the sum of the absorptions along the path of the beam. The more tissue along the path and the denser this tissue, the less intense the x-ray beam becomes. The challenge is to determine the patterns of normal and abnormal tissue from the x-ray patterns.

Aquilion PRIME 160 CT scan showing a brain tumour. The arterial supply and venous drainage can be seen in this 3D image from http://medical.toshiba.com/

The intensity of radiation passing through a homogeneous body of material decreases exponentially with distance,

$\displaystyle I(x) = I_0 \exp(-\mu x)$

where ${I_0}$ is the intensity at ${x=0}$ and ${\mu}$ is the attenuation or absorption coefficient. Differentiating this, we have

$\displaystyle \frac{dI}{dx} = -\mu I$

giving the rate of attenuation with distance. This equation holds more generally where the absorption coefficient varies with distance ${x}$ , that is, for a non-homogeneous body. Integrating in this more general case gives

$\displaystyle I(x) = I_0 \exp\left(-\int\mu\,dx\right)$

If we consider material confined to an interval ${[a,b]}$ with ${0\le a<b\le L}$ , then

$\displaystyle I_L = I(L) = I_0 \exp\left(-\int_a^b \mu\,dx\right)$

So, given ${I_L}$ , we can calculate the total absorption

$\displaystyle M = \int_a^b \mu\,dx = -\log\left(\frac{I_L}{I_0}\right) = \log\left(\frac{I_0}{I_L}\right)$

But from ${I_0}$ and ${I_L}$ we can say nothing more about the distribution of the material in the interval ${[a,b]}$ .

In two dimensions, things change. Given the total absorption for every cross-section through the body, we can construct the absorption coefficient ${\mu(x,y)}$ as a function of position. This was first shown by Johann Radon [Radon, 1917].

The attenuation of an x-ray beam is measured by comparison with the radiation absorbed by water, and is expressed in Houndsfield units, after one of the developers of the CT scanner. Water has a value of 0 Houndsfield units, air is -1000 units and bone is about 1000 units. On an x-ray image, areas with low unit values are black and those with high values are white. Thus, bone shows up as a bright region, the organs are various shades of grey and cavities are black.

The intensity of each x-ray beam gives information about the absorption along a single line through the body. To reconstruct the structure completely, we need x-rays along every line through the body. This is impossible but, from a large set of lines, we can come close to an exact picture.

Any line ${\cal L}$ in the ${xy}$ -plane can be specified by its perpendicular distance ${r}$ from the origin and the angle ${\theta}$ of the perpendicular. Then any point ${P}$ on the line is given by

$\displaystyle P(x,y) = ( r\cos\theta - s\sin\theta , r\sin\theta + s\cos\theta )$

with ${s}$ varying along the line. If we consider the absorption of an x-ray beam along a line ${{\cal L}}$ , we must integrate with respect to ${s}$ . The result will depend on ${r}$ and ${\theta}$ :

$\displaystyle M(r,\theta) = \int_{-\infty}^{+\infty} \mu(r\cos\theta-s\sin\theta,r\sin\theta+s\cos\theta) \,ds$

It is this quantity that is measured directly by an x-ray scanner.

The Radon transform is a function of the polar coordinates ${(r,\theta)}$ . It is a linear operation with respect to the function ${\mu(x,y)}$ being transformed. A graph of ${M(r,\theta)}$ with ${r}$ and ${\theta}$ on orthogonal cartesian axes is called a sinogram. It depicts all the data from a 2D CT scan.

Filtered back projection

How do we recover ${\mu(x,y)}$ from ${M(r,\theta)}$ ? The answer was found in 1917 by the Austrian mathematician Johann Radon (1887–1956) who investigated the properties of the transform and found an inversion method.

The first step is to compute the average of all the line integrals passing through a given point ${(x,y)}$ :

$\displaystyle {\cal B}M(x,y) = \frac{1}{\pi}\int_0^\pi M(x\cos\theta+y\sin\theta,\theta)\,d\theta$

This is called the back projection of ${M}$ . It does not give us ${\mu(x,y)}$ but rather a smoothed version of ${\mu}$ .

The link between the Fourier and Radon transforms is expressed by the Central Slice Theorem [Feeman, 2010]. We denote the 1D and 2D Fourier transforms by ${{\cal F}_1}$ and ${{\cal F}_2}$ and the Radon transform by ${{\cal R}}$ . Then the 1D and 2D Fourier transforms are related by

$\displaystyle {\cal F}_2 \mu(\omega\cos\theta,\omega\sin\theta) = {\cal F}_1({\cal R}\mu)(\omega,\theta)$

To “undo” the smoothing effect of the back projection, the Radon transform is subjected to a filtering procedure in which high frequencies are boosted. This involves a Fourier transform, followed by multiplication by the (absolute value of) frequency, followed by an inverse Fourier transform. The “filtered back projection” then becomes

$\displaystyle \mu(x,y) = \frac{1}{2} {\cal B} \{ {\cal F}^{-1} [|\omega| {\cal F}({\cal R}\mu)(\omega,\theta) ] \} (x,y)$

This is the required formula for inversion of the Radon transform.

Here is the recipe: given the Radon transform ${{\cal R}\mu(r,\theta)}$ , a function of polar coordinates,

Take the 1D Fourier transform with respect to the radial variable ${r}$
Multiply each component by the absolute value of frequency ${|\omega|}$
Take the inverse Fourier transform
Take the back projection.

This leads to the recovery of the original function ${\mu(x,y)}$ .

Filtered back projection is the fundamental method for reconstruction of the x-ray absorption function. However, it assumes complete knowledge of the Radon transform. In practice, we only have a finite number of x-ray cross-sections. Therefore, only an approximation to the function ${\mu(x,y)}$ is possible.

Multiplication by ${|\omega|}$ boosts high frequencies so, if noise is present, it can cause numerical instabilities and spoil the CT image. To avoid this, we replace ${|\omega|}$ by a band-limited function that vanishes for ${|\omega|}$ greater than some cut-off value ${\omega_c}$ . Numerous refinements of the techniques for generating CT scans are described in the book of Feeman [Feeman, 2010].

Further Information

Dean, S. R., 2007: The Radon Transform and some of its Applications. Dover, 295pp. ISBN: 978-0-486-46241-7.

Feeman, T. G., 2010: The Mathematics of Medical Imaging. Springer, 141pp. ISBN: 978-0-387-92711-4.

Radon, J, 1917: On the determination of functions from their integrals along certain manifolds. Translation in [Deans, 2007].

Related

Published by thatsmaths