高斯滤波在图像处理中的应用

た入场券 2022-05-12 04:34 202阅读 0赞

# **卷积：** #

相信很多时候，当我们在看到“卷积”时，总是处于一脸懵逼的状态，不但因为它的本义概念比较难理解，还因为它在不同的应用中发挥出的变幻莫测的作用也时常让人迷糊。但这些应用其实本质上都是同一种东西，理解了卷积的来源，就可以举一反三。其实我个人对于卷积的理解，很长时间都处于似懂非懂的状态，就像傅里叶变换的一些tricky points，只求在应用中不出差错，不求甚解。但是如果想要真正做好learning的东西，我认为在真正的学习理论基础前，必须将概念的本质搞清楚。因此，这篇文章便是为工科，理科以及其他领域对此感兴趣的同学而整理，内容简易通俗，便于理解。

### 什么是卷积 ###

卷积，很多时候都是我们在各种工程领域，信号领域所看到的常用名词，比如系统

![v2-39ea42bd1c601b571ded147b881c3dc6_hd.jpg][]

通俗易懂的说，就是

**输出 = 输入 \* 系统**

虽然它看起来只是个复杂的数学公式，但是却有着重要的物理意义，因为自然界这样的系统无处不在，计算一个系统的输出最好的方法就是运用卷积。更一般的，我们还有很多其他领域的应用：

**统计学**中，加权的滑动平均是一种卷积。

**概率论**中，两个统计独立变量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][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\}][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][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][x_n]是“基”，![a\_n][a_n]是在这个基上的展开系数。两个多项式乘积的在基上展开的系数就是两个多项式各自在基上展开系数的卷积。  
![x^n][x_n]对应着频率不同的![\\exp(ikt)][exp_ikt]，系数对应着其傅里叶变换。自然就是乘积的傅里叶变换等于傅里叶变换的卷积了。

### 卷积的内核（涉及推导过程，可以跳过）： ###

首先我们有这样一个概念：**内积、积分、投影这三者其实从某个角度上讲是一个意思**。

定义一组向量![\\alpha=(\\alpha\_\{1\},\\alpha\_\{2\},...,\\alpha\_\{n\})][alpha_alpha_1_alpha_2_..._alpha_n]，另一组向量![\\beta=(\\beta\_\{1\},\\beta\_\{2\},...,\\beta\_\{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\}\}][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\})][alpha_1_alpha_2_..._alpha_n]，这和一个点在3D欧几里得空间的三轴投影坐标是一个道理。

这样，先来看看Fourier变换在做什么：

![F(\\omega)=\\int\_\{-\\infty \}^\{+\\infty\} f(t)e^\{-j\\omega t\}dt][F_omega_int_-_infty _infty_ f_t_e_-j_omega t_dt]

再引入一个完美的式子，欧拉公式：

![e^\{ix\}=cosx+isinx][e_ix_cosx_isinx]

从Fourier的定义式可以看出是对f(t)和![e^\{-j\\omega t\}][e_-j_omega t]相乘后在无穷域上对其进行积分，那么其实就是将f(t)投影在![e^\{-j\\omega t\}][e_-j_omega 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][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\}][e_-_sigma t]将其进行衰减后再做Fourier变换。到了负无穷的区域这衰减因子可就成了递增因子，所以Laplace变换仅限于大于0的区域，对于小于0的区域用系统初始状态表达就好了。从这点角度上讲，Laplace变换相当于对![f(t)e^\{-\\sigma t\}][f_t_e_-_sigma t]做了一个单边Fourier变换。

然后，分析方法同上，可以看到Laplace把f(t)投影到了![e^\{-st\}][e_-st]空间，这就是s平面。它比Fourier更厉害的地方是不仅可以看到虚轴上![\\omega][omega]的成分，还可以在实轴上看到Fourier看不到的衰减因子![\\sigma][sigma]成分，这是Fourier做不到的。所以Laplace在Fourier的基础上把信号拓展到了衰减因子实轴上，这个衰减因子![\\sigma][sigma]和系统的阻尼![\\zeta][zeta]，自然震荡角频率![\\omega\_\{n\}][omega_n]密切相关，直接影响了系统的调节时间![t\_\{s\}][t_s]。学过自控原理的同学应该知道在频域章节，我们得到系统的频域响应曲线都是通过传递函数来直接转化的，公式就是s=jw。这也就是说Fourier活在一维虚轴空间，Laplace活在二维平面空间，想要得到一维空间上关于w的表达形式，只需要在s平面上做降维处理即可。

