# Basic Physics of Nuclear Medicine/Fourier Methods

Fourier methods are described below in
quite general terms. We will see that any waveform can be broken up
mathematically into sine waves of different
amplitudes and frequencies using the Fourier
transform. We will also see that the
resulting Fourier Spectrum can be modified to enhance and/or suppress
frequencies of interest. The aim of this exercise is to treat the
Modulation Transfer Function (**MTF**) and the filtration stage of the
tomographic image reconstruction process in greater detail than
previously.

## Periodic Functions

A periodic waveform, is a function
which repeats itself regularly over a given interval of time or space. A
good example is a sine wave or a square
wave. When the waveform fluctuates with
respect to time, the wave can be characterised by its **frequency** (see
the figure below), which is defined as the number of cycles passing a
given point each second. Frequency is expressed in units of cycles per
second, or hertz (Hz). A good example of such a periodic waveform in
medicine is the electrocardiogram (ECG).

!A sine wave in the temporal
domain.

When the waveform fluctuates with respect to distance (see the following
figure), the wave is characterised by its **spatial frequency**, which
is defined as the number of cycles per unit distance, e.g. cycles per
mm.

!A sine wave in the spatial
domain.

A medical imaging example of this latter type of fluctuation is the Pb
bar pattern which is widely used to determine the spatial resolution of
an imaging system (see the following figure). In this case, the spatial
frequency is expressed in, for example, **line pairs per mm**, since
each section of Pb and its adjacent void is referred to as a **line
pair**.

!(a) An image of a test object used to determine the limiting spatial
resolution of a gamma camera. (b) An image of a test object used to
determine the limiting spatial resolution of a radiographic imaging
system. (c) A plot of the number of counts (or the radiographic density)
against distance for well resolved regions of image (a) or
(b). An image of a test object used to determine the limiting spatial resolution of a gamma camera. (b) An image of a test object used to determine the limiting spatial resolution of a radiographic imaging system. (c) A plot of the number of counts (or the radiographic density) against distance for well resolved regions of image (a) or (b)."){width="380"}

## Fourier Series

Almost all periodic functions of interest in medical imaging can be
represented by a Fourier Series. This
approach considers that any periodic waveform can be represented by the
sum of a series of sine and cosine waves as follows:

```{=html}
<div class="center">
```
$f(x) = \frac{1}{2}a_0 + a_1 \cos (x) + a_2 \cos (2x)... + a_n \cos (nx) + b_1 \sin (x) + b_2 \sin (2x)... + b_n \sin (nx)$

```{=html}
</div>
```
A square wave, for example, can thus be represented by:

```{=html}
<div class="center">
```
$f(x) = \frac{4h}{\pi} \left ( \sin (x) + \frac{1}{3} \sin (3x) + \frac{1}{5} \sin (5x) + \frac{1}{7} \sin (7x) + ... \right )\ \ (1)$

```{=html}
</div>
```
where *h* is the amplitude of the square wave.

The opposite consideration is also true mathematically, i.e. a square
wave can be constructed by adding together a large number of sine waves
of different frequencies and amplitudes. The addition of the first four
terms of equation (1) is demonstrated in the following figure. The first
term (sin *x*) is shown in panel (a), the addition of the second term to
the first term is shown in panel (b) and so on. Notice that the
fundamental, or first harmonic \[panel
(a)\], has the same frequency as the square wave and that the higher
frequencies progressively build up the shape of the square waveform
\[panels (b)-(d)\]. We can conclude that the higher frequencies
contribute to the sharpness of the sides of the square wave.

!Illustration of the addition of sine waves to approximate a square
wave.{width="380"}

## Fourier Spectrum

The Fourier Series can also be represented as a frequency
spectrum. For example, the amplitudes
of the frequency components for the square wave in equation (1) are
plotted against spatial frequency in the following figure. Note that the
Fourier Spectrum can be used to identify the frequencies and amplitudes
of the sine waves which contribute to make up a given waveform. Note
also that plots of amplitude versus distance are generally referred to
as the **spatial domain** representation, and plots of amplitude versus
spatial frequency as the **frequency domain** representation.

!The Fourier Spectrum for a square
wave.{width="280"}

## Fourier Transform

The Fourier transform is an elegant mathematical technique for
converting data from the spatial domain to the frequency domain (see the
next figure). In other words, the frequencies and amplitudes of the sine
waves which make up any waveform can be readily determined by taking the
Fourier transform of that waveform.

