Transcript Chapter 3: Digital Image Processing

IUST
Chapter 4:
Frequency Domain Processing

Image Transformations
1
Contents








Fourier Transform and DFT
Walsh Transform
Hadamard Transform
Walsh-Hadamard Transform (WHT)
Discrete Cosine Transform (DCT)
Haar Transform
Slant Transform
Comparison of various Transforms
2
Introduction

Although we discuss other transforms in
some detail in this chapter, we emphasize
the Fourier transform because of its wide
range of applications in image processing
problems.
3
Fourier Transform (1-D)
F u  

 f  x  exp  j2  ux  dx 

 F  u   AX
sin   uX
A
u
sin   uX  e
 j  uX

  uX 
4
Fourier Transform (2-D)
 
F(u , v) 

 f  x , y  exp  j2  ux  vy  dx dy


F  u , v   AXY
sin   uX
 uX 

sin   vY

 vY 
5
Discrete Fourier Transform
In the two-variable case the discrete Fourier transform pair is
6
Discrete Fourier Transform
When images are sampled in a squared array, i.e. M=N,
we can write
7
Discrete
Fourier
Transform
Examples
At all of these
examples, the Fourier
spectrum is shifted
from top left corner to
the center of the
frequency square.
8
Discrete
Fourier
Transform
Examples
At all of these
examples, the Fourier
spectrum is shifted
from top left corner to
the center of the
frequency square.
9
Discrete
Fourier
Transform
Display
At all of these
examples, the Fourier
spectrum is shifted
from top left corner to
the center of the
frequency square.
10
Discrete Fourier Transform
Example
Main Image (Gray Level)
DFT of Main image
Logarithmic scaled
(Fourier spectrum)
of the Fourier spectrum
11
Discrete Fourier Transform
MATLAB program page 1 from 3.

% Program written in Matlab for computing FFT of a given gray color image.

% Clear the memory.
clear;
% Getting the name and extension of the image file from the user.
name=input('Please write the name and address of the image : ','s');
% Reading the image file into variable 'a'.
a=imread(name);
% Computing the size of image. Assuming that image is squared.
N=length(a);
% Computing DFT of the image file by using fast Fourier algorithm.
F=fft2(double(a))/N;









12
Discrete Fourier Transform
MATLAB program page 2 from 3.













% Shifting the Fourier spectrum to the center of the frequency square.
for i=1:N/2; for j=1:N/2
G(i+N/2,j+N/2)=F(i,j);
end;end
for i=N/2+1:N; for j=1:N/2
G(i-N/2,j+N/2)=F(i,j);
end;end
for i=1:N/2; for j=N/2+1:N
G(i+N/2,j-N/2)=F(i,j);
end;end
for i=N/2+1:N; for j=N/2+1:N
G(i-N/2,j-N/2)=F(i,j);
end;end
13
Discrete Fourier Transform
MATLAB program page 3 from 3.











% Computing and scaling the logarithmic Fourier spectrum.
H=log(1+abs(G));
for i=1:N
H(i,:)=H(i,:)*255/abs(max(H(i,:)));
end
% Changing the color map to gray scale (8 bits).
colormap(gray(255));
% Showing the main image and its Fourier spectrum.
subplot(2,2,1),image(a),title('Main image');
subplot(2,2,2),image(abs(G)),title('Fourier spectrum');
subplot(2,2,3),image(H),title('Logarithmic scaled Fourier spectrum');
14
Discrete Fourier Transform
(Properties)

Separability
The discrete Fourier transform pair can be expressed in the seperable forms:
F u , v  
f x , y  
1
N
1
N 1
 exp 
x 0
N 1
j2  ux / N   f  x , y  exp  j2  vy / N 
exp  j2  ux

N
u 0
N 1
y0
N 1
/ N   F  u , v  exp  j2  vy / N 
v0
15
Discrete Fourier Transform
(Properties)

