高斯濾波的理解與學習

微信公眾號：幼兒園的學霸

目錄

文章目錄

• 高斯濾波的理解與學習
• • 目錄
  • 前言
  • 高斯函數
  • • 一維高斯函數
    • 二維高斯函數
  • 高斯濾波過程
  • • 高斯核求解
    • 利用高斯核濾波
    • 高斯濾波步驟
  • 高斯濾波實現
  • • 高斯濾波標準差與窗口大小的換算
    • 實現
    • • 常規實現
      • 分離實現高斯濾波
  • 總結
  • • 高斯函數性質
    • 高斯濾波應用
  • 參考資料
  • 附錄
  • • 高斯函數及頻譜繪圖代碼

前言

對一幅圖像而言,低頻部分對應整體灰度級的顯示,高頻部分對應著圖像的細節部分.因此去掉低頻部分(或者增強高頻部分)可以銳化圖像,去掉高頻部分(或者增強低頻部分)可以實現模糊/平滑圖像的作用.

去除低頻或者高頻部分可以通過濾波器來完成,但是在時域空間并不能直觀看出所謂的高頻低頻成分,所以就需要做個傅里葉變換,轉換成頻域,然后直接在頻域上把需要去除的頻率對應的“基”的權值置零即可.

高斯濾波器是一種線性低通濾波器,能夠有效的抑制噪聲,平滑圖像.其本質是帶權值的加權均值濾波,權值大小與濾波器中元素到濾波中心距離有關,卷積核中心權重數值最大,并向四周減小,減小的幅度并不是隨意的,而是要求整個卷積核近似高斯函數的圖像.高斯濾波器對于抑制服從正態分布的噪聲非常有效.

線性濾波器(linear filter) :輸出圖像上每個像素點的值都是由輸入圖像各像素點值加權求和的結果.其原始數據與濾波結果是一種算術運算,即用加減乘除等運算實現,如方框濾波(BoxFilter)、均值濾波(MeanBlur)、高斯濾波(GaussianBlur)等.由于線性濾波器是算術運算,有固定的模板,因此濾波器的轉移函數是可以確定并且是唯一的.

非線性濾波器(non-linear filter):非線性濾波的算子中包含了取絕對值、置零等非線性運算.其原始數據與濾波結果是一種邏輯關系,即用邏輯運算實現,如最大值濾波器、最小值濾波器、中值濾波（medianBlur）和雙邊濾波（bilateralFilter）等,是通過比較一定鄰域內的灰度值大小來實現的,沒有固定的模板,因而也就沒有特定的轉移函數（因為沒有模板作傅里葉變換）,另外,膨脹和腐蝕也是通過最大值、最小值濾波器實現的.

思考1:為什么要進行 加權均值 ?一幅圖像中輪廓和邊緣這種變化比較劇烈的地方高頻信號(噪聲和噪聲周邊的像素相比,其灰度變化同樣比較劇烈),其他的部分(可反應整副圖像的強度)是低頻信號,所以經過低通濾波器后圖像的輪廓和邊緣信號會被濾掉一部分,直觀感受是圖片變模糊了,這是平滑.從另一方面講,加權是將一個像素用鄰域幾個像素的均值代替,這樣會使信號的變化變得平緩(拉大小的值,拉小大的值),信號變化劇烈的地方是高頻信號,所以,這樣就達到了低通的效果,而低通可以濾除高頻噪聲.因此加權的目的在于實現低通.
思考2:如何實現 加權均值 中的求平均目的? 濾波核前面的系數等于矩陣中所有元素之和的倒數

高斯函數

計算高斯濾波的濾波核之前先了解一下高斯函數.

高斯函數在圖像處理中出現的頻率相當高.

一維高斯函數

