闲言碎语

前言：从这篇文章开始，转知乎了，感觉知乎的社区氛围更好，emmm，哈哈哈主要还是想接触一下知乎的图形学大佬们。知乎账号链接：https://www.zhihu.com/people/xuu-27-24

emmm，距离完成作业6已经过去了一个月多的时间（我承认中间有段时间在摆烂==），怎么说呢，感觉作业7真的知识体系非常庞大，虽然实现的代码量不多，但是里面真的很多细节很多trick，对我这种小菜鸡感觉还是有挺大难度的。。。再加上感觉自己挺钻牛角尖的，所以搞了挺久的，也最终算是完成了作业7的全部要求，害，希望后面本科毕设不要太赶吧。

提个建议：如果遇到实在解决不了的问题，可以上games论坛看看，里面有很多共性的问题，上面的问题我几乎都遇到了。。。

最终全部效果展示（均为1024×1024×512ssp）

基础path tracing实现（作业基本要求）

微表面材质

完美镜面反射材质

课程总结与理解（Path Tracing）

课程方面，主要是先讲辐射度量学，然后推导出渲染方程，之后介绍用来求解积分的蒙特卡洛方法，最后讲如何在程序中求解渲染方程计算颜色（包括很多的trick细节，比如对光源进行直接采样提高采样效率，俄罗斯轮盘赌解决无限递归等等），从而实现path tracing。

path tracing流程：整体的渲染流程实际上更多的是在Whitted-Style上的改进，前期添加物体，构建包围盒，这些都是一样的。不一样的地方在于path tracing在最终着色的时候本质上是求解一个渲染方程，需要做大量的采样，近似积分结果，里面小trick非常多。而Whitted-Style并没有做采样求积分，就是简单的光线反射和折射，无法像path tracing一样考虑全局光照。 总体上来说，就是每个像素点发出多条光线，每条光线会通过之前构建的包围盒找到与场景中物体的第一个交点，然后根据渲染方程计算该交点的颜色，并返回，最后将全部发射的光线的颜色求平均，就得到该像素的颜色。

其实这个渲染方程本身的原理还是挺复杂的，可以研究的很深入，尤其辐射度量学那块，课程本身其实讲的并不是很细，如果想细致了解的话可以去翻虎书，虎书可以说对辐射度量学进行了深度的解刨，直接从光子开始讲。。。

框架梳理

和作业6使用的框架其实差不多，可以参考我之前写的https://blog.csdn.net/Xuuuuuuuuuuu/article/details/128556319

作业7的框架实际上主要就是castRay函数有所不同，因为path tracing的最终实现是在castRay中实现。

任务一：迁移相关代码

之前写的在这里：https://blog.csdn.net/Xuuuuuuuuuuu/article/details/128556319，把相关的代码贴进去就行。其中要注意的一个点，IntersectP中的等于号一定要去掉，要不然之后会出现部分物体不可见的情况，去掉的位置在这里，见下图：

为什么这里要去掉等于号，之前闫老师说，图形学编程中很少考虑等于号的情况，几乎是可以忽略的。这是因为这次的场景比较特殊，以右边的绿色墙壁为例，实际上它的Boundbox就是一个与某一轴平行长方形（不是长方体，不信的话可以printBoundbox的那两个点出来看看）所以光线和这个Boundbox相交检测的时候，算出来的t_enter和t_exit是相等的。此时如果相等时仍判断为不相交，就会出错，此时部分物体不可见，如下图：

任务二：实现path tracing

代码量也不是很大，主要就是跟着作业中给的伪代码照着敲，但是想呈现出最终的结果有很多小细节要注意，其中这里面各个向量的方向问题就很容易搞错，一定要仔细琢磨清楚。因为里面要修的细节太多了，我这边直接给出任务二完整的代码（个人建议不要完全照抄，可以自己先跟着作业提供的伪代码敲一下，遇到问题再去自己想办法找资料解决，在最终得到想要的结果后，你会发现这个过程非常值得。），具体的很多细节问题放到后文讨论。很多细节的问题建议直接看代码，可以理解的更清楚。

