一种高效快速的高密度椒盐噪声消除算法（吕宗伟等）

新算法对应的子函数：

%一种高效快速的高密度椒盐噪声滤除算法
%quickdelhighnoise.m
function [R] = quickdelhighnoise(J)

[e,f]=size(J);
R=zeros(e,f);
Nmax=3; %确定最大的滤波半径，那么滤波窗口最大为7
%下面是边界扩展，图像上下左右各增加Nmax像素
z=zeros(e+2*Nmax,f+2*Nmax);%引入一个新矩阵，用以贴上图片
z(Nmax+1:e+Nmax,Nmax+1:f+Nmax)=J;%将图片贴在新矩阵的中心位置

%扩展图像边界
z(1:Nmax,Nmax+1:f+Nmax)=J(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为最后输出的图像

%%%%%1 噪声检测阶段
F=zeros(e+2*Nmax,f+2*Nmax);%噪声标记矩阵，初始所有像素默认为信号点，标记为0

for i=Nmax+1:e+Nmax
  for j=Nmax+1:f+Nmax
    %检测窗口开始用3*3
    if  z(i,j) == 0 | z(i,j) == 255 %疑似噪声点，需要进一步判断
      F(i,j) = 1;
   end
 end
end


%%%%%% 2 噪声滤除阶段
 for i=Nmax+1:e+Nmax
  for j=Nmax+1:f+Nmax
    %初始滤波窗口为3*3
    if F(i,j) == 1 %噪声点，需要进行滤波 
       N=sum(sum(F(i-1:i+1,j-1:j+1) == 0));%统计高密度区域内信号点的个数
       %N=sum(sum( F(i-1:i+1,j-1:j+1) == 0  &  z(i-1:i+1,j-1:j+1) ~= max(z1) &  z(i-1:i+1,j-1:j+1) ~= min(z1) ))
       if  N > 0
         z6 = zeros(1,N);%放置信号点的矩阵
         %存放信号点矩阵
         for i1 = i-1:i+1
           for j1 = j-1:j+1
            if F(i1,j1) == 0
         %if F(i1,j1) == 0 & z(i1,j1) ~= max(z1) & z(i1,j1) ~= min(z1)
              z6(1,N) = z(i1,j1);
              N = N-1;
           end
          end
         end
         %z6内存放着信号点灰度值
         %for i1 = 1:sum(sum(F(i-1:i+1,j-1:j+1) == 0))
         %for i1 = 1:sum(sum( F(i-1:i+1,j-1:j+1) == 0  &  z(i-1:i+1,j-1:j+1) ~= max(z1) &  z(i-1:i+1,j-1:j+1) ~= min(z1) ))
         %    sum1 = sum1 + z6(1,i1)/( 1 + ( z6(1,i1) - median(z6) )^2 );
         %   sum2 = sum2 + 1.0/( 1 + ( z6(1,i1) - median(z6) )^2 );
         %end
         y(i,j) = mean(z6); 
         %medy4 =median( [ y(i-1,j-1),y(i-1,j),y(i-1,j+1),y(i,j-1) ] );%已经处理过的4领域系数
         %y(i,j) = y(i-1,j-1)/( 1 + ( y(i-1,j-1) - medy4  )^2 ) +  y(i-1,j)/( 1 + ( y(i-1,j) - medy4  )^2 ) +  y(i-1,j+1)/( 1 + ( y(i-1,j+1) - medy4  )^2 ) + y(i,j-1)/( 1 + ( y(i,j-1) - medy4  )^2 ) ;
         %y(i,j) =( y(i-1,j-1) + y(i-1,j) + y(i-1,j+1) + y(i,j-1) ) / 4;
         %y(i,j) = median( [ y(i-1,j-1),y(i-1,j),y(i-1,j+1),y(i,j-1) ] );
         %y(i,j) = y(i-1,j-1)/( 1 + ( y(i-1,j-1) - z(i,j)  )^2 ) +  y(i-1,j)/( 1 + ( y(i-1,j) - z(i,j)  )^2 ) +  y(i-1,j+1)/( 1 + ( y(i-1,j+1) - z(i,j)  )^2 ) + y(i,j-1)/( 1 + ( y(i,j-1) - z(i,j)  )^2 ) ;
       end
       %F(i,j) = 0;%处理完了噪声点
       %z(i,j) = y(i,j);
    end
    %处理完该噪声点后，将处理完的新像素用来迭代
  end
 end
 %第一步滤波结束，第二波开始在预处理的图像上进行滤波
 z=y;
 for i=Nmax+1:e+Nmax
  for j=Nmax+1:f+Nmax
    %初始滤波窗口为3*3
    if F(i,j) == 1 %噪声点，需要进行滤波
       N =0;
       for i1= i-1:i+1
         for j1 =j-1:j+1
          if i1 == i & j1 == j%中心像素被排除在外，不进行统计
              continue;
          end
           if z(i1,j1) ~= 0 & z(i1,j1) ~= 250
              N = N + 1;%统计信号点的个数
           end
         end
       end
      if  N > 0
       z6 = zeros(1,N);%放置信号点的矩阵
         %存放信号点矩阵
         for i1 = i-1:i+1
          for j1 = j-1:j+1
            if i1 == i & j1 == j%中心像素被排除在外，不进行统计
              continue;
            end
            if z(i1,j1) ~= 0 & z(i1,j1) ~= 250 
            %if F(i1,j1) == 0 & z(i1,j1) ~= max(z1) & z(i1,j1) ~= min(z1)
              z6(1,N) = z(i1,j1);
              N = N-1;
            end
          end
         end
         y(i,j) = mean(z6); 
         %medy4 =median( [ y(i-1,j-1),y(i-1,j),y(i-1,j+1),y(i,j-1) ] );%已经处理过的4领域系数
         %y(i,j) = y(i-1,j-1)/( 1 + ( y(i-1,j-1) - medy4  )^2 ) +  y(i-1,j)/( 1 + ( y(i-1,j) - medy4  )^2 ) +  y(i-1,j+1)/( 1 + ( y(i-1,j+1) - medy4  )^2 ) + y(i,j-1)/( 1 + ( y(i,j-1) - medy4  )^2 ) ;
         %y(i,j) =( y(i-1,j-1) + y(i-1,j) + y(i-1,j+1) + y(i,j-1) ) / 4;
         %y(i,j) = median( [ y(i-1,j-1),y(i-1,j),y(i-1,j+1),y(i,j-1) ] );
         %y(i,j) = y(i-1,j-1)/( 1 + ( y(i-1,j-1) - z(i,j)  )^2 ) +  y(i-1,j)/( 1 + ( y(i-1,j) - z(i,j)  )^2 ) +  y(i-1,j+1)/( 1 + ( y(i-1,j+1) - z(i,j)  )^2 ) + y(i,j-1)/( 1 + ( y(i,j-1) - z(i,j)  )^2 ) ;
       end
       %F(i,j) = 0;%处理完了噪声点
       %z(i,j) = y(i,j);
    end
    %处理完该噪声点后，将处理完的新像素用来迭代
  end
 end
  
R(1:e,1:f) = y( Nmax+1:e+Nmax,Nmax+1:f+Nmax );%拷出滤波后的图像

主函数：
I=imread('2.png');
I=rgb2gray(I);
[e,f]=size(I);

J=imnoise(I,'salt & pepper',0.7);
%J=imnoise(J,'gaussian',0.002);
R=quickdelhighnoise(J);
tic;R=quickdelhighnoise(J);toc;
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));           %自适应中值算法去噪的峰值信噪比