% filtrace spektra 1
% horni a dolni propust, filtrace periodickych struktur

function sfa6_1

close all

figure(1);
predmet=imread("ctyrka.jpg");
predmet=shift(predmet',30)'; % centrovani
dim=size(predmet);
fxrange=linspace(-1,1,dim(2));
fyrange=linspace(-1,1,dim(1));
[fx,fy]=meshgrid(fxrange,fyrange);  % libovolne jednotky
x=fx; y=fy;

subplot(2,3,1)
imshow(predmet);

subplot(2,3,4);
g=ft2(predmet,1);
si=abs(g);  % spektralni intenzita
imshow(si,[0,max(max(si))/10])


%dolni propust

subplot(2,3,5);
Hd=exp(-200*(fx.*fx+fy.*fy));
gd=g.*Hd;   %potlaceni vysokych frekvenci 
sid=abs(gd);  
imshow(sid,[0,max(max(sid))/10])

subplot(2,3,2);
vystup=abs(ift2(gd,1)); %vysledna intenzita
imshow(vystup,[])

% horni propust

subplot(2,3,6);
Hh=exp(4*(fx.*fx+fy.*fy));
gh=g.*Hh; % zesileni vysokych frekvenci
sih=abs(gh);  
imshow(sih,[0,max(max(sih))/10])

subplot(2,3,3);
vystup=abs(ift2(gh,1)); % vysledna intenzita
imshow(vystup,[])


% prenosove funkce

figure(2)
stred=ceil(dim(1)/2+1);
hd=real(ift2(Hd,1));  hd=hd/norm(hd); 
hh=real(ift2(Hh,1));  hh=hh/norm(hh); 
plot(fx(stred,:),real(hd(stred,:)),"Linewidth",2,fx(stred,:),real(hh(stred,:)),"Linewidth",2);
xlim([-0.15,0.15]); ylim([-0.4,0.5]);
legend("dolni","horni")



% selektivni filtrace


%vytvorim mriz
N=64;
predmet=zeros(N);
for i=1:8:N-2
 predmet(i:i+1,:)=ones(2,N);
 predmet(:,i:i+1)=ones(N,2);
endfor


dim=size(predmet);
fxrange=linspace(-1,1,dim(1));
fyrange=linspace(-1,1,dim(2));
[fx,fy]=meshgrid(fxrange,fyrange);
predmet=predmet.*exp(-fx.*fx-fy.*fy);  % pokles intenzity k okrajum

figure(3)
subplot(2,3,1)
imshow(predmet)

subplot(2,3,4)
g=ft2(predmet,1);
si=abs(g);  % spektralni intenzita
imshow(si,[0,max(max(si))/100])

% ponecham pouze frekvence fx<<1
% tj. odstranim horizontalni detaily

subplot(2,3,5)
gf1=g.*rect(fx,0.2);
sif1=abs(gf1);  
imshow(sif1,[0,max(max(sif1))/10])

subplot(2,3,2);
vystup=abs(ift2(gf1,1));
imshow(vystup,[])

% ponecham frekvence fy<<1
% tj. odstranim vertikalni detaily

subplot(2,3,6)
gf2=g.*rect(fy,0.2);
sif2=abs(gf2);  
imshow(sif2,[0,max(max(sif2))/10])

subplot(2,3,3);
vystup=abs(ift2(gf2,1));
imshow(vystup,[])

end




function G=ft2(g,del)
 G=fftshift(fft2(fftshift(g)))*del*del;
end


function g=ift2(G,delf)
 g=ifftshift(ifft2(ifftshift(G)))*(delf*size(G,1))^2;
end

function z=circ(x,y,D)
 r=sqrt(x.^2+y.^2);
 z=double(r<D/2);
 z(r==D/2)=0.5;
end

function y=rect(x,D)
  if (nargin==1) D=1; end
  x=abs(x);
  y=double(x<D/2);
  y(x==D/2)=0.5;
end


