紋理分割(一)Gabor濾波器學習

reference: http://blog.csdn.net/jorg_zhao?viewmode=contents

第一個項目根據公司那邊提供的學習資料，需要用到Gabor濾波器對圖像進行處理，公司那邊關于項目的說法比較商業化，叫X-Ray Image Auto Judging System，之前找了很久論文都沒有思路，用這個英文查找論文，也是不對路，這讓我在前期浪費不少時間，後來查閱大量論文之後，确定關于目前的項目的學術說法應該是“輪胎X射線圖缺陷檢測”，英文是“X-Ray Defect Detection”，目前放暑假，不在學校檢索論文比較麻煩，隻能根據各路大神網友的部落格在做這個項目。

這個項目涉及的主要學術方向是“複雜多紋理圖像提取”，能查到的最普遍的做法是采用Gabor濾波器對圖像進行處理，可能閱讀論文量不到，也許有更好的圖像處理方法和思路，這裡記錄下我目前的做法吧，有效果，但是不是特别好。

看到的大神寫到的部落格，非常詳細，大部分資訊整理自以下部落格：

1. Gabor濾波器詳細介紹：http://mplab.ucsd.edu/tutorials/gabor.pdf

2. Gabor濾波器學習：http://blog.csdn.net/jinshengtao/article/details/17797641

3. 國外大學對Gabor濾波器的參數的講解：http://matlabserver.cs.rug.nl/edgedetectionweb/web/edgedetection_params.html

4. 備份論文，在校外沒法下載下傳：http://www.researchgate.net/publication/252703921_Designing_multiple_Gabor_filters_for_multitexture_image_segmentation

5. 德國某公司的宣傳網頁，看樣子做的非常棒：http://www.cyxplus.fr/produit/our-industy-activities/cyxpert-automatic-defect-detection

一、Gabor Filter 

一維Gabor濾波器

Gabor濾波器是處理一維信号（比如音頻）最佳的帶通濾波器，一個複雜的Gabor濾波器是一個高斯核函數乘以一個複雜的sin函數（A complex Gabor filter is defined as the product of a Gaussian kernel times a complex sinusoid），比如：

k  theta  fo是濾波器的參數。我們可以把這個複雜的Gabor濾波器想象成兩個相位濾波器分成了一個複雜函數的實部和虛部，

二維Gabor濾波器

但是我要用到的是二維Gabor濾波器，也叫空間Gabor濾波器（The Spatial (2-D) Gabor Filter），但是第一條參考中給出的二維Gabor濾波器講解太過于。。怎麼說。

按照普通的解釋來吧，二維Gabor函數的數學表達式：

如何得到的這個公式呢？大神部落格（http://blog.csdn.net/yanmy2012/article/details/8090400）中給出的詳細求解過程，截圖如下：

根據第三條參考文獻的解釋，關于實部和虛部的說明，隻有一句話，實部可以對圖像進行平滑濾波，虛部可以用來邊緣檢測，具體的用法可以可以參考相關論文，不在學校下論文不友善，暫且寫這點吧。

下面具體看一下Gabor濾波器的參數說明：

二、Gabor濾波器的應用和适應性

根據第二條參考文獻給出的說明，可以總結以下幾點：

1. Gabor濾波器可以很好的近似單細胞的感受野細胞（光強刺激下的傳遞函數），在提取目标的局部空間和頻率域資訊方面具有良好的特性。

2. 雖然Gabor小波本身不能構成正交基，但在特定參數下可構成緊架構。Gabor小波對于圖像的邊緣敏感，能夠提供良好的方向選擇和尺度選擇特性，而且對于光照變化不敏感,能夠提供對光照變化良好的适應性。-------Gabor小波被廣泛應用于視覺資訊了解

3. 二維Gabor小波變換是在時頻域進行信号分析處理的重要工具，其變換系數有着良好的視覺特性和生物學背景。------是以被廣泛應用于圖像處理、模式識别等領域。

4. 與傳統的傅立葉變換相比，Gabor小波變換具有良好的時頻局部化特性。即非常容易地調整Gabor濾波器的方向、基頻帶寬及中心頻率進而能夠最好的兼顧信号在時空域和頻域中的分辨能力.

5. Gabor小波變換具有多分辨率特性即變焦能力。即采用多通道濾波技術，将一組具有不同時頻域特性的Gabor小波應用于圖像變換，每個通道都能夠得到輸入圖像的某種局部特性，這樣可以根據需要在不同粗細粒度上分析圖像。

