文章目录

前言
蒙特卡洛路径追踪
- 蒙特卡洛积分
- 路径追踪
- - Whitted-Style光线追踪的局限性
  - 渲染方程
  - - 直接光照
    - 间接光照
    - 两个问题
    - 优化
  - 总结
Assignment 7
- 准备工作
- 实现过程
- 最终结果
- 提高部分
参考

前言

本篇博客作为Games101关于Ray Tracing部分以及Assignment7路径追踪的复习总结。

蒙特卡洛路径追踪

蒙特卡洛积分

在此之前先复习一下概率密度函数(PDF, Probability Density Functions)，可以理解为取某一随机变量x的概率是多少的函数p(x)p(x)p(x)

PDF具有如下特点

$ p(x) >= 0 \ and \int{p(x)}=1$
E[X]=∫xp(x)dxE[X]=\int{xp(x)}dxE[X]=∫xp(x)dx

对于复杂函数f(x)f(x)f(x)求定积分，用牛顿莱布尼茨公式十分困难，这时候需要用到另一种积分方法，也就是下面的蒙特卡洛积分(Monte Carlo Integration)
∫f(x)dx=1N∑i=1Nf(Xi)p(Xi)Xi∼p(x)\int{f(x)}dx = \frac{1}{N}\sum_{i=1}^{N}{\frac{f(X_i)}{p(X_i)}}\ \ \ \ \ \ \ X_i \sim p(x) ∫f(x)dx=N1i=1∑Np(Xi)f(Xi) Xi∼p(x)

这里的做法可以理解为按p(x)p(x)p(x)的概率取随机变量XiX_iXi，在图形学里也就是随机采样函数值再累加做平均。

一些需要注意的地方：

采样点N越多，方差越小
在x上采样，也是在x上求积分

路径追踪

路径追踪(Path Tracing)根据渲染方程实现全局光照，或者说基于物理的渲染。个人理解叫追踪的原因是因为光线本来是打到物体上再多次弹射到人眼所以人才能看见东西，追踪就是逆着这个路径，从眼睛“射出”光线，“所见即所得”。

Whitted-Style光线追踪的局限性

Whitted-style ray tracing是一种理想化模型，存在一些问题，光线只有反射或折射，在漫反射材质停止弹射，没有漫反射，因此Whitted-style ray tracing渲染出来的效果只有光源和直接光照，如下图左侧所示。我们理应看到的效果是下图右侧这种，就像是照片拍出来的一样，需要实现这样真实的效果就需要用到路径追踪。

渲染方程

根据辐射度量学，我们一通操作最终求得渲染方程，即在物理意义上得出的正确的能够描述我们现实生活所见的方程。
Lo(p,ωo)=Le(p,ωo)+∫ΩLi(p,ωi)fr(p,ωi,ωo),(n⋅ωi)dωiL_o(p,\omega_o)=L_e(p,\omega_o)+\int_{\Omega}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o),(n\cdot\omega_i)d\omega_i Lo(p,ωo)=Le(p,ωo)+∫ΩLi(p,ωi)fr(p,ωi,ωo),(n⋅ωi)dωi

某着色点到相机的radiance由其本身发光产生的部分，和从所有立体角接受到的radiance并反射到相机这个一定角度的部分组成。

为了方便计算和统一，这里将着色点本身产生的Le(p,ωo)L_e(p,\omega_o)Le(p,ωo)省略掉，我们更关注如下这个渲染方程（或者说反射方程）。
Lo(p,ωo)=∫ΩLi(p,ωi)fr(p,ωi,ωo),(n⋅ωi)dωiL_o(p,\omega_o)=\int_{\Omega}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o),(n\cdot\omega_i)d\omega_i Lo(p,ωo)=∫ΩLi(p,ωi)fr(p,ωi,ωo),(n⋅ωi)dωi
根据蒙特卡洛积分可以得到
Lo(p,ωo)≈1N∑i=1NLi(p,ωi)fr(p,ωi,ωo)(n⋅ωi)p(ωi)L_o(p,\omega_o)\approx\frac{1}{N}\sum_{i=1}^{N}\frac{L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)}{p(\omega_i)} Lo(p,ωo)≈N1i=1∑Np(ωi)Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)
这个方程可以理解为按p(ωi)p(\omega_i)p(ωi)的概率随机选取N个方向ωi\omega_iωi采样，再累加求平均。