!The Fourier Transform can be used to generate the Fourier Spectrum of
any waveform, including the square wave
above.{width="380"}

The Fourier transform is used widely throughout medical imaging, where
its applications include:

-   determining the spatial resolution of imaging systems,
-   spatial localisation in magnetic resonance
    imaging,
-   analysis of Doppler
    ultrasound
    signals,
-   image filtering in emission and transmission computed tomography.

The **Inverse Fourier Transform** is the mathematical technique used for
coverting data in the opposite direction, i.e. from the frequency domain
to the spatial domain -- see the following figure:

!The Inverse Fourier Transform can be used to generate any waveform of
interest, including the square wave
above.{width="380"}

In conclusion, the Fourier Transform (**FT**) allows us to identify the
component sine waves of a waveform, and the Inverse Fourier Transform
(**IFT**) allows us to construct a waveform from its component sine
waves.

Finally, it should be noted that the computation of Fourier Transforms
using a digital computer is generally achieved using a special algorithm
called the Fast Fourier Transform (FFT).

## The Dirac Delta Function

An interesting case to consider from a medical imaging perspective is
the delta function:

!The delta function on the left and its Fourier Spectrum on the
right.{width="380"}

The Fourier Transform of the Dirac delta function tells us that it is
composed of an infinite number of sine waves, each of the same
amplitude. If we begin to broaden this function -- as in the following
figure -- we see that the low frequency sine waves have a high amplitude
and that the amplitudes of the sine waves decrease as the spatial
frequency increases:

!The Fourier Transform of a broadened delta
function.{width="380"}

## Modulation Transfer Function

This broadening phenomenon is similar to what happens in medical
imaging. We can consider the amplitude versus distance plot in the
previous figure to be similar to a:

-   density profile through an image of a small hole in a sheet of Pb,
    or a
-   count rate profile through an image of a point source of
    radioactivity.

This type of plot is called a Point Spread Function
(PSF) in medical imaging, and its
Fourier Spectrum is called the Modulation Transfer Function (**MTF**):

!Illustration of the MTF of an ideal and an actual imaging
system.{width="240"}

A comparison of general and imaging terms used in this field is shown in
the following table:

  Domain                  General Terminology   Imaging Terminology
  ----------------------- --------------------- ------------------------------------
  **Spatial**             Input Function        Point Spread Function (PSF)
  **Spatial Frequency**   Fourier Spectrum      Modulation Transfer Function (MTF)

In addition, the MTF may be derived from:

:\* a Line Spread Function (**LSF**)

:\* a differentiated Edge Response Function (**ERF**).

Symbolically, we can write:

```{=html}
<div class="center">
```
$I_{\text{actual}} = I_{\text{ideal}} * PSF\,\!$

```{=html}
</div>
```
where $*$ indicates a convolution operation. In other words, the actual
image is obtained when the ideal image is convolved with the PSF of the
imaging system.

To restore the ideal image, all that is required theoretically is to
remove the effect of the PSF (e.g. attempts to solve imaging problems
encountered by the Hubble space
telescope). This can be achieved readily
using the Fourier transform because the convolution process in the
spatial domain is equivalent to a multiplication process in the spatial
frequency domain, i.e.

```{=html}
<div class="center">
```
$FT(I_{\text{actual}}) = FT(I_{\text{ideal}}) \cdot FT(PSF)\,\!$

```{=html}
</div>
```
and therefore,

```{=html}
<div class="center">
```
$FT(I_{\text{ideal}}) = \frac{FT(I_{\text{actual}})}{FT(PSF)}$

```{=html}
</div>
```
The complete restoration process is referred to as a
deconvolution operation and is given by:

```{=html}
<div class="center">
```
$I_{\text{ideal}} = IFT \left \lbrace \frac{FT(I_{\text{actual}})}{FT(PSF)} \right \rbrace$

```{=html}
</div>
```
## Filtering the Fourier Spectrum

The image restoration process discussed above is an example of Fourier
spectrum filtering. In other words, once a Fourier Spectrum has been
generated for an image, it can be filtered so that certain spatial
frequencies can be modified, e.g. enhanced or suppressed. This filtered
spectrum can then be inverse transformed to generate a filtered image
with, for example, sharpened or smoothed features -- as illustrated in
the following figure:

!The Fourier filtering
process.{width="280"}

