深入剖析高斯濾波

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

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

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

4、DoG的時(shí)候需要做不同尺度的高斯平滑，在 GaussianConvolution_new1函數(shù)中實(shí)現(xiàn)如下：dim = 1 + 2 * (int) (3.0 * sigma);可見(jiàn)autopano也是實(shí)現(xiàn)的 3*sigma 半徑的窗口。在GIMP里，實(shí)現(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;從效果來(lái)看，這個(gè)實(shí)現(xiàn)的窗口半徑是約等于2.2*sigma 。然后是第二個(gè)問(wèn)題：給定窗口大小，怎么計(jì)算 sigma ?OpenCV 的實(shí)現(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沒(méi)有實(shí)現(xiàn)這個(gè)特性。GIMP的實(shí)現(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;注釋講的很清楚了，基本上就是認(rèn)為si

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

9、板大小的增加，計(jì)算量是呈指數(shù)增長(zhǎng)的。 那么有沒(méi)有一種辦法能夠減少計(jì)算量呢？答案是肯定的。由于高斯函數(shù)可以寫(xiě)成可分離的形式，因此可以采用可分離濾波器實(shí)現(xiàn)來(lái)加速。所謂的可分離濾波器，就是可以把一個(gè)多維的卷積化成多個(gè)一維的卷積。具體到二維的高斯濾波，就是指先對(duì)行做一維卷積，再對(duì)列做一維卷 積。這樣就可以將計(jì)算復(fù)雜度從O(M*M*N*N) 降到O(2*M*M*N), M, N分別是圖像和濾波器的窗口大小。問(wèn)題是實(shí)現(xiàn)的時(shí)候怎么計(jì)算一維的卷積核呢？其實(shí)很簡(jiǎn)單，按照前面計(jì)算出來(lái)的窗口大小，將二維的高斯模板合并成一維，計(jì)算所有離散點(diǎn)上一維高斯函數(shù)的權(quán)值，最后將權(quán)值之和歸一化到1。下面是來(lái)自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