用于图像处理的自适应中值滤波（张旭明等）

基于文章“用于图像处理的自适应中值滤波”的matlab代码如下：

%commonfilt2_1.m
%一种自适应调整窗口，一种自适应滤波的算法
function [y]=commonfilt2_1(x)

TD = 9;%判断噪声点所用阈值
N1 = 0; %子块1中噪声点的个数
N2 = 0; %子块2中噪声点的个数
N3 = 0; %子块3中噪声点的个数
N4 = 0; %子块4中噪声点的个数
W1 = 0; %子块1需要滤波的窗口的大小
W2 = 0; %子块2需要滤波的窗口的大小
W3 = 0; %子块3需要滤波的窗口的大小
W4 = 0; %子块4需要滤波的窗口的大小

y=x;
[m,n]=size(x);
F=zeros(m,n);%检查噪声点的矩阵。元素值为1则为噪声点，信号点则为0.初值设为0

%%%%%%%%%%%%%%第一步，检测子块中的噪声，按照文献，我们分图像为4块
%第一步，检测子块1中的噪声点：使用3*3窗口检测
for i=2:255
    for j=2:255
        sum1 = 0; %小窗口内非极值的像素之和
        M1 = 0; %小窗口内非极值的像素个数
        z1=x(i-1:i+1,j-1:j+1);%3*3窗口进行检测
        z1max=max(z1(:));
        z1min=min(z1(:));
       %找出小窗口内与最值不相等的点的均值
       for i1=i-1:i+1
         for j1=j-1:j+1
           if x(i1,j1) ~= z1max & x(i1,j1) ~= z1min
               sum1 = sum1+x(i1,j1);
               M1 = M1+1;
           end
         end
       end
       T1 = sum1/M1;%窗口内非极值的像素的均值
       %下面开始判断噪声点
      if( x(i,j) == z1max | x(i,j) == z1min ) & ( abs( x(i,j) - T1) >= TD )%满足该条件的像素是噪声点
           F(i,j) = 1; %标记相应的点为噪声点
           N1 = N1+1;
      end
    end
end

%第一步，检测子块2中的噪声点：使用3*3窗口检测
for i=2:255
    for j=258:511
        sum2 = 0; %小窗口内非极值的像素之和
        M2 = 0; %小窗口内非极值的像素个数
        z2=x(i-1:i+1,j-1:j+1);%3*3窗口进行检测
        z2max=max(z2(:));
        z2min=min(z2(:));
       %找出小窗口内与最值不相等的点的均值
       for i1=i-1:i+1
         for j1=j-1:j+1
           if x(i1,j1) ~= z2max & x(i1,j1) ~= z2min %非最值
               sum2 = sum2+x(i1,j1);
               M2 = M2+1;
           end
         end
       end
       T2 = sum2/M2;%窗口内非极值的像素的均值
       %下面开始判断噪声点
      if( x(i,j) == z2max | x(i,j) == z2min ) & ( abs( x(i,j) - T2) >= TD )%满足该条件的像素是噪声点
           F(i,j) = 1; %标记相应的点为噪声点
           N2 = N2+1;
      end
    end
end

%第一步，检测子块3中的噪声点：使用3*3窗口检测
for i=258:511
    for j=2:255
        sum3 = 0; %小窗口内非极值的像素之和
        M3 = 0; %小窗口内非极值的像素个数
        z3=x(i-1:i+1,j-1:j+1);%3*3窗口进行检测
        z3max=max(z3(:));
        z3min=min(z3(:));
       %找出小窗口内与最值不相等的点的均值
       for i1=i-1:i+1
         for j1=j-1:j+1
           if x(i1,j1) ~= z3max & x(i1,j1) ~= z3min
               sum3 = sum3+x(i1,j1);
               M3 = M3+1;
           end
         end
       end
       T3 = sum3/M3;%窗口内非极值的像素的均值
       %下面开始判断噪声点
      if( x(i,j) == z3max | x(i,j) == z3min ) & ( abs( x(i,j) - T3) >= TD )%满足该条件的像素是噪声点
           F(i,j) = 1; %标记相应的点为噪声点
           N3 = N3+1;
      end
    end
