【GAMES101】作业7（提高）路径追踪 多线程、Microfacet（全镜面反射）、抗锯齿

1. 作业描述

在之前的练习中，我们实现了 Whitted-Style Ray Tracing 算法，并且用 BVH等加速结构对于求交过程进行了加速。在本次实验中，我们将在上一次实验的基础上实现完整的 Path Tracing 算法。至此，我们已经来到了光线追踪版块的最后一节内容。

1.1 修改的内容

相比上一次实验，本次实验对框架的修改较大，主要在以下几方面：

• 修改了 main.cpp，以适应本次实验的测试模型 CornellBox
• 修改了 Render，以适应 CornellBox 并且支持 Path Tracing 需要的同一 Pixel多次 Sample
• 修改了 Object，Sphere，Triangle，TriangleMesh，BVH，添加了 area 属性与Sample方法，以实现对光源按面积采样，并在 Scene 中添加了采样光源的接口 sampleLight
• 修改了 Material 并在其中实现了 sample, eval, pdf 三个方法用于 Path Tracing 变量的辅助计算

1.2 你需要迁移的内容

你需要从上一次编程练习中直接拷贝以下函数到对应位置：

• Triangle::getIntersection in Triangle.hpp:
  将你的光线-三角形相交函数粘贴到此处，请直接将上次实验中实现的内容粘贴在此。
• IntersectP(const Ray& ray, const Vector3f& invDir,const std::array<int, 3>& dirIsNeg) in the Bounds3.hpp:
  这个函数的作用是判断包围盒 BoundingBox 与光线是否相交，请直接将上次实验中实现的内容粘贴在此处，并且注意检查 t_enter = t_exit的时候的判断是否正确。
• getIntersection(BVHBuildNode* node, const Ray ray)in BVH.cpp:
  BVH查找过程，请直接将上次实验中实现的内容粘贴在此处

1.3 代码框架

在本次实验中，你只需要修改这一个函数: • castRay(const Ray ray, int depth)in Scene.cpp: 在其中实现 Path Tracing 算法
可能用到的函数有：

• 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 的概率

2. 解

2.1 迁移部分（Triangle::getIntersection、Bounds3::IntersectP、BVHAccel::getIntersection）

这部分是将上次作业的代码直接移植过来，没有特别要修改的地方

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
    inter.happened = true;
    inter.obj = this;
    inter.distance = t_tmp;
    inter.normal = normal;
    inter.coords = ray(t_tmp);
    inter.m = this->m;
    return inter;
}

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
    const auto& origin = ray.origin;
	float ten = -std::numeric_limits<float>::infinity();
	float tex = std::numeric_limits<float>::infinity();
	for (int i = 0; i < 3; i++)
	{
		float min = (pMin[i] - origin[i]) * invDir[i];
		float max = (pMax[i] - origin[i]) * invDir[i];
		if (!dirIsNeg[i])
			std::swap(min, max);
		ten = std::max(min, ten);
		tex = std::min(max, tex);
    }
    return ten <= tex && tex >= 0;
}

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    // TODO Traverse the BVH to find intersection
    Intersection intersect, intersectl, intersectr;
    std::array<int, 3> dirIsNeg;
	dirIsNeg[0] = int(ray.direction.x >= 0);
	dirIsNeg[1] = int(ray.direction.y >= 0);
	dirIsNeg[2] = int(ray.direction.z >= 0);
    if(!node->bounds.IntersectP(ray, ray.direction_inv))
        return intersect;
    if(node->left == nullptr && node->right == nullptr){
        intersect = node->object->getIntersection(ray);
        return intersect;
    }
    intersectl = getIntersection(node->left,ray);
    intersectr = getIntersection(node->right,ray);
    return intersectl.distance < intersectr.distance ? intersectl : intersectr;
}

2.2 修改部分（Scene::castRay）

2.2.1 伪代码

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

2.2.2 实现

首先就是利用intersect函数判断光线与场景的交点，如果没有交点自然就不用继续往下求了，如果有交点的话判断是不是打到光源上，是的话也是直接返回光源颜色，因为这里默认光源反射其他方向光线的部分可以忽略不计

