高斯濾波

高斯濾波

第一個問題：高斯函數爲何能做爲圖像處理中的濾波函數？

高斯平滑濾波器不管在空間域仍是在頻率域都是十分有效的低通濾波器，且在實際圖像處理中獲得了工程人員的有效使用．高斯函數具備五個十分重要的性質，它們是：

（1）二維高斯函數具備旋轉對稱性，即濾波器在各個方向上的平滑程度是相同的．通常來講，一幅圖像的邊緣方向是事先不知道的，所以，在濾波前是沒法肯定一個方向上比另外一方向上須要更多的平滑．旋轉對稱性意味着高斯平滑濾波器在後續邊緣檢測中不會偏向任一方向．

（2）高斯函數是單值函數．這代表，高斯濾波器用像素鄰域的加權均值來代替該點的像素值，而每一鄰域像素點權值是隨該點與中心點的距離單調增減的．這一性質是很重要的，由於邊緣是一種圖像局部特徵，若是平滑運算對離算子中心很遠的像素點仍然有很大做用，則平滑運算會使圖像失真．

（3）高斯函數的付立葉變換頻譜是單瓣的．正以下面所示，這一性質是高斯函數付立葉變換等於高斯函數自己這一事實的直接推論．圖像常被不但願的高頻信號所污染(噪聲和細紋理)．而所但願的圖像特徵（如邊緣），既含有低頻份量，又含有高頻份量．高斯函數付立葉變換的單瓣意味着平滑圖像不會被不須要的高頻信號所污染，同時保留了大部分所需信號．

（4）高斯濾波器寬度(決定着平滑程度)是由參數σ表徵的，並且σ和平滑程度的關係是很是簡單的．σ越大，高斯濾波器的頻帶就越寬，平滑程度就越好．經過調節平滑程度參數σ，可在圖像特徵過度模糊(過平滑)與平滑圖像中因爲噪聲和細紋理所引發的過多的不但願突變量(欠平滑)之間取得折衷．

（5）因爲高斯函數的可分離性，大高斯濾波器能夠得以有效地實現．二維高斯函數卷積能夠分兩步來進行，首先將圖像與一維高斯函數進行卷積，而後將卷積結果與方向垂直的相同一維高斯函數卷積．所以，二維高斯濾波的計算量隨濾波模板寬度成線性增加而不是成平方增加．

第二個問題：如何計算高斯函數模板（高斯核）？

其實，只要知道模板的大小和高斯函數的方差sigma,由二維高斯函數的表達式很容易計算出高斯核，只要在歸一化就能夠了。可是由高斯函數的分佈特性可知落在u-3*sigma到u+3*sigma的機率大於百分之九九，因此模板大小的選取每每與sigma的取值是相關的，通常而言：

取dim = 1 + 2 * ((int) (3.0 * sigma));opencv和sift中的源碼也是這麼作的，固然實際中能夠其實沒有這麼嚴格。咱們可使用matlab中的函數直接計算出高斯核，例如3x3的高斯模板：filter=fspecial('gaussian',3,1);其中sigma=1;

sigma的取值決定了高斯函數窗口的大小。在實際中常常看到sigma取值0.8或者1。正常狀況下咱們由高斯函數計算獲得的模板是浮點型數，即double，可是有些狀況咱們爲了加快計算須要將模板處理成整數，對於常見的3x3或者5x5其整數模板以下：

對於浮點數的模板其實每每也都比較固定。不少資料上都有。

第三個問題：可分離的計算

實際上，模板運算（滑動窗口卷積）在數字圖像處理中是一項很是耗時的運算。

以上圖中的3*3高斯模板爲例，每一個像素完成一次模板操做要用9個乘法、8個加法和1個除法。對於一幅n*n的圖像，大約就是9n2個乘法，8n2個加法和n2個除法，這對於比較大的圖像來講，是很是可怕的。並且隨着模板大小的增長，計算量是呈指數增加的。那麼有沒有一種辦法可以減小計算量呢？答案是確定的。因爲高斯函數能夠寫成可分離的形式，所以能夠採用可分離濾波器實現來加速。所謂的可分離濾波器，就是能夠把一個多維的卷積化成多個一維的卷積。具體到二維的高斯濾波，就是指先對行作一維卷積，再對列作一維卷積。這樣就能夠將計算複雜度從O(M*M*N*N)降到O(2*M*M*N)，M, N分別是圖像和濾波器的窗口大小。問題是實現的時候怎麼計算一維的卷積核呢？

其實很簡單，按照前面計算出來的窗口大小，將二維的高斯模板合併成一維，計算全部離散點上一維高斯函數的權值，最後將權值之和歸一化到1。

第四個問題：C代碼實現？

 

 因爲筆者使用的是整數模板，先將浮點數模板代碼列出請參考http://blog.csdn.net/crzy_sparrow/article/details/6998745

