高斯濾波在影象處理中的應用

阿新 • • 發佈：2018-12-09

卷積：

相信很多時候，當我們在看到“卷積”時，總是處於一臉懵逼的狀態，不但因為它的本義概念比較難理解，還因為它在不同的應用中發揮出的變幻莫測的作用也時常讓人迷糊。但這些應用其實本質上都是同一種東西，理解了卷積的來源，就可以舉一反三。其實我個人對於卷積的理解，很長時間都處於似懂非懂的狀態，就像傅立葉變換的一些tricky points，只求在應用中不出差錯，不求甚解。但是如果想要真正做好learning的東西，我認為在真正的學習理論基礎前，必須將概念的本質搞清楚。因此，這篇文章便是為工科，理科以及其他領域對此感興趣的同學而整理，內容簡易通俗，便於理解。

什麼是卷積

卷積，很多時候都是我們在各種工程領域，訊號領域所看到的常用名詞，比如系統

通俗易懂的說，就是

輸出 = 輸入 * 系統

雖然它看起來只是個複雜的數學公式，但是卻有著重要的物理意義，因為自然界這樣的系統無處不在，計算一個系統的輸出最好的方法就是運用卷積。更一般的，我們還有很多其他領域的應用：

統計學中，加權的滑動平均是一種卷積。

概率論中，兩個統計獨立變數X與Y的和的概率密度函式是X與Y的概率密度函式的卷積。

聲學中，回聲可以用源聲與一個反映各種反射效應的函式的卷積表示。

電子工程與訊號處理中，任一個線性系統的輸出都可以通過將輸入訊號與系統函式（系統的衝激響應）做卷積獲得。

物理學中，任何一個線性系統（符合疊加原理）都存在卷積。

電腦科學

中，卷積神經網路（CNN）是深度學習演算法中的一種，近年來被廣泛用到模式識別、影象處理等領域中。

這6個領域中，卷積起到了至關重要的作用。在面對一些複雜情況時，作為一種強有力的處理方法，卷積給出了簡單卻有效的輸出。對於機器學習領域，尤其是深度學習，最著名的CNN卷積神經網路(Convolutional Neural Network, CNN)，在影象領域取得了非常好的實際效果，始一出現便橫掃各類演算法。

那麼，到底什麼是卷積呢？

首先給大家一個果殼上關於卷積的著名暴力講解：

比如說你的老闆命令你幹活，你卻到樓下打檯球去了，後來被老闆發現，他非常氣憤，扇了你一巴掌（注意，這就是輸入訊號，脈衝），於是你的臉上會漸漸地（賤賤地）鼓起來一個包，你的臉就是一個系統，而鼓起來的包就是你的臉對巴掌的響應，好，這樣就和訊號系統建立起來意義對應的聯絡。下面還需要一些假設來保證論證的嚴謹：假定你的臉是線性時不變系統，也就是說，無論什麼時候老闆打你一巴掌，打在你臉的同一位置（這似乎要求你的臉足夠光滑，如果你說你長了很多青春痘，甚至整個臉皮處處連續處處不可導，那難度太大了，我就無話可說了哈哈），你的臉上總是會在相同的時間間隔內鼓起來一個相同高度的包來，並且假定以鼓起來的包的大小作為系統輸出。好了，那麼，下面可以進入核心內容——卷積了！如果你每天都到地下去打檯球，那麼老闆每天都要扇你一巴掌，不過當老闆打你一巴掌後，你5分鐘就消腫了，所以時間長了，你甚至就適應這種生活了……如果有一天，老闆忍無可忍，以0.5秒的間隔開始不間斷的扇你的過程，這樣問題就來了，第一次扇你鼓起來的包還沒消腫，第二個巴掌就來了，你臉上的包就可能鼓起來兩倍高，老闆不斷扇你，脈衝不斷作用在你臉上，效果不斷疊加了，這樣這些效果就可以求和了，結果就是你臉上的包的高度隨時間變化的一個函數了（注意理解）；如果老闆再狠一點，頻率越來越高，以至於你都辨別不清時間間隔了，那麼，求和就變成積分了。可以這樣理解，在這個過程中的某一固定的時刻，你的臉上的包的鼓起程度和什麼有關呢？和之前每次打你都有關！但是各次的貢獻是不一樣的，越早打的巴掌，貢獻越小，所以這就是說，某一時刻的輸出是之前很多次輸入乘以各自的衰減係數之後的疊加而形成某一點的輸出，然後再把不同時刻的輸出點放在一起，形成一個函式，這就是卷積，卷積之後的函式就是你臉上的包的大小隨時間變化的函式。本來你的包幾分鐘就可以消腫，可是如果連續打，幾個小時也消不了腫了，這難道不是一種平滑過程麼？反映到劍橋大學的公式上，f(a)就是第a個巴掌，g(x-a)就是第a個巴掌在x時刻的作用程度，乘起來再疊加就ok了通過這個通俗化的例子我們從基本概念上了解了卷積，那麼更嚴格的定義是怎樣的呢？

