Image in Frequency Domain


By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

Source

  • An Introduction to Digital Image Processing with MATLAB
  • by Alasdair McAndrew at Victoria University of Technology

Table of Contents

1. Fourier Transforms of Images

1.1. Synthetic Images

In [1]:
% generate an image

N = 2^8;
m = [zeros(N,N/2) ones(N,N/2)];
imshow(m)

% just for plotting a boundary
hold on
w = [1 1];
x = [1 N];
y = [N 1];
z = [N N];
boundary = [w;x;z;y;w];  

plot(boundary(:,1),boundary(:,2),'k-');   % draw bounding box
hold off










    


























In [2]:



    


% Fourier Transformation with Image

M = fft2(m);
Ms = fftshift(M);
imshow(mat2gray(log(1+abs(Ms))))



    















    















    


























In [3]:



    


%plot -s 600,400

L = abs(Ms);
plot(1:length(L),L(N/2+1,:)), xlim([N/2-50,N/2+50])



    















    















    


























In [4]:



    


%plot -s 800,400

N = 2^8;
[x,y] = meshgrid(1:N,1:N);

k = 10;
m = cos(2*pi/N*k*x);

M = fft2(m);
Ms = fftshift(M);

subplot(1,2,1), imshow(m,[])
subplot(1,2,2), imshow(mat2gray(log(1+abs(Ms))),[])



    















    















    


























In [5]:



    


%plot -s 800,400

N = 2^8;
[x,y] = meshgrid(1:N,1:N);

k = 5;
m = cos(2*pi/N*k*y);

M = fft2(m);
Ms = fftshift(M);

subplot(1,2,1), imshow(m,[])
subplot(1,2,2), imshow(mat2gray(log(1+abs(Ms))),[])



    















    















    


























In [6]:



    


%plot -s 800,400

N = 2^8;
[x,y] = meshgrid(1:N,1:N);

k = 7;
m = cos(2*pi/N*k*(x-y));

M = fft2(m);
Ms = fftshift(M);

subplot(1,2,1), imshow(m,[])
subplot(1,2,2), imshow(mat2gray(log(1+abs(Ms))),[])



    















    















    


























In [7]:



    


%plot -s 900,300

% square (box) image and its fft

N = 2^8;
m = zeros(N,N);
m(N/2-15:N/2+15,N/2-15:N/2+15) = 1;

Ms = fftshift(fft2(m));

% inverse fft (recovered image)
Mr = ifftshift(Ms);
mr = real(ifft2(Mr));

subplot(1,3,1), imshow(m,[])
subplot(1,3,2), imshow(mat2gray(log(1+abs(Ms))),[])
subplot(1,3,3), imshow(mr,[])



    















    















    


























In [8]:



    


%plot -s 900,300

% a box rotated 45 degree
[x,y] = meshgrid(1:N,1:N);
m = (x+y < 329) & (x+y > 182) & (x-y > -67) & (x-y < 73);
Ms = fftshift(fft2(m));

% inverse fft
Mr = ifftshift(Ms);
mr = ifft2(Mr);

subplot(1,3,1), imshow(m,[])
subplot(1,3,2), imshow(mat2gray(log(1+abs(Ms))),[])
subplot(1,3,3), imshow(real(mr),[])



    















    















    


























In [9]:



    


%plot -s 900,300

% circle
[x,y] = meshgrid(-N/2:N/2-1,-N/2:N/2-1);
z = sqrt(x.^2+y.^2);

% OR
% [x,y] = meshgrid(1:N,1:N);
% z = sqrt((x-N/2).^2+(y-N/2).^2);

m = (z < 15);
Ms = fftshift(fft2(m));

% inverse fft
Mr = ifftshift(Ms);
mr = ifft2(Mr);

subplot(1,3,1), imshow(m,[]), title('original')
subplot(1,3,2), imshow(mat2gray(log(1+abs(Ms)))), title('FFT')
subplot(1,3,3), imshow(real(mr),[]), title('recovered')



    















    















    




























Note "ringing" in the Fourier transform. This artifact is associated with the sharp cutoff of the circle. Think of a sinc function.















1.2. Circle with a Gentle Cutoff














Butterworh




$$ \frac{1}{1+\alpha(x^2+y^2)} $$











In [10]:



    


%plot -s 900,300

