高斯函数的详细分析[通俗易懂]-高斯函数的性质

高斯函数（Gaussian function）的详细分析[通俗易懂]摘要  论文中遇到很重要的一个元素就是高斯核函数，但是必须要分析出高斯函数的各种潜在属性，本文首先参考相关材料给出高斯核函数的基础，然后使用matlab自动保存不同参数下的高斯核函数的变化gif动图，同时分享出源代码，这样也便于后续的论文写作。高斯函数的基础2.1一维高斯函数高斯函数，GaussianFunction，也简称为Gaussian，一维形式如下：对于任…

摘要

论文中遇到很重要的一个元素就是高斯核函数，但是必须要分析出高斯函数的各种潜在属性，本文首先参考相关材料给出高斯核函数的基础，然后使用matlab自动保存不同参数下的高斯核函数的变化gif动图，同时分享出源代码，这样也便于后续的论文写作。

高斯函数的基础

2.1 一维高斯函数

高斯函数，Gaussian Function， 也简称为Gaussian，一维形式如下：

对于任意的实数a,b,c，是以著名数学家Carl Friedrich Gauss的名字命名的。高斯的一维图是特征对称“bell curve”形状，a是曲线尖峰的高度，b是尖峰中心的坐标，c称为标准方差，表征的是bell钟状的宽度。

高斯函数广泛应用于统计学领域，用于表述正态分布，在信号处理领域，用于定义高斯滤波器，在图像处理领域，二维高斯核函数常用于高斯模糊Gaussian Blur，在数学领域，主要是用于解决热力方程和扩散方程，以及定义Weiertrass Transform。

从上图可以看出，高斯函数是一个指数函数，其log函数是对数凹二次函数whose logarithma concave quadratic function。

高斯函数的积分是误差函数error function，尽管如此，其在整个实线上的反常积分能够被精确的计算出来，使用如下的高斯积分

同理可得

当且仅当

上式积分为1，在这种情况下，高斯是正态分布随机变量的概率密度函数，期望值μ=b，方差delta^2 = c^2，即

2.2 二维高斯函数

二维高斯函数，形如

A是幅值，x。y。是中心点坐标，σxσy是方差，图示如下，A= 1,xo= 0,yo= 0, σx= σy= 1

2.3 高斯函数分析

这一节使用matlab直观的查看高斯函数，在实际编程应用中，高斯函数中的参数有

ksize 高斯函数的大小

sigma 高斯函数的方差

center 高斯函数尖峰中心点坐标

bias 高斯函数尖峰中心点的偏移量，用于控制截断高斯函数

为了方便直观的观察高斯函数参数改变而结果也不一样，下面的代码实现了参数的自动递增，并且将所有的结果图保存为gif图像，首先贴出完整代码：

functionmainfunc()%测试高斯函数，递增的方法实现高斯函数参数的改变对整个高斯函数的影响，%并自动保存为gif格式输出。%createdbyzhao.buaa2016.09.28%%保存gif动画item=10;%迭代次数dt=1;%步长大小ksize=20;%高斯大小sigma=2;%方差大小%filename=['ksize-'num2str(ksize)'--'num2str(ksize+dt*item)'-sigma-'num2str(sigma)'.gif'];%必须预先建立gif文件filename=['ksize-'num2str(ksize)'-sigma-'num2str(sigma)'--'num2str(sigma+dt*item)'.gif'];%必须预先建立gif文件%mainloopfori=1:itemcenter=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([xminxmaxyminymaxzminzmax])axis([0ksize0ksize00.008]);view([0,90]);%改变可视角度%ksize递增%ksize=ksize+10;%sigma递增sigma=sigma+dt;%自动保存为gif图像frame=getframe(i);im=frame2im(frame);[I,map]=rgb2ind(im,256);ifi==1imwrite(I,map,filename,'gif','Loopcount',inf,'DelayTime',0.4);elseimwrite(I,map,filename,'gif','WriteMode','append','DelayTime',0.4);endend;closeall;%%截断高斯核函数，截断的程度取决于参数biasfunctionksigma=ksigma(ksize,sigma,center,bias)%ksize=80;sigma=15;ksigma=fspecial('gaussian',ksize,sigma);%构建高斯函数[m,n]=size(ksigma);fori=1:mforj=1:nif((i<center-bias)||(i>center+bias)||(j<center-bias)||(j>center+bias))ksigma(i,j)=0;end;end;end;

是否还在为Ide开发工具频繁失效而烦恼，来吧关注以下公众号获取最新激活方式。亲测可用！

为防止网络爬虫，请关注公众号回复”口令”

激活idea 激活CLion DataGrip DataSpell dotCover dotMemory dotTrace GoLand PhpStorm PyCharm ReSharper ReShaC++ Rider RubyMine WebStorm 全家桶 解锁下载刷新

【正版授权，激活自己账号】：Jetbrains全家桶Ide使用，1年售后保障，每天仅需1毛

【官方授权 正版激活】：官方授权 正版激活 自己使用，支持Jetbrains家族下所有IDE…

窗口尺寸：3*3 sigma = 0.1