從數學上講，卷積只不過是一種運算，對於很多沒有學過訊號處理，自動控制的同學來說各種專業的名詞可以不做了解。我們接著繼續：

本質上卷積是將二元函式 U(x,y) = f(x)g(y) 捲成一元函式 V(t) ，俗稱降維打擊。

怎麼卷？

考慮到函式 f 和 g 應該地位平等，或者說變數 x 和 y 應該地位平等，一種可取的辦法就是沿直線 x+y = t 捲起來：

$V(t) = \int_{x+y=t} U(x,y) \,\mathrm{d}x$

捲了有什麼用？

可以用來做多位數乘法，比如：

$\begin{align}42 \times137 &= (2\times10^0+4\times10^1)(7\times10^0+3\times10^1+1\times10^2) \\&= (2\times7)\times10^0 + (2\times3+4\times7)\times10^1+(2\times1+4\times3)\times10^2 + (4\times1)\times10^3 \\&= 14 + 340+1400+4000 \\&= 5754\end{align}$

注意第二個等號右邊每個括號裡的係數構成的序列 (14,34,14,4)，實際上就是序列 (2,4) 和 (7,3,1) 的卷積

在乘數不大時這麼幹顯得有點蛋疼，不過要計算很長很長的兩個數乘積的話，這種處理方法就能派上用場了，因為你可以用快速傅立葉變換 FFT 來得到卷積，比示例裡的硬乘要快。

這裡有一個不太嚴格的理解： $(\sum_{n=1}^{\infty}{a_nx^n})(\sum_{n=1}^{\infty}{b_nx^n})=\sum_{n=1}^{\infty}(\sum_{k=1}^{n}a_kb_{n-k})x^n$ $x^n$ 是“基”， $a_n$ 是在這個基上的展開係數。兩個多項式乘積的在基上展開的係數就是兩個多項式各自在基上展開係數的卷積。 $x^n$ 對應著頻率不同的 $\exp(ikt)$ ，係數對應著其傅立葉變換。自然就是乘積的傅立葉變換等於傅立葉變換的卷積了。

卷積的核心（涉及推導過程，可以跳過）：

首先我們有這樣一個概念：內積、積分、投影這三者其實從某個角度上講是一個意思。

定義一組向量 $\alpha=(\alpha_{1},\alpha_{2},...,\alpha_{n})$ ，另一組向量 $\beta=(\beta_{1},\beta_{2},...,\beta_{n})$ ，那麼內積可以表達為：

$\alpha \beta = \alpha_{1} \beta_{1}+\alpha_{2} \beta_{2}+...+\alpha_{n} \beta_{n}=\sum_{i=1}^{n}{\alpha_{i} \beta_{i}}$

這即是內積，也是累加（積分）。投影的概念則可以理解為向量a在基向量b上的一組投影，座標為 $(\alpha_{1},\alpha_{2},...,\alpha_{n})$ ，這和一個點在3D歐幾里得空間的三軸投影座標是一個道理。