注意，这里的ωo\omega_oωo和ωi\omega_iωi都是由着色点向外的

直接光照

仅在直接光照下，考虑计算radiance的场景如下

由此可以得到着色方程的shade伪代码

shade(p, wo)Randomly choose N directions wi~pdfLo = 0.0For each wiTrace a ray r(p, wi)If ray r hit the light    // 弹射的光线必须要能打到光源，这个过程中不能被其他物体阻挡Lo += (1 / N) * L_i * f_r * cosine / pdf(wi)Return Lo

其中p是着色点，wo是入射光线，wi是出射光线，r是由眼睛穿过需要渲染的像素到实际着色点的光线。

间接光照

如果路径追踪经过一次弹射后没有打到光源而是打到物体，就构成了间接光照。我们考虑如下这个场景

P接收到的radiance是由Q反射光源得到的，因此为了求点P反射的光，我们需要先求点Q反射出的光。显然，这里用到了递归的思路，将上述shade伪代码稍作修改即可完成

shade(p, wo)Randomly choose N directions wi~pdfLo = 0.0For each wiTrace a ray r(p, wi)If ray r hit the lightLo += (1 / N) * L_i * f_r * cosine / pdf(wi)Else If ray r hit an object at qLo += (1 / N) * shade(q, -wi) * f_r * cosine / pdf(wi)Return Lo

可以看到只是增加了一个判断条件，判断从点p反射出的的光线wi是射到了光源还是其他物体，如果是光源直接加，如果是其他物体则递归先求L_i，也就是shade(q,-wi)。
注意，由于shade的参数要求光线方向向外，所以对于点Q，需要将wi取反。

两个问题

上述shade伪代码还存在两个问题

每次都会采样N个方向，最终导致指数爆炸
递归没有出口

对于第一个问题，可以用下图形象表示

在计算间接光照时，由于每次都要采样N个点，最终导致NxN^xNx的指数增长。但有一个数是可以避免这种情况，那就是N取1。

在采样点每次采样只取1个方向，我们的shade伪代码又有了如下改动

shade(p, wo)Randomly choose ONE direction wi~pdf(w)Trace a ray r(p, wi)If ray r hit the lightReturn L_i * f_r * cosine / pdf(wi)Else If ray r hit an object at qReturn shade(q, -wi) * f_r * cosine / pdf(wi)

由前文的蒙特卡洛积分我们知道，N取越大越等接近原定积分的值，N取越小结果的噪声或者说偏差就会越大。我们这里取N为1，可想而知会有很大偏差。为了解决这个问题，可以采取对同一像素多次采样的方法，根据需要设置对每个像素采样的次数（spp, samples per pixel）。实现的效果如下图所示，可以看到间接光照每次采样的弹射结果虽然都不一样，但多次采样求平均后就可以有效减小噪声。

上述的采样方法用伪代码描述如下，与shade伪代码类似

ray_generation(camPos, pixel)Uniformly choose N sample positions within the pixelpixel_radiance = 0.0For each sample in the pixelShoot a ray r(camPos, cam_to_sample)If ray r hit the scene at ppixel_radiance += 1 / N * shade(p, sample_to_cam)Return pixel_radiance

现在我们来解决第二个问题，这个递归算法没有出口所以会一直持续下去，因此我们需要认为设置一个出口。

