【WebGPU实时光追美少女】降噪与色调映射

少女dtysky

世界Skill

时刻2021.10.04

路标

本系列文章设计的所有代码均已开源，Github仓库在这里：dtysky/webgpu-renderer。

在上一篇文章BRDF与蒙特卡洛积分的最后，我们最终得到了一个充满了噪点的输出，接下来要考虑的就是如何降噪（Denoise），这个分为时间和空间两部分。并且由于整个光照计算都是在高动态范围（HDR）计算的，最后还要做一次色调映射（Tone Mapping）来将其转换到有限的低动态范围（LDR），便于显示器显示。

降噪

实时路径追踪的噪声来自于样本数量不足下的信息量不足，所以有两个方面可以来优化：

提高样本量。
尽可能使用有效的样本。

时间滤波（Temporal Filter）就是用来解决这两个问题的。而由于样本的增加总是需要累计时间，所以会有一个启动时间，所以为了尽可能加快速度达到一个相对可看的程度，我们在空间上也会做一些滤波，或者说是图像处理，这就是空间滤波（Spatial Filter）。

目前工程上的滤波方案往往是混合滤波，即结合时间空间二者，目前最常见的就是SVGF以及其变种，但由于本项目只涉及到静态场景，并未实现Motion Vector这种动态场景的技术，所以只是借鉴其中的一部分。

这里就直接用SVGF的一个流程图来论述整个流程：

首先求得当前帧的1SPP原始图像。
和上一次的结果作时间滤波。
将时间滤波的结果用于空间滤波。
空间滤波的输出反馈到下一帧，同时作为色调映射的来源，做最终输出。

时间滤波

时间滤波顾名思义，是通过不同帧之间的信息来做滤波。路径追踪的时间滤波比较简单：

C_{out} = \alpha C_{i-1} + (1 - \alpha)C_{i}

这个公式本质上是个简单的混合，将上一帧颜色和这一帧颜色通过权重$\alpha$混合起来。这个权重的值是动态变更的，我这里取：

\alpha = \frac{frameCount - 1}{frameCount}

这么做也是有原理可循的，回顾一下上一篇文章所说蒙特卡洛积分的采样相加求平均，我们需要相对保证每帧的结果被均匀平均：