Vector3f Scene::castRay(const Ray& ray, int depth) const
{
    // TO DO Implement Path Tracing Algorithm here
    Vector3f hitColor = this->backgroundColor;
    Intersection shade_point_inter = Scene::intersect(ray);
    if (shade_point_inter.happened)
    {

        Vector3f p = shade_point_inter.coords;
        Vector3f wo = ray.direction;
        Vector3f N = shade_point_inter.normal;
        Vector3f L_dir(0), L_indir(0);

       //sampleLight(inter,pdf_light)
        Intersection light_point_inter;
        float pdf_light;
        sampleLight(light_point_inter, pdf_light);
        //Get x,ws,NN,emit from inter
        Vector3f x = light_point_inter.coords;
        Vector3f ws = normalize(x-p);
        Vector3f NN = light_point_inter.normal;
        Vector3f emit = light_point_inter.emit;
        float distance_pTox = (x - p).norm();
        //Shoot a ray from p to x
        Vector3f p_deviation = (dotProduct(ray.direction, N) < 0) ?
                p + N * EPSILON :
                p - N * EPSILON ;

        Ray ray_pTox(p_deviation, ws);
        //If the ray is not blocked in the middleff
        Intersection blocked_point_inter = Scene::intersect(ray_pTox);
        if (abs(distance_pTox - blocked_point_inter.distance < 0.01 ))
        {
            L_dir = emit * shade_point_inter.m->eval(wo, ws, N) * dotProduct(ws, N) * dotProduct(-ws, NN) / (distance_pTox * distance_pTox * pdf_light);
        }
        //Test Russian Roulette with probability RussianRouolette
        float ksi = get_random_float();
        if (ksi < RussianRoulette)
        {
            //wi=sample(wo,N)
            Vector3f wi = normalize(shade_point_inter.m->sample(wo, N));
            //Trace a ray r(p,wi)
            Ray ray_pTowi(p_deviation, wi);
            //If ray r hit a non-emitting object at q
            Intersection bounce_point_inter = Scene::intersect(ray_pTowi);
            if (bounce_point_inter.happened && !bounce_point_inter.m->hasEmission())
            {
                float pdf = shade_point_inter.m->pdf(wo, wi, N);
				if(pdf> EPSILON)
					L_indir = castRay(ray_pTowi, depth + 1) * shade_point_inter.m->eval(wo, wi, N) * dotProduct(wi, N) / (pdf *RussianRoulette);
            }
        }
        hitColor = shade_point_inter.m->getEmission() + L_dir + L_indir;
    }
    return hitColor;
}

应该能得到下图所示的结果（1024×1024×512ssp）：

任务三：多线程加速（包括其他加速的小trick）

我觉得加速是渲染中一个非常重要的点，也是我个人比较感兴趣的点。tmd没有加速之前，代码的运行速度简直惨不忍睹。。。一开始1024×1024×512ssp的图渲染了七八个小时才渲染了百分之20（渲染到百分之20我就直接放弃了==）。。。最终经过不懈努力，总算把速度优化到了40分钟，可以说快了几乎60-80倍。渲染速度的提高可以说极大地帮助了后面的调试！！！

按照下面三步逐步进行优化，我渲染一张512×512×8ssp的图的速度从13分钟提到了10秒钟。

1.随机数构造器优化

使用C++的性能分析器：https://blog.csdn.net/u011942101/article/details/123656944

不难发现，主要是求交和采样两部分花了很多的时间，求交这部分本来就耗时很长，里面很多递归。但是为啥，采样这部分要这么长时间？

最后发现是这个函数，特别耗时。一分析是因为构建开销过大，实际上把里面三个变量都定义成static就行（只会初始化一次，避免重复构建），没必要重新构建。经测试，改成static后，512×512×8ssp的渲染速度从13分钟提到了5分钟。

参考：https://games-cn.org/forums/topic/zuoyeqidexingnengpingjingshisuijishushengcheng/

2.多线程加速

使用多线程加速之前，问自己一个问题，为什么这里可以使用多线程加速？因为这里满足两个条件：1.这个场景下程序可以写成多线程，2.电脑多核。

实际上实现多线程加速不是很难，因为整个场景太适合写成多线程加速了。。。就每个线程分别独立一定数量计算像素颜色就行。原则上就是要最大化每个线程的计算量，对整个frame计算使用多线程。 具体使用过的是C++里面的thread进行多线程，代码如下：

