close all

% u?itamo sliku
[img,map]=imread('angio0.tif');
figure(1)
imagesc(img)
colormap(gray)

%J = IMNOISE(I,'gaussian',M,V) adds Gaussian white noise of mean M and
%    variance V to the imagesc I. When unspecified, M and V default to 0 and
%    0.01 respectively.

%zasumljna slika
imgGN=imnoise(img,'gaussian'); % dodajemo bijeli šum

figure(2)
imagesc(imgGN)
colormap(gray)
var  = 0.01; % default value

%wavelet transformacija(prva) ulaznog signala
[a1,h1,v1,d1] = swt2(double(img),1,'db1');
[a2,h2,v2,d2] = swt2(a1,2,'db1');
% [a3,h3,v3,d3] = swt2(a2,3,'db1');

%razliite metode raunanja optimalne razine praga
% thr1 = thselect(a, 'rigrsure');
% thr2 = thselect(a, 'heursure');
% thr3 = thselect(a, 'sqtwolog');
% thr4 = thselect(a, 'minimaxi');

h1_thresh = wthresh(h1,'s',0.01);
h2_thresh = wthresh(h2,'s',0.01);
% h3_tresh = wthresh(h3,'s',0.01);

v1_thresh = wthresh(v1,'s',0.01);
v2_thresh = wthresh(v2,'s',0.01);
% v3_tresh = wthresh(v3,'s',0.01);

d1_thresh = wthresh(d1,'s',0.01);
d2_thresh = wthresh(d2,'s',0.01);
% d3_tresh = wthresh(d3,'s',0.01);

a2r = iswt2(a2,h2_thresh,v2_thresh,d2_thresh,'db1');
a1r = iswt2(a2r,h1_thresh,v2_thresh,d2_thresh,'db1'); 
xr = a1r;

%wavelet transformacija(druga) signala na izlazu iz gornje grane i
%ulaznog signala

%gornja grana 
[a1_2,h1_2,v1_2,d1_2] = swt2(xr,1,'sym4');
[a2_2,h2_2,v2_2,d2_2] = swt2(a1_2,2,'sym4');
% [a3_2,h3_2,v3_2,d3_2] = swt2(a2_2,3,'sym4');

%ulazni signal
[a1_3,h1_3,v1_3,d1_3] = swt2(img,1,'sym4');
[a2_3,h2_3,v2_3,d2_3] = swt2(a1_3,2,'sym4');
% [a3_3,h3_3,v3_3,d3_3] = swt2(a2_3,3,'sym4');

%raunanje wienerovog filtra
broj=a2_2.^2;
naz=broj+(var)^2;
hw=broj./naz;
theta2=a2_3.*hw;

theta2r = iswt2(theta2,h2_3,v2_3,d2_3,'db1');
theta2rr = iswt2(theta2r,h1_3,v1_3,d1_3,'db1'); 

%prikaz slike
figure(3)
imagesc(a2r)
colormap(gray)

figure(4)
colormap(gray)
imagesc(a1r)

figure(5)
colormap(gray)
imagesc(theta2r)

figure(6)
colormap(gray)
imagesc(theta2rr)

snri = 10*log10(sum(sum(abs(fft2(img)).^2))/sum(sum(abs(fft2(double(imgGN)-double(img))).^2)));
snro = 10*log10(sum(sum(abs(fft2(img)).^2))/sum(sum(abs(fft2(abs(theta2rr)-double(img))).^2)));