这里采用的方法是俄罗斯轮盘赌。设定一个值Russian Roulette (RR)，光线每次弹射前需先与RR“赌博”，如果成功则继续弹射，否则停止。我们可以设在概率P_RR下，将射出光线返回的渲染方程结果Lo除以P_RR，那么将会有1-P_RR的概率停止弹射直接返回0。在此假设下，我们可以得到期望值E，是与我们期望得到的Lo相同的。
E=P∗(Lo/P)+(1−P)∗0=LoE=P*(Lo/P)+(1-P)*0=Lo E=P∗(Lo/P)+(1−P)∗0=Lo
由此，我们的shade伪代码再次改动如下

shade(p, wo)Manually specify a probability P_RRRandomly select ksi in a uniform dist. in [0, 1]If (ksi > P_RR) return 0.0;Randomly choose ONE direction wi~pdf(w)Trace a ray r(p, wi)If ray r hit the lightReturn L_i * f_r * cosine / pdf(wi) / P_RRElse If ray r hit an object at qReturn shade(q, -wi) * f_r * cosine / pdf(wi) / P_RR

优化

到此，我们终于得到了一个正确的路径追踪算法，的的确确能正确渲染出一个像素，但是目前这种做法的效率其实并不高。在计算直接光照时，如果我们按照上面得伪代码计算直接光照部分，那么每次只会随机射出一根光线，实际上很难打到光源上，计算效率低下。

为了提高成功采样到直接光照的概率，我们可以对直接光照计算的方法进一步修改，在弹射点进行多次采样。与间接光照递归计算相比，因为直接光照只会弹射一次到光源，所以这样并不会造成指数爆炸。

如上图所示，随着着色点采样数的提高，成功采样的几率自然会提高。但是当光源很小时，如果我们继续采用平均采样，那么弹射的数跟光线中可能只有少数几根能打到光源上，其他采样光线可以说都是被浪费掉了，这同样是我们不希望看到的。

如果能对更加确定的东西采样，那么这个问题自然迎刃而解。考虑我们是想要找到弹射点能打到光源的有效光线，所以反过来考虑，我们对光源进行采样。并且实际上，蒙特卡洛积分也是允许我们自定义采样方法的。
Lo(p,ωo)=∫ΩLi(p,ωi)fr(p,ωi,ωo)(n∗ωi)dωiL_o(p,\omega_o)=\int_{\Omega}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n*\omega_i)d\omega_i Lo(p,ωo)=∫ΩLi(p,ωi)fr(p,ωi,ωo)(n∗ωi)dωi
上面是我们在前文我们得到的渲染方程，可以看到我们的积分目标是dωid\omega_idωi，可以解释为我们是从着色点出发考虑这个积分的。如果我们对光源进行采样，积分对象自然也需要变化，所以我们需要找到dωid\omega_idωi和dAdAdA之前的关系。下图很好展现了dωid\omega_idωi和dAdAdA之前的关系。