[x,y] = meshgrid(-N/2:N/2-1,-N/2:N/2-1);
z = sqrt(x.^2+y.^2);
d = 1./(1+(z./15).^2);
Ds = fftshift(fft2(d));

% inverse fft
Dr = ifftshift(Ds);
dr = ifft2(Dr);

subplot(1,3,1), imshow(d,[]), title('original')
subplot(1,3,2), imshow(mat2gray(log(1+abs(Ds)))), title('FFT')
subplot(1,3,3), imshow(real(dr),[]), title('recovered')



    















    















    






























$$e^{\left(-\alpha_1 x^2 - \alpha_2 y^2\right)}$$














In [11]:



    


%plot -s 900,300

[x,y] = meshgrid(-N/2:N/2-1,-N/2:N/2-1);
z = exp(-0.01*x.^2 - 0.05*y.^2);        % Gaussian in time

Cs = fftshift(fft2(z));                 % Also Gaussian in freq.

subplot(1,2,1), imshow(z,[]), title('original')
subplot(1,2,2), imshow(mat2gray(log(1+abs(Cs)))), title('FFT')
%subplot(1,2,2), imshow(abs(Cs))



    















    















    




























1.3. Real Images













In [12]:



    


%plot -s 900,300

% fft with an image
cm = imread('./image_files/cameraman.gif');
imshow(cm)

CMs = fftshift(fft2(cm));

CMr = ifftshift(CMs);
cmr = ifft2(CMr);

subplot(1,3,1), imshow(cm,[]), title('original')
subplot(1,3,2), imshow(mat2gray(log(1+abs(CMs)))), title('FFT')
%subplot(1,3,3), imshow(real(cmr),[])
%subplot(1,3,3), imshow(abs(cmr),[])
subplot(1,3,3), imshow(uint8(cmr),[]), title('recovered')



    















    















    




























2. Image Filtering in Frequency














2.1. Low-Pass Filtering in Frequency














    

  • ideal filter
    • 














In [13]:



    


%plot -s 800,800

[x,y] = meshgrid(-N/2:N/2-1,-N/2:N/2-1);
z = sqrt(x.^2+y.^2);

C = (z < 15);
CMFs = CMs.*C;

CMFr = ifftshift(CMFs);
cmfr = ifft2(CMFr);

subplot(2,2,1), imshow(cm,[]), title('original image')
subplot(2,2,2), imshow(mat2gray(log(1+abs(C))),[]), title('filter in freq.')
subplot(2,2,3), imshow(mat2gray(log(1+abs(CMFs))),[]), title('filtered image in freq.')
subplot(2,2,4), imshow(real(cmfr),[]), title('filtered image')



    















    















    




























    

  • Butterworth $ \frac{1}{1+\alpha(x^2+y^2)} $ type
    • 














In [14]:



    


%plot -s 800,800
 
N = 2^8;
[x,y] = meshgrid(1:N,1:N);
z = sqrt((x-N/2).^2+(y-N/2).^2);
D = 1./(1+(z./10).^2);               % filter in the frequency domain

CMFs = CMs.*D;

CMFr = ifftshift(CMFs);
cmfr = ifft2(CMFr);

subplot(2,2,1), imshow(cm), title('original image')
subplot(2,2,2), imshow(mat2gray(log(1+abs(D))),[]), title('filter in freq.')
subplot(2,2,3), imshow(mat2gray(log(1+abs(CMFs)))), title('filtered image in freq.')
subplot(2,2,4), imshow(real(cmfr),[]), title('filtered image')   
% the second argument is empty square brackets. 
% This maps the minimum value in the image to black and the maximum value in the image to white.



    















    















    




























    

  • Guassian type
    • 














In [15]:



    


%plot -s 800,800

[x,y] = meshgrid(-N/2:N/2-1,-N/2:N/2-1);
Z = exp(-0.0005*x.^2 - 0.0005*y.^2);        % Gaussian 

CMFs = CMs.*Z;

CMFr = ifftshift(CMFs);
cmfr = ifft2(CMFr);

subplot(2,2,1), imshow(cm), title('original image')
subplot(2,2,2), imshow(mat2gray(log(1+abs(Z))),[]), title('filter in freq.')
subplot(2,2,3), imshow(mat2gray(log(1+abs(CMFs)))), title('filtered image in freq.')
subplot(2,2,4), imshow(real(cmfr),[]), title('filtered image')



    















    















    




