這樣，先來看看Fourier變換在做什麼：

$F(\omega)=\int_{-\infty }^{+\infty} f(t)e^{-j\omega t}dt$

再引入一個完美的式子，尤拉公式：

$e^{ix}=cosx+isinx$

從Fourier的定義式可以看出是對f(t)和 $e^{-j\omega t}$ 相乘後在無窮域上對其進行積分，那麼其實就是將f(t)投影在 $e^{-j\omega t}$ 上，如果不理解 $e^{-j\omega t}$ 就變換為兩個正交的三角函式（尤拉公式就在這裡起作用）。所以這就明朗了：Fourier把f(t)投影到了兩個為正交關係的正弦和餘弦空間中。也可以從週期訊號的Fourier級數分解表示式更容易看出這個投影關係。

看完Fourier再看來控制論領域的Laplace變換在做什麼：

$F(s)=\int_{0}^{\infty}f(t)e^{-\sigma t}e^{-j\omega t}dt=\int_{0}^{\infty}f(t)e^{-st}dt$

首先，控制領域裡面經常用到階躍訊號，不幸的是它不滿足狄利克雷第三條件，因此它對Fourier變換免疫，所以聰明的Laplace用了一個衰減因子 $e^{-\sigma t}$ 將其進行衰減後再做Fourier變換。到了負無窮的區域這衰減因子可就成了遞增因子，所以Laplace變換僅限於大於0的區域，對於小於0的區域用系統初始狀態表達就好了。從這點角度上講，Laplace變換相當於對 $f(t)e^{-\sigma t}$ 做了一個單邊Fourier變換。

然後，分析方法同上，可以看到Laplace把f(t)投影到了 $e^{-st}$ 空間，這就是s平面。它比Fourier更厲害的地方是不僅可以看到虛軸上 $\omega$ 的成分，還可以在實軸上看到Fourier看不到的衰減因子 $\sigma$ 成分，這是Fourier做不到的。所以Laplace在Fourier的基礎上把訊號拓展到了衰減因子實軸上，這個衰減因子 $\sigma$ 和系統的阻尼 $\zeta$ ，自然震盪角頻率 $\omega_{n}$ 密切相關，直接影響了系統的調節時間 $t_{s}$ 。學過自控原理的同學應該知道在頻域章節，我們得到系統的頻域響應曲線都是通過傳遞函式來直接轉化的，公式就是s=jw。這也就是說Fourier活在一維虛軸空間，Laplace活在二維平面空間，想要得到一維空間上關於w的表達形式，只需要在s平面上做降維處理即可。

回過頭來再來看看卷積投影：

$f(t)*g(t)=\int_{-\infty }^{+\infty}f(\tau)g(t-\tau)d\tau$

這個投影有點奇怪，它在投影之前先把 $g(\tau)$ 做了一個反對稱，然後再投影。對應到前面推導的系統卷積表示式：

$r(t)*f(t)=\int_{-\infty }^{+\infty}f(\tau)r(t-\tau)d\tau$

相當於在投影之前，先把輸入訊號r(t)在時間軸上翻轉了180°，然後與系統f(t)進行投影。投影的概念我們可以很好理解，無論是向量內積運算相當於線投影，或者空間的一個多面體在三維空間平面上的投影面，這種投影運算就相當於一種重合面積。如果從這個角度去看輸入、系統和輸出三者之間的關係，那麼就可以從圖形角度去理解為什麼一個一階系統在階躍響應輸出下是一條單調上升的曲線了。這裡用一張wikipedia裡關於卷積的一張圖形化解釋，想要了解更多的同學可以自行跳轉：Convolution

（特此感謝知乎學霸王尼莫的幫助）

卷積的應用

