高斯濾波及其原理

一、高斯函數的基礎

1.1 期望、方差與標準差

用來刻畫隨機變量某一方面特征的常數被稱為隨機變量的數字特征，其常用的有：

數學期望：

? 在概率論和統計學中，數學期望（mean）是試驗中每次可能結果的概率乘以其結果的總和，是最基本的數學特征之一。它反映隨機變量平均取值的大小。

? 需要注意的是，期望值并不一定等同于常識中的“期望”——“期望值”也許與每一個結果都不相等。期望值是該變量輸出值的平均數。期望值并不一定包含于變量的輸出值集合里。大數定律表明，隨著重復次數接近無窮大，數值的算術平均值幾乎肯定地收斂于期望值。舉個簡單的例子：擲骰子游戲中，玩家搖到紅色（1和4）可以贏得90元，搖到藍色（2、3、5、6）輸掉60元。而搖到紅色的概率為1/3，藍色的概率為2/3，也就是說玩家有1/3的期望獲得90元，而商家有2/3的 期望獲得60元，客觀來說一次游戲玩家獲得的金錢為30，而商家獲得40元。

? 對于離散型（隨機變量只取得有限個值或無窮能按次序一一列出）的隨機變量的一切可能取值$x$_i與對應的概率$p$($x$_i)乘積之和稱為該離散型隨機變量的數學期望（若該求和絕對收斂），記作$E(x)$或者$\mu$，他是簡繁算數平均的一種推廣， 類似加權平均，即$\mu ={\sum }_{k=1}^{\infty }{x}_k{p}_k$。離散的才有均值。

? 對于連續型隨機變量$x$的概率密度函數為$f(x)$，若積分絕對收斂，則積分的值${\int }_{?\infty }^{\infty }xf(x)dx$為隨機變量的數學期望，即$\mu ={\int }_{?\infty }^{\infty }xf(x)dx$。連續的有數學期望可是沒有均值。

方差：

? 方差是在概率論和統計方差衡量隨機變量或一組數據時離散程度的度量。概率論中方差用來度量隨機變量和其數學期望（即均值）之間的偏離程度。統計中的方差（樣本方差）是每個樣本值與全體樣本值的平均數之差的平方值的平均數。在許多實際問題中，研究方差即偏離程度有著重要意義。簡單來說就是描述一個變量離其期望值的距離，即各個誤差的平方累加，再除以總數，即對于離散型隨機變量${\sigma }^2=\frac{\sum (X?\mu {)}^2}{N}$，其中${\sigma }^2$為總體方差，$X$為變量，$\mu$為總體均值，$N$為總體例數。而對于連續型隨機變量${\sigma }^2=\int (x?\mu {)}^{2}f(x)dx$。

標準差：

? 標準差是方差的算術平方根，用$\sigma$表示，標準差也被稱為標準偏差，反映了一個數據集的離散程度。如下圖（概率密度函數$f(x)$）通過標準差可以直觀得到距離平均值$\mu$到$\mu +\sigma$之間的概率。注意，概率密度函數是對連續隨機變量定義的，且它本身并不是概率，只有對連續隨機變量的概率密度函數在某區間內進行積分后才能得到概率，對于一個連續型隨機變量而言，他取任何固定值的概率都為0，也就是說考察隨機變量在某一點的概率取值是沒有意義的。因此，在考察連續型隨機變量的分布時，討論的是它在某個區間的概率取值。這時就要借助其累積分布函數($F(X)$)。

補充：概率分布函數PDF和累積分布函數CDF

累積分布函數（Cumulative Distribution Function, CDF），簡稱分布函數，對于隨機變量$X$，如下定義的函數F
$$F(X)=P(X\le x)?\infty <x<\infty$$
稱為$X$的累積分布函數，對于任意給定的實數$X$，分布函數等于該隨機變量小于等于$x$的概率。

概率分布函數（Probability Density Function, PDF）:對于連續型隨機變量$X$的累積分布函數$F(x)$,如果存在一個定義在$x$軸上的非負函數$f(x)$，使得對于任意實數$x$，有下式成立
$$F(X)={\int }_{?\infty }^{x}f(t)dt$$
則稱$f(x)$為$X$的概率密度函數，顯然，當概率密度函數存在的時候，累積分布函數是概率密度函數的積分，兩者的圖像如下所示：

1.2 一維高斯函數

