有效去除图像中脉冲噪声的新型滤波算法（曲延锋等）

%imfledfilt.m
%有效去除图像中脉冲噪声的新型滤波算法
%检测窗口5*5，滤波窗口3*3
function [R] = imfledfilt(x,Noiserate)
[e,f]=size(x);
Nmax=2; %确定最大的滤波半径，那么滤波窗口或者检测窗口最大为5
%下面是边界扩展，图像上下左右各增加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为噪声点，开始时默认所有点为信号点（毕竟噪声点很少）

%%%第一步，检测每个像素是否为噪声点
for i=Nmax+1:e+Nmax
  for j=Nmax+1:f+Nmax
    %检测窗口开始用5*5
    z0=z(i-2:i+2,j-2:j+2);
    z1=z0(1:end);
    %if  z(i,j) == max(z1) | z(i,j) == min(z1)%疑似噪声点，需要进一步判断
    if  z(i,j) == 0 | z(i,j) == 255 %疑似噪声点，需要进一步判断
      F(i,j) = 1;
    end
 end
end

%由文献条件
if Noiserate < 0.6
  Wf = 3;
else
  Wf = 5;
end
%%%第二步，噪声滤波
%while F ~= zeros(e+2*Nmax,f+2*Nmax)%当去噪未彻底，那需要反复滤波去噪
%H = 3;
%while H ~= 0
  for i=Nmax+1:e+Nmax
    for j=Nmax+1:f+Nmax
      if F(i,j) == 1%若像素为噪声点
         M = sum(sum( F(i-( (Wf-1)/2 ):i+( (Wf-1)/2 ),j-( (Wf-1)/2 ):j+( (Wf-1)/2 )) == 0 ));%窗口内信号点的个数
         if M >0
           z1 =zeros(1,M);
           for i1=i-( (Wf-1)/2 ):i+( (Wf-1)/2 )
            for j1=j-( (Wf-1)/2 ):j+( (Wf-1)/2 )
              if F(i1,j1) == 0%寻找到信号点
                z1(1,M) = y(i1,j1);
                M = M -1;
              end
            end
           end
          y(i,j) = median(z1);%求取信号点的中值
          F(i,j) = 0;%噪声点处理后标记为信号点
         end
      end 
    end
  end
%  z = y;
% H = H -1;
%end

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

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

Noiserate = 0.9;
J=imnoise(I,'salt & pepper',Noiserate);
R=imfledfilt(J,Noiserate);

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=
MES=sum(sum((I-R).^2))/(e*f);     %自适应中值去噪的均方差
PSNR=20*log10(MAX/sqrt(MES));           %自适应中值算法去噪的峰值信噪比

MES1 = sum(sum((I-R).^2));
MES2 = sum(sum(I.^2));
SNR = 10*log10(MES2/MES1);
tic;R=imfledfilt(J,Noiserate );toc;