严重椒盐噪声污染图像的非线性滤波算法（董继扬）

新算法对应的子函数：

%delseriousnoise.m
%有效去除图像中脉冲噪声的新型滤波算法
%检测窗口5*5，滤波窗口3*3
function [R] = delseriousnoise(x)
[e,f]=size(x);
Nmax=3; %确定最大的滤波半径，那么滤波窗口或者检测窗口最大为7
%下面是边界扩展，图像上下左右各增加Nmax像素
z=zeros(e+2*Nmax,f+2*Nmax);%引入一个新矩阵，用以贴上图片
z(Nmax+1:e+Nmax,Nmax+1:f+Nmax)=x;%将图片贴在新矩阵的中心位置

%扩展图像边界
z(1:Nmax,Nmax+1:f+Nmax)=x(1:Nmax,1:f);                 %扩展上边界
z(1:e+Nmax,f+Nmax+1:f+2*Nmax)=z(1:e+Nmax,f:f+Nmax-1);  %扩展右边界
z(e+Nmax+1:e+2*Nmax,Nmax+1:f+2*Nmax)=z(e:e+Nmax-1,Nmax+1:f+2*Nmax); %扩展下边界
z(1:e+2*Nmax,1:Nmax)=z(1:e+2*Nmax,Nmax+1:2*Nmax);       %扩展左边界
%z即为一个扩展后的图像
y=z;%y复制z,y为最后输出的图像

F=zeros(e+2*Nmax,f+2*Nmax);%判断每个像素点是否是噪声的二值化矩阵.F(i,j)=0为信号点，F(i,j)=1为噪声点，开始时默认所有点为信号点（毕竟噪声点很少）


Q = 5;
T = 30;
%r = floor(5/2);%检测窗口的大小
r = floor(7/2);

%%%第一步，检测每个像素是否为噪声点
for i=Nmax+1:e+Nmax
  for j=Nmax+1:f+Nmax
    z0=z(i-r:i+r,j-r:j+r);
    z1=z0(1:end);
    if  z(i,j) >= min(z1) & z(i,j) <= min(z1) + Q %负脉冲
      F(i,j) = -1;
    end
    if  z(i,j) >= max(z1) - Q  & z(i,j) <= max(z1) %正脉冲
      F(i,j) = 1;
    end
  end
end

%%%第二步，噪声滤波
N = 5;%迭代的次数
while N ~= 0
  for i=Nmax+1:e+Nmax
    for j=Nmax+1:f+Nmax
      if F(i,j) ~= 0 %若像素为噪声点
         z0=z(i-1:i+1,j-1:j+1);
         z1=z0(1:end);
         medz1 = median(z1);%窗口内所有像素的中值
         M = sum(sum( F(i-1:i+1,j-1:j+1) == 0 ));%窗口内信号点的个数
         if M > 0
           z2 =zeros(1,M);%存放信号点的矩阵
           for i1=i-1:i+1
            for j1=j-1:j+1
              if F(i1,j1) == 0%寻找到信号点
                z2(1,M) = z(i1,j1);
                M = M -1;
              end
            end
           end
           averz2 = mean(z2); %信号点的均值
           if abs( z(i,j) - averz2 ) > T
             y(i,j) = averz2;
           end
         else %信号点集合为空，则使方法来用新的方法来判别
          P = sum(sum( F(i-1:i+1,j-1:j+1) == 1 )) - sum(sum( F(i-1:i+1,j-1:j+1) == -1 ))  %正负脉冲数目之差
          if ~( abs(P) > 3  & P*F(i,j) > 0 )
            y(i,j) = medz1;
          end
         end
        %F(i,j) = 0;%噪声点处理后标记为信号点
      end 
    end
  end
  z=y;%用于下次迭代，这样每次在滤波时，当下像素互不影响
  N = N -1;
end


R(1:e,1:f) = y( Nmax+1:e+Nmax,Nmax+1:f+Nmax );%拷出滤波后的图像

主函数：

I=imread('5.png');
%I=rgb2gray(I);
[e,f]=size(I);

J=imnoise(I,'salt & pepper',0.8);
R=delseriousnoise(J);

%subplot(1,2,1),imshow( uint8(I) ),title('原图像');
%subplot(1,2,2),
imshow( uint8(R) ),title('新算法滤波后的图像');
% 计算三种算法的峰值信噪比
B=8;                %编码一个像素用多少二进制位
MAX=2^B-1;          %图像有多少灰度级

I=double(I);
R=double(R);


%%%%%% psnr1=
MES1=sum(sum((I-R).^2))/(e*f);     %自适应中值去噪的均方差
PSNR1=20*log10(MAX/sqrt(MES1));           %自适应中值算法去噪的峰值信噪比
tic;R=delseriousnoise(J);toc;