Games101：作业7（含提高部分）

目录 

作业要求

任务分析

提高部分1：多线程

提高部分2：Microfacet

参考链接

---

作业要求

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

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

        • Triangle::getIntersection in Triangle.hpp: 将你的光线 -三角形相交函数粘贴到此处，请直接将上次实验中实现的内容粘贴在此。

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

        • getIntersection(BVHBuildNode* node, const Ray ray)in BVH.cpp: BVH 查找过程，请直接将上次实验中实现的内容粘贴在此处 .

        

        在本次实验中，你只需要修改这一个函数:

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

任务分析

        本次的任务即完整的实现整个path-tracing，但其实由于前面作业的铺垫，再本次作业中只需要实现castRay，即课上所讲的shade(p,w0)，如下图课件所示。

        但是上面的代码有两个问题：

        ①如果wo直接打到光源，则返回光源信号

        ②如果物体和光源之间有阻碍，则无法得到光源照射，见下图：

        当然，理论和实践之间还是有一定距离的， 本次作业的框架的定义与课件有所不同，具体体现在一些获取变量的函数所需传递的参数上，文档也已给出一个很详尽的伪代码：

         我再借助下面这张图说一下整个流程（图源自水印，侵删，网址在参考链接2）。首先我们已经得到了一条从像素打过来的光线wo：

        1）如果wo没有打到物体（即没有交点），返回（0,0,0）

        2）如果wo打到光源，返回光源信息

        3）如果wo打到物体，则要返回该物体对应的光照下的信息，该光照包括两部分，分别来自直接光源和反射光的间接光源：

                i.直接光源需要对当前环境下的光随机采样，得到一个交点和该采样结果的pdf。但是由于光源是随机采样，所以还需从物体向光源方向发射光线，看这个光线能否打到光源，以此判断是否有物体挡在此光源和物体之间

                ii.间接光源需要先利用俄罗斯轮盘赌判断间接光能否发生以此控制弹射次数，若能发生随机采样一个入射方向wi（这个随机采样实际是把wo当做入射方向，随机旋转一定角度得到“出射方向”wi，因为已知wo求wi），如果这个方向能够打到一个不发光的物体则返回其间接光照

         

         以上整个过程需要用到文档中所需要的各个函数，所以建议先理解其各个函数的意义，这里列出一些我标出注释的：

//判断光线是否与当前场景中的某个包围盒相交
Intersection Scene::intersect(const Ray &ray) const
{
    return this->bvh->Intersect(ray);
}

//对场景中的光源进行随机采样
void Scene::sampleLight(Intersection &pos, float &pdf) const
{
    float emit_area_sum = 0;
    //发光区域面积
    for (uint32_t k = 0; k < objects.size(); ++k) {
        if (objects[k]->hasEmit()){
            emit_area_sum += objects[k]->getArea();
        }
    }
    //对在场景的所有光源上按面积均匀 sample 一个点，并计算该 sample 的概率密度
    float p = get_random_float() * emit_area_sum;
    emit_area_sum = 0;
    for (uint32_t k = 0; k < objects.size(); ++k) {
        if (objects[k]->hasEmit()){
            emit_area_sum += objects[k]->getArea();
            if (p <= emit_area_sum){
                //对光源采样，得到位置和pdf(均匀采样 1/S)
                objects[k]->Sample(pos, pdf);
                break;
            }
            //写不写都行，具体运行之后无意义
            else
                break;
        }
    }
}

    void Sample(Intersection &pos, float &pdf){
        //2*PI
        float theta = 2.0 * M_PI * get_random_float();
        //PI
        float phi = M_PI * get_random_float();
        //cos,sin*cos,sin*sin,空间参数坐标系
        Vector3f dir(std::cos(phi), std::sin(phi)*std::cos(theta), std::sin(phi)*std::sin(theta));
        //圆心+半径*dir
        pos.coords = center + radius * dir;
        pos.normal = dir;
        pos.emit = m->getEmission();
        pdf = 1.0f / area;
    }

Vector3f Material::eval(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
    switch(m_type){
        case DIFFUSE:
        {
            // calculate the contribution of diffuse   model
            //这里认为diffuse均匀反射,每个方向的概率都是kd/PAI
            float cosalpha = dotProduct(N, wo);
            if (cosalpha > 0.0f) {
                Vector3f diffuse = Kd / M_PI;
                return diffuse;
            }
            else
                return Vector3f(0.0f);
            break;
        }
    }
}

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;
        }

        在代码实现中，需要注意浮点数的处理，在判断光线与物体之间是否有障碍物的时候，是比较光物距离①是否等于从物体向光源方向到相交物体的距离②。对于浮点数相等的比较一般是看是否在一个较小区间，又因为①一定大于等于②，所以①-②>= EPSION即可，在我的电脑运行过程中，EPSION取0.001最好。同样的在判断inline bool Bounds3::IntersectP时，需要写成t_enter<=t_exit，对浮点数处理，否则会导致部分地方是黑的：

         最终代码如下：

