%% DEMO: 2D ideal filtering

%% Read and draw a picture (Matlab demo picture)
I  = imread('cameraman.tif');   % type 'uint8' = unsigned integer 8bit
figure(1); clf; 

imshow(I); axis on;
N  = size(I, 1);                % size is (N x N)

%% Compute 2D-DFT and plot spectrum
Id = double(I);                 % convert to 'double' type
IF = fft2(Id);                  % 2-D discrete Fourier transform
figure(2); clf; 
imagesc(log10(abs(IF)));        % axis 0..2pi, 0..2pi like in Figure (c)
%imagesc(log10(abs(fftshift(IF)))); % shifted like in Figure (d)
colormap(gray);                 % grayscale
axis equal                      % interval of x and y equal

%% Filter mask of an ideal lowpass filter
wc = 0.3*pi;                    % cut-off frequency wc, 0 < wc < pi
M  = round(wc*N/(2*pi));        % corresponding index, (M/N) == (wc/2pi)
filterMask = zeros(size(I));    % initialize to zeros
filterMask(1:M,          1:M)         = 1;  % left-top corner, 
filterMask(end-M+2:end,  1:M)         = 1;  % left-bottom corner, ...
filterMask(end-M+2:end,  end-M+2:end) = 1; 
filterMask(1:M,          end-M+2:end) = 1; 
figure(3); clf;
mesh(filterMask); colorbar;
%mesh(fftshift(filterMask)); colorbar  % shifted like in Figure (d)
%colormap(gray);                % grayscale

%% 2-D ideal filtering
IFfilt = IF .* filterMask;      % ideal filtering in transform domain
figure(4); clf; 
imagesc(log10(abs(IFfilt))); 
%imagesc(log10(abs(fftshift(IFfilt)))); % shifted like in Figure (d)
colormap(gray);

%% Inverse 2D-DFT and plot figure
Ifiltered = real(ifft2(IFfilt));  % may contain small complex values
figure(5); clf; 
imshow(Ifiltered, [min(Ifiltered(:)) max(Ifiltered(:))]);
