高斯函数（Gaussian function）的详细分析

转自：https://blog.csdn.net/qinglongzhan/article/details/82348153

摘要

    论文中遇到很重要的一个元素就是高斯核函数，但是必须要分析出高斯函数的各种潜在属性，本文首先参考相关材料给出高斯核函数的基础，然后使用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 logarithm a 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图像，首先贴出完整代码：

function mainfunc()  
% 测试高斯函数，递增的方法实现高斯函数参数的改变对整个高斯函数的影响，  
% 并自动保存为gif格式输出。  
% created by zhao.buaa 2016.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文件  
  
% main loop  
for i = 1:item  
    center  = 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([xmin xmax ymin ymax zmin zmax])  
    axis([0 ksize 0 ksize 0 0.008]);  view([0, 90]);% 改变可视角度     
    % ksize 递增  
%     ksize = ksize + 10;  
    % sigma 递增  
    sigma = sigma + dt;       
      
    % 自动保存为gif图像  
    frame = getframe(i);  
    im = frame2im(frame);  
    [I,map] = rgb2ind(im,256);  
    if i==1  
        imwrite(I,map,filename,'gif','Loopcount',inf, 'DelayTime',0.4);  
    else  
        imwrite(I,map,filename,'gif','WriteMode','append','DelayTime',0.4);  
    end  
end;  
  
close all;  
  
  
%% 截断高斯核函数，截断的程度取决于参数bias  
function ksigma = ksigma(ksize, sigma, center,bias)  
%ksize = 80;    sigma = 15;  
ksigma=fspecial('gaussian',ksize, sigma);   % 构建高斯函数  
[m, n] =size(ksigma);  
for i = 1:m  
    for j = 1:n  
        if(  (i<center-bias)||(i>center+bias)||(j<center-bias)||(j>center+bias)  )  
            ksigma(i,j) = 0;  
        end;  
    end;  
end; 

结果图：

固定ksize为20，sigma从1-9，固定center在高斯中间，并且bias偏移量为整个半径，即原始高斯函数。

随着sigma的增大，整个高斯函数的尖峰逐渐减小，整体也变的更加平缓，则对图像的平滑效果越来越明显。

保持参数不变，对上述高斯函数进行截断，即truncated gaussian function，bias的大小为ksize*3/10，则结果如下图：

truncated gaussian function的作用主要是对超过一定区域的原始图像信息不再考虑，这就保证在更加合理的利用靠近高斯函数中心点的周围像素，同时还可以改变高斯函数的中心坐标，如下图：

为了便于观察截断的效果，改变了可视角度。

高斯核函数卷积

    论文中使用gaussian与feature map做卷积，目前的结果来看，要做到随着到边界的距离改变高斯函数的截断参数，因为图像的边缘如果使用原始高斯函数，就会在边界地方出现特别低的一圈，原因也很简单，高斯函数在与原始图像进行高斯卷积的时候，图像边缘外为0计算的，那么如何解决边缘问题呢？

先看一段代码：

% 截断高斯核函数
ksize = 40; ksigma=fspecial('gaussian',  ksize, 6);
[m, n] =size(ksigma);
for i = 1:m
    for j = 1:n
        if( i<25 )
           ksigma(i,j) = 0;
        end;
    end;
end;
figure, mesh(ksigma);

在i，即row上对高斯核函数进行截断，bias为半径大小，则如图6

并且对下图7进行卷积，

首先卷积核为原始未截断高斯核函数，则结果如图8

可以看出，在图像边缘处的卷积结果出现不想预见的结果，边缘处的值出现大幅度减少的情况，这是高斯核函数在边缘处将图像外的部分当成0计算的结果，因此，需要对高斯核函数进行针对性的截断处理，但是前提是要掌握bias的规律，下面就详细分析。

如果使用图6的高斯核函数与图7做卷积操作，则如图9：

可以看出，相比较于图8，与高斯核函数相对应的部分出现了变化，也就是说：

if( i<25 )
  ksigma(i,j) = 0;
end;

靠近边缘的时候，改变 i 或 j 的值，即可保证边缘处的平滑处理。但是这样改变高斯核函数，使用matlab不是很好解决这个问题，还是使用将待处理图像边缘向外部扩展bias的大小，与标准高斯核函数做卷积，再将超过原始图像大小的部分剪切掉，目前来看在使用matlab中imfilter函数做卷积运算最合适且最简单的处理方法了，先写在这里，此部分并不是论文中的核心部分，只是数值运算的技巧性编程方法。

设计图像处理

深入地探讨高斯滤波与图像处理的关联，包括参数的影响、参数的选取、高斯模板的形成以及自行编程实现高斯滤波的效果与openCV函数实现效果比对。

说道“sigma表示的是标准差，如果标准差比较小，这是就相当于图像点运算，则平滑效果不明显；反之，标准差比较大，则相当于平均模板，比较模糊”，那么这么说可能很多人包括一开始的我并不是很理解，这是为什么呢，那么我们需要从高斯函数谈起：

这样一个高斯函数的概率分布密度如下图所示：

我们要理解好这个图，横轴表示可能得取值x，竖轴表示概率分布密度F(x)，那么不难理解这样一个曲线与x轴围成的图形面积为1。sigma（标准差）决定了这个图形的宽度，我给出下述结论：sigma越大，则图形越宽，尖峰越小，图形较为平缓；sigma越小，则图形越窄，越集中，中间部分也就越尖，图形变化比较剧烈。这其实很好理解，如果sigma也就是标准差越大，则表示该密度分布一定比较分散，由于面积为1，于是尖峰部分减小，宽度越宽（分布越分散）；同理，当sigma越小时，说明密度分布较为集中，于是尖峰越尖，宽度越窄！

理解好上述结论之后，那么（一）中的结论当然也就顺理成章了，sigma越大，分布越分散，各部分比重差别不大，于是生成的模板各元素值差别不大，类似于平均模板；sigma越小，分布越集中，中间部分所占比重远远高于其他部分，反映到高斯模板上就是中心元素值远远大于其他元素值，于是自然而然就相当于中间值得点运算。

程序也可以验证如下：

窗口尺寸：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(int j = 0; j < img->weigth; ++j)

                               {

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

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

                                             {

                                                         ...

                                             }

                               }

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

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

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

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

                               {

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

                                             {

                                                         ...

                                             }

                               }

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

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

                               {

                                       for(int m = 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控制中的选择积分的方法（具体名字忘了···囧）。这个我将在后面进行适当地尝试！

    

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