// Implementation of Path Tracing
Vector3f Scene::castRay(const Ray &ray, int depth) const
{
    // TO DO Implement Path Tracing Algorithm here
    Vector3f L_dir;
    Vector3f L_indir;

    // 从像素发出的光线与物体的交点
    Intersection obj_inter = intersect(ray);
    if(!obj_inter.happened)
        return L_dir;

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

    // 打到物体
    Vector3f p = obj_inter.coords;
    Material* m = obj_inter.m;
    Vector3f N = obj_inter.normal.normalized();
    Vector3f wo = ray.direction; // 像素到物体的向量
    
    // 有交点，对光源采样
    float pdf_L = 1.0; //可以不初始化
    Intersection light_inter ;
    sampleLight(light_inter,pdf_L);    // 得到光源位置和对光源采样的pdf
    
    Vector3f x = light_inter.coords;
    Vector3f ws = (x - p).normalized(); //物体到光源
    Vector3f NN = light_inter.normal.normalized();  
    Vector3f emit = light_inter.emit;
    float d = (x-p).norm();
    
    // 再次从光源发出一条光线，判断是否能打到该物体，即中间是否有阻挡
    Ray Obj2Light(p,ws);
    float d2 = intersect(Obj2Light).distance;
    // 是否阻挡，利用距离判断，需注意浮点数的处理
    if(d2-d > -0.001){
        Vector3f eval = m->eval(wo,ws,N); // wo不会用到
        float cos_theta = dotProduct(N,ws);
        float cos_theta_x = dotProduct(NN,-ws);//ws从物体指向光源，与NN的夹角大于180
        L_dir = emit * eval * cos_theta * cos_theta_x / std::pow(d,2) / pdf_L;
    }
    
    // L_indir
    float P_RR = get_random_float();
    if(P_RR<RussianRoulette){
        Vector3f wi = m->sample(wo,N).normalized();
        Ray r(p,wi);
        Intersection inter = intersect(r);
        // 判断打到的物体是否会发光取决于m
        if(inter.happened && !inter.m->hasEmission()){
            Vector3f eval = m->eval(wo,wi,N);
            float pdf_O = m->pdf(wo,wi,N);
            float cos_theta = dotProduct(wi,N);
            L_indir = castRay(r, depth+1) * eval * cos_theta/ pdf_O / RussianRoulette;
        }
    }
    //4->16min
    return L_dir + L_indir;
}

           下图为sps=16的情况：

提高部分1：多线程

        这部分可以学习一下多线程的使用C++11 多线程（std::thread）详解_sjc_0910的博客-CSDN博客_c++多线程

        需要注意的是对于progress的更新，由于各个线程并行执行访问progress，为了防止修改被吞掉，需要将其设成同步量。我这里使用的是 条形分隔（每一列 或 每一行 为一个线程，更利于编程）

//
// 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 = 16;  //原本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);    
}

sps = 4， 轻薄本，15min

sps = 8， 轻薄本，37min

sps = 16， 轻薄本，65min-70min

sps = 4，轻薄本，利用多线程加速，5min（4核）

sps = 16，轻薄本，利用多线程加速，16min（6核）

        说几点我遇到的问题，首先多线程编译会出现“undefined reference to pthread_create”，解决办法是修改CMakeList文件的配置，修改方式有很多，这里说一种，最终修改结果如下：

cmake_minimum_required(VERSION 3.10)
project(RayTracing)

set(CMAKE_CXX_STANDARD 17)
find_package(Threads) //引入外部依赖包

add_executable(RayTracing main.cpp Object.hpp Vector.cpp Vector.hpp Sphere.hpp global.hpp Triangle.hpp Scene.cpp
        Scene.hpp Light.hpp AreaLight.hpp BVH.cpp BVH.hpp Bounds3.hpp Ray.hpp Material.hpp Intersection.hpp
        Renderer.cpp Renderer.hpp)