影象處理：用一個模板和一幅影象進行卷積，對於影象上的一個點，讓模板的原點和該點重合，然後模板上的點和影象上對應的點相乘，然後各點的積相加，就得到該點的卷積值。對影象上的每個點都這樣處理。由於多數模板都對稱，所以模板不旋轉。卷積是一種積分運算，用來求兩個曲線重疊區域面積。可以看作加權求和，可以用來消除噪聲、特徵增強。把一個點的畫素值用它周圍的點的畫素值的加權平均代替。卷積是一種線性運算,影象處理中常見的mask運算都是卷積，廣泛應用於影象濾波。卷積關係最重要的一種情況，就是在訊號與線性系統或數字訊號處理中的卷積定理。利用該定理，可以將時間域或空間域中的卷積運算等價為頻率域的相乘運算，從而利用FFT等快速演算法，實現有效的計算，節省運算代價。

下面是來自sselssbh部落格的一個例子，非常形象的解釋了卷積在影象領域的作用

有這麼一副影象，可以看到，影象上有很多噪點：這裡寫圖片描述

高頻訊號，就好像平地聳立的山峰：這裡寫圖片描述

看起來很顯眼。

平滑這座山峰的辦法之一就是，把山峰刨掉一些土，填到山峰周圍去。用數學的話來說，就是把山峰周圍的高度平均一下。

平滑後得到：這裡寫圖片描述

4.2 計算

卷積可以幫助實現這個平滑演算法。

有噪點的原圖，可以把它轉為一個矩陣：這裡寫圖片描述

然後用下面這個平均矩陣（說明下，原圖的處理實際上用的是正態分佈矩陣，這裡為了簡單，就用了算術平均矩陣）來平滑影象：

記得剛才說過的演算法，把高頻訊號與周圍的數值平均一下就可以平滑山峰。

比如我要平滑點，就在矩陣中，取出點附近的點組成矩陣 f ,和 g 進行卷積計算後，再填回去這裡寫圖片描述

要注意一點，為了運用卷積， g 雖然和 f 同維度，但下標有點不一樣：這裡寫圖片描述

這裡寫圖片描述

寫成卷積公式就是：

要求，一樣可以套用上面的卷積公式。

這樣相當於實現了 g 這個矩陣在原來影象上的划動（準確來說，下面這幅圖把 g 矩陣旋轉了）：

再比如做饅頭

樓下早點鋪子生意太好了，供不應求，就買了一臺機器，不斷的生產饅頭。

假設饅頭的生產速度是 f(t) ，那麼一天後生產出來的饅頭總量為： ∫240f(t)dt

饅頭生產出來之後，就會慢慢腐敗，假設腐敗函式為 g(t) ，比如，10個饅頭，24小時會腐敗： 10∗g(t)

想想就知道，第一個小時生產出來的饅頭，一天後會經歷24小時的腐敗，第二個小時生產出來的饅頭，一天後會經歷23小時的腐敗。如此，我們可以知道，一天後，饅頭總共腐敗了： ∫240f(t)g(24−t)dt 文章也發表在我的個人部落格中：點選開啟連結，更多與機器學習，數學相關知乎，歡迎訪問~

參考文章：

高斯濾波：

本文主要介紹了高斯濾波器的原理及其實現過程

高斯濾波器是一種線性濾波器，能夠有效的抑制噪聲，平滑影象。其作用原理和均值濾波器類似，都是取濾波器視窗內的畫素的均值作為輸出。其視窗模板的係數和均值濾波器不同，均值濾波器的模板係數都是相同的為1；而高斯濾波器的模板係數，則隨著距離模板中心的增大而係數減小。所以，高斯濾波器相比於均值濾波器對影象個模糊程度較小。

什麼是高斯濾波器

既然名稱為高斯濾波器，那麼其和高斯分佈（正態分佈）是有一定的關係的。一個二維的高斯函式如下：

h(x,y)=e−x2+y22σ2h(x,y)=e−x2+y22σ2