// Implementation of Path Tracing
Vector3f Scene::castRay(const Ray &ray, int depth) const
{
    // TO DO Implement Path Tracing Algorithm here
    Intersection intersec = intersect(ray);
    if (!intersec.happened) {
        return Vector3f();
    }

    // 打到光源
    if (intersec.m->hasEmission()) {
        return intersec.m->getEmission();
    }

    Vector3f l_dir(0,0,0);
    Vector3f l_indir(0,0,0);
    ...

接下来是直接光照的部分，这里采用的是对光源求积分的方式（用蒙特卡洛积分简化了），注意要判断光线是否被遮挡的问题（对光源采样出来的光线做一次求交，如果交点距离小于到光源的距离，说明被遮挡了）：

    ...
    // 直接光照
    Intersection lightInter;
    float lightPdf = 0.0f;
    sampleLight(lightInter, lightPdf);

    Vector3f obj2light = lightInter.coords - intersec.coords;
    Vector3f obj2lightDir = obj2light.normalized();
    float obj2lightPow = obj2light.x * obj2light.x + obj2light.y * obj2light.y + obj2light.z * obj2light.z;

    Ray obj2lightRay(intersec.coords, obj2lightDir);
    Intersection t = intersect(obj2lightRay);
    if (t.distance - obj2light.norm() > -EPSILON)
    {
        l_dir = lightInter.emit * intersec.m->eval(ray.direction, obj2lightDir, intersec.normal) 
            * dotProduct(obj2lightDir, intersec.normal) 
            * dotProduct(-obj2lightDir, lightInter.normal) 
            / obj2lightPow / lightPdf;
    }
    ...

接下来为了保证光线不会无限反射，用俄罗斯轮盘赌的方式决定光线是否继续：

	...
    if (get_random_float() > RussianRoulette) {
        return l_dir;
    }
    ...

若光线存活，这继续对间接光照的求解，因为已经对光源进行积分了，所以这里的间接光照求的是不发光的物体反射的光线：

（因为求的是随机采样的结果，为了保证能量守恒，需要对结果再除一个俄罗斯轮盘赌的概率）

    // 间接光照
    ...
    Vector3f obj2nextobjdir = intersec.m->sample(ray.direction, intersec.normal).normalized();
    Ray obj2nextobjray(intersec.coords, obj2nextobjdir);
    Intersection nextObjInter = intersect(obj2nextobjray);
    if (nextObjInter.happened && !nextObjInter.m->hasEmission())
    {
        float pdf = intersec.m->pdf(ray.direction, obj2nextobjdir, intersec.normal);
        l_indir = castRay(obj2nextobjray, depth + 1) 
            * intersec.m->eval(ray.direction, obj2nextobjdir, intersec.normal) 
            * dotProduct(obj2nextobjdir, intersec.normal)
            / pdf / RussianRoulette;
    }
    return l_dir + l_indir;
}

3. 提高

3.1 多线程

为提升程序执行效率，我们可以考虑利用多线程并发执行，这里我们可以从render函数发射primary ray的时候入手，将屏幕像素分成多块给多个线程执行，比如我们的scene尺寸为784*784，要用32个线程并发执行时，就将每块设置为(784/32) * 784的大小，接下来介绍两种多线程的做法：

3.1.1 std::thread