end

%第一步，检测子块4中的噪声点：使用3*3窗口检测
for i=258:511
    for j=258:511
        sum4 = 0; %小窗口内非极值的像素之和
        M4 = 0; %小窗口内非极值的像素个数
        z4=x(i-1:i+1,j-1:j+1);%3*3窗口进行检测
        z4max=max(z4(:));
        z4min=min(z4(:));
       %找出小窗口内与最值不相等的点的均值
       for i1=i-1:i+1
         for j1=j-1:j+1
           if x(i1,j1) ~= z4max & x(i1,j1) ~= z4min
               sum4 = sum4+x(i1,j1);
               M4 = M4+1;
           end
         end
       end
       T4 = sum4/M4;%窗口内非极值的像素的均值
       %下面开始判断噪声点
      if( x(i,j) == z4max | x(i,j) == z4min ) & ( abs( x(i,j) - T4) >= TD )%满足该条件的像素是噪声点
           F(i,j) = 1; %标记相应的点为噪声点
           N4 = N4+1;
      end
    end
end

%%%%%%%%%第二步，确定各个子块中滤波器窗口的大小
%%每一子块的大小都是相同的，都是256*256
W=[W1,W2,W3,W4];%存放窗口大小的矩阵，初试值为0
p1 = double(N1/(256*256));
p2 = double(N2/(256*256));
p3 = double(N3/(256*256));
p4 = double(N4/(256*256));
p=[p1,p2,p3,p4];
for i =1:4
    if p(1,i) >= 0.002 & p(1,i) < 0.25
      W(1,i) = 3;
    end
    if p(1,i) >= 0.25 & p(1,i) < 0.45
      W(1,i) = 5;
    end
    if p(1,i) >= 0.45
      W(1,i) = 7;
    end
end

%%%%%%%%%%%第三步，噪声滤波
%对第一个子块滤波
k1=floor(W1/2)+1;
for i=k1:257-k1
  for j=k1:257-k1
    if F(i,j) == 1 %若是噪声点
      z1 = x(i-k1+1:i+k1-1,j);
      z11 = median(z1(:));%0°方向上像素的中值
      
      z2 = x(i,j-k1+1:j+k1-1);
      z22 = median(z2(:));%90°方向上像素的中值
      
      %处理45°和135°要麻烦一点
      k=1;z3=zeros(1,2*k1-1);
      for l = 1-k1:k1-1
        z3(1,k) = x(i+l,j+l);
        k = k+1;
      end
      z33 =median(z3(:));%45°方向上像素的中值
      
      k=1;z4=zeros(1,2*k1-1);
      for l = 1-k1:k1-1
        z4(1,k) = x(i+l,j-l);
        k = k+1;
      end
      z44 =median(z4(:));%135°方向上像素的中值
      %根据文献，求出4个方向的中值，再乘以相关的系数累加求和赋给中心点
      y(i,j) = (z11^2 + z22^2 + z33^2 + z44^2)/(z11+z22+z33+z44);
    end
  end
end


%对第二个子块滤波
k2=floor(W2/2)+1;
for i=k2:257-k2
  for j=256+k2:513-k2
    if F(i,j) == 1
      z1 = x(i-k2+1:i+k2-1,j);
      z11 = median(z1(:));%0°方向上像素的中值
      
      z2 = x(i,j-k2+1:j+k2-1);
      z22 = median(z2(:));%90°方向上像素的中值
      
      %处理45°和135°要麻烦一点
      k=1;z3=zeros(1,2*k2-1);
      for l = 1-k2:k2-1
        z3(1,k) = x(i+l,j+l);
        k = k+1;
      end
      z33 =median(z3(:));%45°方向上像素的中值
      
      k=1;z4=zeros(1,2*k2-1);
      for l = 1-k2:k2-1
        z4(1,k) = x(i+l,j-l);
        k = k+1;
      end
      z44 =median(z4(:));%135°方向上像素的中值
      %根据文献，求出4个方向的中值，再乘以相关的系数累加求和赋给中心点
      y(i,j) = (z11^2 + z22^2 + z33^2 + z44^2)/(z11+z22+z33+z44);
    end
  end
