本文转自 taotao1233

本文的目的是用C实现生成Gabor模版，并对图像卷积。并简单提一下，Gabor滤波器在纹理特征提取上的应用。

一、什么是Gabor函数（以下内容含部分翻译自维基百科）

  在图像处理中，Gabor函数是一个用于边缘提取的线性滤波器。Gabor滤波器的频率和方向表达同人类视觉系统类似。研究发现，Gabor滤波器十分适合纹理表达和分离。在空间域中，一个二维Gabor滤波器是一个由正弦平面波调制的高斯核函数。

  还有，生物学实验发现，Gabor滤波器可以很好地近似单细胞的感受野函数（光强刺激下的传递函数），什么视皮层内的超柱，bla...bla，总之是这方面仿生的数学模型。

  另外，网上有一种说法，gabor分为实部和虚部，用实部进行滤波后图像会平滑；虚部滤波后用来检测边缘。【来自百度知道某个大神的回答】，我查了文献，发现的确有人用Gabor的奇函数部分做边缘提取（《基于Gabor滤波器的边缘检测算法》 无线电工程2000年第3卷第30期）。另外，从我的实验结果也有类似的发现。暂且认为这个对的吧。

  Gabor滤波器的脉冲响应，可以定义为一个正弦波（对于二维Gabor滤波器是正弦平面波）乘以高斯函数。由于乘法卷积性质，Gabor滤波器的脉冲响应的傅立叶变换是其调和函数的傅立叶变换和高斯函数傅立叶变换的卷积。该滤波器由实部和虚部组成，二者相互正交。一组不同频率不同方向的Gabor函数数组对于图像特征提取非常有用。

 下面给出二维Gabor函数的数学表达：

复数表达：

实数部分：

虚数部分：

其中：

和

下面介绍公式中各个参数的含义，及参数如何配置问题【都从老外那翻译来的】：

波长（λ）：它的值以像素为单位指定，通常大于等于2.但不能大于输入图像尺寸的五分之一。

方向（θ）：这个参数指定了Gabor函数并行条纹的方向，它的取值为0到360度

相位偏移（φ）：它的取值范围为-180度到180度。其中，0he180度分别对应中心对称的center-on函数和center-off函数，而-90度和90度对应反对称函数。

长宽比（γ）：空间纵横比，决定了Gabor函数形状（support，我翻译为形状）的椭圆率（ellipticity）。当γ= 1时，形状是圆的。当γ< 1时，形状随着平行条纹方向而拉长。通常该值为0.5

带宽（b）：Gabor滤波器的半响应空间频率带宽b和σ/ λ的比率有关，其中σ表示Gabor函数的高斯因子的标准差，如下：

σ的值不能直接设置，它仅随着带宽b变化。带宽值必须是正实数，通常为1，此时，标准差和波长的关系为：σ= 0.56 λ。带宽越小，标准差越大，Gabor形状越大，可见平行兴奋和抑制区条纹数量越多。

下面给出，不同参数配置下的Gabor核函数效果图，大小均100*100：

a.波长对比组【方向为：0，相位偏移量为：0，纵横比率为：0.5，带宽为：1，下图波长分别为5，10，15】

b.方向对比组【波长为：10，相位偏移量为：0，空间纵横比为：0.5，带宽为：1，方向分别为：0，45，90】

c.相位偏移量对比组【波长为：10，方向为：0，空间纵横比：0.5，带宽：1，相位偏移量分别为：0，180，-90，90】

d.空间纵横比对比组【波长：10，相位偏移量：0，方向：0，带宽：1，空间纵横比分别为：0.5，1】

e.带宽对比组【波长：10，方向：0，相位偏移量：0，空间纵横比：0.5，带宽分别为：0.5，1，2】

二、gabor函数实现：

matlab版本，我从pudn上找来的，但他的gabor函数，我没怎么看明白：

gabor函数：

function gabor_k = compute(x,y,f0,theta)

r = 1; g = 1;

x1 = x*cos(theta) + y*sin(theta);

y1 = -x*sin(theta) + y*cos(theta);

gabor_k = f0^2/(pi*r*g)*exp(-(f0^2*x1^2/r^2+f0^2*y1^2/g^2))*exp(i*2*pi*f0*x1);

%绘制一个Gabor滤波器的空域和频域函数图