An example is shown in the following figure in the case of a bone scan
of the hand. Its 2D FFT is shown in the top right panel. The result of
applying a filter which suppresses spatial frequencies below 0.05
cycles/pixel and above 0.3 cycles/pixel, i.e. a bandpass of 3-20 pixels,
is shown in the centre row of the figure, while that of a filter which
suppresses frequencies below 0.01 and above 0.1 cycles/pixels in shown
on the bottom row.

!Illustration of bandpass filtration to selectively suppress high and
low spatial frequencies and preserve intermediate
frequencies.{width="380"}

A feature we need to consider in more detail before proceeding is the
spatial frequency nature of the image data itself. Remember that images
are generally sampled digitally using a square matrix composed of
pixels, the size of which dictate how well a digital image approximates
its analogue counterpart. The resultant digital spatial resolution
places a limit on the maximum spatial frequency that can be
accommodated. The criterion usually applied in digital imaging is based
on the Nyquist-Shannon Sampling Theorem
which implies that:

```{=html}
<div class="center">
```
`<big>`{=html}When an image has spatial frequency components with a
maximum spatial frequency, f, then the image data should be sampled with
a sampling frequency of at least twice f for faithful
reproduction.`</big>`{=html}

```{=html}
</div>
```
This sampling frequency is commonly called the Nyquist
frequency. At lower sampling
frequencies, the resultant digital images can contain artefactual
patterns, called Moiré patterns, and the
phenomenon is sometimes referred to as
aliasing.

## Filtered Back Projection

The streaking inherent in the simple backprojection process makes the
actual image appear as if it had been combined mathematically with a 1/r
function, where `<big>`{=html}r`</big>`{=html} is the radial distance in
the Fourier domain. In filtered backprojection, Fourier filtering can be
used to remove the effect of this **
blurring**.

Symbolically, a measured projection can be considered to be the result
of a convolution with a blurring function:

```{=html}
<div class="center">
```
$P_{\text{measured}} = P \cdot \frac{1}{r}$

```{=html}
</div>
```
The first stage of the filtration process is to calculate the Fourier
transform of the measured projection data, i.e.

```{=html}
<div class="center">
```
$FT(P_{\text{measured}}) = FT(P) \cdot FT \left ( \frac{1}{r} \right )$

```{=html}
</div>
```
The corrected projection, *P*, is then obtained by dividing the FT of
the measured projection by **FT()**, and
taking the inverse transform of the result, i.e.

```{=html}
<div class="center">
```
$P = IFT \left \lbrace \frac{FT(P_{\text{measured}})}{FT \left ( \frac{1}{r}\right )} \right \rbrace$

```{=html}
</div>
```
The **FT()** function is simply that of a
ramp. In other words, when the **FT** of the measured projection is
calculated and the result is divided by this ramp function, then the
corrected projection can be obtained by calculating the **IFT** of this
quotient.

In addition, if some variation is introduced into the
**FT()** function, it is possible to
simultaneously correct for this blurring effect and to enhance or
suppress features in the back-projected image. For example, the blurring
artifact can be removed and, at the same time:

:\* the fine detail can be enhanced (as in the so-called **bone
algorithms** in X-ray CT), or

:\* the image noise can be suppressed (as in the so-called **soft tissue
algorithms** in X-ray CT).

In general, variation in the ramp function can be achieved by
multiplying it by a second function, e.g. a Butterworth
function used in SPECT:

```{=html}
<div class="center">
```
$\text{Amplitude} = \frac{1}{1+(\frac{f}{f_c})^{2n}}$

```{=html}
</div>
```
where:

-   *f*~*c*~: cut-off frequency, which defines the frequency at which
    the amplitude drops by 50%, and
-   *n*: the order of the function.

!A ramp function and a Butterworth function (of variable order and
cut-off frequency) are multiplied to form the Fourier filter used in the
FBP
process. are multiplied to form the Fourier filter used in the FBP process."){width="480"}

Characteristics of a number of filters used in SPECT imaging are shown
in the following table:

