【GAMES101】作業7（提高）路徑追蹤 多線程、Microfacet（全鏡面反射）、抗鋸齒

目錄

• 1. 作業描述
• • 1.1 修改的内容
  • 1.2 你需要遷移的内容
  • 1.3 代碼架構
• 2. 解
• • 2.1 遷移部分（Triangle::getIntersection、Bounds3::IntersectP、BVHAccel::getIntersection）
  • 2.2 修改部分（Scene::castRay）
  • • 2.2.1 僞代碼
    • 2.2.2 實作
• 3. 提高
• • 3.1 多線程
  • • 3.1.1 std::thread
    • 3.1.2 OpenMP （推薦，更快的方法！！！）
    • 3.1.3 效果（32線程，spp=4）
  • 3.2 Microfacet
  • • 3.2.1 全鏡面反射BRDF分析
    • 3.2.2 基本設定
    • 3.2.3 采樣方向
    • 3.2.4 機率密度函數
    • 3.2.5 BRDF
    • 3.2.7 castRay的修改
    • 3.2.8 最終效果
  • 3.3 MSAA抗鋸齒
  • 3.4 其他模型
• 4. 一些可能遇到的問題
• • 4.1 渲染耗時過長？
  • 4.2 渲染結果中光源區域為純黑？
  • 4.3 渲染結果較暗？
  • 4.4 天花闆黑色，牆面沒有影子？
  • 4.5 多線程編譯時出現undefined reference to `pthread_create’？
  • 4.6 渲染出現白色噪點？
• 5. 效果集合
• 6. 附件

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前面，加上

它可以将跟在後面的for循環語句分成多個線程并發執行

然後在cmakelists.txt中加上

重新執行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);
           

同樣也要記住把鎖設定為全局變量：

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材質

然後在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】作業合集