HTI 5720

Digital Imaging and PACS

Fourier Transform

Lawrence WC Chan, Ph.D.

1
 Definition of digital image

• Spatial representation of colors or

gray levels
• 2D array of pixel values (brightness
of 1 or 3 components)
• Pixel value could be integer or
decimal (bits per pixel determine
precision)
2
 Grayscale image
Pixel value at row x and column y
= f (x , y) Origin
(0,0) y

3
 Color image
Pixel values at row x and column y
= [ fR(x , y) , fG(x , y) , fB(x , y) ]T
Origin
(0,0) y

4
1D Representations of digital images
Surface

Horizontal Vertical Diagonal

Line intensity profiles

Jan J. 2006 page 4 5

 Line intensity profile
Pixel value Horizontal
28 30 27 28
35
30
25
28 29 30 1
20
15
27 26 4 6 10
5
0
26 5 2 3 1 2
Pixel
3 4

Pixel value 45 degree

Pixel value Vertical 35
30
30
25 25
20 20
15 15
10 10
5 5

0 0
1 2 3 4 1 2 3 6 4
Pixel Pixel
 Line intensity profile: 1D signal

Any digital Image, f(x,y), can

be regarded as an integration
of 1D signals, f(x), in all
directions.

50

00

50

00

7
50
0 50 100 150 200 250
 What is function “f” in 1D?
Function “f” can be any well-defined function.
e.g. harmonic (cosine) 1D signal, f(x) = cos (x)

f(x)

x
8
Sinusoidal function:

f(x) = A*sin(2fx+b)

Meanings of:
(1) frequency, f
(2) amplitude, A
(3) phase, b
Sinusoidal function
Frequency concept in signal

Low frequency signal High frequency signal

Amplitude concept in signal

Low amplitude signal High amplitude signal

 Phase concept in signal

Zero phase angle 90 degree (1.57 radian) phase angle

Sine wave Cosine wave
Frequency concept in signal
Frequency: 0.05 Hz Frequency: 0 Hz
Amplitude: 50 Amplitude: 50
Phase angle: 0 degree Phase angle: 0 degree

Signal: y=A*sin(2ft) Signal: y=50*sin(0)

Frequency concept in signal
Frequency: 0.05 Hz Frequency: 0 Hz
Amplitude: 50 Amplitude: 50
Phase angle: 90 degree Phase angle: 90 degree

Signal: y=A*cos(2ft) Signal: y=50*cos(0)

 Superimposition of sinusoidal
functions
Example 1: two sine waves and one constant

16
 Superimposition of sinusoidal
functions
Example 2: Seven sine waves and one constant

17
Rectangular wave
approximation
 1 Fundamental

1+2

1+2+3

• Any signal can be decomposed into series of sine and cosine waves
• Using Fourier transform to determine the composition
• The period of the original determines the fundamental frequency
 Fundamental frequency

The smallest non-zero frequency

among the sinusoidal waves
decomposed from the original
signal.
 Superimposition of sinusoidal
functions
Example 3: 1000 sine waves and a constant

21
 What do these examples tell us?

Integration (summation) of constant and harmonic

signals in particular frequencies, phase and
amplitude can yield any signal (well-defined
function).

In the example 3, addition of 1000 different sine

waves and a constant can yield (almost) a
squared wave.

Equivalently, a squared wave is composed of

sinusoidal waves of different frequencies
and a constant (zero frequency signal). 22
 Frequency representation

Any mathematical function can be written as

a combination of sinusoidal (harmonic)
functions with frequencies u (2f), phase
 and amplitudes A.

Laere K.V. et al. 2001 23

 Squared Wave as an Example

Decomposition into sine waves of different

frequencies and amplitudes:

24
 Approximation of squared wave

The number of terms determines the goodness-of-fit to

the desired function, which improves as i grows larger.

Laere K.V., et al. 2001 25

 Any signal can be represented by
superimposing sinusoidal (sine
and cosine) waves.

For example, sawtooth wave

sin(2ft ) + sin(4ft ) + sin(6ft ) + sin(8ft ) + 

1 1 1
2 3 4
Fundamental freq.
Superimposition of sine waves
Fourier Series Representation
Cosine wave Sine wave
Fourier transform decomposes any
signal in a series of sine and cosine
waves with different frequencies,
amplitudes and phases.

Frequency spectrum is then formed.

 Example
Temporal domain
Frequency domain
July Swimming Jogging Freq per day Swimming Jogging

1 30 mins Once per 10 days 0.1

2 20 mins Twice per 10 days 0.2 30 mins

3 3 times per 10 days 0.3
4 20 mins
4 times per 10 days 0.4
5
6 30 mins 20 mins 5 times per 10 days 0.5 20 mins
7
6 times per 10 days 0.6
8 20 mins
9 7 times per 10 days 0.7

10 20 mins
8 times per 10 days 0.8
11 30 mins
9 times per 10 days 0.9
12 20 mins
13 10 times per 10 days 1
14 20 mins
15
16 30 mins 20 mins
 Graphical representation
Temporal domain
Frequency domain
Swimming schedule

