%Please cite the following paper if you use this code for research purpose:

%P. Bao and L. Zhang, “Noise Reduction for Magnetic Resonance Images via Adaptive Multiscale Products Thresholding,” 
%IEEE Trans. on Medical Imaging, vol.22, pp. 1089-1099, Sep. 2003.

clear;

image=imread('house.tif','tif');
image=double(image);
[n,nn]=size(image);%%n must = nn
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1),clf;
imshow(image,[0 255]);
%%%%%%
pimage=mean(mean(image.^2));
v=28;
noi=v*randn(n,n);
snr_o=10*log10(pimage/v^2)  %%the original SNR

sca=3;
imagn=image+noi;
[wr,wc,ss,m]=wt2d(imagn,sca);
[corr,corc]=compcor2d(wr,wc);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rimat=th_w2d(wr,wc,ss,v,m);%%%do the traditional thresholding 

figure(2),clf;
imshow(rimat,[0 255]);

errt=image-rimat;
perrt=mean(mean(errt.^2));
snr_ft=10*log10(pimage/perrt) %% the SNR after traditional thresholding
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% the scale correlation thresholding
[sca,n,nn]=size(wr);

th=getth2d(sca);
c=15;
for i=1:sca-1
   wrt=reshape(wr(i,:,:),n,n);
   corrt=reshape(corr(i,:,:),n,n);
   mask=(corrt>c*v^2*th(i));
   iwr(i,:,:)=wrt.*mask;
   
   wct=reshape(wc(i,:,:),n,n);
   corct=reshape(corc(i,:,:),n,n);
   mask=(corct>c*v^2*th(i));
   iwc(i,:,:)=wct.*mask;
end

nf=normf2d(sca);
wrt=reshape(wr(sca,:,:),n,n);
iwr(sca,:,:)=wrt.*(abs(wrt)>3.*v*nf(sca));
wct=reshape(wc(sca,:,:),n,n);
iwc(sca,:,:)=wct.*(abs(wct)>3.*v*nf(sca));

rima=iwt2d(iwr,iwc,ss,m);

err=image-rima;
perr=mean(mean(err.^2));
snr_f=10*log10(pimage/perr) %% the SNR after scale correlation thresholding

figure(3),clf;
imshow(rima,[0 255]);