Finding Abraham Lincoln















In [16]:



    


%plot -s 900,500

im = imread('./image_files/gala.jpg'); % rgb (color) image 
imgray = rgb2gray(im);

subplot(1,2,1), imshow(im), title('colored original painting')
subplot(1,2,2), imshow(imgray), title('grayed painting')



    















    















    


























In [17]:



    


%plot -s 900,500

IM = fft2(imgray);
IMs = fftshift(IM);

% Gaussian 
[N,M] = size(IMs);
[x,y] = meshgrid(-M/2:M/2-1,-N/2:N/2-1);
Z = exp(-0.005*x.^2 - 0.005*y.^2);        

IMFs = IMs.*Z;

IMFr = ifftshift(IMFs);
imfr = ifft2(IMFr);

subplot(1,2,1), imshow(imgray,[]), title('Grayed Painting')
subplot(1,2,2), imshow(real(imfr),[]), title('Blurred Image')



    















    















    




























2.2. High-Pass Filter in Freq.
    

  • sharpen (or shows the edges of) an image
    • 














In [18]:



    


%plot -s 800,800

N = 2^8;
[x,y] = meshgrid(-N/2:N/2-1,-N/2:N/2-1);
z = sqrt(x.^2+y.^2);

C = (z > 40);
CMFs = CMs.*C;

CMFr = ifftshift(CMFs);
cmfr = ifft2(CMFr);

subplot(2,2,1), imshow(cm,[]), title('original image')
subplot(2,2,2), imshow(mat2gray(log(1+abs(C))),[]), title('filter in freq.')
subplot(2,2,3), imshow(mat2gray(log(1+abs(CMFs))),[]), title('filtered image in freq.')
subplot(2,2,4), imshow(real(cmfr),[]), title('filtered image')



    















    















    






























$$1-e^{\left(-\alpha_1 x^2 - \alpha_2 y^2\right)}$$














In [19]:



    


%plot -s 800,800

[x,y] = meshgrid(-N/2:N/2-1,-N/2:N/2-1);
Z = 1 - exp(-0.0005*x.^2 - 0.0005*y.^2);        % Gaussian 

CMFs = CMs.*Z;

CMFr = ifftshift(CMFs);
cmfr = ifft2(CMFr);

subplot(2,2,1), imshow(cm,[]), title('original image')
subplot(2,2,2), imshow(mat2gray(log(1+abs(Z))),[]), title('filter in freq.')
subplot(2,2,3), imshow(mat2gray(log(1+abs(CMFs))),[]), title('filtered image in freq.')
subplot(2,2,4), imshow(real(cmfr),[]), title('filtered image')



    















    















    






























$$ 1- \frac{1}{1+\alpha(x^2+y^2)} $$














In [20]:



    


%plot -s 800,800
 
N = 2^8;
[x,y] = meshgrid(1:N,1:N);
z = sqrt((x-N/2).^2+(y-N/2).^2);
D = 1 - 1./(1+(z./10).^2);               % filter in the frequency domain

CMFs = CMs.*D;

CMFr = ifftshift(CMFs);
cmfr = ifft2(CMFr);

subplot(2,2,1), imshow(cm,[]), title('original image')
subplot(2,2,2), imshow(mat2gray(log(1+abs(D))),[]), title('filter in freq.')
subplot(2,2,3), imshow(mat2gray(log(1+abs(CMFs))),[]), title('filtered image in freq.')
subplot(2,2,4), imshow(real(cmfr),[]), title('filtered image')  
% the second argument is empty square brackets. 
% This maps the minimum value in the image to black and the maximum value in the image to white.



    















    















    




























2.3. Edge Detection













In [21]:



    


%plot -s 600,500

cm = imread('./image_files/cameraman.gif');
cme = edge(cm,'canny');
imshow(cme)



    















    















    


























In [22]:



    


%plot -s 600,500

cmd = bwmorph(cme,'dilate',1);
imshow(cmd)



    















    















    




























2.4. Notch Filter
    

  • band-stop filter
    • 

  • periodic noise reduction
    • 















In [23]:



    


%plot -s 400,300

% load a test image (with patterned noises)