窗口尺寸：3*3 sigma = 0.8

窗口尺寸：3*3 sigma =2

接着，我们来重点讨论下高斯模板，在初学高斯滤波的时候，用得最多的也是最经典的一个3*3模板就是

[1 2 1 ]

[2 4 2 ]

[1 2 1]

,当时我就很纳闷这个模板是怎么出来的，后来我经过多方查找资料，基本得到了如下的解释：高斯模板实际上也就是模拟高斯函数的特征，具有对称性并且数值由中心向四周不断减小，这个模板刚好符合这样的特性，并且非常简单，容易被大家接受，于是就比较经典！但是这样一个简单的矩阵是远远不能满足我们对图像处理的要求的，我们需要按照自己的要求得到更加精确的模板，那么接下来我们就编程实现自己想要的高斯模板。部分关键函数如下：
double** createG(int iSize, double sigma)
{
double **guass;
double sum = 0;
double x2 = 0;
double y2 = 0;
int center = (iSize – 1) / 2;
guass = new double*[iSize];//注意，double*[k]表示一个有10个元素的指针数组

for (int i = 0; i < iSize; ++i)
{
guass[i] = new double[iSize];
}
for (int i = 0; i<iSize; i++)
{//使用x2，y2降低了运算速度，提高了程序的效率
x2 = pow(double(i – center), 2);
for (int j = 0; j<iSize; j++)
{
y2 = pow(double(j – center), 2);
sum += guass[i][j] = exp(-(x2 + y2) / (2 * sigma*sigma));
}
}
if (sum != 0)
{
//归一化
for (int i = 0; i<iSize; i++)
{
for (int j = 0; j<iSize; j++)
{
guass[i][j] /= sum;
}
}
}
return guass;
}

上述的这个函数就是返回的一个用户想要的高斯模板，其输入参数iSize表示模板大小（这里先只考虑方阵模板），输入参数sigma表示高斯函数标准差，聪明的你应该一眼就看出这个程序实际上就是按照公式根据规定的矩阵大小进行离散化得到相应的数据。纵观整个程序，我觉得最不好理解的是那个参数center，这是个什么参数呢？通过程序可以看出为什么center与iSize呈2center+1的关系呢？而且，为什么x2，y2是那样取值呢？这实际上是模拟的一种距离关系。举个例子来说：

假设我现在想要使用窗口尺寸为2k+1*2k+1的高斯模板对点进行操作，那么假设这个高斯模板的第一个元素地址记为为（1,1），那么高斯模板中心元素很容易算出为（k，k）。假设现在在计算模板中（i，j）元素值，那么公式中的x2模拟为该点到中心点的距离，即i-k-1（为什么要-1，读者不妨自己自己写个3*3的矩阵，看看（1,1）元素到（2,2）元素是不是要走两步）。但是程序中是没有-1的，这是因为程序中的i是从0开始的！于是这个center参数实际上就是代表的中心点！

运行结果：

我的程序运行的结果：3*3，sigma=1

Matlab的程序运行的结果：3*3，sigma=1

可以看出，相差无几，这个程序是成功的！

然后，最重要的部分当然是使用上述高斯模板对图像进行滤波处理得到想要的效果，那么接下来重点论述滤波处理。要理解好高斯滤波，自己写高斯滤波的算法当然是最好的，然后再和openCV的函数进行效果比对，这样，算是对高斯滤波有了比较好的认识和理解了！

在很多资料上，我们都看到高斯函数的这样一个特性，可分离性，意思是一个二维的高斯函数可以分解成相同的一维高斯函数处理两遍，得到的效果是一样的，但是处理时间却大大缩短了。

我们可以想到，二维函数直接处理会是这样的：

for(int i = 0; i < img->height; ++i)

for(intj = 0;j < img->weigth; ++j)

{

for(intm = 0; m < iSize; ++m)

for(intn= 0;n< iSize; ++n)

{

…

}

}

这样的算法复杂度可以看出是为O（height*weigth*iSize^2）

然而使用分解后的一维函数进行两次处理，程序应当如下：

for(int i = 0; i < img->height; ++i)

for(intj = 0;j < img->weigth; ++j)

{

for(intm = 0; m < iSize; ++m)

{

…

}

}

for(intj = 0; j< img->weigth; ++j)

for(inti= 0;i< img->height; ++i)

{

for(intm = 0; m < iSize; ++m)

{

…

}

}

这样的算法复杂度为O（2*height*weigth*iSize），比一维处理要少了很多，所以时间对应来说也会快一点。

用后者，我们的关键函数如下：

/*
生成一维高斯模板，水平的和垂直方向上的模板是一样的
输入参数分别是：模板大小，sigma值
输出：一维数组（高斯模板）
*/
double* CreateMuban(int iSize ,double sigma)
{
double *gauss = new double[iSize];//声明一维模板
int radius = (iSize – 1) / 2;//这是高斯半径
double MySigma = 2 * sigma * sigma;
double value = 0;
for (int i = 0; i < iSize; i++)
{//高斯函数前面的常数项因为在归一化的时候将会消去，故这里不重复计算
gauss[i] = exp(-(i – radius)*(i – radius)/MySigma);
value = value + gauss[i];
}
for (int i = 0; i < iSize; i++)
{//归一化
gauss[i] = gauss[i] / value;
}
return gauss;
}