回顾以前对立体角ω\omegaω的定义
dω=dAr2d\omega = \frac{dA}{r^2} dω=r2dA
单位立体角就是单位面积与距离平方的比值。因此在这里我们可以类比得出相似的关系，类似Bllin-Phong模型，由光源投射出的光线方向与半球面有一个夹角，所以这里的A=Acosθ′A=Acos\theta'A=Acosθ′，而r自然是两点之间的距离r=x′−xr=x'-xr=x′−x。由此，我们得到dωid\omega_idωi和dAdAdA之前的关系
dω=dAcosθ′∣∣x′−x∣∣2cosθ′=nd\omega = \frac{dAcos\theta'}{||x'-x||^2} \ \ \ \ \ cos\theta'=n dω=∣∣x′−x∣∣2dAcosθ′ cosθ′=n
将上式代入渲染方程，我们就可以得到按光源采样的新蒙特卡洛积分形式
Lo(p,ωo)=∫ΩLi(p,ωi)fr(p,ωi,ωo)cosθdωi=∫ΩLi(p,ωi)fr(p,ωi,ωo)cosθcosθ′∣∣x′−x∣∣2dA≈1N∑i=1NLi(p,ωi)fr(p,ωi,ωo)cosθcosθ′∣∣x′−x∣∣2p(A)其中cosθ=n⋅dωcosθ′=n⋅(−dω)\begin{aligned} L_o(p,\omega_o) &=\int_{\Omega}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)cos\theta d\omega_i\\ &=\int_{\Omega}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)cos\theta\frac{cos\theta'}{||x'-x||^2}dA\\ &\approx\frac{1}{N}\sum_{i=1}^{N}\frac{L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)cos\theta\frac{cos\theta'}{||x'-x||^2}}{p(A)}\\ 其中&cos\theta=n\cdot d\omega \ \ \ \ \ cos\theta'=n\cdot(-d\omega)\\ \end{aligned} Lo(p,ωo)其中=∫ΩLi(p,ωi)fr(p,ωi,ωo)cosθdωi=∫ΩLi(p,ωi)fr(p,ωi,ωo)cosθ∣∣x′−x∣∣2cosθ′dA≈N1i=1∑Np(A)Li(p,ωi)fr(p,ωi,ωo)cosθ∣∣x′−x∣∣2cosθ′cosθ=n⋅dω cosθ′=n⋅(−dω)
对于上式，由pdf性质我们知道按p(A)p(A)p(A)对面积积分的结果应该是1（∫p(A)dA=1\int p(A)dA=1∫p(A)dA=1），自然p(A)=1Ap(A)=\frac{1}{A}p(A)=A1。

针对直接光照部分，按面积对光源采样的shade伪代码如下所示

shade(p, wo)# Contribution from the light source.Uniformly sample the light at x’ (pdf_light = 1 / A)L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light# Contribution from other reflectors.L_indir = 0.0Test Russian Roulette with probability P_RRUniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)Trace a ray r(p, wi)If ray r hit a non-emitting object at qL_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RRReturn L_dir + L_indir

除此之外，还有一个遗漏了的问题，那就是在光源上采样的光线必须得能射到着色点上，也就是说着色点ppp和光源x′x'x′之间不能有物体遮挡。

对上述伪代码稍作修改，判断一下是否有物体遮挡

shade(p, wo)# Contribution from the light source.Uniformly sample the light at x’ (pdf_light = 1 / A)Shoot a ray from p to x’If the ray is not blocked in the middle    // 判断是否有物体遮挡L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light# Contribution from other reflectors.L_indir = 0.0Test Russian Roulette with probability P_RRUniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)Trace a ray r(p, wi)If ray r hit a non-emitting object at qL_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RRReturn L_dir + L_indir

自此，我们终于得出了正确的，高效的渲染方式。

总结

为了尽可能展现我们看到的现实世界，实际上也就是渲染一张图片，再根本点就是渲染每一个像素。对每个像素进行多次采样，根据路径追踪得到的渲染方程计算每一次采样结果再累加平均，得到该像素理想的渲染结果。

对于渲染方程，是根据路径追踪得到的基于物理的渲染。与最初我们做的工作相比，之前是简单粗暴假设光线由半球均匀随机采样得到，存在递归指数爆炸，没有出口，效率低下等问题，而现在我们考虑radiance是有两部分组成，分别采用不同的方法计算。

直接光照，对光源按面积进行采样
间接光照，按半球随机采样一次，并且采用RR控制递归深度

Assignment 7

作业7要求我们用代码实现路径追踪算法。在此之前需要用到多个以前作业实现的函数，下面先贴出这些函数的代码，附有详细注解和函数作用讲解。

准备工作

Triangle::getIntersection in Triangle.hpp
判断该光线是否与三角形相交并获得进出时间和相交点的坐标、法线等属性，返回包含这些属性的结构体Intersection

