%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Multimedia
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script IV
%-------------
%
% Contenuti:  - DFT
%             - Trasformata di Fourier
%             - Filtraggio nel dominio della Frequenza

close all
clear all

y=zeros(256,1);
y(1:20)= 2.;
subplot(2,2,1),plot(y),ylim([0,4] );

yy=fft(y,256); % Calcolo della trasformata per un segnale unidimensionale

subplot(2,2,2),plot(abs(yy)),ylim([0 25]);
yy=fftshift(yy);    %Shifting delle frequenze

subplot(2,2,3),plot(abs(yy)),ylim([0 25]);


%Calcolo della Trasformata 2D
%Esempio1
[x,y]= meshgrid( [ -2:0.1:2] ) ;
z= exp( -0.5.* (x.^ 2+ y.^ 2) ./ 0.9) ;
surf(x,y,z) ;
figure; imagesc(z) ;
yy=fft2(z);
I=(fftshift(abs(yy)));
figure;imshow(I)

%Esempio2 
f=imread('barca.gif');
figure; imshow(f)
F=fft2(double(f));
Fc=fftshift(F);
S=abs(Fc);
c=255/ log(max(S(:)));    %Trasformazione logaritmica
S2=c.* log(1.*S);
%S2=uint8(S2);
figure;
%imshow(unit8(S2));
imshow(S2, []);

%Trasformazione inversa

F1=ifftshift(Fc);    %Inverte lo shift
f1=ifft2(F1);       %Inverte la trasformata
f1=real(f1);        %Prende la parte reale

figure;imshow(f1, []);
%figure;imshow(uint8(f1));

%Filtraggio passabasso/passalto 

clear all
close all

f=imread('fishingboat.bmp');
%f=imread('barca.gif');
figure; imshow(f);
[M,N]=size(f);
F=fft2(double(f));

Fc=fftshift(F);
S=abs(Fc);
c=255/ log(max(S(:)));    %Trasformazione logaritmica
S2=c.* log(1.*S);
figure;
imshow(S2, []);


%Il codice che segue permette di organizzare gli indici per i calcoli
%successivi senza la necessità di usare la funzione fftshift

u=0:(M-1);
v=0:(N-1);
    idx=find(u>M/2);
    u(idx)=u(idx)-M;
    idy=find(v>N/2);
    v(idy)=v(idy)-N;
 [V,U]=meshgrid(v,u);
 

D0=10;                                  %Frequenza di taglio
D=sqrt(U.^2+V.^2);
H_LP_Gauss=exp(-(D.^2)./(2*(D0^2)));     %Filtro Gaussiano con frequenza di taglio=dev=Sig
H_LP_Ideale=double(D<=D0);               %Filtro Passabasso ideale

n=1;                                     %Filtro Passabasso di ButterWorth
H_LP_ButterWorth=1./(1+(D./D0).^(2*n));

%H_LP=H_LP_Gauss;
%H_LP=H_LP_Ideale;
H_LP=H_LP_ButterWorth;
figure;imshow(fftshift(H_LP));

H_HP=1-H_LP;                            % Filtro passa alto
figure;imshow(fftshift(H_HP));

G1=H_LP.*F;
G2=H_HP.*F;

g1=real(ifft2(G1));
figure; imshow(g1, []);

g2=real(ifft2(G2));
figure; imshow(g2, []);

% E' possibile realizzare, qualora servisse un padding di zeri nel calcolo
% della trasformata, in analogia a quanto avveniva nel dominio spaziale. 
% Ciò può tornare utile, quando ad esempio le dimensioni del filtro H, sono
% maggiori di quelle della F

[M,N]=size(f);
F=fft2(f, 2*M,2*N);

u=0:(2*M-1);
v=0:(2*N-1);
    idx=find(u>M/2);
    u(idx)=u(idx)-M;
    idy=find(v>N/2);
    v(idy)=v(idy)-N;
 [V,U]=meshgrid(v,u);
 

D0=10;                                  %Frequenza di taglio1
D1=30;                                  %Frequenza di taglio2
D=sqrt(U.^2+V.^2);
H=double((D>=D0)&(D<=D1));                 %Passa Banda...
g=real(ifft2(H.*F));
g=g(1:size(f,1),1:size(f,2));

figure; imshow(g, []);


% Il Padding può essere utilizzato anche per eguagliare le dimensioni dei
% filtri

close all

H=ones(1000,1000);
H(1:100, 1:100)=0;
H(1:100, 900:1000)=0;
H(900:1000, 900:1000)=0;
H(900:1000, 1:100)=0;
figure;imshow(H);
F=fft2(f, size(H,1),size(H,2));

g=real(ifft2(H.*F));
g=g(1:size(f,1),1:size(f,2));

figure; imshow(g, []);


% Relazione tra filtri spaziali e filtri nel dominio della frequenza

h1=fspecial('sobel')';
h2=fspecial('sobel');
h3=fspecial('average');
h4=[-1, 2, -1; 1, 4, 1; -1 2 -1];

h=h1; %Scelta del filtro spaziale da applicare


freqz2(h,100,100)   




f=imread('fishingboat.bmp');
%f=imread('barca.gif');
figure; imshow(f);
[M,N]=size(f);

F=fft2(f, 2*M, 2*N);
S=fftshift(log(1+abs(F)));
figure; imshow(S, []);

H=freqz2(h,2*M, 2*N);
H1=ifftshift(H);                    %L'origine viene riportata nell'angolo in alto a sx
figure; imshow(abs(H),[]);
figure; imshow(abs(H1), []);

%mesh(abs(H));
%surf(abs(H));
%shading interp;

g_spatial=imfilter(double(f), h);

g_freq=real(ifft2(H1.*F));
g_freq=g_freq(1:size(f,1),1:size(f,2));

figure; imshow(g_freq, []);
figure; imshow(g_spatial, []);

d=abs(g_freq - g_spatial);
max(d(:));
