算法对应的子函数代码:
%commonfilt2_1.m
function[y]=commonfilt2_1(x,a,b)
[e f]=size(x);
y=x;
for i=2:e-1
for j=2:f-1
zmed=median([x(i-1,j+1),x(i,j+1),x(i+1,j+1),x(i-1,j),x(i+1,j),x(i-1,j-1),x(i,j-1),x(i+1,j-1)]); %除中心像素外的8个像素的灰度值
m=0;n=0;%窗口内满足开始条件的个数
%3*3的窗口内统计m和n的大小
for i1=i-1:i+1
for j1=j-1:j+1
if i1==i & j1==j
continue;
end
if x(i,j)-x(i1,j1) > b
m=m+1;
end
if x(i1,j1)-x(i,j) > b
n=n+1;
end
end
end
%判断灰度值较大的点
if x(i,j) >= 255-a & x(i,j) <= 255
%if m==0
% y(i,j) = x(i,j);
%end
if m==8
y(i,j) = zmed;
end
if m>0 & m<8
if x(i,j) -zmed < b
y(i,j) = x(i,j);%为了逻辑清楚,这里就不省略了
else%疑是噪声点,需通过相邻像素的相关性进行判断
for i1=i-1:i+1
for j1=j-1:j+1
if i1==1 | i1==n | j1==1 |j1==n %边界点即使满足条件也不能用
continue;
end
if(i1~=1 & i1~=e & j1~=1 &j1~=f & x(i,j)-x(i1,j1) < b)%未越界,找到的新像素满足条件
if x(i1,j1) - median([x(i1-1,j1+1),x(i1,j1+1),x(i1+1,j1+1),x(i1-1,j1),x(i1+1,j1),x(i1-1,j1-1),x(i1,j1-1),x(i1+1,j1-1)]) > b
y(i,j) = zmed;
end
end
end
end
end
end
end
%判断灰度值较小的点
if x(i,j) >=0 & x(i,j) <= a
%if n==0
% y(i,j) = x(i,j);
%end
if n==8
y(i,j) = zmed;
end
if n>0 & n<8
%y(i,j) = zmed;
if zmed-x(i,j) < b
y(i,j) = x(i,j);
end
if zmed-x(i,j) > b%疑是噪声点,需通过相邻像素的相关性进行判断
for i1=i-1:i+1
for j1=j-1:j+1 %寻找灰度值差不多的像素
if i1==1 | i1==n | j1==1 |j1==n %边界点即使满足条件也不能用
continue;
end
if(i1~=1 & i1~=e & j1~=1 &j1~=f & x(i1,j1)- x(i,j) < b)%未越界,找到的新像素满足条件
if median([x(i1-1,j1+1),x(i1,j1+1),x(i1+1,j1+1),x(i1-1,j1),x(i1+1,j1),x(i1-1,j1-1),x(i1,j1-1),x(i1+1,j1-1)])-x(i1,j1) > b
y(i,j) = zmed;
end
end
end
end
end
end
end
end
end
主函数
I=imread('2.png');
I=rgb2gray(I);
J=imnoise(I,'salt & pepper',0.8);
%关于a和b的取值 也是一个研究的热点
a=15;
b=10;
%%%%迭代滤波,反复滤波直至噪声被去除
k1=commonfilt2_1(J,a,b);%a and b 是阈值
%k1=commonfilt2_1(k1,a,b);%a and b 是阈值
%k1=commonfilt2_1(k1,a,b);%a and b 是阈值
%%%%
%k2=commonfilt2_2(J,3,3);%库函数中值滤波方法
subplot(121),imshow( uint8(k1) ),title('一种消除椒盐噪声的有效中值滤波');
%subplot(122),imshow( uint8(k2)),title('库函数中值滤波');
% 计算三种算法的峰值信噪比
[e f]=size(I);
B=8; %编码一个像素用多少二进制位
MAX=2^B-1; %图像有多少灰度级
I=double(I);
k1=double(k1);
%k2=double(k2);
%%%%%% psnr1=
MES1=sum(sum((I-k1).^2))/(e*f);
PSNR1=20*log10(MAX/sqrt(MES1));
%%%%%% psnr2=
%MES2=sum(sum((I-k2).^2))/(e*f);
%PSNR2=20*log10(MAX/sqrt(MES2))