40
Exercise Pattern
30
minutes

20 35
10
30
0
25
0 2 4 6 8 10 12 14 16 18

minutes
20
day
15
10
Jogging schedule
5
25
20 0
minutes

15 0 0.2 0.4 0.6 0.8 1 1.2

10 Freq (times per day)
5
0
0 2 4 6 8 10 12 14 16 18
day
 Weight variation
Temporal domain
Frequency domain
Weight variation due to swimming
Weight variation pattern

64 70
63
62
Weight (kg)

61 60
60
59 50
58

Weight (kg)
57 40
56
0 5 10 15 20 30

Day
20

10

Swimming schedule
0
40 0 0.05 0.1 0.15 0.2 0.25
30 Freq (cycle per day)
minutes

20
10
0
0 2 4 6 8 10 12 14 16 18 Two frequency components with
day different amplitudes
Break

33
1D Fourier Transform (FT)

A mathematical technique to transform

spatial to frequency information

F(u) = FT [ f(x) ]

34
 FT of a constant

f(x) |F(u)|

x u

35
 FT of 5 components

f(x)

x
|F(u)|

f(x)
u

x 36
 FT of 15 components

f(x)

x
|F(u)|

f(x)
u

x 37
Inverse 1D Fourier Transform
A mathematical technique to transform
frequency to spatial information

f(x) = FT-1 [ F(u) ]

= F(u)

38
 Information in spatial domain (mm)

f(x)
x → signal at x

Information in frequency domain (mm-1)

F(u)
u → (A, ) at u

where F(u) gives a complex number, A=|F(u)| and =F(u)

39
 2D Fourier Transform

Jan J. 2006 page 9

Inverse 2D Fourier Transform

Jan J. 2006 page 13

40
2D Digital Fourier Transform

Inverse 2D Digital Fourier

Transform

Jan J. 2006 pages 94,95 41

 Information in spatial domain (mm)

f ( x, y )
(x,y) → pixel value at (x,y)

Information in frequency domain (mm-1)

F ( u, v )
(u,v) → (A, ) at (u,v)

where F(u,v) gives a complex number, A=|F(u,v)| and

=F(u,v)
42
Frequency concept in image

Low frequency pattern High frequency pattern

Waveform Waveform
propagating propagating
horizontally vertically
 Fourier Transform

Fourier Transform

Spatial domain Frequency domain

(0,0) Frequency in ky
Column y direction

kx
(0,0)
Row Frequency in
represents the origin (0,0) x direction
 Frequency domain:
What are meant by high freq and low freq?
Frequency in
vertical direction v
Relatively high freq
(10,10)

Relatively low freq (-5,5) (5,5) Relatively low freq

u
Frequency in
Relatively low freq (-4,-4)
horizontal direction

(8,-8) Relatively high freq

Determined by the
closeness to origin 45
 X-ray image and
its amplitude and phase spectra

f(x,y)

|F(u,v)| F(u,v)

Jan J. 2006 page 12 46

Ripple direction and freq.

47
Ripple direction
f(x,y) |F(u,v)|

Jan J. 2006 page 102 48

f(x,y) |F(u,v)|

49
Remove ripple

50
f(x,y) |F(u,v)|

“Brick-wall” filter

51
Rectangular Window
f(x,y) |F(u,v)|

Jan J. 2006 page 102 52

f(x,y) |F(u,v)|

Jan J. 2006 page 11 53

Ripple in Rectangular Window

f(x,y) |F(u,v)|

Jan J. 2006 page 102 54

f(x,y) |F(u,v)|

55
 Zero-frequency component
f(x,y) F(u,v)

Relatively low freq

56
Rel. high freq (other pixels) (white pixel at origin)
 Horizontal frequency component
f(x,y) F(u,v)

57
Rel. high freq Rel. low freq
 Vertical frequency component
f(x,y) F(u,v)

58
Rel. high freq Rel. low freq
 3x3 centered at (0,0)
f(x,y) F(u,v)

59
Rel. high freq Rel. low freq
 5x5 centered at (0,0)
f(x,y) F(u,v)

60
Rel. high freq Rel. low freq
 9x9 centered at (0,0)
f(x,y) F(u,v)

61
Rel. high freq Rel. low freq
 17x17 centered at (0,0)
f(x,y) F(u,v)

62
Rel. high freq Rel. low freq
 33x33 centered at (0,0)
f(x,y) F(u,v)

63
Rel. high freq Rel. low freq
 65x65 centered at (0,0)
f(x,y) F(u,v)

64
Rel. high freq Rel. low freq
Excluding high freq.

“Brick-wall”
filter

65
Adjustment of |F(u,v)|

66
Break

67
Fourier Transform of convolution integral
G=HxF

What do you think about the meaning of

convolution filter h?
g=h*f

68
 FT of Convolution Integral

Imply |G(u,v)| = |H(u,v)| |F(u,v)|

|H(u,v)| specifies the degree of amplification or attenuation of the

particular harmonic component by the system.

Jan J. 2006 pages 25,26 69

 Low-pass Filtering

f(x,y) |F(u,v)|

