1 简介

D.Gabor 1946年提出

窗口Fourier变换，为了由信号的Fourier变换提取局部信息，引入了时间局部化的窗函数。

由于窗口Fourier变换只依赖于部分时间的信号，所以，现在窗口Fourier变换又称为短时Fourier变换，这个变换又称为Gabor变换。

1) Gabor优点

Gabor小波与人类视觉系统中简单细胞的视觉刺激响应非常相似。它在提取目标的局部空间和频率域信息方面具有良好的特性。虽然Gabor小波本身并不能构成正交基，但在特定参数下可构成紧框架。Gabor小波对于图像的边缘敏感，能够提供良好的方向选择和尺度选择特性，而且对于光照变化不敏感,能够提供对光照变化良好的适应性。上述特点使Gabor小波被广泛应用于视觉信息理解。

Gabor滤波器和脊椎动物视觉皮层感受野响应的比较：第一行代表脊椎动物的视觉皮层感受野，第二行是Gabor滤波器，第三行是两者的残差。可见两者相差极小。Gabor滤波器的这一性质，使得其在视觉领域中经常被用来作图像的预处理。

2) Gabor定义

① 具体窗函数――Gaussaion的 Gabor变换定义式

Gabor变换的基本思想：把信号划分成许多小的时间间隔，用傅里叶变换分析每一个时间间隔，以便确定信号在该时间间隔存在的频率。其处理方法是对f(t)加一个滑动窗，再作傅里叶变换。

② 窗口的宽高关系

经理论推导可以得出：高斯窗函数条件下的窗口宽度与高度，且积为一固定值。

3) 离散Gabor变换的一般求法

① 首先选取核函数

可根据实际需要选取适当的核函数。如，如高斯窗函数；

② 离散Gabor变换的表达式

4) Gabor变换的解析理论

对偶函数可以使计算更为简洁方便。

5) 适用条件

① 临界采样Gabor展开要求条件：TΩ=2π；

② 过采样展开要求条件：TΩ≤2π；

当TΩ＞2π时，欠采样Gabor展开，已证明会导致数值上的不稳定。

6) 应用

① 暂态信号检测

如果对信号波形有一定的先验知识且可以据此选取合适的基函数，可以用Gabor变换对信号作精确的检测统计计量。

② 图象分析与压缩

二维Gabor变换可以应用到图象分析与压缩中。

3. 二维Gabor滤波器

用Gabor 函数形成的二维Gabor 滤波器具有在空间域和频率域同时取得最优局部化的特性，因此能够很好地描述对应于空间频率(尺度)、空间位置及方向选择性的局部结构信息。

Gabor滤波器的频率和方向表示接近人类视觉系统对于频率和方向的表示，并且它们常备用于纹理表示和描述。

在图像处理领域，Gabor滤波器是一个用于边缘检测的线性滤波器。

在空域，一个2维的Gabor滤波器是一个正弦平面波和高斯核函数的乘积。

Gabor滤波器是自相似的，也就是说，所有Gabor滤波器都可以从一个母小波经过膨胀和旋转产生。

实际应用中，Gabor滤波器可以在频域的不同尺度，不同方向上提取相关特征。

 

Gabor滤波器的傅里叶变换：峰值响应在复正弦的空域频率(u0,v0)：

Gabor滤波器示意图，3种角度5种方向：

​2 部分代码

b1=[0 1 0

    1 1 1

    0 1 0];

b2=[0 0 1 0 0

    0 1 1 1 0

    1 1 1 1 1

    0 1 1 1 0

    0 0 1 0 0];

b31=[0 0 0

     1 1 1

     0 0 0];

b32=[0 0 1

     0 1 0

     1 0 0];

b33=[0 1 0

     0 1 0

     0 1 0];

b34=[1 0 0

     0 1 0

     0 0 1];

I = imread('C:\Documents and Settings\Administrator\桌面\matlab code\邵海峰\matlab代码\原始指纹图像2.jpg');

