1. Fourier Transforms of Images¶

1.1. Synthetic Images¶

% 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

% Fourier Transformation with Image

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

%plot -s 600,400

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

%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))),[])

%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))),[])

%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))),[])

%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,[])

%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),[])

%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)} $$

%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)}$$

%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¶

Download the image

%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

%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

%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.