// 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);
    int m = 0;

    // change the spp value to change sample ammount
    int spp = 512;//16
    int thread_num = 8;//我的电脑有8核，所以开8个线程。注：屏幕的高度一定要是线程数的倍数
    int thread_height = scene.height / thread_num;
    std::vector<std::thread> threads(thread_num);
    std::cout << "SPP: " << spp << "\n";
    
    //多线程实现
    std::mutex mtx;
    float process=0;
    float Reciprocal_Scene_height=1.f/ (float)scene.height;
    auto castRay = [&](int thread_index) 
    {
        int height = thread_height * (thread_index + 1);
        for (uint32_t j = height - thread_height; j < height; 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;

                //eye的位置对结果有影响
                Vector3f dir = normalize(Vector3f(-x, y, 1));
                for (int k = 0; k < spp; k++)
                {
                    framebuffer[j*scene.width+i] += scene.castRay(Ray(eye_pos, dir), 0) / spp;  
                }
            }
            mtx.lock();
            process = process + Reciprocal_Scene_height;
            UpdateProgress(process);
            mtx.unlock();
        }
    };

    for (int k = 0; k < thread_num; k++)
    {
        threads[k] = std::thread(castRay,k);
    }
    for (int k = 0; k < thread_num; k++)
    {
        threads[k].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);    
}

512×512×8ssp的渲染速度从5分钟提到了1分钟！！！
真的很担心cpu被烧坏。。。这种情况下说明多线程编程成功，CPU被疯狂地调用。

3.debug改成release版本

无心之举，512×512×8ssp的渲染速度从1分钟提到了10秒钟。。。 只能说release比debug版本快了真的不止一点半点。

任务四：微表面模型

101讲的不是很细，因此我就花了不少时间仔细地看了下202和一些其他资料。

参考1：https://zhuanlan.zhihu.com/p/434964126（写得贼nb，推导地很深入）

参考2：https://zhuanlan.zhihu.com/p/152226698（这个更加全面一点，涉及的更广，但1更深）

本质上微表面材质是描述高光项，真正的brdf是由漫反射项＋高光项组成的，如下（参考2）：

各自的比例是由菲涅尔项计算的，而菲涅尔项是折射和反射的比例，这里漫反射等价于折射（因为漫反射从微观角度是折射进入物体，然后多次反弹弹出表面，产生漫反射现象）。

说白了就是一部分算折射然后反弹出来的能量，一部分算反射出来的能量。

1. 微表面材质定义

说白了，就是把一个物体表面看作有多个微表面组成的表面，从而展示更加真实的材质细节。

下图的虚线就是物体表面，而对其着色的时候，微表面模型会将物体表面近似成下图实线尖尖样子进行着色，具体就是通过定义微表面模型的prdf实现的。

思想和凹凸贴图那边有点想，实际上没有改变物体的几何模型，但是通过改变了prdf，使物体有了多个微表面组成的感觉。所以有了这句话：远处看是材质和外观，近处看是几何。

以太空图为例，远处看一个物体的时候看不到很微小的东西，看到的是一个全局效果。

ok，上面主要大致讲微表面模型大概是个什么东西，下面就是具体深入的理解，

微表面模型主要分成两个部分（也就是说设计一个微表面模型就必须要考虑到这两部分）：

①微表面分布的表达（全局，比如法线分布、几何遮挡等）

②光在每个微表面散射的BRDF（细节，单个微表面的材质，即prdf，可以是完美镜面反射，又可以是理想漫反射）

接下来依次介绍以上两部分，然后再介绍最常见的Torrance-Sparrow模型。

2. 微表面分布的表达

包括D(h)和G(i,o,h)两块，即法线分布和几何遮挡。

2.1. D(h)：法线分布函数（Normal Distribution Fuction，NDF）

不同的法线分布会产生完全不一样的材质。

D(h)本质上是面积微元dA上的法向为h的全部微表面的面积和与面积微元dA的比例。

注：dA是宏观表面上的面积微元。

不难发现，D(h)实际上是针对物体宏观表面上的一个面积微元的，而不是整个物体。

D(h)的推导公式。

①beckmann NDF

类似高斯分布函数，里面就是有那两个量，α应该是超参数，控制整个分布的宽瘦，θ就是h代进去计算得到的。为什么这么去设计这个函数，是有它的意义在里面。

②GGX NDF

相较于Beckmann NDF，有了个长尾巴，意味着会有更好的过渡。

2.1. G(i,o,h)：几何遮挡函数（Shadowing and Making term）