Translation
The translation properties of the
Fourier transform pair are :
f  x , y  exp  j2   u 0 x  v 0 y  / N   F  u  u 0 , v  v 0 
and
f  x  x 0 , y  y 0   F  u , v  exp  j2   ux 0  vy 0  / N 
16
Discrete Fourier Transform
(Properties)

Periodicity
The discrete Fourier transform and its
inverse are periodic with period N; that is,
F(u,v)=F(u+N,v)=F(u,v+N)=F(u+N,v+N)
If f(x,y) is real, the Fourier transform also
exhibits conjugate symmetry:
F(u,v)=F*(-u,-v)
Or, more interestingly:
|F(u,v)|=|F(-u,-v)|
17
Discrete Fourier Transform
(Properties)

Rotation
If we introduce the polar coordinates
x  r cos 
u   cos 
y  r sin 
v   sin 
Then we can write:
f  r ,    0   F  ,    0 
In other words, rotating F(u,v)
rotates f(x,y) by the same angle.
18
Discrete Fourier Transform
(Properties)

Convolution
The convolution theorem in
two dimensions is expressed
by the relations :
f  x , y  * g ( x , y )  F u , v  G u , v 
and
f  x , y  g  x , y   F u , v  * G u , v 
Note :
f x , y  * g x , y  
 
  f  ,   g  x   , y    d  d 

19
Discrete Fourier Transform
(Properties)

Correlation
The correlation of two continuous
functions f(x) and g(x) is defined
by the relation
f x   g x  



f
*
  g  x    d 
So we can write:
f  x , y   g  x , y   F u , v  G u , v 
*
and
f
*
 x , y  g  x , y   F u , v   G u , v 
20
Discrete
Fourier
Transform

Sampling
(Properties)
1-D
The Fourier transform and the convolution theorem provide the
tools for a deeper analytic study of sampling problem. In particular, we
want to look at the question of how many samples should be taken so that
no information is lost in the sampling process. Expressed differently, the
problem is one of the establishing the sampling conditions under which a
continuous image can be recovered fully from a set of sampled values. We
begin the analysis with the 1-D case.
As a result, a function which is band-limited in frequency domain
must extend from negative infinity to positive infinity in time domain (or x
domain).
21
Discrete
Fourier
Transform
(Properties)

Sampling
1-D
f(x) : a given function
F(u): Fourier Transform of f(x)
which is band-limited
s(x) : sampling function
Recovered f(x) from sampled data
S(u): Fourier Transform of s(x)
G(u): window for recovery of the
main function F(u) and f(x).
22
Discrete
Fourier
Transform
(Properties)

Sampling
1-D
f(x) : a given function
F(u): Fourier Transform of f(x)
which is band-limited
s(x) : sampling function
S(u): Fourier Transform of s(x)
h(x): window for making f(x)
finited in time.
H(u): Fourier Transform of h(x)
23
Discrete
Fourier
Transform
(Properties)

Sampling
1-D
s(x)*f(x) (convolution
function) is periodic, with
period 1/Δu. If N samples of
f(x) and F(u) are taken and the
spacings between samples are
selected so that a period in
each domain is covered by N
uniformly spaced samples,
then NΔx=X in the x domain
and NΔu=1/Δx in the
frequency domain.
24
Discrete
Fourier
Transform
(Properties)

Sampling
2-D
The sampling process for 2-D
functions can be formulated
mathematically by making use
of the 2-D impulse function
δ(x,y), which is defined as
 
  f x , y   x  x
0
, y  y 0  dx dy  f  x 0 , y 0 

A 2-D sampling function is
consisted of a train of impulses
separated Δx units in the x
direction and Δy units in the y
direction as shown in the figure.
25
Discrete
Fourier
Transform
(Properties)
 Sampling
2-D
If f(x,y) is band limited (that is, its
Fourier transform vanishes outside
some finite region R) the result of
covolving S(u,v) and F(u,v) might
look like the case in the case shown
in the figure. The function shown is
periodic in two dimensions.



