%===========================================================================================
% Filename: denoise.m
%
% Operation: denoisinig greyscale images using spatially adaptive ProbShrink [1] 
% (called GenLik in [2]) implemented with a NON-DECIMATED wavelet transform and under 
% the generalized Laplacian marginal prior.
%
%[1] A. Pizurica and W. Philips, "Estimating probability of presence of a signal of interest
%    in multiresolution single- and multiband image denoising," submitted to IEEE Trans. Image
%    Processing, May 2004.
%
%[2] A. Pizurica, Image denoising using wavelets and spatial context modeling, Ph.D. thesis, 
%    Ghent University, June 2002.
%
% Author: Aleksandra Pizurica, TELIN/Ghent University <Aleksandra.Pizurica@telin.UGent.be>
% Last modified: 5. May 2004.
%=============================================================================================

clear all;
close all;

Ref=input('Reference noise-free image available? (y/n) ','s');

% Read reference, noise-free image
if (Ref=='y')
    InputType=0;
    while (InputType~=1 & InputType~=2 & InputType~=3 & InputType~=4)  
     InputType=input('\nReference image format:\n1. "*.tif"\n2. "*.sg2"\n3. "*.mat"\n4. "*.bin"\n>>  ');
    end
    imload;
    Xclean=SCAN;
    clear SCAN;
    [M,N]=size(Xclean);
  
    Choice=input('Choose 1 or 2:\n1. An arbitrary noise variance\n2. A stored noisy image\n ');

    if(Choice==1)
       NoiseStd=input('Standard Deviation of added noise: ');
       Xnoisy=Xclean+randn(M,N)*NoiseStd;
       clean_image=input('Reference image name: ','s');
       image_name=strcat(clean_image,'_std',num2str(NoiseStd));
    else
       InputType=0;
       while (InputType~=1 & InputType~=2 & InputType~=3 & InputType~=4)  
         InputType=input('\nNoisy image format:\n1. "*.tif"\n2. "*.sg2"\n3. "*.mat"\n4. "*.bin"\n>>  ');
       end
       imload;
       Xnoisy=SCAN;
       clear SCAN;
  
       %Extract the image name from the filename
       image_name=filename(1:(length(filename)-4))
       SigmaInfo=input('Noise standard deviation known? (y/n) ','s');
       if(SigmaInfo=='y')
           NoiseStd=input('Standard Deviation of added noise: ');
       else
          [clp,chor,cver,cdiag]=dwt2(Xnoisy,'db1');
          NoiseStd=Thr_univ(cdiag);
      end
    end
  
else % No reference noise-free image image available

    % Read noisy image
    InputType=0;
    while (InputType~=1 & InputType~=2 & InputType~=3 & InputType~=4)  
        InputType=input('\nNoisy image format:\n1. "*.tif"\n2. "*.sg2"\n3. "*.mat"\n4. "*.bin"\n>>  ');
    end
    imload;
    Xnoisy=SCAN;
    clear SCAN;
    imagesc(Xnoisy);
    colormap(gray);

    %Extract the image name from the filename
    image_name=filename(1:(length(filename)-4))
    [clp,chor,cver,cdiag]=dwt2(Xnoisy,'db1');
    NoiseStd=Thr_univ(cdiag);
    EstimatedNoiseStd=NoiseStd
end
 
[M,N]=size(Xnoisy); 

figure
imagesc(Xnoisy);
colormap(gray);
title('Input image');


%======================================================
% Parameters

WindowSize=7;
WindowSizeOPT=input('Window size is set to default 7. Choose other (y/n)? ','s');

if (WindowSizeOPT=='y')
  WindowSize=0;
  while(WindowSize~=1 & WindowSize~=3 & WindowSize~=5 & WindowSize~=7)
     WindowSize=input('window size: 1, 3, 5 or 7: ');
  end
end

W=(WindowSize-1)/2;

SigIntDEF=1;

SigIntOPT=input('DEFAULT/other signal of interest (any key/other) ','s');

if(SigIntOPT=='other')
   SigIntDEF=input('Multiplication factor for the threshold that define the signal of interest (0.5-1): ');
end

