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

一、Gabor变换概述

 　　Gabor变换是一种加窗短时Fourier变换（Window Fourier transform or Short Time Fourier Transform）。Fourier变换是整体上将信号分解为不同的频率分量（任何信号都可分解为复正弦信号之和），对确定性信号及平稳信号使用。其缺点为缺乏时间的局部性信息，并且对时变信号、非平稳信号的分析存在严重不足，(1)无法告知某些频率成分发生在哪些时间内；(2)无法告知某个时刻信号频谱的分布情况。

　　Gabor函数可以在频域不同尺度、不同方向上提取相关的特征。另外Gabor函数与人眼的生物作用相仿，所以经常用作纹理识别上，并取得了较好的效果，Gabor变换是短时Fourier变换中当窗函数取为高斯函数时的一种特殊情况。Gabor变换的本质实际上还是对二维图像求卷积。因此二维卷积运算的效率就直接决定了Gabor变换的效率。

　　下面简要说一下二维卷积运算的计算过程：

　　A可以理解成是待处理的笔迹纹理，B可以理解成Gabor变换的核函数，现在要求A与B的离散二维叠加卷积，我们首先对A的右边界和下边界填充0（zero padding)，然后将B进行水平翻转和垂直翻转，如下图：

　　

然后用B中的每个值依次乘以A中相对位置处的值并进行累加，结果填入相应位置处（注意红圈位置）。通常二维卷积的结果比A、B的尺寸要大，如下图：

二、Gabor函数表达式

　　通过频率参数和高斯函数参数的选取，Gabor变换可以选取很多纹理特征，但是Gabor是非正交的，不同特征分量之间有冗余。 Gabor是有高斯核函数与复正弦函数调制而成，如下图所示：

   图(a)为偏移x轴30°的正弦函数，图(b)为高斯核，图(c)为对应的Gabor filter。可以看出正弦函数是如何在空间上具有局部性的。

        二维Gabor函数的第一种表示形式：

 

       为正弦函数的波长，核函数方向，为相位偏移，为高斯标准差，为x、y两个方向的纵横比（指定了Gabor函数的椭圆率）。

       二维Gabor函数的第二种形式：

 

   v的取值决定了Gabor滤波的波长，u的取值表示Gabor核函数的方向，K表示总的方向数。参数决定了高斯窗口的大小，这里取

三、Gabor特征提取

　　提取Gabor特征之前，我们可以先对需要处理的图像I(x,y)进行实数形式的Gabor变换，得到处理后的图像。举个例子，我们需要处理的图像大小为128x128，可以先通过变换得到处理后的图像，大小也是128x128.但是我们不方便直接提取特征，因为这样提取出来的特征维数太高，不利于我们进行后续处理。这里我们对图像分块，分别水平和垂直方向取16等分，这样就可以将整个图像分成64个16x16大小的子图像块。如下图所示:

 

接着计算每一块对应的能量，第k块的能量定义如下：

 

这样计算之后就可以形成如下图所示的联合空间频率能量矩阵Energy。

 

然后，我们将能量矩阵降维成1x64的行向量，作为我们的原始图像在某一方向和尺度变换后的特征向量，如下：

 

接着就可以利用得到的特征向量进行下一步的处理，比如相似度判别或者聚类等等。当然一般情况下，需要提取多个方向和尺度特征，变换参数重复上述过程就可以了。

四、基于OpenCV 2.x的实现

　　下面的代码是根据上面Gabor变换的第二种形式来的，原始代码来源于这里，由于原始代码使用c语言实现，里面很多函数是基于OpenCV1.x的，比如使用了大量的IplImage类型，这种是比较老的图片处理类型，而且要手动进行内存管理，加之不方便进行矩阵运算，所以我在代码的基础上进行了修改，使用了Mat类型代替了原始的函数参数，并对原函数中的一些不必要的函数进行了精简。整个函数基于C++实现，已经调试通过。

　　头文件(cvgabor.h)：

#ifndef CVGABOR_H
#define CVGABOR_H

#include <iostream>


#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>

