Path Tracing

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

• 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：在场景的所有光源上按面积 统一采样一个点，并计算该采样点的概率密度
• 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：给定一对入射、出射方向与法向量，计算这种情况下的 值

可能用到的变量有：

• RussianRoulette in Scene.cpp： P_RR

Path Tracing 伪代码：

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

Monte Carlo Integration

概率相关知识

设是一个随机变量，是任意实数，则称 为的累计分布函数（Cumulative Distribution Function，CDF）。如果对于随机变量 的累计分布函数 ，存在非负函数，使对任意实数，有

则称 为连续性随机变量，其中函数 称为 的概率密度函数，简称概率密度（Probability Distribution Function，PDF）。概率密度具有以下几个性质：

• 对于任意实数，有
• 若 在点 处连续，则有

设连续性随机变量 的概率密度函数为 ，若积分 绝对收敛，则称积分 的值为随机变量 的数学期望，记为 ，简称为期望：

设是一个随机变量，若存 在，则称它的值为 的方差，记为 ，即：

同时 的标准差或均方差为 ，记为，方差表示的是随机变量与其均值的偏移程度，随机变量的方差 可按下列公式计算：

几个定理

切比雪夫 不等式 ：设随机变量 具有数学期望 ，方差，则对于任意正数，不等式

成立。这一不等式称为切比雪夫不等式

辛钦大数定理 ：设 是相互独立，服从同一分布随机变量序列，且具有数学期望 ，作前 个变量的算术平均 ，则对于任意，有：

伯努利大数定理 ：设 是次独立重复试验中事件 A 发生的次数，是事件 A 在每次试验中发生的概率，则对于任意正数 ，有：

辛钦大数定理解释了：在大量重复试验下，样本的平均值约等式总体的平均值。伯努利大数定理解释了：在大量重复试验下，样本的频率收敛于其概率。

蒙特卡洛法积分

采用蒙特卡洛方法来计算函数积分，一般的定义为：设 是相互独立的样本且服从同一分布，概率密度函数表示为 ，则函数的积分可以表示为：

这就是蒙特卡洛法积分的一般等式

重要性采样

重要性采样 是已知被积函数的一些分布信息而采用的一种缩减方差的策略，还有别的策略像 俄罗斯轮盘切割 ， 分层采样 ， 拉丁超立方体采样 等，是通过控制采样的策略达到缩减方差的目的

考虑一个简单的例子，设一个函数为

采用蒙特卡洛法估计它的积分，选择区间 之间的均匀分布作为它的随机数，那么存在绝大部分的采样点在区间 之间，但是它对积分估计的贡献只有 0.01，小部分采样点在区间 之间，它对积分估计的贡献却非常的大，这种现象就会导致误差非常的大，简单提高采样数对估计量收敛的影响较小。那么，如果简单的选择一个采样策略：多采样区间 的样本点，少采样 的样本点，这就违背了蒙特卡洛法的本质，产生的统计结果就没有任何意义。蒙特卡洛法的核心，是根据某一概率分布来随机采样。

GGX（Trowbridge-Reitz）分布

GGX 即 Trowbridge-Reitz 分布，最初由 Trowbridge 和 Reitz^{[Trowbridge 1975]} 推导出，在 Blinn 1977 年的论文 ^{[Blinn 1977]}中也有推荐此分布函数，但一直没有受到图形学界的太多关注。30 多年后，Trowbridge-Reitz 分布被 Walter 等人独立重新发现 ^{[Walter 2007]}，并将其命名为 GGX 分布。之后，GGX 分布采用风潮开始在电影^{[Burley 2012]} 和游戏 ^{[Karis 2013]，[Lagarde 2014]} 行业中广泛传播，成为了如今游戏行业和电影行业中最常用的法线分布函数。

GGX 分布的公式为：

Path Tracing

A Simple Monte Carlo Solution

在直接光照的环境下，渲染下图中一个点

计算 点的在相机方向上的 radiance 有以下渲染方程

对于 蒙特卡洛 ：， 就是 ， 概率密度函数 可以是 （在半球面上进行均匀采样）。所以 点的 radiance 可以表示为：

Introducing Global Illumination

可以写出如下伪代码：

shade(p, wo)
    Randomly choose N directions wi~pdf
    Lo = 0.0
    For each wi
        Trace a ray r(p, wi)
        If ray r hit the light
            Lo += (1 / N) * L_i * f_r * cosine / pdf(wi)
        Else If ray r hit an object at q
            Lo += (1 / N) * shade(q, -wi) * f_r * cosine / pdf(wi)
    Return Lo

Problem 1: Explosion of #rays as #bounces go up

但是光线的数量随着弹射的次数增加，以指数的量级增长 #rays = N^#bounces，所以假设每一个着色点只有一条光线，随之而来的就是会存在大量的噪音，通过增加经过每个像素的 paths，再做radiance 的均值便可以减少噪音。

生成 paths 的伪代码如下：

ray_generation(camPos, pixel)
    Uniformly choose N sample positions within the pixel
    pixel_radiance = 0.0
    For each sample in the pixel
        Shoot a ray r(camPos, cam_to_sample)
        If ray r hit the scene at p
            pixel_radiance += 1 / N * shade(p, sample_to_cam)
    Return pixel_radiance

Problem 2: The recursive algorithm will never stop

可以通过 Russian Roulette（RR）（即以一个固定的概率 决定 路径追踪 是否可以继续——是否发从相机处射射线），来解决递归算法无法结束的问题。如果采用 RR 的方法需要更改每一次 路径追踪 后radiance的值 ，使得该点radiance 与实际相差不大。将 除以概率 ，即。那么可以求出 的期望，可以证明 RR 的合理性：