描述微表面之间互相遮挡的几何现象的数学模型，本质上是法向为h的微表面在w方向上可见的比例（因为有部分被遮挡了，见下面第一张图）。

遮挡项，因为光线或者视线非常平（第二张图的球边界，此时角度称为grazing angle），f的分母会非常小，因此f就会非常大，边界就会非常亮，此时就要引入遮挡项解决这个问题，进行衰减，因为这种情况下，是最容易出现shadowing-masking现象的。

具体推导可以看上述的链接，讲的非常清楚！！！

3.光在每个微表面散射的BRDF

这里又可以有很多文章可以做，可以把每个微表面看作理想镜面（完美反射），也可以看成理想玻璃表面（又有折射又有反射），还可以看成漫反射表面。

结合上述两个微表面的分布项，再定义每个微表面的材质，就能推导出相关的微表面模型。

下文以Torrance-Sparrow Model为例。

4.Torrance-Sparrow Model

把每个微表面看成理想镜面。

推导过程看不懂，我太菜了，就这样吧，该章节一开始的链接里有。

5.代码实现

实现的是上述的Torrance-Sparrow Model，F项是代码中已经写好的菲涅尔项，D项和G项选择GGX，如下图所示：

其中各参数的含义可以查看原文https://zhuanlan.zhihu.com/p/434964126，或者可以直接看代码，应该也能够理解。代码实现方面，先在定义里面添加Microfacet材质，如下：

enum MaterialType { DIFFUSE, MICROFACET}

然后在main函数中添加物体，并对其赋予Microfacet材质，如下：

    Material* microfacet = new Material(MICROFACET, Vector3f(0.0f));
    microfacet->Ks = Vector3f(0.45, 0.45, 0.45);
    microfacet->Kd = Vector3f(0.3, 0.3, 0.25);
    microfacet->ior = 12.85;
    Sphere sphere1(Vector3f(150, 100, 200), 100, microfacet);

    scene.Add(&sphere1);

sample和pdf直接照抄原本的就行。最后修改eval，在其中实现Torrance-Sparrow Model模型即可，如下：

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 MICROFACET:
        {
            float cosalpha = dotProduct(N, wo);
            if (cosalpha > 0.0f) 
            {
                // calculate the contribution of Microfacet model
                float F, G, D;

                fresnel(wi, N, ior, F);

                float Roughness = 1;//超参数，控制粗糙度
                auto G_function = [&](const float& Roughness, const Vector3f& wi, const Vector3f& wo, const Vector3f& N)
                {
                    float A_wi, A_wo;
                    A_wi = (-1 + sqrt(1 + Roughness * Roughness * pow(tan(acos(dotProduct(wi, N))), 2))) / 2;
                    A_wo = (-1 + sqrt(1 + Roughness * Roughness * pow(tan(acos(dotProduct(wo, N))), 2))) / 2;
                    float divisor = (1 + A_wi + A_wo);
                    if (divisor < 0.001)
                        return 1.f;
                    else
                        return 1.0f / divisor;
                };
                G = G_function(Roughness, -wi, wo, N);

                auto D_function = [&](const float& Roughness, const Vector3f& h, const Vector3f& N)
                {
                    float cos_sita = dotProduct(h, N);
                    float divisor = (M_PI * pow(1.0 + cos_sita * cos_sita * (Roughness * Roughness - 1), 2));
                    if (divisor < 0.001)
                        return 1.f;
                    else 
                        return (Roughness * Roughness) / divisor;
                };
                Vector3f h = normalize(-wi + wo);
                D = D_function(Roughness, h, N);

                // energy balance
                Vector3f diffuse = (Vector3f(1.0f) - F) * Kd / M_PI;
                Vector3f specular;
                float divisor= ((4 * (dotProduct(N, -wi)) * (dotProduct(N, wo))));
                if (divisor < 0.001)
                    specular= Vector3f(1);
                else
                    specular = F *G * D / divisor;
                //std::cout << "F:"<<F << "\n";
                //std::cout << "diffuse:"<<diffuse<<"\n";
                //std::cout << "specular:" << specular << "\n";
                return diffuse+specular;
            }
            else
                return Vector3f(0.0f);
            break;
        }
    }
}

最终效果（均为1024×1024×512ssp）

