图像分析：Gabor滤波器解析与编程-蒲公英云

我现在的研究方向主要在于图像分析这一块，下面借整理过去学习的技术资料重新复习一遍。

图像分析：SIFT特征解析与编程

图像分析：Gabor滤波器解析与编程

图像分析：积分图像与代码实现

图像分析：LBP特征解析与代码

图像分析：分水岭算法解析

Blogger： LinJM 微博（ LinJM-机器视觉）

在图像分析里面，Gabor滤波器是应用非常广泛的一种工具，也算是我本科第一次接触图像分析领域时用到的一种工具，刚开始学Gabor的时候，因为找不到门路，看到它的几种形式特别伤脑筋，不知道这几种之间到底是个什么关系，为它的复杂形式纠结了好久。后来随着学习的深入，发现Gabor其实是一种很“傻瓜式的万用工具”，可以应用在很多场合，就是参数设置麻烦了点。在应用的时候，只要把它当成一个滤波器来对待就行，当然，若是有时间那是最好去理解一下它公式的含义。

Gabor滤波器的历史过程

十九世纪初，法国科学家傅里叶在向巴黎科学院呈交的关于热传导的著名论文中提出了Fourier级数，其意义是无法估量的。今天，Fourier分析方法已成为了各种信号数据处理中最基本的数学工具。学过《信号与系统》和《数字信号处理》的同学都知道，Fourier变换是在整体上将信号分解为不同的频率分量，但是它只适用于统计量不随时间变化的平稳信号，因为它并不能告之某种频率发生在哪些时间内，而这对于非平稳信号是十分重要的。

为了研究信号在局部范围内的频率特性，Dennis Gabor于1946年在“Theory of communication”一文中提出了著名的“窗口”傅里叶变换（也叫短时Fourier变换，STFT），即Gabor变换。

Gabor滤波器解析

Gabor滤波器实质上一种加了高斯窗的傅里叶变换，对于图像来说，窗函数决定了它在空域的局部性，所以可以通过移动窗口的中心來获得不同位置的空间域信息。此外，由于髙斯函数在经过傅里叶变换以后仍然是高斯函数，这就使得Gabor变换在频域上仍然是局部的。
1985年，Daugman将Gabor函数扩展为2维形式，并在此基础上构造了2D Gabor滤波器。人们惊奇地发现，2D Gabor滤波器不仅同样可以同时获取时间和频率域的最小不确定性，而且与哺乳动物视网膜神经细胞的接收模型相吻合。
Daugman认为，尽管Gabor滤波器的基函数不能构成一个完备的正交集，Gabor滤波器也可以看作是一种小波滤波器。他给出了2D Gabor滤波器基函数的一般形式：

SouthEast

从上式可以看出，Gabor函数是一个被复正弦函数调制的Gaussian函数。其中λ和θk分别是正弦波的波长和方向；θk的定义如下：

SouthEast 1

k决定了滤波器方向的个数；σ_x和σ_y分别为高斯包络在x方向和y方向上的标准差，它们决定了高斯包络的空间扩展；

SouthEast 2

Gabor滤波器的参数设置：目前主流的方法是基于一定的约束条件，通过实验进行Gabor滤波器的参数选择。

Gabor滤波器的缺点：计算量过大。因为（1）Gabor滤波器的实现上需要耗费时间的卷积操作，然而却没有行之有效的快速算法；（2）Gabor特征维数较高，Gabor滤波器的优点是多通道、多分辨率分析，通过Gabor滤波器抽取的特征通常能够取得较高的识别率，但也产生了高维数的Gabor特征矢量，由此带来了较大的计算量和存储负担。例如，对于96×48的样本图像来说，当选择3尺度、8方向时，Gabor特征矢量会达到110，592（86×48×3×8）。

Gabor滤波器程序实现

Matlab实现：

