Matlab图像处理之Refined Lee滤波器 [精致Lee滤波器]

一、介绍

  Refined Lee滤波器是在Lee滤波器的基础上，改进的滤波算法，与Lee滤波器最主要的区别是采用了非正方形的滤波窗口，使其在对图像降噪的同时，更好地保留图像的细节信息，因此在雷达和卫星图像处理中得到广泛应用。

二、原理

  Refined Lee滤波器的算法实现主要包括两步：模板窗口选择和局部统计滤 波，具体步骤为：
  （1）模板窗口选择
  将滤波窗口分解为3x3的9个子窗格， 计算各子窗格内的均值，然后采用如下图所示 4个边缘梯度算子进行检测，每个边缘梯度算子对应图中每列的2个边缘方向窗口 。选择梯度最大的算子所对应的两个边缘方向窗口中均值与中心像素均值比较接近的那个。

  （2）局部统计滤波
  在边缘方向窗口中，利用滤波器进行滤波。

三、matlab代码

clear;clc;close all;

fid=fopen("C:\Users\Administrator\Desktop\3.tif",'r');
A=fread(fid,[3080 2816],'float32');
fclose(fid);
A=A';
complexA=complex(A(:,1:2:end),A(:,2:2:end));
I=abs(complexA);               %读取图像
I0=log(I)*40;

[m,n]=size(I);
figure;
imshow(uint8(I0));

 VarV=1;
 u(1:7,1:7)=0;
 ver=[-1 0 1;-1 0 1;-1 0 1];         %4个梯度模版
 deg_135=[0 1 1;-1 0 1;-1 -1 0];
 leve=[1 1 1;0 0 0;-1 -1 -1];
 deg_45=[1 1 0;1 0 -1;0 -1 -1];                 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 for i=1:m-6  
     for k=1:n-6   
        u=I(i:i+6,k:k+6);   %滑动窗口7*7
        h=1;
        for t=1:2:5          %子窗口3*3
            
            for f=1:2:5
                u_1=u(t:t+2,f:f+2);
                u_1_mean=sum(sum(u_1))/9;       %子窗口均值
                mean(h)=u_1_mean;
                h=h+1;
            end
        end
        
        me=[mean(1:3);mean(4:6);mean(7:9)];    %3*3子窗口均值矩阵
        
        mo1=ver.*me;                      %利用梯度模板判别邻域窗口
        mo1_1=abs(sum(sum(mo1)));  
        mo2=deg_135.*me;
        mo2_1=abs(sum(sum(mo2)));  
        mo3=leve.*me;
        mo3_1=abs(sum(sum(mo3)));
        mo4=deg_45.*me;
        mo4_1=abs(sum(sum(mo4)));
        
        max_mo=max([mo1_1 mo2_1 mo3_1 mo4_1]);   %判断最大值
        
        switch max_mo  
            
            case  mo1_1           %垂直情况
                  model=ver;
                  if  abs(me(2,3)-me(2,2))>abs(me(2,1)-me(2,2)) 
                      u(1:7,5:7)=0;      %%%%%%%%%%%%
                  else abs(me(2,3)-me(2,2))<=abs(me(2,1)-me(2,2))
                       u(1:7,1:3)=0;     %%%%%%%%%%%%%
                  end
                   u_mean=sum(sum(u))/49;  
                   Vary=std2(u)^2;
                   Varx=(Vary-u_mean^2)/(1+VarV);
                   if Varx<0
                       Varx=0;
                   end
                   b=Varx/Vary;
                   X(i,k)=u_mean+b*(u(4,4)-u_mean);
                   
            case mo2_1
                 model=deg_135;     %135度情况
                 if abs(me(2,3)-me(2,2))>abs(me(2,1)-me(2,2))
                    L2_1=1;
                    for L2=2:7
                    u(L2_1,L2:7)=0;
                    L2_1=L2_1+1;
                    end
                else abs(me(2,3)-me(2,2))<abs(me(2,1)-me(2,2))
                     L2_2=1;
                     for L2=2:7
                     u(L2:7,L2_2)=0;        
                     L2_2=L2_2+1;
                     end
                 end
                    u_mean=sum(sum(u))/49;
                    Vary=std2(u)^2;
                    Varx=(Vary-u_mean^2)/(1+VarV);
                    if Varx<0
                       Varx=0;
                   end
                    b=Varx/Vary;
                    X(i,k)=u_mean+b*(u(4,4)-u_mean);
                
            case mo3_1
                 model=leve;               %水平情况
                 if abs(me(3,3)-me(2,2))>abs(me(1,1)-me(2,2))
                    u(5:7,1:7)=0;
                 else abs(me(3,3)-me(2,2))<=abs(me(1,1)-me(2,2))
                   u(1:3,1:7)=0;
                 end
                 u_mean=sum(sum(u))/49;
                 Vary=std2(u)^2;
                 Varx=(Vary-u_mean^2)/(1+VarV);
                 if Varx<0
                       Varx=0;
                   end
                 b=Varx/Vary;
                 X(i,k)=u_mean+b*(u(4,4)-u_mean);
               
            case mo4_1
                 model=deg_45;        %45度情况
                 if abs(me(3,3)-me(2,2))>abs(me(1,1)-me(2,2))
                    L4_1=2;
                    for q=7:-1:2
                    u(L4_1,q:7)=0;
                    L4_1=L4_1+1;
                    end
                else abs(me(3,3)-me(2,2))<=abs(me(1,1)-me(2,2))
                     L4_1=1;
                     for  q=6:-1:1
                     u(L4_1,1:q)=0;
                     L4_1=L4_1+1; 
                     end
                 end
                 u_mean=sum(sum(u))/49;
                 Vary=std2(u)^2;
                 Varx=(Vary-u_mean^2)/(1+VarV);
                 if Varx<0
                       Varx=0;
                   end
                 b=Varx/Vary;
                 X(i,k)=u_mean+b*(u(4,4)-u_mean);                                           
        end   
     end
 end
 
  im=X/max(max(X));            %拉伸作图
  im=imadjust(im,stretchlim(im),[],1);
  figure;
  imshow(im);