LII: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).

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:

$$

 
   
     
       f
       (
       x
       )
       =
       
         
           1
           2
         
       
       
         a
         
           0
         
       
       +
       
         a
         
           1
         
       
       cos
       ⁡
       (
       x
       )
       +
       
         a
         
           2
         
       
       cos
       ⁡
       (
       2
       x
       )
       .
       .
       .
       +
       
         a
         
           n
         
       
       cos
       ⁡
       (
       n
       x
       )
       +
       
         b
         
           1
         
       
       sin
       ⁡
       (
       x
       )
       +
       
         b
         
           2
         
       
       sin
       ⁡
       (
       2
       x
       )
       .
       .
       .
       +
       
         b
         
           n
         
       
       sin
       ⁡
       (
       n
       x
       )
     
   
   {\displaystyle 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)}
 

A square wave, for example, can thus be represented by:

$$

 
   
     
       f
       (
       x
       )
       =
       
         
           
             4
             h
           
           π
         
       
       
         (
         
           sin
           ⁡
           (
           x
           )
           +
           
             
               1
               3
             
           
           sin
           ⁡
           (
           3
           x
           )
           +
           
             
               1
               5
             
           
           sin
           ⁡
           (
           5
           x
           )
           +
           
             
               1
               7
             
           
           sin
           ⁡
           (
           7
           x
           )
           +
           .
           .
           .
         
         )
       
        
        
       (
       1
       )
     
   
   {\displaystyle 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)}
 

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.

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.

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.

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.

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.

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.

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.

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:

$$

 
   
     
       
         I
         
           actual
         
       
       =
       
         I
         
           ideal
         
       
       ∗
       P
       S
       F
       
       
     
   
   {\displaystyle I_{\text{actual}}=I_{\text{ideal}}*PSF\,\!}
 

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.

$$

 
   
     
       F
       T
       (
       
         I
         
           actual
         
       
       )
       =
       F
       T
       (
       
         I
         
           ideal
         
       
       )
       ⋅
       F
       T
       (
       P
       S
       F
       )
       
       
     
   
   {\displaystyle FT(I_{\text{actual}})=FT(I_{\text{ideal}})\cdot FT(PSF)\,\!}
 

and therefore,

$$

 
   
     
       F
       T
       (
       
         I
         
           ideal
         
       
       )
       =
       
         
           
             F
             T
             (
             
               I
               
                 actual
               
             
             )
           
           
             F
             T
             (
             P
             S
             F
             )
           
         
       
     
   
   {\displaystyle FT(I_{\text{ideal}})={\frac {FT(I_{\text{actual}})}{FT(PSF)}}}
 

The complete restoration process is referred to as a deconvolution operation and is given by:

$$

 
   
     
       
         I
         
           ideal
         
       
       =
       I
       F
       T
       
         {
         
           
             
               F
               T
               (
               
                 I
                 
                   actual
                 
               
               )
             
             
               F
               T
               (
               P
               S
               F
               )
             
           
         
         }
       
     
   
   {\displaystyle I_{\text{ideal}}=IFT\left\lbrace {\frac {FT(I_{\text{actual}})}{FT(PSF)}}\right\rbrace }
 

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.

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.

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:

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.

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 r is the radial distance in the Fourier domain. In filtered backprojection, Fourier filtering can be used to remove the effect of this ¹⁄_r blurring.

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

$$

 
   
     
       
         P
         
           measured
         
       
       =
       P
       ⋅
       
         
           1
           r
         
       
     
   
   {\displaystyle P_{\text{measured}}=P\cdot {\frac {1}{r}}}
 

The first stage of the filtration process is to calculate the Fourier transform of the measured projection data, i.e.

$$

 
   
     
       F
       T
       (
       
         P
         
           measured
         
       
       )
       =
       F
       T
       (
       P
       )
       ⋅
       F
       T
       
         (
         
           
             1
             r
           
         
         )
       
     
   
   {\displaystyle FT(P_{\text{measured}})=FT(P)\cdot FT\left({\frac {1}{r}}\right)}
 

The corrected projection, P, is then obtained by dividing the FT of the measured projection by FT(¹⁄_r), and taking the inverse transform of the result, i.e.

$$

 
   
     
       P
       =
       I
       F
       T
       
         {
         
           
             
               F
               T
               (
               
                 P
                 
                   measured
                 
               
               )
             
             
               F
               T
               
                 (
                 
                   
                     1
                     r
                   
                 
                 )
               
             
           
         
         }
       
     
   
   {\displaystyle P=IFT\left\lbrace {\frac {FT(P_{\text{measured}})}{FT\left({\frac {1}{r}}\right)}}\right\rbrace }
 

The FT(¹⁄_r) 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(¹⁄_r) 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:

