function [dX, SNRul, SNRizl]=wavelet_hard(X,Xsum,Nsigma, level, wname)
% Ovo je funkcija koja vrsi potiskivanje suma uz pomoc wavelet
% razlaganja i metode praga (HARD threshloding)
%
% Ulazne varijable:
% X            -- originalna slika
% Xsum         -- Slika sa bijelim sumom
% Nsigma       -- Broj koji odreduje iznos praga
% level        -- broj razina razlaganja u wavelet stablu
% wname        -- ime wavelet funkcije
% Izlazne varijable:
% dX           -- rekonstrukcijska slika s uklonjenim sumom
% SNRul        -- Odnos signal-sum na ulazu (prije rekonstrukcije)
% SNRizl       -- Odnos signal-sum na izlazu (poslije rekonstrukcije)



x=size(X,1);
y=size(X,2);



% Izracunavanje efektivne vrijednosti suma
Sum=abs(Xsum-X);
eN=sum(sum(sum(Sum.*Sum)));
MS=eN/(x*y*3);
RMSul=sqrt(MS);

eX=sum(sum(sum(X.*X)));
D1=sqrt(eN/eX);

% Odnos signal/sum za sliku s bijelim sumom
SNRul=-20*log10(D1);


% Potiskivanje suma za crvenu komponentu

% razlaganje na Wavelet filtarske slogove
% dwtmode('per');

[C_R,S]=wavedec2(Xsum(:,:,1),level,wname);
[H0,H1,F0,F1] = wfilters(wname);

diagR = detcoef2('d',C_R,S,1);   % dijagonalni koeficijenti

% Izracunavanje standardne devijacije dijagonalnih komponenata
sigma_R=std(diagR(:));

% vrijednosti praga
prag=Nsigma.*sigma_R;

new_C_R = wthresh(C_R, 'h', prag);     % HARD thresholding


% Rekonstrukcija 
% dwtmode('per');
r_X_R=waverec2(new_C_R,S,wname);

% Limitiranje koeficijenata na [0,255]
for i=1:x
    for j=1:y
        if (r_X_R(i,j)<0)
            r_X_R(i,j)=0;
        elseif (r_X_R(i,j)>255)
            r_X_R(i,j)=255;
        end
    end
end



% Potiskivanje suma za zelenu komponentu

% razlaganje na Wavelet filtarske slogove
% dwtmode('per');
[C_G,S]=wavedec2(Xsum(:,:,2),level,wname);
diagG = detcoef2('d',C_G,S,1);   % dijagonalni koeficijenti

% Izracunavanje standardne devijacije dijagonalnih komponenata
sigma_G=std(diagG(:));

% vrijednosti praga
prag=Nsigma.*sigma_G;

new_C_G = wthresh(C_G, 'h', prag);     % HARD thresholding


% Rekonstrukcija 
% dwtmode('per');
r_X_G=waverec2(new_C_G,S,wname);

% Limitiranje koeficijenata na [0,255]
for i=1:x
    for j=1:y
        if (r_X_G(i,j)<0)
            r_X_G(i,j)=0;
        elseif (r_X_G(i,j)>255)
            r_X_G(i,j)=255;
        end
    end
end



% Potiskivanje suma za plavu komponentu

% razlaganje na Wavelet filtarske slogove
% dwtmode('per');
[C_B,S]=wavedec2(Xsum(:,:,3),level,wname);
diagB = detcoef2('d',C_B,S,1);   % dijagonalni koeficijenti

% Izracunavanje standardne devijacije dijagonalnih komponenata
sigma_B=std(diagB(:));

% vrijednosti praga
prag=Nsigma.*sigma_B;

new_C_B = wthresh(C_B, 'h', prag);     % HARD thresholding


% Rekonstrukcija 
% dwtmode('per');
r_X_B=waverec2(new_C_B,S,wname);

% Limitiranje koeficijenata na [0,255]
for i=1:x
    for j=1:y
        if (r_X_B(i,j)<0)
            r_X_B(i,j)=0;
        elseif (r_X_B(i,j)>255)
            r_X_B(i,j)=255;
        end
    end
end




% Rekonstrukcijska slika

dX(:,:,1)=r_X_R;
dX(:,:,2)=r_X_G;
dX(:,:,3)=r_X_B;


% Izracunavanje efektivne vrijednosti suma na izlazu 
Sum=abs(dX-X);
eN=sum(sum(sum(Sum.*Sum)));
MS=eN/(x*y*3);
RMSizl=sqrt(MS);

% Odnos signal sum na izlazu
D2=sqrt(eN/eX);
SNRizl=-20*log10(D2);


imwrite(uint8(dX),'wavelet_hard_rezultat_db10.bmp','bmp');