$$f(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{\frac{?(x?\mu {)}^2}{2{\sigma }^2}}$$
其中,μ是x的均值,σ是x的方差.本質上,x,μ都是空間中的坐標,x是卷積核內任一點的坐標,μ是卷積核中心的坐標.
由于每次計算時,都是以當前計算點為原點(中心點就是原點),因此均值μ等于0.所以公式進一步簡化為:
$$f(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{\frac{?x^2}{2{\sigma }^2}}$$
其函數曲線如下圖所示:

可以看到其是一鐘形曲線,σ決定了分布的形狀或者說表征了鐘的寬度,σ越小形狀越瘦高,σ越大越矮胖

二維高斯函數

在計算機視覺中,圖像是二維的,因此引入下面的二維高斯函數.
二維高斯函數為X,Y兩個方向的一維高斯函數的乘積:
$$f(x,y)=f(x)f(y)=\frac{1}{\sqrt{2\pi }{\sigma }_x}{e}^{?\frac{{\left(x?{\mu }_x\right)}^2}{2{\sigma }_x^2}}\frac{1}{\sqrt{2\pi }{\sigma }_y}{e}^{?\frac{{\left(x?{\mu }_y\right)}^2}{2{\sigma }_y^2}}$$
同樣,在圖像濾波中,一般情況下${\mu }_x={\mu }_y=0$,因此二維高斯函數可表示如下:
$$f(x,y)=f(x)f(y)=\frac{1}{2\pi {\sigma }^2}{e}^{?\frac{x^2+y^2}{2{\sigma }^2}}$$

其中,(x,y)為點的坐標,在圖像處理中,其為整數.
其函數曲線如下圖所示:

從函數曲線可以看到,在整個定義域內,高斯函數值都是非負的,$f(x,y)>0$成立,這也決定了高斯核中所有元素的值均為正值.

高斯濾波過程

高斯核求解

要想得到一個高斯核,可以對高斯函數進行離散化,得到的高斯函數值作為高斯核的元素.例如:要產生一個3x3的高斯核,以高斯核的中心位置為坐標原點進行離散取樣, 那么高斯核中心的8鄰域坐標為:

將各個位置的坐標代入到高斯函數中,得到的值就是初步的高斯核:

假設σ=1,則上面的高斯核求解結果如下:

上面初步計算3x3高斯核的權重和為0.7794836.考慮這樣一個問題,假若某一鄰域內所有像素的灰度值為255,那么通過該高斯核進行卷積之后,模板中心像素的灰度值為0.7794836×255=198.77 < 255,偏離了實際的灰度值,使得圖像亮度相比原圖偏暗,產生了誤差.因此有必要對該高斯核進行歸一化處理,保證權重和為1以保證圖像的均勻灰度區域不受影響.

歸一化的原因更加正確的表述應該為:對灰度級為常數的圖像區域,高斯模板的響應該為1.這樣才能確保灰度級為常數的圖像區域經過高斯模板處理之后,依舊為其本身,而不是被改變為其他灰度級.

將求得的高斯核中各權重值除以權重和進行歸一化后如下:

作為對比,在matlab中,z = fspecial('gaussian', [3 3], 1);得到的結果和上面結果一致.

上面歸一化后求解的高斯核為小數形式的.

除了小數形式,高斯核也可以寫為整數形式.只用對求得的初始高斯核進行取整處理,就能得到整數形式的高斯核.從高斯函數圖像可以看到:1)從曲線中心位置到兩端,高斯函數值是單調遞減的;2)曲線位于x軸上方,其值總是正的.因此高斯核中距離中心最遠的位置(如高斯核左上角)的元素總是最小的,如果進行取整,那么左上角元素的最小值也要保證≥1.所以對求解整數形式的高斯核過程如下:
1). 高斯核左上角的值歸一化為1.即高斯核中所有元素均除以左上角的元素.在本例中,左上角的元素為0.058549833,歸一化后高斯核為:

2).取整(對于取整方式,未查找到具體的原則.有文章選擇了選擇向下取整的方式,也許有向上取整,或者四舍五入等方式,但是從不同取整方式的結果來看,四舍五入的取整方式更接近于小數形式的結果).
3).取整后,高斯核所有元素的和明顯大于1,為了保證圖像的均勻灰度區域不受影響,因此需要進行權重歸一化,給高斯核一系數,以保證高斯和所有元素的和與系數的乘積為1.顯然,該系數為高斯核和的倒數,該例中取整、歸一化后高斯核為:

總結高斯核求解過程圖示如下:

利用高斯核濾波

濾波計算過程和CNN中的單通道卷積過程一致,只是此處濾波時的步長為1,濾波核為單通道,而CNN中卷積步長可以指定為1,2等值.

有了高斯核,就可以對圖像進行高斯濾波了.具體濾波過程和其他濾波方法一致.以一個3x3數據為例,利用上面求解的高斯核進行濾波.
其中心點的高斯濾波過程如下:

圖中中心像素值為50,經過濾波后其值仍為50,這只是由于數字選擇的原因導致其湊巧在數值上相等而已.如果將中心像素值修改為52,其濾波后值也為50.

在邊界附近的點,其鄰域內沒有足夠的點進行濾波處理,因此在濾波之前需要進行上下、左右邊界擴充操作,各邊界擴充的寬度為濾波半徑大小.具體的邊界擴充類型可以參考OpenCV中BorderTypes指定的類型進行.
考慮此,對上述3x3輸入進行濾波的過程如下:

易計算得濾波前輸入像素的和為450,濾波后輸出矩陣各元素相加后的和仍為450,即圖像的亮度在經過該高斯核濾波前后沒有發生變化.

高斯濾波步驟

按照上面的流程,可以歸納高斯濾波的步驟如下:
1.生成高斯核.根據選擇的濾波半徑以及標準差,確定高斯核鄰域內各點到中心點的距離,代入二維高斯函數,求解高斯核各位置值,然后進行歸一化等處理;
2.邊界擴充.根據高斯核大小對輸入圖像進行擴充邊界操作;
3.卷積.針對圖像矩陣中的某個元素P,把高斯核的中心和待處理元素P空間對齊,將高斯核中各權值和P的鄰域內元素相乘后相加,作為P的輸出值;
4.對圖像中的每個元素(不包含邊界擴充的元素)進行步驟3的操作,得到的結果就是高斯濾波的結果.

