flt='cl';fltap='clap'; %使用的多小波和预滤波器
deimg='c:\denoise\cameraman512.bmp'; %要处理的图像
img=imread(deimg); %从c:\noise目录下读入图像
%cc=0.0384;img1=imnoise(img,'gaussian',0,cc);
%在某些论文中,都说加入了sigma=??的噪声,例如sigma=10等等,
%见论文<基于小波变换的自适应多阀值图像去噪,中国图形图像学报,2005年第5期,表1和表2>,
%我认为,应该用下面一行代码就可在MATLAB中实现,这是因为,
%若数列Xn的标准差是1,则由D(cc*Xn)=cc^2*D(cc*Xn)=cc^2,得cc*Xn标准差是cc,
%其中函数randn(m,n)生成均值是0,标准差是1的正态分布的mxn矩阵,
%那么,它与使用imnoise(img,'gaussian',0,sigma_square)有什么关系?
%imnoise中的sigma_square是方差,0=sigma_square=1,而cc是标准差,0=cc=255,
%两者的关系是: sigma_square=(cc/255)^2 或者 sigma=cc/255,
%例如:cc=10对应sigma=10/255, 也对应sigma_square=(10/255)^2
%另外,下面的这行程序已经变成一个此目录下的子程序,如下:
% noise_img=addnoise(img,sigma)
%上面的函数对于给定的图像矩阵img,加入均值为0,标准差为sigma的符合正态分布的噪声,
%即高斯噪声,它与imnoise(img,'gaussian',0,sigma_square)区别在于,imnoise中的
%sigma_square是规范化后的方差,它在[0,1]之间,而addnoise的sigma是没有规范化的标准差,
%它在[0,255]之间,其两者的转换关系见上面的说明,你可以这样理解这两个函数的不同,
%对于取值于[0,255]的灰度图像矩阵img,函数imnoise先将img中的每个元素变化到[0,1]之间,
%然后加入噪声,而函数addnoise则在img中直接加入噪声.
cc=50;
[m,n]=size(img); img1=double(img); img1=img1+cc*randn(m,n);
%disp(rmse(img,img1));
pre=prep2D_appe(img1,fltap); %预滤波
trans_num=3;w=dec2D_pe(pre,flt,trans_num);
%进行trans_num次多小波变换,第1次图像分成4x4块,每块(m/4)x(n/4)
%第2次变换,左上角的4个图像块,被变换成4x4的块,每块(m/16)x(n/16)
%下面对每次变换后的细节系数使用阀值算法修改块内的小波系数
au=0;bu=0;num=1;t=0.6;rrr=4;
%disp(max(max(abs(w))));
while num=trans_num+1
for i=1:m/4:m-m/4+1
for j=n/2+1:n/4:n-n/4+1
mat=w(i:i+m/4-1,j:j+n/4-1);dd=dim(mat);
xx=diffx(diffx(mat));yy=diffy(diffy(mat));
%fprintf('(%f,%f),(%f,%f)\n',i,j,i+m/4-1,j+n/4-1);
[mm,nn]=size(mat); %ee=energy(mat);disp(ee);
u0=(max(max(mat))-min(min(mat)))/rrr; qq=1/(mm*nn);
for ii=1:mm
for jj=1:nn
%u=u0*exp(-sqrt(log10(1+t*abs(xx(ii,jj)+yy(ii,jj)))));
u=u0*dd*exp(-qq*t*abs(xx(ii,jj)+yy(ii,jj))^2);
%if num=10 fprintf('%f*',u);num=num+1;else fprintf('%f\n',u);;num=0;end
if mat(ii,jj)>u
mat(ii,jj)=mat(ii,jj)-u*t;au=au+1;
elseif mat(ii,jj)<-u
mat(ii,jj)=mat(ii,jj)+u*t;au=au+1;
else
mat(ii,jj)=0;bu=bu+1;
end
end
end
w(i:i+m/4-1,j:j+n/4-1)=mat; %pause;
end
end
for i=m/2+1:m/4:m-m/4+1
for j=1:n/4:n/2-n/4+1
mat=w(i:i+m/4-1,j:j+n/4-1);dd=dim(mat);
xx=diffx(diffx(mat));yy=diffy(diffy(mat));
%fprintf('(%f,%f),(%f,%f)\n',i,j,i+m/4-1,j+n/4-1);
[mm,nn]=size(mat); %ee=energy(mat);
u0=(max(max(mat))-min(min(mat)))/rrr; qq=1/(mm*nn);
for ii=1:mm
for jj=1:nn
%u=u0*exp(-sqrt(log10(1+t*abs(xx(ii,jj)+yy(ii,jj)))));
u=u0*dd*exp(-qq*t*abs(xx(ii,jj)+yy(ii,jj))^2);
%if num=10 fprintf('%f*',u);num=num+1;else fprintf('%f\n',u);;num=0;end
if mat(ii,jj)>u
mat(ii,jj)=mat(ii,jj)-u*t;au=au+1;
elseif mat(ii,jj)<-u
mat(ii,jj)=mat(ii,jj)+u*t;au=au+1;
else
mat(ii,jj)=0;bu=bu+1;
end
end
end
w(i:i+m/4-1,j:j+n/4-1)=mat; %pause;
end
end
m=m/2;n=n/2;num=num+1;
end
ww=rec2D_pe(w,flt,trans_num);
www=postp2D_appe(ww,fltap);
%www=ordfilt2(www,5,ones(3)); %3x3中值滤波
%www=www-0.01*t*(diffx(diffx(www))+diffy(diffy(www)));
fprintf('%% 滤波前:snr(img,img1)=%f, snr(img,www)=%f\n',snr(img,img1),snr(img,www));
www=ordfilt2(www,5,ones(3));
%www=wiener2(www,[3,3]);
fprintf('%% 下面是图像%s当t=%f,sigma=%f时的处理结果:\n',deimg,t,cc);
fprintf('%% multiwavelets->%s, frefilter->%s, u=0, sigma=%f\n',flt,fltap,cc);
fprintf('%% mse(img,img1)=%f, mse(img,www)=%f\n',mse(img,img1),mse(img,www));
fprintf('%% rmse(img,img1)=%f, rmse(img,www)=%f\n',rmse(img,img1),rmse(img,www));
fprintf('%% psnr(img,img1)=%f, psnr(img,www)=%f\n',psnr(img,img1),psnr(img,www));
fprintf('%% snr(img,img1)=%f, snr(img,www)=%f\n',snr(img,img1),snr(img,www));
fprintf('%% w(i,j>u:%f%%, w(i,j)=u:%f%%\n\n',au/(au+bu)*100,bu/(au+bu)*100);
subplot(121);colormap(gray(256));image(img1);axis('off');
subplot(122);colormap(gray(256));image(www);axis('off');;