6. 在特征提取方面，Gabor小波變換與其它方法相比：一方面其處理的資料量較少，能滿足系統的實時性要求；另一方面，小波變換對光照變化不敏感，且能容忍一定程度的圖像旋轉和變形，當采用基于歐氏距離進行識别時，特征模式與待測特征不需要嚴格的對應，故能提高系統的魯棒性。

為什麼能夠有這些應用呢？我們需要深入的了解下Gabor的特性。根據大神部落格（http://blog.csdn.net/yanmy2012/article/details/8090400）中的詳細解釋。在基本視覺皮層裡的簡單細胞的感受野局限在很小的空域範圍内，并且高度結構化。

1. Gabor變換所采用的核（Kernels）與哺乳動物視覺皮層簡單細胞2D感受野剖面（Profile）非常相似，具有優良的空間局部性和方向選擇性，能夠抓住圖像局部區域内多個方向的空間頻率（尺度）和局部性結構特征。這樣，Gabor分解可以看作一個對方向和尺度敏感的有方向性的顯微鏡。

2. 二維Gabor函數也類似于增強邊緣以及峰、谷、脊輪廓等底層圖像特征，這相當于增強了被認為是面部關鍵部件的眼睛、鼻子、嘴巴等資訊，同時也增強了諸于黑痣、酒窩、傷疤等局部特征，進而使得在保留總體人臉資訊的同時增強局部特性成為可能。它的小波特性說明了Gabor濾波結果是描述圖像局部灰階分布的有力工具,是以,可以使用Gabor濾波來抽取圖像的紋理資訊。