從上面的步驟可以看到,高斯濾波僅考慮了圖像中各元素空間位置之間的關系,沒有將像素值之間的關系考慮進去. 而之前講到的雙邊濾波就不僅考慮像素空間之間的關系,同時考慮像素值之間的關系,它是空間域和像素值域之間的綜合結果.并且空間域和像素值域的求解都是套用高斯函數公式得到的.只是f(x,y)中x,y,μ的含義不同:針對空間域,它和這里的高斯核是一樣的,x,y,μ表示的是空間位置,而μ=0;針對值域,x,y表示每一點的像素值,μ表示中心的像素值.所以它不光考慮了像素在空間中位置遠近程度的影響,還考慮了像素亮度相近程度的影響.從這個角度來看,雙邊濾波就很好理解了.

高斯濾波實現

根據上面的濾波步驟,可以比較方便的實現高斯濾波代碼.

高斯濾波標準差與窗口大小的換算

如下圖所示為高斯函數的分布特點.一般3σ外的數值已接近于0,可忽略.數值分布在(μ—3σ,μ+3σ)中的概率為0.997,或者說以均值為中心,半徑為3σ的范圍內已經包含了高斯函數99.7%以上的信息。如果高斯核的尺寸大,那么相應的標準差σ也更大;相對應,如果σ大,那么意味著高斯核的尺寸也要相應的進行增加.因此當我們知道高斯核的大小后,可以推斷出標準差的值.

我們大概可以感受到高斯核的半徑大約為3σ的大小.據此,在OpenCV中,有下面公式:
$$\sigma =0.3?(\frac{ksize}{2}?1)+0.8$$
推導標準差的取值.公式中0.3比較好理解,但0.8的來歷并不是很清楚.

當知道標準差而不知道高斯核尺寸時,同樣可以由上式進行反向求解,得到高斯核尺寸,公式如下:
$$ksize=[2?(3\sigma )+1]|1$$
公式為先取整,然后和1進行按位或運算,以保證高斯核尺寸為奇數.
OpenCV中該部分代碼如下:

代碼中半徑為σ的3或4倍.

4倍的來源不明,無資料說明

實現

常規實現

常規實現過程按照上面的公式計算即可.其中高斯核或標準差的計算按上節公式實現.
代碼如下:

