Gabor滤波进行目标图像纹理特征的提取

1.傅里叶变换

1) 简介

数字图像处理的方法主要分成两大部分：空域分析法和频域分析法。空域分析法就是对图像矩阵进行处理；频域分析法是通过图像变换将图像从空域变换到频域，从另外一个角度来分析图像的特征并进行处理。频域分析法在图像增强、图像复原、图像编码压缩及特征编码压缩方面有着广泛应用。

如果一个信号f(t)在上满足：

① f(t)在任一有限区间上满足狄氏条件；

② f(t)在上绝对可积即

就可以通过傅里叶变换把时域信号f(t)转化到频域进行处理：

然后再通过傅里叶反变换把频域信号转化到时域：

傅里叶变换是线性系统分析的有力工具，提供了一种把时域信号转换到频域进行分析的途径，时域和频域之间是一对一的映射关系。图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对 于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。

傅立叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅立叶变换就表示f的 谱。从纯粹的数学意义上看，傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将 图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为 灰度分布函数。

2) 不足之处

经典Fourier变换只能反映信号的整体特性（时域，频域）。对傅里叶谱中的某一频率，无法知道这个频率是在什么时候产生的。从傅里叶变换的定义也可看出，傅里叶变换是信号在整个时域内的积分，因此反映的是信号频率的统计特性，没有局部化分析信号的功能。另外，要求信号满足平稳条件。傅里叶变换时域和频域是完全分割开来的。

l 由式可知，要用Fourier变换研究时域信号频谱特性，必须要获得时域中的全部信息；

l 信号在某时刻的一个小的邻域内发生变化，那么信号的整个频谱都要受到影响，而频谱的变化从根本上来说无法标定发生变化的时间位置和发生变化的剧烈程度。也就是说，Fourier变换对信号的齐性不敏感。不能给出在各个局部时间范围内部频谱上的谱信息描述。然而在实际应用中齐性正是我们所关心的信号局部范围内的特性。如，音乐，语言信号等。即：局部化时间分析，图形边缘检，地震勘探反射波的位置等信息极重要。

l 为解决傅里叶变换的局限性，产生了Gabor变换和小波变换。

2.Gabor变换

Gabor变换是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)加一个滑动窗，再作傅里叶变换。

设函数f为具体的函数，且，则Gabor变换定义为

其中，，是高斯函数，称为窗函数。其中a>0,b>0.

是一个时间局部化的“窗函数”。其中，参数b用于平行移动窗口，以便于覆盖整个时域。对参数b积分，则有

信号的重构表达式为

Gabor取g(t)为一个高斯函数有两个原因：一是高斯函数的Fourier变换仍为高斯函数，这使得Fourier逆变换也是用窗函数局部化，同时体现了频域的局部化；二是Gabor变换是最优的窗口Fourier变换。其意义在于Gabor变换出现之后，才有了真正意义上的时间－频率分析。即Gabor变换可以达到时频局部化的目的：它能够在整体上提供信号的全部信息而又能提供在任一局部时间内信号变化剧烈程度的信息。简言之，可以同时提供时域和频域局部化的信息。

 

② 窗口的宽高关系

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

 矩形时间――频率窗：宽为，高。

由此，可以看出Gabor变换的局限性：时间频率的宽度对所有频率是固定不变的。实际要求是：窗口的大小应随频率而变化，频率高窗口应愈小，这才符合实际问题中的高频信号的分辨率应比低频信号的分辨率要低。

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

① 首先选取核函数

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

则其对偶函数为

② 离散Gabor变换的表达式

其中，

是的对偶函数，二者之间有如下双正交关系。

4) Gabor变换的解析理论

Gabor变换的解析理论就是由g(t)求对偶函数的方法。

定义g(t)的Zak变换为

可以证明对偶函数可由下式求出：

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

5) 适用条件

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

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

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

6) 应用

① 暂态信号检测

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

② 图象分析与压缩

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

3. 二维Gabor滤波器

用Gabor 函数形成的二维Gabor 滤波器具有在空间域和频率域同时取得最优局部化的特性，因此能够很好地描述对应于空间频率(尺度)、空间位置及方向选择性的局部结构信息。Gabor滤波器的频率和方向表示接近人类视觉系统对于频率和方向的表示，并且它们常备用于纹理表示和描述。在图像处理领域，Gabor滤波器是一个用于边缘检测的线性滤波器。，在空域，一个2维的Gabor滤波器是一个正弦平面波和高斯核函数的乘积。Gabor滤波器是自相似的，也就是说，所有Gabor滤波器都可以从一个母小波经过膨胀和旋转产生。实际应用中，Gabor滤波器可以在频域的不同尺度，不同方向上提取相关特征。

 