3.  由于Gabor特征具有良好的空間局部性和方向選擇性，而且對光照、姿态具有一定的魯棒性，是以在人臉識别中獲得了成功的應用。然而，大部分基于Gabor特征的人臉識别算法中，隻應用了Gabor幅值資訊，而沒有應用相位資訊，主要原因是Gabor相位資訊随着空間位置呈周期性變化，而幅值的變化相對平滑而穩定，幅值反映了圖像的能量譜，Gabor幅值特征通常稱為Gabor 能量特征(Gabor Energy Features）。Gabor小波可像放大鏡一樣放大灰階的變化，人臉的一些關鍵功能區域(眼睛、鼻子、嘴、眉毛等)的局部特征被強化，進而有利于區分不同的人臉圖像。

4. Gabor小波核函數具有與哺育動物大腦皮層簡單細胞的二維反射區相同的特性，即具有較強的空間位置和方向選擇性，并且能夠捕捉對應于空間和頻率的局部結構資訊；Gabor濾波器對于圖像的亮度和對比度變化以及人臉姿态變化具有較強的健壯性，并且它表達的是對人臉識别最為有用的局部特征。Gabor 小波是對進階脊椎動物視覺皮層中的神經元的良好逼近，是時域和頻域精确度的一種折中。

三、Gabor濾波器的程式設計實作

這裡我用opencv做的，因為對matlab不是很熟，二是工業應用一般都是C或C++做的，代碼給出吧，我用mfc做的界面，之前有好幾個版本，驗證了幾個算法，最後這個就是為了看效果的，比較簡單。因為用的opencv3.0，是以顯示圖像時遇到了麻煩，查找相關資料，解決辦法是自己建立CvvImage.h CvvImage.cpp檔案，添加到工程中即可。還有一個要說明的是，Gabor濾波器的實作，是用到的網上搜到的Gabor.h Gabor.cpp實作的，隻是進行了簡單的注釋和文檔整理。原作者。。。額，sorry，找的資料太多，不知道出處了，見諒啊！

Gabor.h頭檔案

[cpp]  view plain  copy  print ?

1. #ifndef _GABOR_H  
2. #define _GABOR_H  
3. #include <stdio.h>  
4. #include <iostream>  
5. #include <cv.h>  
6. #define PI 3.14159  
7. #define GAMMA  0.5   //The default value of γ,which is the spatial aspect ratio (sigma_x/sigma_y)  
8. #define RATIO_S2L 0.56  //The default value of σ/λ  
9. #define THETA 45  
10. class Gabor  
11. {  
12. public:  
13.     //@construct 構造函數  
14.     Gabor(float dLambda, float dTheta, float dRatio_S2L = RATIO_S2L, float dGamma = GAMMA, float dPhi = 0);  
15.     //@abolish 析構函數  
16.     ~Gabor();  
17.     //@init 初始化函數  
18.     void init(float dLambda, float dTheta, float dPhi, float dGamma = GAMMA);  
19.     //@init 初始化函數  
20.     void init(float dSigma, float dTheta, float dPhi);  
21.     //@init 初始化函數  
22.     void init();  
23.     //判斷Gabor核心是否建立成功[email protected] find whether the Gabor kernel is created  
24.     bool is_kernel(){ return bKernel; }  
25.     //判斷是否初始化成功[email protected] find whether the parameters is inited  
26.     bool is_init() { return bInit; }  
27.     //判斷初始化的參數是否足夠[email protected] find whether the parameters inited is enough  
28.     bool is_param() { return bParam; }  
29.     //得到核心函數所在矩陣[email protected] the kernel in matrix form  
30.     CvMat* get_Mat() { return pGaborfilter; }  
31.     //得到歸一化圖像[email protected] the kernel in image form  
32.     IplImage* get_NormImage();  
33.     //使用Gabor核函數對輸入圖像進行處理[email protected] the filtering operation to input image with Gabor kernel  
34.     IplImage* do_Filter(const IplImage *src);  
35. protected:  
36.     bool bParam;        //初始化的參數--bool , if the parameters inited are enough  
37.     bool bKernel;        //bool  
38.     bool bInit;            //bool  
39.     float Lambda;      //餘弦函數波長--Wavelength of the cosine factor, which represent the central frequency of Gabor filter   
40.     float Theta;          //核函數的方向--Orientation of the Gabor function, the axis x'  
41.     float Sigma;         // 标準差--The standard deviation of x, and for y , it is Sigma/Gamma;  
42.     float Gamma;       // 空間方向率，指定Gabor函數支援的橢圓率--The spatial aspect ratio  
43.     float Phi;              //Gabor的相位偏移--The phase offset of Gabor  
44.     CvSize GaborWindow;    //Gabor視窗的寬度--The width of  window  
45.     CvMat *pGaborfilter;     //The kernel of Gabor filter  
46. private:  
47.     void create_kernel();  
48. };  
49. #endif  

Gabor.cpp源檔案

[cpp]  view plain  copy  print ?

1. #include "stdafx.h"  
2. #include <iostream>  
3. #include <cv.h>  
4. #include <highgui.h>  
5. #include <cstdlib>  
6. #include "Gabor.h"  
7. Gabor::Gabor(float dLambda, float dTheta, float dRatio_S2L, float dGamma, float dPhi)  
8. {  
9.     Lambda = dLambda;  
10.     Theta = dTheta;  
11.     Sigma = dLambda*dRatio_S2L;  
12.     Gamma = dGamma;  
13.     Phi = dPhi;  
14.     pGaborfilter = NULL;  
15.     bParam = 1;  
16. }  
17. Gabor::~Gabor()  
18. {  
19.     cvReleaseMat(&pGaborfilter);  
20. }  
21. void Gabor::init()  
22. {  
23.     float dtmp;  
24.     int itmp;  
25.     if (is_param() == 0)  
26.     {  
27.         AfxMessageBox("The parameters are not enough!");  
28.         return;  
29.     }  
30.     //沒明白這裡是什麼意思？  
31.     dtmp = sqrt(48 * pow(Sigma, 2) + 1);//根号下( 48*Sigma^2 +1 )  
32.     itmp = cvRound(dtmp);//對一個double型的數進行四舍五入，并傳回一個整型數！  
33.     if (itmp % 2 == 0)  
34.         itmp++;  
35.     GaborWindow.height = GaborWindow.width = itmp;//建立itmp*itmp的Gabor窗函數  
36.     bInit = 1;  
37.     create_kernel();  
38. }  
39. void Gabor::init(float dSigma, float dTheta, float dPhi)  
40. {  
41.     float dtmp;  
42.     int itmp;  
43.     Sigma = dSigma;  
44.     Theta = dTheta;  
45.     Phi = dPhi;  
46.     Gamma = GAMMA;  
47.     Lambda = Sigma / RATIO_S2L;  
48.     bParam = 1;  
49.     dtmp = sqrt(24 * pow(Sigma, 2));  
50.     itmp = cvRound(dtmp);  
51.     if (itmp % 2 == 0)  
52.         itmp++;  
53.     GaborWindow.height = GaborWindow.width = itmp;  
54.     bInit = 1;  
55.     create_kernel();  
56. }  
57. void Gabor::init(float dLambda, float dTheta, float dPhi, float dGamma)  
58. {  
59.     float dtmp;  
60.     int itmp;  
61.     Lambda = dLambda;  
62.     Theta = dTheta;  
63.     Phi = dPhi;  
64.     Gamma = dGamma;  
65.     Sigma = Lambda * RATIO_S2L;  
66.     bParam = 1;  
67.     dtmp = sqrt(24 * pow(Sigma, 2));  
68.     itmp = cvRound(dtmp);  
69.     if (itmp % 2 == 0)  
70.         itmp++;  
71.     GaborWindow.height = GaborWindow.width = itmp;  
72.     bInit = 1;  
73.     create_kernel();  
74. }  
75. void Gabor::create_kernel()  
76. {  
77.     float tmp1, tmp2, xtmp, ytmp, re;  
78.     int i, j, x, y;  
79.     if (is_init() == 0)  
80.     {  
81.         AfxMessageBox("The paremeters haven't been initialed!");  
82.     }  
83.     pGaborfilter = cvCreateMat(GaborWindow.height, GaborWindow.width, CV_32FC1);  
84.     for (i = 0; i < GaborWindow.height; i++)  
85.     {  
86.         for (j = 0; j < GaborWindow.width; j++)  
87.         {  
88.             x = j - GaborWindow.width / 2;  
89.             y = i - GaborWindow.height / 2;  
90.             //源代碼此處的計算公式有誤  
91.             //xtmp = (float)x*cos(Theta) - (float)y*sin(Theta);  
92.             //ytmp = (float)x*sin(Theta) + (float)y*cos(Theta);  
93.             xtmp = (float)x*cos(Theta) + (float)y*sin(Theta);  
94.             ytmp = -(float)x*sin(Theta) + (float)y*cos(Theta);  
95.             tmp1 = exp(-(pow(xtmp, 2) + pow(ytmp*Gamma, 2)) / (2 * pow(Sigma, 2)));  
96.             tmp2 = cos(2 * PI*xtmp / Lambda + Phi);  
97.             // int p=sizeof(float);  
98.             re = tmp1*tmp2;  
99.             cvSetReal2D((CvMat*)pGaborfilter, i, j, re);  
100.         }  
101.     }  
102.     bKernel = 1;  
103. }  
104. IplImage* Gabor::get_NormImage()  
105. {  
106.     if (is_kernel() == 0)  
107.     {  
108.         AfxMessageBox("The filter hasn't bee created!");  
109.         return NULL;  
110.     }  
111.     IplImage *pImg  = cvCreateImage(GaborWindow, IPL_DEPTH_32F, 1);  
112.     IplImage *pImgU8 = cvCreateImage(GaborWindow, IPL_DEPTH_8U, 1);  
113.     CvMat * pMat = cvCreateMat(GaborWindow.height, GaborWindow.width, CV_32FC1);  
114.     cvCopy(pGaborfilter, pImg);  
115.     //歸一化，數組的數值被平移或縮放到一個指定的範圍  
116.     cvNormalize((IplImage*)pImg, (IplImage*)pImg, 0, 255, CV_MINMAX, NULL);  
117.     //使用線性變換轉換輸入數組元素成8位無符号整型   
118.     cvConvertScaleAbs(pImg, pImgU8, 1, 0);  
119.     return pImgU8;  
120.     //return pImg;  
121. }  
122. IplImage * Gabor::do_Filter(const IplImage *src)  
123. {  
124.     //判斷核函數時候建立成功  
125.     if (is_kernel() == false)  
126.     {  
127.         printf("The Gabor Kernel has not been created!");  
128.         return NULL;  
129.     }  
130.     IplImage *pDestImage = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_8U, 1);//建立一個同原圖像相同的目标圖像  
131.     IplImage *tmpImg = cvCloneImage(src);//同cvCreateImage一樣，隻是不需要開辟記憶體空間，直接将原圖像複制到目标圖像  
132.     IplImage *tmpGrayImg = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_8U, 1);  
133.     //判斷是不是有色圖像  
134.     if (tmpImg->nChannels != 1)  
135.     {  
136.         cvCvtColor(tmpImg, tmpGrayImg, CV_BGR2GRAY);//從有色圖轉到灰階圖  
137.     }  
138.     else  
139.     {  
140.         cvReleaseImage(&tmpGrayImg);  
141.         tmpGrayImg = tmpImg;  
142.     }  
143.     CvMat *pGaborKernel = get_Mat();//其實就是CvMat *pGaborKernel=pGaborfilter;  
144.     //Gabor核函數與原圖像進行卷積計算  
145.     cvFilter2D(tmpGrayImg, pDestImage, pGaborKernel, cvPoint((GaborWindow.width - 1) / 2, (GaborWindow.height - 1) / 2));  
146.     cvReleaseImage(&tmpImg);  
147.     cvReleaseImage(&tmpGrayImg);  
148.     return  pDestImage;  
149. }  