const double PI = 3.14159265;
const int CV_GABOR_REAL = 1;
const int CV_GABOR_IMAG = 2;
const int CV_GABOR_MAG  = 3;
const int CV_GABOR_PHASE = 4;

using namespace cv;
using namespace std;

class CvGabor{
public:
    CvGabor();
    ~CvGabor();

    CvGabor(int iMu, int iNu);
    CvGabor(int iMu, int iNu, double dSigma);
    CvGabor(int iMu, int iNu, double dSigma, double dF);
    CvGabor(double dPhi, int iNu);
    CvGabor(double dPhi, int iNu, double dSigma);
    CvGabor(double dPhi, int iNu, double dSigma, double dF);

    void Init(int iMu, int iNu, double dSigma, double dF);
    void Init(double dPhi, int iNu, double dSigma, double dF);

    bool IsInit();
    bool IsKernelCreate();
    int mask_width();
    int get_mask_width();

    void get_image(int Type, Mat& image);
    void get_matrix(int Type, Mat& matrix);

    void conv_img(Mat& src, Mat& dst, int Type);

protected:
    double Sigma;
    double F;
    double Kmax;
    double K;
    double Phi;
    bool bInitialised;
    bool bKernel;
    int Width;
    Mat Imag;
    Mat Real;
  
private:
    void creat_kernel();
    
};

#endif

　　函数实现(cvgabor.cpp)：

#include "cvgabor.h"

CvGabor::CvGabor()
{
}


CvGabor::~CvGabor()
{
}


/*!

Parameters:
iMu        The orientation iMu*PI/8,
iNu         The scale,
dSigma         The sigma value of Gabor,
dPhi        The orientation in arc
dF        The spatial frequency

*/

CvGabor::CvGabor(int iMu, int iNu)
{
    double dSigma = 2*PI;
    F = sqrt(2.0);
    Init(iMu, iNu, dSigma, F);
}

CvGabor::CvGabor(int iMu, int iNu, double dSigma)
{ 
    F = sqrt(2.0);
    Init(iMu, iNu, dSigma, F);
}

CvGabor::CvGabor(int iMu, int iNu, double dSigma, double dF)
{
    Init(iMu, iNu, dSigma, dF);
}

CvGabor::CvGabor(double dPhi, int iNu)
{
    Sigma = 2*PI;
    F = sqrt(2.0);
    Init(dPhi, iNu, Sigma, F);
}

CvGabor::CvGabor(double dPhi, int iNu, double dSigma)
{
    F = sqrt(2.0);
    Init(dPhi, iNu, dSigma, F);
}

CvGabor::CvGabor(double dPhi, int iNu, double dSigma, double dF)
{
    Init(dPhi, iNu, dSigma,dF);
}

/*!
Parameters:
iMu     The orientations which is iMu*PI.8
iNu     The scale can be from -5 to infinit
dSigma     The Sigma value of gabor, Normally set to 2*PI
dF     The spatial frequence , normally is sqrt(2)

Initilize the.gabor with the orientation iMu, the scale iNu, the sigma dSigma, the frequency dF, it will call the function creat_kernel(); So a gabor is created.
*/
void CvGabor::Init(int iMu, int iNu, double dSigma, double dF)
{
    //Initilise the parameters
    bInitialised = false;
    bKernel = false;

    Sigma = dSigma;
    F = dF;

    Kmax = PI/2;

    //Absolute value of K
    K = Kmax / pow(F, (double)iNu);
    Phi = PI*iMu/8;
    bInitialised = true;

    Width = mask_width();
    creat_kernel();
}

/*!
Parameters:
dPhi     The orientations
iNu     The scale can be from -5 to infinit
dSigma     The Sigma value of gabor, Normally set to 2*PI
dF     The spatial frequence , normally is sqrt(2)

Initilize the.gabor with the orientation dPhi, the scale iNu, the sigma dSigma, the frequency dF, it will call the function creat_kernel(); So a gabor is created.filename     The name of the image file
file_format     The format of the file,
*/
void CvGabor::Init(double dPhi, int iNu, double dSigma, double dF)
{

    bInitialised = false;
    bKernel = false;
    Sigma = dSigma;
    F = dF;

    Kmax = PI/2;

    // Absolute value of K
    K = Kmax / pow(F, (double)iNu);
    Phi = dPhi;
    bInitialised = true;

    Width = mask_width();
    creat_kernel();
}