clear;

x = 0;

theta = 0;

f0 = 0.2;

for i = linspace(-15,15,50)

    x = x + 1;

    y = 0;

   for j = linspace(-15,15,50)

       y = y + 1;

       z(y,x)=compute(i,j,f0,theta);

   end

end

x = linspace(-15,15,50);

y = linspace(-15,15,50);

surf(x,y,real(z))

title('Gabor filter:real component');

xlabel('x');

ylabel('y');

zlabel('z');

figure(2);

surf(x,y,imag(z))

title('Gabor filter:imaginary component');

xlabel('x');

ylabel('y');

zlabel('z');

Z = fft2(z);

u = linspace(-0.5,0.5,50);

v = linspace(-0.5,0.5,50);

figure(3);

surf(u,v,abs(fftshift(Z)))

title('Gabor filter:frequency component');

xlabel('u');

ylabel('v');

zlabel('Z');

运行结果：

%4个方向的Gabo滤波器通过图像显示

clear;

x = 0;

theta = pi*3/4;%用弧度0,pi/4,pi/2,pi*3/4

f0 = 0.2; 

for i = linspace(-15,15,50)

    x = x + 1;

    y = 0;

   for j = linspace(-15,15,50)

       y = y + 1;

       z(y,x)=compute(i,j,f0,theta);

   end

end

z_real = real(z);

m = min(z_real(:));

z_real = z_real+abs(m);

M = max(z_real(:));

imshow(1/M*z_real);

figure(2)

z_imag = imag(z);

m = min(z_imag(:));

z_imag = z_imag+abs(m);

M = max(z_imag(:));

imshow(1/M*z_imag);

运行效果：

实数部分：

虚数部分：

%4个方向的Gabor滤波器对lena进行滤波

clear;

I = imread('.\pic\lena.bmp');

f0 = 0.2; 

count = 0;

for theta = [0,pi/4,pi/2,pi*3/4];%用弧度0,pi/4,pi/2,pi*3/4

   count = count + 1;

    x = 0;

   for i = linspace(-8,8,11)

       x = x + 1;

       y = 0;

       for j = linspace(-8,8,11)

           y = y + 1;

           z(y,x)=compute(i,j,f0,theta);

       end

   end

   figure(count);

   filtered = filter2(z,I,'valid');

    f = abs(filtered);

   imshow(f/max(f(:)))

end

运行效果：

好吧，不管他了。大概感受一下吧。由于我没看明白他的gabor函数怎么定义的，参数设置也不一样，实验结果很不相同，我希望我是对的，天地良心呐！！我只能按照维基百科给出的函数，编写了以下C代码：

// my_gabor.cpp : 定义控制台应用程序的入口点。

//

#include "stdafx.h"

#include<iostream>

#include <opencv2/core/core.hpp>  

#include <opencv2/highgui/highgui.hpp>  

#include <opencv2/imgproc/imgproc.hpp>

#include "math.h"

#define PI 3.1415926

#define N 4

using namespace std;

using namespace cv;

void m_filer(double *src,int height,int width,double *mask_rel,double *mask_img,int mW,int mH,int k)

{

       IplImage *tmp;

       double a,b,c;

       char res[20];         //保存的图像名称

       tmp = cvCreateImage(cvSize(width,height),IPL_DEPTH_8U,1);

       for (int i = 0;i < height;i++)

       {

              for (int j = 0;j < width;j++)

              {

                     a = 0.0;

                     b = 0.0;

                     c = 0.0;

                     //去掉靠近边界的行列，无法滤波，超出范围

                     if (i > int(mH/2) && i < height - int(mH/2) && j > int(mW) && j < width - int(mW/2))

                     {

                            for (int m = 0;m < mH;m++)

                            {

                                   for(int n = 0;n < mW;n++)

                                   {

                                          //printf("%f\n",src[(i+m-int(mH/2))*width+(j+n-int(mW))]);

                                          a += src[(i+m-int(mH/2))*width+(j+n-int(mW))]*mask_rel[m*mW+n];

                                          b += src[(i+m-int(mH/2))*width+(j+n-int(mW))]*mask_img[m*mW+n];

                                          //printf("%f,%f\n",a,b);

                                   }

                            }

                     }

                     c = sqrt(a*a+b*b);

                     c /= mW*mH;

                     tmp->imageData[i*width+j] = (unsigned char)c;

              }

       }

       sprintf(res,"result%d.jpg",k);

       cvSaveImage(res,tmp);

       cvReleaseImage(&tmp);

}