可以写出如下伪代码：

shade(p, wo)
    Manually specify a probability P_RR
    Randomly select ksi in a uniform dist. in [0, 1]
    If (ksi > P_RR)
        return 0.0
    
    Randomly choose N directions wi~pdf
    Lo = 0.0
    For each wi
        Trace a ray r(p, wi)
        If ray r hit the light
            Lo += (1 / N) * L_i * f_r * cosine / pdf(wi) / P_RR
        Else If ray r hit an object at q
            Lo += (1 / N) * shade(q, -wi) * f_r * cosine / pdf(wi) / P_RR
    Return Lo

Sampling the Light

光源离物体越远，需要从相机发射的射线越多（在半球上采样次数越多）才能到达光源，这样会造成性能上的问题，可以通过对光源的采样来解决在半球上采样次数越多的问题，即对于直接光照，原来是从 点所在的半球进行采样，变成对 所在的光源平面采样。

此时 概率密度函数 为（），但之前的渲染方程是对立体角进行积分，现在需要对面积进行积分。由上图和立体角的定义（Projected area on the unit sphere）可得

此时渲染方程也可以写为：

现在的渲染迭代流程主要考虑以下两个部分：

• 光源能够直接到达的点，通过对光源采样后计算 radiance 的（需要判断是否存在遮挡），同时不需要考虑 RR
• 递归的间接光照，需要考虑 RR

最终的伪代码如下：

shade(p, wo)

    # Contribution from the light source.
    Uniformly sample the light at x’ (pdf_light = 1 / A)
    Shoot a ray from p to x’
    If the ray is not blocked in the middle
        L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light
        
    # Contribution from other reflectors.
    L_indir = 0.0
    Test Russian Roulette with probability P_RR
    Uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)
    Trace a ray r(p, wi)
    If ray r hit a non-emitting object at q
        L_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RR
      
    Return L_dir + L_indir

相关代码

// Implementation of Path Tracing
Vector3f Scene::castRay(const Ray &ray, int depth) const
{
    // TO DO Implement Path Tracing Algorithm here
    Vector3f ldir = {0, 0, 0 };
    Vector3f lindir = {0, 0, 0 };

    Intersection objectInter = intersect(ray);
    if (!objectInter.happened)
    {
        return {};
    }   

    if (objectInter.m->hasEmission())
    {
        return objectInter.m->getEmission();
    }

    Intersection lightInter;
    float lightPdf = 0.0f;
    sampleLight(lightInter, lightPdf);

    Vector3f obj2light = lightInter.coords - objectInter.coords;
    Vector3f obj2lightdir = obj2light.normalized();
    float distancePow2 = obj2light.x * obj2light.x + obj2light.y * obj2light.y + obj2light.z * obj2light.z;

    Ray obj2lightray = {objectInter.coords, obj2lightdir};
    Intersection t = intersect(obj2lightray);
    if (t.distance - obj2light.norm() > -EPSILON)
    {
        ldir = lightInter.emit * objectInter.m->eval(ray.direction, obj2lightdir, objectInter.normal) * dotProduct(obj2lightdir, objectInter.normal) * dotProduct(-obj2lightdir, lightInter.normal) / distancePow2 / lightPdf;
    }

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

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

    return ldir + lindir;
}

相关函数

Scene::sampleLight→MeshTriangle::Sample→BVHAccel::Sample

最终结果是通过遍历 BVH 树到子节点，计算子节点光源 Object 的三角形面积，存在疑问？

Material::sample

• 首先在局部坐标系下的半球面上进行采样得到localRay，其中，

  服从上图分布，限制了上半球。

• 然后通过该点的法向量 将localRay转换成世界坐标系中的向量

  Coordinate system from a vector

  We can use the fact that the cross product gives a vector orthogonal to the two vectors to write a function that takes one vector and returns two new vectors so that the three of them form an orthonormal coordiante system. Specifically, all three of the vectors will be perpendicular to each other. Note that the other two vectors returned are only unique up to a rotation about the given vector. This function assumes that the vector passed in, , has already been normalized. We first construct a perpendicular vector by zeroing one of the two components of the original vector and permuting the remaining two. Inspection of the two cases should make clear that will be normalized and that the dot product will be equal to zero. Given these two perpendicular vectors, one more cross product wraps things up to give us the third, which by definition of the cross product will be be perpendicular to the first two^[1].

  inline void CoordinateSystem(const Vector &v1, Vector *v2, Vector *v3) {
      if (fabsf(v1.x) > fabsf(v1.y)) {
          Float invLen = 1.f / sqrtf(v1.x*v1.x + v1.z*v1.z);
          *v2 = Vector(-v1.z * invLen, 0.f, v1.x * invLen);
      }
      else {
          Float invLen = 1.f / sqrtf(v1.y*v1.y + v1.z*v1.z);
          *v2 = Vector(0.f, v1.z * invLen, -v1.y * invLen);
      }
      *v3 = Cross(v1, *v2);
  }

Scene::intersect

通过构建的 BVH 寻找射线与场景中物体的交点

Material::pdf

计算交点处的 概率密度函数，在半球内的返回常数0.5f / M_PI

Material::eval

计算交点处的漫反射系数，在半球内的返回常数Kd / M_PI

结果

---

• [1] Physically Based Image Synthesis: Design and Implementation of A Rendering Sysytem