# wieneruv filtr

close all

obr=double(imread("m33cb.jpg")); # musim pretypovat, jinak by obr byl int8!!!
sizeobr=size(obr);
sizex=sizeobr(1);
sizey=sizeobr(2);


figure(1)
imshow(obr,[0,255]);
max(max(obr))
min(min(obr))


#odezvova fce
# Gaussovska PSF, simuluje sferickou aberaci
[y,x]=meshgrid(linspace(-10,10,sizey),linspace(-10,10,sizex));
aber=0.015;
odezva=exp(-(x.*x+y.*y)/aber);
figure(2)
imshow(odezva,[]) 



#prenosova fce
prenos=ifftshift(fft2(fftshift(odezva)));
# regularizace, abychom nedelili nulou
prenos=prenos./max(max(abs(prenos)))+1e-9;
figure(3)
imshow(abs(prenos),[])



#rozmazani obrazu
spektrum=ifftshift(fft2(fftshift(obr)));
rozmazf=prenos.*spektrum;
rozmaz=ifftshift(ifft2(fftshift(rozmazf)));
figure(4)
imshow(real(rozmaz),[])


# nejprve inverze bez sumu

# transformuji obraz bez sumu
dataf=ifftshift(fft2(fftshift(rozmaz)));
# inverzni filtr
invfilt=1./prenos;
# dekonvoluce
dekonv1=dataf.*invfilt;
obr2=ifftshift(ifft2(fftshift(dekonv1)));
figure(5)
imshow(abs(obr2),[])


# inverze se sumem
# pridam sum k obrazu
rozmaz2=rozmaz+rand(sizeobr)*max(max(rozmaz))*1e-9;
dataf=ifftshift(fft2(fftshift(rozmaz2)));
# dekonvoluce
dekonv2=dataf.*invfilt;
obr2=abs(ifftshift(ifft2(fftshift(dekonv2))));
figure(6)
imshow(abs(obr2),[0,max(max(obr2))])



# Wieneruv filtr, optimalni zaostreni
nsr=1e-8;  # N/S pomer bileho sumu
# Wieneruv filtr
wiener=conj(prenos)./(abs(prenos).^2+nsr);
# dekonvoluce
dekonv3=dataf.*wiener;
obr3=abs(ifftshift(ifft2(fftshift(dekonv3))));
figure(7)
imshow(abs(obr3),[0,max(max(obr3))/1.3])

# kreslim rezy inverznim a Wienerovym filtrem a rez 
# horni inset je horizontalne zvetseny 
# dolni inset ukazuje rez prenosovou funkci
idx=((x==0)&(y>=0));
rezinv=invfilt(idx);
rezwig=wiener(idx);
rezprenos=prenos(idx);
yy=y(idx);
figure(8)
plot(yy,abs(rezinv),"r","LineWidth",2,yy,abs(rezwig),"b","LineWidth",2)
xlim([0,5]); ylim([0,1e4]);
legend("inverzni","Wieneruv");
hax=axes();
set(hax,"position",[0.2,0.57,0.3,0.3]);
plot(yy,abs(rezwig),"b","LineWidth",2)
legend("Wieneruv");
xlim([0,6]); ylim([0,10]);
hax2=axes();
set(hax2,"position",[0.2,0.17,0.3,0.3]);
plot(yy,abs(rezprenos),"g","LineWidth",2)
legend("S(f)")
xlim([0,6]); ylim([-0.01,1]);