G u , v   



1
(u,v) inside one of the rectangles
enclosing R
0
elsewhere
The inverse Fourier transform of
G(u,v)[S(u,v)*F(u,v)] yields f(x,y).
26
The Fast Fourier Transform (FFT) Algorithm
27
Other Seperable Image Transforms
For 2-D square arrays the forward and inverse transforms are
T u , v  
N 1 N 1
  f x , y  g x , y , u , v 
g(x,y,u,v) : forward transformation kernel
x 0 y0
and
f x , y  
N 1 N 1
  T u , v  h  x , y , u , v 
h(x,y,u,v) : inverse transformation kernel
u 0 v0
The forward kernel is said to be seperable if
g(x,y,u,v)=g1(x,u)g2(y,v)
In addition, the kernel is symmetric if g1 is functionally equal to g2. In this case
we can write:
g(x,y,u,v)=g1(x,u)g1(y,v)
28
Other Seperable Image Transforms
If the kernel g(x,y,u,v) is seperable and symmetric,
T u , v  
N 1 N 1
  f x , y  g x , y , u , v 
x 0 y0
also may be expressed in matrix form:
T  AFA
Where F is the N×N image matrix,
A is an N×N symmetric transformation matrix
T is the resulting N×N transform.
And for inverse transform we have:
BTB  BAFAB
B  A
1
 F  BTB
29
Walsh Transform
When N=2n, the 2-D forward and inverse Walsh kernels are given by the relations
g x , y , u , v  
1
N
n 1
 b i  x  b n 1 i  u  b i  y  b n 1 i  v  
  1
and
h x , y , u , v  
i0
1
N
n 1
b  x b



1

i
n 1 i
 u  b i  y  b n 1 i  v  
i0
Where bk(z) is the kth bit in the binary representation of z.
So the forward and inverse Walsh transforms are equal in form; that is:
30
Walsh Transform
In 1-D case we have :
g x , u  
1
N
n 1
b



1

i
 x  b n 1 i  u 
1-D kernel
i0
In the following table N=8 so n=3 (23=8).
“+” denotes for +1 and “-” denotes for -1.
31
Walsh Transform
This figure shows the basis functions (kernels) as
a function of u and v (excluding the 1/N constant
term) for computing the Walsh transform when
N=4. Each block corresponds to varying x and y
form 0 to 3 (that is, 0 to N-1), while keeping u
and v fixed at the values corresponding to that
block. Thus each block consists of an array of
4×4 binary elements (White means “+1” and
Black means “-1”). To use these basis functions
to compute the Walsh transform of an image of
size 4×4 simply requires obtaining W(0,0) by
multiplying the image array point-by-point with
the 4×4 basis block corresponding to u=0 and
v=0, adding the results, and dividing by 4, and
continue for other values of u and v.
32
Walsh Transform
Example
Main Image (Gray Level)
WT of Main image
Logarithmic scaled
(Walsh spectrum)
of the Walsh spectrum
33
Walsh Transform (WT)
MATLAB program page 1 from 3.













% Program written in Matlab for computing WT of a given gray color image.
clear;
% Getting the name and extension of the image file from the user.
name=input('Please write the name and address of the image : ','s');
a=imread(name);
N=length(a);
% Computing Walsh Transform of the image file.
n=log2(N);n=1+fix(n);f=ones(N,N);
for x=1:N; for u=1:N
p=dec2bin(x-1,n); q=dec2bin(u-1,n);
for i=1:n; f(x,u)=f(x,u)*((-1)^(p(n+1-i)*q(i)));
end;end;end
F=(1/N)*f*double(a)*f;
34
Walsh Transform (WT)
MATLAB program page 2 from 3.