int _tmain(int argc, _TCHAR* argv[])

{

       IplImage *src;

       double *rel,*img,*src_data,xtmp,ytmp,tmp1,tmp2,tmp3,re,im;

       double Theta,sigma,Gamma,Lambda,Phi;        //公式中的5个参数

       int gabor_height,gabor_width,x,y;

       src = cvLoadImage("test.jpg",CV_LOAD_IMAGE_GRAYSCALE);

       gabor_height = 10;

       gabor_width = 10;

       Gamma = 1.0;

       Lambda = 10.0;

       sigma = 100;

       Phi = 0;

       rel = (double *)malloc(sizeof(double)*gabor_width*gabor_height);//实数部分

       img = (double *)malloc(sizeof(double)*gabor_width*gabor_height);//虚数部分

       src_data = (double *)malloc(sizeof(double)*src->widthStep*src->height);  //图像数据

       for (int i=0;i<src->height;i++)

       {

              for (int j=0;j<src->widthStep;j++)

              {

                     src_data[i*src->widthStep+j]=(unsigned char)src->imageData[i*src->widthStep+j];

                     //printf("%f\n",src_data[i*src->widthStep+j]);

              }

       }

       //构造gabor函数

       for (int k = 0;k < N;k++)                                  //定义N方向

       {

              Theta = PI*((double)k/N);

              for (int i = 0;i < gabor_height;i++)    //定义模版大小

              {

                     for (int j = 0;j < gabor_width;j++)

                     {

                            x = j - gabor_width/2;

                            y = i - gabor_height/2;

                            xtmp = (double)x*cos(Theta) + (double)y*sin(Theta);

                            ytmp = (double)y*cos(Theta) - (double)x*sin(Theta);

                            tmp1 = exp(-(pow(xtmp,2)+pow(ytmp*Gamma,2))/(2*pow(sigma,2)));

                            tmp2 = cos(2*PI*xtmp/Lambda + Phi);

                            tmp3 = sin(2*PI*xtmp/Lambda + Phi);

                            re = tmp1*tmp2;

                            im = tmp1*tmp3;

                            rel[i*gabor_width+j] = re;

                            img[i*gabor_width+j] = im;

                            //printf("%f,%f\n",re,im);

                     }

              }

              //用不同方向的GABOR函数对图像滤波并保存图片

              m_filer(src_data,src->height,src->width,rel,img,10,10,k);

       }

       free(rel);free(img);free(src_data);

       return 0;

}

运行效果：

大概就这样凑活吧。我这边实数部分和虚数部分的处理是采用求模的方式。有问题的，请广大人民群众提出来啊。

三、用gabor提取纹理特征的思路【抄别人的论文】

  Gabor滤波方法的主要思想是：不同纹理一般具有不同的中心频率及带宽，根据这些频率和带宽可以设计一组Gabor滤波器对纹理图像进行滤波，每个Gabor滤波器只允许与其频率相对应的纹理顺利通过，而使其他纹理的能量受到抑制，从各滤波器的输出结果中分析和提取纹理特征，用于之后的分类或分割任务。Gabor滤波器提取纹理特征主要包括两个过程：①设计滤波器(例如函数、数目、方向和间隔)；②从滤波器的输出结果中提取有效纹理特征集。Gabor滤波器是带通滤波器，它的单位冲激响应函数(Gabor函数)是高斯函数与复指数函的乘积。它是达到时频测不准关系下界的函数，具有最好地兼顾信号在时频域的分辨能力。

  实现步骤：

（1）将输入图像分为3×3（9块）和4×4（16块）的图像块；

（2）建立Gabor滤波器组：选择4个尺度，6个方向，这样组成了24个Gabor滤波器；

（3）Gabor滤波器组与每个图像块在空域卷积，每个图像块可以得到24个滤波器输出，这      些输出是图像块大小的图像，如果直接将其作为特征向量，特征空间的维数会很大，      所以需要“浓缩”；