/*!
Returns:
a boolean value, TRUE is created or FALSE is non-created.

Determine whether a gabor kernel is created.
*/

bool CvGabor::IsInit()
{
    return bInitialised;
}

bool CvGabor::IsKernelCreate()
{
    return bKernel;
}

/*!
Returns:
The long type show the width.

Return the width of mask (should be NxN) by the value of Sigma and iNu.
*/
int CvGabor::mask_width()
{
    int lWidth;
    if (IsInit() == false)
    {
        cerr << "Error: The Object has not been initilised in mask_width()!
" << endl;
        return 0;
    }
    else
    {
        //determine the width of Mask
        double dModSigma = Sigma/K;
        int dWidth = cvRound(dModSigma*6 + 1);

        //test whether dWidth is an odd.
        if((dWidth % 2) == 0)
        {
            lWidth = dWidth + 1;
        }
        else
        {
            lWidth = dWidth;
        }
        return lWidth;
    }
}

/*!

Returns:
Pointer to long type width of mask.

*/
int CvGabor::get_mask_width()
{
    return Width;
}

/*!
fn CvGabor::creat_kernel()
Create gabor kernel

Create 2 gabor kernels - REAL and IMAG, with an orientation and a scale 
*/
void CvGabor::creat_kernel()
{

    if (IsInit() == false)
    {
        cerr << "Error: The Object has not been initilised in creat_kernel()!" << endl;
    }
    else
    {
        Mat mReal(Width, Width, CV_32FC1);
        Mat mImag(Width, Width, CV_32FC1);

        /**************************** Gabor Function ****************************/ 
        int x, y;
        double dReal;
        double dImag;
        double dTemp1, dTemp2, dTemp3;

        for (int i = 0; i < Width; i++)
        {
            for (int j = 0; j < Width; j++)
            {
                x = i-(Width-1)/2;
                y = j-(Width-1)/2;
                dTemp1 = (pow(K,2)/pow(Sigma,2))*exp(-(pow((double)x,2)+pow((double)y,2))*pow(K,2)/(2*pow(Sigma,2)));
                dTemp2 = cos(K*cos(Phi)*x + K*sin(Phi)*y) - exp(-(pow(Sigma,2)/2));
                dTemp3 = sin(K*cos(Phi)*x + K*sin(Phi)*y);
                dReal = dTemp1*dTemp2;
                dImag = dTemp1*dTemp3;

                mReal.row(i).col(j) = dReal;
                mImag.row(i).col(j) = dImag;
            }
        }
        /**************************** Gabor Function ****************************/
        bKernel = true;

        mReal.copyTo(Real);
        mImag.copyTo(Imag);
    }
}


/*!
fn CvGabor::get_image(int Type)
Get the speific type of image of Gabor

Parameters:
Type        The Type of gabor kernel, e.g. REAL, IMAG, MAG, PHASE   

Returns:
Pointer to image structure, or NULL on failure    

Return an Image (gandalf image class) with a specific Type   "REAL"    "IMAG" "MAG" "PHASE"  
*/
void CvGabor::get_image(int Type, Mat& image)
{
    if(IsKernelCreate() == false)
    { 
        cerr << "Error: the Gabor kernel has not been created in get_image()!" << endl;
        return;
    }
    else
    {  
        Mat re(Width, Width, CV_32FC1);
        Mat im(Width, Width, CV_32FC1);
        Mat temp;

        switch(Type)
        {
        case 1:  //Real
            temp = Real.clone();
            normalize(temp, temp, 255.0, 0.0, NORM_MINMAX);
            break;
        case 2:  //Imag
            temp = Imag.clone();
            break; 
        case 3:  //Magnitude
            re = Real.clone();
            im = Imag.clone();

            pow(re, 2, re);
            pow(im, 2, im);
            add(im, re, temp);
            pow(temp, 0.5, temp);
            break;
        case 4:  //Phase
            ///@todo
            break;
        }

        convertScaleAbs(temp, image, 1, 0);
    }
}

