高斯滤波

转载地址：http://blog.csdn.NET/sunmc1204953974/article/details/50634652

高斯滤波

图像滤波之高斯滤波（Gauss filter）

概述：

高斯滤波：

高斯滤波在图像处理概念下，将图像频域处理和时域处理相联系，作为低通滤波器使用，可以将低频能量（比如噪声）滤去，起到图像平滑作用。

高斯滤波是一种线性平滑滤波，适用于消除高斯噪声，广泛应用于图像处理的减噪过程。通俗的讲，高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到。高斯滤波的具体操作是：用一个模板（或称卷积、掩模）扫描图像中的每一个像素，用模板确定的邻域内像素的加权平均灰度值去替代模板中心像素点的值用。高斯平滑滤波器对于抑制服从正态分布的噪声非常有效。

高斯模糊：

我们常说的高斯模糊就是使用高斯滤波器完成的，高斯模糊是低通滤波的一种，也就是滤波函数是低通高斯函数，但是高斯滤波是指用高斯函数作为滤波函数，至于是不是模糊，要看是高斯低通还是高斯高通，低通就是模糊，高通就是锐化。

算法步骤：

在图像处理中，高斯滤波一般有两种实现方式，一是用离散化窗口滑窗卷积，另一种通过傅里叶变换。最常见的就是第一种滑窗实现，只有当离散化的窗口非常大，用滑窗计算量非常大（即使用可分离滤波器的实现）的情况下，可能会考虑基于傅里叶变化的实现方法。

由于高斯函数可以写成可分离的形式，因此可以采用可分离滤波器实现来加速。所谓的可分离滤波器，就是可以把多维的卷积化成多个一维卷积。具体到二维的高斯滤波，就是指先对行做一维卷积，再对列做一维卷积。这样就可以将计算复杂度从O(M*M*N*N)降到O(2*M*M*N)，M，N分别是图像和滤波器的窗口大小。

高斯模糊是一个非常典型的图像卷积例子，本质上，高斯模糊就是将(灰度)图像 I 和一个高斯核进行卷积操作:

其中 * 表示卷积操作; Gσ 是标准差为σ 的二维高斯核,定义为:

这里补充以下卷积的知识： 
卷积是分析数学中一种重要的运算。 
设:f(x),g(x)是R1上的两个可积函数，作积分：

可以证明，关于几乎所有的实数x，上述积分是存在的。这样，随着x的不同取值，这个积分就定义了一个新函数h(x)，称为函数f与g的卷积，记为h(x)=(f*g)(x)。

卷积是一个单纯的定义，本身没有什么意义可言，但是其在各个领域的应用是十分广泛的，在滤波中可以理解为一个加权平均过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到,而如何加权则是依据核函数高斯函数。

平均的过程:

对于图像来说，进行平滑和模糊，就是利用周边像素的平均值

“中间点”取”周围点”的平均值，就会变成1。在数值上，这是一种”平滑化”。在图形上，就相当于产生”模糊”效果，”中间点”失去细节。

显然，计算平均值时，取值范围越大，”模糊效果”越强烈。

上面分别是原图、模糊半径3像素、模糊半径10像素的效果。模糊半径越大，图像就越模糊。从数值角度看，就是数值越平滑。

如果使用简单平均，显然不是很合理，因为图像都是连续的，越靠近的点关系越密切，越远离的点关系越疏远。因此，加权平均更合理，距离越近的点权重越大，距离越远的点权重越小。

而正态分布显然是一种可取的权重分配模式。由于图像是二维的，所以需要使用二维的高斯函数：

二维高斯函数图像:

一维高斯函数:

中心点就是原点，μ等于0:

二维高斯函数（中心为原点）:

计算平均值的时候，我们只需要将”中心点”作为原点，其他点按照其在正态曲线上的位置，分配权重，就可以得到一个加权平均值，而这就是上述的与二维高斯核进行卷积的过程。

计算权重矩阵的过程:

假定中心点的坐标是（0,0），那么距离它最近的8个点的坐标如下：

假定σ=1.5，则模糊半径为1的权重矩阵如下：

