202 系列完结撒花 ✿✿ヽ(°▽°)ノ✿！
0x00 To begin with
这篇文章将会包含以下内容：
- 部分课程内容回顾
- GAMES202 作业 5
For reference👇：
- 📺B 站视频： GAMES202-高质量实时渲染
- 📦代码仓库： congyuxiaoyoudao/GAMES202_Homework at working
- 作业5发布公告及场景下载
- 📖 GAMES202-作业5：实时光线追踪降噪
master 分支上是 202 全部作业汇总，working 分支用于提交代码。需要原始作业可下载 master 分支的包
0x01 问题修复 TroubleShooting
做作业时遇到一些出现问题的地方，这里一并提一下：
- 运行 build.bat 后 build 目录下没有出现可执行文件
ReadFloatImage方法报错- mathutil.h 中
uint32_twas not declared in this scope
第一个问题纯属对 Cmake 不熟练，build.bat 提供的命令只是生成工程文件，要运行构建还需要如下命令（或者直接用 VS 等 IDE 运行）：
1cmake --build build
这里为了简化运行操作，新增了一个 run.bat：
1cmake --build build
2cd build
3Denoise.exe
4cd ..
第二个问题为路径错误，下载完场景后我的根目录文件结构如下：
1assignment5/
2├── examples/
3| └── ...
4├── src/
5| ├── main.cpp
6│ └── ...
7└── ...
所以需要在 main.cpp 中更改一下路径（或者直接把 examples 放到 src 下）：
1filesystem::path inputDir("../examples/box/input");
2filesystem::path outputDir("../examples/box/output");
第三个问题，缺少头文件，在 mathutil.h 中包含 <cstdint> 即可。
0x02 联合双边滤波 Joint Bilateral Filter
对于单帧图像的降噪，实际上就是低通滤波的过程，但是常见的低通滤波核（例如高斯核）会无条件地抹除输入图像中的高频部分，包括图像中不同物体的边界。而我们希望在降噪的过程中能够保留这些边缘，避免 Overblur，就需要寻找别的方法。
双边滤波提供了一个不错的思路：双边滤波核的权重不仅考虑了像素间几何距离的差异，还考虑了图像值域上的差异：
$$
w(i, j, k, l) = \exp\left( -\frac{(i - k)^2 + (j - l)^2}{2\sigma_d^2} - \frac{\|I(i, j) - I(k, l)\|^2}{2\sigma_r^2} \right)
$$
当两个像素在值（颜色、深度等）上差异过大时，采样像素对中心像素的贡献就减少，这样就可以维持住边缘信息。
所以我们可以在双边滤波的基础上，为核加入更多控制以应对多种复杂边缘情况，这就是联合双边滤波的思想。
使用联合双边滤波对单帧图像进行降噪，作业中给出的核为：
$$
J(\mathbf{i}, \mathbf{j}) = \exp\left( -\frac{\|\mathbf{i} - \mathbf{j}\|^2}{2\sigma_p^2} - \frac{\|\tilde{C}[i] - \tilde{C}[j]\|^2}{2\sigma_c^2} - \frac{D_{\text{normal}}(\mathbf{i}, \mathbf{j})^2}{2\sigma_n^2} - \frac{D_{\text{plane}}(\mathbf{i}, \mathbf{j})^2}{2\sigma_d^2} \right)
$$
其中：
- $\lVert \mathbf{i} - \mathbf{j} \rVert$：像素坐标差值
- $\lVert \tilde{C}[i] - \tilde{C}[j] \rVert$：颜色差值
- $D_{\text{normal}}(\mathbf{i}, \mathbf{j}) = \arccos\left( \text{Normal}[\mathbf{i}] \cdot \text{{Normal}}[\mathbf{j}] \right)$，描述面朝向差异
- $D_{\text{plane}}(\text{i}, \text{j}) = \text{Normal}[\text{i}] \cdot \frac{\text{Position}[\text{j}] - \text{Position}[\text{i}]}{\lVert \text{Position}[\text{j}] - \text{Position}[\text{i}] \rVert}$，描述距离差异（沿法向变化的距离）
一直有点疑惑为什么很多作业实现都是把第一项当作世界位置差异，试了一下效果也的确差不多
各个 $\sigma$ 为控制不同项的平滑程度的标准差，$\sigma$ 越大表示对上述距离差异越不敏感，即更能容忍权重差异。
在 denoise.cpp 中完成 Filter 函数，按照上述公式完成联合双边滤波逻辑：
1Buffer2D<Float3> Denoiser::Filter(const FrameInfo &frameInfo) {
2 int height = frameInfo.m_beauty.m_height;
3 int width = frameInfo.m_beauty.m_width;
4 Buffer2D<Float3> filteredImage = CreateBuffer2D<Float3>(width, height);
5 int kernelRadius = 16;
6#pragma omp parallel for
7 for (int y = 0; y < height; y++) {
8 for (int x = 0; x < width; x++) {
9 // Joint bilateral filter
10 // get current pixel vals
11 Float3 color = frameInfo.m_beauty(x, y);
12 Float3 normal = frameInfo.m_normal(x, y);
13 Float3 pos = frameInfo.m_position(x, y);
14
15 float weight_sum = 0.f;
16 Float3 color_sum = Float3(0.f);
17 for(int i = -kernelRadius; i <= kernelRadius; i++){
18 for(int j = -kernelRadius; j <= kernelRadius; j++){
19 // clamp coordinate
20 int sampleCoord_x = std::min(std::max(0,x+i),width-1);
21 int sampleCoord_y = std::min(std::max(0,y+j),height-1);
22
23 // sample sampled pixel vals
24 Float3 sampledColor = frameInfo.m_beauty(sampleCoord_x,sampleCoord_y);
25 Float3 sampledNormal = frameInfo.m_normal(sampleCoord_x,sampleCoord_y);
26 Float3 sampledPos = frameInfo.m_position(sampleCoord_x,sampleCoord_y);
27
28 // calculate weight
29 float weight = 0.0;
30 float d_coord2 = i*i+j*j;
31
32 float d_color2 = SqrDistance(color,sampledColor);
33
34 float cosTheta = Dot(normal,sampledNormal);
35 float d_normal2 = SafeAcos(cosTheta)*SafeAcos(cosTheta);
36
37 float d_plane = 0.f;
38 float d_pos = Distance(pos, sampledPos);
39 if (d_pos > 0)
40 d_plane = Dot(normal,Normalize(sampledPos-pos));
41 float d_plane2 = d_plane*d_plane;
42
43 float w_coord = - d_coord2/(2*m_sigmaCoord*m_sigmaCoord);
44 float w_color = - d_color2/(2*m_sigmaColor*m_sigmaColor);
45 float w_normal = - d_normal2/(2*m_sigmaNormal*m_sigmaNormal);
46 float w_plane = - d_plane2/(2*m_sigmaPlane*m_sigmaPlane);
47
48 weight = exp(w_coord+w_color+w_normal+w_plane);
49 // accumulate weight and color buffer
50 weight_sum += weight;
51 color_sum += sampledColor * weight;
52 }
53 }
54 // protect 0 weight
55 if(weight_sum == 0.0f){
56 filteredImage(x, y) = color;
57 }
58 else{
59 filteredImage(x, y) = color_sum / weight_sum;
60 }
61 }
62 }
63 return filteredImage;
64}
0x03 反投影 Reprojection
反投影利用当前帧的信息找到上一帧的对应像素位置，通过反投影可以在时间上复用已经滤波的上一帧信息。
已知当前帧的像素坐标，可以在 GBuffer 上获取世界位置和各类变换矩阵，使用下式即可通过像素世界位置找到上一帧该像素的屏幕坐标：
$$
\text{Screen}_{i-1} = P_{i-1} V_{i-1} M_{i-1} M_i^{-1} \, \text{World}_i
$$
投影到屏幕坐标下会存在一些屏幕空间固有的问题，导致反投影后的像素不合法，若强行复用会导致 lagging ，对于这些情况我们要手动进行判断：
- 反投影后的像素在屏幕外，即当前帧像素对应新进入屏幕的物体
- 反投影后的像素和当前帧像素对应物体的 id 不同
第一种情况可以通过判断反投影后的像素坐标是否合法屏幕范围内（$[0,w-1],[0,h-1]$）；第二种情况可以通过判断 GBuffer Id 是否相等。仅对上述判断后合法的像素进行时间累积，更新 m_valid 和 m_misc。
按照上述思路完成 Reprojection 函数：
1void Denoiser::Reprojection(const FrameInfo &frameInfo) {
2 int height = m_accColor.m_height;
3 int width = m_accColor.m_width;
4 Matrix4x4 preWorldToScreen =
5 m_preFrameInfo.m_matrix[m_preFrameInfo.m_matrix.size() - 1];
6 Matrix4x4 preWorldToCamera =
7 m_preFrameInfo.m_matrix[m_preFrameInfo.m_matrix.size() - 2];
8#pragma omp parallel for
9 for (int y = 0; y < height; y++) {
10 for (int x = 0; x < width; x++) {
11 // Reproject
12 m_valid(x, y) = false;
13 m_misc(x, y) = Float3(0.f);
14
15 float id = frameInfo.m_id(x, y);
16 // skip ivalid id
17 if(id == -1)
18 continue;
19
20 Float3 pos = frameInfo.m_position(x, y);
21 Matrix4x4 WorldToLocal = Inverse(frameInfo.m_matrix[id]);
22 Matrix4x4 preLocalToWorld = m_preFrameInfo.m_matrix[id];
23 // P_{i-1}V_{i-1}M_{-1}M_{i}^{-1}
24 Matrix4x4 ReprojectionMatrix = preWorldToScreen*preLocalToWorld*WorldToLocal;
25
26 Float3 prePosScreen = ReprojectionMatrix(pos, Float3::EType::Point);
27 float pre_x = prePosScreen.x;
28 float pre_y = prePosScreen.y;
29
30 // is in viewport
31 if(pre_x<0 || pre_x>width-1 || pre_y<0 || pre_y>height-1)
32 continue;
33
34 float preId = m_preFrameInfo.m_id(pre_x, pre_y);
35
36 if(preId == id){
37 m_valid(x, y) = true;
38 m_misc(x, y) = m_accColor(pre_x, pre_y);
39 }
40 }
41 }
42 std::swap(m_misc, m_accColor);
43}
0x04 时间累积 Temporal Accumulation
反投影将当前帧所有像素对应的上一帧反投影的像素颜色（若合法）存入 m_accColor，以便时间累积复用。
如果直接使用上一帧的信息，当其颜色与当前颜色差异过大时可能导致 ghosting，为了避免这个情况，需要对上一帧信息和当前帧信息做钳制，即：
$$
C^{(i-1)}=\text{clamp}(C^{(i-1)}, \mu - k\sigma, \mu + k\sigma)
$$
其中，$\mu,\sigma$ 分别为当前帧像素某一邻域的均值和方差。钳制后即可使用下式混合颜色：
$$
C^{(i)} = \alpha \bar{C}^{(i)} + (1 - \alpha) C^{(i-1)}
$$
其中，$\bar{C}^{(i)}$ 为当前帧已降噪像素信息，$C^{(i-1)}$ 为上一帧已降噪反投影像素信息，$\alpha$ 通常取 0.2，表示更多地使用上一帧信息。
在 denoise.cpp 中完成 TemporalAccumulation 函数，通过遍历当前像素的 7×7 邻域计算均值和方差，并将上一帧信息限制到 $(\mu-k\sigma,\mu+k\sigma)$ 范围内。为了避免求方差再做一遍循环，可以利用方差的另一个计算公式：
$$
Var(X)=E(X^{2})-E(X)^{2}
$$
1void Denoiser::TemporalAccumulation(const Buffer2D<Float3> &curFilteredColor) {
2 int height = m_accColor.m_height;
3 int width = m_accColor.m_width;
4 int kernelRadius = 3;
5 int kernelSize = 2*kernelRadius+1;
6 int sampleNum = kernelSize * kernelSize;
7#pragma omp parallel for
8 for (int y = 0; y < height; y++) {
9 for (int x = 0; x < width; x++) {
10 // Temporal clamp
11 float alpha = 1.0f;
12 // previous frame color
13 Float3 color = m_accColor(x, y);
14
15 if(m_valid(x, y)){
16 alpha = m_alpha;
17 // calculate mu and sigma in 7*7 block
18 Float3 mu, mu2 = Float3(0.f);
19 Float3 sigma = Float3(0.f);
20 for(int i = -kernelRadius; i <= kernelRadius; i++){
21 for(int j = -kernelRadius; j <= kernelRadius; j++){
22 int sampleCoord_x = std::min(std::max(0,x+i),width-1);
23 int sampleCoord_y = std::min(std::max(0,y+j),height-1);
24
25 Float3 curSampledColor = curFilteredColor(sampleCoord_x, sampleCoord_y);
26 mu+=curSampledColor;
27 mu2+=Sqr(curSampledColor);
28 }
29 }
30 mu /= sampleNum;
31 mu2 /= sampleNum;
32 sigma = SafeSqrt(mu2 - mu*mu);
33
34 color = Clamp(color, mu-sigma*m_colorBoxK, mu+sigma*m_colorBoxK);
35 }
36
37 // Exponential moving average
38 m_misc(x, y) = Lerp(color, curFilteredColor(x, y), alpha);
39 }
40 }
41 std::swap(m_misc, m_accColor);
42}
0x05 A-Torus Wavelet 加速
A-Torus Wavelet，即进行多次滤波，每次使用逐渐增加的滤波核尺寸和步长，以保证每次滤波的采样点数不变。总共进行 5 次滤波，每次滤波输入图像为上一次滤波结果，第 $i$ 次滤波的步长为 $2^{i}$，核尺寸为 $2^{i+2}+1$ 。通过 A-Torus Wavelet 的加速，可以将原本 $64^{2}$ 的采样数减少为 $125$，大幅提升单帧降噪速度。
在 Denoise 类中新增一个 ATWFilter 函数，实现由 A-Trous Wavelet 加速的滤波，在原本的两层 for 循环外部再嵌套一层 for：
1Buffer2D<Float3> Denoiser::ATWFilter(const FrameInfo &frameInfo) {
2 int height = frameInfo.m_beauty.m_height;
3 int width = frameInfo.m_beauty.m_width;
4 Buffer2D<Float3> filteredImage = CreateBuffer2D<Float3>(width, height);
5 int passes = 5;
6#pragma omp parallel for
7 // Accelerated by A-Trous Wavelet
8 for (int y = 0; y < height; y++) {
9 for (int x = 0; x < width; x++) {
10 // Joint bilateral filter
11 // get current pixel vals
12 Float3 color = frameInfo.m_beauty(x, y);
13 Float3 normal = frameInfo.m_normal(x, y);
14 Float3 pos = frameInfo.m_position(x, y);
15
16 float weight_sum = 0.f;
17 Float3 color_sum = Float3(0.f);
18
19 for(int pass = 0; pass < passes; pass++){
20 int stride = pow(2, pass);
21 int kernelRadius = stride * 2;
22 for(int i = -kernelRadius; i <= kernelRadius; i+=stride){
23 for(int j = -kernelRadius; j <= kernelRadius; j+=stride){
24 // clamp coordinate
25 int sampleCoord_x = std::min(std::max(0,x+i),width-1);
26 int sampleCoord_y = std::min(std::max(0,y+j),height-1);
27
28 // sample sampled pixel vals
29 Float3 sampledColor = frameInfo.m_beauty(sampleCoord_x,sampleCoord_y);
30 Float3 sampledNormal = frameInfo.m_normal(sampleCoord_x,sampleCoord_y);
31 Float3 sampledPos = frameInfo.m_position(sampleCoord_x,sampleCoord_y);
32
33 // calculate weight
34 float weight = 0.0;
35 float d_coord2 = i*i+j*j;
36
37 float d_color2 = SqrDistance(color,sampledColor);
38
39 float cosTheta = Dot(normal,sampledNormal);
40 float d_normal2 = SafeAcos(cosTheta)*SafeAcos(cosTheta);
41
42 float d_plane = 0.f;
43 float d_pos = Distance(pos, sampledPos);
44 if (d_pos > 0)
45 d_plane = Dot(normal,Normalize(sampledPos-pos));
46 float d_plane2 = d_plane*d_plane;
47
48 float w_coord = - d_coord2/(2*m_sigmaCoord*m_sigmaCoord);
49 float w_color = - d_color2/(2*m_sigmaColor*m_sigmaColor);
50 float w_normal = - d_normal2/(2*m_sigmaNormal*m_sigmaNormal);
51 float w_plane = - d_plane2/(2*m_sigmaPlane*m_sigmaPlane);
52
53 weight = exp(w_coord+w_color+w_normal+w_plane);
54 // accumulate weight and color buffer
55 weight_sum += weight;
56 color_sum += sampledColor * weight;
57 }
58 }
59 // protect 0 weight
60 if(weight_sum == 0.0f){
61 filteredImage(x, y) = color;
62 }
63 else{
64 filteredImage(x, y) = color_sum / weight_sum;
65 }
66 }
67 }
68 }
69 return filteredImage;
70}
加速效果也是十分明显的，ATWFilter 比 Filter 几乎快了两分钟，但是就单帧降噪质量来说是不如 Filter 的。
0x06 Afterwords
所以闫大的离线渲染课什么时候开 :p