target_link_libraries (${PROJECT_NAME} ${CMAKE_THREAD_LIBS_INIT}) //链接 Thread 库

        其次也是最重要的问题，编译正常之后，我最初几次运行时间都没有什么变化，最后和一个大佬交流之后，觉得应该是因为自己的虚拟机是单核的，虽然可以执行多线程，但是单核情况下多线程的切换频率最终在速度上没有收益。所以需要修改虚拟机的CPU分配，对于VirtualBox而言，只需在关机之后的主界面修改：

提高部分2：Microfacet

        这部分其实并不难，关键是理解Microfacet的BRDF的公式，因为课上没有详细介绍各个参数应该怎么计算，我强推下面这个链接，看不懂英文也问题不大，只需要看公式就行，

LearnOpenGL - Theory

        或者也可以看博客：

Games101,作业7（微表面模型）_Elsa的迷弟的博客-CSDN博客_games101 作业7

         原理部分我就不多搬运啦，我补充一下如何使用兔子模型，首先，如果你只是在场景中加个兔子，会发现完全没有任何效果。如果你在作业3的时候尝试过换模型就会知道是因为兔子模型太小了，并且对于cornell盒子这个模型部分身体甚至不在场景中。打印一下兔子部分的坐标（左下）和原模型中的地板坐标（右下），可以知道需要对兔子进行缩放和平移。这里也要感谢一下论坛中大佬们的讨论，详见最后一个参考链接。

         我选择的参数和论坛中所说的一致，translate用（300，0，300）， scale用（2000，2000，2000），最终代码如下：

   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))
    {
        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 = Scale*vert+Trans;
                face_vertices[j] = vert;

                min_vert = Vector3f(std::min(min_vert.x, vert.x),
                                    std::min(min_vert.y, vert.y),
                                    std::min(min_vert.z, vert.z));
                max_vert = Vector3f(std::max(max_vert.x, vert.x),
                                    std::max(max_vert.y, vert.y),
                                    std::max(max_vert.z, vert.z));
            }

            triangles.emplace_back(face_vertices[0], face_vertices[1],
                                   face_vertices[2], mt);
        }

        bounding_box = Bounds3(min_vert, max_vert);

        std::vector<Object*> ptrs;
        for (auto& tri : triangles){
            ptrs.push_back(&tri);
            area += tri.area;
        }
        bvh = new BVHAccel(ptrs);
    }

    ....
    Material* whiteM = new Material(Microfacet, Vector3f(0.0f));
    whiteM->Ks = Vector3f(0.45, 0.45, 0.45);
	whiteM->Kd = Vector3f(0.3, 0.3, 0.25);
    Material* light = new Material(DIFFUSE, (8.0f * Vector3f(0.747f+0.058f, 0.747f+0.258f, 0.747f) + 15.6f * Vector3f(0.740f+0.287f,0.740f+0.160f,0.740f) + 18.4f *Vector3f(0.737f+0.642f,0.737f+0.159f,0.737f)));
    light->Kd = Vector3f(0.65f);

    MeshTriangle floor("../models/cornellbox/floor.obj", white);
    // MeshTriangle shortbox("../models/cornellbox/shortbox.obj", white);
    // MeshTriangle tallbox("../models/cornellbox/tallbox.obj", white);
    MeshTriangle bunny("../models/bunny/bunny.obj", whiteM, Vector3f(300,0,300), 
        Vector3f(2000,2000,2000));
    MeshTriangle left("../models/cornellbox/left.obj", red);
    MeshTriangle right("../models/cornellbox/right.obj", green);
    MeshTriangle light_("../models/cornellbox/light.obj", light);

    scene.Add(&floor);
    // scene.Add(&shortbox);
    // scene.Add(&tallbox);
    scene.Add(&bunny);
    scene.Add(&left);
    scene.Add(&right);
    scene.Add(&light_);
    ...

        渲染结果如下，左下为Microfacet，右下为diffuse，都是sps=4的结果，所以diffuse的噪点有点多。

一些疑惑

参考链接

C++11 多线程（std::thread）详解_sjc_0910的博客-CSDN博客_c++多线程

Games 101 | 作业7 + 路径追踪 Path Tracing + 多线程 - 知乎

cmake undefined reference to `pthread_create‘_早睡的叶子的博客-CSDN博客

Cmake编译pthreads报错：undefined reference to pthread_create_v俊逸的博客-CSDN博客_cmake pthread_create

LearnOpenGL - Theory

Games101,作业7（微表面模型）_Elsa的迷弟的博客-CSDN博客_games101 作业7

作业7 请问如何添加其他模型 – 计算机图形学与混合现实在线平台