其中(x,y)(x,y)為點座標，在影象處理中可認為是整數；σσ是標準差。要想得到一個高斯濾波器的模板，可以對高斯函式進行離散化，得到的高斯函式值作為模板的係數。例如：要產生一個3×33×3的高斯濾波器模板，以模板的中心位置為座標原點進行取樣。模板在各個位置的座標，如下所示（x軸水平向右，y軸豎直向下）

這樣，將各個位置的座標帶入到高斯函式中，得到的值就是模板的係數。對於視窗模板的大小為 (2k+1)×(2k+1)(2k+1)×(2k+1)，模板中各個元素值的計算公式如下：

Hi,j=12πσ2e−(i−k−1)2+(j−k−1)22σ2Hi,j=12πσ2e−(i−k−1)2+(j−k−1)22σ2

這樣計算出來的模板有兩種形式：小數和整數。

小數形式的模板，就是直接計算得到的值，沒有經過任何的處理；
整數形式的，則需要進行歸一化處理，將模板左上角的值歸一化為1，下面會具體介紹。使用整數的模板時，需要在模板的前面加一個係數，係數為1∑(i,j)∈wwi,j1∑(i,j)∈wwi,j，也就是模板係數和的倒數。

高斯模板的生成

知道模板生成的原理，實現起來也就不困難了

void generateGaussianTemplate(double window[][11], int ksize, double sigma)
{
    static const double pi = 3.1415926;
    int center = ksize / 2; // 模板的中心位置，也就是座標的原點
    double x2, y2;
    for (int i = 0; i < ksize; i++)
    {
        x2 = pow(i - center, 2);
        for (int j = 0; j < ksize; j++)
        {
            y2 = pow(j - center, 2);
            double g = exp(-(x2 + y2) / (2 * sigma * sigma));
            g /= 2 * pi * sigma;
            window[i][j] = g;
        }
    }
    double k = 1 / window[0][0]; // 將左上角的係數歸一化為1
    for (int i = 0; i < ksize; i++)
    {
        for (int j = 0; j < ksize; j++)
        {
            window[i][j] *= k;
        }
    }
}

需要一個二維陣列，存放生成的係數（這裡假設模板的最大尺寸不會超過11）；第二個引數是模板的大小（不要超過11）；第三個引數就比較重要了，是高斯分佈的標準差。生成的過程，首先根據模板的大小，找到模板的中心位置ksize/2。然後就是遍歷，根據高斯分佈的函式，計算模板中每個係數的值。需要注意的是，最後歸一化的過程，使用模板左上角的係數的倒數作為歸一化的係數（左上角的係數值被歸一化為1），模板中的每個係數都乘以該值（左上角係數的倒數），然後將得到的值取整，就得到了整數型的高斯濾波器模板。下面截圖生成的是，大小為3×3,σ=0.83×3,σ=0.8的模板

對上述解結果取整後得到如下模板：

116⎡⎣⎢121242121⎤⎦⎥116[121242121]

這個模板就比較熟悉了，其就是根據σ=0.8σ=0.8的高斯函式生成的模板。

至於小數形式的生成也比較簡單，去掉歸一化的過程，並且在求解過程後，模板的每個係數要除以所有係數的和。具體程式碼如下：

void generateGaussianTemplate(double window[][11], int ksize, double sigma)
{
    static const double pi = 3.1415926;
    int center = ksize / 2; // 模板的中心位置，也就是座標的原點
    double x2, y2;
    double sum = 0;
    for (int i = 0; i < ksize; i++)
    {
        x2 = pow(i - center, 2);
        for (int j = 0; j < ksize; j++)
        {
            y2 = pow(j - center, 2);
            double g = exp(-(x2 + y2) / (2 * sigma * sigma));
            g /= 2 * pi * sigma;
            sum += g;
            window[i][j] = g;
        }
    }
    //double k = 1 / window[0][0]; // 將左上角的係數歸一化為1
    for (int i = 0; i < ksize; i++)
    {
        for (int j = 0; j < ksize; j++)
        {
            window[i][j] /= sum;
        }
    }
}

3×3,σ=0.83×3,σ=0.8的小數型模板。

