快速高斯模糊算法的原理以及webgl工程实现

[forthealllight]

快速高斯模糊算法的原理以及webgl工程实现

---

    在日常生活中，高斯模糊很常见，相比与马赛克的方式，高斯模糊给人的感觉更加的自然。高斯模糊顾名思义，就是卷积(掩膜)的采样服从高斯分布，本文主要聊一聊，如何对一张图片进行高斯模糊，以及如何减少高斯模糊的算法计算量。此外，在前端对图片的渲染中，使用webgl可以充分利用GPU来进行计算和渲染，可以起到加速的效果，本文也会在快速高斯模糊算法的基础上，介绍用webgl如何实现快速高斯模糊。

• 滤波和卷积(掩膜)
• 高斯模糊
• 快速高斯模糊原理

---

一、滤波和卷积

    首先介绍一下滤波和卷积，滤波是信号处理中的一个概念，对于信号，可以由很多不同频率的波组成，滤波的本质就是增强或者减弱某一个频率的波，信号经过滤波处理，得到的新的信号。

    简单来说滤波就是：从信号中得到指定波长的波

(1)图像的时域和频域

    那么什么是图片或者说图像中的滤波呢，首先我们要明白滤波其实是一个频域的概念，频域顾名思义就是频率，那么什么是图像的频率呢？

• 时域

    要了解图像的频域首先先介绍一下图像的时域。

    一幅数字图像可定义为一个二维函数f(x,y)，其中x和y是空间（平面）坐标，而在任何一对空间坐标（x,y）出的幅值f称为图像在该点处的灰度或强度。

    上述图片中我们肉眼看到的灰度或者说颜色，就是时域，也就是说我们平常说的某个像素的颜色，就是从时域观测到的。

• 频域

    图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。简单来说：频域就是灰度或者颜色变换的趋势(梯度)。

    同样我们以上面那幅图为例，我们用rgb颜色画出图片每一个点的颜色。可以得到如下图：

    颜色波动较大点就是频率较大的地方，相对的没有颜色波动的地方频率较小。

(2)卷积(掩膜)

    我们可以通过傅立叶变换将时域转换成频域，时域的卷积等于频的乘积。如果滤波可以在频域做，那么只需要做乘积。但是从时域到频域的转换并不是很直观，是相对复杂的，因此，在工程上对于图像的滤波大部分集中在时域，也就是对颜色值做卷积。

    我们将原始图像，用一个窗口做卷积，得到的新的图像就是滤波后的图像，这个窗口就称为掩膜。这个窗口每个元素的值一般与滤波函数有关，我们接着举例来说明滤波的过程：

上述就是一个简单的原始图像，图像上面的值表示灰度值，我们现在来假设一个掩膜。

我们将图像的灰度与窗口最卷积，这里为了计算方便，只计算高亮的一个像素点，其原始像素值为98。卷积后计算得到滤波后的图像为如下所示：

    我们发现像素98灰度变成了82。这是一个像素，将图片的所有像素都与窗口做卷积，得到的新的灰度，就是图片经过滤波后得到的新图片。

注意：当然为了保证灰度或者说颜色值的范围不会超过阈值，窗口的值必须是归一化的。

二、高斯模糊

(1)什么是高斯模糊

    模糊就是一种特殊的滤波，经过这种滤波后图像变得不清晰。我们知道滤波 = 原始图像和掩膜的卷积，当掩膜(窗口)服从高斯分布时，此时我们称这种滤波为高斯滤波，也称为高斯模糊。

    首先来看一维的高斯分布：

其概率图像为：

在图像中，我们是在一个2D的坐标系内，因此应该服从的是一个二维的概率密度函数：

    在二维的高斯分布函数中:

• x 表示的是一个二维的坐标(x,y)。
• μ表示的是中心点的一般会假设称(0,0)。

    在计算掩膜的时候，因为高斯模糊的本质是用周围的点来估计中心点，因此我们将需要被估计的那个中心点的坐标设置为(0,0)。其周边N个点的值，用高斯分布计算出来。就得到了滤波的窗口函数。我们以3x3滤波窗口的计算为例：我们取一个像素点为原点，然后采样其周边8个点，总共9个点。

    我们以中心点为原点，共9个点，用二维高斯分布计算其概率，取方差 σ为1。得到9个点的高斯分布的值为：

    归一化这9个点，可以得到完整的滤波窗口(掩膜)：

    这就是最终的掩膜或者说滤波窗口，将原始图像的每一个像素与该掩膜卷积，就得到了一个模糊的图片，这个过程就称之高斯模糊或者说高斯滤波。

    概括的讲，高斯模糊跟高斯分布不是一回事，当滤波的窗口服从高斯分布，此时的滤波称为高斯模糊。前置的数学知识就告一段路，接着我们来看工程上如何实现和优化高斯模糊。

