% Homomorphic filtering in the frequency domain

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

IMG_FILE = 'PET.png';

% flag
ZEROPAD = 1;

% filters
gamma_L = .3; gamma_H = 2;

COLOR_FLAG=0;
a = imread(IMG_FILE);
if (size(a, 3) > 1),
    % if the image is a color image, use the intensity
    hsv_a = rgb2hsv(a); 
    a = hsv_a(:,:,3);
    COLOR_FLAG=1;
end;

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');

% filter construction
H = zeros(size(fA));
orig = size(H)/2;
%D0 = 40; % cut-off frequency
D0 = min(size(fA)*.2); % 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 computation
H = (gamma_H-gamma_L)*(1-exp(-(D.^2)/(2*D0.^2)))+gamma_L;

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

% inverse transform
B = real(ifft2(fB));
% crop
B = B(1:size(A,1), 1:size(A,2));

if (COLOR_FLAG==1),
    hsv_b(:,:,1) = hsv_a(:,:,1); 
    hsv_b(:,:,2) = hsv_a(:,:,2); 
    hsv_b(:,:,3) = B; 
    B = hsv2rgb(hsv_b);
    % remapping B (that is in double) in 0..255 
    % it is an optional step (depending on the study case)
    if (1),
        minB = min(min(min(B(:,:,:))));
        maxB = max(max(max(B(:,:,:))));
        B = ((B-minB)/(maxB-minB))*255;
    end;
end;
% integer representation of the filtered image
b = uint8(B);
figure; imshow(b); title('Filtered image');