1. //  一維高斯分佈函數，用於平滑函數中生成的高斯濾波係數  
2. /*  
3.  *  @parameter sigma:     高斯函數參數  
4.  *  @parameter pdKernel:    高斯核函數模板  
5.  *  @parameter pnWidowSize:  高斯模板大小  
6.  */  
7. void CreatGauss(double sigma, double **pdKernel, int *pnWidowSize)  
8. {  
9.     LONG i;  
10.       
11.     int nCenter;//數組中心點  
12.     double dDis;//數組中一點到中心點距離  
13.   
14.     //中間變量  
15.     double dValue;  
16.     double dSum;  
17.     dSum = 0;  
18.       
19.     *pnWidowSize = 1+ 2*ceil(3*sigma);// [-3*sigma,3*sigma] 之內數據，會覆蓋絕大部分濾波係數  
20.   
21.     nCenter = (*pnWidowSize)/2;  
22.     //生成高斯數據  
23.     for(i=0;i<(*pnWidowSize);i++)  
24.     {  
25.         dDis = (double)(i - nCenter);  
26.         dValue = exp(-(1/2)*dDis*dDis/(sigma*sigma))/(sqrt(2*3.1415926)*sigma);  
27.         (*pdKernel)[i] = dValue;  
28.         dSum+=dValue;  
29.     }  
30.     //歸一化  
31.     for(i=0;i<(*pnWidowSize);i++)  
32.     {  
33.         (*pdKernel)[i]/=dSum;  
34.     }  
35. }  

 

1. //用高斯濾波器平滑原圖像  
2. /* 
3.  *  @parameter sz   :  圖像尺寸 
4.  *  @parameter pGray   :   圖像灰度值 
5.  *  @parameter pResult:  圖像 
6.  *  @parameter sigma:     高斯函數參數 
7.  */  
8. void GaussianSmooth(SIZE sz, LPBYTE pGray, LPBYTE pResult, double sigma)  
9. {  
10.     LONG x, y;  
11.     LONG i;  
12.       
13.     int nWindowSize;//高斯濾波器長度  
14.       
15.     int nLen;//窗口長度  
16.       
17.     double *pdKernel;//一維高斯濾波器  
18.       
19.     double dDotMul;//高斯係數與圖像數據的點乘  
20.   
21.     double dWeightSum;//濾波係數總和  
22.   
23.     double *pdTemp;  
24.   
25.     nWindowSize = 1+ 2*ceil(3*sigma);// [-3*sigma,3*sigma] 之內數據，會覆蓋絕大部分濾波係數  
26.   
27.   
28.     if ((pdTemp = (double *)malloc(sz.cx*sz.cy*sizeof(double)))==NULL)  
29.     {  
30.         printf("melloc memory for pdTemp failed!!");  
31.         exit(0);  
32.     }  
33.     if ((pdKernel = (double *)malloc(nWindowSize*sizeof(double)))==NULL)  
34.     {  
35.         printf("malloc memory for pdKernel,failed!!");  
36.         exit(0);  
37.     }  
38.   
39.     //產生一維高斯數據  
40.     CreatGauss(sigma, &pdKernel, &nWindowSize);  
41.   
42.     nLen = nWindowSize/2;  
43.   
44.     //x方向濾波  
45.     for(y=0;y<sz.cy;y++)  
46.     {  
47.         for(x=0;x<sz.cx;x++)  
48.         {  
49.             dDotMul = 0;  
50.             dWeightSum = 0;  
51.             for(i=(-nLen);i<=nLen;i++)  
52.             {  
53.                 //判斷是否在圖像內部  
54.                 if((i+x)>=0 && (i+x)<sz.cx)  
55.                 {  
56.                     dDotMul+=(double)(pGray[y*sz.cx+(i+x)] * pdKernel[nLen+i]);  
57.                     dWeightSum += pdKernel[nLen+i];  
58.                 }  
59.             }  
60.             pdTemp[y*sz.cx+x] = dDotMul/dWeightSum;     
61.         }  
62.     }  
63.   
64.     //y方向濾波  
65.     for(x=0; x<sz.cx;x++)  
66.     {  
67.         for(y=0; y<sz.cy; y++)  
68.         {  
69.             dDotMul = 0;  
70.             dWeightSum = 0;  
71.             for(i=(-nLen);i<=nLen;i++)  
72.             {  
73.                 if((i+y)>=0 && (i+y)< sz.cy)  
74.                 {  
75.                     dDotMul += (double)pdTemp[(y+i)*sz.cx+x]*pdKernel[nLen+i];  
76.                     dWeightSum += pdKernel[nLen+i];  
77.                 }  
78.             }  
79.             pResult[y*sz.cx+x] = (unsigned char)dDotMul/dWeightSum;  
80.         }  
81.     }  
82.       
83.     free(pdTemp);//釋放內存  
84.     free(pdKernel);