消除椒盐噪声的改进滤波算法（李双全等）

算法对应的子函数代码：

%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))