Digital Image Processing

Image Transforms
Unitary Transforms and the 2D Discrete Fourier Transform
DR TANIA STATHAKI
READER (ASSOCIATE PROFFESOR) IN SIGNAL PROCESSING
IMPERIAL COLLEGE LONDON
 What is this lecture about?

• Welcome back to the Digital Image Processing lecture!

• In this lecture we will learn about the Discrete Fourier Transform in
images. This is a two-dimensional transform (2D DFT).
• We will briefly recall the one-dimensional Discrete Fourier Transform (1D
DFT) and show how this transform can be extended into two dimensions.
• The 1D DFT is a member of a large family of transforms called unitary
transforms.
One-dimensional transforms in Science and
Engineering.
One dimensional unitary
transforms
Two- One-
dimensional dimensional
DFT DFT
 The generic form of a one-dimensional signal transform

• Consider a one-dimensional discrete signal with 𝑁 samples.

• Call the signal 𝑓(𝑥), 0 ≤ 𝑥 ≤ 𝑁 − 1.
• A transform of the signal 𝑓(𝑥) will “convert” this signal into a new signal
which let us call 𝑔(𝑢).
• It makes sense that if 𝑓(𝑥) has 𝑁 samples, then 𝑔(𝑢) has 𝑁 samples.
• In engineering and science most of the transforms we deal with have the
following generic form:
𝑔 𝑢 = 𝑁−1 𝑥=0 𝑇 𝑢, 𝑥 𝑓(𝑥), 0 ≤ 𝑢 ≤ 𝑁 − 1, 𝑇 𝑢, 𝑥 is a function of 𝑢, 𝑥 called
the forward transformation kernel.
• Note that all values of 𝑓(𝑥) are required to produce a single value of 𝑔 𝑢 .
• We can put the signals 𝑓(𝑥) and 𝑔 𝑢 in column vectors. Denote these
with 𝑓 and 𝑔, respectively .
• Therefore, the transform can be written in the form 𝑔 = 𝑇 ∙ 𝑓.
• The matrix 𝑇 is of dimension 𝑁 × 𝑁 and contains the values 𝑇 𝑢, 𝑥 for
different 𝑢, 𝑥 .
 Example of a one-dimensional signal transform: DFT

• Recall the generic form:

𝑔 𝑢 = 𝑁−1𝑥=0 𝑇 𝑢, 𝑥 𝑓(𝑥), 0 ≤ 𝑢 ≤ 𝑁 − 1, 𝑇 𝑢, 𝑥 is a function of 𝑢, 𝑥.
• Now recall the one-dimensional Discrete Fourier Transform (DFT)
𝑁−1 1 −𝑗2𝜋𝑢𝑥
𝐹 𝑢 = 𝑥=0 𝑓(𝑥) 𝑁 𝑒 𝑁,0≤𝑢 ≤𝑁−1
Very often 𝑥 represents time. We can write:
𝑓(0)
1 𝑢 𝑢(𝑁−1) 𝑓(1)
𝐹 𝑢 = −𝑗2𝜋𝑁 −𝑗2𝜋 ∙
𝑁 1 𝑒 … 𝑒 𝑁 ⋮
𝑓(𝑁 − 1)
• We see that for each 𝑢 the signal is projected (observe the inner product)
onto complex sinusoids which depend on 𝑢 (the frequency of the complex
sinusoid).
• Try to observe the analogy with the generic form: for the DFT we have
1 −𝑗2𝜋𝑢𝑥
𝑇 𝑢, 𝑥 = 𝑒 𝑁.
𝑁
 The generic form of the inverse
