基於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()!\n" << 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