im = imread('./image_files/ClownOrig.jpg');
im = rgb2gray(im);
imshow(im,[])



    















    















    


























In [24]:



    


%plot -s 400,300

% transform to Fourier Domain

IM = fft2(im);
IMs = fftshift(IM);

imshow(mat2gray(log(1+abs(IMs))))
% recognize unwnated spots that are related to periodic noises



    















    















    


























In [25]:



    


%plot -s 800,800

%run it in Matlab for interactive demo

N = length(IMs);
[x,y] = meshgrid(1:N,1:N);

% notch filter generation (need to understand)
a1 = 0.008;
a2 = 0.008;
NF1 = 1 - exp(-a1*(x-190).^2 - a2*(y-123).^2);        % Gaussian 
NF2 = 1 - exp(-a1*(x-104).^2 - a2*(y-172).^2);        % Gaussian 
NF3 = 1 - exp(-a1*(x-126).^2 - a2*(y-135).^2);        % Gaussian 
NF4 = 1 - exp(-a1*(x-168).^2 - a2*(y-161).^2);        % Gaussian 

Z = NF1.*NF2.*NF3.*NF4;

IMFs = IMs.*Z;

IMFr = ifftshift(IMFs);
imfr = ifft2(IMFr);

subplot(2,2,1), imshow(im,[])
subplot(2,2,2), imshow(mat2gray(log(1+abs(Z))),[])
subplot(2,2,3), imshow(mat2gray(log(1+abs(IMFs))),[])
subplot(2,2,4), imshow(real(imfr),[])



    















    















    








































In [26]:



    


%plot -s 800,400

% load a test image (with patterned noises)
im = imread('./image_files/lena2.jpg');
im = rgb2gray(im);
IM = fft2(im);
IMs = fftshift(IM);

subplot(1,2,1), imshow(im,[])
subplot(1,2,2), imshow(mat2gray(log(1+abs(IMs))))



    















    















    








































In [27]:



    


%plot -s 800,400

% load a test image (with patterned noises)
im = imread('./image_files/script.jpg');
im = rgb2gray(im);
IM = fft2(im);
IMs = fftshift(IM);

subplot(1,2,1), imshow(im,[])
subplot(1,2,2), imshow(mat2gray(log(1+abs(IMs))))



    















    















    








































In [28]:



    


%plot -s 600,600

im = imread('./image_files/halftone.png');
imshow(im,[])



    















    















    


























In [29]:



    


IM = fft2(im);
IMs = fftshift(IM);
imshow(mat2gray(log(1+abs(IMs))))



    















    















    


























In [30]:



    


%plot -s 800,800

N = length(IMs);
[x,y] = meshgrid(1:N,1:N);

% notch filter generation (need to understand)
a1 = 0.005;
a2 = 0.005;
NF1 = 1 - exp(-a1*(x-400).^2 - a2*(y-134).^2);        % Gaussian 
NF2 = 1 - exp(-a1*(x-400).^2 - a2*(y-267).^2);        % Gaussian 
NF3 = 1 - exp(-a1*(x-400).^2 - a2*(y-400).^2);        % Gaussian 
NF4 = 1 - exp(-a1*(x-400).^2 - a2*(y-533).^2);        % Gaussian 
NF5 = 1 - exp(-a1*(x-266).^2 - a2*(y-134).^2);        % Gaussian 
NF6 = 1 - exp(-a1*(x-266).^2 - a2*(y-267).^2);        % Gaussian 
NF7 = 1 - exp(-a1*(x-266).^2 - a2*(y-400).^2);        % Gaussian 
NF8 = 1 - exp(-a1*(x-266).^2 - a2*(y-533).^2);        % Gaussian 

NF9  = 1 - exp(-a1*(x-533).^2 - a2*(y-134).^2);        % Gaussian 
NF10 = 1 - exp(-a1*(x-533).^2 - a2*(y-267).^2);        % Gaussian 
NF11 = 1 - exp(-a1*(x-533).^2 - a2*(y-400).^2);        % Gaussian 
NF12 = 1 - exp(-a1*(x-533).^2 - a2*(y-533).^2);        % Gaussian 
NF13 = 1 - exp(-a1*(x-133).^2 - a2*(y-134).^2);        % Gaussian 
NF14 = 1 - exp(-a1*(x-133).^2 - a2*(y-267).^2);        % Gaussian 
NF15 = 1 - exp(-a1*(x-133).^2 - a2*(y-400).^2);        % Gaussian 
NF16 = 1 - exp(-a1*(x-133).^2 - a2*(y-533).^2);        % Gaussian 