#include <iostream> #include <chrono> #include <opencv2/opencv.hpp> #include <numeric>/*!* 生成小數形式的二維高斯核* @param ksize 高斯核大小,奇數.當ksize<=0,并且sigma>0時,自動根據sigma計算ksize* @param sigma 高斯濾波標準差.當sigma<=0,并且ksize>0為奇數時,自動根據ksize計算sigma* @return 小數形式的高斯核.CV_32FC1.尺寸:ksize * ksize* @author liheng*/ cv::Mat GetGaussianKernel( int ksize, float sigma) {//根據sigma自動計算ksizeif( ksize <= 0 && sigma > 0 )ksize = cvRound((sigma*3)*2 + 1)|1;// | 1操作以保證ksize為奇數//校驗ksizeCV_Assert( ksize > 0 && ksize % 2 == 1);//未設定sigma值時,根據窗口大小進行計算sigma = (sigma>0 ? sigma : 0.3f*((ksize-1)*0.5f - 1) + 0.8f);auto sigma2_inverse = 1.0/(2*sigma*sigma);// 1/(2*sigma^2)cv::Mat kernel(ksize, ksize,CV_32FC1);const int center = ksize/2;//高斯核中心坐標for(int row=0;row<ksize;++row){int y=row-center; //以center為中心,第row行的 y坐標auto y2 = y*y;//y^2for(int col=0;col<ksize;++col){int x=col-center; //以center為中心,第col列的x坐標//float f_xy = 1.0/(2*CV_PI*sigma*sigma) * std::exp( -(x*x+y*y)/(2*sigma*sigma) );float f_xy = std::exp(-(x*x+y2)*sigma2_inverse );//減少不必要計算;系數在后面歸一化時會被消去;kernel.at<float>(row,col) = f_xy;}}//小數形式高斯核auto sum = cv::sum(kernel)[0];kernel /= sum;整數形式高斯核左上角系數歸一化為1//kernel /= kernel.at<float>(0,0);向下取整//kernel.convertTo(kernel,CV_32SC1);//向上取整?//kernel.convertTo(kernel,CV_32FC1);//系數未計算return kernel; }/*!* 利用卷積核對輸入圖像進行二維卷積(濾波)* @param src 輸入.CV_8UC1* @param dst 輸出.CV_8UC1* @param kernel 卷積核(相關核).CV_32FC1. kernel.size()需要為奇數* @note 該函數在OpenCV中的對應函數:cv::filter2D(...)* @author liheng*/ void Filter2D(const cv::Mat& src,cv::Mat& dst,const cv::Mat& kernel) {CV_Assert( kernel.cols % 2 == 1 && kernel.rows % 2 == 1);//邊界填充int border_lr = kernel.cols/2;//左右填充寬度int border_tb = kernel.rows/2;//上下填充寬度cv::Mat borderedSrc;cv::copyMakeBorder(src,borderedSrc,border_tb,border_tb,border_lr,border_lr,cv::BORDER_REFLECT);dst.create(src.size(),src.type());//對圖像中各元素(不包含填充元素)分別進行濾波int rows = borderedSrc.rows;int cols = borderedSrc.cols;for(int row=border_tb; row<rows-border_tb; ++row){for(int col=border_lr; col<cols-border_lr; ++col){//遍歷卷積核鄰域求加權和作為該點輸出值float sum = 0;for(int y=-border_tb; y<=border_tb; ++y)//鄰域高度:2*border_tb+1{for(int x=-border_lr;x<=border_lr; ++x)//鄰域寬度:2*border_lr+1sum += kernel.at<float>(y+border_tb,x+border_lr)*borderedSrc.at<uchar>(row+y,col+x);//權重與元素相乘}//輸出dst.at<uchar>(row-border_tb,col-border_lr) = cv::saturate_cast<uchar>(sum);}} }/*!* 自定義實現高斯濾波.原始版本,未采用分離卷積進行加速* @param src 輸入.CV_8UC1* @param dst 輸出.CV_8UC1* @param ksize 高斯核大小,奇數.當ksize<=0,并且sigma>0時,自動根據sigma計算ksize* @param sigma 高斯濾波標準差.當sigma<=0,并且ksize>0為奇數時,自動根據ksize計算sigma* @author liheng*/ void MyGaussianFilter(const cv::Mat& src,cv::Mat& dst,int ksize,float sigma) {//生成二維高斯核cv::Mat kernel = GetGaussianKernel(ksize,sigma);//濾波Filter2D(src,dst,kernel); }int main(int argc, char *argv[]) {cv::Mat kernel = GetGaussianKernel(3,1);std::cout<<"二維高斯核:"<<std::endl;std::cout<<kernel<<std::endl;cv::Mat data = (cv::Mat_<uchar>(3,3)<<10,20,30,40,50,60,70,80,90);//cv::Mat dst;MyGaussianFilter(data,dst,3,1);std::cout<<"自定義實現結果::"<<std::endl;std::cout<<dst<<std::endl;cv::Mat dst2;cv::filter2D(data,dst2,CV_8UC1,kernel,cv::Point(-1,-1),0,cv::BORDER_REFLECT);std::cout<<"filter2D結果::"<<std::endl;std::cout<<dst2<<std::endl;cv::GaussianBlur(data,data,cv::Size(3,3),1,1,cv::BORDER_REFLECT);std::cout<<"cv::GaussianBlur結果::"<<std::endl;std::cout<<data<<std::endl;//輸出結果如下:/**二維高斯核:[0.075113609, 0.1238414, 0.075113609;0.1238414, 0.20417996, 0.1238414;0.075113609, 0.1238414, 0.075113609]自定義實現結果::[ 21, 28, 35;43, 50, 57;65, 72, 79]filter2D結果::[ 21, 28, 35;43, 50, 57;65, 72, 79]cv::GaussianBlur結果::[ 21, 28, 35;43, 50, 57;65, 72, 79]**/return 0; }

可以參考上面圖片中流程以輔助理解.

上面的濾波過程,其循環次數為row*col*(ksize*ksize),其中row,col為圖像的尺寸.其時間復雜度為$O(N^2)$,隨著高斯核尺寸的增加而呈平方增長,當高斯核尺寸較大時,其運算效率比較低.為了提高運算速度,可以將二維的高斯核進行拆分,然后分別進行濾波.

分離實現高斯濾波

從二維高斯函數的公式可以看到,其為兩個一維高斯函數的乘積.因此二維高斯核具有可分離性,可以將高斯核分為兩個一維高斯核的乘積,如下圖所示:

同樣,利用高斯函數進行卷積(高斯濾波)的過程也具有可分離性. 其示意圖如下:

證明過程如下.
高斯濾波函數如上面所示的$f(x,y)$,輸入圖像為$I(x,y)$,兩者進行卷積運算,記邊界索引為$a=\frac{n?1}{2}$,則兩者卷積過程表示如下:

兩者之間為 離散二維卷積.下面公式即為離散二維卷積的定義.