(2)高斯模糊的实现

    从二维高斯模糊的公式中我们可以看出来，对于高斯模糊而言，在计算窗口掩膜的时候只有一个 变量就是模糊半径和方差。对于方差我们都会设为1，模糊半径R决定了采样多少个点，前面的例子中模糊半径为3，所以采样了9个点。在二维中模糊半径（也就是取多少个点周围）可以为：

采样半径	采样点个数
3	9
5	25
7	49
9	81
11	121
13	169

    我们以半径为3的9点采样为例，来看看高斯模糊的webgl代码实现：

      vec4 orColor=texture2D(u_image,v_texCoord);
      const float PI = 3.1415926535897932384626433832795; 
      float w = 1./degree;
      float temWeight[11];
      temWeight[0] = (1.0/(2.0 * PI * w * w));
      temWeight[1] = (1.0/(2.0 * PI * w * w)) * exp(-1 * w * w);
      temWeight[2] = (1.0/(2.0 * PI * w * w)) * exp(-2 * w * w);
      float total;
      total = temWeight[0] + temWeight[1] * 4. + temWeight[2] * 4. 
      float orAlpha=orColor.a;
      float weight[3];
      weight[0] = temWeight[0]/total; 
      weight[1] = temWeight[1]/total; 
      weight[2] = temWeight[2]/total;
      vec4 color= texture2D(u_image,v_texCoord + tex_offset * vec2( 0,0)) * weight[0] 
      + texture2D(u_image,v_texCoord + tex_offset * vec2( 0,-1)) * weight[1] 
      + texture2D(u_image,v_texCoord + tex_offset * vec2( 1,0)) * weight[1] 
      + exture2D(u_image,v_texCoord + tex_offset * vec2( -1,0)) * weight[1] 
      + texture2D(u_image,v_texCoord + tex_offset * vec2( 1, 0)) * weight[1] 
      + texture2D(u_image,v_texCoord + tex_offset * vec2( 1, 1)) * weight[2] 
      + texture2D(u_image,v_texCoord + tex_offset * vec2(-1, 1)) * weight[2] 
      + texture2D(u_image,v_texCoord + tex_offset * vec2( 1, -1)) * weight[2] 
      + texture2D(u_image,v_texCoord + tex_offset * vec2( -1, -1)) * weight[2] 

     gl_FragColor=vec4(color.rgb,orAlpha);

    上述就是webgl版的9点采样高斯模糊实现方式，上述着色器的写法 ，就可以对一张图片做高斯模糊。效果图如下所示：

    在9点采样的时候，对于图片中的每一个像素点，做了9次乘法，8次加法（计算total的不算和归一化权重的部分不算，是通用部分）。*这样对于一张400300像素的图片，一共需要做400300（8+9） 共近20万次计算**。虽然我们用webgl利用gpu提高了运算效率，但是某些场景特别是实时性较高或者FPS较高的动画中，如此庞大的计算量是不能接受的，因此我们需要通过一定的方法来减少计算量。

三、快速高斯模糊原理

    高斯模糊的计算量减化基本上围绕一下两个假设进行的：

1. 二维的高斯模糊可以看成两个一维高斯模糊的乘积
2. 标准差2σ的高斯模糊等于2个标准差为σ的高斯模糊之和

（1） 杨辉三角(二项式系数)来代替权重求值

    我们前面在计算权重的时候是用二维高斯分布来计算的，这个计算的过程可以简化成用二项式 系数来表示。原理就是：

高斯分布是正太分布，而离散的高斯分布可以用二项式分布 的二项式系数的值来近似离散高斯分布的采样权重。

    上述就是一个二项式系数的图片，同样的当半径为3，进行9点采样的时候，我们取N = R+ 1 = 4.此时的二项式系数就 为1 4 6 4 1。权重就一目了然，不需要跟之前一样经过计算高斯分布的系数来求权重。

（2）二维高斯滤波拆成两个一维高斯滤波

    我们可以将二维的高斯滤波拆成两个一维的高斯滤波。同样对于半径为3的高斯模糊，此时我们需要采样的不是9个点，而是 5个点。

    我们分别求x方向和y方向的一维的高斯滤波，这样9次乘法和8次 加法就变成了5次乘法和4次加法，几乎减少了一半的计算量。