（出现undefined reference to `pthread_create’或其他编译问题的可以划到下面的第四部分看解决方法哦）

使用的头文件：

#include < thread >
#include < mutex >

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);
    int m = 0;
    int thread_num = 32;
    int thread_step = scene.height / thread_num;
    std::vector<std::thread> rays;
    // change the spp value to change sample ammount
    int spp = 1024;
    std::cout << "SPP: " << spp << "\n";
    for (int i = 0; i < thread_num; i++) 
        rays.push_back(std::thread(para, eye_pos, std::ref(framebuffer), std::ref(scene), spp, 
                    imageAspectRatio, scale, i * thread_step, (i + 1) * thread_step));
    for (int i = 0; i < thread_num; i++) 
        rays[i].join();
    UpdateProgress(1.f);
    ...

这里我们定义一个para函数作为多线程的入口，别忘记定义一个mutex锁用来锁住进度条更新相关语句，用来避免多线程同时访问全局变量的时候出现冲突

int prog = 0;
std::mutex lock;

void para(Vector3f eye_pos, std::vector<Vector3f> &framebuffer, const Scene& scene, int spp, float imageAspectRatio, float scale, int start, int end){
    int width, height;
    width = height = sqrt(spp);
    float step = 1.0f / width;
    for (uint32_t j = start; j < end; ++j) {
        for (uint32_t i = 0; i < scene.width; ++i) {
            // generate primary ray direction   
            for (int k = 0; k < spp; k++){
                float x = (2 * (i + step / 2 + step * (k % width)) / (float)scene.width - 1) *
                        imageAspectRatio * scale;
                float y = (1 - 2 * (j + step / 2 + step * (k / height)) / (float)scene.height) * scale;
                Vector3f dir = normalize(Vector3f(-x, y, 1));
                framebuffer[j * scene.width + i] += scene.castRay(Ray(eye_pos, dir), 0) / spp;  
            }
        }
        lock.lock();
        prog++;
        UpdateProgress(prog / (float)scene.height);
        lock.unlock();
    }
}

同时别忘记用join方法等待所有线程结束，防止有的线程还没结束主程序就结束了

for (int k = 0; k < spp; k++)
        rays[k].join();

3.1.2 OpenMP （推荐，更快的方法！！！）

使用的头文件：
#include <omp.h>

操作非常简单：
在你需要并行化的for前面，加上

#pragma omp parallel for

它可以将跟在后面的for循环语句分成多个线程并发执行

然后在cmakelists.txt中加上

set(CMAKE_CXX_FLAGS "${CAMKE_CXX_FLAGS} -O3 -fopenmp")

重新执行cmake …然后编译即可。
上面的-O3也可以提速。

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);
    int m = 0;
    int thread_num = 32;
    int thread_step = scene.height / thread_num;
    std::vector<std::thread> rays;
    // change the spp value to change sample ammount
    int spp = 4;
    std::cout << "SPP: " << spp << "\n";
    #pragma omp parallel for
        for (int i = 0; i < thread_num; i++) 
            para(eye_pos, std::ref(framebuffer), std::ref(scene), spp, 
                    imageAspectRatio, scale, i * thread_step, (i + 1) * thread_step);
    UpdateProgress(1.f);

更新进度条时加锁：

		omp_set_lock(&lock1);
        prog++;
        UpdateProgress(prog / (float)scene.height);
        omp_unset_lock(&lock1);

同样也要记住把锁设置为全局变量：

omp_lock_t lock1;

3.1.3 效果（32线程，spp=4）

不用多线程：

std::thread（32线程）：

openMP：

可见，openMP简直就是为本次作业量身定做的方法

3.2 Microfacet

3.2.1 全镜面反射BRDF分析

对于一个全镜面反射来说, 只需要考虑菲涅反射的系数, 也就是说

这里 wr 和 wo 是关于表面法线对称的.

现在的问题是, 我们如何描述这里的 fr BRDF项? 答案是使用狄拉克 函数, 狄拉克 函数满足:

也就是只能接收到唯一一个方向的光线

这样, 将狄拉克函数带入到方程中, 得到:

得到BRDF的表示为:

3.2.2 基本设置

首先在Material枚举类中添加MIRROR材质

enum MaterialType { DIFFUSE, MIRROR };

然后在MAIN函数里设置MIRROR参数，并且将两个长方体设置为MIRROR材质

Material* white1 = new Material(MIRROR, Vector3f(0.0f));
white1->Kd = Vector3f(0.0f, 0.0f, 0.0f);
white1->ior = 40.0f;
MeshTriangle shortbox("../models/cornellbox/shortbox.obj", white1);
MeshTriangle tallbox("../models/cornellbox/tallbox.obj", white1);

3.2.3 采样方向

全镜面反射只有一个方向的光线能被眼睛接收，所以采样函数就只返回反射方向

Vector3f Material::sample(const Vector3f &wi, const Vector3f &N){
    switch(m_type){
        case DIFFUSE:
        {
            // uniform sample on the hemisphere
            float x_1 = get_random_float(), x_2 = get_random_float();
            float z = std::fabs(1.0f - 2.0f * x_1);
            float r = std::sqrt(1.0f - z * z), phi = 2 * M_PI * x_2;
            Vector3f localRay(r*std::cos(phi), r*std::sin(phi), z);
            return toWorld(localRay, N);
            break;
        }
        case MIRROR:
        {
            Vector3f localRay = reflect(wi, N);
            return localRay;
            break;
        }
    }
}

Vector3f reflect(const Vector3f &I, const Vector3f &N) const
{
	return I - 2 * dotProduct(I, N) * N;
}

3.2.4 概率密度函数

因为全镜面反射只有一个方向的光线能被眼睛接收，所以pdf就设置为1

float Material::pdf(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
    switch(m_type){
        case DIFFUSE:
        {
            // uniform sample probability 1 / (2 * PI)
            if (dotProduct(wo, N) > 0.0f)
                return 0.5f / M_PI;
            else
                return 0.0f;
            break;
        }
        case MIRROR:
        {
            if (dotProduct(wo, N) > EPSILON)
                return 1.0f;
            else
                return 0.0f;
            break;
        }
    }
}

3.2.5 BRDF

为了保证最终结果只和菲涅尔项和反射光线有关，brdf里还要抵消掉cosθi的影响：

Vector3f Material::eval(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
    switch(m_type){
        case DIFFUSE:
        {
            // calculate the contribution of diffuse   model
            float cosalpha = dotProduct(N, wo);
            if (cosalpha > 0.0f) {
                Vector3f diffuse = Kd / M_PI;
                return diffuse;
            }
            else
                return Vector3f(0.0f);
            break;
        }
        case MIRROR:
        {
            float cosalpha = dotProduct(N, wo);
            float kr;
            if (cosalpha > EPSILON) {
                fresnel(wi, N, ior, kr);
                Vector3f mirror = 1 / cosalpha;
                return kr * mirror;        
            }
            else
                return Vector3f(0.0f);
            break;
        }
    }
}

void fresnel(const Vector3f &I, const Vector3f &N, const float &ior, float &kr) const
{
	float cosi = clamp(-1, 1, dotProduct(I, N));
	float etai = 1, etat = ior;
	if (cosi > 0) {  std::swap(etai, etat); }
	// Compute sini using Snell's law
	float sint = etai / etat * sqrtf(std::max(0.f, 1 - cosi * cosi));
	// Total internal reflection
	if (sint >= 1) {
		kr = 1;
	}
	else {
		float cost = sqrtf(std::max(0.f, 1 - sint * sint));
		cosi = fabsf(cosi);
		float Rs = ((etat * cosi) - (etai * cost)) / ((etat * cosi) + (etai * cost));
		float Rp = ((etai * cosi) - (etat * cost)) / ((etai * cosi) + (etat * cost));
		kr = (Rs * Rs + Rp * Rp) / 2;
	}
}

3.2.7 castRay的修改

当然，完成了镜面反射的BRDF，我们其实还有一件事要做，还记得之前的渲染方程里面我们计算了直接光照和间接光照吗，记得直接光照是对光源积分吗？没错，我们这里不能再对光源积分了，因为镜面材质不像漫反射那样会把各个方向光源的光都反射过来，加上对光源积分可能会导致面向光源的面过曝全白（感谢评论区 @木木是小呆呆 同学反馈的图片）：

所以我们在castRay里的直接和间接光照计算需要分Diffuse和Mirror两个情况来讨论，对于Mirror材质，我们直接把直接光照设置为0，但是间接光照要注意不能照搬，之前我们对于Diffuse材质的间接光照是只对不发光的物体采样积分，这是因为我们在直接光照里已经对光源积分了，但是Mirror我们并没有计算直接光照，所以别忘了把非发光物体的判断条件给去掉，我们接收所有物体入射的光！

// Implementation of Path Tracing
Vector3f Scene::castRay(const Ray &ray, int depth) const
{
    // TO DO Implement Path Tracing Algorithm here
    Intersection intersec = intersect(ray);
    if (!intersec.happened) {
        return Vector3f();
    }

    // 打到光源
    if (intersec.m->hasEmission()) {
        return intersec.m->getEmission();
    }


    Vector3f l_dir(0,0,0);
    Vector3f l_indir(0,0,0);
    switch(intersec.m->getType()){
        case DIFFUSE:{
            // 对光源积分
            Intersection lightInter;
            float lightPdf = 0.0f;

            sampleLight(lightInter, lightPdf);

            Vector3f obj2light = lightInter.coords - intersec.coords;
            Vector3f obj2lightDir = obj2light.normalized();
            float obj2lightPow = obj2light.x * obj2light.x + obj2light.y * obj2light.y + obj2light.z * obj2light.z;

            Ray obj2lightRay(intersec.coords, obj2lightDir);
            Intersection t = intersect(obj2lightRay);
            if (t.distance - obj2light.norm() > -EPSILON)
            {
                l_dir = lightInter.emit * intersec.m->eval(ray.direction, obj2lightDir, intersec.normal) 
                    * dotProduct(obj2lightDir, intersec.normal) 
                    * dotProduct(-obj2lightDir, lightInter.normal) 
                    / obj2lightPow / lightPdf;
            }

            if (get_random_float() > RussianRoulette) {
                return l_dir;
            }

            // 对其他方向积分
            Vector3f obj2nextobjdir = intersec.m->sample(ray.direction, intersec.normal).normalized();
            Ray obj2nextobjray(intersec.coords, obj2nextobjdir);
            Intersection nextObjInter = intersect(obj2nextobjray);
            if (nextObjInter.happened && !nextObjInter.m->hasEmission())
            {
                float pdf = intersec.m->pdf(ray.direction, obj2nextobjdir, intersec.normal);
                if (pdf > EPSILON)
                {
                    l_indir = castRay(obj2nextobjray, depth + 1) 
                        * intersec.m->eval(ray.direction, obj2nextobjdir, intersec.normal) 
                        * dotProduct(obj2nextobjdir, intersec.normal)
                        / pdf / RussianRoulette;
                }
            }           
            break;
        }
        case MIRROR:{
            if (get_random_float() > RussianRoulette) {
                return l_dir;
            }
            Vector3f obj2nextobjdir = intersec.m->sample(ray.direction, intersec.normal).normalized();
            Ray obj2nextobjray(intersec.coords, obj2nextobjdir);
            Intersection nextObjInter = intersect(obj2nextobjray);
            if (nextObjInter.happened)
            {
                float pdf = intersec.m->pdf(ray.direction, obj2nextobjdir, intersec.normal);
                if (pdf > EPSILON)
                {
                    l_indir = castRay(obj2nextobjray, depth + 1) 
                        * intersec.m->eval(ray.direction, obj2nextobjdir, intersec.normal) 
                        * dotProduct(obj2nextobjdir, intersec.normal)
                        / pdf / RussianRoulette;
                }
            }
            break;
        }
    }
    return l_dir + l_indir;
}

3.2.8 最终效果

镜面反射到墙上的高光区域的效果不太好，看起来像是一堆噪点，是因为该区域在进行采样的时候直接光照只对光源积分，但是对于全镜面反射过来的光源也应该看做一种光源，然后进行专门的重要性采样，而不是当作漫反射光源进行采样，这样的话可能只有较少的概率采样到镜面反射过来的光源，也就是图中很明显的离散的高光区域，这个问题可以通过提升spp进行解决，如下图，但是效率真的比较差：

除了针对镜面反射的光源的重要性采样之外，为了进一步优化视觉效果，还可以加入伽马矫正以符合人眼色彩观测经验，大概效果如下，因为时间问题我没有继续做了，有兴趣的朋友可以自己去尝试一下：

3.3 MSAA抗锯齿

原代码中只是重复计算spp次从一个像素发出的光线，最终取平均而已，但是这样就只是得到该点像素中心比较接近现实光追的颜色，但是对于计算机显示来说，他并没有解决该点像素周围的平滑过渡问题，比如图中两个长方体的边界：

为了解决这个问题，我们在对一个像素进行spp次采样的同时将这个像素分为spp个小像素，并从这些像素的中心发出primary ray，这样每个像素的颜色就可以实现和周围像素的平滑过渡了：

for (int k = 0; k < spp; k++){
	float x = (2 * (i + step / 2 + step * (k % width)) / (float)scene.width - 1) *
	        imageAspectRatio * scale;
	float y = (1 - 2 * (j + step / 2 + step * (k / height)) / (float)scene.height) * scale;
	Vector3f dir = normalize(Vector3f(-x, y, 1));
	framebuffer[m] += scene.castRay(Ray(eye_pos, dir), 0) / spp;  
}

3.4 其他模型

添加作业自带的bunny模型时，记得要对模型进行缩放，平移和旋转，否则模型将无法在镜头里面显示，这里给出我用到的参数：平移是(200,-60,150)，缩放是Vector3f(1500,1500,1500)，旋转是绕y轴180°

我们直接对MeshTriangle的构造函数修改一下就行了：

MeshTriangle bunny("../models/bunny/bunny.obj", white1, Vector3f(200,-60,150), 
        Vector3f(1500,1500,1500), Vector3f(-1,0,0), Vector3f(0,1,0), Vector3f(0,0,-1));

MeshTriangle(const std::string& filename, Material *mt = new Material(),
                Vector3f Trans = Vector3f(0.0,0.0,0.0), Vector3f Scale = Vector3f(1.0,1.0,1.0), 
                Vector3f xr = Vector3f(1.0,0,0), Vector3f yr = Vector3f(0,1.0,0),  Vector3f zr = Vector3f(0,0,1))
    {
        objl::Loader loader;
        loader.LoadFile(filename);
        area = 0;
        m = mt;
        assert(loader.LoadedMeshes.size() == 1);
        auto mesh = loader.LoadedMeshes[0];

        Vector3f min_vert = Vector3f{std::numeric_limits<float>::infinity(),
                                     std::numeric_limits<float>::infinity(),
                                     std::numeric_limits<float>::infinity()};
        Vector3f max_vert = Vector3f{-std::numeric_limits<float>::infinity(),
                                     -std::numeric_limits<float>::infinity(),
                                     -std::numeric_limits<float>::infinity()};
        for (int i = 0; i < mesh.Vertices.size(); i += 3) {
            std::array<Vector3f, 3> face_vertices;

            for (int j = 0; j < 3; j++) {
                auto vert = Vector3f(mesh.Vertices[i + j].Position.X,
                                     mesh.Vertices[i + j].Position.Y,
                                     mesh.Vertices[i + j].Position.Z);
                                     
                vert.x = dotProduct(vert, xr);
                vert.y = dotProduct(vert, yr);
                vert.z = dotProduct(vert, zr);//旋转
                vert = Scale*vert+Trans;//平移，缩放
                ...

4. 一些可能遇到的问题

4.1 渲染耗时过长？

使用了多线程，但是一次spp为1000以上的渲染还是要几个小时？可能是global.hpp下的get_random_float()随机数生成函数存在问题，它会导致在重复调用该函数时，返回同一个值。

解决方法：把其中定义的dev，rng，dist变量定义为静态变量（加个static修饰），这样最后的时间消耗就缩短了大概几十倍

4.2 渲染结果中光源区域为纯黑？

作业给的伪代码中缺少了光线直接与光源相交的部分，所以是纯黑，记得加上这部分

4.3 渲染结果较暗？

如果有渲染结果较暗，出现横向黑色条纹的情况，那么，很可能是因为直接光部分由于精度问题，被错误判断为遮挡，可以试着通过精度调整放宽相交限制（将EPSILON变量增大），同时可能因为老师课上说浮点数相等可能性很低，所以在判断条件中设置为 t_enter < t_exit而忘了等号，所以可能会有一定的出入，这次作业一定不要忘记加上

4.4 天花板黑色，墙面没有影子？

同样的，因为cornell box是由墙壁面组成的，使得包围盒高度为0，所以 t_exit >= 0也不要忘了等号，否则会导致天花板黑色，且盒子在地板和墙面没有影子

4.5 多线程编译时出现undefined reference to `pthread_create’？