對話框Dlg.cpp，這裡隻給出了實作按鈕的實作函數代碼，完整的工程代碼，看文章最後的附錄：

[cpp]  view plain  copy  print ?

1. void CTyreXDlg::OnBnClickedBtndstimg()  
2. {  
3.     // TODO:  在此添加控件通知處理程式代碼  
4.     UpdateData(TRUE);  
5.     CSliderCtrl   *pSlidCtrlTheta = (CSliderCtrl*)GetDlgItem(IDC_SLIDER_THETA);  
6.     int tmptheta = pSlidCtrlTheta->GetPos();//取得目前位置值    
7.     CString sValueTheta = "";  
8.     sValueTheta.Format("%d", tmptheta);  
9.     SetDlgItemText(IDC_StaticTheta, sValueTheta);  
10.     //顯示lambda值  
11.     CSliderCtrl   *pSlidCtrlLambda = (CSliderCtrl*)GetDlgItem(IDC_SLIDER_LAMBDA);  
12.     int tmplambda = pSlidCtrlLambda->GetPos();//取得目前位置值    
13.     CString sValueLambda = "";  
14.     sValueLambda.Format("%d", tmplambda);  
15.     SetDlgItemText(IDC_StaticLambda, sValueLambda);  
16.     //圖像是否加載成功  
17.     if (SrcImg == NULL)  
18.     {  
19.         AfxMessageBox("沒有可處理圖像！！！");  
20.         return;  
21.     }  
22.     //構造函數  
23.     Gabor GaborFilter(tmplambda *0.1, PI*tmptheta / 180, RATIO_S2L, GAMMA, 0);  
24.     //初始化  
25.     GaborFilter.init();  
26.     //擷取歸一化圖像  
27.     IplImage* pGaborNormImg = GaborFilter.get_NormImage();  
28.     //用Gabor核函數對輸入圖像處理，傳回目标圖像  
29.     IplImage* poutGaborimg;  
30.     poutGaborimg = GaborFilter.do_Filter(SrcImg);  
31.     IplImage *pGrayImage;  
32.     // 轉為灰階圖  
33.     pGrayImage = cvCreateImage(cvGetSize(poutGaborimg), IPL_DEPTH_8U, 1);  
34.     //判斷是不是彩色圖像  
35.     if (poutGaborimg->nChannels != 1)  
36.         cvCvtColor(poutGaborimg, pGrayImage, CV_BGR2GRAY);//從彩色圖轉到灰階圖  
37.     else  
38.     {  
39.         cvReleaseImage(&pGrayImage);  
40.         pGrayImage = poutGaborimg;  
41.     }  
42.     //二值化  
43.     IplImage *pBinaryImage;  
44.     pBinaryImage = cvCreateImage(cvGetSize(pGrayImage), IPL_DEPTH_8U, 1);  
45.     cvThreshold(pGrayImage, pBinaryImage, 0, 255, CV_THRESH_BINARY);  
46.     //顯示到picture控件  
47.     ShowImgFunc(pBinaryImage, IDC_ShowDstImg);          //二值化圖像  
48.     //ShowImgFunc(pGaborNormImg, IDC_ShowSrcImg);   //歸一化圖像  
49.     //ShowImgFunc(pGrayImage, IDC_ShowSrcImg);          //二值化之前的灰階圖  
50.     cvReleaseImage(&DstImg);  
51.     cvReleaseImage(&DstGaborImg);  
52.     cvReleaseImage(&DstBinaryImg);  
53. }  

附錄：

1. 上面代碼的完整的工程檔案，可以到這裡下載下傳，不需要積分。我運作成功，但是出現錯誤，或者有不合理的地方，希望各位看到能告訴我一聲，共同進步。

工程檔案完整下載下傳位址：http://download.csdn.net/detail/jorg_zhao/8949247