%This code can only be used for research purpose only. 
%Please cite the following paper if you use it:
%L. Zhang et al, “Multiscle LMMSE-based image denoising with optimal wavelet selection,” 
%IEEE Trans. on Circuits and Systems for Video Technology, vol. 15, pp. 469-481, April 2005.


%Note:
%1. You need Matlab Wavelet Toolbox to run the code.
%2. In this code, we suppose the input image is a nxn square matrix (just for some convenience:).
%   You may easily revise the code for non-square nxm images.

clear;
ima=imread('lena.tif','tif');
n=length(ima);%% We suppose ima is a nxn square image. You may easily revise the code for nxm images.
ima=double(ima);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v=25;
noi=v*randn(n,n);
iman=ima+noi;
figure(1);clf;
imshow(iman,[0 255]);

pima=mean(mean(ima.^2));
snro=10*log10(pima/(v^2))

%%%%%%%%%%%%%%%%%%%%%%%%%

[ld,hd,lr,hr]=wfilters('bior1.3');
wbase=2;
%%%%wbase:
%1-bior1.1; 2-bior1.3; 3-bior2.2; 4-bior2.4
%5-bior3.3; 6-db2; 7-db3; 8-db4

J=3;
[S,HW,WH,WW,etl] = ocwt2d(iman,ld,hd,J);
%%%%%%%%denss%%%%%%%%%%%%%%%%%%%
s=[4 6 8 10];
rWW=denss(WW,v,wbase,1,s);
rWH=denss(WH,v,wbase,2,s);
rHW=denss(HW,v,wbase,3,s);
rima=iocwt2d(S,rHW,rWH,rWW,etl,lr,hr);
err=ima-rima;
perr=mean(mean(err.^2));

snrss=10*log10(pima/perr)
figure(2);clf;
imshow(rima,[0 255]);