$$\begin{array}{rl}f(x,y)?I(x,y)& =\sum _{i=?a}^a\sum _{j=?a}^{a}f(i,j)I(x?i,y?j)\\ & =\sum _{i=?a}^a\sum _{j=?a}^{a}I(x?i,y?j)\frac{1}{2\pi {\sigma }^2}{e}^{?\frac{i^2+j^2}{2{\sigma }^2}}\\ & =\sum _{i=?a}^a\left[\frac{1}{\sqrt{2\pi }\sigma }{e}^{\frac{?i^2}{2{\sigma }^2}}?\left[\sum _{j=?a}^{a}I(x?i,y?j)?\frac{1}{\sqrt{2\pi }\sigma }{e}^{\frac{?j^2}{2{\sigma }^2}}\right]\right]\end{array}\text{\hspace{1em}(1)}$$
記$f(y)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{\frac{?y^2}{2{\sigma }^2}}$,則上面第2個中括號部分整理如下:
$$\begin{array}{r}\sum _{j=?a}^{a}I(x?i,y?j)?\frac{1}{\sqrt{2\pi }\sigma }{e}^{\frac{?j^2}{2{\sigma }^2}}=f(y)?I(x?i,y)\end{array}\text{\hspace{1em}(2)}$$
可視為圖像I與f(y)之間的卷積運算.

記$f(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{\frac{?x^2}{2{\sigma }^2}}$,同時將公式(2)代入公式(1)整理:
$$\begin{array}{rl}f(x,y)?I(x,y)& =\sum _{i=?a}^a\left[\frac{1}{\sqrt{2\pi }\sigma }{e}^{\frac{?i^2}{2{\sigma }^2}}?\left[f(y)?I(x?i,y)\right]\right]\\ & =f(x)?\left(f(y)?I(x,y)\right)\end{array}\text{\hspace{1em}(3)}$$
在上面推導過程中,f(y),f(x)并不代表具體的行列順序,f(x)可以表示行方向，也可以代表列方向,f(y)同理.因此公式(3)可整理為:
$$\begin{array}{rl}f(x,y)?I(x,y)& =f(x)?\left(f(y)?I(x,y)\right)\\ & =I(x,y)?f(y)?f(x)\\ & =I(x,y)?f(x)?f(y)\end{array}\text{\hspace{1em}(4)}$$

從而:高斯函數進行卷積(高斯濾波)的過程同樣具有可分離性得證.

上面結論可以進行推廣.如果卷積核kernel是可分離的,并且$kernel=kernel1?kernel2$,則對圖像I進行卷積有:
$$\begin{array}{rl}I?kernel& =I?(kernel1?kernel2)\\ & =(I?kernel1)?kernel2\end{array}$$
需要注意的是,上面的卷積核kernel為一維水平方向和一維垂直方向共2個卷積核的全卷積,這2個卷積核是滿足交換律的,但是對于多個卷積核的全卷積,其并不滿足交換律.
對于均值濾波,其卷積核也是可分離的,如下:
$$\frac{1}{9}?\begin{array}{lll}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}?=\frac{1}{3}?\begin{array}{l}1\\ 1\\ 1\end{array}??\frac{1}{3}\left[\begin{array}{lll}1& 1& 1\end{array}\right]$$
用于邊緣檢測的Sobel卷積核也是可分離的,如下:
$$?\begin{array}{lll}?1& 0& 1\\ ?2& 0& 2\\ ?1& 0& 1\end{array}?=?\begin{array}{l}1\\ 2\\ 1\end{array}?\times \left[\begin{array}{lll}?1& 0& 1\end{array}\right]$$

CNN中非對稱卷積(Asymmetric Convolution)/空間可分離卷積(Spatial Separable Convolution) 的思想和這里有點相似.

按照公式的推導結論,可對原始高斯濾波進行優化,改寫為可分離卷積的形式.
為方便,上一小節代碼同樣拷貝到該代碼中.