（4）每个图像块经过Gabor滤波器组的24个输出，要“浓缩”（文中提到“average filter        responses within the block”我的理解是取灰度均值）为一个24×1的列向量作为该图像       块的纹理特征。查阅相关文献，发现也可以用方差。

    利用一幅真实图像，按照文献原文所说，利用4scales*6orientations的Gabor滤波器组进行纹理特征提取，可以有效获得图像纹理信息。其中，单独拿出某组相同scale的结果，展示如下所示。【别人的实验结果，也没给代码，我也没去做】

本文转自 sunlinju

在数字图像处理领域，Gabor滤波器是以Dennis Gabor命名的，Gabor滤波器是用作边缘检测的线性滤波器。Gabor滤波器的频率和方向的表达与人类的视觉系统很相似。研究发现，Gabor滤波器非常适合纹理表达和分离。在空间域中，一个二维Gabor滤波器是一个由正弦平面波调制的高斯核函数。  

Figure 1. Example of a two-dimensional Gabor filter 

Gabor滤波器的脉冲响应，可以定义为一个正弦波（对于二维Gabor滤波器是正弦平面波）乘以高斯函数。由于乘法卷积性质，Gabor滤波器的脉冲响应的傅立叶变换是其调和函数的傅立叶变换和高斯函数傅立叶变换的卷积。该滤波器由实部和虚部组成，二者相互正交。 

复数表达： 

g(x,y;λ,θ,ψ,σ,γ)=exp(−x2+γ2y′22σ2)exp(i(2πx′λ+ψ))

实数部分： 

g(x,y;λ,θ,ψ,σ,γ)=exp(−x2+γ2y′22σ2)cos((2πx′λ+ψ))

虚数部分： 

g(x,y;λ,θ,ψ,σ,γ)=exp(−x2+γ2y′22σ2)sin((2πx′λ+ψ))

其中: 

x′=xcosθ+ysinθy′=−xsinθ+ycosθ

λ：代表正弦因子的波长，它的值以像素为单位指定，通常大于等于2.但不能大于输入图像尺寸的五分之一； 

θ:指定了Gabor函数并行条纹的方向，它的取值为0到360度 

ψ：代表相位偏移，它的取值范围为-180度到180度。其中，0和180度分别对应中心对称的center-on函数和center-off函数，而-90度和90度对应反对称函数。 

γ:空间纵横比，决定了Gabor函数形状的椭圆率。当γ= 1时，形状是圆的。当γ< 1时，形状随着平行条纹方向而拉长。通常该值为0.5

带宽（b）：Gabor滤波器的半响应空间频率带宽b和σ/ λ的比率有关，其中σ表示Gabor函数的高斯因子的标准差，如下：

参考： 

1.https://www.mathworks.com/matlabcentral/fileexchange/23253-gabor-filter《程序》 

2、http://www.cs.utah.edu/~arul/report/node13.html《说明》 

3、https://en.wikipedia.org/wiki/Gabor_filter《维基》 

4、http://blog.csdn.net/jinshengtao/article/details/17797641《参考的翻译》

本文转自网页

拉普拉斯金字塔

一般来说，由于PyrDown的过程会损失部分图像信息，因此通常情况下：

Gi≠PyrUp(PyrDown(Gi))Gi≠PyrUp(PyrDown(Gi))

为了使PyrUp和PyrDown可逆，这里引入拉普拉斯金字塔的概念。其定义如下：

Li=Gi−PyrUp(Gi+1)=Gi−Up(Gi+1)⊗H5×5Li=Gi−PyrUp(Gi+1)=Gi−Up(Gi+1)⊗H5×5

其中UP操作是将源图像中位置为(x,y)的像素映射到目标图像的(2x+1,2y+1)位置，H5×5H5×5表示5x5的高斯核。

整个拉普拉斯金字塔运算过程可以通过下图来概括：

图中最左列和最右列都是高斯金字塔，中间一列是拉普拉斯金字塔。

注：严格来说，这样的金字塔应该叫做高斯差分金字塔，只不过高斯差分恰好是拉普拉斯算子的近似解，故名。

Steerable金字塔

将拉普拉斯金字塔中的高斯滤波函数，换成Steerable滤波函数即可。

双边滤波（Bilateral filter）

双边滤波（Bilateral filter）是一种可以保边去噪的滤波器。其输出像素的值依赖于邻域像素的值的加权组合，即：

也就是：