回过头来再来看看卷积投影：

![f(t)\*g(t)=\\int\_\{-\\infty \}^\{+\\infty\}f(\\tau)g(t-\\tau)d\\tau][f_t_g_t_int_-_infty _infty_f_tau_g_t-_tau_d_tau]

这个投影有点奇怪，它在投影之前先把![g(\\tau)][g_tau]做了一个反对称，然后再投影。对应到前面推导的系统卷积表达式：

![r(t)\*f(t)=\\int\_\{-\\infty \}^\{+\\infty\}f(\\tau)r(t-\\tau)d\\tau][r_t_f_t_int_-_infty _infty_f_tau_r_t-_tau_d_tau]

相当于在投影之前，先把输入信号r(t)在时间轴上翻转了180°，然后与系统f(t)进行投影。投影的概念我们可以很好理解，无论是向量内积运算相当于线投影，或者空间的一个多面体在三维空间平面上的投影面，这种投影运算就相当于一种重合面积。如果从这个角度去看输入、系统和输出三者之间的关系，那么就可以从图形角度去理解为什么一个一阶系统在阶跃响应输出下是一条单调上升的曲线了。这里用一张wikipedia里关于卷积的一张图形化解释，想要了解更多的同学可以自行跳转：[Convolution][]

（特此感谢知乎学霸王尼莫的帮助）

![0df97130e37577d4e176620ab432a7fe_hd.jpg][]

### 卷积的应用 ###

图像处理：用一个模板和一幅图像进行卷积，对于图像上的一个点，让模板的原点和该点重合，然后模板上的点和图像上对应的点相乘，然后各点的积相加，就得到该点的卷积值。对图像上的每个点都这样处理。由于多数模板都对称，所以模板不旋转。  
卷积是一种积分运算，用来求两个曲线重叠区域面积。可以看作加权求和，可以用来消除噪声、特征增强。 把一个点的像素值用它周围的点的像素值的加权平均代替。  
  
卷积是一种线性运算,图像处理中常见的mask运算都是卷积，广泛应用于图像滤波。   
卷积关系最重要的一种情况，就是在信号与线性系统或数字信号处理中的卷积定理。利用该定理，可以将时间域或空间域中的卷积运算等价为频率域的相乘运算，从而利用FFT等快速算法，实现有效的计算，节省运算代价。

下面是来自sselssbh博客的一个例子，非常形象的解释了卷积在图像领域的作用

有这么一副图像，可以看到，图像上有很多噪点：   
![这里写图片描述][SouthEast]

高频信号，就好像平地耸立的山峰：   
![这里写图片描述][SouthEast 1]

看起来很显眼。

平滑这座山峰的办法之一就是，把山峰刨掉一些土，填到山峰周围去。用数学的话来说，就是把山峰周围的高度平均一下。

平滑后得到：   
![这里写图片描述][SouthEast 2]

4.2 计算

卷积可以帮助实现这个平滑算法。

有噪点的原图，可以把它转为一个矩阵：   
![这里写图片描述][SouthEast 3]

然后用下面这个平均矩阵（说明下，原图的处理实际上用的是正态分布矩阵，这里为了简单，就用了算术平均矩阵）来平滑图像：

记得刚才说过的算法，把高频信号与周围的数值平均一下就可以平滑山峰。

比如我要平滑 点，就在矩阵中，取出点附近的点组成矩阵 f ,和 g 进行卷积计算后，再填回去   
![这里写图片描述][SouthEast 4]

要注意一点，为了运用卷积， g 虽然和 f 同维度，但下标有点不一样：   
![这里写图片描述][SouthEast 5]

![这里写图片描述][SouthEast 6]

写成卷积公式就是：

要求，一样可以套用上面的卷积公式。

这样相当于实现了 g 这个矩阵在原来图像上的划动（准确来说，下面这幅图把 g 矩阵旋转了 ）：

再比如做馒头

楼下早点铺子生意太好了，供不应求，就买了一台机器，不断的生产馒头。

假设馒头的生产速度是 f(t) ，那么一天后生产出来的馒头总量为： ∫240f(t)dt