function [G,gabout] = gaborfilter1(I,Sx,Sy,f,theta)
%The Gabor filter is basically a Gaussian (with variances sx and sy along x and y-axes respectively)
%modulated by a complex sinusoid (with centre frequencies U and V along x and y-axes respectively) 
%described by the following equation
%%%%%%%%%%%%%%%%%%%%%%
%                                  -1     x' ^      y' ^             
%%% G(x,y,theta,f) =  exp ([----{(----) 2+(----) 2}])*cos(2*pi*f*x'); cos代表实部
%                                   2     sx'       sy'
%%% x' = x*cos(theta)+y*sin(theta);
%%% y' = y*cos(theta)-x*sin(theta);
%% Describtion :
%% I : Input image
%% Sx & Sy : Variances along x and y-axes respectively 方差
%% f : The frequency of the sinusoidal function
%% theta : The orientation of Gabor filter
%% G : The output filter as described above
%% gabout : The output filtered image
% %%isa判断输入参量是否为指定类型的对象
if isa(I,'double')~=1 
    I = double(I);
end
%%%%Sx,Sy在公式里分别表示Guass函数沿着x,y轴的标准差，相当于其他的gabor函数中的sigma. 
%%同时也用Sx,Sy指定了gabor滤波器的大小。（滤波器矩阵的大小）
%%这里没有考虑到相位偏移.fix(n)是取小于n的整数（往零的方向靠）
for x = -fix(Sx):fix(Sx)
    for y = -fix(Sy):fix(Sy)
        xPrime = x * cos(theta) + y * sin(theta);
        yPrime = y * cos(theta) - x * sin(theta);
        G(fix(Sx)+x+1,fix(Sy)+y+1) = exp(-.5*((xPrime/Sx)^2+(yPrime/Sy)^2))*cos(2*pi*f*xPrime);
    end
end
Imgabout = conv2(I,double(imag(G)),'same');
Regabout = conv2(I,double(real(G)),'same');
gabout = sqrt(Imgabout.*Imgabout + Regabout.*Regabout);

OpenCV中的实现：

#include "precomp.hpp"
/*
 Gabor filters and such. To be greatly extended to have full texture analysis.
 For the formulas and the explanation of the parameters see:
 http://en.wikipedia.org/wiki/Gabor_filter
*/
cv::Mat cv::getGaborKernel( Size ksize, double sigma, double theta,
                            double lambd, double gamma, double psi, int ktype )
{
    double sigma_x = sigma;
    double sigma_y = sigma/gamma;
    int nstds = 3;
    int xmin, xmax, ymin, ymax;
    double c = cos(theta), s = sin(theta);
    if( ksize.width > 0 )
        xmax = ksize.width/2;
    else
        // cvRound返回和参数最接近的整数值；fabs返回浮点数的绝对值
        xmax = cvRound(std::max(fabs(nstds*sigma_x*c), fabs(nstds*sigma_y*s)));
    if( ksize.height > 0 )
        ymax = ksize.height/2;
    else
        ymax = cvRound(std::max(fabs(nstds*sigma_x*s), fabs(nstds*sigma_y*c)));
    xmin = -xmax;
    ymin = -ymax;
    //CV_Assert（）若括号中的表达式值为false，则返回一个错误信息。
    CV_Assert( ktype == CV_32F || ktype == CV_64F );
    //初始化kernel矩阵
    Mat kernel(ymax - ymin + 1, xmax - xmin + 1, ktype);
    double scale = 1;
    double ex = -0.5/(sigma_x*sigma_x);
    double ey = -0.5/(sigma_y*sigma_y);
    double cscale = CV_PI*2/lambd;
    for( int y = ymin; y <= ymax; y++ )
        for( int x = xmin; x <= xmax; x++ )
        {
            double xr = x*c + y*s;
            double yr = -x*s + y*c;
            // real gabor
            double v = scale*exp(ex*xr*xr + ey*yr*yr)*cos(cscale*xr + psi);
            if( ktype == CV_32F )
                kernel.at<float>(ymax - y, xmax - x) = (float)v;
            else
                kernel.at<double>(ymax - y, xmax - x) = v;
        }
    return kernel;
}