基于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 <wbr>变换（1）决定了高斯窗口的大小，这里取Gabor <wbr>变换（1）

三、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)：

 1 #ifndef CVGABOR_H
 2 #define CVGABOR_H
 3 
 4 #include <iostream>
 5 
 6 
 7 #include <opencv2/core/core.hpp>
 8 #include <opencv2/highgui/highgui.hpp>
 9 #include <opencv2/imgproc/imgproc.hpp>
10 
11 const double PI = 3.14159265;
12 const int CV_GABOR_REAL = 1;
13 const int CV_GABOR_IMAG = 2;
14 const int CV_GABOR_MAG  = 3;
15 const int CV_GABOR_PHASE = 4;
16 
17 using namespace cv;
18 using namespace std;
19 
20 class CvGabor{
21 public:
22     CvGabor();
23     ~CvGabor();
24 
25     CvGabor(int iMu, int iNu);
26     CvGabor(int iMu, int iNu, double dSigma);
27     CvGabor(int iMu, int iNu, double dSigma, double dF);
28     CvGabor(double dPhi, int iNu);
29     CvGabor(double dPhi, int iNu, double dSigma);
30     CvGabor(double dPhi, int iNu, double dSigma, double dF);
31 
32     void Init(int iMu, int iNu, double dSigma, double dF);
33     void Init(double dPhi, int iNu, double dSigma, double dF);
34 
35     bool IsInit();
36     bool IsKernelCreate();
37     int mask_width();
38     int get_mask_width();
39 
40     void get_image(int Type, Mat& image);
41     void get_matrix(int Type, Mat& matrix);
42 
43     void conv_img(Mat& src, Mat& dst, int Type);
44 
45 protected:
46     double Sigma;
47     double F;
48     double Kmax;
49     double K;
50     double Phi;
51     bool bInitialised;
52     bool bKernel;
53     int Width;
54     Mat Imag;
55     Mat Real;
56   
57 private:
58     void creat_kernel();
59     
60 };
61 
62 #endif

 

　　函数实现(cvgabor.cpp)：

  1 #include "cvgabor.h"
  2 
  3 CvGabor::CvGabor()
  4 {
  5 }
  6 
  7 
  8 CvGabor::~CvGabor()
  9 {
 10 }
 11 
 12 
 13 /*!
 14 
 15 Parameters:
 16 iMu        The orientation iMu*PI/8,
 17 iNu         The scale,
 18 dSigma         The sigma value of Gabor,
 19 dPhi        The orientation in arc
 20 dF        The spatial frequency
 21 
 22 */
 23 
 24 CvGabor::CvGabor(int iMu, int iNu)
 25 {
 26     double dSigma = 2*PI;
 27     F = sqrt(2.0);
 28     Init(iMu, iNu, dSigma, F);
 29 }
 30 
 31 CvGabor::CvGabor(int iMu, int iNu, double dSigma)
 32 { 
 33     F = sqrt(2.0);
 34     Init(iMu, iNu, dSigma, F);
 35 }
 36 
 37 CvGabor::CvGabor(int iMu, int iNu, double dSigma, double dF)
 38 {
 39     Init(iMu, iNu, dSigma, dF);
 40 }
 41 
 42 CvGabor::CvGabor(double dPhi, int iNu)
 43 {
 44     Sigma = 2*PI;
 45     F = sqrt(2.0);
 46     Init(dPhi, iNu, Sigma, F);
 47 }
 48 
 49 CvGabor::CvGabor(double dPhi, int iNu, double dSigma)
 50 {
 51     F = sqrt(2.0);
 52     Init(dPhi, iNu, dSigma, F);
 53 }
 54 
 55 CvGabor::CvGabor(double dPhi, int iNu, double dSigma, double dF)
 56 {
 57     Init(dPhi, iNu, dSigma,dF);
 58 }
 59 
 60 /*!
 61 Parameters:
 62 iMu     The orientations which is iMu*PI.8
 63 iNu     The scale can be from -5 to infinit
 64 dSigma     The Sigma value of gabor, Normally set to 2*PI
 65 dF     The spatial frequence , normally is sqrt(2)
 66 
 67 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.
 68 */
 69 void CvGabor::Init(int iMu, int iNu, double dSigma, double dF)
 70 {
 71     //Initilise the parameters
 72     bInitialised = false;
 73     bKernel = false;
 74 
 75     Sigma = dSigma;
 76     F = dF;
 77 
 78     Kmax = PI/2;
 79 
 80     //Absolute value of K
 81     K = Kmax / pow(F, (double)iNu);
 82     Phi = PI*iMu/8;
 83     bInitialised = true;
 84 
 85     Width = mask_width();
 86     creat_kernel();
 87 }
 88 
 89 /*!
 90 Parameters:
 91 dPhi     The orientations
 92 iNu     The scale can be from -5 to infinit
 93 dSigma     The Sigma value of gabor, Normally set to 2*PI
 94 dF     The spatial frequence , normally is sqrt(2)
 95 
 96 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
 97 file_format     The format of the file,
 98 */
 99 void CvGabor::Init(double dPhi, int iNu, double dSigma, double dF)