g(x,y) |H(u,v)|

Jan J. 2006 page 27 70

 High-pass Filtering

f(x,y) |F(u,v)|

g(x,y) |H(u,v)|

Jan J. 2006 page 27 71

 Narrowband Filtering

f(x,y) |F(u,v)|

g(x,y) |H(u,v)|

Jan J. 2006 page 27 72

Which filter is applied?
 Design H in frequency domain
(1) Inverse FT

Ideal h in continuous spatial domain

(2) Discretization

Truncated h in discrete spatial domain

(3) Convolution with image f

Filtered image, g
75
 Fourier filter, H

(1) Inverse FT (2) Truncate to smaller matrix

Spatial filter, h

Jan J. 2006 page 598 76

 Example: Ramp filter
Ramp filter in frequency domain
1
Three samples only if the range
0.8 is quantized into levels of 0.1.
1
0.6

0.4

0.2
0.5
0
-pi -pi/2 0 pi/2 pi
Spatial Freq
FT-1
Ramp filter in spatial domain
1 0
0.8

0.6

0.4

0.2 -0.5
54 56 58 60 62 64 66 68
0

-0.2

-0.4
20 40 60 80 100 120
77
Spatial location
 Example: 2D Ramp filter
H: 201x201 matrix h: 201x201 matrix

FT-1

21x21 matrix at the centre

Spatial ramp filter 78
 Two Approaches
Original image, f

(1) FT

FT of original image, F

(2) G = H x F g=h*f

FT of filtered image, G

(3) Inverse FT

Fourier filtered image, g Spatial filtered image, g 79

Fourier Filtering

80
 Fourier Transform
f: 201x201 matrix F: 201x201 matrix

FT

81
 Fourier filtering: G = H x F
H: 201x201 matrix F: 201x201 matrix

G: 201x201 matrix

82
 Inverse FT
G: 201x201 matrix g: 201x201 matrix

FT-1

83
Spatial Filtering

84
 Spatial filtering: g = h* f
f: 201x201 matrix h: 21x21 matrix

g: 201x201 matrix

85
Comparison of two approaches

Fourier filtering Spatial filtering

86
 Four filters

w1 = 3x3 impulse

w2 = 3x3 averaging filter

w3 = 3x3 high-pass filter = w1 - w2

w4 = 3x3 sharpening filter = w1 - 0.9 w2

 3x3 Filter

w1 = 3x3 impulse

0 0 0
0 1 0
0 0 0

f = f * w1
 3x3 Filter

w2 = 3x3 averaging filter

1/9 1/9 1/9

1/9 1/9 1/9
1/9 1/9 1/9

g = f * w2
Low frequency components retained
High frequency components removed
 3x3 Filter

w3 = 3x3 high-pass filter = w1 - w2

0 0 0 1/9 1/9 1/9

0 1 0 - 1/9 1/9 1/9
0 0 0 1/9 1/9 1/9

g = f * w3 = f * (w1 – w2)
= f * w1 – f * w2 = f – f * w2
High frequency components retained
Low frequency components removed
 3x3 Filter

w4 = 3x3 sharpening filter = w1 - 0.9 w2

0 0 0 1/9 1/9 1/9

0 1 0 - 0.9 1/9 1/9 1/9
0 0 0 1/9 1/9 1/9

g = f * w4 = f * (w1 – 0.9 w2)

= f * w1 – 0.9 f * w2 = f – 0.9 f * w2
High frequency components retained
90% Low frequency components removed
Aliasing

92
(pixels/cycle)

93
 Nyquist frequency (Nq)
The maximum frequency of signal, can be
measured using sampling, is 0.5 cycles/pixel.
(fixed sampling frequency)

Laere K.V., et al. 2001

This value corresponds to the minimum

sampling rate, 2 pixels/cycle, which is the
critical value for aliasing.
(fixed signal frequency)

94
Sine wave

95
Aliasing Artifacts

96
 Insufficient Sampling

(a) original (well sampled image, (b) half of density of sampling, (c) quarter
density of sampling, and (d) quarter density of sampling after image
smoothing (i.e., after simple low-pass filtering).
97
Jan J. 2006 page 60
 Prevent Aliasing

The sampling frequency should be greater

than two times of the frequency of signal.

Nyquist frequency is half of the sampling

frequency, determining which frequency
components can be reconstructed in
digitization process.

Jan J. 2006 pages 58, 59 98

 Line pairs
Squared wave is the intensity profile of line pairs.
One cycle of squared wave represents one bright and one dark lines.

Intensity profile: squared wave of frequency u

Sine wave component of lowest frequency u

sin(2πu x)

1 cycle = 1 lp
99
Line pairs and sine wave

Source: [Link] 100

iBlur: 19,128,314,262
Modulation Transfer Function (MTF)

101
Source: [Link]
 Frequency components
For a fixed sampling frequency,
which component can be reconstructed?

|F(u)|

u (cycle/mm)
Nyquist frequency (cycle/mm)
= 0.5 (cycle/pixel) * Sampling frequency (pixel/mm)
102
Thank You!

Good Luck!

103