I=rgb2gray(I);

I=double(I);

figure,imshow(I,[])

xs1=size(I,1);

xs2=size(I,2);

K1=sum(sum(I));

M1=K1/(xs1*xs2);

M2=ones(xs1,xs2);

M0=M2*M1;

M3=zeros(xs1,xs2);

for i=1:xs1

    for j=1:xs2

        M3(i,j)=(I(i,j)-M0(i,j))^2;   %求一个平均值 

    end

end

K2=sum(sum(M3));

V0=K2/(xs1*xs2);

P1=ones(fix(xs1/15),fix(xs2/15));

P2=ones(15,15);

P3=ones(fix(xs1/15),fix(xs2/15));

for i=1:xs1

    for j=1:xs2

        a=mod(i,15);

        b=mod(j,15);

        if a==0

            a=15;

        end

        if b==0

            b=15;

        end

        P2(a,b)=I(i,j);

        if a==15&&b==15

            PK1=sum(sum(P2));

            PM1=PK1/(15*15);

            PM2=ones(15,15);

            PM0=PM2*PM1;

            PM3=zeros(15,15);

            for m=1:15

                for n=1:15

                    PM3(m,n)=(PM2(m,n)-PM0(m,n))^2;

                end

            end

            PK2=sum(sum(PM3));

            PV0=PK2/(15*15);

            P1(fix(i/15),fix(j/15))=PM1;

            P3(fix(i/15),fix(j/15))=PV0;

        end        

    end

end

I_n=zeros(xs1,xs2);

for i=1:xs1

    for j=1:xs2

        m=ceil(i/15);

        n=ceil(j/15);

        if I(i,j)>M0(i,j)&&m~=ceil(xs1/15)&&n~=ceil(xs2/15)

            I_n(i,j)=M0(i,j)+sqrt(V0*(I(i,j)-P1(m,n))^2/P3(m,n));

        end

        if I(i,j)<=M0(i,j)&&m~=ceil(xs1/15)&&n~=ceil(xs2/15)

            I_n(i,j)=M0(i,j)-sqrt(V0*(P1(m,n)-I(i,j))^2/P3(m,n));

        end

    end

end

bw1= I;

vertical=[ 1  2  1

           0  0  0

          -1  -2  -1];

horizontal=[ -1  0  1

             -2  0  2

             -1  0  1];

fx=conv2(bw1,vertical,'same');

fy=conv2(bw1,horizontal,'same');

mubanver=15;mubanhor=15;  %模板的垂直和水平大小

middle_mu1=size(bw1,1)/mubanver;

middle_mu2=size(bw1,2)/mubanhor;

middle_mu1=fix(middle_mu1);

middle_mu2=fix(middle_mu2);

thea= zeros(middle_mu1, middle_mu2);

thea1= zeros(middle_mu1, middle_mu2);

sum1= zeros(middle_mu1, middle_mu2);

sum2= zeros(middle_mu1, middle_mu2);

sum3= zeros(middle_mu1, middle_mu2);

for i=1:middle_mu1

    for j=1:middle_mu2

        for y=1:mubanver

            for x=1:mubanhor

                sum1(i,j)=sum1(i,j)+2*fx(y+(i-1)*(mubanver),x+(j-1)*(mubanhor))*fy(y+(i-1)*(mubanver),x+(j-1)*(mubanhor));

                %2Gxy

                sum2(i,j)=sum2(i,j)+fx(y+(i-1)*(mubanver),x+(j-1)*(mubanhor))*fx(y+(i-1)*(mubanver),x+(j-1)*(mubanhor)) - fy(y+(i-1)*(mubanver),x+(j-1)*(mubanhor))*fy(y+(i-1)*(mubanver),x+(j-1)*(mubanhor));

                %Gxx-Gyy

                sum3(i,j)=sum3(i,j)+fx(y+(i-1)*(mubanver),x+(j-1)*(mubanhor))*fx(y+(i-1)*(mubanver),x+(j-1)*(mubanhor)) + fy(y+(i-1)*(mubanver),x+(j-1)*(mubanhor))*fy(y+(i-1)*(mubanver),x+(j-1)*(mubanhor));

                %Gxx+Gyy

            end

        end

        if sum2(i,j)>=0

           fi=1/2*atan(sum1(i,j)/sum2(i,j));

        end

        if sum2(i,j)<0 && sum1(i,j)>=0

            fi=1/2*(atan(sum1(i,j)/sum2(i,j))+pi);

        end

        if sum2(i,j)<0 && sum1(i,j)<0

            fi=1/2*(atan(sum1(i,j)/sum2(i,j))-pi);

        end

           thea(i,j)=pi-fi;   %方向角度 

    end