1)定义

空域来看：是高斯核函数调制正弦平面波

s(x,y)是复杂的正弦函数，相当于载波；w(x,y)是2维高斯函数包迹。

（u0,v0）定义了正弦平面波的时域频率，在极坐标中可用f和Θ来表示。

a,b 为x和y方向的椭圆高斯的方差

K=1/ab 为高斯包迹的参数

r 为角度旋转的下标

Θ为旋转角度

（x0,y0）为函数峰值，也是接受域的中心

        

f(x,y) f(x',y')

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

  

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

  
2.

二．Gabor 滤波器简介（部分资料来自维基百科）

   在图像处理领域，Gabor滤波器是一个用于边缘检测的线性滤波器。Gabor滤波器的频率和方向表示接近人类视觉系统对于频率和方向的表示，并且它们常备用于纹理表示和描述。在空域，一个2维的Gabor滤波器是一个正弦平面波和高斯核函数的乘积。Gabor滤波器是自相似的，也就是说，所有Gabor滤波器都可以从一个母小波经过膨胀和旋转产生。实际应用中，Gabor滤波器可以在频域的不同尺度，不同方向上提取相关特征。

三．Gabor滤波器公式化定义 

公式中：

λ：正弦函数波长；

θ：Gabor核函数的方向 

ψ：相位偏移

σ：高斯函数的标准差 

γ： 空间的宽高比（这个没太理解）