%=====================================================  
 % Decomposition 
 
 % Choose a wavelet
 wtypeN=0;
  while((wtypeN~=1) & (wtypeN~=2) & (wtypeN~=3) & (wtypeN~=4) & (wtypeN~=5) & (wtypeN~=6))
    wtypeN=input('Wavelet selection:\n 1 (db2)\n 2 (db4)\n 3 (db8)\n 4 (sym2)\n 5 (sym4)\n 6 (sym8)\n enter your choice 1-6: ');
 end
 
 if(wtypeN==1)
    wtype='db2'
    h = daubcqf(4,'min');  
 elseif(wtypeN==2)
    wtype='db4'
    h = daubcqf(8,'min');
 elseif(wtypeN==3)
    wtype='db8'
    h = daubcqf(16,'min');
 elseif(wtypeN==4)
    wtype='sym2'
    [LO_D,HI_D,LO_R,HI_R] = wfilters('sym2');
    h=LO_R;
 elseif(wtypeN==5)
    wtype='sym4'
    [LO_D,HI_D,LO_R,HI_R] = wfilters('sym4');
    h=LO_R;
 elseif(wtypeN==6)
    wtype='sym8'
    [LO_D,HI_D,LO_R,HI_R] = wfilters('sym8');
    h=LO_R;
 end
 
 
 t=clock;
 
 L=4;
 [yl,yh,L] = mrdwt(Xnoisy/4,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);
 
 
% Due to rescaling of gray values

sigma=NoiseStd/4;
T=SigIntDEF*1.5*sigma;

%====================================
% Processing scale 2^4

Scale=4


HL4p=rem_noise_GenLap(HL4,W,sigma,T);
LH4p=rem_noise_GenLap(LH4,W,sigma,T);
HH4p=rem_noise_GenLap(HH4,W,sigma,T);


%====================================
% Processing scale 2^3


Scale=3


HL3p=rem_noise_GenLap(HL3,W,sigma,T);
LH3p=rem_noise_GenLap(LH3,W,sigma,T);
HH3p=rem_noise_GenLap(HH3,W,sigma,T);


%====================================
%====================================
% Processing scale 2^2

Scale=2


HL2p=rem_noise_GenLap(HL2,W,sigma,T);
LH2p=rem_noise_GenLap(LH2,W,sigma,T);
HH2p=rem_noise_GenLap(HH2,W,sigma,T);

%====================================
% Processing scale 2^1

Scale=1


HL1p=rem_noise_GenLap(HL1,W,sigma,T);
LH1p=rem_noise_GenLap(LH1,W,sigma,T);
HH1p=rem_noise_GenLap(HH1,W,sigma,T);



%================================================
% 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);
Xden=Xden*4;

Time=etime(clock,t)  



figure
imagesc(Xden);
colormap(gray);
title('Denoised image');
   
   
if(Ref=='y')
  Ps=std2(Xclean);
  Pn=std2(Xclean-Xnoisy);
  Prn=std2(Xclean-Xden);

  snr_input=20*log10(Ps/Pn)
  snr_res=20*log10(Ps/Prn)
  
  MSE=mean(mean((Xclean-Xden).^2));
  MSEin=mean(mean((Xclean-Xnoisy).^2));


PSNR=10*log10(255^2/MSE)
PSNRin=10*log10(255^2/MSEin)
end

   
%============================================
% Save result

Save_image=input('Save in mat format? y/n ','s')

if (Save_image=='y');
   Filename=strcat(image_name,'_genlikGL','_',num2str(wtype),'.mat');
   X=Xden;
   WindowSize=W;
   if(Ref=='y')
      save(Filename,'X','snr_input','snr_res','PSNR','PSNRin','NoiseStd', 'MSE','MSEin','WindowSize','wtype','SigIntDEF');
   else
      save(Filename,'X','WindowSize');
   end 
end

Save_tif=input('Save in tif format? y/n ','s')

if (Save_tif=='y');
    Xden(Xden<0)=0;
    Xden(Xden>255)=255;
    Filename=strcat(image_name,'_genlikGL','_',num2str(wtype),'.tif');
    imwrite(uint8(Xden),gray(256),Filename,'tif');
end