inline Intersection Triangle::getIntersection(Ray ray)
{Intersection inter;if (dotProduct(ray.direction, normal) > 0)return inter;double u, v, t_tmp = 0;Vector3f pvec = crossProduct(ray.direction, e2);double det = dotProduct(e1, pvec);if (fabs(det) < EPSILON)return inter;double det_inv = 1. / det;Vector3f tvec = ray.origin - v0;u = dotProduct(tvec, pvec) * det_inv;if (u < 0 || u > 1)return inter;Vector3f qvec = crossProduct(tvec, e1);v = dotProduct(ray.direction, qvec) * det_inv;if (v < 0 || u + v > 1)return inter;t_tmp = dotProduct(e2, qvec) * det_inv;// TODO find ray triangle intersectionif (t_tmp < 0)return inter;inter.happened = true;inter.coords = ray(t_tmp);inter.normal = normal;inter.distance = t_tmp;inter.obj = this;inter.m = m;return inter;
}

IntersectP(const Ray& ray, const Vector3f& invDir, const std::array<int, 3>& dirIsNeg) in the Bounds3.hpp
判断光线ray是否与某一包围盒是否与相交，用于bvh递归对每个划分节点进行判断

inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,const std::array<int, 3>& dirIsNeg) const
{// invDir: ray direction(x,y,z), invDir=(1.0/x,1.0/y,1.0/z), use this because Multiply is faster that Division// dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x>0),int(y>0),int(z>0)], use this to simplify your logic// TODO test if ray bound intersectsfloat t_Min_x = (pMin.x - ray.origin.x) * invDir[0];float t_Min_y = (pMin.y - ray.origin.y) * invDir[1];float t_Min_z = (pMin.z - ray.origin.z) * invDir[2];float t_Max_x = (pMax.x - ray.origin.x) * invDir[0];float t_Max_y = (pMax.y - ray.origin.y) * invDir[1];float t_Max_z = (pMax.z - ray.origin.z) * invDir[2];if (dirIsNeg[0]) {float t = t_Min_x;t_Min_x = t_Max_x;t_Max_x = t;}if (dirIsNeg[1]) {float t = t_Min_y;t_Min_y = t_Max_y;t_Max_y = t;}if (dirIsNeg[2]) {float t = t_Min_z;t_Min_z = t_Max_z;t_Max_z = t;}float t_enter = std::max(t_Min_x, std::max(t_Min_y, t_Min_z));float t_exit = std::min(t_Max_x, std::min(t_Max_y, t_Max_z));// 这里需要注意判断t_enter == t_exit的情况// 若相等时返回false，由于light的y值一样所以厚度为0，则光线与灯相交时也为false，导致全黑if (t_enter <= t_exit && t_exit >= 0)return true;elsereturn false;
}

getIntersection(BVHBuildNode* node, const Ray ray)in BVH.cpp
加速结构，用于计算光线ray与物体相交的结果，可只对部分物体判断而不用遍历所有物体。通过bvh划分objects并递归判断，获得最小的distance也就是时间，得到最近的射线求交结果

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{// TODO Traverse the BVH to find intersectionIntersection inter;Vector3f invdir(1 / ray.direction.x, 1 / ray.direction.y, 1 / ray.direction.z);std::array<int, 3> dirIsNeg;dirIsNeg[0] = ray.direction.x < 0;dirIsNeg[1] = ray.direction.y < 0;dirIsNeg[2] = ray.direction.z < 0;// 若不与任何包围盒相交，则返回空也就是默认的interif (!node->bounds.IntersectP(ray, invdir, dirIsNeg))return inter;// 若为叶子节点则获得ray与该节点内triangle相交的结果if (node->left == nullptr && node->right == nullptr)return node->object->getIntersection(ray);// 对孩子节点进行递归，返回最近交点Intersection hit1 = getIntersection(node->left, ray);Intersection hit2 = getIntersection(node->right, ray);return hit1.distance < hit2.distance ? hit1 : hit2;
}

实现过程

在深入理解并正确实现上述代码以后，我们再来看一下本次作业框架提供的用于实现路径追踪算法的函数。