NF17 = 1 - exp(-a1*(x-66).^2 - a2*(y-65).^2);        % Gaussian 
NF18 = 1 - exp(-a1*(x-66).^2 - a2*(y-200).^2);        % Gaussian 
NF19 = 1 - exp(-a1*(x-66).^2 - a2*(y-334).^2);        % Gaussian 
NF20 = 1 - exp(-a1*(x-66).^2 - a2*(y-468).^2);        % Gaussian 
NF21 = 1 - exp(-a1*(x-66).^2 - a2*(y-602).^2);        % Gaussian 

NF22 = 1 - exp(-a1*(x-200).^2 - a2*(y-65).^2);        % Gaussian 
NF23 = 1 - exp(-a1*(x-200).^2 - a2*(y-200).^2);        % Gaussian 
NF24 = 1 - exp(-a1*(x-200).^2 - a2*(y-334).^2);        % Gaussian 
NF25 = 1 - exp(-a1*(x-200).^2 - a2*(y-468).^2);        % Gaussian 
NF26 = 1 - exp(-a1*(x-200).^2 - a2*(y-602).^2);  

NF27 = 1 - exp(-a1*(x-334).^2 - a2*(y-65).^2);        % Gaussian 
NF28 = 1 - exp(-a1*(x-334).^2 - a2*(y-200).^2);        % Gaussian 
NF29 = 1 - exp(-a1*(x-334).^2 - a2*(y-334).^2);        % Gaussian 
NF30 = 1 - exp(-a1*(x-334).^2 - a2*(y-468).^2);        % Gaussian 
NF31 = 1 - exp(-a1*(x-334).^2 - a2*(y-602).^2);  

NF32 = 1 - exp(-a1*(x-466).^2 - a2*(y-65).^2);        % Gaussian 
NF33 = 1 - exp(-a1*(x-466).^2 - a2*(y-200).^2);        % Gaussian 
NF34 = 1 - exp(-a1*(x-466).^2 - a2*(y-334).^2);        % Gaussian 
NF35 = 1 - exp(-a1*(x-466).^2 - a2*(y-468).^2);        % Gaussian 
NF36 = 1 - exp(-a1*(x-466).^2 - a2*(y-602).^2);  

NF37 = 1 - exp(-a1*(x-600).^2 - a2*(y-65).^2);        % Gaussian 
NF38 = 1 - exp(-a1*(x-600).^2 - a2*(y-200).^2);        % Gaussian 
NF39 = 1 - exp(-a1*(x-600).^2 - a2*(y-334).^2);        % Gaussian 
NF40 = 1 - exp(-a1*(x-600).^2 - a2*(y-468).^2);        % Gaussian 
NF41 = 1 - exp(-a1*(x-600).^2 - a2*(y-602).^2);  

%Z = NF1.*NF2.*NF3.*NF4.*NF5.*NF6.*NF7.*NF8;
% Z = Z.*NF9.*NF10.*NF11.*NF12.*NF13.*NF14.*NF15.*NF16;

% Z = Z.*NF17.*NF18.*NF19.*NF20.*NF21.*NF22.*NF23.*NF24.*NF25.*NF26;
% Z = Z.*NF27.*NF28.*NF29.*NF30.*NF31.*NF32.*NF33.*NF34.*NF35.*NF36;
% Z = Z.*NF37.*NF38.*NF39.*NF40.*NF41;

Z = NF2.*NF3.*NF6.*NF7;
IMFs = IMs.*Z;

IMFr = ifftshift(IMFs);
imfr = ifft2(IMFr);

subplot(2,2,1), imshow(im,[])
subplot(2,2,2), imshow(mat2gray(log(1+abs(Z))),[])
subplot(2,2,3), imshow(real(imfr),[]) 
subplot(2,2,4), imshow(mat2gray(log(1+abs(IMFs))),[])



    















    















    


























In [31]:



    


%plot -s 600,600

imshow(real(imfr),[])



    















    















    


























In [32]:



    


%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')