pthread 库不是 Linux 系统默认的库，连接时需要使用静态库 libpthread.a.

所以在使用pthread_create()创建线程时，需要链接该库。

解决方法是打开CMakelists.txt，添加下列语句：

cmake_minimum_required (VERSION 2.6)
find_package (Threads)
add_executable (myapp main.cpp …)
target_link_libraries (myapp ${CMAKE_THREAD_LIBS_INIT})

本作业中为：

同时还要注意，往多线程的启动函数里传入参数的时候，需要引用的参数一定要用std::ref来拷贝，否则构析的时候会出错，如：

rays.push_back(std::thread(para, eye_pos, std::ref(framebuffer), std::ref(scene), spp, 
                    		imageAspectRatio, scale, i * thread_step, (i + 1) * thread_step));

4.6 渲染出现白色噪点？

原因应该是pdf接近于0时，除以它计算得到的颜色会偏向极限值，体现在图上也就是白色

要解决这个问题，对于pdf接近于0的情况直接将它的radience算作0就行：

float pdf = intersec.m->pdf(ray.direction, obj2nextobjdir, intersec.normal);
if (pdf > EPSILON)
{
     l_indir = castRay(obj2nextobjray, depth + 1) 
     * intersec.m->eval(ray.direction, obj2nextobjdir, intersec.normal) 
     * dotProduct(obj2nextobjdir, intersec.normal)
     / pdf / RussianRoulette;
}

优化后效果：

5. 效果集合

6. 附件

附上源代码，有兴趣的朋友可以自己尝试一下效果：
CSDN:【GAMES101】作业7
GITHUB:【GAMES101】作业合集