高斯滤波

Posted 2020-10-19 心要让你听见

高斯滤波

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

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

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

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

第四个问题：C代码实现？

 

 由于笔者使用的是整数模板，先将浮点数模板代码列出请参考http://blog.csdn.net/crzy_sparrow/article/details/6998745

1. //  一维高斯分布函数，用于平滑函数中生成的高斯滤波系数  
2. /*  
3.  *  @parameter sigma:     高斯函数参数  
4.  *  @parameter pdKernel:    高斯核函数模板  
5.  *  @parameter pnWidowSize:  高斯模板大小  
6.  */  
7. void CreatGauss(double sigma, double **pdKernel, int *pnWidowSize)  
8. {  
9.     LONG i;  
10.       
11.     int nCenter;//数组中心点  
12.     double dDis;//数组中一点到中心点距离  
13.   
14.     //中间变量  
15.     double dValue;  
16.     double dSum;  
17.     dSum = 0;  
18.       
19.     *pnWidowSize = 1+ 2*ceil(3*sigma);// [-3*sigma,3*sigma] 以内数据，会覆盖绝大部分滤波系数  
20.   
21.     nCenter = (*pnWidowSize)/2;  
22.     //生成高斯数据  
23.     for(i=0;i<(*pnWidowSize);i++)  
24.     {  
25.         dDis = (double)(i - nCenter);  
26.         dValue = exp(-(1/2)*dDis*dDis/(sigma*sigma))/(sqrt(2*3.1415926)*sigma);  
27.         (*pdKernel)[i] = dValue;  
28.         dSum+=dValue;  
29.     }  
30.     //归一化  
31.     for(i=0;i<(*pnWidowSize);i++)  
32.     {  
33.         (*pdKernel)[i]/=dSum;  
34.     }  
35. }  

 

1. //用高斯滤波器平滑原图像  
2. /* 
3.  *  @parameter sz   :  图像尺寸 
4.  *  @parameter pGray   :   图像灰度值 
5.  *  @parameter pResult:  图像 
6.  *  @parameter sigma:     高斯函数参数 
7.  */  
8. void GaussianSmooth(SIZE sz, LPBYTE pGray, LPBYTE pResult, double sigma)  
9. {  
10.     LONG x, y;  
11.     LONG i;  
12.       
13.     int nWindowSize;//高斯滤波器长度  
14.       
15.     int nLen;//窗口长度  
16.       
17.     double *pdKernel;//一维高斯滤波器  
18.       
19.     double dDotMul;//高斯系数与图像数据的点乘  
20.   
21.     double dWeightSum;//滤波系数总和  
22.   
23.     double *pdTemp;  
24.   
25.     nWindowSize = 1+ 2*ceil(3*sigma);// [-3*sigma,3*sigma] 以内数据，会覆盖绝大部分滤波系数  
26.   
27.   
28.     if ((pdTemp = (double *)malloc(sz.cx*sz.cy*sizeof(double)))==NULL)  
29.     {  
30.         printf("melloc memory for pdTemp failed!!");  
31.         exit(0);  
32.     }  
33.     if ((pdKernel = (double *)malloc(nWindowSize*sizeof(double)))==NULL)  
34.     {  
35.         printf("malloc memory for pdKernel,failed!!");  
36.         exit(0);  
37.     }  
38.   
39.     //产生一维高斯数据  
40.     CreatGauss(sigma, &pdKernel, &nWindowSize);  
41.   
42.     nLen = nWindowSize/2;  
43.   
44.     //x方向滤波  
45.     for(y=0;y<sz.cy;y++)  
46.     {  
47.         for(x=0;x<sz.cx;x++)  
48.         {  
49.             dDotMul = 0;  
50.             dWeightSum = 0;  
51.             for(i=(-nLen);i<=nLen;i++)  
52.             {  
53.                 //判断是否在图像内部  
54.                 if((i+x)>=0 && (i+x)<sz.cx)  
55.                 {  
56.                     dDotMul+=(double)(pGray[y*sz.cx+(i+x)] * pdKernel[nLen+i]);  
57.                     dWeightSum += pdKernel[nLen+i];  
58.                 }  
59.             }  
60.             pdTemp[y*sz.cx+x] = dDotMul/dWeightSum;     
61.         }  
62.     }  
63.   
64.     //y方向滤波  
65.     for(x=0; x<sz.cx;x++)  
66.     {  
67.         for(y=0; y<sz.cy; y++)  
68.         {  
69.             dDotMul = 0;  
70.             dWeightSum = 0;  
71.             for(i=(-nLen);i<=nLen;i++)  
72.             {  
73.                 if((i+y)>=0 && (i+y)< sz.cy)  
74.                 {  
75.                     dDotMul += (double)pdTemp[(y+i)*sz.cx+x]*pdKernel[nLen+i];  
76.                     dWeightSum += pdKernel[nLen+i];  
77.                 }  
78.             }  
79.             pResult[y*sz.cx+x] = (unsigned char)dDotMul/dWeightSum;  
80.         }  
81.     }  
82.       
83.     free(pdTemp);//释放内存  
84.     free(pdKernel);