//对像素进行操作
IplImage* operatorImage(IplImage* img, double* Muban, int iSize)
{
//创建一张新的图片来进行滤波操作
IplImage* NewImage = cvCreateImage(cvSize(img->width, img->height), 8, 3);
int radius = (iSize – 1) / 2;
int r = 0;
int g = 0;
int b = 0;
CvScalar cs;
//复制图片
cvCopy(img, NewImage);
//先对I，也就是垂直方向进行操作
for (int j = 0; j < NewImage->width; ++j)
{
for (int i = 0; i < NewImage->height; ++i)
{
//先判断是否是边缘,不是则操作，是则跳过不处理，保持原样
if (!JudgeEdge(i, NewImage->height, radius))
{
for (int k = 0; k < iSize; ++k)
{
/*b = b + (int)((double)(NewImage->imageData + (i – radius + k) * NewImage->widthStep)[j * (int)NewImage->nChannels + 0] * Muban[k]);
g = g + (int)((double)(NewImage->imageData + (i – radius + k) * NewImage->widthStep)[j * (int)NewImage->nChannels + 1] * Muban[k]);
r = r + (int)((double)(NewImage->imageData + (i – radius + k) * NewImage->widthStep)[j * (int)NewImage->nChannels + 2] * Muban[k]);*/

cs = cvGet2D(NewImage, i – radius + k, j); //获取像素
b = b + (int)((double)cs.val[0] * Muban[k]);
g = g + (int)((double)cs.val[1] * Muban[k]);
r = r + (int)((double)cs.val[2] * Muban[k]);

}
/*((uchar *)(NewImage->imageData + i * NewImage->widthStep))[j * NewImage->nChannels + 0] = b; //改变该像素B的颜色分量
((uchar *)(NewImage->imageData + i * NewImage->widthStep))[j * NewImage->nChannels + 1] = g; //改变该像素G的颜色分量
((uchar *)(NewImage->imageData + i * NewImage->widthStep))[j * NewImage->nChannels + 2] = r; //改变该像素R的颜色分量 */
cs = cvGet2D(NewImage, i, j);
cs.val[0] = b;
cs.val[1] = g;
cs.val[2] = r;
cvSet2D(NewImage, i, j, cs);
b = 0;
g = 0;
r = 0;
}
}
}
//在对J，也就是水平方向进行操作
for (int i = 0; i < NewImage->height; ++i)
{
for (int j = 0; j < NewImage->width; ++j)
{
//先判断是否是边缘,不是则操作，是则跳过不处理，保持原样
if (!JudgeEdge(j, NewImage->width, radius))
{
for (int k = 0; k < iSize; ++k)
{
/*b = b + (int)((double)(NewImage->imageData + i * NewImage->widthStep)[(j – radius + k) * (int)NewImage->nChannels + 0] * Muban[k]);
g = g + (int)((double)(NewImage->imageData + i * NewImage->widthStep)[(j – radius + k) * (int)NewImage->nChannels + 1] * Muban[k]);
r = r + (int)((double)(NewImage->imageData + i * NewImage->widthStep)[(j – radius + k) * (int)NewImage->nChannels + 2] * Muban[k]);*/

cs = cvGet2D(NewImage, i, j – radius + k); //获取像素
b = b + (int)((double)cs.val[0] * Muban[k]);
g = g + (int)((double)cs.val[1] * Muban[k]);
r = r + (int)((double)cs.val[2] * Muban[k]);

}
/*((uchar *)(NewImage->imageData + i * NewImage->widthStep))[j * NewImage->nChannels + 0] = b; //改变该像素B的颜色分量
((uchar *)(NewImage->imageData + i * NewImage->widthStep))[j * NewImage->nChannels + 1] = g; //改变该像素G的颜色分量
((uchar *)(NewImage->imageData + i * NewImage->widthStep))[j * NewImage->nChannels + 2] = r; //改变该像素R的颜色分量*/

cs = cvGet2D(NewImage, i, j);
cs.val[0] = b;
cs.val[1] = g;
cs.val[2] = r;
cvSet2D(NewImage, i, j, cs);
b = 0;
g = 0;
r = 0;
//cout << r << ” ” << g << ” ” << b << endl;
}
}
}

return NewImage;
}

实现效果：

自己编的程序与openCV函数处理效果并无差别，成功！

最后，对高斯滤波部分进行扩展——自适应高斯滤波。为了达到更好的滤波效果，我们的手法需要更加灵活，往往可以加上一些限制条件进行判断是否需要处理。这一点很想PID控制中的选择积分的方法（具体名字忘了···囧）。这个我将在后面进行适当地尝试！

学习高斯滤波只是一个开始，但是学习其他的图像处理方法诸如此类，我要学会举一反三，而且要实践与理论紧密结合，真正做到知其甚解才好啊！