% Shifting the Fourier spectrum to the center of the frequency square.
for i=1:N/2; for j=1:N/2
G(i+N/2,j+N/2)=F(i,j);
end;end
for i=N/2+1:N; for j=1:N/2
G(i-N/2,j+N/2)=F(i,j);
end;end
for i=1:N/2; for j=N/2+1:N
G(i+N/2,j-N/2)=F(i,j);
end;end
for i=N/2+1:N; for j=N/2+1:N
G(i-N/2,j-N/2)=F(i,j);
end;end
35
Walsh Transform (WT)
MATLAB program page 3 from 3.











% Computing and scaling the logarithmic Walsh spectrum.
H=log(1+abs(G));
for i=1:N
H(i,:)=H(i,:)*255/abs(max(H(i,:)));
end
% Changing the color map to gray scale (8 bits).
colormap(gray(255));
% Showing the main image and its Walsh spectrum.
subplot(2,2,1),image(a),title('Main image');
subplot(2,2,2),image(abs(G)),title('Walsh spectrum');
subplot(2,2,3),image(H),title('Logarithmic scaled Walsh spectrum');
36
Hadamard Transform
When N=2n, the 2-D forward and inverse Hadamard kernels are given by the relations
g x , y , u , v  
1
N
n 1
 b i  x  b i  u  b i  y  b i  v  
  1
and
h x , y , u , v  
i0
1
N
n 1
 b  x  b  u  b  y  b  v  



1

i
i
i
i
i0
Where bk(z) is the kth bit in the binary representation of z.
So the forward and inverse Hadamard transforms are equal in form; that is:
37
Hadamard Transform
In 1-D case we have :
g x , u  
1
N
n 1
b



1

i
 x b i  u 
1-D kernel
i0
In the following table N=8 so n=3 (23=8).
“+” denotes for +1 and “-” denotes for -1.
38
Hadamard Transform
This figure shows the basis functions
(kernels) as a function of u and v (excluding
the 1/N constant term) for computing the
Hadamard transform when N=4. Each block
corresponds to varying x and y form 0 to 3
(that is, 0 to N-1), while keeping u and v
fixed at the values corresponding to that
block. Thus each block consists of an array
of 4×4 binary elements (White means “+1”
and Black means “-1”) like Walsh transform.
If we compare these two transforms we can
see that they only differ in the sense that the
functions in Hadamard transform are ordered
in increasing sequency and thus are more
“natural” to interpret.
39
Hadamard Transform
Example
Main Image (Gray Level)
HT of Main image
Logarithmic scaled
(Hadamard spectrum)
of the Hadamard spectrum
40
Hadamard Transform (HT)
MATLAB program page 1 from 3.













% Program written in Matlab for computing HT of a given gray color image.
clear;
% Getting the name and extension of the image file from the user.
name=input('Please write the name and address of the image : ','s');
a=imread(name);
N=length(a);
% Computing Hadamard Transform of the image file.
n=log2(N);n=1+fix(n);f=ones(N,N);
for x=1:N; for u=1:N
p=dec2bin(x-1,n); q=dec2bin(u-1,n);
for i=1:n; f(x,u)=f(x,u)*((-1)^(p(n+1-i)*q(n+1-i)));
end;end;end
F=(1/N)*f*double(a)*f;
41
Hadamard Transform (HT)
MATLAB program page 2 from 3.













% Shifting the Fourier spectrum to the center of the frequency square.
for i=1:N/2; for j=1:N/2
G(i+N/2,j+N/2)=F(i,j);
end;end
for i=N/2+1:N; for j=1:N/2
G(i-N/2,j+N/2)=F(i,j);
end;end
for i=1:N/2; for j=N/2+1:N
G(i+N/2,j-N/2)=F(i,j);
end;end
for i=N/2+1:N; for j=N/2+1:N
G(i-N/2,j-N/2)=F(i,j);
end;end
42
Hadamard Transform (HT)
MATLAB program page 3 from 3.