#include <iostream> #include <chrono> #include <opencv2/opencv.hpp> #include <numeric>/*!* 生成小數形式的二維高斯核* @param ksize 高斯核大小,奇數.當ksize<=0,并且sigma>0時,自動根據sigma計算ksize* @param sigma 高斯濾波標準差.當sigma<=0,并且ksize>0為奇數時,自動根據ksize計算sigma* @return 小數形式的高斯核.CV_32FC1.尺寸:ksize * ksize* @author liheng*/ cv::Mat GetGaussianKernel( int ksize, float sigma) {//根據sigma自動計算ksizeif( ksize <= 0 && sigma > 0 )ksize = cvRound((sigma*3)*2 + 1)|1;// | 1操作以保證ksize為奇數//校驗ksizeCV_Assert( ksize > 0 && ksize % 2 == 1);//未設定sigma值時,根據窗口大小進行計算sigma = (sigma>0 ? sigma : 0.3f*((ksize-1)*0.5f - 1) + 0.8f);auto sigma2_inverse = 1.0/(2*sigma*sigma);// 1/(2*sigma^2)cv::Mat kernel(ksize, ksize,CV_32FC1);const int center = ksize/2;//高斯核中心坐標for(int row=0;row<ksize;++row){int y=row-center; //以center為中心,第row行的 y坐標auto y2 = y*y;//y^2for(int col=0;col<ksize;++col){int x=col-center; //以center為中心,第col列的x坐標//float f_xy = 1.0/(2*CV_PI*sigma*sigma) * std::exp( -(x*x+y*y)/(2*sigma*sigma) );float f_xy = std::exp(-(x*x+y2)*sigma2_inverse );//減少不必要計算;系數在后面歸一化時會被消去;kernel.at<float>(row,col) = f_xy;}}//小數形式高斯核auto sum = cv::sum(kernel)[0];kernel /= sum;整數形式高斯核左上角系數歸一化為1//kernel /= kernel.at<float>(0,0);向下取整//kernel.convertTo(kernel,CV_32SC1);//向上取整?//kernel.convertTo(kernel,CV_32FC1);//系數未計算return kernel; }/*!* 利用卷積核對輸入圖像進行卷積(濾波)* @param src 輸入.CV_8UC1* @param dst 輸出.CV_8UC1* @param kernel 卷積核(相關核).CV_32FC1. kernel.size()需要為奇數* @note 該函數在OpenCV中的對應函數:cv::filter2D(...)* @author liheng*/ void Filter2D(const cv::Mat& src,cv::Mat& dst,const cv::Mat& kernel) {CV_Assert( kernel.cols % 2 == 1 && kernel.rows % 2 == 1);//邊界填充int border_lr = kernel.cols/2;//左右填充寬度int border_tb = kernel.rows/2;//上下填充寬度cv::Mat borderedSrc;cv::copyMakeBorder(src,borderedSrc,border_tb,border_tb,border_lr,border_lr,cv::BORDER_REFLECT);dst.create(src.size(),src.type());//對圖像中各元素(不包含填充元素)分別進行濾波int rows = borderedSrc.rows;int cols = borderedSrc.cols;for(int row=border_tb; row<rows-border_tb; ++row){for(int col=border_lr; col<cols-border_lr; ++col){//遍歷卷積核鄰域求加權和作為該點輸出值float sum = 0;for(int y=-border_tb; y<=border_tb; ++y)//鄰域高度:2*border_tb+1{for(int x=-border_lr;x<=border_lr; ++x)//鄰域寬度:2*border_lr+1sum += kernel.at<float>(y+border_tb,x+border_lr)*borderedSrc.at<uchar>(row+y,col+x);//權重與元素相乘}//輸出dst.at<uchar>(row-border_tb,col-border_lr) = cv::saturate_cast<uchar>(sum);}} }/*!* 自定義實現高斯濾波.原始版本,未采用分離卷積進行加速* @param src 輸入.CV_8UC1* @param dst 輸出.CV_8UC1* @param ksize 高斯核大小,奇數.當ksize<=0,并且sigma>0時,自動根據sigma計算ksize* @param sigma 高斯濾波標準差.當sigma<=0,并且ksize>0為奇數時,自動根據ksize計算sigma* @author liheng*/ void MyGaussianFilter(const cv::Mat& src,cv::Mat& dst,int ksize,float sigma) {//生成二維高斯核cv::Mat kernel = GetGaussianKernel(ksize,sigma);//濾波Filter2D(src,dst,kernel); }///*!* 生成小數形式的 "一維" 高斯核* @param ksize 高斯核大小,奇數.當ksize<=0,并且sigma>0時,自動根據sigma計算ksize* @param sigma 高斯濾波標準差.當sigma<=0,并且ksize>0為奇數時,自動根據ksize計算sigma* @return 小數形式的高斯核.CV_32FC1.尺寸: 1 * ksize* @author liheng*/ cv::Mat GetSingleDimGaussianKernel( int ksize, double sigma) {//根據sigma自動計算ksizeif( ksize <= 0 && sigma > 0 )ksize = cvRound((sigma*3)*2 + 1)|1;// | 1操作以保證ksize為奇數//校驗ksizeCV_Assert( ksize > 0 && ksize % 2 == 1);//未設定sigma值時,根據窗口大小進行計算sigma = (sigma>0 ? sigma : 0.3f*((ksize-1)*0.5f - 1) + 0.8f);auto sigma2_inverse = 1.0/(2*sigma*sigma);// 1/(2*sigma^2)cv::Mat kernel(1, ksize,CV_32FC1);const int center = ksize/2;//高斯核中心坐標for(int col=0;col<ksize;++col){int x=col-center; //以center為中心,第col列的 x坐標//float f_x = 1.0/(std::sqrt(2*CV_PI)*sigma) * std::exp( -(x*x)/(2*sigma*sigma) );float f_x = std::exp( -(x*x)*sigma2_inverse );//減少不必要計算;系數在后面歸一化時會被消去;kernel.at<float>(0,col) = f_x;}auto sum = cv::sum(kernel)[0];kernel /= sum;return kernel; }/*!* 自定義實現高斯濾波.改進版.利用高斯函數的可分離性以加速卷積過程* @param src 輸入.CV_8UC1* @param dst 輸出.CV_8UC1* @param ksize 高斯核大小,奇數.當ksize<=0,并且sigma>0時,自動根據sigma計算ksize* @param sigma 高斯濾波標準差.當sigma<=0,并且ksize>0為奇數時,自動根據ksize計算sigma* @author liheng*/ void separateGaussianFilter(const cv::Mat& src,cv::Mat& dst,int ksize,float sigma) {//獲取水平和豎直高斯核cv::Mat kernelX = GetSingleDimGaussianKernel(ksize,sigma);cv::Mat kernelY = kernelX.t();cv::Mat Z;//中間變量.存儲水平方向高斯濾波后結果.//水平方向高斯濾波Filter2D(src,Z,kernelX);//對水平方向濾波后結果進行豎直方向高斯濾波Filter2D(Z,dst,kernelY);}int main(int argc, char *argv[]) {cv::Mat kernel = GetGaussianKernel(3,1);std::cout<<"二維高斯核:"<<std::endl;std::cout<<kernel<<std::endl;cv::Mat data = (cv::Mat_<uchar>(3,3)<<10,20,30,40,50,60,70,80,90);//cv::Mat dst;MyGaussianFilter(data,dst,3,1);std::cout<<"自定義實現結果::"<<std::endl;std::cout<<dst<<std::endl;cv::Mat dst1;separateGaussianFilter(data,dst1,3,1);std::cout<<"自定義分離實現結果::"<<std::endl;std::cout<<dst1<<std::endl;cv::Mat dst2;cv::filter2D(data,dst2,CV_8UC1,kernel,cv::Point(-1,-1),0,cv::BORDER_REFLECT);std::cout<<"filter2D結果::"<<std::endl;std::cout<<dst2<<std::endl;cv::GaussianBlur(data,data,cv::Size(3,3),1,1,cv::BORDER_REFLECT);std::cout<<"cv::GaussianBlur結果::"<<std::endl;std::cout<<data<<std::endl;//輸出結果如下:/**二維高斯核:[0.075113609, 0.1238414, 0.075113609;0.1238414, 0.20417996, 0.1238414;0.075113609, 0.1238414, 0.075113609]自定義實現結果::[ 21, 28, 35;43, 50, 57;65, 72, 79]自定義分離實現結果::[ 21, 28, 35;43, 50, 57;65, 72, 79]filter2D結果::[ 21, 28, 35;43, 50, 57;65, 72, 79]cv::GaussianBlur結果::[ 21, 28, 35;43, 50, 57;65, 72, 79]**/return 0; }

