蒙特卡洛路径跟踪详解：原理、改进与实战应用-CSDN博客

本文链接：https://blog.csdn.net/Listoree/article/details/124081645

蒙特卡洛路径跟踪

原理
实现
结果

github代码链接: https://github.com/Arieys/MonteCarloPathTracing.

原理

蒙特卡洛路径跟踪的基本目标为求解rendering equation。而rendering equation中最难求解的积分项，可以通过蒙特卡洛积分方法求解。下面我们给出在半球面积分的渲染方程：
在这里插入图片描述
其中，Lo为x点到w0方向的出射radiance，Le为x点到w0的自发光radiance，Lr为x到w0的反射radiance，由于物体的发光属性一般是已知的，我们默认Le为已知，所以求解rendering equation即为求解Lr中的积分项。该积分项可以使用蒙特卡洛积分的方法进行求解。

蒙特卡洛积分

蒙特卡洛积分是通过统计和采样的方式求解一个定积分，假设我们要求解如下一个定积分I
在这里插入图片描述
给定一个均匀采样概率密度函数p(x) = 1，可以得到该积分的蒙特卡洛估计子F

理论证明，该蒙特卡洛估计子的期望最终会收敛到积分值，即我们可以通过采样的方式计算出一个定积分，而该定积分最终能够收敛到真值，并且采样数越多，方差越小，且选择的采样概率密度函数也会影响积分的收敛速度，因为在实际中，我们需要尽可能地选用比较好的概率密度函数。

蒙特卡洛路径跟踪