castRay(const Ray ray, int depth)in Scene.cpp:
在其中实现 Path Trac-ing 算法，可能用到的函数有：
intersect(const Ray ray)in Scene.cpp:
求一条光线与场景的交点
sampleLight(Intersection pos, float pdf) in Scene.cpp:
在场景的所有光源上按面积 uniform 地 sample 一个点，并计算该 sample 的概率密度
sample(const Vector3f wi, const Vector3f N) in Material.cpp:
按照该材质的性质，给定入射方向与法向量，用某种分布采样一个出射方向
pdf(const Vector3f wi, const Vector3f wo, const Vector3f N) in Material.cpp:
给定一对入射、出射方向与法向量，计算 sample 方法得到该出射方向的概率密度
eval(const Vector3f wi, const Vector3f wo, const Vector3f N) in Material.cpp:
给定一对入射、出射方向与法向量，计算这种情况下的 f_r 值

除此之外，可能用到的变量有：
RussianRoulette in Scene.cpp: P_RR, Russian Roulette 的概率

除此之外，本次作业由于框架原因所以对原伪代码做了一些修改（主要是wo定义与课程内容相反，由像素指向着色点，而不再是由着色点向外）

shade(p, wo)sampleLight(inter, pdf_light)Get x, ws, NN, emit from interShoot a ray from p to xIf the ray is not blocked in the middleL_dir = emit * eval(wo, ws, N) * dot(ws, N) * dot(ws, NN) / |x-p|^2 / pdf_lightL_indir = 0.0Test Russian Roulette with probability RussianRoulettewi = sample(wo, N)Trace a ray r(p, wi)If ray r hit a non-emitting object at qL_indir = shade(q, wi) * eval(wo, wi, N) * dot(wi, N) / pdf(wo, wi, N) / RussianRouletteReturn L_dir + L_indir

根据上述伪代码我们可以得到下面这个图，方便帮助我们理解并写代码

由于在前文学习过程中我们简化了渲染方程，忽略了物体本身发光的部分，所以如果按上述伪代码实现得到的渲染图中光源部分会是黑色的。很简单，如果弹射光线遇到发光的物体则直接返回，我们补上该部分即可得到正确的渲染图。

最终实现的castRay函数如下，附有详细注解

// Implementation of Path Tracing
Vector3f Scene::castRay(const Ray &ray, int depth) const
{// TO DO Implement Path Tracing Algorithm here// 计算直接光照部分Vector3f L_dir (0, 0, 0);// 从像素射出光线作为着色点的入射光线，获得渲染方程中入射光线与着色点相交的各项参数Intersection inter_obj = intersect(ray);if (!inter_obj.happened) // 需先判断从像素采样的该光纤是否有射到物体return L_dir;if (inter_obj.m->hasEmission()) // 如果追踪到光源则直接返回光源，保证光源是亮的return inter_obj.m->getEmission();Vector3f p = inter_obj.coords;Vector3f N = inter_obj.normal.normalized(); // 法线需要计算cos，所以需要单位化Vector3f wo = ray.direction;   // 由于框架原因，这里的wo为像素到着色点的方向，与着色点向外的方向相反// 计算反射光线与光源相交，为提高采样效率所以从光源出发采样Intersection inter_light;float pdf_light;sampleLight(inter_light, pdf_light); // 对光源按面积采样,直接填入默认参数，得到着色点与光源的交点Vector3f x = inter_light.coords;Vector3f ws = (x - p).normalized(); // 从物体指向光源// 从着色点射出光线，与着色点反射的光线比较，判断着色点反射的光线是否能射到光源Ray objTolight(p, ws);float d1 = (x - p).norm(); // 着色点与光源距离float d2 = intersect(objTolight).distance; // 着色点与相交物体距离// 计算L_dir直接光照if (d2 - d1 > -0.001) { // 做差判断浮点数是否相等Vector3f eval = inter_obj.m->eval(wo, ws, N); // 注意参数的方向，需要入射方向、出射方向Vector3f emit = inter_light.emit;   // 光源投射到着色点的光线的radianceVector3f NN = inter_light.normal.normalized();float cos_theta_obj = dotProduct(N, ws);float cos_theta_light = dotProduct(NN, -ws); // 这里要取负号，因为ws是从着色点指向光源，与NN夹角大于90度为负L_dir = emit * eval * cos_theta_obj * cos_theta_light / std::pow(d1, 2) / pdf_light;}// 递归计算间接光照部分Vector3f L_indir (0, 0, 0);float P_RR = get_random_float(); // 采用俄罗斯轮盘赌，避免无穷递归，最终数学期望一样if (P_RR < RussianRoulette) {Vector3f wi = inter_obj.m->sample(wo, N).normalized(); // 对着色点随机采样出射光线Ray ray_objToobj(p, wi); // 根据采样获得从着色点射向其他物体的光线Intersection inter = intersect(ray_objToobj);if (inter.happened && !inter.m->hasEmission()) { // 如果有交点且交点材质不发光Vector3f eval = inter_obj.m->eval(wo, wi, N); // 计算着色点的eval,而不是反射光线与其他物体交点的float pdf = inter_obj.m->pdf(wo, wi, N); // 计算着色点采样的pdffloat cos_theta = dotProduct(wi, N);L_indir = castRay(ray_objToobj, depth + 1) * eval * cos_theta / pdf / RussianRoulette;}}return L_dir + L_indir;
}

