[图像处理] 高斯模糊的C++实现（Gaussian Blur）

程序完成时间：  2017/2/21

      事实上，写下这篇博文的时候，距离完成这个小程序已经两个月了。现在才跑过来填个坑儿。

参考文献：

          [1] Wiki百科：https://en.wikipedia.org/wiki/Gaussian_blur

          [2] 博客文章： http://www.ruanyifeng.com/blog/2012/11/gaussian_blur.html

实现效果截图：（图一为原图、图二为模糊结果。分割线用以调整处理区域）

          我的二女神！！西野七濑Nishino Nanase！！(づ￣ 3￣)づ

                   

                   

高斯模糊原理

       The Gaussian blur is a type of image-blurring filter that uses a Gaussian function (which also expresses the normal distribution in statistics) for calculating the transformation to apply to each pixel in the image. The equation of a Gaussian function in one dimension is:

       高斯模糊是一种图像模糊过滤器，通过使用高斯函数（在统计学中也被称为正态分布）来计算应用到图像中每个像素的变换关系（事实上是矩阵Matrix）。一维的高斯函数的表达式为：

       

         in two dimensions, it is the product of two such Gaussians, one in each dimension:

    在二维情况下，他是两个这样的高斯函数的积分，表达式如下：

    

      where x is the distance from the origin in the horizontal axis, y is the distance from the origin in the vertical axis, and σ is the standard deviation of the Gaussian distribution. When applied in two dimensions, this formula produces a surface whose contours are concentric circles with a Gaussian distribution from the center point. Values from this distribution are used to build a convolutionmatrix which is applied to the original image. Each pixel's new value is set to a weighted average of that pixel's neighborhood. The original pixel's value receives the heaviest weight (having the highest Gaussian value) and neighboring pixels receive smaller weights as their distance to the original pixel increases.

      这里，x是在水平方向上与图像原点（图像左上角）的距离，y是在竖直方向上与图像原点的距离，σ是高斯分布的标准差。当在二维情况下时，这个公式表现为一个轮廓由多个同心圆围绕一个中轴线组成的曲面（如下图）。这个分布的数值被用于构建作用在原始图像上的卷积矩阵。每个像素的新值被设置为临近的矩形区域像素的加权平均（权值来自于前面提到的的卷积矩阵）。某像素的原始值具有最高的权重（拥有一个最大的高斯函数值）并且临近像素的权值随着与该中心像素的距离的增加而减小。

    二维高斯分布的示意图如下（来源：WikiPedia）

    

高斯模糊的C++实现：

1. 设定初始值

        高斯模糊的卷积矩阵半径：r（即每个点的新值是由一个宽高均为2r+1的矩形区域的加权平均得出）

        标准差σ ： sigma（取卷积矩阵半径r的1/3为宜）

        下图为一个r=1的卷积矩阵，该矩阵可以通过二维的高斯函数带入，并将权值归一化求得。

        

2. 边界处理函数： int edge(int, int, int)

        之所以引入这个函数，是因为在高斯模糊中，每个点的新值都是由临近的矩形区域的平均值得出，但是若矩形区域超出图像范围，则会产生黑边。为了消除这一影响，我们采用对称的方法，将矩形的超出部分用内部的对应点代替。

int Image::edge(int x, int i, int w)
{
    // x为中心位置, i为偏移量,w为宽度
    int inx = x+i;
    if(inx<0 || inx>=w)
        return x-i;
}

3. 高斯函数处理

       这种方法非常简单，不做赘述，因为我想把更多的篇幅放在后面的高斯函数的优化版本中。

       在这里，我们只需要对每一个点都进行一次计算：用上述的卷积矩阵的每个数值分别与该点周围的相应矩形区域的像素值相乘，求出加权平均即可。但是每个点都需要进行(2r+1)*(2r+1)次运算，对于一个较大的图像来说，这种平方级的运算时间是难以忍受的。

4. 我们实际采用的高斯模糊的优化版本！！

    下面，我想重点谈一下关于高斯函数的优化！！

    在Wiki百科中有下面一段话：

      In practice, it is best to take advantage of the Gaussian blur’s separable property by dividing the process into two passes. In the first pass, a one-dimensional kernel is used to blur the image in only the horizontal or vertical direction. In the second pass, the same one-dimensional kernel is used to blur in the remaining direction. The resulting effect is the same as convolving with a two-dimensional kernel in a single pass, but requires fewer calculations.     

       ”在实际应用中，一种更好的方式为：对图像在两个独立的一维空间分别进行计算。在第一次处理中，应用一个一维的高斯函数得到的线性空间对图像在水平方向进行一次模糊。在第二次处理中，同样的一维线性空间对图像在竖直方向进行一次模糊。通过这种方式得到的结果可以达到与使用二维卷积空间处理的同样效果，但是却大大节省了计算时间”

         显然，通过这种方式，每个点只需要进行2*(2r+1)次计算即可，时间复杂度降低了一个幂。

初始化线性权重空间的代码：

// 一维高斯函数
double Image::gaussFunc1D(int x)
{
    double A = 1/(sigma*sqrt(2*3.141592653));
    double index = -1.0*(x*x)/(2*sigma*sigma);

    return A*exp(index);
}

// 获取线性的权值空间
QImage Image::getKernal(double* weight)
{ 
    weight= new double[2*r+1];
    double sum=0;

    // 获取权值空间weight[]
    for(int i=0; i<2*r+1;i++)
    {
        weight[i] = gaussFunc1D(i-radius);
        sum+=weight[i];
    }
    // 归一化
    for(int i=0; i<2*radius+1;i++)
    {
        weight[i]/=sum;
    }
}

进行两个方向的高斯模糊获得最终图像代码（Qt）:

// 获取最终的模糊图像
QImage Image::getBlurImage()
{ 
    // 原图为 originImg
    QImage originImg(......);
    // 第一个方向处理的图像为tmpImg
    QImage tmpImg(.......);
    // 经过第二次处理的最终结果newImg
    QImage newImg(.......);
    
    // 在横向进行一次相加
    for(int y=0;y<originImg.height();y++)
    {
        for(int x=0;x<originImg.width();x++)
        {
            double red=0,green=0,blue=0;
            for(int i=-r; i<=r; i++)
            {
                // 边界处理后的对应的权值矩阵实际值
                int inx = edge(x,i,width());

                QColor rgb = QColor(originImgpixel(inx,y));
                red+=rgb.red()*weight[r+i];
                green+=rgb.green()*weight[r+i];
                blue+=rgb.blue()*weight[r+i];
            }
            tmpImg.setPixel(x,y,qRgb(red,green,blue));
        }
    }
    // 在纵方向对第一次的结果重新进行一次
    for(int y=0;y<tmpImg.height();y++)
    {
        for(int x=0;x<tmpImg.width();x++)
        {
            double red=0,green=0,blue=0;
            for(int i=-r; i<=r; i++)
            {
                int iny = edge(y,i,height());

                QColor rgb = QColor(tmpImg.pixel(x,iny));
                red+=rgb.red()*weight[r+i];
                green+=rgb.green()*weight[r+i];
                blue+=rgb.blue()*weight[r+i];
            }
            newImg.setPixel(x,y,qRgb(red,green,blue));
        }
    }
    return newImg;
}