從代碼可以看到,其實現原理還是比較簡單的.首先得到一維的高斯核,在卷積（濾波）的過程中,保持行不變,列變化,在水平方向上做卷積運算;接著在上述得到的結果上,保持列不邊,行變化,在豎直方向上做卷積運算. 這樣分解開來,算法的時間復雜度為$O(N)$,運算量和高斯核的尺寸呈線性增長.
利用該分離性,卷積復雜度大大降低,并且我們也不需要生成二維卷積核,而只要生成一維核即可,這也減少了計算量.在實際應用中,甚至還可以直接存儲常用維數的高斯核數值來提升性能.在OpenCV的源碼中就存儲了ksize=1,3,5,7這4中高斯核的權值,當ksize<=7時直接取存儲的值而不必單獨計算一遍高斯核.

總結

高斯函數性質

高斯函數具有五個重要的性質,這些性質使得它在早期圖像處理中特別有用．這些性質表明,高斯平滑濾波器無論在空間域還是在頻率域都是十分有效的低通濾波器,且在實際圖像處理中得到了工程人員的有效使用．高斯函數具有五個十分重要的性質,它們是:
1).二維高斯函數具有旋轉對稱性,即濾波器在各個方向上的平滑程度是相同的．一般來說,一幅圖像的邊緣方向是事先不知道的,因此,在濾波前是無法確定一個方向上比另一方向上需要更多的平滑．旋轉對稱性意味著高斯平滑濾波器在后續邊緣檢測中不會偏向任一方向．
2).高斯函數是單值函數．這表明,高斯濾波器用像素鄰域的加權均值來代替該點的像素值,而每一鄰域像素點權值是隨該點與中心點的距離單調增減的．這一性質是很重要的,因為邊緣是一種圖像局部特征,如果平滑運算對離算子中心很遠的像素點仍然有很大作用,則平滑運算會使圖像失真．
3).高斯函數的傅立葉變換頻譜是單瓣的．如下圖所示,這一性質是高斯函數傅立葉變換等于高斯函數本身這一事實的直接推論．圖像常被不希望的高頻信號所污染(噪聲和細紋理)．而所希望的圖像特征（如邊緣）,既含有低頻分量,又含有高頻分量．高斯函數傅立葉變換的單瓣意味著平滑圖像不會被不需要的高頻信號所污染,同時保留了大部分所需信號．

高斯函數傅立葉變換等于高斯函數本身的證明可自行查找資料

