%=================================================================================
%
% Interactive image denoisinig demo using 
% standard wavelet soft and hard thresholding.
% 
% 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=5;
 [y,L] = mdwt(Xnoisy,h,L);
 
  HH1 = y(N/2+1:N,N/2+1:N);
  HL1 = y(1:N/2,N/2+1:N); 
  LH1 = y(N/2+1:N,1:N/2); 
  
  HH2 = y(N/4+1:N/2,N/4+1:N/2);
  HL2 = y(1:N/4,N/4+1:N/2); 
  LH2 = y(N/4+1:N/2,1:N/4); 
  
  HH3 = y(N/8+1:N/4,N/8+1:N/4);
  HL3 = y(1:N/8,N/8+1:N/4); 
  LH3 = y(N/8+1:N/4,1:N/8); 
  
  HH4 = y(N/16+1:N/8,N/16+1:N/8);
  HL4 = y(1:N/16,N/16+1:N/8); 
  LH4 = y(N/16+1:N/8,1:N/16); 
  
  HH5 = y(N/32+1:N/16,N/32+1:N/16);
  HL5 = y(1:N/32,N/32+1:N/16); 
  LH5 = y(N/32+1:N/16,1:N/32); 

  sigma=Thr_univ(HH1);

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

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

    MultipFac=input('Multiplication factor for the universal threshold (default 1): ');
    
    HL5p=rem_noise_WTHresh(HL5,sigma,MultipFac,THRrule);
    LH5p=rem_noise_WTHresh(LH5,sigma,MultipFac,THRrule);
    HH5p=rem_noise_WTHresh(HH5,sigma,MultipFac,THRrule);

    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

    yp=y;

    yp(N/2+1:N,N/2+1:N)=HH1p;
    yp(1:N/2,N/2+1:N)=HL1p; 
    yp(N/2+1:N,1:N/2)=LH1p; 
  
    yp(N/4+1:N/2,N/4+1:N/2)=HH2p;
    yp(1:N/4,N/4+1:N/2)=HL2p; 
    yp(N/4+1:N/2,1:N/4)=LH2p; 
  
    yp(N/8+1:N/4,N/8+1:N/4)=HH3p;
    yp(1:N/8,N/8+1:N/4)=HL3p; 
    yp(N/8+1:N/4,1:N/8)=LH3p; 
     
    yp(N/16+1:N/8,N/16+1:N/8)=HH4p;
    yp(1:N/16,N/16+1:N/8)=HL4p; 
    yp(N/16+1:N/8,1:N/16)=LH4p; 
 
    y(N/32+1:N/16,N/32+1:N/16)=HH5p;
    y(1:N/32,N/32+1:N/16)=HL5p; 
    y(N/32+1:N/16,1:N/32)=LH5p; 
 
    [Xden,L] = midwt(yp,h,L);

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