/*!
fn CvGabor::get_matrix(int Type)
Get a matrix by the type of kernel

Parameters:
Type        The type of kernel, e.g. REAL, IMAG, MAG, PHASE

Returns:
Pointer to matrix structure, or NULL on failure.

Return the gabor kernel.
*/
void CvGabor::get_matrix(int Type, Mat& matrix)
{
    if (!IsKernelCreate())
    {
        cerr << "Error: the gabor kernel has not been created!" << endl;
        return;
    }
    switch (Type)
    {
    case CV_GABOR_REAL:
        matrix = Real.clone();
        break;
    case CV_GABOR_IMAG:
        matrix = Imag.clone();
        break;
    case CV_GABOR_MAG:
        break;
    case CV_GABOR_PHASE:
        break;
    }
}

/*!
fn CvGabor::conv_img_a(IplImage *src, IplImage *dst, int Type)
*/
void CvGabor::conv_img(Mat &src, Mat &dst, int Type)
{
    Mat mat = src.clone();

    Mat rmat(src.rows, src.cols, CV_32FC1);
    Mat imat(src.rows, src.cols, CV_32FC1);

    switch (Type)
    {
    case CV_GABOR_REAL:
        filter2D(mat, mat, 1, Real, Point((Width-1)/2, (Width-1)/2));
        break;

    case CV_GABOR_IMAG:
        filter2D(mat, mat, 1, Imag, Point( (Width-1)/2, (Width-1)/2));
        break;

    case CV_GABOR_MAG:
        /* Real Response */
        filter2D(mat, rmat, 1, Real, Point((Width-1)/2, (Width-1)/2));

        /* Imag Response */
        filter2D(mat, imat, 1, Imag, Point( (Width-1)/2, (Width-1)/2));

        /* Magnitude response is the square root of the sum of the square of real response and imaginary response */
        pow(rmat, 2, rmat);
        pow(imat, 2, imat);
        add(rmat, imat, mat);
        pow(mat, 0.5, mat);
        break;

    case CV_GABOR_PHASE:
        break;
    }

//    cvNormalize(mat, mat, 0, 255, CV_MINMAX, NULL);
    mat.copyTo(dst);
}

 函数的使用方法如下：

首先显示核函数图像：

    //创建一个方向是PI/4而尺度是3的gabor
    double Sigma = 2*PI;
    double F = sqrt(2.0);
    CvGabor gabor(PI/4, 3, Sigma, F);

    //获得实部并显示它
    Mat kernel(gabor.get_mask_width(), gabor.get_mask_width(), CV_8UC1);
    gabor.get_image(CV_GABOR_REAL, kernel);
    imshow("Kernel", kernel);
    cout << kernel.rows << endl;
    cout << kernel.cols << endl;
    cvWaitKey(0);

 显示效果如下：

接着，对输入的图像进行处理：

    //载入一个图像并显示
    Mat img = imread( "test.jpg", CV_LOAD_IMAGE_GRAYSCALE );
    imshow("Original Image", img);
    cvWaitKey(0);

    //获取载入图像的gabor滤波响应的实部并且显示
    Mat reimg(img.rows,img.cols, CV_32FC1);
    gabor.conv_img(img, reimg, CV_GABOR_REAL);
    imshow("After Image", reimg);
    cvWaitKey(0);

下面是显示效果:

接着就可以对得到的图像，按照上面体征提取部分的步骤进行特征提取。

参考链接：http://blog.csdn.net/renjinr/article/details/13768655

　　　　　http://blog.sina.com.cn/s/blog_75e063c10101455s.html

　　　　　http://blog.163.com/huai_jing@126/blog/static/171861983201172091718341/

参考文献：基于Gabor变换的特征提取及其应用

本文为原创内容，转载请注明出处！http://www.cnblogs.com/Jack-Lee/p/3649114.html