4).高斯濾波器寬度(決定著平滑程度)是由參數σ表征的,而且σ和平滑程度的關系是非常簡單的．σ越大,高斯濾波器的頻帶就越寬,平滑程度就越好．通過調節平滑程度參數σ,可在圖像特征過分模糊(過平滑)與平滑圖像中由于噪聲和細紋理所引起的過多的不希望突變量(欠平滑)之間取得折衷．
5).由于高斯函數的可分離性,較大尺寸的高斯濾波器可以得以有效地實現．二維高斯函數卷積可以分兩步來進行,首先將圖像與一維高斯函數進行卷積,然后將卷積結果與方向垂直的相同一維高斯函數卷積．因此,二維高斯濾波的計算量隨濾波模板寬度成線性增長而不是成平方增長．

高斯濾波應用

高斯濾波后圖像被平滑的程度取決于標準差.它的輸出是領域像素的加權平均,同時離中心越近的像素權重越高.因此,相對于均值濾波(mean filter)它的平滑效果更柔和,而且邊緣保留的也更好.

高斯濾波被用作為平滑濾波器的本質原因是因為它是一個低通濾波器(其頻譜圖如上節所示).而且,大部份基于卷積平滑濾波器都是低通濾波器.

參考資料

1.圖像濾波之高斯濾波介紹
2.圖像處理基礎(4):高斯濾波器詳解
3.簡單易懂的高斯濾波
4.線性濾波、非線性濾波區別
5.如何確定高斯濾波的標準差和窗口大小
6.二維圖像處理中的可分離卷積核
7.可分離卷積的運算量比較與分析（一維、二維卷積）
8.關于高斯濾波的一些理解

附錄

高斯函數及頻譜繪圖代碼

繪制高斯函數曲線代碼:
一維高斯函數:

from matplotlib import pyplot as plt import numpy as npdef gaussian(x, mu, sig):a = 1 / (np.sqrt(2 * np.pi) * sig)return a*np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))x_values = np.linspace(-11, 11, 120) for mu, sig in [(0, 1), (0, 2), (0, 3)]:plt.plot(x_values, gaussian(x_values, mu, sig),label='σ={}'.format(sig))plt.xlabel('x') plt.ylabel('y') plt.legend(loc='best') plt.grid(linestyle='--') plt.show()

二維高斯函數:

import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3Dx,y = np.meshgrid(np.arange(-5,5,0.1),np.arange(-5,5,0.1)) sigma = 1 z = 1/(2 * np.pi * (sigma**2)) * np.exp(-(x**2+y**2)/(2 * sigma**2))fig = plt.figure() ax = Axes3D(fig) ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap='rainbow',alpha = 0.9)ax.set_xlabel("x") ax.set_ylabel("y") ax.set_zlabel("z") # plt.title("二維高斯分布") plt.show()

繪制高斯傅里葉變換頻譜:

# Reference:https://www.cnblogs.com/jingsupo/p/9989559.html import numpy as np#導入一個數據處理模塊 import pylab as pl#導入一個繪圖模塊,matplotlib下的模塊def gaussian(x, mu, sig):a = 1 / (np.sqrt(2 * np.pi) * sig)return a*np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))sampling_rate = 8000#采樣頻率為8KHz fft_size = 8000 #FFT處理的取樣長度 t = np.arange(0, 4.0, 1.0/sampling_rate)#np.arange(起點,終點,間隔)產生4s長的取樣時間 x = gaussian(t,0,1) xs = x[:fft_size]# 從波形數據中取樣fft_size個點進行運算 xf = np.fft.rfft(xs)# 利用np.fft.rfft()進行FFT計算,rfft()是為了更方便對實數信號進行變換,由公式可知/fft_size為了正確顯示波形能量 xf = xf / len(xf) # rfft函數的返回值是N/2+1個復數,分別表示從0(Hz)到sampling_rate/2(Hz)的分. #于是可以通過下面的np.linspace計算出返回值中每個下標對應的真正的頻率: freqs = np.linspace(0, sampling_rate/2, fft_size//2+1)#頻率 xfp = np.abs(xf)#幅值#繪圖顯示結果 pl.figure(figsize=(8,4)) pl.subplot(211) pl.plot(t, x) pl.xlim(0,) pl.ylim(0,) pl.xlabel('t') pl.title('Gaussias\'s WaveForm And Freq')pl.subplot(212) freqs = freqs[:100] xfp = xfp[:100] _range = np.max(xfp) - np.min(xfp) xfp=(xfp - np.min(xfp)) / _rangepl.plot(freqs, xfp) pl.xlim(0,) pl.ylim(0,) pl.xlabel('Freq/Hz') pl.ylabel('Amplitude') pl.subplots_adjust(hspace=0.4) pl.show()

---

下面的是我的公眾號二維碼圖片，按需關注。

總結

以上是生活随笔為你收集整理的高斯滤波的理解与学习的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。