作业要求
本次作业，我们需要实现一个简单的实时光线追踪的降噪方法。光线追踪的渲染结果，G-Buffer 及其他相关信息会以文件的形式提供给大家。本次作业的重点在于对渲染结果的降噪，而不是如何进行渲染。为了在作业框架中实现上述目标，我们将作业分为三个部分：单帧图像的降噪，计算 motion vector，累积帧间信息。
根据上述作业要求，我们把具体实现分为三个步骤
- 单帧降噪
- 累加上一帧
- 用A-Trous Wavelet加速单帧降噪（提高项）
单帧降噪
首先我们来看一下单帧降噪的效果。
未降噪
降噪后
观察上面两幅图我们发现，第二张图的效果明显要比第一张图的好，但是同样的，第二张图片球的边缘相对于第一张要模糊了。由此我们很容易想到，降噪可以通过滤波来做，但是滤波存在一个很明显的问题就是，普通的滤波方式得到的结果是全局模糊的。我们想保留一些细节，在保留边缘的情况下减少噪声，我们不得不提一种特殊的滤波方式：联合双边滤波
联合双边滤波
首先了解一下滤波的原理：高斯滤波以采样距离的远近来决定该采样点对主像素的贡献 ，距离主像素点越近的采样点，贡献到主像素点的颜色的比例就越大。
我们来看看普通滤波的结果，普通高斯滤波会丢失细节
我们再看看双边滤波的结果
上面的结果在保留了边界细节的同时模糊了内容，双边滤波正是做这个的，我们可以通过加入一些限定条件，当达到某些条件的时候，减少某些像素对中心像素的影响，就可以保留想要的细节。
举个例子，我们看下图的山峰边缘，蓝色突变到黑色，他们的颜色变换十分剧烈
在通过距离判断贡献的基础上，我们加一个条件：颜色差值越大，那么该采样点对中心像素的贡献就越小，这样就可以把距离和颜色一起考虑进去，得到一张颜色相近的位置模糊，边缘处清晰的图片了，公式参考下图
i、j、k、l分别是中心像素点和采样点的x、y坐标
I指的是颜色
要得到我们期望的结果，只通过距离和颜色的限制是不够的，我们通过多方因素来限制滤波，因此我们引入联合双边滤波。
联合双边滤波可以参考G-buffers进行滤波，就比如通过对比像素点的深度图、Albedo图和法线贴图来对滤波进行限制，这样可以使得滤波的结果保留更多的细节
作业中给出了公式，在该公式中我们根据distance（距离）、color（颜色）、normal（法线）、plane（综合考虑法线和平面的夹角）四个因素来限制滤波。
在得知了如何计算每个像素点对中心像素点的贡献后，我们参考以下公式，把各个贡献累加起来并且归一化。
作业中给出了公式的详细计算方法，我们直接套入公式计算
//denoiser.cpp
Buffer2D<Float3> Denoiser::Filter(const FrameInfo &frameInfo) {
int height = frameInfo.m_beauty.m_height;
int width = frameInfo.m_beauty.m_width;
Buffer2D<Float3> filteredImage = CreateBuffer2D<Float3>(width, height);
int kernelRadius = 32;
#pragma omp parallel for
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
// TODO: Joint bilateral filter
//filteredImage(x, y) = frameInfo.m_beauty(x, y);
//
//定义Filter范围
int x_start = std::max(0, x - kernelRadius);
int x_end = std::min(width-1,x+kernelRadius);
int y_start = std::max(0,y-kernelRadius);
int y_end = std::min(height - 1, y + kernelRadius);
//获取Filter中心点的数据
auto center_position = frameInfo.m_position(x, y);
auto center_normal = frameInfo.m_normal(x, y);
auto center_color = frameInfo.m_beauty(x, y);
//先定义 最终颜色值和权重总和方便归一化
Float3 result_color;
auto total_weight = 0.0f;
for (int i = x_start; i < x_end; i++) {
for (int j = y_start; j < y_end; j++) {
//对采样点进行分析
auto current_position = frameInfo.m_position(i, j);
auto current_normal = frameInfo.m_normal(i, j);
auto current_color = frameInfo.m_beauty(i, j);
//套公式算
// Distance
auto D = Dot((center_position - current_position),
(center_position - current_position)) /
(2.0f * m_sigmaCoord * m_sigmaCoord);
// Color
auto C = Dot((center_color - current_color),
(center_color - current_color)) /
(2.0f * m_sigmaColor * m_sigmaColor);
// Normal
////SafeAcos和acos(x)
///一样，不过如果x超出范围的话，它会夹紧最近的可用值。
//和标准c函数acos(x)一致，这个值返回0到pi。
auto N = SafeAcos(Dot(center_normal, current_normal));
N *= N;
//这里遇到一个非常奇怪的问题，下面直接除结果是正确的，但是我改成N/=(2 * m_sigmaNormal * m_sigmaNormal)就会报错，希望有大佬能解答
N/(2 * m_sigmaNormal * m_sigmaNormal);
// Plane
auto P = Dot(center_normal,(current_position - center_position)) *
Dot(center_normal,
(current_position - center_position)) /
(2 * m_sigmaPlane * m_sigmaPlane);
float weight = std::exp(-D - C - N - P);
total_weight += weight;
result_color += current_color * weight;
}
}
//filteredImage(x, y) = frameInfo.m_beauty(x, y);
//filteredImage(x, y) =(result_color.x / total_weight, result_color.y / total_weight,result_color.z / total_weight);
//归一化
filteredImage(x, y) = result_color / total_weight;
}
}
return filteredImage;
}
上面通过采样周围的像素并计算各个像素对中心像素的贡献来滤波得出单帧也就是当前帧的颜色值，相比于1spp的光追结果，降噪后的结果清晰了很多，但是时间是原来的好几倍。
累加上一帧
单帧降噪的结果非常的不错，但是只经过单帧降噪不能消除所有的早点，降噪说白了还是把早点的颜色按一定比例分给旁边的像素，所以一次滤波会出现“不干净”的结果，由此我们可以再对画面进行一次滤波，不过这次滤波我们不再是求周围像素点对中心点的贡献了，我们计算上一帧对应位置对这一帧某个像素点的贡献。
假如上一帧是没有噪声的，我们令最终颜色为80%上一帧颜色+20%当前帧滤波后颜色，这样可以得到一个非常好的效果，这便是基于时间上的滤波Temporal。
对时间的滤波Temporal
Temporal的核心思想是运用了递归，当前帧的前一帧一定是已经滤波好了的。我们假设每一帧的上一帧都是没有噪声的，那么得到的这一帧结果会比只做一次当前帧滤波的结果好很多，接着这一帧的结果作为下一帧的前一帧继续递归，会得到一个非常好的结果。
我们详细讲解Temporal的原理。
motion vector
在Temporal里，我们需要用motion vector 来记录上一帧对应物体的位置。
motion vector是用来记录某个物体这一帧相对上一帧的运动情况的
我们可以从motion vector里找到该像素上一帧的已经滤波好了的对应位置的像素，然后把上一帧已经滤波过的结果加权作用到这一帧上。
这被称为时间上的复用。
要用motion vector做到当前帧某个位置找到前一帧的某个位置，就需要用到G-Buffers中的一些数据
我们可以用G-Buffer存储世界坐标coord也可以用MVPE(E是Viewport)矩阵来逆变换位置来找到该点的世界坐标，通过乘一个Translate的逆矩阵来得到上一帧的世界坐标，再通过记录上一帧的MVP矩阵得到上一帧的对应屏幕位置。
该方法和光流法有一定区别，该方法可以用G-Buffers的信息获取一个准确的结果，但光流法只通过两帧的结果进行分析，不使用G-Buffers。
下面得出利用motion vector的公式
~代表没有滤波的结果
-代表已滤波的结果
- 首先要对当前帧进行滤波处理
- 然后最后该帧的结果其实是前一帧的已经滤波的结果和这一帧已经滤波的结果的混合
我们用Temporal降噪后的结果
这里补充一点：
第二张图是降噪后的图片，看起来比第一张亮很多，其实这种理解是不对的，因为降噪过程并不会导致能量增加，导致这种问题的原因是第一张图的噪点的能量实际上是非常高的，往往他的亮度超过了1，也就是说第二张图的能量来自于噪点能量的均分，所以能量并没有增加。
计算当前像素在上一帧的对应点
我们要把上一帧的结果累加到这一帧，第一步首先得先找出我们要算的点在上一帧的对应点。
我们可以根据下面这个公式来找这一帧某个像素的对应点
首先把我们要算的当前帧的某个像素的世界位置通过乘Model的逆矩阵变回到局部空间，再乘上一帧的MVP矩阵，便可以知道，该像素点上一帧的屏幕位置了。
我们可以通过motion vector来存储各个物体或者场景上一帧的MVPl矩阵来分别计算上一帧的对应点。
由于作业里发生移动的只有小球本身，所以我们默认只需要一个motion vector记录小球上一帧的Model矩阵便可。
Matrix4x4 M =frameInfo.m_matrix[id];
auto world_position = frameInfo.m_position(x, y);
Matrix4x4 preWorldToScreen =
m_preFrameInfo.m_matrix[m_preFrameInfo.m_matrix.size() - 1];
我们可以直接从作业框架调取当前帧的Model矩阵、当前帧的WorldPosition已经Pi-1Vi1矩阵，所以我们要计算的只有上一帧的Model矩阵和这一帧Model的逆矩阵。
void Denoiser::Reprojection(const FrameInfo &frameInfo) {
int height = m_accColor.m_height;
int width = m_accColor.m_width;
Matrix4x4 preWorldToScreen =
m_preFrameInfo.m_matrix[m_preFrameInfo.m_matrix.size() - 1];
Matrix4x4 preWorldToCamera =
m_preFrameInfo.m_matrix[m_preFrameInfo.m_matrix.size() - 2];
#pragma omp parallel for
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
// TODO: Reproject
m_valid(x, y) = false;
m_misc(x, y) = Float3(0.f);
//判断id是否合理
int id = frameInfo.m_id(x, y);
if (id==-1) {
continue;
}
//已知Pi-1Vi-1也就是上面的preWorldToScreen
//我们只需要求Mi-1和Mi^(-1)即可
//Mi^(-1)
Matrix4x4 M_i_inv = Inverse(frameInfo.m_matrix[id]);
//Mi-1
Matrix4x4 M_pre_i = m_preFrameInfo.m_matrix[id];
auto world_position = frameInfo.m_position(x, y);
auto local_position = M_i_inv(world_position, Float3::EType::Point);
auto pre_world_postion = M_pre_i(local_position, Float3::EType::Point);
auto pre_screen_position =
preWorldToScreen(pre_world_postion, Float3::EType::Point);
//一定要加上这个判断，因为不确定该点对应上一帧的位置是否在屏幕外
if (pre_screen_position.x<0||pre_screen_position.x>=width||pre_screen_position.y<0||pre_screen_position.y>=height) {
continue;
}
else {
int pre_id = m_preFrameInfo.m_id(pre_screen_position.x,pre_screen_position.y);
if (id == pre_id) {
m_valid(x, y) = true;
//m_misc是临时颜色缓冲区、m_accColor是上一帧的颜色缓冲区
m_misc(x, y) =
m_accColor(pre_screen_position.x, pre_screen_position.y);
}
}
}
}
//需要交换缓冲区，不然下一帧需要同时对m_accColor进行读写，这可能会导致数据的不一致或竞争
std::swap(m_misc, m_accColor);
}
有几点需要注意的是：
我们需要在找上一帧对应位置的时候考虑几个问题
- 如果该像素点的内容在上一帧对应位置被遮挡了
- 如果该像素点的内容在上一帧在屏幕外
我们先看第一个问题， 如果该像素点的内容在上一帧对应位置被遮挡了，最终颜色有80%的贡献还来自上一帧的话，会出现残影
就比如下图，这一帧墙壁的位置出现了箱子的残影，因为墙壁没有动，他上一帧对应像素点的内容是移动过程中的箱子，我们仍取80%上一帧的贡献就会出现箱子的残影
在实时渲染里面，相比于噪声，其实残影给玩家带来的体验会更差，所以我们选择对像素点进行处理，如果上一帧的对应点和这一帧不是一个物体，那我们不再累加上一帧的颜色，上面代码直接通过id来判断，如果不是通过物体则跳过该循环。
对于第二个问题，我的措施是只有对应点在屏幕空间内才会进行颜色的累加，否则跳过该循环。
计算累加后的颜色
在实现上面Reprojection函数来找到上一帧对应点后，我们开始计算累加上一帧的当前帧的最终颜色，根据下面的流程来计算
我们来看看课程给出的公式
在判断这一帧和前一帧对应点都是同个物体的时候，需要把前一帧对应位置的颜色Clamp到这一帧来，简单来说，就是我们计算这一帧对应点附近一个范围的均值和方差，用此范围来限制上一帧对应点的颜色值，这是降低比如这一帧颜色突变了或者是上一帧有很亮的噪点这种问题的影响
void Denoiser::TemporalAccumulation(const Buffer2D<Float3> &curFilteredColor) {
int height = m_accColor.m_height;
int width = m_accColor.m_width;
int kernelRadius = 3;
#pragma omp parallel for
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
// TODO: Temporal clamp
Float3 color = m_accColor(x, y);
// TODO: Exponential moving average
float alpha = 1.0f;
if (m_valid(x,y)) {
alpha = m_alpha;
//根据采样范围，计算均值和标准差，目的是把上一帧颜色clamp到这一帧
int x_start = std::max(0, x - kernelRadius);
int x_end = std::min(width - 1, x + kernelRadius);
int y_start = std::max(0, y - kernelRadius);
int y_end = std::min(height - 1, y + kernelRadius);
Float3 mean_value(0.0f);
Float3 variance(0.0f);
for (int i=x_start;i<x_end;i++) {
for (int j=y_start;j<y_end;j++) {
mean_value += curFilteredColor(i, j);
variance += Sqr(curFilteredColor(i, j) - curFilteredColor(x, y));
}
}
int sum = (kernelRadius * 2 + 1) * (kernelRadius * 2 + 1);
//均值
mean_value /= (float)(sum);
//标准差
variance = SafeSqrt(variance/(float)(sum));
color = Clamp(color, mean_value - variance * m_colorBoxK,
mean_value + variance * m_colorBoxK);
}
m_misc(x, y) = Lerp(color, curFilteredColor(x, y), alpha);
}
}
std::swap(m_misc, m_accColor);
}
这里我们采样7×7范围的点进行均值和方差的计算，再用这个结果约束上一帧对应点的颜色。
在上面Reprojection和TemporalAccumulation代码里面其实有一个小细节不易于理解，就是框架里给出的是在代码末尾用swap来交换缓冲区，其实我们也可以不这么做。
//denoiser.cpp
Reprojection(const FrameInfo &frameInfo)
{
//..
//删除std::swap(m_misc, m_accColor);
}
TemporalAccumulation(const Buffer2D<Float3> &curFilteredColor)
{
//..
Float3 color = m_misc(x, y);
//..
m_accColor(x, y) = Lerp(color, curFilteredColor(x, y), alpha);
//..
//删除std::swap(m_misc, m_accColor);
}
其实可以这么理解，原框架中m_misc是临时颜色缓冲区，m_accColor是上一帧的颜色缓冲区，按我们的思想，用m_misc存储当前点需要累加的上一帧颜色，用color获取这个值，最后把当前帧累加后的最终结果更新到m_accColor作为下一帧的上一帧颜色缓冲，以此递归循环。
但是原框架用了swap来交换缓冲区，第一次交换缓冲区后m_misc存储的是上一帧的颜色，m_accColor则存当前点需要累加的上一帧颜色。因此在TemporalAccumulation里，color读取m_accColor的值其实也是在读取当前点需要累加的上一帧颜色，m_misc存的值是最终颜色。最后交换缓冲区后和我们写法结果是一样的。
交换缓冲区的目的其实是为了避免同时对m_accColor缓冲区的读写，同时读写可能会导致数据的不一致或竞争。
接下来我们来对比结果。
未经过Temporal的第19帧
经过Temporal的第19帧
上图中，经过Temporal的第19帧很明显会比未经过累加处理的第19帧效果要好，但出现了一个问题，就是有很明显的残。
这里残影出现的问题原因很明显，当前帧出现残影的点内容是墙壁，该点对应上一帧的点内容其实是球体， 按理来说我们判断过，如果不是同一物体颜色应该不会累加。
出现这个问题说明通过id判断这一帧在上一帧对应点是否为同一物体是不准确的。
我们来看看id到底是什么
//main.cpp
Buffer2D<float> id =
ReadFloatImage((inputDir / ("ID_" + std::to_string(idx) + ".exr")).str());
在main里存在上面代码，说明该场景id是由这个inputDir路径下的文件导入的，我们看看这个文件
我们打开这张图可以发现，图中颜色并没有区分球和墙体，也就是说导入的球和墙体的id是一样的
这就是导致了不管前一帧对应位置内容是球还是墙体，该值都会被累加到最终结果中，最终导致残影的出现。
用A-Trous Wavelet加速单帧降噪
作业中的单帧降噪中每个点都会进行一个64×64的采样操作，这无疑是巨费时间的，我们通过A-Trous Wavelet的方式来减小这个采样数量级。
A-Trous Wavelet采取的是多趟间隔采样的思想。该方法总共进行5次采样，每次都只采样一5×5的区域，每次采样间隔是2的i次方，i是采样趟数，参考下图
该方法进行采样总次数为555也就是125次采样，和64*64的4096次采样不在一个数量级，得到的效果却相差无几。
参照第一次Filter的结果，我们写出FIlter的ATrousWaveletFilter重载
Buffer2D<Float3> Denoiser::ATrousWaveletFilter(const FrameInfo &frameInfo) {
int height = frameInfo.m_beauty.m_height;
int width = frameInfo.m_beauty.m_width;
Buffer2D<Float3> filteredImage = CreateBuffer2D<Float3>(width, height);
int kernelRadius = 32;
#pragma omp parallel for
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
// TODO: Joint bilateral filter
// filteredImage(x, y) = frameInfo.m_beauty(x, y);
//
//定义Filter范围
//获取Filter中心店的数据
auto center_position = frameInfo.m_position(x, y);
auto center_normal = frameInfo.m_normal(x, y);
auto center_color = frameInfo.m_beauty(x, y);
Float3 result_color;
auto total_weight = 0.0f;
for (int pass= 0; pass < 5; pass++) {
for (int x_start = -2; x_start <= 2;x_start++) {
for (int y_start = -2; y_start <= 2;y_start++) {
int i = x + std::pow(2, pass) * x_start;
int j = y + std::pow(2, pass) * y_start;
auto current_position = frameInfo.m_position(i, j);
auto current_normal = frameInfo.m_normal(i, j);
auto current_color = frameInfo.m_beauty(i, j);
// Distance
auto D = Dot((center_position - current_position),
(center_position - current_position)) /
(2.0f * m_sigmaCoord * m_sigmaCoord);
// Color
auto C = Dot((center_color - current_color),
(center_color - current_color)) /
(2.0f * m_sigmaColor * m_sigmaColor);
// Normal
////SafeAcos和acos(x)
///一样，不过如果x超出范围的话，它会夹紧最近的可用值。
//和标准c函数acos(x)一致，这个值返回0到pi。
auto N = SafeAcos(Dot(center_normal, current_normal));
N *= N;
N / (2 * m_sigmaNormal * m_sigmaNormal);
// Plane
auto P =
Dot(center_normal, (current_position - center_position)) *
Dot(center_normal, (current_position - center_position)) /
(2 * m_sigmaPlane * m_sigmaPlane);
float weight = std::exp(-D - C - N - P);
total_weight += weight;
result_color += current_color * weight;
}
}
}
// filteredImage(x, y) = frameInfo.m_beauty(x, y);
// filteredImage(x, y) =(result_color.x / total_weight, result_color.y /
// total_weight,result_color.z / total_weight);
filteredImage(x, y) = result_color / total_weight;
}
}
return filteredImage;
}
最后来解释一下为什么可以跳着采样（沿用202和101课上的说法，没学过信号处理，不太清楚本质的原因）
首先我们要知道连个概念
- 越大的采样面积会溢出越多的低频信息
- 采样其实是在搬移频谱
蓝色部分是第一趟采样，去除了高频信息，而隔一定距离进行采样可以理解为移动了原理的频谱，使得频谱之间恰好不想交，所以刚好没有走样，可以参考下面这幅图理解
参考文章：
GAMES202-作业5：实时光线追踪降噪