? 高斯函數（Gaussian Function），也簡稱為Gaussian，是以數學家Carl Friedrich Gaussian的名字命名的。其一維形式如下：
$$\begin{gathered} f(x)=a{e}^{?\frac{(x?b{)}^2}{2c^2}} \\ (a>0) \end{gathered}$$

? 高斯函數的一維圖像是對稱的鐘形（bell curve），以上公式中各個字符在圖像中的含義為：a是曲線尖峰的高度，b是尖峰中心的坐標，c稱為標準方差，表征的是bell鐘狀的寬度。高斯函數廣泛應用于模糊、正態分布統計、卷積等方面，具體來說可以做熒光筆效果、毛玻璃效果、邊緣平滑、抗鋸齒等。

? 由高斯積分公式${\int }_{?\infty }^{\infty }{e}^{?x^2}dx=\sqrt{\pi }$，令$y=x?b$，可得$dx=dy$，令$z^2=\frac{y^2}{2c^2}$，可得$dy=\sqrt{2}cdz$。

? 則
$${\int }_{?\infty }^{\infty }a{e}^{\frac{?(x?b{)}^2}{2c^2}}dx$$
? 當且僅當$a=\frac{1}{c\sqrt{2\pi }}$時上式積分值為1，在這種情況下，高斯是正態分布隨機變量的概率密度函數。而正態分布就是使$b=\mu ,c=\sigma$，此時:
$$f(x)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{?\frac{(x?\mu {)}^2}{{\sigma }^2}}$$
? 當$\mu =0,{\sigma }^2=1$時為標準正態分布，此時：
$$f(x)=\frac{1}{\sqrt{2\pi }}{e}^{?\frac{x^2}{{\sigma }^2}}$$

1.2二維高斯函數

二維高斯函數形式為：
$$f(x,y)=A{e}^{?(\frac{(x?x_0{)}^2}{2{\sigma }_x^2}+\frac{(y?y_0{)}^2}{2{\sigma }_y^2})}$$
其中A是幅值，$x_0$,$y_0$是中心點坐標，${\sigma }_x^2和{\sigma }_y^2$是方差,圖示如下圖所示，$A=1$, $x_0=y_0=0$，${\sigma }_x^2={\sigma }_y^2=1$

同樣，類似于一維高斯函數，得到的高斯函數如下：
$$G(x,y)=\frac{1}{2\pi {\sigma }^2}{e}^{?\frac{x^2+y^2}{2{\sigma }^2}}$$

二、高斯濾波原理

2.1高斯圖像分析

使用matlab直觀的分析高斯圖像，在實際的matlab編程中，高斯函數用到的參數有：

· ksize 高斯函數的大小

· sigma 高斯函數的方差

· center 高斯函數尖峰中心點坐標

· **bias ** 高斯函數尖峰中心點的偏移量，用于控制階段高斯函數

為了方便直觀的觀察高斯函數不同參屬下的結果圖像，使用代碼實現參數的自動遞增，并且生成gif圖像，完整代碼如下：

function mainfunc() % 測試高斯函數，遞增的方法實現高斯函數參數的改變對整個高斯函數的影響， % 并自動保存為gif格式輸出。 % created by zhao.buaa 2016.09.28 % 引用https://blog.csdn.net/qinglongzhan/article/details/82348153中的代碼，修改了視角，輸入mainfunc();即可運行%% 保存gif動畫 item = 10; % 迭代次數 dt = 1; % 步長大小 ksize =20; % 高斯大小 sigma = 2; % 方差大小 filename = ['ksize-' num2str(ksize) '-sigma-' num2str(sigma) '--' num2str(sigma+dt*item) '.gif']; %必須預先建立gif文件 % main loop for i = 1:item center = round(ksize/2); % 中心點 bias = ksize*10/10; % 偏移中心點量 ksigma = ksigma(ksize, sigma, center, bias); % 構建核函數 tname = ['ksize-' num2str(ksize) '-sigma-' num2str(sigma) '-center-' num2str(center)]; figure(i), mesh(ksigma), title(tname); %設置固定的x-y-z坐標范圍，便于觀察，axis([xmin xmax ymin ymax zmin zmax]) axis([0 ksize 0 ksize 0 0.008]); view([30, 30]);% 改變可視角度 sigma = sigma + dt; % 自動保存為gif圖像 frame = getframe(i); im = frame2im(frame); [I,map] = rgb2ind(im,256); if i==1 imwrite(I,map,filename,'gif','Loopcount',inf, 'DelayTime',0.4); else imwrite(I,map,filename,'gif','WriteMode','append','DelayTime',0.4); end end; close all; %% 截斷高斯核函數，截斷的程度取決于參數bias function ksigma = ksigma(ksize, sigma, center,bias) %ksize = 80; sigma = 15; ksigma=fspecial('gaussian',ksize, sigma); % 構建高斯函數 [m, n] =size(ksigma); for i = 1:m for j = 1:n if( (i<center-bias)||(i>center+bias)||(j<center-bias)||(j>center+bias) ) ksigma(i,j) = 0; end; end; end;

