蒙特卡洛路径跟踪 C++实现 图形学冬学期作业

蒙特卡洛路径跟踪

• 原理
• • 蒙特卡洛积分
  • 蒙特卡洛路径跟踪
  • 俄罗斯轮盘赌
  • 改写rendering equation
  • 改进的蒙特卡洛路径跟踪算法
  • 对每个像素发射光线
• 实现
• • 场景管理
  • BVH
  • • 光线和AABB求交（axis aligned bounding box）
    • 光线和三角形面片求交
    • 光线和基于BVH的场景求交
  • PathTracer
  • • 相机配置与光线生成
    • 光线跟踪
    • 随机采样
    • 处理折射
• 结果

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的构建步骤为：

1. 根据空间morton码将物体在空间中排序
   sort(scene.f.begin(), scene.f.end(), compare);
2. 计算和存储物体的最小包围盒，只在叶节点上存储物体信息
3. 自底向上构建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;
	}
	if (tymin > tymax) {
		double temp = tymin;
		tymin = tymax, tymax = temp;
	}
	if (tzmin > tzmax) {
		double temp = tzmin;
		tzmin = tzmax, tzmax = temp;
	}

	if (txmax < 0 || tymax < 0 || tzmax < 0) return false;
	if (txmin <= 0 && tymin <= 0 && tzmin <= 0) return true;
	if (dmax(txmin, tymin, tzmin) <= dmin(txmax, tymax, tzmax)) return true;
	else return false;
}

光线和三角形面片求交

先计算光线和三角形面片所在平面的交点，再判断该交点是否在三角形内部，即将问题分解为两步。

1. 对于光线和三角形面片所在平面的交点，设交点为P，则P点既在光线上，也在平面上，由此得到平面方程n·(v1-P) = 0,将射线方程代入平面方程中，解得 t = ((triangle.v1 - ray.start) * norm) / (norm * ray.direction); 将t代入光线方程得到光线与平面的交点。
   

   
   

2. 判断该交点是否在三角形内部，连接该交点和三角形的三个顶点得到三个向量AP,BP,CP，三角形三条边的向量为AB,BC,CA，若点在三角形内部，则边与交点与三角形顶点的连线向量的叉乘一定为同号，即AB×AP与BC×BP与CA×CP得到的向量与法向量的点乘为同号。
   

   
   

代码如下：
intersect函数判断ray和三角形面片是否存在交点，若存在交点，将交点存储至ret中并返回true

bool intersect(Ray &ray, Face &triangle, Vertex &ret)
{
	double t;
	Vertex norm = triangle.norm;
	t = ((triangle.v1 - ray.start) * norm) / (norm * ray.direction);

	Vertex p = ray.start + (ray.direction * t);
	Vertex ap = p - triangle.v1, bp = p - triangle.v2, cp = p - triangle.v3;
	Vertex ab = triangle.v2 - triangle.v1, bc = triangle.v3 - triangle.v2, ca = triangle.v1 - triangle.v3;
	Vertex cross1 = ab.cross(ap), cross2 = bc.cross(bp), cross3 = ca.cross(cp);
	double dir1 = cross1 * norm, dir2 = cross2 * norm, dir3 = cross3 * norm; //叉乘法线
	double j1, j2, j3;
	j1 = dir1 * dir2, j2 = dir1 * dir3, j3 = dir2 * dir3;
	bool judge = false;
	if (j1 >= 0 && j2 >= 0 && j3 >= 0) judge = true; //叉乘符号相同
	ret = p;
	return judge;
}

光线和基于BVH的场景求交

光线和BVH求交的步骤如下：

1. 光线和BVH root求交，如果相交，分别对其左右子树进行相交测试，如果不相交，就终止
2. 对左右子树的相交测试如同步骤1一样进行
3. 当进行到叶子节点时，如果光线和叶子节点的包围盒相交，将光线与叶子节点中的物体做相交测试