σσ值的意義及選取

通過上述的實現過程，不難發現，高斯濾波器模板的生成最重要的引數就是高斯分佈的標準差σσ。標準差代表著資料的離散程度，如果σσ較小，那麼生成的模板的中心繫數較大，而周圍的係數較小，這樣對影象的平滑效果就不是很明顯；反之，σσ較大，則生成的模板的各個係數相差就不是很大，比較類似均值模板，對影象的平滑效果比較明顯。

來看下一維高斯分佈的概率分佈密度圖：

橫軸表示可能得取值x，豎軸表示概率分佈密度F(x)，那麼不難理解這樣一個曲線與x軸圍成的圖形面積為1。σσ（標準差）決定了這個圖形的寬度，可以得出這樣的結論：σσ越大，則圖形越寬，尖峰越小，圖形較為平緩；σσ越小，則圖形越窄，越集中，中間部分也就越尖，圖形變化比較劇烈。這其實很好理解，如果sigma也就是標準差越大，則表示該密度分佈一定比較分散，由於面積為1，於是尖峰部分減小，寬度越寬（分佈越分散）；同理，當σσ越小時，說明密度分佈較為集中，於是尖峰越尖，寬度越窄！於是可以得到如下結論： σσ越大，分佈越分散，各部分比重差別不大，於是生成的模板各元素值差別不大，類似於平均模板； σσ越小，分佈越集中，中間部分所佔比重遠遠高於其他部分，反映到高斯模板上就是中心元素值遠遠大於其他元素值，於是自然而然就相當於中間值得點運算。

基於OpenCV的實現

在生成高斯模板好，其簡單的實現和其他的空間濾波器沒有區別，具體程式碼如下：

void GaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)
{
    CV_Assert(src.channels() || src.channels() == 3); // 只處理單通道或者三通道影象
    const static double pi = 3.1415926;
    // 根據視窗大小和sigma生成高斯濾波器模板
    // 申請一個二維陣列，存放生成的高斯模板矩陣
    double **templateMatrix = new double*[ksize];
    for (int i = 0; i < ksize; i++)
        templateMatrix[i] = new double[ksize];
    int origin = ksize / 2; // 以模板的中心為原點
    double x2, y2;
    double sum = 0;
    for (int i = 0; i < ksize; i++)
    {
        x2 = pow(i - origin, 2);
        for (int j = 0; j < ksize; j++)
        {
            y2 = pow(j - origin, 2);
            // 高斯函式前的常數可以不用計算，會在歸一化的過程中給消去
            double g = exp(-(x2 + y2) / (2 * sigma * sigma));
            sum += g;
            templateMatrix[i][j] = g;
        }
    }
    for (int i = 0; i < ksize; i++)
    {
        for (int j = 0; j < ksize; j++)
        {
            templateMatrix[i][j] /= sum;
            cout << templateMatrix[i][j] << " ";
        }
        cout << endl;
    }
    // 將模板應用到影象中
    int border = ksize / 2;
    copyMakeBorder(src, dst, border, border, border, border, BorderTypes::BORDER_REFLECT);
    int channels = dst.channels();
    int rows = dst.rows - border;
    int cols = dst.cols - border;
    for (int i = border; i < rows; i++)
    {
        for (int j = border; j < cols; j++)
        {
            double sum[3] = { 0 };
            for (int a = -border; a <= border; a++)
            {
                for (int b = -border; b <= border; b++)
                {
                    if (channels == 1)
                    {
                        sum[0] += templateMatrix[border + a][border + b] * dst.at<uchar>(i + a, j + b);
                    }
                    else if (channels == 3)
                    {
                        Vec3b rgb = dst.at<Vec3b>(i + a, j + b);
                        auto k = templateMatrix[border + a][border + b];
                        sum[0] += k * rgb[0];
                        sum[1] += k * rgb[1];
                        sum[2] += k * rgb[2];
                    }
                }
            }
            for (int k = 0; k < channels; k++)
            {
                if (sum[k] < 0)
                    sum[k] = 0;
                else if (sum[k] > 255)
                    sum[k] = 255;
            }
            if (channels == 1)
                dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
            else if (channels == 3)
            {
                Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
                dst.at<Vec3b>(i, j) = rgb;
            }
        }
    }
    // 釋放模板陣列
    for (int i = 0; i < ksize; i++)
        delete[] templateMatrix[i];
    delete[] templateMatrix;
}

