% Filtraggio nel dominio della frequenza

% Effettua il filtraggio nel dominio della frequenza 
% costruendo il filtro dalla definizione del Gonzalez Woods

IMG_FILE = 'bianco_nero.png';

% flag
ZEROPAD = 1;

% filtri disponibili
FILTRO_LOWPASS_IDEALE = 1
FILTRO_LOWPASS_BUTTERWORTH = 3
FILTRO_HIGHPASS_GAUSS = 10

% selezione del filtro da applicare
%filtro = FILTRO_LOWPASS_IDEALE
filtro = FILTRO_LOWPASS_BUTTERWORTH; butt_order = 1;
%filtro = FILTRO_HIGHPASS_GAUSS

a = imread(IMG_FILE);
a = a(:,:,1); % nel caso sia a colori

figure; imshow(a); title('Immagine originale');

%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 dell''immagine');

% filtro ideale
H = zeros(size(fA));
orig = size(H)/2;
D0 = 40; % frequenza di riferimento

% calcolo della distanza rispetto al centro della matrice
[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);

% generazione del filtro
if filtro == FILTRO_LOWPASS_IDEALE,
  [i j] = find(D<D0);
  H(i,j) = 1;
elseif filtro == FILTRO_LOWPASS_BUTTERWORTH, 
    H = 1./(1+((D/D0).^(2*butt_order)));
elseif filtro == FILTRO_HIGHPASS_GAUSS, 
  H = 1-exp(-(D.^2)/(2*D0.^2));
end;

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

% scambia i quadranti del filtro
H = fftshift(H);

disp('pause'); pause;

% filtraggio nel dominio della frequenza
fB = fA.*H;

% antitrasformata
B = real(ifft2(fB));

% rappresentazione intera dell'immagine filtrata
b = uint8(B(1:size(A,1), 1:size(A,2)));
figure; imshow(b); title('Immagine filtrata');