void PathTracer::bvh_intersect(Ray ray, BVH &bvh, intersection &v, int current_node, int current_level, bool &flag)
{
	int index = bvh.findIndex(current_node, current_level); //找到当前节点对应的bvh的物理内存号
	if (bvh.bvh[index].level == bvh.Level && intersect(ray,bvh.bvh[index].b)) { //如果为叶节点
		Vertex ret;
		if (intersect(ray, *bvh.bvh[index].object, ret)) { //将光线与三角形面片求交
			intersection i;
			i.p = ret;
			i.ray = ray;
			i.t = (i.p.x - ray.start.x) / ray.direction.x;
			i.f = *bvh.bvh[index].object;

			Vertex garcov = findGarCor(i.f, i.p); //得到重心插值坐标
			i.pn = (i.f.vn1 * garcov.x) + (i.f.vn2 * garcov.y) + (i.f.vn3 * garcov.z);
	
			if (flag == false) {
				if (i.t > 0) {
					v = i;
					flag = true;
				}
			}
			else if (i.t > 0 && i.t < v.t) v = i; //保留最小正值t，确保交点为最近交点
		}

		return;
	}
	else { //如果不是叶节点，分别对bvh的左右子树做相交测试
		if (intersect(ray, bvh.bvh[index].b)) {
			int next_node1 = 2 * current_node + 1;
			int next_node2 = 2 * current_node + 2;
			bvh_intersect(ray, bvh, v, next_node1, current_level+1,flag);
			bvh_intersect(ray, bvh, v, next_node2, current_level+1,flag);
		}
	}
	
}

PathTracer

在pathtracer中，我们定义了几个重要的函数：

相机配置与光线生成

void generateImg(scene_data &scene, BVH &bvh, image &img, int N_ray_per_pixel);

该函数向每个像素发射N根光线，并与场景求交得到radiance。

从.camera文件中读取的相机参数有eye（相机在世界坐标系下的位置），look_at（相机在世界坐标系下看向的位置），up（相机的上方向），fovy（相机的视角 field of view），屏幕的分辨率。

由于需要进行从视点到屏幕的每个像素发射光线的 *** 作，定义屏幕空间在世界坐标系下的x、y单位向量。

Y方向的单位向量即为up归一化，由于相机看向的方向为向量look_at – eye，而屏幕空间的X方向与Y方向核相机看向的Z方向垂直，所以X方向的单位向量等于Y叉乘Z向量，得到如下代码

	scene.camera.up.normalize();
	Vertex dir = scene.camera.look_at - scene.camera.eye;
	double l = dir.norm();
	double dy = tan(scene.camera.fovy/2/180 * pi) * l; //根据fovy求世界坐标系中半屏幕的x,y大小
	double dx = dy / scene.camera.height * scene.camera.width;

	Vertex screen_center = scene.camera.look_at; //屏幕在世界坐标系下中心点的坐标
	double pdx = 2 * dx / scene.camera.width, pdy = 2 * dy / scene.camera.height; //每个像素在世界坐标系下的大小
	
	Vertex screen_x_dir = dir.cross(scene.camera.up);
	screen_x_dir.normalize();
	Vertex screen_y_dir = scene.camera.up;

	Vertex screen_pdy = screen_y_dir * pdy;
	Vertex screen_pdx = screen_x_dir * pdx;

	Vertex start_point = screen_center - (screen_x_dir * dx) + (scene.camera.up * dy); //逐像素发射光线 像素的中心位置

从像素（0，0）到像素（width，height）分别打出N根光线，并对每根光线应用光线跟踪算法，得到当前像素的辐射度值（radiance）

光线跟踪

Vertex shade(intersection &p, Vertex dir, scene_data &data, BVH &bvh);

shade函数接收四个参数，P为交点数据，包含交点坐标、交点所在的平面、平面的法向量和材质信息，dir为从交点p到视线的方向，data和bvh分别存储了场景信息和BVH信息。

该算法分别计算光源对点的直接光照和来自其他半球面的间接光照，再将两个积分叠加。

随机采样

对于直接光照部分：
按照点对于光源面计算，随机采样光源点。

先根据光源三角形面片的面积大小以不同概率随机选择三角形面片进行采样，实现面积越大选择的概率越大。

static std::default_random_engine e; //c++11随机数生成引擎
static std::uniform_real_distribution<double>u1(0, total_aera); //均匀采样
double rnd = u1(e);

按照离散采样求逆的方法得到选择的三角形面片，再根据三角形重心坐标随机采样，即P = aA+bB+cC
随机均匀采样a,b,c并且归一化得到三角形面片中一个点以及该点的法向量