CGaborFilter::CGaborFilter(float dLambda, float dTheta, float dRatio_S2L, float dGamma, float dPhi)  2{  3    Lambda = dLambda;  4    Theta = dTheta;  5    sigma = dLambda*dRatio_S2L;  6    Gamma = dGamma;  7    Phi = dPhi;  8    m_pGaborFilter = NULL;  9    bParam = 1; 10} 11 12 13CGaborFilter::~CGaborFilter(void) 14{ 15    cvReleaseMat(&m_pGaborFilter); 16} 17 18void CGaborFilter::Init() 19{ 20    float dtmp; 21    int itmp; 22    if(is_param() == 0) 23    { 24        printf("The parameters are not enough!"); 25    } 26    else 27    { 28        dtmp = sqrt(48*pow(sigma,2)+1); 29        itmp = cvRound(dtmp); 30        if(itmp%2 == 0) 31            itmp ++; 32        GaborWindow.height = GaborWindow.width = 16; 33        bInit = 1; 34 35        create_kernel(); 36    } 37} 38 39void CGaborFilter::Init(float dSigma,float dTheta,float dPhi) 40{ 41    float dtmp; 42    int itmp; 43 44    sigma = dSigma; 45    Theta = dTheta; 46    Phi = dPhi; 47    Gamma = GAMMA; 48    Lambda = sigma/RATIO_S2L; 49    bParam = 1; 50 51    dtmp = sqrt(24*pow(sigma,2)); 52    itmp = cvRound(dtmp); 53    if(itmp%2 == 0) 54        itmp ++; 55    GaborWindow.height = GaborWindow.width = itmp; 56    bInit = 1;  57 58    create_kernel(); 59} 60 61void CGaborFilter::Init(float dLambda,float dTheta, float dPhi,float dGamma) 62{ 63    float dtmp; 64    int itmp; 65 66    Lambda = dLambda; 67    Theta  = dTheta; 68    Phi    = dPhi; 69    Gamma  = dGamma; 70    sigma  = Lambda * RATIO_S2L; 71    bParam = 1; 72 73    dtmp = sqrt(24*pow(sigma,2)); 74    itmp = cvRound(dtmp); 75    if(itmp%2 == 0) 76        itmp ++; 77    GaborWindow.height = GaborWindow.width = itmp; 78    bInit = 1;  79 80    create_kernel(); 81} 82 83void CGaborFilter::create_kernel() 84{ 85    float tmp1,tmp2,xtmp,ytmp,re; 86    int i,j,x,y; 87 88    if(is_init() == 0) 89        printf("The parameters haven't been initialed!"); 90 91 92    else{ 93 94 95        m_pGaborFilter = cvCreateMat(GaborWindow.height,GaborWindow.width,CV_32FC1); 96        for(i= 0; i< GaborWindow.height; i++) 97            for(j = 0; j< GaborWindow.width; j++) 98            { 99                x = j - GaborWindow.width/2;100                y = i - GaborWindow.height/2;101102103                xtmp = (float)x*cos(Theta) - (float)y*sin(Theta);104                ytmp = (float)x*sin(Theta) + (float)y*cos(Theta);105106                tmp1 = exp(-(pow(xtmp,2)+pow(ytmp*Gamma,2))/(2*pow(sigma,2)));107                tmp2 = cos(2*PI*xtmp/Lambda + Phi);108                // int p=sizeof(float);109                re   = tmp1*tmp2;110                cvSetReal2D((CvMat*)m_pGaborFilter,i,j,re);111112            }113            bKernel = 1;114    }115}116117IplImage * CGaborFilter::get_Image()118{119    if(is_kernel() == 0)120    {121        printf("The filter hasn't bee created!");122    }123    else124    {125        IplImage *pImg = cvCreateImage(GaborWindow,IPL_DEPTH_32F,1);126        IplImage *pImgU8 = cvCreateImage(GaborWindow,IPL_DEPTH_8U,1);127        CvMat * pMat = cvCreateMat(GaborWindow.height,GaborWindow.width,CV_32FC1);128129        cvCopy(m_pGaborFilter,pImg);130        //pImg->imageData = (char *)pMat->data;131        cvNormalize((IplImage*)pImg, (IplImage*)pImg,0,255,CV_MINMAX,NULL);132        cvConvertScaleAbs(pImg,pImgU8,1,0);133        return pImgU8;134    }135}136137IplImage * CGaborFilter::do_filter(const IplImage * src)138{139    if(is_kernel()==false)140    {141        printf("The Gabor Kernel has not been created!");142    }143    else{144145        IplImage *pDestImage = cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);146        // IplImage * pGaborImage = get_Imge();147        // CvMat GaborKernel = cvMat(pGaborImage->height,pGaborImage->width,CV_8U,pGaborImage->imageData);148        IplImage *tmpImg = cvCloneImage(src);149        IplImage *tmpGrayImg = cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);150151        if(tmpImg->nChannels != 1)152        {153            cvCvtColor(tmpImg,tmpGrayImg,CV_BGR2GRAY);154        }155        else 156        {157            cvReleaseImage(&tmpGrayImg);158            tmpGrayImg = tmpImg;159        }160        CvMat * pGaborKernel = get_Mat();161162        cvFilter2D(tmpGrayImg,pDestImage,pGaborKernel,cvPoint((GaborWindow.width-1)/2,(GaborWindow.height-1)/2));163164        cvReleaseImage(&tmpImg);165        return  pDestImage;  166    }167}

3.

人脸光照调整一直是人脸识别问题中的难点，作者就不同处理方法，并结合其在实际应用中的表现，在此分章节谈一些个人看法，有不当之处，还望各园友指正。

光照调整主要分在频率域和空间域的处理，频率域中有DCT变换、小波变换等。而在空间域有直方图均衡化、Gamma校正等。作者在此只介绍一些主流并有较广适应范围的方法。首先，我们来谈谈DCT变换吧：

DCT（离散余弦变换）对高相关性的数据（信号），具有非常好的能量聚焦性，经过变换，信号能量的绝大部分被集中到变换域的少数系数上。因此，对于受光照影响的图像，我们只需要修改很少的频域系数，就可以对图像的光照做出较好调整，避免了需要调节多个参数以适合不同图像的问题，操作简便易行。

1.       一维DCT变换的实现步骤：

1)        计算DCT变换的点数，并对时域空间进行延拓；

2)        调用一维傅里叶变换；

3)        调整系数并存储；

2.       二维DCT变换的实现步骤：

1)        计算进行二维图像DCT变换的高度和宽度，如果不是2的整数次幂则要进行调整，并计算在水平和垂直方向上变换时迭代的次数；

2)        用一维DCT变换进行水平方向上的变换

3)        用一维DCT变换进行垂直方向上的变换；

4)        得到二维离散余弦变换系数并存储。

3.       试验结果：

图1是对高曝光图像的二维DCT变换结果，其中（b）是DCT变换后的频谱图像，可以看出图像的低频能量都集中在左上角区域，而向右下角方向，频率越来越高。图（c）和原始图像相比，脸上（左脸）的高光照部分得到了一定的抑制。

     　　　　 　　　　　　　　　　

