RealTime-Rendering11-Gaussian Blur

Gaussian Blur

Gaussian Function

高斯函数（或正态分布函数）在各个领域中被广泛应用。标准差 σ 决定了数据在均值周围的分散程度。较大的标准差会使数据分布得更分散。（体现在高斯模糊上就是图片会变的更糊）

高斯函数的函数图像类似钟形如下图所示：

In Compute Grahpic

在图像处理中，高斯函数通过与输入图像和高斯滤波器（核）进行二维卷积，用于平滑或模糊输入图像，或从源图像中去除噪声。在计算机图形学中，高斯模糊用于柔化阴影贴图的硬边缘，或创建泛光（bloom）效果。

Generate Gaussian Kernel

一般来说，二维卷积需要一个二维核（滤波器）。但如果核是可分离的，可以通过使用一维核进行两次一维卷积来优化。（使用 M×N 核的二维卷积需要 M×N 次乘法，而使用 M×1 和 1×N 可分离核的一维卷积只需要 M+N 次乘法。）
高斯滤波器在垂直和水平方向上都是对称函数，它是一个可分离核。因此，高斯模糊可以通过沿垂直和水平方向进行一维卷积来优化，而不是二维卷积。
基于标准差 σ 和kernelSize生成一维高斯卷积核：

function generateGaussianKernel(sigma, kernelSize)
{
  let kernel = new Float32Array(kernelSize);

  // compute kernel elements of Gaussian
  // do only half(positive side) and mirror to negative side
  // because Gaussian is even function, symmetric to Y-axis.
  let center = Math.floor(kernelSize / 2);   // center value of n-array(0 ~ n-1)

  let result = 0;
  let sum = 0;
  if(sigma == 0)
  {
    kernel.fill(0);
    kernel[center] = 1.0;
  }
  else
  {
    const SS2 = sigma * sigma * 2;
    kernel[center] = 1;
    sum = 1;
    for(let x = 1; x <= center; ++x)
      {
        // dividing (sqrt(2*PI)*sigma) is not needed because normalizing result later
        result = Math.exp(-(x*x)/SS2);
        kernel[center+x] = kernel[center-x] = result;
        sum += result;
        sum += result;
      }

    // normalize kernel
    // make sum of all elements in kernel to 1
    for(let i = 0; i <= center; ++i)
      kernel[center+i] = kernel[center-i] /= sum;
  }
  return kernel;
}

由于高斯核沿Y轴对称（偶函数），我们可以将数组大小缩减为仅存储正半轴的值（包括中心值）。然后，对负Y轴侧的数值进行镜像。

// generate half-size 1D seperable gaussian kernel
function generateHalfGaussianKernel(sigma, halfKernelSize)
{
    let kernel = new Float32Array(halfKernelSize);
    let result = 0;
    let sum = 0;
    if(sigma == 0)
    {
        kernel.fill(0);
        kernel[0] = 1.0;
    }
    else
    {
        const SS2 = sigma * sigma * 2;
        kernel[0] = 1;
        sum = 1;
        for(let x = 1; x < halfKernelSize; ++x)
        {
            result = Math.exp(-(x*x)/SS2);
            kernel[x] = result;
            sum += result * 2;
        }

        // normalize kernel
        for(let i = 0; i <= halfKernelSize; ++i)
            kernel[i] /= sum;
    }
    return kernel;
}

GaussianBlur Shader

// Gaussian Blur Fragment Shader
// constants
const int MAX_KERNEL = 21;              // max half kernel size including the center
const float ZERO = 0.0;
const float ONE  = 1.0;

// uniforms
uniform float kernel[MAX_KERNEL];       // half gaussian kernel from center
uniform vec2  imageSize;
uniform vec2  direction;                 // horizontal=(1,0) or vertical=(0,1)
uniform sampler2D map0;                 // input image

// input varying variables
varying vec2 texCoord0;

void main(void)
{
    // compute the center first
    vec3 color = texture2D(map0, texCoord0).rgb * kernel[0];
    vec2 offset;
    vec2 texelSize = 1.0f / imageSize;

    // compute with other kernel elements
    for(int i = 1; i < MAX_KERNEL; ++i)
    {
        offset = direction * float(i) * texelSize;
        color += texture2D(map0, texCoord0 + offset).rgb * kernel[i]; // positive side
        color += texture2D(map0, texCoord0 - offset).rgb * kernel[i]; // negative side
    }

    gl_FragColor = vec4(color, ONE);
}

Optimization

GLSL着色器中的卷积可以通过GPU的纹理线性过滤进一步优化。之前的卷积使用N×1的核对纹理采样进行N次乘法运算，其中Pi是一个纹素(texel)，Hi是一个核权重。
我们可以利用两个纹素之间的线性插值，这样只需要对一半的样本进行乘法运算。例如，openGL的双线性线性插值会对纹理样本P₀和P₁进行插值，返回P₀,₁：

因此，加权的2个样本：H₀P₀ + H₁P₁ 可以简化为 H₀,₁P₀,₁，其中：

使用线性插值优化后的高斯模糊Shader：

// Gaussian Blur Fragment Shader with texture filtering
// constants
const int MAX_KERNEL = 21;              // max half kernel size including the center
const float ZERO = 0.0;
const float ONE  = 1.0;

// uniforms
uniform float kernel[MAX_KERNEL];       // half gaussian kernel from center
uniform float imageDimension;
uniform vec2 direction;                 // horizontal=(1,0) or vertical=(0,1)
uniform sampler2D map0;                 // input image

// input varying variables
varying vec2 texCoord0;

void main(void)
{
    // compute the center texel first
    vec3 color = texture2D(map0, texCoord0).rgb * kernel[0];
    vec2 offset;

    // optimize using linear texture filtering (with half texels)
    float k;    // interpolated kernel, H = H1 + H2
    float t;    // interpolated alpha, t = H2 / H
    for(int i = 1; i < MAX_KERNEL; i += 2)
    {
        k = kernel[i] + kernel[i+1];
        t = kernel[i+1] / k;
        offset = direction * (float(i) + t) / imageDimension;
        color += texture2D(map0, texCoord0 + offset).rgb * k;
        color += texture2D(map0, texCoord0 - offset).rgb * k;
    }

    gl_FragColor = vec4(color, ONE);
}

Conclussion

优化技术	原理	效果
可分离滤波	将 2D 卷积拆分为两个 1D 卷积（水平+垂直）	复杂度从 $O(n^2)$ 降到 $O(2n)$
线性采样	利用 GPU 双线性过滤，一次采样覆盖两个像素	采样次数减半
降采样	先在低分辨率进行模糊，再放大	大幅减少计算量
限制核大小	使用较小的 σ 和核半径	权衡质量与性能