如果黑色噪点很多，尝试修改一下sphere的求交判断阈值，有可能是精度的问题，如下：

    Intersection getIntersection(Ray ray){
        Intersection result;
        result.happened = false;
        Vector3f L = ray.origin - center;
        float a = dotProduct(ray.direction, ray.direction);
        float b = 2 * dotProduct(ray.direction, L);
        float c = dotProduct(L, L) - radius2;
        float t0, t1;
        if (!solveQuadratic(a, b, c, t0, t1)) return result;
        if (t0 < 0) t0 = t1;
        if (t0 < 0) return result;
        if (t0 > 0.5)
        {
            result.happened = true;
            result.coords = Vector3f(ray.origin + ray.direction * t0);
            result.normal = normalize(Vector3f(result.coords - center));
            result.m = this->m;
            result.obj = this;
            result.distance = t0;
        }

        return result;

    }

如果白色噪点很多，在castRay函数的最后，限制一下hitColor的范围，避免出现特别大的值，影响采样质量。如下：

        hitColor = shade_point_inter.m->getEmission()+L_dir + L_indir;
        hitColor.x = (clamp(0, 1, hitColor.x));
        hitColor.y = (clamp(0, 1, hitColor.y));
        hitColor.z = (clamp(0, 1, hitColor.z));

任务五：完美反射模型（看到顺便实现了，作业没要求）

核心是重要性采样，参考这篇文章实现的：https://blog.csdn.net/ycrsw/article/details/124408789，各种细节这里面已经讲得非常清楚了，因此就不再赘述，直接给出代码，和Microfacet材质的实现流程差不多，首先添加材质，如下：

enum MaterialType { DIFFUSE, MICROFACET,MIRROR};

main函数中添加物体，赋予材质，

    Material* mirror = new Material(MIRROR, Vector3f(0.0f));
    mirror->Ks = Vector3f(0.45, 0.45, 0.45);
    mirror->Kd = Vector3f(0.3, 0.3, 0.25);
    mirror->ior = 12.85;
    
    MeshTriangle floor("../models/cornellbox/floor.obj", white);
    MeshTriangle shortbox("../models/cornellbox/shortbox.obj", white);
    MeshTriangle tallbox("../models/cornellbox/tallbox.obj", mirror);
    MeshTriangle left("../models/cornellbox/left.obj", red);
    MeshTriangle right("../models/cornellbox/right.obj", green);
    MeshTriangle light_("../models/cornellbox/light.obj", light);

最后修改eval，在eval中实现完美镜面反射（Mirror），如下：

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);
            if (cosalpha > 0.0f) 
            {
                float divisor = cosalpha;
                if (divisor < 0.001) return 0;
                Vector3f mirror = 1 / divisor;
                float F;
                fresnel(wi, N, ior, F);
                return F * mirror;
            }
            else
                return Vector3f(0.0f);
            break;
        }
        case MICROFACET:
        {
            float cosalpha = dotProduct(N, wo);
            if (cosalpha > 0.0f) 
            {
                // calculate the contribution of Microfacet model
                float F, G, D;

                fresnel(wi, N, ior, F);

                float Roughness = 1;//超参数，控制粗糙度
                auto G_function = [&](const float& Roughness, const Vector3f& wi, const Vector3f& wo, const Vector3f& N)
                {
                    float A_wi, A_wo;
                    A_wi = (-1 + sqrt(1 + Roughness * Roughness * pow(tan(acos(dotProduct(wi, N))), 2))) / 2;
                    A_wo = (-1 + sqrt(1 + Roughness * Roughness * pow(tan(acos(dotProduct(wo, N))), 2))) / 2;
                    float divisor = (1 + A_wi + A_wo);
                    if (divisor < 0.001)
                        return 1.f;
                    else
                        return 1.0f / divisor;
                };
                G = G_function(Roughness, -wi, wo, N);

                auto D_function = [&](const float& Roughness, const Vector3f& h, const Vector3f& N)
                {
                    float cos_sita = dotProduct(h, N);
                    float divisor = (M_PI * pow(1.0 + cos_sita * cos_sita * (Roughness * Roughness - 1), 2));
                    if (divisor < 0.001)
                        return 1.f;
                    else 
                        return (Roughness * Roughness) / divisor;
                };
                Vector3f h = normalize(-wi + wo);
                D = D_function(Roughness, h, N);

                // energy balance
                Vector3f diffuse = (Vector3f(1.0f) - F) * Kd / M_PI;
                Vector3f specular;
                float divisor= ((4 * (dotProduct(N, -wi)) * (dotProduct(N, wo))));
                if (divisor < 0.001)
                    specular= Vector3f(1);
                else
                    specular = F *G * D / divisor;
                //std::cout << "F:"<<F << "\n";
                //std::cout << "diffuse:"<<diffuse<<"\n";
                //std::cout << "specular:" << specular << "\n";
                return diffuse+specular;
            }
            else
                return Vector3f(0.0f);
            break;
        }
    }
}