在实现过程中有几点值得注意

在判断直接光照弹射的光线是否有遮挡的时候，需要用浮点数比较着色点到光源的距离d1和着色点弹射光线与最近相交物体的距离d2是否相等。判断浮点数是否相等我们可以采用做差的方法，由于d1>=d2，所以需要d2-d1>=EPSION，我在这里的EPSION取值为-0.001，能够满足判断要求。
对于判断光线是否与某包围盒相交的函数IntersectP，在上课时老师讲到图形学里并不会太关注边界情况，所以t_enter<或者<=t_eixt都是没问题的，但是在本次实验中我们需要将边界条件设为小于等于，否则会渲染出的图片有部分区域是黑色的。

最终结果

设置不同的spp，能够得到不同的渲染效果。ssp越大，自然渲染的效果越好，花费的时间也更长。下面给出不同spp渲染的结果

spp=1

spp=16

spp=64

spp=256

提高部分

这里实现了多线程加速部分的内容，可以明显加快渲染过程。具体实现过程不再赘述，这里仅贴一下代码和不同线程下的渲染时长对比。

渲染时间	1	8	32	64	256
单线程	2m	21m	87m
多线程	32s	4m	16m	33m	145m

Render.cpp:

//
// Created by goksu on 2/25/20.
//
#include <fstream>#include "Scene.hpp"
#include "Renderer.hpp"
#include <atomic>
#include <thread>// 多线程加速渲染
inline float deg2rad(const float& deg) { return deg * M_PI / 180.0; }const float EPSILON = 0.00001;
std::atomic_int progress = 0;// The main render function. This where we iterate over all pixels in the image,
// generate primary rays and cast these rays into the scene. The content of the
// framebuffer is saved to a file.
void Renderer::Render(const Scene& scene)
{std::vector<Vector3f> framebuffer(scene.width * scene.height);float scale = tan(deg2rad(scene.fov * 0.5));float imageAspectRatio = scene.width / (float)scene.height;Vector3f eye_pos(278, 273, -800);// change the spp value to change sample ammount// spp： sample per pixelint spp = 256;  //原本16std::cout << "SPP: " << spp << "\n";int thred = 24;int per = scene.height/thred;  // 960/24=40std::thread th[24]; //多线程auto renderRow = [&](uint32_t lrow, uint32_t hrow){for (uint32_t j = lrow; j < hrow; ++j) {for (uint32_t i = 0; i < scene.width; ++i) {// generate primary ray directionfloat x = (2 * (i + 0.5) / (float)scene.width - 1) *imageAspectRatio * scale;float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;Vector3f dir = normalize(Vector3f(-x, y, 1));  // ??? (x,y,-1) ???for (int k = 0; k < spp; k++){framebuffer[(int)(j*scene.width+i)] += scene.castRay(Ray(eye_pos, dir), 0) / spp;  }}progress += 1;UpdateProgress(progress / (float)scene.height);}};for(int i=0;i<thred;i++){th[i] = std::thread(renderRow,i*per,(i+1)*per);}for(int i=0;i<thred;i++){th[i].join();}UpdateProgress(1.f);// save framebuffer to fileFILE* fp = fopen("binary.ppm", "wb");(void)fprintf(fp, "P6\n%d %d\n255\n", scene.width, scene.height);for (auto i = 0; i < scene.height * scene.width; ++i) {static unsigned char color[3];color[0] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].x), 0.6f));color[1] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].y), 0.6f));color[2] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].z), 0.6f));fwrite(color, 1, 3, fp);}fclose(fp);
}