end

for i=1:size(sum3,1)

   for j=1:size(sum3,2)

       if sum3(i,j)==0

           coh(i,j)=1;

       else

           coh(i,j)=sqrt(sum2(i,j)*sum2(i,j)+sum1(i,j)*sum1(i,j))/sum3(i,j);

       end

   end

end

figure,imshow(coh,[]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%根据置信度改变方向%%%%%%%%%%%%%%%%%%%%%%%%

[X,Y] = meshgrid(size(thea,1):-1:1,1:1:size(thea,2));

xx=X';yy=Y';

for i=1:size(thea,1)

   for j=1:size(thea,2)

       sum2(i,j)=cos(thea(i,j));

       sum1(i,j)=sin(thea(i,j));

    end     

end

figure,quiver(yy,xx,sum2,sum1);

    xs1=size(I,1);

    xs2=size(I,2);

    thea_jz = thea;

    I=double(I);

    bw1= I;

mubanver=15;mubanhor=15;  %模板的垂直和水平大小

middle_mu1=size(bw1,1)/mubanver;

middle_mu2=size(bw1,2)/mubanhor;

middle_mu1=fix(middle_mu1);

middle_mu2=fix(middle_mu2);

for i=1:middle_mu1

       for j=1:middle_mu2

           if thea_jz(i,j)>=7/16*pi&&thea_jz(i,j)<9/16*pi

               thea_jz(i,j)=1/2*pi;

           end

           if thea_jz(i,j)>=9/16*pi&&thea_jz(i,j)<11/16*pi

               thea_jz(i,j)=5/8*pi;

           end

           if thea_jz(i,j)>=11/16*pi&&thea_jz(i,j)<13/16*pi

               thea_jz(i,j)=3/4*pi;

           end

           if thea_jz(i,j)>=13/16*pi&&thea_jz(i,j)<15/16*pi

               thea_jz(i,j)=7/8*pi;

           end

           if thea_jz(i,j)>=15/16*pi&&thea_jz(i,j)<17/16*pi

               thea_jz(i,j)=pi;

           end

           if thea_jz(i,j)>=17/16*pi&&thea_jz(i,j)<19/16*pi

               thea_jz(i,j)=9/8*pi;

           end

           if thea_jz(i,j)>=19/16*pi&&thea_jz(i,j)<21/16*pi

              thea_jz(i,j)=5/4*pi;

           end

           if thea_jz(i,j)>=21/16*pi&&thea_jz(i,j)<23/16*pi

              thea_jz(i,j)=11/8*pi;

           end

           if thea_jz(i,j)>=23/16*pi&&thea_jz(i,j)<25/16*pi

               thea_jz(i,j)=1/2*pi;

           end

       end

end

    for i=1:middle_mu1

       for j=1:middle_mu2

           for a=1:15

               for b=1:15

                   c=a+(i-1)*15;

                   d=b+(j-1)*15;

                   theta_all(c,d)=thea_jz(i,j);

               end

           end

       end

    end

  %%%%%%%%%%%%%%%%%%%%%thea=1/2*pi

    Sx1=5;

    Sy1=5;

    theta1=1/2*pi;

    f1=0.125;

    if isa(I,'double')~=1

        I = double(I);

    end

    for x = -fix(Sx1):fix(Sx1)

    for y = -fix(Sy1):fix(Sy1)

        xPrime = x * cos(theta1) + y * sin(theta1);

        yPrime = y * cos(theta1) - x * sin(theta1);

        G1(fix(Sx1)+x+1,fix(Sy1)+y+1) = exp(-.5*((xPrime/Sx1)^2+(yPrime/Sy1)^2))*cos(2*pi*f1*xPrime);

    end

    end

    Imgabout = conv2(I,double(imag(G1)),'same');

    Regabout = conv2(I,double(real(G1)),'same');

    Gabout1 = sqrt(Imgabout.*Imgabout + Regabout.*Regabout)/10;

    for i=1:xs1-8

        for j=1:xs2-1

            if theta_all(i,j)~=1/2*pi&&theta_all(i,j)~=11/8*pi

                 Gabout1(i,j)=0;

            end

            if Gabout1(i,j)>160

                Gabout1(i,j)=256;

            end

        end

    end

%  Gabout1=uint8(Gabout1);

%   F1=im2bw(Gabout1,85/255);1

%   h=ones(3,3)/15;

%   Gabout1=imfilter(F1,h);

%figure,imshow(Gabout1,[]);

% %%%%%%%%%%%%%%%%%%%% thea=5/8*pi

    Sx2=7;

    Sy2=7;

    theta2=5/8*pi;

    f2=0.1;

    if isa(I,'double')~=1

        I = double(I);

    end

    for x = -fix(Sx2):fix(Sx2)

    for y = -fix(Sy2):fix(Sy2)

        xPrime = x * cos(theta2) + y * sin(theta2);

        yPrime = y * cos(theta2) - x * sin(theta2);

        G2(fix(Sx2)+x+1,fix(Sy2)+y+1) = exp(-.5*((xPrime/Sx2)^2+(yPrime/Sy2)^2))*cos(2*pi*f2*xPrime);

    end

    end

    Imgabout = conv2(I,double(imag(G2)),'same');

    Regabout = conv2(I,double(real(G2)),'same');

    Gabout2 = sqrt(Imgabout.*Imgabout + Regabout.*Regabout)/10;

    for i=1:xs1-8

        for j=1:xs2-1

            if theta_all(i,j)~=5/8*pi

                 Gabout2(i,j)=0;

            end

            if Gabout2(i,j)>220

                Gabout2(i,j)=256;

            end

        end

    end

%   Gabout2=uint8(Gabout2);

%   F1=im2bw(Gabout2,20/255);

%   h=ones(3,3)/15;

%   Gabout2=imfilter(F1,h);

%figure,imshow(Gabout2,[]);

    %%%%%%%%%%%%%%%%%%%% thea=3/4*pi  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

%     Sx3=13;

%     Sy3=13;

%     theta3=3/4*pi;

%     f3=0.09;

%     if isa(I,'double')~=1   

%         I = double(I);

%     end

%     for x = -fix(Sx3):fix(Sx3)

%     for y = -fix(Sy3):fix(Sy3)

%         xPrime = x * cos(theta3) + y * sin(theta3);

%         yPrime = y * cos(theta3) - x * sin(theta3);

%         G3(fix(Sx3)+x+1,fix(Sy3)+y+1) = exp(-.5*((xPrime/Sx3)^2+(yPrime/Sy3)^2))*cos(2*pi*f3*xPrime);

%     end

%     end

%     Imgabout = conv2(I,double(imag(G3)),'same');

%     Regabout = conv2(I,double(real(G3)),'same');

%     Gabout3 = sqrt(Imgabout.*Imgabout + Regabout.*Regabout)/10;

%     for i=1:xs1-1

%         for j=1:xs2-1

%             if theta_all(i,j)~=3/4*pi

%                  Gabout3(i,j)=0;

%             end

% %              if Gabout3(i,j)>150

% %                 Gabout3(i,j)=256;

% %             end

%         end

%      end

% Gabout3=uint8(Gabout3);

% F1=im2bw(Gabout3,88/255);

% h=ones(3,3)/15;

% Gabout3=imfilter(F1,h);

% figure,imshow(Gabout3,[]);  

    %%%%%%%%%%%%%%%%%%%% thea=7/8*pi

    Sx4=7;

    Sy4=7;

    theta4=7/8*pi;

    f4=0.12;

    if isa(I,'double')~=1

        I = double(I);

    end

    for x = -fix(Sx4):fix(Sx4)

    for y = -fix(Sy4):fix(Sy4)

        xPrime = x * cos(theta4) + y * sin(theta4);

        yPrime = y * cos(theta4) - x * sin(theta4);

        G4(fix(Sx4)+x+1,fix(Sy4)+y+1) = exp(-.5*((xPrime/Sx4)^2+(yPrime/Sy4)^2))*cos(2*pi*f4*xPrime);

    end

    end

    Imgabout = conv2(I,double(imag(G4)),'same');

    Regabout = conv2(I,double(real(G4)),'same');

    Gabout4 = sqrt(Imgabout.*Imgabout + Regabout.*Regabout)/10;

    for i=1:xs1-8

        for j=1:xs2-1

            if theta_all(i,j)~=7/8*pi&&theta_all(i,j)~=3/4*pi&&theta_all(i,j)~=pi

                 Gabout4(i,j)=0;

            end

            if Gabout4(i,j)>150

                Gabout4(i,j)=256;

            end

        end

    end

%     Gabout4=uint8(Gabout4);

%     F1=im2bw(Gabout4,30/255);

%     h=ones(3,3)/15;

%     Gabout4=imfilter(F1,h);

%figure,imshow(Gabout4,[]);

%     

%     %%%%%%%%%%%%%%%%%%%% thea=pi  xxxxxxxxxxxxxxxxxxxxxxxxxxx

%     Sx5=6;

%     Sy5=6;

%     theta5=pi;

%     f5=0.08;

%     if isa(I,'double')~=1   

%         I = double(I);

%     end

%     for x = -fix(Sx5):fix(Sx5)

%     for y = -fix(Sy5):fix(Sy5)

%         xPrime = x * cos(theta5) + y * sin(theta5);

%         yPrime = y * cos(theta5) - x * sin(theta5);

%         G5(fix(Sx5)+x+1,fix(Sy5)+y+1) = exp(-.5*((xPrime/Sx5)^2+(yPrime/Sy5)^2))*cos(2*pi*f5*xPrime);

%     end

%     end

%     

%     Imgabout = conv2(I,double(imag(G5)),'same');

%     Regabout = conv2(I,double(real(G5)),'same');

%     Gabout5 = sqrt(Imgabout.*Imgabout + Regabout.*Regabout)/10;

%         

%     for i=1:xs1-1

%         for j=1:xs2-1

%             if theta_all(i,j)~=pi

%                  Gabout5(i,j)=0;

%             end

%         end

%          end

%     Gabout5=uint8(Gabout5);

%     F1=im2bw(Gabout5,230/255);

%     h=ones(3,3)/15;

%     Gabout5=imfilter(F1,h);

%     figure,imshow(Gabout5,[]);

    for i=1:xs1

        for j=1:xs2

            gabout_18=[Gabout1(i,j);Gabout2(i,j);Gabout4(i,j);Gabout6(i,j)];

            gabout_last(i,j)=max(gabout_18); 

            if gabout_last(i,j)~=256

                gabout_last(i,j)=0;

            end

           if i<20||j<30||i>xs1-20||j>xs2-20

                gabout_last(i,j)=256;

            end

        end

    end

figure,imshow(gabout_last,[]);

3 仿真结果

4 参考文献

[1]贡玉南, 华建兴, 黄秀宝. 基于匹配GabOr滤波器的规则纹理缺陷检测方法[J]. 中国图象图形学报, 2001, 006(007):624-628.

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。