100 {
101 
102     bInitialised = false;
103     bKernel = false;
104     Sigma = dSigma;
105     F = dF;
106 
107     Kmax = PI/2;
108 
109     // Absolute value of K
110     K = Kmax / pow(F, (double)iNu);
111     Phi = dPhi;
112     bInitialised = true;
113 
114     Width = mask_width();
115     creat_kernel();
116 }
117 
118 /*!
119 Returns:
120 a boolean value, TRUE is created or FALSE is non-created.
121 
122 Determine whether a gabor kernel is created.
123 */
124 
125 bool CvGabor::IsInit()
126 {
127     return bInitialised;
128 }
129 
130 bool CvGabor::IsKernelCreate()
131 {
132     return bKernel;
133 }
134 
135 /*!
136 Returns:
137 The long type show the width.
138 
139 Return the width of mask (should be NxN) by the value of Sigma and iNu.
140 */
141 int CvGabor::mask_width()
142 {
143     int lWidth;
144     if (IsInit() == false)
145     {
146         cerr << "Error: The Object has not been initilised in mask_width()!\n" << endl;
147         return 0;
148     }
149     else
150     {
151         //determine the width of Mask
152         double dModSigma = Sigma/K;
153         int dWidth = cvRound(dModSigma*6 + 1);
154 
155         //test whether dWidth is an odd.
156         if((dWidth % 2) == 0)
157         {
158             lWidth = dWidth + 1;
159         }
160         else
161         {
162             lWidth = dWidth;
163         }
164         return lWidth;
165     }
166 }
167 
168 /*!
169 
170 Returns:
171 Pointer to long type width of mask.
172 
173 */
174 int CvGabor::get_mask_width()
175 {
176     return Width;
177 }
178 
179 /*!
180 \fn CvGabor::creat_kernel()
181 Create gabor kernel
182 
183 Create 2 gabor kernels - REAL and IMAG, with an orientation and a scale 
184 */
185 void CvGabor::creat_kernel()
186 {
187 
188     if (IsInit() == false)
189     {
190         cerr << "Error: The Object has not been initilised in creat_kernel()!" << endl;
191     }
192     else
193     {
194         Mat mReal(Width, Width, CV_32FC1);
195         Mat mImag(Width, Width, CV_32FC1);
196 
197         /**************************** Gabor Function ****************************/ 
198         int x, y;
199         double dReal;
200         double dImag;
201         double dTemp1, dTemp2, dTemp3;
202 
203         for (int i = 0; i < Width; i++)
204         {
205             for (int j = 0; j < Width; j++)
206             {
207                 x = i-(Width-1)/2;
208                 y = j-(Width-1)/2;
209                 dTemp1 = (pow(K,2)/pow(Sigma,2))*exp(-(pow((double)x,2)+pow((double)y,2))*pow(K,2)/(2*pow(Sigma,2)));
210                 dTemp2 = cos(K*cos(Phi)*x + K*sin(Phi)*y) - exp(-(pow(Sigma,2)/2));
211                 dTemp3 = sin(K*cos(Phi)*x + K*sin(Phi)*y);
212                 dReal = dTemp1*dTemp2;
213                 dImag = dTemp1*dTemp3;
214 
215                 mReal.row(i).col(j) = dReal;
216                 mImag.row(i).col(j) = dImag;
217             }
218         }
219         /**************************** Gabor Function ****************************/
220         bKernel = true;
221 
222         mReal.copyTo(Real);
223         mImag.copyTo(Imag);
224     }
225 }
226 
227 
228 /*!
229 \fn CvGabor::get_image(int Type)
230 Get the speific type of image of Gabor
231 
232 Parameters:
233 Type        The Type of gabor kernel, e.g. REAL, IMAG, MAG, PHASE   
234 
235 Returns:
236 Pointer to image structure, or NULL on failure    
237 
238 Return an Image (gandalf image class) with a specific Type   "REAL"    "IMAG" "MAG" "PHASE"  
239 */
240 void CvGabor::get_image(int Type, Mat& image)
241 {
242     if(IsKernelCreate() == false)
243     { 
244         cerr << "Error: the Gabor kernel has not been created in get_image()!" << endl;
245         return;
246     }
247     else
248     {  
249         Mat re(Width, Width, CV_32FC1);
250         Mat im(Width, Width, CV_32FC1);
251         Mat temp;
252 
253         switch(Type)
254         {
255         case 1:  //Real
256             temp = Real.clone();
257             normalize(temp, temp, 255.0, 0.0, NORM_MINMAX);
258             break;
259         case 2:  //Imag
260             temp = Imag.clone();
261             break; 
262         case 3:  //Magnitude
263             re = Real.clone();
264             im = Imag.clone();
265 
266             pow(re, 2, re);
267             pow(im, 2, im);
268             add(im, re, temp);
269             pow(temp, 0.5, temp);
270             break;
271         case 4:  //Phase
272             ///@todo
273             break;
274         }
275 
276         convertScaleAbs(temp, image, 1, 0);
277     }
278 }
279 
280 /*!
281 \fn CvGabor::get_matrix(int Type)
282 Get a matrix by the type of kernel
283 
284 Parameters:
285 Type        The type of kernel, e.g. REAL, IMAG, MAG, PHASE
286 
287 Returns:
288 Pointer to matrix structure, or NULL on failure.
289 
290 Return the gabor kernel.
291 */
292 void CvGabor::get_matrix(int Type, Mat& matrix)
293 {
294     if (!IsKernelCreate())
295     {
296         cerr << "Error: the gabor kernel has not been created!" << endl;
297         return;
298     }
299     switch (Type)
300     {
301     case CV_GABOR_REAL:
302         matrix = Real.clone();
303         break;
304     case CV_GABOR_IMAG:
305         matrix = Imag.clone();
306         break;
307     case CV_GABOR_MAG:
308         break;
309     case CV_GABOR_PHASE:
310         break;
311     }
312 }
313 
314 /*!
315 \fn CvGabor::conv_img_a(IplImage *src, IplImage *dst, int Type)
316 */
317 void CvGabor::conv_img(Mat &src, Mat &dst, int Type)
318 {
319     Mat mat = src.clone();
320 
321     Mat rmat(src.rows, src.cols, CV_32FC1);
322     Mat imat(src.rows, src.cols, CV_32FC1);
323 
324     switch (Type)
325     {
326     case CV_GABOR_REAL:
327         filter2D(mat, mat, 1, Real, Point((Width-1)/2, (Width-1)/2));
328         break;
329 
330     case CV_GABOR_IMAG:
331         filter2D(mat, mat, 1, Imag, Point( (Width-1)/2, (Width-1)/2));
332         break;
333 
334     case CV_GABOR_MAG:
335         /* Real Response */
336         filter2D(mat, rmat, 1, Real, Point((Width-1)/2, (Width-1)/2));
337 
338         /* Imag Response */
339         filter2D(mat, imat, 1, Imag, Point( (Width-1)/2, (Width-1)/2));
340 
341         /* Magnitude response is the square root of the sum of the square of real response and imaginary response */
342         pow(rmat, 2, rmat);
343         pow(imat, 2, imat);
344         add(rmat, imat, mat);
345         pow(mat, 0.5, mat);
346         break;
347 
348     case CV_GABOR_PHASE:
349         break;
350     }
351 
352 //    cvNormalize(mat, mat, 0, 255, CV_MINMAX, NULL);
353     mat.copyTo(dst);
354 }

 

 函数的使用方法如下：

首先显示核函数图像：

 1     //创建一个方向是PI/4而尺度是3的gabor
 2     double Sigma = 2*PI;
 3     double F = sqrt(2.0);
 4     CvGabor gabor(PI/4, 3, Sigma, F);
 5 
 6     //获得实部并显示它
 7     Mat kernel(gabor.get_mask_width(), gabor.get_mask_width(), CV_8UC1);
 8     gabor.get_image(CV_GABOR_REAL, kernel);
 9     imshow("Kernel", kernel);
10     cout << kernel.rows << endl;
11     cout << kernel.cols << endl;
12     cvWaitKey(0);

 显示效果如下：

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

 1     //载入一个图像并显示
 2     Mat img = imread( "test.jpg", CV_LOAD_IMAGE_GRAYSCALE );
 3     imshow("Original Image", img);
 4     cvWaitKey(0);
 5 
 6     //获取载入图像的gabor滤波响应的实部并且显示
 7     Mat reimg(img.rows,img.cols, CV_32FC1);
 8     gabor.conv_img(img, reimg, CV_GABOR_REAL);
 9     imshow("After Image", reimg);
10     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

　　　　

posted @ 2014-04-06 20:08  Teckee  阅读(17511)  评论(10)  编辑  收藏  举报