馒头生产出来之后，就会慢慢腐败，假设腐败函数为 g(t) ，比如，10个馒头，24小时会腐败： 10∗g(t)

想想就知道，第一个小时生产出来的馒头，一天后会经历24小时的腐败，第二个小时生产出来的馒头，一天后会经历23小时的腐败。 如此，我们可以知道，一天后，馒头总共腐败了： ∫240f(t)g(24−t)dt   
  
文章也发表在我的个人博客中：[点击打开链接][Link 1] ，更多与机器学习，数学相关知乎，欢迎访问~

### 参考文章： ###

1. CSDN：最容易理解的对卷积(convolution)的解释[http://blog.csdn.net/bitcarmanlee/article/details/54729807][http_blog.csdn.net_bitcarmanlee_article_details_54729807]

2. 知乎：如何通俗易懂的解释卷积？[https://www.zhihu.com/question/22298352][http_blog.csdn.net_bitcarmanlee_article_details_54729807]

3. 卷积为什么叫卷子？[https://www.zhihu.com/question/54677157][https_www.zhihu.com_question_54677157]

4. 果壳：关于卷积的一个血腥讲解 [https://www.guokr.com/post/342476/][https_www.guokr.com_post_342476]

5. 小木虫：[https://zh.wikipedia.org/wiki/%E5%8D%B7%E7%A7%AF%E5%AE%9A%E7%90%86 ][https_zh.wikipedia.org_wiki_E5_8D_B7_E7_A7_AF_E5_AE_9A_E7_90_86]

原文转自:[https://blog.csdn.net/qq\_39521554/article/details/79083864][https_blog.csdn.net_qq_39521554_article_details_79083864]

## 高斯滤波： ##

本文主要介绍了高斯滤波器的原理及其实现过程

高斯滤波器是一种线性滤波器，能够有效的抑制噪声，平滑图像。其作用原理和均值滤波器类似，都是取滤波器窗口内的像素的均值作为输出。其窗口模板的系数和均值滤波器不同，均值滤波器的模板系数都是相同的为1；而高斯滤波器的模板系数，则随着距离模板中心的增大而系数减小。所以，高斯滤波器相比于均值滤波器对图像个模糊程度较小。

## 什么是高斯滤波器 ##

既然名称为高斯滤波器，那么其和高斯分布（正态分布）是有一定的关系的。一个二维的高斯函数如下：

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

其中(x,y)(x,y)为点坐标，在图像处理中可认为是整数；σσ是标准差。要想得到一个高斯滤波器的模板，可以对高斯函数进行离散化，得到的高斯函数值作为模板的系数。例如：要产生一个3×33×3的高斯滤波器模板，以模板的中心位置为坐标原点进行取样。模板在各个位置的坐标，如下所示（x轴水平向右，y轴竖直向下）  
![439761-20170216214155582-1833428662.png][]

这样，将各个位置的坐标带入到高斯函数中，得到的值就是模板的系数。  
对于窗口模板的大小为 (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的模板  
![439761-20170216214219035-2030191739.png][]

对上述解结果取整后得到如下模板：

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的小数型模板。

![439761-20170216214250660-931432778.png][]

## σσ值的意义及选取 ##

通过上述的实现过程，不难发现，高斯滤波器模板的生成最重要的参数就是高斯分布的标准差σσ。标准差代表着数据的离散程度，如果σσ较小，那么生成的模板的中心系数较大，而周围的系数较小，这样对图像的平滑效果就不是很明显；反之，σσ较大，则生成的模板的各个系数相差就不是很大，比较类似均值模板，对图像的平滑效果比较明显。

来看下一维高斯分布的概率分布密度图：  
![439761-20170216214529738-2113566517.png][]

横轴表示可能得取值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`的结果对比  
![439761-20170216214326675-853983120.png][]

上图是5×5,σ=0.85×5,σ=0.8的高斯滤波器，可以看出两个实现得到的结果没有很大的区别。

## 总结 ##

高斯滤波器是一种线性平滑滤波器，其滤波器的模板是对二维高斯函数离散得到。由于高斯模板的中心值最大，四周逐渐减小，其滤波后的结果相对于均值滤波器来说更好。  
高斯滤波器最重要的参数就是高斯分布的标准差σσ，标准差和高斯滤波器的平滑能力有很大的能力，σσ越大，高斯滤波器的频带就较宽，对图像的平滑程度就越好。通过调节σσ参数，可以平衡对图像的噪声的抑制和对图像的模糊。

原文转自：[https://www.cnblogs.com/wangguchangqing/p/6407717.html][https_www.cnblogs.com_wangguchangqing_p_6407717.html]

[v2-39ea42bd1c601b571ded147b881c3dc6_hd.jpg]: /images/20220512/444828cbcbbc4a4096fb77052ddad2f3.png
[V_t_ _ _int_x_y_t_ U_x_y_ _mathrm_d_x]: https://www.zhihu.com/equation?tex=V%28t%29+%3D+%5Cint_%7Bx%2By%3Dt%7D+U%28x%2Cy%29+%5C%2C%5Cmathrm%7Bd%7Dx
[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]: https://www.zhihu.com/equation?tex=%5Cbegin%7Balign%7D%0A42+%5Ctimes137+%26%3D+%282%5Ctimes10%5E0%2B4%5Ctimes10%5E1%29%287%5Ctimes10%5E0%2B3%5Ctimes10%5E1%2B1%5Ctimes10%5E2%29+%5C%5C%0A%26%3D+%282%5Ctimes7%29%5Ctimes10%5E0+%2B+%282%5Ctimes3%2B4%5Ctimes7%29%5Ctimes10%5E1%2B%282%5Ctimes1%2B4%5Ctimes3%29%5Ctimes10%5E2+%2B+%284%5Ctimes1%29%5Ctimes10%5E3+%5C%5C%0A%26%3D+14+%2B+340%2B1400%2B4000+%5C%5C%0A%26%3D+5754%0A%5Cend%7Balign%7D
[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]: /images/20220512/8ab2cbbe5e154bcaaac9303cdee8884b.png
[x_n]: https://www.zhihu.com/equation?tex=x%5En
[a_n]: /images/20220512/f6c07d12ecc74588a4e833dd216c0bad.png
[exp_ikt]: /images/20220512/ce3aeeafff1547e896fd8a0d75f014fe.png
[alpha_alpha_1_alpha_2_..._alpha_n]: https://www.zhihu.com/equation?tex=%5Calpha%3D%28%5Calpha_%7B1%7D%2C%5Calpha_%7B2%7D%2C...%2C%5Calpha_%7Bn%7D%29
[beta_beta_1_beta_2_..._beta_n]: /images/20220512/9abad7d797d042ba9181c25499602b8e.png
[alpha _beta _ _alpha_1_ _beta_1_alpha_2_ _beta_2_..._alpha_n_ _beta_n_sum_i_1_n_alpha_i_ _beta_i]: https://www.zhihu.com/equation?tex=%5Calpha+%5Cbeta+%3D+%5Calpha_%7B1%7D+%5Cbeta_%7B1%7D%2B%5Calpha_%7B2%7D+%5Cbeta_%7B2%7D%2B...%2B%5Calpha_%7Bn%7D+%5Cbeta_%7Bn%7D%3D%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7B%5Calpha_%7Bi%7D+%5Cbeta_%7Bi%7D%7D+
[alpha_1_alpha_2_..._alpha_n]: /images/20220512/90c188909e1340f4a9ecab478e9821a5.png
[F_omega_int_-_infty _infty_ f_t_e_-j_omega t_dt]: https://www.zhihu.com/equation?tex=F%28%5Comega%29%3D%5Cint_%7B-%5Cinfty+%7D%5E%7B%2B%5Cinfty%7D+f%28t%29e%5E%7B-j%5Comega+t%7Ddt+
[e_ix_cosx_isinx]: /images/20220512/a0779d697f724b60a2acf4bb9355c2f1.png+
[e_-j_omega t]: https://www.zhihu.com/equation?tex=e%5E%7B-j%5Comega+t%7D+
[F_s_int_0_infty_f_t_e_-_sigma t_e_-j_omega t_dt_int_0_infty_f_t_e_-st_dt]: https://www.zhihu.com/equation?tex=F%28s%29%3D%5Cint_%7B0%7D%5E%7B%5Cinfty%7Df%28t%29e%5E%7B-%5Csigma+t%7De%5E%7B-j%5Comega+t%7Ddt%3D%5Cint_%7B0%7D%5E%7B%5Cinfty%7Df%28t%29e%5E%7B-st%7Ddt
[e_-_sigma t]: https://www.zhihu.com/equation?tex=e%5E%7B-%5Csigma+t%7D+
[f_t_e_-_sigma t]: https://www.zhihu.com/equation?tex=f%28t%29e%5E%7B-%5Csigma+t%7D+
[e_-st]: https://www.zhihu.com/equation?tex=e%5E%7B-st%7D+
[omega]: https://www.zhihu.com/equation?tex=%5Comega+
[sigma]: https://www.zhihu.com/equation?tex=%5Csigma+
[zeta]: https://www.zhihu.com/equation?tex=%5Czeta+
[omega_n]: https://www.zhihu.com/equation?tex=%5Comega_%7Bn%7D+
[t_s]: /images/20220512/5d67e9bb9e4646c794642f0daf674059.png+
[f_t_g_t_int_-_infty _infty_f_tau_g_t-_tau_d_tau]: https://www.zhihu.com/equation?tex=f%28t%29%2Ag%28t%29%3D%5Cint_%7B-%5Cinfty+%7D%5E%7B%2B%5Cinfty%7Df%28%5Ctau%29g%28t-%5Ctau%29d%5Ctau
[g_tau]: /images/20220512/5d83fe01ed0b41f78d3ff549ba462908.png
[r_t_f_t_int_-_infty _infty_f_tau_r_t-_tau_d_tau]: https://www.zhihu.com/equation?tex=r%28t%29%2Af%28t%29%3D%5Cint_%7B-%5Cinfty+%7D%5E%7B%2B%5Cinfty%7Df%28%5Ctau%29r%28t-%5Ctau%29d%5Ctau
[Convolution]: https://link.zhihu.com/?target=https%3A//en.wikipedia.org/wiki/Convolution
[0df97130e37577d4e176620ab432a7fe_hd.jpg]: /images/20220512/00d8f770a0674763a4e016bfba56640a.png
[SouthEast]: /images/20220512/c00b8f3d591d4e57af483ffcf5d8af14.png
[SouthEast 1]: /images/20220512/9432b418138c45338e8429fb3875a4cb.png
[SouthEast 2]: /images/20220512/f1750837e4864a4f97dc54618f176dcf.png
[SouthEast 3]: /images/20220512/9e921ac58ed44a7bb82bd3561f17ce0f.png
[SouthEast 4]: /images/20220512/a87296a08a15451aade4c0271a7a976e.png
[SouthEast 5]: /images/20220512/b699d8ae79ed4a60a7a772e093cb0b13.png
[SouthEast 6]: /images/20220512/4fe471279fde4132aab7913584a1f23e.png
[Link 1]: http://www.yiiran.wang/index.php/2018/01/17/%E5%88%AB%E6%80%95%EF%BC%8C%E5%8D%B7%E7%A7%AF%E5%85%B6%E5%AE%9E%E5%BE%88%E7%AE%80%E5%8D%95/
[http_blog.csdn.net_bitcarmanlee_article_details_54729807]: http://blog.csdn.net/bitcarmanlee/article/details/54729807
[https_www.zhihu.com_question_54677157]: https://www.zhihu.com/question/54677157
[https_www.guokr.com_post_342476]: https://www.guokr.com/post/342476/
[https_zh.wikipedia.org_wiki_E5_8D_B7_E7_A7_AF_E5_AE_9A_E7_90_86]: https://zh.wikipedia.org/wiki/%E5%8D%B7%E7%A7%AF%E5%AE%9A%E7%90%86
[https_blog.csdn.net_qq_39521554_article_details_79083864]: https://blog.csdn.net/qq_39521554/article/details/79083864
[439761-20170216214155582-1833428662.png]: /images/20220512/b7270fa9abfe4bdf95c19655bf5f1fbb.png
[439761-20170216214219035-2030191739.png]: /images/20220512/77a24dc8bffe4ce186a2e95206623641.png
[439761-20170216214250660-931432778.png]: /images/20220512/5f3ce4d2555340d5af7a452bc9cbc5f5.png
[439761-20170216214529738-2113566517.png]: /images/20220512/62165ad932e14fc687cb33540d725487.png
[439761-20170216214326675-853983120.png]: /images/20220512/b5d6dc9a85994fae91823ab76b663194.png
[https_www.cnblogs.com_wangguchangqing_p_6407717.html]: https://www.cnblogs.com/wangguchangqing/p/6407717.html