只處理單通道或者三通道影象，模板生成後，其濾波（卷積過程）就比較簡單了。不過，這樣的高斯濾波過程，其迴圈運算次數為m×n×ksize2m×n×ksize2，其中m，n為影象的尺寸；ksize為高斯濾波器的尺寸。這樣其時間複雜度為O(ksize2)O(ksize2)，隨濾波器的模板的尺寸呈平方增長，當高斯濾波器的尺寸較大時，其運算效率是極低的。為了，提高濾波的運算速度，可以將二維的高斯濾波過程分解開來。

分離實現高斯濾波

由於高斯函式的可分離性，尺寸較大的高斯濾波器可以分成兩步進行：首先將影象在水平（豎直）方向與一維高斯函式進行卷積；然後將卷積後的結果在豎直（水平）方向使用相同的一維高斯函式得到的模板進行卷積運算。具體實現程式碼如下：

// 分離的計算
void separateGaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)
{
    CV_Assert(src.channels()==1 || src.channels() == 3); // 只處理單通道或者三通道影象
    // 生成一維的高斯濾波模板
    double *matrix = new double[ksize];
    double sum = 0;
    int origin = ksize / 2;
    for (int i = 0; i < ksize; i++)
    {
        // 高斯函式前的常數可以不用計算，會在歸一化的過程中給消去
        double g = exp(-(i - origin) * (i - origin) / (2 * sigma * sigma));
        sum += g;
        matrix[i] = g;
    }
    // 歸一化
    for (int i = 0; i < ksize; i++)
        matrix[i] /= sum;
    // 將模板應用到影象中
    int border = ksize / 2;
    copyMakeBorder(src, dst, border, border, border, border, BorderTypes::BORDER_REFLECT);
    int channels = dst.channels();
    int rows = dst.rows - border;
    int cols = dst.cols - border;
    // 水平方向
    for (int i = border; i < rows; i++)
    {
        for (int j = border; j < cols; j++)
        {
            double sum[3] = { 0 };
            for (int k = -border; k <= border; k++)
            {
                if (channels == 1)
                {
                    sum[0] += matrix[border + k] * dst.at<uchar>(i, j + k); // 行不變，列變化；先做水平方向的卷積
                }
                else if (channels == 3)
                {
                    Vec3b rgb = dst.at<Vec3b>(i, j + k);
                    sum[0] += matrix[border + k] * rgb[0];
                    sum[1] += matrix[border + k] * rgb[1];
                    sum[2] += matrix[border + k] * rgb[2];
                }
            }
            for (int k = 0; k < channels; k++)
            {
                if (sum[k] < 0)
                    sum[k] = 0;
                else if (sum[k] > 255)
                    sum[k] = 255;
            }
            if (channels == 1)
                dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
            else if (channels == 3)
            {
                Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
                dst.at<Vec3b>(i, j) = rgb;
            }
        }
    }
    // 豎直方向
    for (int i = border; i < rows; i++)
    {
        for (int j = border; j < cols; j++)
        {
            double sum[3] = { 0 };
            for (int k = -border; k <= border; k++)
            {
                if (channels == 1)
                {
                    sum[0] += matrix[border + k] * dst.at<uchar>(i + k, j); // 列不變，行變化；豎直方向的卷積
                }
                else if (channels == 3)
                {
                    Vec3b rgb = dst.at<Vec3b>(i + k, j);
                    sum[0] += matrix[border + k] * rgb[0];
                    sum[1] += matrix[border + k] * rgb[1];
                    sum[2] += matrix[border + k] * rgb[2];
                }
            }
            for (int k = 0; k < channels; k++)
            {
                if (sum[k] < 0)
                    sum[k] = 0;
                else if (sum[k] > 255)
                    sum[k] = 255;
            }
            if (channels == 1)
                dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
            else if (channels == 3)
            {
                Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
                dst.at<Vec3b>(i, j) = rgb;
            }
        }
    }
    delete[] matrix;
}