　　（a）原图　　　　　　　　（b）DCT变换的频谱图 　　（c）将频谱图中的低频减去并反DCT变换的结果

　　　　　　　　　　　　　　　　　　　　图1 高曝光图像进行DCT变换的结果

图2是光照不足的二维DCT变换结果，原始图像中左边人脸基本看不到什么信息，经过处理后可以看到左眼信息。

                                             

        （a）原图　　　　　　　　（b）DCT变换的频谱图 　　（c）将频谱图中的低频减去并反DCT变换的结果

图2 “阴阳脸”图像进行DCT变换的结果

4.       DCT方法小结：

DCT变换的方法，只需要对频域图像做极少的调整（如将最左上角的频率置0），就可以达到对整体光照的调整，不需要像在空间域中那样，不断的调整参数，修改阈值等步骤，这是它的优点。但是，DCT变换的时间稍长，而且对于光照复杂度比较大的图像，其调整的能力也比较有限。当然，也有不少人对其进行了改进，比如在Log域的DCT变换等。网上也有不少关于其改进的文献，在此就不一一赘述。作者在此介绍的方法，代表了频率域调整光照的一般思路，希望对大家有作用。

 4.

在图像处理、模式识别以及计算机视觉等领域中，Gabor 滤波器得到了广泛的应用。 用Gabor 函数形成的二维Gabor 滤波器具有在空间域和频率域同时取得最优局部化的特性，与人类生物视觉特性很相似，因此能够很好地描述对应于空间频率(尺度)、空间位置及方向选择性的局部结构信息。

　　Gabor变换是一种短时傅里叶变换方法，其实质是在傅里叶变换中加入一个窗函数，通过窗函数来实现信号的时频分析。当选取高斯函数作为窗函数时，短时傅里叶变换称为Gabor变换。

　　常用的偶对称二维Gabor滤波器可表示为：

1. 不同方向下的Gabor滤波器：

　　　　　　　　　　　　　　　　　　　　图1 不同方向上的滤波器

　　在实际应用时，可以根据检测对象的方向趋势，选择合适的方向参数进行滤波。如在检测人脸的五官时，可以根据人脸的偏转角度进行滤波，可以使特征点的定位更加准确。

　　2.  不同频率下的滤波器：

　　　　　　　　　　　　　　　　　　　　　　图2 不同频率下的滤波器

　　从图2可以看出随着的变化，Gabor滤波器中出现了很多宽窄与纹理不同的明暗条纹。当滤波器纹理与图像作用时，滤波器覆盖下的局部纹理频率与滤波器的频率越接近响应就越大，反之越小。

　　3.   试验结果:

　　在“人脸光照调整之DCT变换”随笔中，原始图像经过DCT变换处理后，并不能完全去除光照在人脸上分布不均的影响，而且人脸的本真信息也难以被全部表达。为此在DCT变换的基础上，用Gabor滤波对其进行再处理，可以达到更好的结果。

图三 基于DCT变换的Gabor滤波

图四 基于DCT变换的Gabor滤波

　　图三（c）是在（b）图基础上做的Gabor滤波，效果显示已基本完全消除了高曝光对图像的影响。同理，图四（c）的右边脸的光照也被抑制下来。图四（d）是对原始图像直接做Gabor滤波，虽然局部效果较（c）图更清晰，但整体纹理没有（c）图平滑，这样会给后续特征点定位的收敛性带来影响，因此定位效果欠稳定。

　　这两种方法合在一起使用，时间开销还是挺大的，在人脸识别等实时系统中，需要优化或精简。一般情况下，就单比处理效果和稳定性，Gabor要由于DCT变换。因此，在容许情况下，我们可以只取Gabor对图像进行处理。比如，作者在“眼睛定位”随笔中，就只用Gabor滤波对人脸处理，以提高眼睛定位精度。

　　下面，作者再贴几张图，看看这两种方法合在一起时，对AAM的帮助。

图五 光照调整对AAM定位的帮助

　　图五中的（a）图是AAM对原始图像直接定位的结果，（b）图是在去光照后的定位效果。比较两组图像，可以很明显的看到（b）图的定位精度有了大幅度的提高。

参考文章：

Gabor滤波简介和实现(Matlab，OpenCV)

2D-Gabor Matlab实现