这9个点的权重总和等于0.4787147，如果只计算这9个点的加权平均，还必须让它们的权重之和等于1，因此上面9个值还要分别除以0.4787147，得到最终的权重矩阵：

有了权重矩阵，就可以计算高斯模糊的值了。 
中心点以及周边n个点，每个点乘以自己的权重值并将这些值相加，就是中心点的高斯模糊的值。对所有点重复这个过程，就得到了高斯模糊后的图像。

如果原图是彩色图片，可以对RGB三个通道分别做高斯模糊。

对于边界点来说，周边没有足够的点，一个变通方法是把已有的点拷贝到另一面的对应位置，模拟出完整的矩阵

第一个问题：高斯函数为什么能作为图像处理中的滤波函数？

高斯平滑滤波器无论在空间域还是在频率域都是十分有效的低通滤波器，且在实际图像处理中得到了工程人员的有效使用．高斯函数具有五个十分重要的性质，它们是：

（1）二维高斯函数具有旋转对称性，即滤波器在各个方向上的平滑程度是相同的．一般来说，一幅图像的边缘方向是事先不知道的，因此，在滤波前是无法确定一个方向上比另一方向上需要更多的平滑．旋转对称性意味着高斯平滑滤波器在后续边缘检测中不会偏向任一方向．

（2）高斯函数是单值函数．这表明，高斯滤波器用像素邻域的加权均值来代替该点的像素值，而每一邻域像素点权值是随该点与中心点的距离单调增减的．这一性质是很重要的，因为边缘是一种图像局部特征，如果平滑运算对离算子中心很远的像素点仍然有很大作用，则平滑运算会使图像失真．

（3）高斯函数的付立叶变换频谱是单瓣的．正如下面所示，这一性质是高斯函数付立叶变换等于高斯函数本身这一事实的直接推论．图像常被不希望的高频信号所污染(噪声和细纹理)．而所希望的图像特征（如边缘），既含有低频分量，又含有高频分量．高斯函数付立叶变换的单瓣意味着平滑图像不会被不需要的高频信号所污染，同时保留了大部分所需信号．

（4）高斯滤波器宽度(决定着平滑程度)是由参数σ表征的，而且σ和平滑程度的关系是非常简单的．σ越大，高斯滤波器的频带就越宽，平滑程度就越好．通过调节平滑程度参数σ，可在图像特征过分模糊(过平滑)与平滑图像中由于噪声和细纹理所引起的过多的不希望突变量(欠平滑)之间取得折衷．

（5）由于高斯函数的可分离性，大高斯滤波器可以得以有效地实现．二维高斯函数卷积可以分两步来进行，首先将图像与一维高斯函数进行卷积，然后将卷积结果与方向垂直的相同一维高斯函数卷积．因此，二维高斯滤波的计算量随滤波模板宽度成线性增长而不是成平方增长．

第二个问题：如何计算高斯函数模板（高斯核）？

其实，只要知道模板的大小和高斯函数的方差sigma,由二维高斯函数的表达式很容易计算出高斯核，只要在归一化就可以了。但是由高斯函数的分布特性可知落在u-3*sigma到u+3*sigma的概率大于百分之九九，所以模板大小的选取往往与sigma的取值是相关的，一般而言：

取dim = 1 + 2 * ((int) (3.0 * sigma));OpenCV和sift中的源码也是这么做的，当然实际中可以其实没有这么严格。我们可以使用matlab中的函数直接计算出高斯核，例如3x3的高斯模板：filter=fspecial('gaussian',3,1);其中sigma=1;

sigma的取值决定了高斯函数窗口的大小。在实际中经常看到sigma取值0.8或者1。正常情况下我们由高斯函数计算得到的模板是浮点型数，即double，但是有些情况我们为了加快计算需要将模板处理成整数，对于常见的3x3或者5x5其整数模板如下：

对于浮点数的模板其实往往也都比较固定。很多资料上都有。

第三个问题：可分离的计算

实际上，模板运算（滑动窗口卷积）在数字图像处理中是一项非常耗时的运算。