哦对，由于要重要性采样（采样方向只需要采完美反射的方向即可），所以采样的那两个函数也需要改变，不能像Microfacet一样不变。如下：

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;
        }
        case MICROFACET:
        {
            // 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;
        }
    }
}

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) > 0.0f)
                return 1.0f;
            else
                return 0.0f;
            break;
        }
        case MICROFACET:
        {
            // uniform sample probability 1 / (2 * PI)
            if (dotProduct(wo, N) > 0.0f)
                return 0.5f / M_PI;
            else
                return 0.0f;
            break;
        }
    }
}

哦对了，最后的最后，因为要做重要性采样，castRay函数也要进行修改，否则会有过曝现象，mirror材质不再需要直接光L_dir采样了，就直接算L_indir就行，而且在算L_indir的时候记得那个把里面的避免光源采样的判断条件去了。如下：

Vector3f Scene::castRay(const Ray& ray, int depth) const
{
    // TO DO Implement Path Tracing Algorithm here
    Vector3f hitColor = this->backgroundColor;
    Intersection shade_point_inter = Scene::intersect(ray);
    if (shade_point_inter.happened)
    {
        Vector3f p = shade_point_inter.coords;
        Vector3f wo = ray.direction;
        Vector3f N = shade_point_inter.normal;
        Vector3f L_dir(0), L_indir(0);
        Vector3f p_deviation = (dotProduct(ray.direction, N) < 0) ?
            p + N * EPSILON :
            p - N * EPSILON;

        switch (shade_point_inter.m->getType())
        {
            case MIRROR:
            {
                //Test Russian Roulette with probability RussianRouolette
                float ksi = get_random_float();
                if (ksi < RussianRoulette)
                {
                    //wi=sample(wo,N)
                    Vector3f wi = normalize(shade_point_inter.m->sample(wo, N));
                    //Trace a ray r(p,wi)
                    Ray ray_pTowi(p_deviation, wi);
                    //If ray r hit a object at q
                    Intersection bounce_point_inter = Scene::intersect(ray_pTowi);
                    if (bounce_point_inter.happened)
                    {
                        float pdf = shade_point_inter.m->pdf(wo, wi, N);
                        if (pdf > EPSILON)
                            L_indir = castRay(ray_pTowi, depth + 1) * shade_point_inter.m->eval(wo, wi, N) * dotProduct(wi, N) / (pdf * RussianRoulette);
                    }
                }
                break;
            }
            default:
            {
                //sampleLight(inter,pdf_light)
                Intersection light_point_inter;
                float pdf_light;
                sampleLight(light_point_inter, pdf_light);
                //Get x,ws,NN,emit from inter
                Vector3f x = light_point_inter.coords;
                Vector3f ws = normalize(x - p);
                Vector3f NN = light_point_inter.normal;
                Vector3f emit = light_point_inter.emit;
                float distance_pTox = (x - p).norm();
                //Shoot a ray from p to x
                Ray ray_pTox(p_deviation, ws);
                //If the ray is not blocked in the middleff
                Intersection blocked_point_inter = Scene::intersect(ray_pTox);
                if (abs(distance_pTox - blocked_point_inter.distance < 0.01))
                {
                    L_dir = emit * shade_point_inter.m->eval(wo, ws, N) * dotProduct(ws, N) * dotProduct(-ws, NN) / (distance_pTox * distance_pTox * pdf_light);
                }
                //Test Russian Roulette with probability RussianRouolette
                float ksi = get_random_float();
                if (ksi < RussianRoulette)
                {
                    //wi=sample(wo,N)
                    Vector3f wi = normalize(shade_point_inter.m->sample(wo, N));
                    //Trace a ray r(p,wi)
                    Ray ray_pTowi(p_deviation, wi);
                    //If ray r hit a non-emitting object at q
                    Intersection bounce_point_inter = Scene::intersect(ray_pTowi);
                    if (bounce_point_inter.happened && !bounce_point_inter.m->hasEmission())
                    {
                        float pdf = shade_point_inter.m->pdf(wo, wi, N);
                        if (pdf > EPSILON)
                            L_indir = castRay(ray_pTowi, depth + 1) * shade_point_inter.m->eval(wo, wi, N) * dotProduct(wi, N) / (pdf * RussianRoulette);
                    }
                }
                break;
            }
        }
        hitColor = shade_point_inter.m->getEmission()+L_dir + L_indir;
        hitColor.x = (clamp(0, 1, hitColor.x));
        hitColor.y = (clamp(0, 1, hitColor.y));
        hitColor.z = (clamp(0, 1, hitColor.z));

    }
    return hitColor;
}