end

%对第三个子块滤波
k3=floor(W3/2)+1;
for i=256+k3:513-k3
  for j=k3:257-k3
    if F(i,j) == 1
      z1 = x(i-k3+1:i+k3-1,j);
      z11 = median(z1(:));%0°方向上像素的中值
      
      z2 = x(i,j-k3+1:j+k3-1);
      z22 = median(z2(:));%90°方向上像素的中值
      
      %处理45°和135°要麻烦一点
      k=1;z3=zeros(1,2*k3-1);
      for l = 1-k3:k3-1
        z3(1,k) = x(i+l,j+l);
        k = k+1;
      end
      z33 =median(z3(:));%45°方向上像素的中值
      
      k=1;z4=zeros(1,2*k3-1);
      for l = 1-k3:k3-1
        z4(1,k) = x(i+l,j-l);
        k = k+1;
      end
      z44 =median(z4(:));%135°方向上像素的中值
      %根据文献，求出4个方向的中值，再乘以相关的系数累加求和赋给中心点
      y(i,j) = (z11^2 + z22^2 + z33^2 + z44^2)/(z11+z22+z33+z44);
    end
  end
end

%对第四个子块滤波
k4=floor(W4/2)+1;
for i=256+k4:513-k4
  for j=256+k4:513-k4
    if F(i,j) == 1
      z1 = x(i-k4+1:i+k4-1,j);
      z11 = median(z1(:));%0°方向上像素的中值
      
      z2 = x(i,j-k4+1:j+k4-1);
      z22 = median(z2(:));%90°方向上像素的中值
      
      %处理45°和135°要麻烦一点
      k=1;z3=zeros(1,2*k4-1);
      for l = 1-k4:k4-1
        z3(1,k) = x(i+l,j+l);
        k = k+1;
      end
      z33 =median(z3(:));%45°方向上像素的中值
      
      k=1;z4=zeros(1,2*k4-1);
      for l = 1-k4:k4-1
        z4(1,k) = x(i+l,j-l);
        k = k+1;
      end
      z44 =median(z4(:));%135°方向上像素的中值
      %根据文献，求出4个方向的中值，再乘以相关的系数累加求和赋给中心点
      y(i,j) = (z11^2 + z22^2 + z33^2 + z44^2)/(z11+z22+z33+z44);
    end
  end
end

%commonfilt2_2.m
%库函数中值滤波
function [y]=commonfilt2_2(x,a,b)
y=x;
[m,n]=size(x);
k=floor(a/2);
for i=k+1:m-k
    for j=k+1:n-k
        z=x(i-k:i+k,j-k:j+k);
        y(i,j)=median(z(:));
    end
end

I=imread('1.png');
I=rgb2gray(I);
J=imnoise(I,'salt & pepper',0.4);

k1=commonfilt2_1(J);%一种自适应调整窗口，一种自适应滤波的算法
k2=commonfilt2_2(J,3,3);%库函数中值滤波方法
k3=commonfilt2_2(J,5,5);%库函数中值滤波方法

subplot(222),imshow( uint8(k2)),title('库函数中值滤波3*3');
subplot(224),imshow( uint8(k1) ),title('自适应调整窗口，自适应滤波');
subplot(223),imshow( uint8(k3)),title('库函数中值滤波5*5');
subplot(221),imshow( uint8(J)),title('噪声图像');

% 计算三种算法的峰值信噪比
[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));