Games101 Path Tracing笔记以及作业7代码实现

Path tracing笔记以及Games101作业7代码相关解读以及实现。笔记包括

1. 辐射度量学相关概念
2. 渲染方程
3. 蒙特卡洛积分
4. Path Tracing逻辑

辐射度量学 (Radiometry)

• 辐射通量 (Radiant flux (power)) $\Phi =\frac{\mathrm{d}Q}{\mathrm{d}t}$ Energy per unit time，表示单位时间内光穿过一个截面的光能。

• 辐射强度 (Radiant intensity)，$I(\omega )=\frac{\mathrm{d}\Phi }{\mathrm{d}\omega }$ Power per unit solid angle，表示单位立体角上的辐射通量。

• 辐射率 (Radiance)，$L(\mathrm{p},\omega )=\frac{{\mathrm{d}}^2\Phi (p,\omega )}{\mathrm{d}\omega \cdot \mathrm{d}A\cos \theta }$ Power emitted, reflected, transmitted or received by a surface, per unit solid angle, per projected unit area，表示单位立体角，单位投影面积上的辐射通量。这里多出来一个$\cos \theta$是因为Lamber’s Cosine Law，即阳光在斜照的时候强度会由于倾角变成原来的$\cos \theta$倍。

• 辐照度 (Irradiance)，$E(\mathrm{x})=\frac{\mathrm{d}\Phi (\mathrm{x})}{\mathrm{d}A}$ Power per unit area incident on a surface point，表示单位面积上的辐射通量。这里有一种记忆方法，在Games104上学到的，将Irradiance拆分成Ir(In)和Radiance记为进入到表面的光照强度。

• 立体角 (Solid Angles)

  

Bidirectional Reflectance Distribution Function (BRDF)

• 渲染方程 (Rendering Equation)
  $$L_o(p,{\omega }_0)=L_e(p,{\omega }_0)+{\int }_{\Omega +}{L}_i(p,{\omega }_i)f_r(p,{\omega }_i,{\omega }_0)(n\cdot {\omega }_i)\mathrm{d}{\omega }_i$$
  其中第一项表示物体本身的发出的光照强度，第二项的积分所包括的内容就是我们的反射函数 (Reflection Function)，即
  $$L_r(\mathrm{p},{\omega }_r)={\int }_{\Omega +}{L}_i(p,{\omega }_i)f_r(p,{\omega }_i,{\omega }_0)\cos {\theta }_i\mathrm{d}{\omega }_i$$
  其中第一项为入射方向的Radiance，代表从单位立体角方向入射到单位面积上的光照强度，第二项就是双向反射分布函数，他表示以下比值
  $$f_r({\omega }_i\to {\omega }_r)=\frac{\mathrm{d}{L}_r({\omega }_r)}{\mathrm{d}{E}_i({\omega }_i)}=\frac{\mathrm{d}{L}_r({\omega }_r)}{L_i({\omega }_i)\cos {\theta }_i\mathrm{d}{\omega }_i}[\frac{1}{\text{sr}}]$$
  也就是出射光Radiance与入射光Irradiance的比值，也就是代表了在各个方向上多少光照会被反射出去。这里的积分是对入射光的Radiance进行积分，得到的就是入射光在半球面上的Irradiance，然后和BRDF相乘得到的就是反射出的在单位立体角上的Radiance。

蒙特卡洛积分 (Monte Carlo Integration)

• Monte Carlo estimator $F_N=\frac{1}{N}{\sum }_{i=1}^N\frac{f(X_i)}{p(X_i)}$, Example, Basic (Uniform) Monte Carlo estimator $F_N=\frac{b-a}{N}{\sum }_{i=1}^{N}f(X_i)$
• Monte Carlo Integration $\int f(x)\mathrm{d}x=\frac{1}{N}{\sum }_{i=1}^N\frac{f(X_i)}{p(X_i)}$ The more samples, the less variance.

路径追踪 (Path Tracing)

Whitted-style 光线追踪始终采用镜面反射，在光线与diffuse表面相交时停止反射，这一假设会造成一些问题。首先，比如说诸如Glossy (感觉类似磨砂)的效果就无法表现出来。其次，由于不在diffuse表面进行反射，物体之间由于光线反射形成的现象比如说康奈尔Box中的墙面会映射出左右两面墙的红绿两种颜色(如下图)就无法表现出来。

Path tracing使用Rendering Equation进行光照的计算。在计算Rendering Equation中的积分部分时，我们就可以使用Monte Carlo Integration计算。一个简单的Monte Carlo的方法如下图所示 (Direct Illumination)，其中的pdf最直观的办法我们可以视为uniform的，也就是$\mathrm{pdf}=\frac{1}{2\pi }$。

以上是最简单的直接光照的蒙特卡洛积分，然后我们再加上递归的间接光照以及控制递归深度就完成了最终的Path Tracing。为了控制深度，我们使用一种叫做俄罗斯轮盘赌 (Russian Roulette)的方法。非常简单，这种方法就是在每次光线进行反射的时候进行一次随机选择，p的概率继续递归，(1-p)的概率停止，如果继续反射，则我们计算反射的光照强度为$\frac{L_r(p,{\omega }_r)}{p}$，这样我们最后根据数学期望得到的光照强度仍然是$L_r(p,{\omega }_r)$。使用后RR的Path Tracing伪代码如下所示。