最终结果（均为1024×1024×512）

这时候有人会问了，tmd，重要性采样怎么这么麻烦，那我不用，把ssp调高点不就行了？emmm，如果头铁，不做重要性采样，我也实验了一下，就是下图的结果，ssp已经是2048了，但是噪点还是非常多，就是采样效率太低了，做了重要性采样后，ssp为512的结果可以完美爆杀这张图。

代码实现可能遇到的问题

1.光源融入到整个天花板，没有出现黑色

光源融入到整个天花板了，==tmd，感觉网上没有人和我这个情况是一样的。。。

光源尼玛的，就是融入整个天花板了，==tmd，感觉网上没有人和我这个情况是一样的。。。

测试的时候发现，其实这个光源是能被检测出来的，只是着色的时候变成这样了，然后分析了好几个小时的源码，艹，一直找不到原因，最后debug了半天，总算知道为什么了。

就是当对光源着色时，目前是分为两个部分，一个是直接光，一个是间接光。

直接光的话应该是不会提供颜色的，因为直接光同样是在光源上采样，相当于ws和N是几乎垂直的，cos(ws,N)项就会使得整个直接光的值接近0。代码中点偏移的处理只会让两者夹角大于90，此时cos算出来甚至还是小于0的。综上，直接光不会给其提供颜色。

间接光会对其提供颜色，并且也是导致光源颜色和天花板颜色接近的直接原因，因为计算间接光时，光源本身是 漫反射材质，所以这边代码会让光源应该会有light->Kd差不多的颜色，而这个值又和天花板的颜色其实差不多，所以看上去才只能看到天花板。

这个地方我看其他博主的解决方法都很简单粗暴，检测出是光源就直接返回光源本身颜色，但实际上光源也是会吸收反射其他光的，如果场景中有其他光源，那这种处理有可能就会出问题。**因此就需要保证光源是可以正常接收直接光和间接光的，不能直接返回光源本身颜色。**此时观察之前的公式，少了一个自身发光项（下面第二张图，之前一直不知道这个有啥用，此时起作用了，醍醐灌顶），最后把这个加上去就可以正常显示了。

emm，之前好几个小时一直纠结在直接光上面，想方设法要避免直接光，其实并不是这个原因导致的，害==，也是尼玛的运气好，多实验了一下，结果和想的不一样，就再思考分析了下，就想出来原因了。。。

也明白了个道理，发现问题并打算编写代码去解决的时候，先看看能不能设计实验用代码验证一下这个问题的提出是不是正确的，如果这个问题本身就是错误的，那就没有必要花很多时间去写代码解决。我这次就花了很多时间去思考怎么避免直接光，实际上这个问题本身就是不成立的==我应该试一下没有直接光的话，光源区域会不会变黑，如果变黑了，就说明我之前提出的问题是成立的。


我草，事后去games论坛上，居然真的找到了个一样的了==
https://games-cn.org/forums/topic/zuoye7wenti/

2.箱子背阴处有很多黑色噪点

没对好伪代码中的那些向量方向定义，说白了就是没完全按照伪代码来，改了之后就好了。

修改后的结果

3.有白色噪点（ssp调高仍然有）

画面上有少部分白色噪点。

用ps工具放大，发现是某个像素为白色，此时SPP已经是512，所以应该不是采样不够的问题。而且如果是采样不够，应该是会在其周围同时出现大量噪点，而不是只出现这么一个这么突出的，例如下面第三张图。


