一种用于抑制椒盐噪声的多窗口中值滤波器（邢藏菊等）

几种算法的比较

1）递进开关中值滤波
%commonfilt2_1.m
%递进开关中值滤波
function [R]=commonfilt2_1(x,NR0)
[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为噪声点，开始时默认所有点为信号点（毕竟噪声点很少）
%文献系数
ND1 = 3;%迭代次数
ND2 = 3;%迭代次数


%%%第一步，检测每个像素是否为噪声点
while ND1 ~= 0%要求每一次循环都影响下一次
 for i=Nmax+1:e+Nmax
  for j=Nmax+1:f+Nmax
    if ND1 == 3 %第一次检测，用5*5窗口
       NR = NR0;
       TD = 65 - 50*NR;
       z0=z(i-2:i+2,j-2:j+2);
       zmed=median(z0(:));
       if abs( z(i,j)-zmed ) > TD 
          y(i,j) = zmed;%该点为噪声
          F(i,j) = 1;%标定该点为噪声点
       end 
    else
       NR = sum(sum( F == 1 )) /(e*f);%当前窗口下噪声的密度
       if NR <= 0.25
          WD = 3;
       else
          WD = 5;
       end
       TD = 65 - 50*NR;
       z0=z(i-( (WD-1)/2 ):i+( (WD-1)/2 ),j-( (WD-1)/2 ):j+( (WD-1)/2 ));
       zmed=median(z0(:));
       if abs( z(i,j)-zmed ) > TD 
          y(i,j) = zmed;%该点为噪声
          F(i,j) = 1;%标定该点为噪声点
       end 
    end
  end
 end
 z=y;%再进行下一个大循环之前，将处理好的图像y作为下一次要处理的对象，即y赋给z
 ND1=ND1-1;
end
%此时 x即为处理后图像，下面程序的处理对象
z = y;
%%%第二步，噪声滤波
 while ND2 ~=0
  for i=Nmax+1:e+Nmax
   for j=Nmax+1:f+Nmax
     if F(i,j) == 1%若像素为噪声点
      M=sum(sum( F(i-1:i+1,j-1:j+1) == 0 ));%统计窗口内信号点的数目
      if M > 0
        z1=zeros(1,M);%放置窗口内的信号点，一共M个，M位置后面的数都为0
        for i1=i-1:i+1
         for j1=j-1:j+1
           if F(i1,j1) == 0%寻找到信号点
            z1(1,M) = z(i1,j1);
            M = M-1;
           end
         end
        end
        y(i,j)= median(z1);%信号点的中值 
      end 
      if ND2 == 1 
        F(i,j) = 0;%最后一次迭代了，噪声点标记为信号点
      end
    end
   end
  end
  z=y;
  ND2=ND2-1;
 end
R(1:e,1:f) = y( Nmax+1:e+Nmax,Nmax+1:f+Nmax );%拷出滤波后的图像

2）中值滤波
%commonfilt2_1.m
%自己编写的中值滤波
function [y]=commonfilt2_1(x,a)
y=x;
k=floor(a/2);%中心点旁边对称点的数目
[m,n]=size(x);


for i=k+1:m-k
    for j=k+1:n-k
        z=x(i-k:i+k,j-k:j+k);
        z1=z(1:end);%矩阵z转化为行向量z1
        %a*a个数冒泡排序法
        for i1=1:1:(a^2-1)
            for j1=1:1:(a^2-i1)
                if z1(1,j1) > z1(1,j1+1)
                    z0=z1(1,j1);
                    z1(1,j1)=z1(1,j1+1);
                    z1(1,j1+1)=z0;
                end
            end
        end
        y(i,j)=z1(1,k);%z1中的中值
    end
end

3）MWA算法
%commonfilt2_3.m
%MWA算法（利用阈值确定信号点）
function [y]=commonfilt2_3(x)
y=x;
[m,n]=size(x);
W=[5 3 5 3];%备选有3个窗口，其中数字为窗口的宽度
TD = 10;%文献阈值，由此判断某一点的为中心的窗口内的点是否是信号点
%k=floor( W(1,k0)/2 );
for i=3:m-2
   for j=3:n-2
      %(i,j)初始窗口都是5*5
      z=x(i-2:i+2,j-2:j+2);
      zmin=min(z(:));
      zmax=max(z(:));
      zmed=median(z(:));
      if x(i,j) == zmin | x(i,j) == zmax %灰度值极端的点，则为疑似噪声点
         M=0;%当前窗口中信号点的数目
         sum=0;%信号点的和
         k0 = 1;%问题就在这，怪不得每次调程序没调好：每次对一个新的像素，窗口都要从W(1,1)开始
        while M == 0
         k0 = k0+1;%窗口换做下一个,此时窗口大小为W(1,k0)
         k=floor( W(1,k0)/2 );
         %判断新窗口内是否有信号点
         for i1=i-k:i+k
          for j1=j-k:j+k
           %判断x(i1,j1)是不是信号点：
           %判断x（i1,j1）是否为信号点
            if abs( x(i1,j1) - zmed ) < TD
              M = M+1;
              sum = sum+x(i1,j1);
            end
          end
         end
        %循环结束后，M的个数也确定了
       end
       if M ~= 0
         y(i,j) = double(sum)/M;%信号点的均值赋给中心点
       end
    end
  end
end

4）主函数
I=imread('1.png');
I=rgb2gray(I);
NR0 = 0.6;%噪声密度
J=imnoise(I,'salt & pepper',NR0);

k1=commonfilt2_1(J,NR0);%递进开关中值滤波方法（PSM算法）
tic;k1=commonfilt2_1(J,NR0);toc;
%k1=commonfilt2_1_1(J);%NMA算法（跟最大值最小值比较，确定信号点）
%k2=commonfilt2_2(J,5);%库函数中值滤波方法
%k3=commonfilt2_3(J);%MWA算法（利用阈值，确定信号点）

%subplot(221),imshow( uint8(J) ),title('噪声图像');
%subplot(222),imshow( uint8(k1) ),title('MWA算法中值滤波(最大最小值确定信号点)');
subplot(222),imshow( uint8(k1) ),title('递进开关中值滤波');
%subplot(223),imshow( uint8(k2)),title('库函数中值滤波');
%subplot(224),imshow( uint8(k3)),title('MWA算法中值滤波（阈值确定信号点）');


% 计算三种算法的峰值信噪比
[e f]=size(I);
B=8;                %编码一个像素用多少二进制位
MAX=2^B-1;          %图像有多少灰度级


I=double(I);
k1=double(k1);
%k2=double(k2);
%k3=double(k3);


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


%%%%%% psnr3=
%MES3=sum(sum((I-k3).^2))/(e*f);    
%PSNR3=20*log10(MAX/sqrt(MES3));

特别的是PSMF算法为（转载）：

%% Implementation of Progressive Switching Median Filter
%% Base Paper : Zhou Wang and David Zhang, "Progressive Switching Median
%% Filter for the Removal of Impulse Noise from Highly Corrupted Images", 
%% IEEE Trans. on Cir. and Sys., vol. 46, no. 1, Jan. 1999.
%% Function Y = PSMF(x)
%% input    x = Image is corrupted by Salt & Pepper Noise
%%          
%%  Example: Y = PSMF(x);
%%      Posted date   : 16 - 10 - 2008
%%      Modified date : 
%%                  
%% Developed By : K.Kannan (kannan.keizer@gmail.com) 
%%                & Jeny Rajan (jenyrajan@gmail.com)
%%                  Medical Imaging Research Group (MIRG), NeST,
%%                  Trivandrum.
%% Progressive Switching Median Filter
function Y = PSMF(x)
x = double(x);
WF = 3; ND = 3;T = 40;a = 65;b = -50;
M = medfilt2(x,[3 3]);
N = abs(x - M);
N(N>T)=0;
N = N ~= 0;
N = double(N);
R = sum(N(:))/(size(x,1) * size(x,2));
if R <= 0.25
    WD = 3;
else
    WD = 5;
end
TD = a + (b * R);
z = IMPDET(x,ND,WD,TD);
Y = NF(x,z,WF);

%% Impulse Detection
function F1 = IMPDET(x,ND,WD,TD)
X = x;
M = medfilt2(X,[WD WD]);
D = abs(X - M);
F = zeros(size(x));
F(D>=TD)=1;
F1 = F;
X(F1==F)=X(F1==F);
X(F1~=F)=M(F1~=F);
for i = 1:ND-1
    M = medfilt2(X,[WD WD]);
    F1(abs(X - M)<TD)=F(abs(X - M)<TD);
    F1((X - M)>=TD)=1;
    X(F1==F)=X(F1==F);
    X(F1~=F)=M(F1~=F);
    F = F1;
end
return;

%% Noise Filtering
function Y = NF(x,f,WF)
g = f;
Y = x;
Y1 = Y;
g1 = g;
s = sum(g(:));
while s ~= 0
    M = medfilt2(Y,[WF WF]);
    Y1(g==1)=M(g==1);
    g1(Y~=Y1)=0;
    Y = Y1;
    g = g1;
    s1 = sum(g(:));
    if s1 ~= s
        s = s1;
    else
        s = 0;
    end
end
return;