%=================================================================================
%
% Interactive image denoisinig demo using 
% wavelet thresholding in a non-decimated wavelet representation.
% 
% Aleksandra Pizurica, TELIN/UGent, 5. November 2005.
%=================================================================================

clear all;
%close all;

global InputType SCAN;

InputType=1;
imload;
Xnoisy=SCAN;
clear SCAN;
figure('units','Normalized', 'Position',[0.1 0.1 0.4 0.7]);
subplot(2,1,1),imagesc(Xnoisy); title('Input image');
colormap(gray);

THRrule=0;
while((THRrule~=1) & (THRrule~=2))
   THRrule=input('Choose thresholding rule:\n1- soft\n2 - hard: ');
end
 
%=====================================================  
 % Decomposition 
 
 [M,N]=size(Xnoisy);
 wtype='db2'
[LO_D,HI_D,LO_R,HI_R] = wfilters(wtype); 
 h=LO_R;
 L=4;
 [yl,yh,L] = mrdwt(Xnoisy,h,L);
 
 LH1 = yh(:,1:N); 
 HL1 = yh(:,N+1:2*N); 
 HH1 = yh(:,2*N+1:3*N);
 
 LH2 = yh(:,3*N+1:4*N); 
 HL2 = yh(:,4*N+1:5*N); 
 HH2 = yh(:,5*N+1:6*N);
 
 LH3 = yh(:,6*N+1:7*N); 
 HL3 = yh(:,7*N+1:8*N); 
 HH3 = yh(:,8*N+1:9*N);
 
 LH4 = yh(:,9*N+1:10*N); 
 HL4 = yh(:,10*N+1:11*N); 
 HH4 = yh(:,11*N+1:12*N);
 

  sigma=Thr_univ(HH1);

  %====================================
   % Processing

  
 Repeat='y';
 while(Repeat=='y')

    MultipFac=input('Multiplication factor for the universal threshold (default 1): ');

    HL4p=rem_noise_WTHresh(HL4,sigma,MultipFac,THRrule);
    LH4p=rem_noise_WTHresh(LH4,sigma,MultipFac,THRrule);
    HH4p=rem_noise_WTHresh(HH4,sigma,MultipFac,THRrule);

    HL3p=rem_noise_WTHresh(HL3,sigma,MultipFac,THRrule);
    LH3p=rem_noise_WTHresh(LH3,sigma,MultipFac,THRrule);
    HH3p=rem_noise_WTHresh(HH3,sigma,MultipFac,THRrule);

    HL2p=rem_noise_WTHresh(HL2,sigma,MultipFac,THRrule);
    LH2p=rem_noise_WTHresh(LH2,sigma,MultipFac,THRrule);
    HH2p=rem_noise_WTHresh(HH2,sigma,MultipFac,THRrule);

    HL1p=rem_noise_WTHresh(HL1,sigma,MultipFac,THRrule);
    LH1p=rem_noise_WTHresh(LH1,sigma,MultipFac,THRrule);
    HH1p=rem_noise_WTHresh(HH1,sigma,MultipFac,THRrule);

    %================================================
    % Reconstruction

    yhp=yh;

    yhp(:,1:N)=LH1p; 
    yhp(:,N+1:2*N)=HL1p; 
    yhp(:,2*N+1:3*N)=HH1p;
 
    yhp(:,3*N+1:4*N)=LH2p; 
    yhp(:,4*N+1:5*N)=HL2p; 
    yhp(:,5*N+1:6*N)=HH2p;
 
    yhp(:,6*N+1:7*N)=LH3p; 
    yhp(:,7*N+1:8*N)=HL3p; 
    yhp(:,8*N+1:9*N)=HH3p;
 
    yhp(:,9*N+1:10*N)=LH4p; 
    yhp(:,10*N+1:11*N)=HL4p; 
    yhp(:,11*N+1:12*N)=HH4p;
 
    [Xden,L] = mirdwt(yl,yhp,h,L);

    subplot(2,1,2),imagesc(Xden);
    title('Reconstructed image');
    Repeat=input('Repeat with different thresholds (y/n)? ', 's');
    
end
    
   
   