static std::uniform_real_distribution<double>u2(0, 1);
double rnd1 = u2(e), rnd2 = u2(e), rnd3 = u2(e);
double p1 = rnd1 / (rnd1 + rnd2 + rnd3), p2 = rnd2 / (rnd1 + rnd2 + rnd3), p3 = rnd3 / (rnd1 + rnd2 + rnd3); //归一化
xl = sample_face.v1 * p1 + sample_face.v2 * p2 + sample_face.v3 * p3; //采样点
vn = sample_face.vn1 * p1 + sample_face.vn2 * p2 + sample_face.vn3 * p3; //法向量

得到采样的光源点后，根据计算蒙特卡洛估计子，f(x)/p(x)。

并且需要判断着色点到光源到中间是否存在遮挡，即打出一根shadow ray。

		double visibility = 1;
		Ray rl;
		rl.start = p.p + (direction * 0.01);//add micro turbulance
		rl.direction = direction;
		intersection inter;

		ray_intersect(rl, data, bvh, inter);
		if (inter.f.material != sample_face.material) { //如果打到的第一个物体不是光源，即不可见
			visibility = 0;
		}

计算着色点的radiance

double pdf_light = double(1)/total_aera; //pdf_light = 1/A where A is the aera of light source
double cos_theta = abs(direction * vn / direction.norm() / vn.norm());
double cos_theta_hat = abs(direction * p.pn / direction.norm() / p.pn.norm());
Vertex intensity = l.radiance  *  cos_theta * cos_theta_hat / pow(distance(xl, p.p), 2) / pdf_light * visibility;
double kd_dots = direction * p.pn; //cos between light to intersection and face normal	
//only add diffuse
if (kd_dots > 0) {
	L_dir.x += kd.x * intensity.x * kd_dots / pi;
	L_dir.y += kd.y * intensity.y * kd_dots / pi;
	L_dir.z += kd.z * intensity.z * kd_dots / pi;
}

由于反射项需要按照BRDF重要性采样，在采样直接光照时反射项的贡献暂时先忽略不计，全部放到间接光照中求。

在间接光照中求漫反射的贡献时是需要判断打到的物体是否是光源的，如果是光源，那么这根光线的贡献是需要舍弃的，而间接光照求镜面反射的贡献时，如果打到光源，这根光线的贡献仍然有效，这保证了在计算直接光照时不计算镜面反射项，而放到后续一起算，结果也是正确的。

对于间接光照部分：
对于需要计算间接光照的情况，按照BRDF重要性采样。

根据前面我们介绍的，可以用任意一个概率密度函数采样，得到蒙特卡洛估计子，这个蒙特卡洛估计子只要除以采样的概率密度函数就可以了。

而蒙特卡洛估计子的分子是有BRDF函数的，如果我们能够用BRDF作为概率密度函数进行采样，那么就可以直接将分子分母约分，去掉BRDF项，有了这种想法，我们有了以下的代码思路。

1. 先进行俄罗斯轮盘赌判断，如果需要进行采样，则进行步骤2，否则返回0
2. 根据BRDF采样下一根光线（需要得到光线的方向即可）
   我们根据物体的材质信息，按照不同的概率进行漫反射和镜面反射。
   

   
   

   概率为材质的ks和kd的模的比值

double kd_norm = kd.norm(), ks_norm = m->ks.norm();
	type ray_type;
	if (ks_norm != 0 && kd_norm / ks_norm < u2(e)) {
		//sample specular
		Vertex reflect; //reflect ray direction
		Vertex incoming = dir.negative();
		reflect = incoming - p.pn * (incoming * p.pn) * 2;
		direction = BRDFImportanceSampling(reflect, SPECULAR,p.f, m->Ns);
		ray_type = SPECULAR;
	}
	else {
		//sample diffuse
		direction = BRDFImportanceSampling(p.pn, DIFFUSE, p.f, m->Ns);
		ray_type = DIFFUSE;
	}

	Ray r(p.p + (direction * 0.01), direction, ray_type);
	return r;