以上主要用到了fspecial這個matlab函數，該函數的官方解釋為：

H = fspecial(‘gaussian’,HSIZE,SIGMA) returns a rotationally symmetric Gaussian lowpass filter of size HSIZE with standard deviation SIGMA (positive). HSIZE can be a vector specifying the number of rows and columns in H or a scalar, in which case H is a square matrix.

即高斯低通濾波，有兩個參數，HSIZE表示模板尺寸，sigma為濾波器的標準差，可以看到$\sigma$越大，鐘形越胖（同一維高斯函數），同時整個高斯函數的尖峰逐漸減小，整體也更加平滑，也就是對圖像的平滑效果也越明顯。

保持其他參數不變，對上述高斯函數進行截斷，即將上述代碼中的$bias=ksize?10/10$;改為$bias=ksize?3/10$，則結果如下圖，也就是說對超過一定區域的原始圖像信息不再考慮（可能類似于高斯模板有大有小），從而保證在更加合理的利用靠近高斯函數中心點的周圍像素。

2.2高斯核與高斯函數的關系以及計算

以上關于高斯圖像的討論，并沒有解釋高斯模板的由來以及和高斯濾波的關系，其實高斯濾波用的就是高斯函數圖像的特性進行的。簡單來說就是高斯濾波器是一類根據高斯函數的形狀來選擇權值的線性平滑濾波器。在對圖像進行卷積的時候其實是對圖像進行點操作，用所求像素點周圍的像素通過模板的權值計算得到新的像素值。而圖像形狀的不同也就代表了不同的參數，從而得到不同的模板。還是以上面程序中提到的fspecial函數為例，看看怎么計算高斯模板。

常用零均值離散高斯濾波器函數進行計算，上面提到二維高斯函數的形式為$G(x,y)=\frac{1}{2\pi {\sigma }^2}{e}^{?\frac{x^2+y^2}{2{\sigma }^2}}$(曲線積分面積為1)，而在圖像上是離散的點，將該函數離散化得到：
$$H(i,j)=\frac{1}{2\pi {\sigma }^2}{e}^{?\frac{(i?k?1{)}^2+(j?k?1{)}^2}{2{\sigma }^2}}$$
高斯卷積核H:(2k+1)×(2k+1),$\sigma$為方差。

驗證fspecial的計算：在matlab中輸入$k=fspecial{(}^{\prime }gaussia{n}^{\prime },3,1);$，得到的kernel如下：

