一种改进的中值滤波算法（张恒等）

主函数：

%1 三种算法比较
I=imread('11.png');
I=rgb2gray(I);
%加入混合噪声
%var=100/(255^2);%imnoise中的方差等于实际方差除以灰度级的平方
%按照文献设置参数
J=imnoise(I,'salt & pepper',0.8);
%J=imnoise(J,'gaussian',0,var);%var=0.0015


%k1=commonfilt2_1(J,3);%普通的中值滤波
%k2=commonfilt2_2(J,3);%普通的均值滤波
%k3=commonfilt2_3(J,3);%MTM算法
k4=commonfilt2_4(J,3);%改进的中值滤波IMF


%subplot(221),imshow( uint8(k1)),title('普通的中值滤波');
%subplot(222),imshow( uint8(k2)),title('普通的均值滤波');
%subplot(223),imshow( uint8(k3)),title('MTM算法滤波');
%subplot(224),
imshow( uint8(k4)),title('改进的中值滤波算法IMF');

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

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

%%%%% psnr4=
MES4=sum(sum((I-k4).^2))/(e*f);    
PSNR4=20*log10(MAX/sqrt(MES4));

%S3=PIF(I,J,k3);%信噪比改善因子
%S4=PIF(I,J,k4);%信噪比改善因子

子函数：

%commonfilt2_1.m
%中值滤波,这里我未添加边界，所以相比于算法3（p =0）效果略差
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);
        y(i,j)=median(z(:));
    end
end

%commonfilt2_2.m
%均值滤波
function [y]=commonfilt2_2(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);
        y(i,j)=mean(z(:));
    end
end

%commonfilt2_3.m
%处理（i，j）这个像素时，确定一个3*3窗口，再确定一个区间：[z0-p,z0+p]，对落在其中的点求其均值，取代中心点
function [R]=commonfilt2_3(x,a) %p为阈值
[e,f]=size(x);
Nmax=1; %确定最大的滤波半径，那么滤波窗口最大为3
%下面是边界扩展，图像上下左右各增加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为最后输出的图像

p = 3;%参数
k = floor(a/2);

for i=Nmax+1:e+Nmax
  for j=Nmax+1:f+Nmax
     zsum=0;
     num=0;
     z0=z(i-k:i+k,j-k:j+k);%确定一个窗口
     zmedian=median(z0(:));%小窗口内的中值
     for i1=i-k:i+k
        for j1=j-k:j+k
          if ( z(i1,j1) >= (zmedian-p)  ) & ( z(i1,j1) <= (zmedian + p)  )
             zsum = zsum + z (i1,j1);%满足该条件，则将该灰度值相加，最后求取其均值
             num = num+1;
          end 
        end
     end
     y(i,j)=zsum/num;
   end
end


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

%commonfilt2_4.m
%改进的中值滤波（IMF）
function [R]=commonfilt2_4(x,a)
[e,f]=size(x);
k=floor(a/2);
Nmax=1; %确定最大的滤波半径，那么滤波窗口最大为3
%下面是边界扩展，图像上下左右各增加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);%标记矩阵

for i=Nmax+1:e+Nmax
  for j=Nmax+1:f+Nmax
    %检测窗口开始用3*3
    %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


for i=Nmax+1:e+Nmax
  for j=Nmax+1:f+Nmax
      if F(i,j) == 1
        z0=z(i-k:i+k,j-k:j+k);
        zmedian=median(z0(:));%3*3窗口内的中值
        sum=0;r = 0;%加权系数乘以系数之和 以及 加权系数之和
        for i1=i-k:i+k
            for j1=j-k:j+k
              sum=sum + z(i1,j1) / ( 1+( z(i1,j1) - zmedian )^2 );
              r = r + 1.0 / ( 1+( z(i1,j1) - zmedian )^2 );
            end
        end
        y(i,j) = sum / r;%加权系数的均值
      end
    end
end


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