文章目录

前言
蒙特卡洛路径追踪
Assignment 7
参考

前言

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

蒙特卡洛路径追踪

蒙特卡洛积分

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

PDF具有如下特点

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

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

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

一些需要注意的地方：

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

路径追踪

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

Whitted-Style光线追踪的局限性

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

渲染方程

根据辐射度量学，我们一通操作最终求得渲染方程，即在物理意义上得出的正确的能够描述我们现实生活所见的方程。
$L_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$

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

为了方便计算和统一，这里将着色点本身产生的 $L_e(p,\omega_o)$ 省略掉，我们更关注如下这个渲染方程（或者说反射方程）。
$L_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$
根据蒙特卡洛积分可以得到
$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)}$
这个方程可以理解为按 $p(\omega_i)$ 的概率随机选取N个方向 $\omega_i$ 采样，再累加求平均。

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

直接光照

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

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

shade(p, wo)
	Randomly choose N directions wi~pdf
	Lo = 0.0
	For each wi
		Trace 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~pdf
	Lo = 0.0
	For each wi
		Trace a ray r(p, wi)
		If ray r hit the light
			Lo += (1 / N) * L_i * f_r * cosine / pdf(wi)
		Else If ray r hit an object at q
			Lo += (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个点，最终导致 $N^x$ 的指数增长。但有一个数是可以避免这种情况，那就是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 light
		Return L_i * f_r * cosine / pdf(wi)
	Else If ray r hit an object at q
		Return 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 pixel
    pixel_radiance = 0.0
    For each sample in the pixel
    	Shoot a ray r(camPos, cam_to_sample)
    	If ray r hit the scene at p
    		pixel_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 * (L o / P) + (1 - P) * 0 = L o$
由此，我们的shade伪代码再次改动如下

shade(p, wo)
    Manually specify a probability P_RR
    Randomly 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 light
    	Return L_i * f_r * cosine / pdf(wi) / P_RR
    Else If ray r hit an object at q
   		Return shade(q, -wi) * f_r * cosine / pdf(wi) / P_RR

优化

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

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

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

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

回顾以前对立体角 $\omega$ 的定义
$d\omega = \frac{dA}{r^2}$
单位立体角就是单位面积与距离平方的比值。因此在这里我们可以类比得出相似的关系，类似Bllin-Phong模型，由光源投射出的光线方向与半球面有一个夹角，所以这里的 $A=Acos\theta'$ ，而r自然是两点之间的距离 $r = x^{'} - x$ 。由此，我们得到 $d\omega_i$ 和 $d A$ 之前的关系
$d\omega = \frac{dAcos\theta'}{||x'-x||^2} \ \ \ \ \ cos\theta'=n$
将上式代入渲染方程，我们就可以得到按光源采样的新蒙特卡洛积分形式
$\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}$
对于上式，由pdf性质我们知道按 $p (A)$ 对面积积分的结果应该是1（ $\int p(A)dA=1$ ），自然 $p(A)=\frac{1}{A}$ 。

针对直接光照部分，按面积对光源采样的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.0
    Test Russian Roulette with probability P_RR
    Uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)
    Trace a ray r(p, wi)
    If ray r hit a non-emitting object at q
    	L_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RR
    
    Return L_dir + L_indir

除此之外，还有一个遗漏了的问题，那就是在光源上采样的光线必须得能射到着色点上，也就是说着色点 $p$ 和光源 $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.0
    Test Russian Roulette with probability P_RR
    Uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)
    Trace a ray r(p, wi)
    If ray r hit a non-emitting object at q
    	L_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RR
    
    Return 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 intersection
    if (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 intersects
    float 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;
    else
        return 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 intersection
    Intersection 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;
    // 若不与任何包围盒相交，则返回空也就是默认的inter
    if (!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 inter
    Shoot a ray from p to x
    If the ray is not blocked in the middle
    	L_dir = emit * eval(wo, ws, N) * dot(ws, N) * dot(ws, NN) / |x-p|^2 / pdf_light

    L_indir = 0.0
    Test Russian Roulette with probability RussianRoulette
    wi = sample(wo, N)
    Trace a ray r(p, wi)
    If ray r hit a non-emitting object at q
    	L_indir = shade(q, wi) * eval(wo, wi, N) * dot(wi, N) / pdf(wo, wi, N) / RussianRoulette

    Return 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;   // 光源投射到着色点的光线的radiance
        Vector3f 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); // 计算着色点采样的pdf
            float 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 pixel
    int spp = 256;  //原本16
    std::cout << "SPP: " << spp << "\n";
    int thred = 24;
    int per = scene.height/thred;  // 960/24=40
    std::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 direction
                float 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 file
    FILE* 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