+----------------------+----------------------+----------------------+
| Filter               | Equation             | Comment              |
+======================+======================+======================+
| Ram-Lak              | $\text{Amp           | star artifact        |
|                      | litude} = 1\ \,\text | removal; noise       |
|                      | {for}\ f < f_c\,\!$\ | sensitivity          |
|                      | \                    |                      |
|                      | $\text{Am            |                      |
|                      | plitude} = 0\ \,\tex |                      |
|                      | t{for}\ f > f_c\,\!$ |                      |
+----------------------+----------------------+----------------------+
| Butterworth          | $\text{Amplit        | noise reduction      |
|                      | ude} = \frac{1}{1+(\ |                      |
|                      | frac{f}{f_c})^{2n}}$ |                      |
+----------------------+----------------------+----------------------+
| Metz                 | $\text{Amplitude} =  | noise reduction;     |
|                      | \frac{1-\big(1-MTF(f | contrast enhancement |
|                      | )^2\big)^x}{MTF(f)}$ |                      |
+----------------------+----------------------+----------------------+
| Wiener | de} = \frac{1}{MTF(f | contrast enhancement |
|                      | )} \left \lbrace \fr |                      |
|                      | ac{MTF(f)^2}{S+MTF(f |                      |
|                      | )^2} \right \rbrace$ |                      |
+----------------------+----------------------+----------------------+
| Scramp               | $\text{Am            | noise reduction;     |
|                      | plitude} = 1 + \frac | contrast enhancement |
|                      | {f}{f_c} \lbrace MTF |                      |
|                      | (f_c)-1 \rbrace\ \,\ |                      |
|                      | text{for}\ f < f_c$\ |                      |
|                      | \                    |                      |
|                      | $\text{Am            |                      |
|                      | plitude} = 0\ \,\tex |                      |
|                      | t{for}\ f > f_c\,\!$ |                      |
+----------------------+----------------------+----------------------+
| Inverse MTF          | $                    | noise reduction;     |
|                      | \text{Amplitude} = \ | contrast suppression |
|                      | frac{1}{MTF(f)}\ \,\ |                      |
|                      | text{for}\ f < f_c$\ |                      |
|                      | \                    |                      |
|                      | $\text{Am            |                      |
|                      | plitude} = 0\ \,\tex |                      |
|                      | t{for}\ f > f_c\,\!$ |                      |
+----------------------+----------------------+----------------------+
| Hamming              | $\text{Amplitude     | noise reduction      |
|                      | } = 0.54 + 0.46 \cos |                      |
|                      |  \left ( \pi \frac{f |                      |
|                      | }{f_c} \right )\ \,\ |                      |
|                      | text{for}\ f < f_c$\ |                      |
|                      | \                    |                      |
|                      | $\text{Am            |                      |
|                      | plitude} = 0\ \,\tex |                      |
|                      | t{for}\ f > f_c\,\!$ |                      |
+----------------------+----------------------+----------------------+
| Parzen               | $\text{Ampl          | noise reduction      |
|                      | itude} = 2 \left \lb |                      |
|                      | race 1 - \left (\fra |                      |
|                      | c{f}{f_c} \right )^3 |                      |
|                      |  \right \rbrace\ \,\ |                      |
|                      | text{for}\ f < f_c$\ |                      |
|                      | \                    |                      |
|                      | $\text{Am            |                      |
|                      | plitude} = 0\ \,\tex |                      |
|                      | t{for}\ f > f_c\,\!$ |                      |
+----------------------+----------------------+----------------------+
| Shepp-Logan          | $\                   | noise reduction      |
|                      | text{Amplitude} = \f |                      |
|                      | rac{2f_c}{\pi f} \si |                      |
|                      | n \left (\pi \frac{f |                      |
|                      | }{f_c} \right )\ \,\ |                      |
|                      | text{for}\ f < f_c$\ |                      |
|                      | \                    |                      |
|                      | $\text{Am            |                      |
|                      | plitude} = 0\ \,\tex |                      |
|                      | t{for}\ f > f_c\,\!$ |                      |
+----------------------+----------------------+----------------------+
| Hann                 | $\text{Amplitude}    | noise reduction      |
|                      |  = \frac{f}{2} \left |                      |
|                      |  (1+ \cos \pi \frac{ |                      |
|                      | f}{f_c} \right )\ \, |                      |
|                      | \text{for} f < f_c$\ |                      |
|                      | \                    |                      |
|                      | $\text{Am            |                      |
|                      | plitude} = 0\ \,\tex |                      |
|                      | t{for}\ f > f_c\,\!$ |                      |
+----------------------+----------------------+----------------------+

The choice of filter for a given image reconstruction task is generally
a compromise between the extent of noise reduction and fine detail
suppression (and of contrast enhancement in some cases) as well as the
spatial frequency pattern of the image data of interest.

## Further reading

-   Engineering Analysis/Multi-Dimensional Fourier
    Series
-   The Feynman lectures Chapter 50
    Harmonics