这意味着计算光照时，出现了一个非常大的值。排查代码发现原因很可能在于除数过小（无限接近于0，或等于0，这里就是等于0），结果接近无穷大。使这个像素的全部采样都收到这个无限大值的影响，直接变成白色。这显然是我们不想要的。这就对应作业里的提示：pdf接近0。解决方案就是pdf低于某个指定阈值时，L_indir为0，代码如下：

float pdf = shade_point_inter.m->pdf(wo, wi, N);
if(pdf> EPSILON)
	L_indir = castRay(ray_pTowi, depth + 1) * shade_point_inter.m->eval(wo, wi, N) * dotProduct(wi, N) / (pdf *RussianRoulette);

修改后的结果

通过这个案例，分析一下为什么要考虑数值精度，emmm，如果不考虑数值精度，很有可能一次数值的不正确就会产生难以想象的结果，例如上述除数为0导致无限大inf的出现，会使得当次像素颜色计算的其他采样结果无效，相当于一个极偏值对整体结果产生了巨大影响，而我们想要的是能够避免极偏值的影响，类似机器学习里面的回归分类问题。所以这个时候就要对其特殊处理，降低他对其他结果的影响。

理论思考

1.偏移为什么要分两种情况这么偏移？

考虑了背光情况（因为世间万物都有厚度，如果光在前面，从背面看就是会被挡住），如下面第二张图。

2.为什么path tracing有软阴影？

阴影部分应该都是没有直接光的，那么导致软阴影的原因就是环境光，接受的间接光是不一样的，离得近的会被挡住更多的间接光，离得远的间接光挡住的比较少，所以就会有软阴影。

而witted style就没有软阴影，这是因为没有考虑全局光，就是只考虑了那一个方向。

3.把witted style和path tracing的结果对比一下，两者的实现的不同点有哪些？

区别很大，path tracing是基于witted style的改进，很容易看到这两者共同的地方。明显，path tracing是合理的，并且可以基于物理，还能考虑全局光照，能量守恒保证结果更加精确。具体的可以看 Ray Tracing笔记那一章。

反正我觉得两者最大的区别，就是path tracing通过渲染方程考虑了全局光照，而witted style无法考虑全局光照（没有全局积分，一个点的颜色就由那么几根光线决定，实现那种完美镜面效果倒是好的很。。。），它因此很难实现很多的细节，大部分材质表现也没有path tracing好。

4.为什么程序中，渲染方程的入射radiance是光源的emit，是一个常量？

emmm，如果要深究，非常复杂，后面的光学理论看得我头痛，大致就下面这么理解就够了。

反正radiance在程序中就是一根光线。

5.为什么说光栅化比光追快？

光栅化没有大量的采样，光追里面因为要求解积分，需要大量的采样。

个人感悟

emmm，可以说这是我从小到大，除一些项目以外，做某门课作业花的最长的时间了，很感谢闫老师的101课程，提供了优质的作业，通过作业7对光线追踪技术有了一个基础的基本认识，也在做的过程中不断发现问题，解决问题，提高自己的coding和思考能力，真的很感谢闫老师，我有一种预感，这门课将会为中国图形学持续输送数以万计的人才。这次作业做完后差不多打好了离线渲染的一些基本的底子，打算明天开始正式做本科的毕设了（之前10月到12月一直在摆烂，1月到2月在补渲染的基础知识，所以说本科毕设正式开始时间就差不多是现在2月，希望来得及==害），加油加油！

也希望这个博客可以对做这个作业的人有一定帮助！

参考资料

1.https://blog.csdn.net/u011942101/article/details/123656944（代码性能分析器）

2.https://blog.csdn.net/ycrsw/article/details/124408789

3.https://blog.csdn.net/qq_41765657/article/details/121942469

4.https://blog.csdn.net/ycrsw/article/details/124565054

5.https://blog.csdn.net/weixin_44491423/article/details/127552276

6.https://games-cn.org/forums/topic/guanyuguangxianzhuizong3zhongradiancedeshizi/

7.https://www.zhihu.com/question/28476602/answer/41003204

8.https://games-cn.org/forums/topic/guanyuxuanranfangchengdiguidingyideyiwen/

9.https://www.jianshu.com/p/0cfc3204af77

10.https://zhuanlan.zhihu.com/p/434964126

11.https://zhuanlan.zhihu.com/p/152226698