参考

Games101：作业7（含提高部分）_Q_pril的博客-CSDN博客_games101 作业7

GAMES101: 现代计算机图形学入门

GAMES101蒙特卡洛光线追踪及Assignment7相关推荐

蒙特卡洛树搜索_蒙特卡洛树是什么算法？
点击上方"MLNLP",选择"星标"公众号重磅干货,第一时间送达编辑:忆臻 https://www.zhihu.com/question/39916945 ...
《一周学完光线追踪》学习十一点五离焦模糊代码原理分析
蒙特卡洛光线追踪技术系列见蒙特卡洛光线追踪技术首先分析一下生成随机Ray的程序: vec3 random_in_unit_disk() {vec3 p;do {p = 2.0*vec3(rand ...
计算机图形学——光线追踪（RayTracing）算法
转自:https://blog.csdn.net/hmbxsy/article/details/80509876?depth_1-utm_source=distribute.pc_relevant.n ...
《学一辈子光线追踪》三光散射
蒙特卡洛光线追踪技术系列见蒙特卡洛光线追踪技术这一小节全都是文字,但是比较重要.虽然之前在三维渲染光能辐射基础三维渲染 BRDF双向反射分布函数三维渲染相位函数三维渲染体照明模型 ...
《学一辈子光线追踪》十你的余生+极简后记
蒙特卡洛光线追踪技术系列见蒙特卡洛光线追踪技术这本书的目的是显示一种组织基于物理渲染器的采样方法的所有数学i的细节.现在你可以探索很多不同的潜在途径. 如果你想探索蒙特卡罗方法,可以研究双向和路 ...
《学一辈子光线追踪》十完成pdf管理
蒙特卡洛光线追踪技术系列见蒙特卡洛光线追踪技术到目前为止,我已经使用color函数创建了两个硬编码的pdf: 1.与光的形状有关的p0 2.与法向量和曲面类型相关的p1 我们可以将关于光的信息( ...
《一周学完光线追踪》学习三光线相机和背景
蒙特卡洛光线追踪技术系列见蒙特卡洛光线追踪技术所有光线跟踪器都有一个光线类,以及计算沿光线看到的颜色.让我们把射线看作一个函数,p(t)=A+t*B,这里p是3D中沿直线的3D位置,a是射线的原 ...
蒙特卡洛树是什么算法？
点击上方"小白学视觉",选择加"星标"或"置顶" 重磅干货,第一时间送达编辑:忆臻 https://www.zhihu.com/quest ...
《再学一周光线追踪》学习四 Perlin噪声
蒙特卡洛光线追踪技术系列见蒙特卡洛光线追踪技术为了获得酷酷的固体纹理,大多数人使用某种形式的柏林噪音.这些是以他们的发明家肯·佩林命名的.Perlin纹理不会返回这样的白噪声: 相反,它返回类似 ...

GAMES101蒙特卡洛光线追踪及Assignment7