以上图中的3*3高斯模板为例，每个像素完成一次模板操作要用9个乘法、8个加法和1个除法。对于一幅n*n的图像，大约就是9n²个乘法，8n²个加法和n²个除法，这对于比较大的图像来说，是非常可怕的。而且随着模板大小的增加，计算量是呈指数增长的。那么有没有一种办法能够减少计算量呢？答案是肯定的。由于高斯函数可以写成可分离的形式，因此可以采用可分离滤波器实现来加速。所谓的可分离滤波器，就是可以把一个多维的卷积化成多个一维的卷积。具体到二维的高斯滤波，就是指先对行做一维卷积，再对列做一维卷积。这样就可以将计算复杂度从O(M*M*N*N)降到O(2*M*M*N)，M, N分别是图像和滤波器的窗口大小。问题是实现的时候怎么计算一维的卷积核呢？

其实很简单，按照前面计算出来的窗口大小，将二维的高斯模板合并成一维，计算所有离散点上一维高斯函数的权值，最后将权值之和归一化到1。

编程实例：

void CreateGaussKernel(double sigma, int size,double *pKernel) {double dSum=0,dValue,dDis;for (int i = 0; i < size; i++) {dDis = i - size / 2;dValue = exp(-dDis*dDis / (sigma*sigma) / 2) / sqrt(2 * 3.1415926) / sigma;pKernel[i] = dValue;dSum += dValue;}for (int i = 0; i < size; i++)pKernel[i] /= dSum;}void GaussianBlu(InputArray _src, OutputArray _dst, int size,double sigma) {Mat src = _src.getMat();_dst.create(src.size(), src.type());//imshow("123", src);Mat dst = _dst.getMat();src.copyTo(dst);if (size == 1){src.copyTo(dst);return;}//奇数核size = (size <= 0) ? 1 + 2 *(int) ceil(3 * sigma) : size;CV_Assert(size % 2 == 1);double *pKernel = new double[size];CreateGaussKernel(sigma, size, pKernel);CV_Assert(src.isContinuous());int nLen = size / 2,rows=src.rows,cols=src.cols;double dMul1, dMul2, dMul3,dWeight1=0, dWeight2=0, dWeight3=0;for (int k = -nLen; k <= nLen; k++) {dWeight1 += pKernel[k + nLen];dWeight2 += pKernel[k + nLen];dWeight3 += pKernel[k + nLen];}if (src.type() == CV_8UC3) {for(int i=nLen;i<rows-nLen;i++){uchar *pData = src.ptr<uchar>(i);uchar *pDstData = dst.ptr<uchar>(i);//pData += 3*nLen; //pDstData += 3 * nLen;for (int j = nLen; j < cols - nLen; j++){dMul1 = dMul2 = dMul3 = 0;for (int k = -nLen; k <= nLen; k++) {dMul1 += pKernel[k + nLen] * (double)pData[(j + k) * 3];dMul2 += pKernel[k + nLen] * (double)pData[(j + k) * 3+1];dMul3 += pKernel[k + nLen] * (double)pData[(j + k) * 3+2];}pDstData[j*3] = (uchar)dMul1 / dWeight1;pDstData[j * 3+1] = (uchar)dMul2/ dWeight2;pDstData[j * 3+2] = (uchar)dMul3 / dWeight3;}}//xfor (int j = nLen; j<cols - nLen; j++) {uchar *pData = src.col(j).data;uchar *pDstData = dst.col(j).data;for (int i = nLen; i < rows - nLen; i++){dMul1 = dMul2 = dMul3 = 0;for (int k = -nLen; k <= nLen; k++) {dMul1 += pKernel[k + nLen] * (double)pDstData[(i + k) *cols * 3];dMul2 += pKernel[k + nLen] * (double)pDstData[(i + k) *cols * 3 + 1];dMul3 += pKernel[k + nLen] * (double)pDstData[(i + k) *cols * 3 + 2];}pDstData[i *cols * 3] = (uchar)(dMul1 / dWeight1);pDstData[i *cols * 3 + 1] = (uchar)(dMul2 / dWeight2);pDstData[i *cols * 3 + 2] = (uchar)(dMul3 / dWeight3);}}//y}delete[]pKernel;}