最后还有一个问题需要解决，使用以上方法进行采样的效率实际上是不高的，因为很多光线实际上是反射不到光源上的，那样就会造成很多浪费。因此，我们可以使用另外一种采样方式，Sampling the Light，从光源处采样。为了实现光源采样，我们可以重写反射方程，将对入射角度$\omega$的积分转化为对光源面积$A$的积分，由于$\mathrm{d}\omega =\frac{\mathrm{d}A\cos {\theta }^{\prime }}{\Vert x^{\prime }-x{\Vert }^2}$ (由立体角的定义得出)，因此
$$\begin{array}{rl}{L}_r(\mathrm{p},{\omega }_r)& ={\int }_{\Omega +}{L}_i(x,{\omega }_i)f_r(x,{\omega }_i,{\omega }_0)\cos {\theta }_i\mathrm{d}{\omega }_i\\ & ={\int }_A{L}_i(x,{\omega }_i)f_r(p,{\omega }_i,{\omega }_0)\frac{\cos {\theta }_i\cos {\theta }_i^{\prime }}{\Vert x^{\prime }-x{\Vert }^2}\mathrm{d}{\omega }_i\end{array}$$

最终版的Path Tracing伪代码如下

作业7

首先，我们需要将上一次作业的几个函数的实现复制粘贴到这里来，包括Triangle::getIntersection() in Triangle.hpp，IntersectP(const Ray& ray, const Vector3f& invDir, const std::array<int, 3>& dirIsNeg) in Bounds3.hpp以及getIntersection(BVHBuildNode* node, const Ray ray) in BVH.cpp。代码如下所示

// Triangle.hpp
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;

    // find ray triangle intersection
    if (t_tmp < 0)
        return inter;
    inter.happened = true;
    inter.distance = t_tmp;
    inter.coords = ray(t_tmp);
    inter.normal = normal;
    inter.obj = this;
    inter.m = m;
    return inter;
}

// Bounds3.hpp
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
    // test if ray bound intersects
    Vector3f t_min = (this->pMin - ray.origin) * invDir;
    Vector3f t_max = (this->pMax - ray.origin) * invDir;
    if (!dirIsNeg[0]) std::swap(t_min.x, t_max.x);
    if (!dirIsNeg[1]) std::swap(t_min.y, t_max.y);
    if (!dirIsNeg[2]) std::swap(t_min.z, t_max.z);
    float t_enter = fmax(t_min.x, fmax(t_min.y, t_min.z));
    float t_exit = fmin(t_max.x, fmin(t_max.y, t_max.z));
    return t_enter <= t_exit && t_exit >= 0;
}

//BVH.cpp
Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    Intersection inter;
    if (node == nullptr)
        return inter;
    std::array<int, 3> dirIsNeg = {int(ray.direction.x > 0), int(ray.direction.y > 0), int(ray.direction.z > 0)};
    if (!node->bounds.IntersectP(ray, ray.direction_inv, dirIsNeg))
    {
        return inter; // No intersection in the current node
    }
    if (node->left == nullptr && node->right == nullptr)
    {
        if (node->object == nullptr)
            throw std::runtime_error("NULL Object Error");
        inter = node->object->getIntersection(ray);
        return inter;
    }
    Intersection left_inter = getIntersection(node->left, ray);
    Intersection right_inter = getIntersection(node->right, ray);
    return left_inter.distance < right_inter.distance ? left_inter : right_inter;
}

完成后确保编译运行都没有问题后，我们再正式开始这次作业的实现，我们只需要对castRay(const Ray ray, int depth) in Scene.cpp进行实现，代码如下所示

Vector3f Scene::castRay(const Ray &ray, int depth) const
{
    Vector3f L_dir = {0, 0, 0}, L_indir = {0, 0, 0};
    Intersection intersection = Scene::intersect(ray);
    if (!intersection.happened)
        return {};
    if (intersection.m->hasEmission())
        return intersection.m->getEmission();

    Intersection light_pos;
    float light_pdf = 0.0f;
    sampleLight(light_pos, light_pdf);
    Vector3f collision_light = light_pos.coords - intersection.coords;
    float dis = dotProduct(collision_light, collision_light);
    Vector3f collision_light_dir = collision_light.normalized();
    Ray light_to_object_ray(intersection.coords, collision_light_dir);
    Intersection light_ray_inter = Scene::intersect(light_to_object_ray);
    
    auto f_r = intersection.m -> eval(ray.direction, collision_light_dir, intersection.normal);
    
    // Pay attention to the precision here
    if (light_ray_inter.distance - collision_light.norm() > -0.005){
        L_dir = light_pos.emit * f_r * dotProduct(collision_light_dir, intersection.normal)
                * dotProduct(-collision_light_dir, light_pos.normal) / dis / light_pdf;
    }

    if (get_random_float() > RussianRoulette)
        return L_dir;

    Vector3f w0 = intersection.m -> sample(ray.direction, intersection.normal).normalized();
    Ray object_to_object_ray(intersection.coords, w0);
    Intersection islight = Scene::intersect(object_to_object_ray);
    if (islight.happened && !islight.m->hasEmission())
    {
        float pdf = intersection.m->pdf(ray.direction, w0, intersection.normal);
        f_r = intersection.m->eval(ray.direction, w0, intersection.normal);
        L_indir = castRay(object_to_object_ray, depth + 1) * f_r * dotProduct(w0, intersection.normal) / pdf / RussianRoulette;
    }
    return L_dir + L_indir;
}

最终生成的图片，这里为了效果好一点我跑了128的SPP，大概需要跑几个小时，如下图

未完待续

生成一张这样的图片几个小时确实有点太慢了，提高篇要求我们使用多线程来提高渲染的速度，并且采用微表面模型的材质。这也是十分重要的概念，不过还是没有深入理解，准备等到以后学习完Games202之后再回来看看如何实现。

参考资料

GAMES101:作业7_南酒猫的博客-CSDN博客_games101作业7