3. 根据不同的反射类型计算下一根光线的方向

	static std::default_random_engine e(time(NULL));
	static std::uniform_real_distribution<double>u1(0, 1);

	double phi = u1(e) * 2 * pi;
	double theta;
	if (type == DIFFUSE) { //diffuse
		// for cosine-weighted Lambertian
		theta = asin(sqrt(u1(e)));
	}
	else if (type == SPECULAR) {
		// for sampling specular term
		theta = acos(pow(u1(e), (double)1 / (Ns + 1)));
	}

采样的theta和phi为球坐标系下的，需要转换成直接坐标系，并且将着色点的法向作为其中一个正方向，转换的公式为：

Vertex sample(sin(theta)*cos(phi), cos(theta), sin(theta)*sin(phi));

4. 得到下一根光线后将下一根光线与场景求交，得到返回的radiance

if (russian_Roulette(P_RR)) {
		Ray r = nextRay(p, dir, kd, data); //sample next ray based on BRDF
		//find the first object ray hits
		intersection ret;
		if (ray_intersect(r, data, bvh, ret)) {
			Vertex intensity = shade(ret, r.direction.negative(),data,bvh) / P_RR;
			if (r.ray_type == DIFFUSE) {
				if (data.light_map.find(ret.f.material) == data.light_map.end()) { //hit none emitting object
					L_indir.x += kd.x * intensity.x;
					L_indir.y += kd.y * intensity.y;
					L_indir.z += kd.z * intensity.z;
				}
			}
			else if (r.ray_type == SPECULAR) {
				L_indir.x += m->ks.x *  intensity.x;
				L_indir.y += m->ks.y *  intensity.y;
				L_indir.z += m->ks.z *  intensity.z;
			}
			else {
				L_indir = L_indir + intensity;
			}
		}
	}

处理折射

光线的折射方向只有一个，因此需要对折射物体做特殊的处理，即不需要进行随机采样，下一根光线必定是上一根光线折射的方向。

下面的函数判断该光线在p点能否产生折射（是否会发生全反射）

bool PathTracer::Refract(Vertex dir, Vertex &normal, double eta, Vertex &refract_dir) {
	// Ref: https://www.cnblogs.com/night-ride-depart/p/7429618.html
	Vertex &i = dir;
	Vertex &n = normal;
	float cosi = i * n;
	float cost2 = 1.0f - eta * eta * (1.0f - cosi * cosi);
	if (cost2 >= 0.0f) {
		refract_dir = i * eta - n * (eta * cosi + std::sqrt(cost2));
		return true;
	}
	else {
		// total internal reflection
		return false;
	}
}

对于折射的物体，产生下一根光线的代码如下：先判断是折射进入物体还是折射出物体，然后根据物体的折射率和菲尼尔公式近似出折射的方向，如果不能发射折射，就按照全反射处理。

	//for fraction
	if (m->Ni > 1) {
		double n1, n2;
		double cos_in = dir.negative() * p.pn;
		Vertex normal;
		if (cos_in > 0) {
			// out of glass
			normal = p.pn.negative();
			n1 = m->Ni;
			n2 = 1.0;
		}
		else {
			// in to the glass
			normal = p.pn;
			n1 = 1.0;
			n2 = m->Ni;
		}

		// Ref : https://en.wikipedia.org/wiki/Schlick%27s_approximation
		double rf0 = pow((n1 - n2) / (n1 + n2),2);
		double fresnel = rf0 + (1.0f - rf0) * pow(1.0f - std::abs(cos_in),5);
		if (fresnel < u2(e) ){
			if (Refract(dir.negative(), normal, n1 / n2, direction)) {
				return Ray(p.p, direction, TRANSMISSION);
			}
			else {
				// only specular for refraction material
				Vertex reflect; //reflect ray direction
				Vertex incoming = dir.negative();
				reflect = incoming - normal * (incoming * normal) * 2;
				return Ray(p.p, reflect, SPECULAR);
			}

		}
	}

结果

至此，我们已经实现了一个基本的蒙特卡洛路径跟踪算法，可以对场景进行渲染。

渲染时每个像素的采样数会影响到最终结果的质量和渲染的时间，下面展示在不同采样数下不同场景
Cornell-box-spp25

Cornell-box-spp100

bedroom-spp10

bedroom-spp25

bedroom-spp100

bedroom-spp256

veach-mis-spp25

veach-mis-spp100