$$

 
   
     
       
         Amplitude
       
       =
       
         
           1
           
             1
             +
             (
             
               
                 f
                 
                   f
                   
                     c
                   
                 
               
             
             
               )
               
                 2
                 n
               
             
           
         
       
     
   
   {\displaystyle {\text{Amplitude}}={\frac {1}{1+({\frac {f}{f_{c}}})^{2n}}}}
 

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.

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

Filter	Equation	Comment
Ram-Lak	Amplitude = 1   for   f < f c {\displaystyle {\text{Amplitude}}=1\ \,{\text{for}}\ f<f_{c}\,\!} $$ Amplitude = 0   for   f > f c {\displaystyle {\text{Amplitude}}=0\ \,{\text{for}}\ f>f_{c}\,\!}	star artifact removal; noise sensitivity
Butterworth	$$ Amplitude = 1 1 + ( f f c ) 2 n {\displaystyle {\text{Amplitude}}={\frac {1}{1+({\frac {f}{f_{c}}})^{2n}}}}	noise reduction
Metz	$$ Amplitude = 1 − ( 1 − M T F ( f ) 2 ) x M T F ( f ) {\displaystyle {\text{Amplitude}}={\frac {1-{\big (}1-MTF(f)^{2}{\big )}^{x}}{MTF(f)}}}	noise reduction; contrast enhancement
Wiener	$$ Amplitude = 1 M T F ( f ) { M T F ( f ) 2 S + M T F ( f ) 2 } {\displaystyle {\text{Amplitude}}={\frac {1}{MTF(f)}}\left\lbrace {\frac {MTF(f)^{2}}{S+MTF(f)^{2}}}\right\rbrace }	noise reduction; contrast enhancement
Scramp	Amplitude = 1 + f f c { M T F ( f c ) − 1 }   for   f < f c {\displaystyle {\text{Amplitude}}=1+{\frac {f}{f_{c}}}\lbrace MTF(f_{c})-1\rbrace \ \,{\text{for}}\ f<f_{c}} $$ Amplitude = 0   for   f > f c {\displaystyle {\text{Amplitude}}=0\ \,{\text{for}}\ f>f_{c}\,\!}	noise reduction; contrast enhancement
Inverse MTF	Amplitude = 1 M T F ( f )   for   f < f c {\displaystyle {\text{Amplitude}}={\frac {1}{MTF(f)}}\ \,{\text{for}}\ f<f_{c}} $$ Amplitude = 0   for   f > f c {\displaystyle {\text{Amplitude}}=0\ \,{\text{for}}\ f>f_{c}\,\!}	noise reduction; contrast suppression
Hamming	Amplitude = 0.54 + 0.46 cos ⁡ ( π f f c )   for   f < f c {\displaystyle {\text{Amplitude}}=0.54+0.46\cos \left(\pi {\frac {f}{f_{c}}}\right)\ \,{\text{for}}\ f<f_{c}} $$ Amplitude = 0   for   f > f c {\displaystyle {\text{Amplitude}}=0\ \,{\text{for}}\ f>f_{c}\,\!}	noise reduction
Parzen	Amplitude = 2 { 1 − ( f f c ) 3 }   for   f < f c {\displaystyle {\text{Amplitude}}=2\left\lbrace 1-\left({\frac {f}{f_{c}}}\right)^{3}\right\rbrace \ \,{\text{for}}\ f<f_{c}} $$ Amplitude = 0   for   f > f c {\displaystyle {\text{Amplitude}}=0\ \,{\text{for}}\ f>f_{c}\,\!}	noise reduction
Shepp-Logan	Amplitude = 2 f c π f sin ⁡ ( π f f c )   for   f < f c {\displaystyle {\text{Amplitude}}={\frac {2f_{c}}{\pi f}}\sin \left(\pi {\frac {f}{f_{c}}}\right)\ \,{\text{for}}\ f<f_{c}} $$ Amplitude = 0   for   f > f c {\displaystyle {\text{Amplitude}}=0\ \,{\text{for}}\ f>f_{c}\,\!}	noise reduction
Hann	Amplitude = f 2 ( 1 + cos ⁡ π f f c )   for f < f c {\displaystyle {\text{Amplitude}}={\frac {f}{2}}\left(1+\cos \pi {\frac {f}{f_{c}}}\right)\ \,{\text{for}}f<f_{c}} $$ Amplitude = 0   for   f > f c {\displaystyle {\text{Amplitude}}=0\ \,{\text{for}}\ f>f_{c}\,\!}	noise reduction

Note that x in the equation for the Metz filter is called the power factor, and that S in the Wiener filter equation is called the shape factor.

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

Notes

This article is a direct transclusion of the Wikibooks article and therefore may not meet the same editing standards as LIMSwiki.