of a one-dimensional signal transform
• Transforming a signal 𝑓(𝑥) into a signal 𝑔(𝑢) might be handy but we
must have a method to get back to 𝑓(𝑥) from 𝑔(𝑢).
• Therefore, we need an inverse transform.
• The inverse transform has the generic form:
𝑓 𝑥 = 𝑁−1 𝑢=0 𝐼 𝑢, 𝑥 𝑔(𝑢), 0 ≤ 𝑥 ≤ 𝑁 − 1, 𝐼 𝑢, 𝑥 is a function of 𝑢, 𝑥 called
the inverse transformation kernel.
• Note again that all values of 𝑔(𝑢) are required to produce a single value
of 𝑓 𝑥 .
• The inverse transform can be written in the form 𝑓 = 𝐼 ∙ 𝑔.
• The matrix 𝐼 is of dimension 𝑁 × 𝑁 and contains the values 𝐼 𝑢, 𝑥 for
different 𝑢, 𝑥 .
• Note that 𝑓 = 𝐼 ∙ 𝑔 = 𝐼 ∙ 𝑇 ∙ 𝑓 ⇒ 𝐼 = 𝑇 −1 . We must invert 𝑇 to obtain 𝐼.
• Requiring to invert a matrix is a tedious task. Think that if 𝑁 = 1000, i.e.,
𝑓(𝑥) is a signal with 1000 samples, then 𝑇 is of dimension 1000 × 1000.
• There must be a way to do this easily, otherwise using a transform would
not be efficient.
 The one-dimensional inverse DFT (IDFT)
• Recall the generic form:
𝑓 𝑥 = 𝑁−1𝑢=0 𝐼 𝑢, 𝑥 𝑔(𝑢), 0 ≤ 𝑥 ≤ 𝑁 − 1, 𝐼 𝑢, 𝑥 is a function of 𝑢, 𝑥.
• Now recall the one-dimensional Inverse Discrete Fourier Transform (IDFT).
𝑢𝑥
𝑁−1 𝑗2𝜋 𝑁
𝑓 𝑥 = 𝑢=0 𝑒 𝐹(𝑢),0≤𝑥 ≤𝑁−1
• Observe that for the IDFT we have
𝑢𝑥
𝑗2𝜋 1 𝑗2𝜋𝑢𝑥 1 −𝑗2𝜋
𝑢𝑥
1 −𝑗2𝜋𝑢𝑥 ∗
𝐼 𝑢, 𝑥 = 𝑒 𝑁 =𝑁 𝑒 𝑁 = 𝑁 (𝑒 𝑁 )∗ = 𝑁( 𝑒 𝑁𝑇 ∗ (𝑢, 𝑥).
𝑁) =
𝑁 𝑁 𝑁
Therefore, 𝐼 = 𝑇 = 𝑁𝑇 ∗ .
−1

• The above relationship reveals that the inverse of matrix 𝑇 is simply its
complex conjugate multiplied with 𝑁, i.e., the size of the signal.
1
• The inverse of matrix 𝑇 would be identical to 𝑇 if we put in both forward
𝑁
1
and inverse transform instead of and 1.
𝑁
• Therefore, we don’t really have to calculate any inverse matrix in the case
of DFT. This is a universal observation: in all transforms we use in
engineering and science, the inverse of the transformation matrix is easily
obtained from the forward transformation matrix.
 Unitary matrices: some brief revision of linear algebra
• A square matrix 𝑇 is unitary if 𝑇 𝐻 = 𝑇 −1 , where 𝑇 𝐻 denotes the Hermitian
𝑇
(conjugate transpose) of 𝑇, 𝑇 𝐻 = 𝑇 ∗ and 𝑇 −1 is the matrix inverse.
• Unitary matrices leave the length of a complex vector unchanged, i.e.,
𝑇𝑥 = 𝑥 .
• For real matrices, 𝑇 𝐻 = 𝑇 𝑇 = 𝑇 −1 and the term unitary matrix becomes the
same as the term orthogonal matrix.
• The term orthogonal is quite obvious since 𝑇 𝑇 = 𝑇 −1 implies that 𝑇 𝑇 ∙ 𝑇 = 𝐈
which means that the rows/columns of 𝑇 are orthogonal (and more
accurately orthonormal).
• In fact, the rows 𝑟𝑖 of a unitary matrix are a unitary basis. That is, each row
𝑇
has length one, and their Hermitian inner product is zero, i.e., 𝑟𝑖 ∙ 𝑟𝑗 ∗ = 0, 𝑖 ≠
𝑇
𝑗 and 𝑟𝑖 ∙ 𝑟𝑖 ∗ = 1.
• Similarly, the columns are also a unitary basis.
• Given any unitary basis, the matrix whose rows are that basis is a unitary
matrix. It is shown that the columns are another unitary basis.
• The eigenvalues of a unitary matrix have the property 𝜆 = 1.
 Let us go back to the one-dimensional DFT
