基於OpenCV的Gabor變換及特徵提取
阿新 • • 發佈：2018-11-11
一、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函式的第一種表示形式：
二維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)：
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/[email protected]/blog/static/171861983201172091718341/
參考文獻：基於Gabor變換的特徵提取及其應用