matlab椒盐噪声中值滤波

主要参照了其他网站的一些信息，结合得到以下结果：

a=imread('mengnalisha.jpg');
b=rgb2gray(a);
subplot(241),imshow(a),title('原图');
subplot(242),imshow(b),title('灰度图');
c=imnoise(b,'salt & pepper',0.02);
subplot(243),imshow(c),title('添加椒盐噪声图像');
x=ones(3,3)/9;    %定义一个3*3的均值滤波器
d=imfilter(c,x);  %滤波
subplot(244),imshow(d),title('均值滤波后的图像');
e=medfilt2(c);    %中值滤波
subplot(245),imshow(e),title('中值滤波后的图像');
f=wiener2(c,[3 3]);  %维纳滤波，[3 3]是设置的邻域大小
subplot(246),imshow(f),title('维纳滤波后的图像');

[C1,S1] = wavedec2(c,2,'coif3');  %提取小波分解中第一层的低频图像，即实现了低通滤波去噪
n = [1,2]; %设置尺度向量n
p = [10.12,23.28]; %设置阈值向量p
nc1 = wthcoef2('h',C1,S1,n,p,'s');%对三个方向h,v,d(水平、垂直、对角线）高频系数进行软阈值处理
nc1 = wthcoef2('v',nc1,S1,n,p,'s');
nc1 = wthcoef2('d',nc1,S1,n,p,'s');
x1 = waverec2(nc1,S1,'coif3');  %对新的小波分解结构[c,s]进行重构
Z1=uint8(x1);                  %将double型改为unit8型显示图像
subplot(247),imshow(Z1),title('小波软阈值处理后的图像');
nc2 = wthcoef2('h',C1,S1,n,p,'h');%三个方向硬阈值处理
nc2 = wthcoef2('v',nc2,S1,n,p,'h');
nc2= wthcoef2('d',nc2,S1,n,p,'h');
x2 = waverec2(nc2,S1,'coif3'); 
Z2=uint8(x2);
subplot(248),imshow(Z2),title('小波硬阈值处理后的图像');

[a1,a2]=psnr(Z1,b);  %求软阈值去噪的峰值信噪比a1和信噪比a2
[a3,a4]=psnr(Z2,b);  %求硬阈值去噪的峰值信噪比a1和信噪比a2
% RMSE=sqrt((sum((c-b).^2))./2);

%软阈值结果评价
B=double(b);
A=B-x1;
MSE = sum(A(:).*A(:))/numel(B);  %均方误差MSE，numel()函数返回矩阵元素个数
RMSE=sqrt(MSE);
SNR = 10*log10(sum(B(:).*B(:))/MSE/numel(B));%信噪比SNR
PSNR = 10*log10(255^2/MSE);  %峰值信噪比PSNR

%硬阈值结果评价
C=B-x2;
MSE2 = sum(C(:).*C(:))/numel(B);  
RMSE2=sqrt(MSE2);
SNR2 = 10*log10(sum(B(:).*B(:))/MSE2/numel(B));
PSNR2 = 10*log10(255^2/MSE2); 

%均值滤波结果评价
D=double(d);
D1=B-D;
MSE3 = sum(D1(:).*D1(:))/numel(B);  
RMSE3=sqrt(MSE3);
SNR3 = 10*log10(sum(B(:).*B(:))/MSE3/numel(B));
PSNR3 = 10*log10(255^2/MSE3); 

%中值滤波结果评价
E=double(e);
E1=B-E;
MSE4 = sum(E1(:).*E1(:))/numel(B);  
RMSE4=sqrt(MSE4);
SNR4 = 10*log10(sum(B(:).*B(:))/MSE4/numel(B));
PSNR4 = 10*log10(255^2/MSE4); 

%维纳滤波结果评价
F=double(f);
F1=B-F;
MSE5 = sum(F1(:).*F1(:))/numel(B);  
RMSE5=sqrt(MSE5);
SNR5 = 10*log10(sum(B(:).*B(:))/MSE5/numel(B));
PSNR5 = 10*log10(255^2/MSE5);

结果：

关于以上方法，设置的尺度向量和阈值向量不太理解，于是我又在网上找到了另一种软阈值去噪的方法：

%方法2（该方法得到的软阈值的信噪比低）
[C1,S1] = wavedec2(c,3,'db3');  
%求取阈值  
N = numel(c);                                                  %求图像的像素个数
[chd1,cvd1,cdd1] = detcoef2('all',C1,S1,1);    %[H，V，D]=detcoef2（'all'，C，S，N）返回N级的水平H、垂直V和对角线D细节系数。
cdd1=cdd1(:)';                                                 %m=[1 2;3 4]，m=m(:)，变为[1;3;2;4]，4x1的矩阵.后面的'表示求矩阵转置，m变为 [1,3,2,4]，1x4的矩阵
% sigma = median(abs(cdd1))/0.6745;           %提取细节系数求中值并除以0.6745,abs求取数据的绝对值
 sigma=std(cdd1,1);
thr = sigma*sqrt(2*log(N));                              %求阈值

X1 = wthresh(C1,'h',thr);%硬阈值去噪  
X2 = wthresh(C1,'s',thr);%软阈值去噪 

%小波重建  
xchard = waverec2(X1,S1,'db3');  
xcsoft = waverec2(X2,S1,'db3'); 

y1=uint8(xcsoft);
subplot(247),imshow(y1),title('小波软阈值处理后的图像');
y2=uint8(xchard);
subplot(248),imshow(y2),title('小波硬阈值处理后的图像');

将原图加入高斯噪声：

c=imnoise(b,'gaussian',0.02);
subplot(243),imshow(c),title('添加高斯噪声图像');

结果：

各种方法得到的精度评价结果此处不再展示。