% Computing and scaling the logarithmic Hadamard spectrum.
H=log(1+abs(G));
for i=1:N
H(i,:)=H(i,:)*255/abs(max(H(i,:)));
end
% Changing the color map to gray scale (8 bits).
colormap(gray(255));
% Showing the main image and its Hadamard spectrum.
subplot(2,2,1),image(a),title('Main image');
subplot(2,2,2),image(abs(G)),title('Hadamard spectrum');
subplot(2,2,3),image(H),title('Logarithmic scaled Hadamard
spectrum');
43
Walsh-Hadamard Transform (WHT)
44
1-D WHT Kernel Functions
45
2-D WHT Kernel Functions
46
WHT and Fourier Transform
47
WHT Example
48
Discrete Cosine Transform
49
Discrete Cosine Transform (DCT)
Each block consists
of 4×4 elements,
corresponding to x
and y varying from 0
to 3. The highest
value is shown in
white. Other values
are shown in grays,
with darker meaning
smaller.
50
Discrete Cosine Transform
Example
Main Image (Gray Level)
DCT of Main image
Logarithmic scaled
(Cosine spectrum)
of the Cosine spectrum
51
Discrete Cosine Transform
MATLAB program page 1 from 3.

% Program written in Matlab for computing DCT of a given gray color image.

% Clear the memory.
clear;
% Getting the name and extension of the image file from the user.
name=input('Please write the name and address of the image : ','s');
% Reading the image file into variable 'a'.
a=imread(name);
% Computing the size of image. Assuming that image is squared.
N=length(a);
% Computing DCT of the image file.
F=dct2(double(a));









52
Discrete Cosine Transform
MATLAB program page 2 from 3.













% Shifting the Fourier spectrum to the center of the frequency square.
for i=1:N/2; for j=1:N/2
G(i+N/2,j+N/2)=F(i,j);
end;end
for i=N/2+1:N; for j=1:N/2
G(i-N/2,j+N/2)=F(i,j);
end;end
for i=1:N/2; for j=N/2+1:N
G(i+N/2,j-N/2)=F(i,j);
end;end
for i=N/2+1:N; for j=N/2+1:N
G(i-N/2,j-N/2)=F(i,j);
end;end
53
Discrete Cosine Transform
MATLAB program page 3 from 3.











% Computing and scaling the logarithmic Cosine spectrum.
H=log(1+abs(G));
for i=1:N
H(i,:)=H(i,:)*255/abs(max(H(i,:)));
end
% Changing the color map to gray scale (8 bits).
colormap(gray(255));
% Showing the main image and its Cosine spectrum.
subplot(2,2,1),image(a),title('Main image');
subplot(2,2,2),image(abs(G)),title('Cosine spectrum');
subplot(2,2,3),image(H),title('Logarithmic scaled Cosine spectrum');
54
DCT and Fourier Transform
55
DCT Example
56
Blockwise DCT Example
57
Fast DCT Algorithm
58
Fast DCT Algorithm
59
Haar Transform
The Haar transform is based on the Haar functions, hk(z), which are defined over the
continuous, closed interval [0,1] for z, and for k=0,1,2,…,N-1, where N=2n. The first
step in generating the Haar transform is to note that the integer k can be decomposed
uniquely as
k=2p+q-1
where 0≤p≤n-1, q=0 or 1 for p=0, and 1≤q≤2p for p≠0.
With this background, the Haar functions are defined as

h 0  z   h 00  z  
1
for z  0 ,1
N
and

h k  z   h 00  z  
1
N
q 1
q 1/ 2
 p/2
2

z

p
p

2
2

q 1/ 2
q

p/2

z

 2
p
p
2
2

otherwise for z  0 ,1
0


60
Haar Transform

Haar transform matrix for sizes N=2,4,8
Hr 2 
1 1

2 1
1 

 1
1

1 1
Hr 4 
4 1

1
1
1
2

2
1
0
1
0
0 

0 
2 

 2 
1

1
1

1 1
Hr 8 
8 1

1

1

1
1
2
0
2
0
0
1
2
0
2
0
0
1

2
0
0
2
0
1

2
0
0
2
0
1
0
2
0
0
2
1
0
2
0
0
2
1
0

2
0
0
0
1
0

2
0
0
0
0 

0 
0 

0 
0 

0 

2

 2 
Can be computed by taking sums and differences.
Fast algorithms by recursively applying Hr2.
61
Haar Transform Example
62
Haar Transform Example
63
Haar Transform
MATLAB program page 1 from 1.
% Program written in Matlab for computing HT of a given gray color image.
% Getting the name and extension of the image file from the user.
name=input('Please write the name and address of the image : ','s');
a=imread(name); N=length(a);
% Computing Haar Transform of the image file.
for i=1:N p=fix(log2(i)); q=i-(2^p);
for j=1:N
z=(j-1)/N;
if (z>=(q-1)/(2^p))&&(z<(q-1/2)/2^p) f(i,j)=(1/(sqrt(N)))*(2^(p/2));
elseif (z>=(q-1/2)/(2^p))&&(z<(q/2^p)) f(i,j)=(1/(sqrt(N)))*(-(2^(p/2)));
else f(i,j)=0;
end;end;end
F=f*double(a)*f;
% Changing the color map to gray scale (8 bits).
colormap(gray(255));
% Showing the main image and its Hadamard spectrum.
subplot(2,2,1),image(a),title('Main image');
subplot(2,2,2),image(abs(F)),title('Haar spectrum');
64
Slant Transform

The Slant Transform matrix of order N*N is the recursive expression
65
Slant Transform
Where I is the identity matrix, and
1 1
S2 

2 1
1 

 1
aN
 3N

 

2
4
N

1


bN
 N 4 
 

2
4
N

1


2

2

2

1
1
2

66
Slant Transform
Example
Main Image (Gray Level)
Slant Transform of Main image
(Slant spectrum)
67
Slant Transform
MATLAB program page 1 from 1.
% Program written in Matlab for computing ST of a given gray color image.
% Getting the name and extension of the image file from the user.
name=input('Please write the name and address of the image : ','s');
a=imread(name); N=length(a);
% Computing Slant Transform of the image file.
A=[ 1/(2^0.5) 1/(2^0.5);1/(2^0.5) -1/(2^0.5)];
for i=2:log2(N) N=2^i; aN=((3*(N^2))/(4*((N^2)-1)))^0.5;
bN=(((N^2)-4)/(4*((N^2)-1)))^0.5;
m=1/(2^0.5)*[1 0 zeros(1,(N/2)-2) 1 0 zeros(1,(N/2)-2)
aN bN zeros(1,(N/2)-2) -aN bN zeros(1,(N/2)-2)
zeros((N/2)-2,2) eye((N/2)-2) zeros((N/2)-2,2) eye((N/2)-2)
0 1
zeros(1,(N/2)-2) 0 -1 zeros(1,(N/2)-2)
-bN aN zeros(1,(N/2)-2) bN aN zeros(1,(N/2)-2)
zeros((N/2)-2,2) eye((N/2)-2) zeros((N/2)-2,2) -eye((N/2)-2)];
n=[A A-A;A-A A]; A=m*n;
end
F=A*double(a)*A;
% Changing the color map to gray scale (8 bits).
colormap(gray(255));
% Showing the main image and its Hadamard spectrum.
subplot(2,2,1),image(a),title('Main image');
subplot(2,2,2),image(abs(F)),title('Slant spectrum');
68
Comparison Of Various Transforms
69
Comparison Of Various Transforms
70
Comparison Of Various Transforms
71
Comparison Of Various Transforms
72
References




“Digital image processing” by Rafael C.
Gonzalez and Richard E. Woods
“Digital image processing” by Jain
Image communication I by Bernd Girod
Lecture 3, DCS339/AMCM053 by Pengwei
Hao, University of London
73