在蒙特卡洛路径跟踪里，我们的最终目的就是求解渲染方程，计算出rendering equation中积分的蒙特卡洛估计值，并且尽量得到较小方差的结果。通过蒙特卡洛积分的方法求解渲染方程自然要涉及到采样，因为蒙特卡洛积分方法的原理就是采样，求均值。而不同的蒙特卡洛路径积分算法即是对不同采样方法的研究，采样方法不同，得到的结果的效果也会呈现比较大的差异。
我们需要明确的一点是，蒙特卡洛路径跟踪，目的是为了解决渲染方程，即求解从点x到方向w0的radiance。明确了这一点，我们开始思考求解渲染方程需要哪几步。首先，我们需要一个概率密度函数来对光线进行采样，其次，我们需要用这个概率密度函数来对光线进行实际的采样，最后，根据f(x）和p(x)得到蒙特卡洛估计子的值，并且我们将这种随即采用重复N次，在实际的应用中，N越大，得到的效果越可能收敛，需要的渲染时间也越长。到这里，我们已经可以写出一个蒙特卡洛路径跟踪的基础版本了。
(图片来自games101课件）
在这里插入图片描述
这里的p和w0即为我们求解渲染方程的p点和方向w0，随机在p点的法向半球面采样一根光线，计算来自该跟光线的radiance，将它代入渲染方程积分项里的函数f(x)，再除以采样概率密度函数p(x)得到蒙特卡洛估计子。在这里，由于是在半球面均匀采样，概率密度函数就是半球面的面积2pi。
然而，基础版本的蒙特卡洛路径跟踪算法存在两个问题，第一个问题为算法里面嵌套了递归，即shade函数调用了自身，而递归函数需要递归出口，该函数存在永远无法停下的问题。第二个问题为一般场景中的光源面积不会特别大，而要通过在半球面随机采样，可能采样的光线很难打到光源上，造成了很多光线对点p实际上并没有贡献，导致收敛的速度很缓慢，可能需要采集几千几万跟光线才能收敛到正确结果。
接下来我们分别解决这两个问题，对于第一个问题，我们采用俄罗斯轮盘赌的方法，对于第二个问题，我们采用分别计算直接光照和间接光照的方法。

俄罗斯轮盘赌

在之前的算法里，我们在着色点p上，会一直发射一根光线，根据光线返回的radiance计算蒙特卡洛估计子的值。假设手动设置一个概率p0，在着色点处，我们以概率p0打出一根光线，像之前一样，以概率1-p0停止，如此以来，便有了递归出口。即由1-p0的概率会停下，即使1-p0的概率再小，只要不是0，就会有停下的时候。并且，我们在计算着色点的radiance贡献时除以一个p0，那么在p0的概率下我们得到原有的radiance除以p0，在1-p0的概率下我们得到0，总体的数学期望即为

p0 * radiance / p0 + （1-p0）* 0 = radiance

即我们这么做了之后，在着色点p处得到的radiance的期望依然为原有的radiance，并且现在的算法有了递归出口。

改写rendering equation

在这里，我们可以将积分的求解分为两部分，一部分为直接光照，另一部分为间接光照。即x点到w0的出射radiance由两部分组成，一部分是来自光源的直接光照，另一部分是来自其他物体的间接光照。
对于直接光照的部分，我们可以先将渲染方程改写。最开始的渲染方程的积分域是在半球面积分，我们可以将其改写为对光源的面积分。
在这里插入图片描述
其中，x’为出射方向，即改写后的渲染方程计算的是x点对于x’点的出射radiance，x’‘点为对面的采样点，Fr为BRDF函数的点形式，其本质仍为BRDF函数，G为改变积分域产生的几何项
（图片来自浙大图形学课件）
在这里插入图片描述
V为可见性项，即x和x’‘之间是否存在遮挡，G的分母为x到x’'的距离平方，分子中的cos分别为入射光线和点法线的夹角以及入射光线和对面光源采样点处夹角的cos值，积分中的L为x‘’到x的radiance。由此，我们得到了两种渲染方程的形式，第一种的积分域为半球面，而第二种的积分域则为一个面。
改写rendering equation后，在着色点p的积分值可以分为两部分求解，一部分是该积分来自光源的贡献，另一部分为该积分来自其他非光源物体的贡献。对于来自光源的贡献部分，我们通过对光源面积分，得到radiance，对于来自非光源物体的部分，我们仍旧通过原有的在半球面采样的方式得到radiance。

改进的蒙特卡洛路径跟踪算法

由此，我们已经将前面的两个问题都得到了解决，已经可以得到一个基础的蒙特卡洛路径跟踪算法的改进版本了。
（图片来自games101课件）
在这里插入图片描述
到目前为止，我们的算法已经可以得到任意一点着色点p对于w0方向的radiance了。

对每个像素发射光线

有了着色点的radiance，我们可以得到一个蒙特卡洛路径跟踪的简单版本了。将场景中的相机和每个像素的位置连接，对于每个不同的像素，相机和像素中心点连接成为一根光线，发射到场景中，并且计算打到场景中的点p的radiance。发射N根光线，将每次的radiance相加求平均，即为将蒙特卡洛估计子求平均，根据蒙特卡洛积分方法，最终收敛到积分的值。由此我们得到了每个像素的颜色值（radiance），整个过程如下图所示
在这里插入图片描述
对于每个像素，也可以去像素中的任意一点发射光线，即对于同一个像素，每次发射的光线也可以是不同的，如下图所示（图片来自real time rendering）

实现

github代码链接: https://github.com/Arieys/MonteCarloPathTracing.

有了原理上的支持，下面我们就用C++实现一个简单的蒙特卡洛路径跟踪算法，我将代码分为三个模块，分别为场景管理、PathTracer、BVH。

场景管理

在场景管理模块，我们将读入场景的几何信息（obj文件），纹理和材质信息（mtl文件）以及相机信息（camera文件），我们定义各种类存储各种不同的信息。

class VertexT {
   //贴图坐标
public:
	double vtx, vty;
};

class Face {
   //面对应的顶点索引
public:
	Vertex v1, v2, v3;
	Vertex vn1, vn2, vn3;
	VertexT vt1, vt2, vt3;
	Vertex norm;
	std::string material;
	unsigned int morton_code;
};

class Material {
   //材质信息
public:
	std::string name;
	Vertex kd, ks;
	double Ns;
	double Ni;
	vector <Face> f; //存放所有对应该材质的三角形面片
	std::string map_Kd_filename;
	bool mapping_flag;
	cv::Mat img;
	int map_height, map_width;
};

class Camera {
   //场景相机信息
public:
	Vertex eye, look_at, up;
	double fovy;
	int width, height;
};

class Light {
   //光源信息
public:
	std::string name;
	Vertex radiance;
};

class Ray {
   //光线类
public:
	Vertex start;
	Vertex direction;
	enum type ray_type;
};

BVH

在光线跟踪的过程中，需要实现光线和场景中物体的求交，由此涉及到大量的求交工作。而在较大的场景中，采用遍历的方法，将光线与每个物体求交，再排序找出一个最近的物体为最近交点，这种方法在时间上几乎是不可行的，因此我们需要用到加速结构。在这里，我实现的是最常用的加速结构BVH中比较新的Binary Ostensibly-Implicit Trees，根据Binary Ostensibly-Implicit Trees for Fast Collision Detection实现隐式BVH的构建和使用，BVH的构建步骤为：

根据空间morton码将物体在空间中排序
sort(scene.f.begin(), scene.f.end(), compare);
计算和存储物体的最小包围盒，只在叶节点上存储物体信息
自底向上构建BVH

关于BVH实现的细节，可以参考文献《Binary Ostensibly-Implicit Trees for Fast Collision Detection》
有了BVH之后，我们就可以将光线快速与场景求交。基本的逻辑为，如果光线和BVH中的包围盒没有交点，那么它和包围盒中的任何一个物体或包围盒都没有交点。下面先介绍光线和AABB以及光线和三角形求交的算法。

光线和AABB求交（axis aligned bounding box）

将光线表示为Light = Start_position + t * direction的形式
应用SLABS算法，分别计算光线和XYZ平面的交点的坐标，得到光线进入和离开XYZ轴的时间t，若光线进入XYZ轴的最晚时间t大于光线离开XYZ轴的最早时间t，则说明在光线完全进入包围盒前，某一个轴光线已经离开，即判断光线与AABB没有交点，如下图所示
在这里插入图片描述
以下为光线和AABB有交点的情况

代码如下：

bool intersect(Ray &ray, boundingBox &b)
{
   
	double txmin, txmax, tymin, tymax, tzmin, tzmax;
	txmin = (b.min_x - ray.start.x) / ray.direction.x;
	txmax = (b.max_x - ray.start.x) / ray.direction.x;
	tymin = (b.min_y - ray.start.y) / ray.direction.y;
	tymax = (b.max_y - ray.start.y) / ray.direction.y;
	tzmin = (b.min_z - ray.start.z) / ray.direction.z;
	tzmax = (b.max_z - ray.start.z) / ray.direction.z;
	if (txmin > txmax) {
   
		double temp = txmin;
		txmin = txmax, txmax = temp;
	}