對于上述公式，取k=1,$\sigma =1$,即可得到3×3的高斯卷積核如下：
KaTeX parse error: Undefined control sequence: \matrix at position 13: H = \left[ \?m?a?t?r?i?x?{ \frac{1}{2\…
然而計算出的高斯卷積核還應該進一步歸一化，也就是計算出來的高斯模板中各個數值其和必須為1，因為不歸一化計算出來系的像素值會偏離原來的像素范圍。歸一化如下所示：
KaTeX parse error: Undefined control sequence: \matrix at position 13: H = \left[ \?m?a?t?r?i?x?{ 0.0585 & 0.…

2.3 高斯核與高斯模板的關系

為了區別，暫且把以上計算出來的稱為高斯核吧，而我們經常會用到一些形如以下的高斯模板，那么這個高斯模板和高斯核有什么關系呢，高斯模板實際上也就是模擬高斯函數的特征，具有對稱性并且數值由中心向四周不斷減小，這個模板剛好符合這樣的特性，并且非常簡單，容易被大家接受，于是就比較經典！
KaTeX parse error: Undefined control sequence: \matrix at position 21: …c{1}{16}\left[ \?m?a?t?r?i?x?{ 1 & 2 & 1\\…

2.4 OpenCV中cvSmooth函數的用法

? 但是以上這樣一個簡單的矩陣是遠遠不能滿足我們對圖像處理的要求的，我們需要按照自己的要求得到更加精確的模板，那么接下來我們先看看OpenCV中高斯核的計算，緊接著再進一步編程實現自己想要的高斯模板。

OpenCV中GaussianBlur的原型為：

void cv::GaussianBlur ( InputArray src, OutputArray dst, Size ksize, double sigmaX, double sigmaY = 0, int borderType = BORDER_DEFAULT )

關于該函數的官方文檔如下：

https://docs.opencv.org/4.x/d4/d86/group__imgproc__filter.html#gaabe8c836e97159a9193fb0b11ac52cf1

Blurs an image using a Gaussian filter.

The function convolves the source image with the specified Gaussian kernel. In-place filtering is supported.

Parameters
src input image; the image can have any number of channels, which are processed independently, but the depth should be CV_8U, CV_16U, CV_16S, CV_32F or CV_64F.
dst output image of the same size and type as src.
**ksize ** Gaussian kernel size. ksize.width and ksize.height can differ but they both must be positive and odd. Or, they can be zero’s and then they are computed from sigma.
sigmaX Gaussian kernel standard deviation in X direction.
sigmaY Gaussian kernel standard deviation in Y direction; if sigmaY is zero, it is set to be equal to sigmaX, if both sigmas are zeros, they are computed from ksize.width and ksize.height, respectively (see getGaussianKernel for details); to fully control the result regardless of possible future modifications of all this semantics, it is recommended to specify all of ksize, sigmaX, and sigmaY.
borderType pixel extrapolation method, see BorderTypes. BORDER_WRAP is not supported.

也就是說，高斯核的Size必須是正奇數或者也可以是0，當為0時，可以通過sigma進行計算，另外sigmaX和sigmaY分別代表X、Y方向的均方差。當sigmaY=0時，其值與sigmaY相同，而當兩者皆沒有給出時，可以通過Ksize進行計算。其中用到了getGaussianKernel()函數，而getGaussianKernel函數原型為：

Mat cv::getGaussianKernel ( int ksize, double sigma, int ktype = CV_64F )

同樣分析其官方文檔：

Returns Gaussian filter coefficients.

The function computes and returns the (\texttt{ksize} \times 1) matrix of Gaussian filter coefficients:

[G_i= \alpha e^{-(i-( \texttt{ksize} -1)/2)^2/(2 \texttt{sigma}^2)},]

where (i=0…\texttt{ksize}-1) and (\alpha) is the scale factor chosen so that (\sum_i G_i=1).

Two of such generated kernels can be passed to sepFilter2D. Those functions automatically recognize smoothing kernels (a symmetrical kernel with sum of weights equal to 1) and handle them accordingly. You may also use the higher-level GaussianBlur.

• Parameters

  ksizeAperture size. It should be odd ( (\texttt{ksize} \mod 2 = 1) ) and positive.sigmaGaussian standard deviation. If it is non-positive, it is computed from ksize as sigma = 0.3*((ksize-1)*0.5 - 1) + 0.8.ktypeType of filter coefficients. It can be CV_32F or CV_64F .

簡單總結一下就是根據ksize去計算$\sigma$的大小，根據公式$\sigma =0.3((ksize?1)\times 0.5?1)+0.8$，反過來也可以根據$\sigma$去計算核的大小。這里根據的是高斯分布的特點（如圖所示，數值分布在（μ—3σ,μ+3σ)中的概率為0.9974），如果核矩陣更大，那么相應的Sigma也更大，相反，如果Sigma更大，那么核矩陣覆蓋范圍也更大。

2.5 編寫計算高斯模板的自定義函數

編寫了一個自定義的計算高斯模板的函數如下，在$\sigma$缺省（等于0）時，通過OpenCV中的公式計算$\sigma$值，同樣，在高斯模板尺寸缺省時，通過該公式由$\sigma$反過來計算模板大小，完整代碼如下：

double ** createGuassion(int , double); double ** createGuassion(int kSize, double sigma) { if (kSize%2==0 || kSize<0 || (kSize==0&&sigma==0)){return 0;//當模板大小為偶數、或者小于0、或者模板大小與sigma同時為0時返回}else if (kSize==0){kSize = 2 * (sigma - 0.8) / 0.3 + 3;kSize = round(kSize);//if (kSize % 2 == 0){kSize += 1;}}else if (sigma == 0){sigma = 0.3*((kSize - 1)*0.5 - 1) + 0.8;//sigma缺省時，利用該公式計算sigma}double **guass;//生成的高斯矩陣double sum = 0;//用于歸一化求和double x2 = 0;double y2 = 0;int center = (kSize - 1) / 2;//高斯核中心點的行標和列標（從0開始計），所以必須保證kSize為奇數guass = new double*[kSize];//注意，double*[k]表示一個有10個元素的指針數組for (int i = 0; i < kSize; ++i){//先創建一個kSize×kSize的模板矩陣guass[i] = new double[kSize];}for (int i = 0; i < kSize; i++){x2 = pow(double(i - center), 2); //高斯函數中的x平方for (int j = 0; j < kSize; j++){y2 = pow(double(j - center), 2);//高斯函數中的y平方sum += guass[i][j] = exp(-(x2 + y2) / (2 * sigma*sigma));//賦值并累加}}//歸一化處理if (sum != 0){for (int i = 0; i < kSize; i++){for (int j = 0; j < kSize; j++){guass[i][j] /= sum;}}}return guass; }int main() {double **p = createGuassion(3, 1);for (int i = 0; i < 3; ++i){for (int j = 0; j < 3; ++j){printf("%-9f\x20", *(*(p + i) + j));//p+i表示a[i]的地址//"%-9f中, '-'表示左對齊(無'-'默認為右對齊)；9表示這個元素輸出時占兩個空格的空間//"\x20"表示以十六進制數20所代表的字符，即空格（space的ASCII為32）}printf("\n");}printf("\n");for (int i = 0; i < 3; ++i){for (int j = 0; j < 3; ++j){printf("%7.4f\x20", *(*(p + i) + j)); //p+i表示a[i]的地址//"%-7.4f中, '-'表示左對齊(無'-'默認為右對齊)；//7.4表示這個元素輸出時占7列，其中包括4個小數}printf("\n");}printf("\n");return 0;

運行結果如下，上面是原輸出，當取計算的4位小數時，結果與matlab中計算的一模一樣，由此也驗證了該函數的正確性。

0.075114 0.123841 0.075114 0.123841 0.204180 0.123841 0.075114 0.123841 0.0751140.0751 0.1238 0.07510.1238 0.2042 0.12380.0751 0.1238 0.0751

另外分析在使用創建的高斯模板進行卷積的時候，二維函數處理應該是這樣：

for(i=0;i<img.height;i++)for(j=0;j<img.width;j++)for(m=0;m<kSize;kSize++)for(n=0;n<kSize;kSize++)......

這樣卷積算法的復雜度應該是O（height*width*kSize*kSize）

然而高斯函數具有可分離性，也就是說一個二維的高斯函數可以分解成相同的一維高斯函數，用其來處理兩遍圖像，效果一樣，但是算法復雜度就變成了O（height*width*kSize*2），這樣算法復雜度就會減小不少。

相應的模板函數變為：

//生成一維高斯模板，水平的和垂直方向上的模板是一樣的,與上面的二維類似，只不過少計算一維 double* CreateMuban(int kSize ,double sigma) {double *gauss = new double[kSize];//聲明一維模板int radius = (kSize - 1) / 2;//這是高斯中心double sum = 0;for (int i = 0; i < kSize; i++){//高斯函數前面的常數項因為在歸一化的時候將會消去，故這里不重復計算gauss[i] = exp(-(i - radius)*(i - radius)/(2*sigma*sigma));sum = sum + gauss[i];}for (int i = 0; i < kSize; i++){//歸一化gauss[i] = gauss[i] / sum;}return gauss; } //這里得到得guassian模板就是一維的，用來處理圖像時，要分別進行兩次操作，即對水平和垂直方向分別進行卷積。

2.6 邊界點的處理

應該注意到在OpenCV中的GuassianBlur中有一個參數為BorderType，這就涉及到邊界點的處理，因為卷積中如果一個像素點處于邊界，周邊沒有足夠的點，這個時候應該怎么處理？常見的方法要么是補0，要么就是把已有的點拷貝到另一點的對應位置，模擬出完整的矩陣。邊界點的處理以及高斯的分布處理將放到卷積中總結。

參考文獻：

百度百科

https://zhuanlan.zhihu.com/p/40060966

https://www.pianshen.com/article/1236562010/

https://blog.csdn.net/qinglongzhan/article/details/82348153

https://www.cnblogs.com/ping-36/articles/2810203.html

https://blog.csdn.net/qinglongzhan/article/details/82348153

https://docs.opencv.org/4.x/d4/d86/group__imgproc__filter.html#gaabe8c836e97159a9193fb0b11ac52cf1

https://www.cnblogs.com/bingdaocaihong/p/7007346.html

https://blog.csdn.net/qinglongzhan/article/details/82348153

總結

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

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