程式碼沒有重構較長，不過其實現原理是比較簡單的。首先得到一維高斯函式的模板，在卷積（濾波）的過程中，保持行不變，列變化，在水平方向上做卷積運算；接著在上述得到的結果上，保持列不邊，行變化，在豎直方向上做卷積運算。這樣分解開來，演算法的時間複雜度為O(ksize)O(ksize)，運算量和濾波器的模板尺寸呈線性增長。

在OpenCV也有對高斯濾波器的封裝GaussianBlur，其宣告如下：

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

二維高斯函式的標準差在x和y方向上應該分別有一個標準差，在上面的程式碼中一直設其在x和y方向的標準是相等的，在OpenCV中的高斯濾波器中，可以在x和y方向上設定不同的標準差。下圖是自己實現的高斯濾波器和OpenCV中的GaussianBlur的結果對比

上圖是5×5,σ=0.85×5,σ=0.8的高斯濾波器，可以看出兩個實現得到的結果沒有很大的區別。

總結

高斯濾波器是一種線性平滑濾波器，其濾波器的模板是對二維高斯函式離散得到。由於高斯模板的中心值最大，四周逐漸減小，其濾波後的結果相對於均值濾波器來說更好。高斯濾波器最重要的引數就是高斯分佈的標準差σσ，標準差和高斯濾波器的平滑能力有很大的能力，σσ越大，高斯濾波器的頻帶就較寬，對影象的平滑程度就越好。通過調節σσ引數，可以平衡對影象的噪聲的抑制和對影象的模糊。

高斯濾波在影象處理中的應用

卷積：

什麼是卷積

卷積的核心（涉及推導過程，可以跳過）：

卷積的應用

參考文章：

高斯濾波：

什麼是高斯濾波器

高斯模板的生成

σσ值的意義及選取

基於OpenCV的實現

分離實現高斯濾波

總結

高斯濾波——影象處理

高斯濾波在影象處理中的應用

基於MATLAB影象處理的中值濾波、均值濾波以及高斯濾波的實現與對比

[Python影象處理] 四.影象平滑之均值濾波、方框濾波、高斯濾波及中值濾波

VC++高斯濾波\中值濾波實現影象模糊處理

OpenCV計算機視覺學習（4）——影象平滑處理（均值濾波，高斯濾波，中值濾波，雙邊濾波）

影象處理（三）——高斯濾波

opencv影象處理-------高斯濾波

opencv實現影象鄰域均值濾波、中值濾波、高斯濾波

影象識別與處理之Opencv——高斯濾波GaussianBlur() 11月2日暫存

基於opencv下對視頻的灰度變換，高斯濾波，canny邊緣檢測處理，同窗體顯示並保存

高斯濾波邊界優化處理

opencv3 圖片模糊操作-均值濾波高斯濾波中值濾波雙邊濾波

模糊集理論在影象處理中的應用

數字影象處理——中值濾波

方框濾波/均值濾波/高斯濾波/中值濾波/雙邊濾波

數字影象處理中的頻率域濾波

【機器學習】EM演算法在高斯混合模型學習中的應用

高斯濾波對影象方差有什麼影響

數字影象處理中的空間濾波

高斯濾波在影象處理中的應用

卷積：

什麼是卷積

卷積的核心（涉及推導過程，可以跳過）：

卷積的應用

參考文章：

高斯濾波：

什麼是高斯濾波器

高斯模板的生成

σσ值的意義及選取

基於OpenCV的實現

分離實現高斯濾波

總結

相關推薦