% Filtering in the frequency domain

% Operate the filtering in the frequency domain,
% using the process defined in the Gonzalez Woods textbook.

%IMG_FILE = 'white_black.png';
%IMG_FILE = 'car.png';
IMG_FILE = 'PET.png';

% flag
ZEROPAD = 1;

% filters
IDEAL_LOWPASS_FILTER = 1;
BUTTERWORTH_LOWPASS_FILTER = 3;
GAUSS_HIGHPASS_FILTER = 10;
GAUSS_HOMOMORPHIC_FILTER = 20;
NOTCH_FILTER = 30;

% filter selection
%filter = IDEAL_LOWPASS_FILTER
filter = BUTTERWORTH_LOWPASS_FILTER; butt_order = 1;
%filter = GAUSS_HIGHPASS_FILTER
%filter = NOTCH_FILTER; notch_x = 60; notch_y = 80; sigma = 30; ...
%         attenuation = .99;

a = imread(IMG_FILE);
a = a(:,:,1); % if the image is a color image

figure; imshow(a); title('Original image');

%disp('pause'); pause;

A = double(a);

if ZEROPAD,
        fA = fft2(A, 2*size(A,1), 2*size(A,2)); % zero padding
else
        fA = fft2(A);
end;

%figure;
%mesh(log(abs(fA)));

figure;
mesh(fftshift(log(abs(fA)+eps)));
view(0,90);
axis equal;
colorbar;
colormap gray;
title('FFT of the image');

% ideal filter
H = zeros(size(fA));
orig = size(H)/2;
D0 = 40; % cut-off frequency

% distance computation with respect to the center of the matrix
[mu mv] = meshgrid(1:size(H,2), 1:size(H,1));
mu = mu-orig(2);
mv = mv-orig(1);
D = sqrt(mu.^2+mv.^2);

% filter construction
if filter == IDEAL_LOWPASS_FILTER,
  [i j] = find(D<D0);
  H(i,j) = 1;
elseif filter == BUTTERWORTH_LOWPASS_FILTER, 
    H = 1./(1+((D/D0).^(2*butt_order)));
elseif filter == GAUSS_HIGHPASS_FILTER, 
  H = 1-exp(-(D.^2)/(2*D0.^2));
elseif filter == NOTCH_FILTER,
    % notch filter for car.png
    % decrease the components near (+/-notch_x, +/-notch_y) 
    H = ones(size(H));
    %(+notch_x, +notch_y)
    H = H-attenuation*exp(-((mu-notch_x).^2+(mv-notch_y).^2)./sigma^2);
    %(+notch_x, -notch_y)
    H = H-attenuation*exp(-((mu-notch_x).^2+(mv+notch_y).^2)./sigma^2);
    %(-notch_x, -notch_y)
    H = H-attenuation*exp(-((mu+notch_x).^2+(mv+notch_y).^2)./sigma^2);
    %(-notch_x, +notch_y)
    H = H-attenuation*exp(-((mu+notch_x).^2+(mv-notch_y).^2)./sigma^2);
end;

figure;
mesh(mu,mv,H);
title('Filter');

% switch the filter quadrants
H = fftshift(H);

disp('pause'); pause;

% filtering in the frequency domain
fB = fA.*H;

figure;
mesh(fftshift(log(abs(fB)+eps)));
view(0,90);
axis equal;
colorbar;
colormap gray;
title('FFT of the filtered image');

% inverse transform
B = real(ifft2(fB));

% integer representation of the filtered image
b = uint8(B(1:size(A,1), 1:size(A,2)));
figure; imshow(b); title('Filtered image');