\frac {C_{n}} n + \frac {n-1}n(\frac {C_{n-1}} {n-1} + \frac {n-2} {n-1}(\frac {C_{n-2}} {n-2} + ...) = \frac 1 n(C_{n} + C_{n - 1} + C_{n - 2} + ...)

所以只要在实现的Shader内接受完结传来的这个权重u_preWeight，然后在JS里不断更新权重即可：

[[stage(compute), workgroup_size(16, 16, 1)]]
fn main(
  [[builtin(workgroup_id)]] workGroupID : vec3<u32>,
  [[builtin(local_invocation_id)]] localInvocationID : vec3<u32>
) {
  let size: vec2<i32> = textureDimensions(u_current);
  let groupOffset: vec2<i32> = vec2<i32>(workGroupID.xy) * 16;
  let baseIndex: vec2<i32> = groupOffset + vec2<i32>(localInvocationID.xy);

  let pre: vec4<f32> = textureLoad(u_pre, baseIndex, 0);
  let current: vec4<f32> = textureLoad(u_current, baseIndex, 0);
  let mixed: vec4<f32> = vec4<f32>(mix(current.rgb, pre.rgb, vec3<f32>(material.u_preWeight)), current.a);

  textureStore(u_output, baseIndex, mixed);
}

空间滤波

空间滤波本质上就是图像处理，最简单的空间滤波就是高斯模糊。

高斯模糊

exp(-\frac {(x-x_c)^2 + (y-y_c)^2} {2\sigma_d^2})

其以某个像素为中心开一个固定大小“窗口”，根据周边像素到中心像素的半径计算权重，然后将窗口中的所有像素颜色加权平均，最终得到一个模糊的效果：

单纯的高斯模糊其实是一个低通滤波器，其本质是在频域将高频噪声过滤，留下低频信号，这种频域的处理就会得到一个如图的看起来很糊的效果，这个糊来源于边界的消失，而边界的消失则是基于以下两个事实：

高频信息并非都是噪声。
低频信息也并非不都是噪声。

所以我们需要又一个更好的策略，来在降噪的同时还能保留有效的高频信息，这就要提到联合双边滤波（Joint Bilateral Filter）。

联合双边滤波

exp({-\frac {(x-x_c)^2 + (y-y_c)^2} {2\sigma_d^2}} - {\frac {||color_{x,y} - color_{x_c,y_c}||^2} {2\sigma_c^2}})

变换后为：

exp({-\frac {(x-x_c)^2 + (y-y_c)^2} {2\sigma_d^2}}) * exp({-\frac {||color_{x,y} - color_{x_c,y_c}||^2} {2\sigma_c^2}})

如公式所示，在方才高斯模糊的基础上，我们加了一个和颜色相关的项，来共同决定窗口中某个像素的权重，同时还可以控制二者的方差来决定像素位置和颜色分别决定权重的比例。那么接下来就可以很自然想到——既然如此，是否可以再加上别的项，通过别的信息来共同决定这个权重呢？当然是可以的，这就又要用到Gbuffer了。

已知Gbuffer中提供了世界坐标、法线等等信息，我们可以考虑对于当前图像“边界”到底意味着什么？很自然可以想到就是一些合理可控的突变的点，比如在一个立方体的转折处（法线发生了较大变化），再比如一前一后的两个物体（深度发生了较大变化）。利用这些特性，就可以构造一个比较合理的滤波器：

exp({-\frac {(x-x_c)^2 + (y-y_c)^2} {2\sigma_d^2}}) * exp({-\frac {||color_{x,y} - color_{x_c,y_c}||^2} {2\sigma_c^2}}) * exp({-\frac {||normal_{x,y} - normal_{x_c,y_c}||^2} {2\sigma_n^2}}) * exp({-\frac {||z_{x,y} - z_{x_c,y_c}||^2} {2\sigma_n^2}})

然后接下来就是调整系数（方差）的工作了，我在CPU中通过方差事先算好了系数，来供Shader中使用：

export function genFilterParams(sigmas: Float32Array): Float32Array {
  const res = new Float32Array(sigmas.length);

  for (let i: number = 0; i < sigmas.length; i += 1) {
    const s = sigmas[i];
    res[i] = -0.5 * s * s;
  }

  return res;
}

到此就可以实现最终的滤波器了，联合双边滤波的结果如下，可见解决了边界被模糊的问题：

Outlier Removal

滤波滤波到此结束，但还有另一个问题，如下图：

图上有一些亮点，这些亮点的亮度明显高于周边。这是由于采样率不足，导致没有足够的样本做平均来达到真正的效果，换言之如果样本足够这些亮点是可以被平均掉的。但对于实时应用，我们等不起这个过程，所以必须用方法将它们“抹除”掉。这确实会造成能量的不守恒，但却是一个必要的权衡。

这个过程称为Outlier Removal，采用的是统计的方法：

在滤波操作之前进行。
计算所有点亮度的均值和标准差。
根据均值和方差算出上界，如果某像素亮度大于上界，则判定为Outlier。
针对Outlier，将其拉到上界范围内。

代码实现

至此，所有的要素都已集齐，就可以进行代码实现了（当然写的有点糙，凑合看吧）：

fn calcWeightNumber(factor: f32, a: f32, b: f32) -> f32 {
  return exp(factor * (a - b) * (a - b));
}

fn calcWeightVec2(factor: f32, a: vec2<i32>, b: vec2<i32>) -> f32 {
  let diff: vec2<f32> = vec2<f32>(a - b);
  return exp(factor * dot(diff, diff));
}

fn calcWeightVec3(factor: f32, a: vec3<f32>, b: vec3<f32>) -> f32 {
  let diff: vec3<f32> = a - b;
  return exp(factor * dot(diff, diff));
}

fn calcLum(color: vec3<f32>) -> f32 {
  return dot(color, vec3<f32>(0.2125, 0.7154, 0.0721));
}

fn blur(center: vec2<i32>, size: vec2<i32>) -> vec4<f32> {
  let radius: i32 = WINDOW_SIZE / 2;
  let sigmaD: f32 = material.u_filterFactors.x;
  let sigmaC: f32 = material.u_filterFactors.y;
  let sigmaZ: f32 = material.u_filterFactors.z;
  let sigmaN: f32 = material.u_filterFactors.w;
  var centerColor: vec4<f32> = textureLoad(u_preFilter, center, 0);
  let alpha: f32 = centerColor.a;
  let centerPosition = textureLoad(u_gbPositionMetalOrSpec, center, 0).xyz;
  let centerNormal = textureLoad(u_gbNormalGlass, center, 0).xyz;
  var colors: array<array<vec3<f32>, ${WINDOW_SIZE}>, ${WINDOW_SIZE}>;
  var lums: array<array<f32, ${WINDOW_SIZE}>, ${WINDOW_SIZE}>;

  var minUV: vec2<i32> = max(center - vec2<i32>(radius, radius), vec2<i32>(0));
  var maxUV: vec2<i32> = min(center + vec2<i32>(radius, radius), size);
  var localUV: vec2<i32> = vec2<i32>(0, 0);
  var sumLum: f32 = 0.;
  var count: f32 = 0.;

  // 计算每个像素的颜色亮度存起来
  for (var r: i32 = minUV.x; r <= maxUV.x; r = r + 1) {
    localUV.y = 0;
    for (var c: i32 = minUV.y; c <= maxUV.y; c = c + 1) {
      let iuv: vec2<i32> = vec2<i32>(r, c);
      let color: vec3<f32> = textureLoad(u_preFilter, iuv, 0).rgb;
      let lum: f32 = calcLum(color);
      colors[localUV.x][localUV.y] = color;
      lums[localUV.x][localUV.y] = lum;

      sumLum = sumLum + lum;
      count = count + 1.;
      localUV.y = localUV.y + 1;
    }
    localUV.x = localUV.x + 1;
  }

  // 计算亮度均值
  let meanLum: f32 = sumLum / count;

  // 计算亮度方差
  var std: f32 = 0.;
  for (var r: i32 = 0; r < localUV.x; r = r + 1) {
    for (var c: i32 = 0; c < localUV.y; c = c + 1) {
      let lum: f32 = lums[r][c];
      std = std + (lum - meanLum) * (lum - meanLum);
    }
  }
  std = sqrt(std / (count - 1.));

  // 确定亮度最大边界
  let largestLum: f32 = max(meanLum + std * 2., 1.);

  var lum: f32 = calcLum(centerColor.rgb);
  if (lum > largestLum) {
    centerColor = centerColor * largestLum / lum;
  }

  localUV = vec2<i32>(0, 0);
  var weightsSum: f32 = 0.;
  var res: vec3<f32> = vec3<f32>(0., 0., 0.);

  for (var r: i32 = minUV.x; r <= maxUV.x; r = r + 1) {
    localUV.y = 0;
    for (var c: i32 = minUV.y; c <= maxUV.y; c = c + 1) {
      var color: vec3<f32> = colors[localUV.x][localUV.y];
      lum = lums[localUV.x][localUV.y];

      if (lum > largestLum) {
        // 当前像素亮度大于最大边界，将其拉回来
        color = color * largestLum / lum;
      }

      let iuv: vec2<i32> = vec2<i32>(r, c);
      let position: vec4<f32> = textureLoad(u_gbPositionMetalOrSpec, iuv, 0);
      let normal: vec4<f32> = textureLoad(u_gbNormalGlass, iuv, 0);
      // 联合双边滤波
      let weight: f32 = calcWeightVec2(sigmaD, iuv, center)
        * calcWeightVec3(sigmaC, color.rgb, centerColor.rgb)
        * calcWeightVec3(sigmaN, normal.xyz, centerNormal.xyz)
        * calcWeightNumber(sigmaZ, position.z, centerPosition.z);
      weightsSum = weightsSum + weight;
      res = res + weight * color;

      localUV.y = localUV.y + 1;
    }
    localUV.x = localUV.x + 1;
  }

  res = res / weightsSum;

  return vec4<f32>(res, alpha);
}

色调映射

在完成了滤波后，如果不加任何处理，我们会得到这样的一个结果：

可见其有明显的过曝，所以就需要色调映射。评估一个色调映射质量最重要的维度是对比度的保留，历史上出现过很多种策略，像是Reinhard曲线、指数曲线等，但都不足够好，会有灰蒙蒙的感觉。直到后期，人们最终找到了高次曲线拟合的方法，最终实现了目前最通用的拟合策略，其计算高效，实现简单：

let A: f32 = 2.51;
let B: f32 = 0.03;
let C: f32 = 2.43;
let D: f32 = 0.59;
let E: f32 = 0.14;  

fn acesToneMapping(color: vec3<f32>) -> vec3<f32> {
  return (color * (A * color + B)) / (color * (C * color + D) + E); 
}

色调映射后的最终的结果如下：

结语

本篇文章所言都是针对静态场景的，在这种实现下，物体和相机都不能移动，这在实际的应用中显然是不够的。一般在真正可用的实现中，我们必须引入Motion Vector，来计算当前一个像素上一帧所对应的像素，然后分别存储权重，最终达到动态场景的滤波。但这也会引入更多其他的问题，由于篇幅和时间确实没空完成了。

当然如果只是玩票，理论上加上Motion Vector，对SVGF之类的进行实现也并非难事，这篇文章的很多知识是可以用到的，有兴趣的读者可以自行实现，或者等有生之年我进行优化（逃

点击查看评论