其中：

这里的d(i,j,k,l)由于只与定义域有关，通常叫做“定义域核”。实际上，这就是一个高斯滤波核。而r(i,j,k,l)

由于和像素值的差有关（像素差越大，权重越小），也被叫做“值域核”。

从效果来说，双边滤波可产生类似美肤的效果。皮肤上的皱纹和斑，与正常皮肤的差异，远小于黑白眼珠之间的差异，因此前者被平滑，而后者被保留。

为了体现效果，这里来张大叔的照片。

参考：

https://mp.weixin.qq.com/s/JbrZRTwkI7yeC5hOr8g_Mw

Bilateral Solver双边算子原理及应用

Steerable滤波

高斯滤波是一种各向同性滤波，如果想要对特定方向进行滤波的话，可使用Steerable滤波。

=e−(x2+y

这就是两个轴向上的1阶Steerable滤波函数。

任意角度的1阶Steerable滤波函数为：

如果对高斯函数求2阶偏导,还可得到2阶Steerable滤波函数。进一步的讨论详见参考论文。

参考：

1991年IEEE论文：《The Design and Use of Steerable Filters》

作者：William T. Freeman，斯坦福大学本科+斯坦福/康奈尔大学双料硕士+麻省理工学院博士，麻省理工学院教授。1987年，曾做为访问学者在太原理工大学待了一学年。不知道爱不爱吃刀削面(^ω^)

Edward H. Adelson，密歇根大学博士，麻省理工学院教授。

Gabor滤波

基、线性无关、正交

一般的函数可以展开为幂级数或者Fourier级数。这些级数中的幂函数或者正弦函数，被称作“基(basis)函数”。

基的属性主要涉及“线性无关”和“正交”这两个名词。

线性无关的几何含义：在R³

R³(3维空间)中，如果三个向量不共面，则它们相互线性无关。

基如果线性无关，则其函数的级数展开式是唯一的。由于线性相关基使用的比较少，以下如无特指，基均为线性无关基。

正交的几何含义：两个向量正交，则它们是相互垂直的。

正交基一定线性无关，反之则不成立。一般采用施密特正交化方法，将线性无关基，转换为正交基。

幂级数是线性无关基，而Fourier级数是正交基。

Gabor wavelet

除了以上两种常用的基函数外，其他函数也可以作为基函数。其中使用最多的基函数是小波(wavelet)函数，其变换也被称作小波变换。

需要指出的是，小波函数不是一个函数，而是一类函数。Gabor函数就是小波函数的其中一种，其定义如下：

这里的^a,b a,b为常数，^gg为L^2（R）L²（R）（立方可积函数），且^||g||=1  ||g||=1。

注：Dennis Gabor（1900~1979），全息学创始人，1971年获诺贝尔物理学奖，著有《Theory of Communication》（1946）。

当g为高斯函数时，可得到Gabor wavelet：

Gabor wavelet的性质：

1.Gabor wavelet的Fourier变换还是Gabor wavelet：

2.从物理上来说，Gabor wavelet等效于在一个正弦载波（频域）上，调制一个高斯函数（时空域）。这也是Dennis Gabor最早提出它的时候的用途。

3.Fourier变换是信号在整个时域内的积分，因此反映的是信号频率的统计特性，没有局部化分析信号的功能。而Gabor变换是一种短时Fourier变换，具有良好的时频局部化特性，即非常容易地调整Gabor滤波器的方向、基频带宽及中心频率，从而能够最好的兼顾信号在时空域和频域中的分辨能力。

Gabor filter

将Gabor wavelet扩展到2维，可得到Gabor filter(图像实际上就是一种2维信号)：

其中，

λ ：正弦函数波长；θ

θ：Gabor核函数的方向；ψ

ψ：相位偏移；σ

σ：高斯函数的标准差；γ

γ： 空间的宽高比。

可以看出Gabor filter是一个复函数，其实部为：

其虚部为：

此外，还有对数Gabor函数：

Gabor滤波的效果

参考文献3，给出了Gabor滤波的效果图，如下所示：

图1

从效果来看，该滤波可获得美术上的浮雕效果。但实际上，大多数的边缘检测算法都可得到类似效果，这并不是Gabor滤波的主流用法。

以下对参考文献3做一个补充说明：

1.Gabor滤波是复数域的，这点和之前提到的滤波算法有很大的不同。因此，Gabor滤波计算核的方法有3种：复数、实部和虚部。参考文献3采用的是实部法。1987年，J.P. Jones和L.A. Palmer发现Gabor变换所采用的核（Kernels）的实部与哺乳动物视觉皮层简单细胞2D感受野剖面（Profile）非常相似。

2.实部计算的结果有正有负。参考文献3给出的归一化算法，很有通用性，摘录如下：

其中，F为源图像所有像素的集合，D为总的灰度级数。

Gabor滤波采样方式与图像压缩

Gabor滤波和之前的滤波算法的另一大差异是：Gabor滤波核不是一个，而是由若干不同参数组合而成的一组核，其中的每一个参数组合被称为一个采样点。

从Gabor filter的计算公式亦可看出，组成采样点的参数，既有时空域参数，也有频域参数。这些采样点在时空域和频域中如何分布，才能达到最终效果呢？

由于Gabor filter不是正交基，因此针对采样点分布提出了Tight Frame的概念。参考文献1给出了满足Tight Frame要求的采样点分布方式（简称采样方式）的条件。这里的推导非常复杂，但从概念上可以类比信号处理中的奈奎斯特采样定理。

Tight Frame有个重要特性：

如果采样方式满足Tight Frame条件，且 图像集B=Gabor(图像A),图像集C=Gabor^—1(图像B)，那么图像A≈图像CC，其中的Gabor$Gabor$和Gabor−1$Gabo{r}^{-1}$分别表示Gabor变换及其逆变换。

参考文献2给出了采用上述方法对Lena图进行压缩并还原的例子。这也是Gabor滤波在图像处理领域的早期典型应用。

Gabor滤波与模式识别

图2

2000年以后，科学界对Gabor滤波的研究，主要集中在模式识别方面。比如图2就是参考文献4中给出的人脸识别方面的Gabor滤波效果图。其中，左边是原图，而右边是40组不同参数的Gabor滤波器所得到的滤波效果图。

注：1幅原图变成40幅滤波效果图的过程，在数学上是个升维过程。在后处理阶段为了处理的方便，往往会进行数据降维，如参考文献5所示。

从中还可以看出，虽然图1显示出一定的艺术处理效果，但大多数情况下，Gabor滤波所得的图像是如图2所示的极度扭曲而无明显意义的图片。Gabor滤波的真正用途，并不是给人看，而是给机器看。

从上面的讨论可知，Gabor滤波是一种带通滤波，使用不同的时空域或频域参数，可以过滤出不同的时空域或频域特征。这些特征正是模式识别所需要的。

图3

图3是参考文献4给出的一种Gabor滤波器的使用场景图，从中可以看出Gabor滤波效果图是如何应用到人脸识别技术中的。

必须指出的是：Gabor滤波效果图的后处理方法有很多种，而图3仅是其中一种而已。

参考

1.1996年IEEE论文：《Image Representation Using 2D Gabor Wavelets》

作者：Tai Sing Lee，哈佛大学博士，卡内基梅隆大学教授。

2.1988年IEEE论文：《Complete Discrete 2-D Gabor Transforms by Neural Networks for Image Analysis and Compression》

作者：JohnG. Daugman，哈佛大学博士，剑桥大学教授。

3.http://blog.csdn.net/xiaowei_cqu/article/details/24745945

Gabor滤波器

4.《Face recognition using Ada-Boosted Gabor features》

作者：Peng Yang，Shiguang Shan，Wen Gao，Stan Z. Li，Dong Zhang，中科院计算所和微软亚洲研究院的几个小牛。

高文，1956年生，哈尔滨工业大学博士（1988）+东京大学博士（1991）。中科院计算所所长，中国工程院院士。

山世光，1975年生，哈尔滨工业大学本硕（1997，1999）+中科院博士（2004）。中科视拓CEO。

5.http://www.cnblogs.com/Jack-Lee/p/3649114.html

基于OpenCV的Gabor变换及特征提取

https://mp.weixin.qq.com/s/RS0Zw6iLFQBUlpOdkuoFbg

从傅立叶变换到Gabor滤波器

https://mp.weixin.qq.com/s/UDJAYB79KaPB0pgnnCIPTw

不确定性原理的前世今生