• Now recall again the one-dimensional IDFT
𝑢𝑥
𝑁−1 𝑗2𝜋 𝑁
𝑓 𝑥 = 𝑢=0 𝑒 𝐹(𝑢),0≤𝑥 ≤𝑁−1
and the result we have shown, i.e., 𝐼 = 𝑇 −1 = 𝑁𝑇 ∗ .
𝑇
• In that case𝑇 −1
∙𝑇=𝐈= 𝑁𝑇 ∗
∙𝑇 = 𝑁𝑇 ∗ ∙ 𝑇 (observe that 𝑇 and 𝑇 ∗ are
∗ ∗ 𝑇
symmetric and therefore, 𝑇 = 𝑇 .
𝑇 1
• Therefore, 𝑇 ∗ ∙ 𝑇 = 𝐈. We see that the rows/columns of 𝑇 are not
𝑁
orthonormal but they are orthogonal. This property is equally good to
make the transformation matrix easily invertible; the inverse is simply the
conjugate, multiplied by a constant factor.
• However, why are we still discussing one-dimensional DFT (1D-DFT) in
an Image Processing course?
• The 1D-DFT is the building block for the implementation of the two-
dimensional DFT introduced next.
 Why do we use image transforms?

• Before we embark into the two-dimensional DFT (2D DFT) let us comment
in general on Image Transforms. Why do we use them?
• Often, image processing tasks are best performed in a domain other than
the spatial domain. The key steps are:
o Transform the image.
o Carry the task(s) of interest in the transformed domain.
o Apply inverse transform to return to the spatial domain.
 The generic form of a two-dimensional transform
• Consider a two-dimensional signal 𝑓(𝑥, 𝑦) 𝑥 = 0, … , 𝑀 − 1 and 𝑦 =
0, … , 𝑁 − 1. We can assume that 𝑓(𝑥, 𝑦) is an image.
• The generic form of a two-dimensional (image) transform is a
straightforward extension from the one-dimensional case. It is defined as
follows:
𝑀−1 𝑁−1

𝑔 𝑢, 𝑣 = 𝑇 𝑢, 𝑥, 𝑣, 𝑦 𝑓(𝑥, 𝑦)
𝑥=0 𝑦=0
• A 2D transform is separable if 𝑇 𝑢, 𝑥, 𝑣, 𝑦 = 𝑇1 (𝑢, 𝑥) ∙ 𝑇2 (𝑣, 𝑦).
• A 2D transform is symmetric if 𝑇1 𝑢, 𝑥 = 𝑇2 (𝑢, 𝑥)
• The inverse transform is defined as follows:
𝑀−1 𝑁−1

𝑓 𝑥, 𝑦 = 𝐼 𝑢, 𝑥, 𝑣, 𝑦 𝑔(𝑢, 𝑣)
𝑢=0 𝑣=0
• The same definitions for separability and symmetry apply to the inverse
transform.
 Separable and symmetric image transforms
• Consider the generic image transform below. Assume that 𝑥 indicates row
and 𝑦 indicates column.
𝑀−1 𝑁−1

𝑔 𝑢, 𝑣 = 𝑇 𝑢, 𝑥, 𝑣, 𝑦 𝑓(𝑥, 𝑦)
𝑥=0 𝑦=0
• Suppose that it is separable and symmetric:
𝑀−1 𝑁−1

𝑔 𝑢, 𝑣 = 𝑇 𝑢, 𝑥 𝑇(𝑣, 𝑦)𝑓(𝑥, 𝑦)
𝑥=0 𝑦=0
𝑀−1 𝑁−1
• We can write 𝑔 𝑢, 𝑣 = 𝑥=0 𝑇 𝑢, 𝑥 𝑦=0 𝑇(𝑣, 𝑦)𝑓(𝑥, 𝑦)
• The sum 𝑁−1𝑦=0 𝑇(𝑣, 𝑦)𝑓(𝑥, 𝑦), 𝑣 = 0, … , 𝑁 − 1 is basically the corresponding
one-dimensional transform applied along the 𝑥 −th row of the image.
• By applying the one-dimensional
transform in each row 𝑥 = 0, … , 𝑀 − 1
we obtain an intermediate image 𝐹 𝑥, 𝑣 .
 Separable and symmetric image transforms cont.

• We can now write 𝑔 𝑢, 𝑣 = 𝑀−1 𝑥=0 𝑇 𝑢, 𝑥 𝐹(𝑥, 𝑣)

• The above sum is now the corresponding one-dimensional transform
applied along the 𝑣 −th column of the intermediate image 𝐹(𝑥, 𝑣).
• We see, therefore, that the implementation of a separable and symmetric
transform in an image requires the sequential implementation of the
corresponding one-dimensional transform row-by-row and then column-by-
column (or the inverse).
• The above process is depicted in the figure below.
 Two crucial properties of most signal and image transforms
• Energy preservation. A transform possesses the property of energy
preservation if both the original and the transformed signal have the same
2 2
energy, i.e., 𝑔 = 𝑓 . For images this property can be written as
𝑀−1 𝑁−1 2 𝑀−1 𝑁−1
𝑥=0 𝑦=0 𝑓(𝑥, 𝑦) = 𝑢=0 𝑣=0 𝑔(𝑢, 𝑣) 2 .
• Energy compaction. Most of the energy of the original data is
concentrated in only a few transform coefficients, which are placed close to
the origin; remaining coefficients have small values.
• The property of energy compaction facilitates the compression of the
original image. A generic compression scheme based on this idea can be
summarized in the following steps.
o Transform the image.
o Keep a small fraction of the transformed image values close to the origin (the
(0,0) point). That way you save space.
o For reconstruction, replace the values that you discarded with zeros and use
the inverse transform.
 Two-dimensional Discrete Fourier Transform (2D DFT)
• Consider an image 𝑓(𝑥, 𝑦) of size 𝑀 × 𝑁. The 2D DFT is defined as
follows:
𝑀−1 𝑁−1
1 𝑢𝑥 𝑣𝑦
−𝑗2𝜋 𝑀 + 𝑁
𝐹 𝑢, 𝑣 = 𝑓(𝑥, 𝑦) 𝑒
𝑀𝑁
𝑥=0 𝑦=0
• The 2D Inverse DFT (2D IDFT) is defined as:
𝑀−1 𝑁−1
𝑢𝑥 𝑣𝑦
𝑗2𝜋 +
𝑓(𝑥, 𝑦) = 𝐹 𝑢, 𝑣 𝑒 𝑀 𝑁
𝑢=0 𝑣=0
• It is separable and symmetric, with its one-dimensional version being
(almost) unitary.
• Therefore, the one-dimensional DFT and more specifically the FFT, can be
used for the implementation of the two-dimensional DFT ( 2D FFT.
 Two-dimensional Discrete Fourier Transform (2D DFT) cont.
1
• The location of the is not important.
𝑀𝑁
• Sometimes it is located in front of the inverse transform.
1
• Other times it is found split into two equal terms of multiplying the
𝑀𝑁
transform and its inverse.
• In the cases of square images 𝑀 = 𝑁.
• Very often in Image Processing we work with square images whose size is
a power of 2. Powers of 2 facilitate easier implementation of DSP
algorithms as you already know.
• The amplitude spectrum of 2D DFT is 𝐹(𝑢, 𝑣) = 𝑅2 𝑢, 𝑣 + 𝐼2 𝑢, 𝑣 1/2 ,
𝑅(𝑢, 𝑣) and 𝐼(𝑢, 𝑣) are the real and imaginary parts of 𝐹(𝑢, 𝑣), respectively.
𝐼(𝑢,𝑣)
• The phase spectrum is 𝜙 𝑢, 𝑣 = tan−1
𝑅(𝑢,𝑣)
• The power spectrum is 𝑃 𝑢, 𝑣 = 𝐹(𝑢, 𝑣) 2 = 𝑅2 𝑢, 𝑣 + 𝐼2 𝑢, 𝑣
 The visualisation of the range of values of 2D DFT
• The range of values of 𝐹(𝑢, 𝑣) is typically very large.
• Due to quantization of 𝐹(𝑢, 𝑣) small values are not distinguishable when
we attempt to display the amplitude of 𝐹 𝑢, 𝑣 .
• We, therefore, apply a logarithmic transformation to enhance small values.
𝐷 𝑢, 𝑣 = 𝑐log 1 + 𝐹(𝑢, 𝑣)
• The parameter 𝑐 is chosen so that
the range of 𝐷(𝑢, 𝑣) is [0, 255], i.e.,
255
𝑐= .
log 1+max{ 𝐹 𝑢,𝑣 }

Original image Display of amplitude Display of the logarithmic

of DFT amplitude of DFT
 The visualisation of the range of values of 2D DFT cont.

Display of amplitude of DFT Display of the logarithmic amplitude of DFT

The visualisation of the range of values of the 2D DFT cont.

Display of amplitude of DFT Display of the logarithmic amplitude of DFT

 Properties of 2D DFT
• Periodicity: The 2D DFT and its inverse are periodic
𝐹 𝑢, 𝑣 = 𝐹 𝑢 + 𝑀, 𝑣 = 𝐹 𝑢, 𝑣 + 𝑁 = 𝐹(𝑢 + 𝑀, 𝑣 + 𝑁)
Therefore, use of DFT implies virtual periodicity in space.
• Conjugate symmetry: 𝐹 𝑢, 𝑣 = 𝐹 ∗ −𝑢 + 𝑝𝑀, −𝑣 + 𝑞𝑁 with 𝑝, 𝑞 any
integers. This property also implies that 𝐹(𝑢, 𝑣) = 𝐹(−𝑢, −𝑣) .
• If 𝑓 𝑥, 𝑦 is real and even then 𝐹 𝑢, 𝑣 is real and even.
• If 𝑓 𝑥, 𝑦 is real and odd then 𝐹 𝑢, 𝑣 is imaginary and odd.
• ℱ 𝑓 𝑥, 𝑦 + 𝑔 𝑥, 𝑦 = ℱ 𝑓 𝑥, 𝑦 + ℱ 𝑔 𝑥, 𝑦 where ℱ{∙} indicates the
Discrete Fourier Transform operator.
• ℱ 𝑓 𝑥, 𝑦 ∙ 𝑔 𝑥, 𝑦 ≠ ℱ 𝑓 𝑥, 𝑦 ∙ ℱ 𝑔 𝑥, 𝑦
• Translation in spatial and frequency domain:
𝑢𝑥0 𝑣𝑦0
−𝑗2𝜋 +
𝑓 𝑥 − 𝑥0 , 𝑦 − 𝑦0 ↔ 𝐹(𝑢, 𝑣)𝑒 𝑀 𝑁
𝑢0 𝑥 𝑣0 𝑦
𝑗2𝜋 𝑀 + 𝑁
𝑓 𝑥, 𝑦𝑒 ↔ 𝐹(𝑢 − 𝑢0 , 𝑣 − 𝑣0 )
1 𝑀−1 𝑁−1
• Average value of the signal 𝑓 (𝑥, 𝑦) = 𝑥=0 𝑦=0 𝑓(𝑥, 𝑦)
𝑀𝑁
1 𝑀−1 𝑁−1
𝐹 0,0 = 𝑥=0 𝑦=0 𝑓(𝑥, 𝑦) ⇒ 𝑓 (𝑥, 𝑦) = 𝐹(0,0)
𝑀𝑁
 Properties of 2D DFT cont.

• Rotation: Rotating 𝑓(𝑥, 𝑦) by 𝜃 rotates 𝐹 𝑢, 𝑣 by 𝜃.

 Visualisation of a full period of 2D DFT

• In order to display a full period of the 2D DFT in the center of the image,
we need to translate the origin of the transform at 𝑢, 𝑣 = (𝑀/2, 𝑁/2).
• To move 𝐹 𝑢, 𝑣 at (𝑀/2, 𝑁/2) we use the second translation property in
the previous slide with 𝑢0 = 𝑀/2 and 𝑣0 = 𝑁/2.
𝑢 𝑥 𝑣 𝑦
0 + 0
𝑗2𝜋 𝑀
In that case ℱ 𝑓 𝑥, 𝑦 𝑒 𝑁 = 𝐹(𝑢 − 𝑢0 , 𝑣 − 𝑣0 ) becomes
ℱ 𝑓 𝑥, 𝑦 𝑒 𝑗𝜋 𝑥+𝑦 = ℱ 𝑓 𝑥, 𝑦 (−1) 𝑥+𝑦 = 𝐹(𝑢 − 𝑀/2, 𝑣 − 𝑁/2)

Original image Original DFT amplitude Translated DFT amplitude

 Applications:
Exploiting the property of energy compaction of 2D DFT
 Applications: Low- and high- pass filtering:
How do frequencies show up in an image?

• Low frequencies correspond to slowly varying pixel intensities (e.g.,

smooth surfaces).
• High frequencies correspond to quickly varying pixel intensities (e.g.,
edges)

original image low-pass filtered image

 Low- and high- pass filtering using amplitude of 2D DFT

original image low pass filtered high pass filtered

amplitude of DFT amplitude of DFT with amplitude of DFT with

eliminated high frequencies eliminated low frequencies
 Amplitude and phase of 2D DFT

original image amplitude of DFT phase of DFT

amplitude of DFT amplitude of DFT phase of DFT

Exercise: use IDFT to reconstruct the image using magnitude or phase only information
Amplitude and phase of 2D DFT cont.

Reconstructed image using the

amplitude of DFT only and zero phase.
The magnitude determines the
strength of each frequency component but
not its location.

Reconstructed image using the phase

of DFT only and amplitude equal to 1.
• The phase determines the locations of
individual frequencies within the image.
• High frequencies → abrupt changes →
edges → object silhouettes → image
content.
Reconstruction from the amplitude of one image
and the phase of another

phase of cameraman phase of grasshopper

amplitude of grasshopper amplitude of cameraman
Reconstruction from the amplitude of one image
and the phase of another cont.

phase of buffalo phase of rocks

amplitude of rocks amplitude of buffalo
What follows in the next slides are experiments done by students of the
MSc in Communications and Signal Processing
 Experiment run by Zheng Zhang
MSc Student in Communications and Signal Processing, 2018-2019.
 Experiment run by Zheng Zhang,
MSc Student in Communications and Signal Processing, 2018-2019.
 Experiment run by Zheng Zhang,
MSc Student in Communications and Signal Processing, 2018-2019.
 Experiment run by Ximo Chen
MSc Student in Communications and Signal Processing, 2018-2019.
 Experiment run by Ximo Chen
MSc Student in Communications and Signal Processing, 2018-2019.
 Experiment run by Yunping Zhang
MSc Student in Communications and Signal Processing, 2018-2019.

Original image Reconstruction using Reconstruction using

amplitude only phase only
 Experiment run by Yunping Zhang
MSc Student in Communications and Signal Processing, 2018-2019.

Original image “Tea” Original image “Drawing”

Phase of “Tea”, amplitude of “Drawing” Phase of “Drawing”, amplitude of “Tea”

 Experiment run by Ai BO
MSc Student in Communications and Signal Processing, 2018-2019.
 Experiment run by Yunyi Guang
MSc Student in Communications and Signal Processing, 2018-2019.

Original image Reconstruction using Reconstruction using

amplitude only phase only
 Experiment run by Yunyi Guang
MSc Student in Communications and Signal Processing, 2018-2019.

Original image “Girl” Original image “Sticker”

Phase of “Girl”, amplitude of “Sticker” Phase of “Sticker”, amplitude of “Girl”