深入剖析高斯濾波

1、高斯濾波（高斯平滑）是數(shù)字圖像處理和計算機視覺里面最常見的操作。平時，我們都是用Matlab 或者OpenCV 的函數(shù)調(diào)用：imfilter 或者cvSmooth ,并不關(guān)心底層的實現(xiàn)。然 而當開發(fā)者要自己實現(xiàn)高斯濾波的時候，往往就會很迷惘，經(jīng)常會被下面幾個問題所困擾：1、給定sigma，即標準偏差，怎么確定離散化后的濾波器窗口大?。?、 給定窗口大小，怎么計算高斯核的sigma ,即標準偏差？3、怎么實現(xiàn)可分離濾波器？在網(wǎng)上搜了一下，還真沒幾個人把實現(xiàn)的細節(jié)講清楚了。這里，我嘗試結(jié)合三份源碼，做個小小的總結(jié)。三份源碼分別是：1、OpenCV 中的 cvfilter.cpp2、autopano

2、-sift-c 中的 GaussianConvolution.c3、GIMP 中的 blur-gauss.c 和 unsharp-mask.c在圖像處理中，高斯濾波一般有兩種實現(xiàn)方式，一種是用離散化的滑動窗口進行卷積運算，另一種則是通過傅里葉變換來實現(xiàn)。最常見的就是第一種卷積核實現(xiàn)，只有當離散化的窗口非常大，計算非常耗時（這時可使用可分離濾波器）的情況下，可能會考慮基于傅里葉變換 的實現(xiàn)方法。這里我們只討論第一種方法。二維高斯函數(shù)的形式是這樣的：I (J-Jo)2 , (v-yp2 y) = Ae ' 氣 2try有著如下的形狀:基本上，離散化的主旨就是保留高斯函數(shù)中心能量最集中的中間

3、部分，忽略四周能量很小的平坦區(qū)域。這只是個很感性的描述， 具體實現(xiàn)起來，就會出現(xiàn)千奇百怪的版本。下面結(jié)合三 份源碼，看看現(xiàn)實世界里的高斯平滑到底長的什么樣子。首先是第一個問題：給定 sigma ,怎么計算窗口大??？直接上OpenCV 的源碼，在cvFilter函數(shù)中：paraml = cvRound(sigma1*(depth = CV_8U ? 3 : 4)*2 + 1)|1;OpenCV認為半徑為3*sigma的窗口就是高斯函數(shù)的能量最集中的區(qū)域。(在圖像位深度不是8的時候，使用4*sigma 半徑的窗口？)autopan0-sift-c是圖像拼接軟件 hugin里面的sift實現(xiàn)，在實現(xiàn)

4、DoG的時候需要做不同尺度的高斯平滑，在 GaussianConvolution_new1函數(shù)中實現(xiàn)如下：dim = 1 + 2 * (int) (3.0 * sigma);可見autopano也是實現(xiàn)的 3*sigma 半徑的窗口。在GIMP里，實現(xiàn)比較奇特，在 blur_gauss.c 的make_rle_curve函數(shù)里面，const gdouble sigma2 = 2 * sigma * sigma;const gdouble l = sqrt (-sigma2 * log (1.0 / 255.0);int n = ceil (l) * 2;if (n % 2) = 0)n +=

5、1;從效果來看，這個實現(xiàn)的窗口半徑是約等于2.2*sigma 。然后是第二個問題：給定窗口大小，怎么計算 sigma ?OpenCV 的實現(xiàn)，在 cvFilter.cpp 的 init_gaussian_kernel函數(shù)中:sigmaX = sigma > 0 ? sigma : (n/2 1)*0.3 + 0.8;sigma的大小約為窗口半徑的0.3倍再加上0.8 。autopano沒有實現(xiàn)這個特性。GIMP的實現(xiàn)：/* we want to generate a matrix that goes out a certain radius* from the center, so we

6、 have to go out ceil(rad-0.5) pixels,* inlcuding the center pixel. Of course, thatn 6 nes otiniction,* so we have to go the same amount in the other direction, but not count* the center pixel again. So we double the previous result and subtract* one.* The radius parameter that is passed to this func

7、tion is used ass a little confusing.* the standard deviation, and the radius of effect is the* standard deviation * 2. It*/ radius = fabs (radius) + 1.0;std_dev = radius;radius = std_dev * 2;/* go out ' radius ' in each direction */matrix_length = 2 * ceil (radius 0.5) + 1;注釋講的很清楚了，基本上就是認為si

8、gma應(yīng)該等于窗口半徑的一半?？赐赀@三份源碼，結(jié)論就是，關(guān)于 sigma和窗口的大小，你愛怎么算都可以，這個由你的 實際需要決定，別太離譜就行。(根據(jù)概率理論，高斯分布的樣本絕大部分集中在3倍標準差以內(nèi)的區(qū)域。)第三個問題是可分離濾波器：首先說明為什么要使用可分離濾波器。實際上，模板運算(滑動窗口卷積)在數(shù)字圖像處理中是一項非常耗時的運算。1 X 16。2 1'2 42121I以上圖中的3*3高斯模板為例，每個像素完成一次模板操作要用9個乘法、8個加法和1個除法。對于一幅n*n的圖像，大約就是 9n2個乘法，8n 2個加法和n2個除法，這對于比 較大的圖像來說，是非?？膳碌?。而且隨著模

9、板大小的增加，計算量是呈指數(shù)增長的。 那么有沒有一種辦法能夠減少計算量呢？答案是肯定的。由于高斯函數(shù)可以寫成可分離的形式，因此可以采用可分離濾波器實現(xiàn)來加速。所謂的可分離濾波器，就是可以把一個多維的卷積化成多個一維的卷積。具體到二維的高斯濾波，就是指先對行做一維卷積，再對列做一維卷 積。這樣就可以將計算復(fù)雜度從O(M*M*N*N) 降到O(2*M*M*N), M, N分別是圖像和濾波器的窗口大小。問題是實現(xiàn)的時候怎么計算一維的卷積核呢？其實很簡單，按照前面計算出來的窗口大小，將二維的高斯模板合并成一維，計算所有離散點上一維高斯函數(shù)的權(quán)值，最后將權(quán)值之和歸一化到1。下面是來自O(shè)penCV 的源碼：for( i = 0; i <= n/2; i+ )(double t = fixed_kernel ? (double)fixed_kerneli : exp(scale2X*i*i);if( type = CV_32FC1 )(cf(n/2+i)*step = (float)t;sum += cf(n/2+i)*step*2;)else(cd(n/2+i)*step = t;sum += cd(n/2+i)*step*2;)sum = 1./sum;for( i = 0; i &l