（3）线性滤波进一步减少计算

    我们以半径为R的25点采样为例，通过将二维高斯模糊拆成2个一维高斯模糊，实现的代码如下：

uniform sampler2D image;

out vec4 FragmentColor;

uniform float offset[5] = float[](0.0, 1.0, 2.0, 3.0, 4.0);
uniform float weight[5] = float[](0.2270270270, 0.1945945946, 0.1216216216,
                                  0.0540540541, 0.0162162162);

void main(void) {
    FragmentColor = texture2D(image, vec2(gl_FragCoord) / 1024.0) * weight[0];
    for (int i=1; i<5; i++) {
        FragmentColor +=
            texture2D(image, (vec2(gl_FragCoord) + vec2(0.0, offset[i])) / 1024.0)
                * weight[i];
        FragmentColor +=
            texture2D(image, (vec2(gl_FragCoord) - vec2(0.0, offset[i])) / 1024.0)
                * weight[i];
    }
}

    因为是半径为5，25个点的采样，降维后 还需要实现9次乘法，8次加法。此时我们可以通过线性滤波的方法，用一个点来表示两个点。

    线性滤波、线性插值是数学上一个常用的方法。我们简单来介绍一下线性插值：

    上图中我们已知(x0,y0)和(x1,y1)。当我们如果已知x,如何求y,这就需要用到线性插值。

*y = y0 + [(y1-y0)/(x1-x0)] x**

    也就是说我们可以通过2个点，线性差值的方法得到一个新的点。线性插值在图像中广泛使用，比如

    上述一个五彩的正方形，其实除了顶点的4个颜色，其他点的颜色都是通过双线性插值得到的。

    我们可以由于2个点得到一个新点，这就是线性插值。这个新点可以表示原来的两个 点。

    在高斯模糊中我们可以用上述公式，将2个点的权重和距离用一个点来表示。通过此公式，我们半径为R = 5，25点的采样可以进一步简化。从线性的 5个不同权重，变成 3个权重。

uniform sampler2D image;

out vec4 FragmentColor;

uniform float offset[3] = float[](0.0, 1.3846153846, 3.2307692308);
uniform float weight[3] = float[](0.2270270270, 0.3162162162, 0.0702702703);

void main(void) {
    FragmentColor = texture2D(image, vec2(gl_FragCoord) / 1024.0) * weight[0];
    for (int i=1; i<3; i++) {
        FragmentColor +=
            texture2D(image, (vec2(gl_FragCoord) + vec2(0.0, offset[i])) / 1024.0)
                * weight[i];
        FragmentColor +=
            texture2D(image, (vec2(gl_FragCoord) - vec2(0.0, offset[i])) / 1024.0)
                * weight[i];
    }
}

    就简化成了5次 乘法，4次 加法。

（4）用空间换时间

    除此之外，我们也可以进一步减少乘法的次数，对于半径为3，通过拆成两个一维的高斯滤波后，由9次乘法、8次加法变成了5次乘法和4次加法，我们可以通过查表的方法，将乘法给简化。

• 对于半径为3的高斯模糊，窗口掩膜的系数是确定的。简化后只有2个不同的值。
• 对于像素而言，其像素值的范围也是确定的，0-255.

    这样我们就可以 创建一个256 * 3的 map对象，保存在数组中，这样就不用去计算乘法，可以直接拿到掩膜后的值。这样就省去了5次乘法，进一步变成了只有4次加法。

（5）最后FBO多进行几次

    通过快速高斯滤波方法得到的高斯模糊可能效果不是很好，我们此时可以在第一次滤波的结果上，再一次滤波，多次滤波的方法来提升高斯模糊的效果。

[forthealllight]

σ $\overrightarrow{μ} $

​

[MrQinYQ]

大佬请教几个问题： temWeight[0] = (1.0/(2.0 PI w w)); temWeight[1] = (1.0/(2.0 PI w w)) exp(-1 w w); temWeight[2] = (1.0/(2.0 PI w w)) exp(-2 w w); 根据公式 下面两个是不是应该分别是 exp(-0.5 w w);和exp(-1 w w);呢？ 我还有一个疑问我看exp(0) exp(-0.5) exp(-1)就是图片里对应的1，0.606，0.368。公式前面的(1.0/(2.0 PI w w))是什么意义呢？ w就是标准差吧？degree的含义是什么呢？

[yangnaiyue]

自动回复邮件。您好，我已收到您的来信，会尽快给您回复。