Physically Based Rendering

Ray differentials 包含了射线的相关的信息，其实就是在当前像素点上分别在 x 和 y 方向上偏移一像素生成的射线.

Li 函数用于计算给定射线的入射辐射率

• ray 指定被计算的射线
  
• scene 指定被渲染的场景
  
• sampler 样本生成器，使用 MonteCarlo 积分来计算光照传输方程
  
• arena MemoryArena 对象，用于积分器高效地分配临时内存空间
  
• depth 射线从摄像机发出，直到当前调用 Li()，弹射的次数
  

- 依据胶片大小划分tile, 每个tile并行进行渲染
  - 为tile分配内存
  - 为 tile 克隆sampler实例
  - 计算tile的采样边界
  - 依据tile的采样边界获取对应的FileTile
  - 遍历tile内的每个像素，对每个像素进行渲染
  	- 对每个像素的每个样本进行处理
      - 根据像素坐标获取对应的CameraSample
      - 为当前的CameraSample 生成摄像机射线
			- 计算沿着摄像机射线的辐射率 L
      - 记录当前样本贡献的辐射率，以及对应的权重

参考项目 unitycatlikecoding\RenderingAndAdvancedR\Assets\MyTest\08BoundingBoxAlgorithm 中计算。

template <typename T>
inline bool Bounds3<T>::IntersectP(const Ray &ray, Float *hitt0,
                                   Float *hitt1) const {
    Float t0 = 0, t1 = ray.tMax;
    for (int i = 0; i < 3; ++i) {
        // Update interval for _i_th bounding box slab
        Float invRayDir = 1 / ray.d[i];
        Float tNear = (pMin[i] - ray.o[i]) * invRayDir;
        Float tFar = (pMax[i] - ray.o[i]) * invRayDir;

        // Update parametric interval from slab intersection $t$ values
        if (tNear > tFar) std::swap(tNear, tFar);

        // Update _tFar_ to ensure robust ray--bounds intersection
        tFar *= 1 + 2 * gamma(3);
        t0 = tNear > t0 ? tNear : t0;
        // tFar为无穷大时，t1=t1=ray.tMax
        t1 = tFar < t1 ? tFar : t1;
        if (t0 > t1) return false;
    }
    if (hitt0) *hitt0 = t0;
    if (hitt1) *hitt1 = t1;
    return true;
}

template <typename T>
inline bool Bounds3<T>::IntersectP(const Ray &ray, const Vector3f &invDir, const int dirIsNeg[3]) const
{
    const Bounds3f &bounds = *this;
    // Check for ray intersection against $x$ and $y$ slabs
    // 求xAxis yAxis方向ray和BoundBox的相交区域
    Float tMin = (bounds[dirIsNeg[0]].x - ray.o.x) * invDir.x;
    Float tMax = (bounds[1 - dirIsNeg[0]].x - ray.o.x) * invDir.x;
    Float tyMin = (bounds[dirIsNeg[1]].y - ray.o.y) * invDir.y;
    Float tyMax = (bounds[1 - dirIsNeg[1]].y - ray.o.y) * invDir.y;

    // Update _tMax_ and _tyMax_ to ensure robust bounds intersection
    tMax *= 1 + 2 * gamma(3);
    tyMax *= 1 + 2 * gamma(3);
    // xAxis方向上的相交区域向yAxis方向延伸 yAxis方向上的相交区域向xAxis方向延伸 如果两个延伸后的区域不相交 说明ray和BoundingBox不相交
    if (tMin > tyMax || tyMin > tMax) return false;
    if (tyMin > tMin) tMin = tyMin;
    if (tyMax < tMax) tMax = tyMax;

    // Check for ray intersection against $z$ slab
    Float tzMin = (bounds[dirIsNeg[2]].z - ray.o.z) * invDir.z;
    Float tzMax = (bounds[1 - dirIsNeg[2]].z - ray.o.z) * invDir.z;

    // Update _tzMax_ to ensure robust bounds intersection
    tzMax *= 1 + 2 * gamma(3);
    if (tMin > tzMax || tzMin > tMax) return false;
    if (tzMin > tMin) tMin = tzMin;
    if (tzMax < tMax) tMax = tzMax;
    return (tMin < ray.tMax) && (tMax > 0);
}

下面文章对算法解释的很清楚：

• 射线与包围盒的相交测试 https://zhuanlan.zhihu.com/p/138259656
  

n 阶的贝塞尔曲线，是两个 n-1 阶贝塞尔曲线之间的插值

下图演示了从一阶 Bezier 曲线导出二阶 Bezier 曲线:

./PhysicallyBasedRendering/00_03_bezier_curve.ggb

将一条 BezierCurve 分为多段，可以通过 Blossoming 来确定每一段 BezierCurve 的控制点

// p[4] 为BezierCurve的4个控制点
static Point3f BlossomBezier(const Point3f p[4], Float u0, Float u1, Float u2)
{
    Point3f a[3] = {Lerp(u0, p[0], p[1]), Lerp(u0, p[1], p[2]),
                    Lerp(u0, p[2], p[3])};
    Point3f b[2] = {Lerp(u1, a[0], a[1]), Lerp(u1, a[1], a[2])};
    return Lerp(u2, b[0], b[1]);
}

// 下面代码求出了当前子段Curve对应的控制点 TODO Why?
Point3f cpObj[4];
cpObj[0] = BlossomBezier(common->cpObj, uMin, uMin, uMin);
cpObj[1] = BlossomBezier(common->cpObj, uMin, uMin, uMax);
cpObj[2] = BlossomBezier(common->cpObj, uMin, uMax, uMax);
cpObj[3] = BlossomBezier(common->cpObj, uMax, uMax, uMax);

• 贝塞尔曲线 https://zh.wikipedia.org/wiki/%E8%B2%9D%E8%8C%B2%E6%9B%B2%E7%B7%9A
  
• Bezier Blossoms https://courses.engr.illinois.edu/cs418/fa2018/slides/10-5%20Bezier%20Blossoms.pdf
  
• Blossom https://en.wikipedia.org/wiki/Blossom_(functional)
  
• A Primer on Bézier Curves https://pomax.github.io/bezierinfo/
  

• Subdivision Surface 比多边形 Mesh 要平滑
  
• 大多数现存的基础设施建模系统可以被重定向为 Subdivision，经典的多边形建模工具可用于构建 Subdivision 控制模型
  
• Subdivision 方法通常是基于样条的表面表示的一般化。
  
• 通过在局部增加控制表面就可以很容易地为局部区域增加细节，这对于样条表现形式则比较困难
  

SDFace 表示 subdivision mesh 的三角面
SDVertex 表示 subdivision mesh 的顶点

subdivision mesh 可以是封闭的 mesh 也可以是开放的 mesh

• 封闭的 mesh 没有边界，所有的面在每个边上都有邻接面，即每个面都有 3 个邻接面
  
• 开放的 mesh 有边界
  

在 mesh 内部的，大部分顶点有 6 个邻接面，有 6 个邻接顶点
在开放 mesh 的边界上，大多数顶点有 3 个邻接面，有 4 个邻接顶点

顶点的邻接顶点数量被称为该顶点的 valence
valence 值大于 6 的内部顶点，和 valence 值大于 4 的边界顶点被称为 Extraordinary Vertices；否则为 Regular Vertices.

Loop Subdivision 表面在所有地方都是平滑的，除了在 Extraordinary Vertices 处。

原来在 Mesh 中已经存在的点被称为 Even Vertices（SDVertex.child 指向的 vertex），分割边的新的顶点被称为 Odd Vertices(SDFace.children[].v[]数组中的某个 vertex);

下图表示了 SDFace 中顶点和邻接面数据的拓扑结构

下图表示了判断顶点是否在边界上

对控制 Mesh 进行 Subdivision 后，每个 SDVertex 有一个 children SDVertex，每个 SDFace 有 4 个 children SDFace

void SDVertex::oneRing(Point3f *p) {
    if (!boundary) {
        // Get one-ring vertices for interior vertex
        SDFace *face = startFace;
        do {
            *p++ = face->nextVert(this)->p;
            face = face->nextFace(this);
        } while (face != startFace);
    } else {
        // Get one-ring vertices for boundary vertex
        SDFace *face = startFace, *f2;
        while ((f2 = face->nextFace(this)) != nullptr) face = f2;
        *p++ = face->nextVert(this)->p;
        do {
            *p++ = face->prevVert(this)->p;
            face = face->prevFace(this);
        } while (face != nullptr);
    }
}

下图展示了计算非边界顶点周围一圈顶点的过程：

无法在有限的内存中表示无限的实数。

在 Ray-BoundingBox 相交计算中，我们沿着 Ray 方向找到进入 BoundingBox 的交点处的 t_min，然后找到离开 BoundingBox 的交点处的 t_max，如果 t_min<t_max说明有交点，否则说明无交点。
使用浮点数计算，可能会导致错误，使得相交测试失败。

三角形和射线相交计算中，将三角形转化到求交坐标系下，使用重心坐标判断射线和三角形是否相交，避免了上面计算方式导致的误差。尽管错误可以从三角形转化到求交坐标系时产生，但是该错误不会影响相交测试，因为相同的变换被应用于所有的三角形(这个误差非常小，不会明显影响交点计算的精确性)。

该算法的关键在于使用浮点算术运算，edge 函数可以保证有正确的符号。

下图分析了按照普通方法求 Ray 和 Shape 交点的误差：

下图展示了使用投影法求交点：

下图展示了三角形和射线交点的求解：

下图展示了其他 Shape 交点求解：

下图展示了矩阵变换引入的误差：

下图展示了 SpawnRay 时，Ray Origin 的计算

  inline Ray Transform::operator()(const Ray &r, Vector3f *oError, Vector3f *dError) const
  {
      Point3f o = (*this)(r.o, oError);
      Vector3f d = (*this)(r.d, dError);
      Float tMax = r.tMax;
      Float lengthSquared = d.LengthSquared();
      if (lengthSquared > 0) {
          Float dt = Dot(Abs(d), *oError) / lengthSquared;
          o += d * dt;
          //        tMax -= dt;
      }
      return Ray(o, d, tMax, r.time, r.medium);
  }

下图展示了 RayBoundingBox 交点计算中，舍入操作没有影响到交点计算精确性：

下图展示了 Ray Triangle 计算中，如何避免在 Ray 原点后产生交点：

下图展示了划分 primitive 时，如何选择划分的方向：

Surface Area Heuristic 提供了有充分依据的消耗模型，可以回答以下这类问题：

1. 将几何单元划分为多少可以得到一个更好的 BVH 用于射线交点测试？
   
2. 将空间划分为多少子空间可以得到一个更好的加速结构？
   

射线通过一个区域，执行相交测试需要的消耗为下面公式 1；
另一种情况是一个区域被划分为了两个区域，则相交测试需要的消耗为下面公式 2；

t_trav 表示遍历 interior 节点，确定射线经过了哪个子区域所消耗的时间。
pA 和 pB 分别为射线通过 A 区域和 B 区域的概率
Na 和 Nb 分别为 A 区域和 B 区域中的几何单元数量

下图为计算射线通过区域的概率的原理：

不是考虑所有可能的划分方法，然后计算每种划分的消耗，最后选择消耗最少的划分。这样的方式消耗太大了。
沿划分的轴，将整个区域划分为多个等距的桶，只需要考虑在桶边界上的划分。这种方式更加高效，而且依然可以获得有效的划分。下图展示了这种方式：

在构建 Bounding Volume Hierachies 时使用 Surface Area Heuristic 方法可以给出很好的结果，但是，这种方式有两个缺点：

1. 在树的每一层都需要对 SAH 消耗进行计算。
   
2. 从上至下的构建 BVH 对于并行非常不友好。最明显的并行方案(并行构建子树)受制于只有上一级树构建完成后才能构建下一级。
   

LBVH 就是用来解决以上问题的。
LBVH 背后关键的理念是将 BVH 构建转化为排序问题。

• Introduction to Acceleration Structure https://www.scratchapixel.com/lessons/advanced-rendering/introduction-acceleration-structure/introduction
  

将树的节点所占用内存大小压缩到 8 字节，可以保证 8 个节点刚好可以放到 64 字节的 cache（8 字节的 Node 比 16 字节的 Node 速度提高了 20%）。

  struct KdAccelNode
  {
      union {
          // 表示分割位置
          Float split;                 // Interior
          // 如果只有一个图元，表示几何体的偏移
          int onePrimitive;            // Leaf
          // 如果由多个图元，表示第一个几何体的偏移
          int primitiveIndicesOffset;  // Leaf
      };

  private:
      union {
          // Leaf     第两位为11表示为Leaf   其他位表示几何体数量
          // Interior 第两位表示分割面所沿的轴  其他位表示分割后上面子节点的偏移
          int flags;       // Both
          int nPrims;      // Leaf
          int aboveChild;  // Interior
      };
  };

kd-trees 构建使用的 SAH 和 BVH trees 不同之处在于，前者在选择分割面时，如果该分割面生成的两个区域中，有一个区域是空区域则这种分割更优，因为空区域和射线求交消耗最小。
kd-trees 使用的 SAH 消耗模型如下：

kd-trees 将所有几何单元的 BoundingBox 的各个面做为分割面，选出消耗最小的分割面。通过将 BoundingBox 投影到各个轴上，可以方便地进行处理。
对 BoudingBox 在轴上的投影进行大小排序，小于分割面的图元为 below，大于分割面的图元为 above

下图展示了遍历树的操作：

光是电磁辐射

• 不同颜色的光对应不同的波长
  
• 经典电磁学里，光波的能量由振幅确定，与频率或波长无关。在量子力学的情形，单光子的能量由频率或波长确定，而光波是大量光子的集合，所以光波的能量或振幅对应的是光子的数量。
  

光子能量: E=hc/λ=hf=hv
h 为普朗克常量，f和 v 为光的频率，c为光速，λ为光的波长

可见光的波长范围为 380nm - 780nm
下图展示了电磁波的波长分布

energy, flux(power), intensity, irradiance, radiance 这些物理量可以用于描述电磁辐射。

SPD(spectral power distribution) 是一个波长函数，用于描述每种波长的光的总量。
Relative SPD 是每种波长的光的总量的相对值。

下图展示了常见光源的 SPD：

不同物体反射光的 SPD 不同，如下图展示了黄油、番茄、生菜反射光的 SPD：

• https://baike.baidu.com/item/%E5%85%89%E8%B0%B1%E5%88%86%E5%B8%83
  
• 光谱功率分布（spectral power distribution） http://mycolordoc.com/archives/741
  
• spectral power distribution https://en.wikipedia.org/wiki/Spectral_power_distribution
  
• 光的波长和振幅对应光子的什么属性？ https://www.zhihu.com/question/309647712
  
• 光子能量 https://baike.baidu.com/item/%E5%85%89%E5%AD%90%E8%83%BD%E9%87%8F
  

• Color and Light http://www.cse.psu.edu/~rtc12/CSE486/lecture26.pdf
  
• 色觉地图的建立（一）：光感受器、色匹配实验与 CIE RGB 坐标系 https://zhuanlan.zhihu.com/p/72803768
  
• HumanSensitivity to color http://graphics.stanford.edu/courses/cs178/applets/locus.html
  
• ColorMatching http://graphics.stanford.edu/courses/cs178/applets/colormatching.html
  

人类视觉系统的特性决定了可以只使用 3 个浮点值就可以表示人类看到的颜色。色彩感知的三原色理论表示，可以使用三个值 xλ,yλ,zλ为人类观察者精确表示所有可见光谱。
给定一个 SPD S(λ)，这三个值可通过下面公式计算得到，X(λ) Y(λ) Z(λ) 分别为光谱匹配曲线，这些曲线由 CIE 标准通过一系列主观试验测试绘制出来。这些匹配曲线和人眼视网膜上的三种颜色感知视锥细胞的响应曲线相似。
下图为光谱匹配曲线：

下图为由 SPD 和 X(λ) Y(λ) Z(λ)计算 xλ,yλ,zλ的方法：

需要注意的是，不同的光谱分布可能会得到相似的 xλ,yλ,zλ值。对于人类视觉来说，这些不同的 SPDs 确有相同的可视效果。这些光谱被称为条件等色。

尽管 XYZ 可以很好的将一个给定的 SPD 表示为人类观察者看到的颜色。但是其也有不足的地方，例如，计算两个颜色各自 XYZ 值的乘积所得的 XYZ 值和将两个 SPD 相乘然后计算结果的 XYZ 值有明显的差别。
XYZ 空间的 y 系数和 luminance 有紧密联系，其用于度量颜色的亮度感觉。

当我们将 RGB 颜色显示到显示器上时，实际上显示的光谱基本由 r、g、b 三个光谱响应曲线的权重和决定。下图显示了 LED 显示器和 LCD 显示器发出的 r、g、b 的分布：

下图显示了 RGB 颜色(0.6,0.3,0.2)在 LED、LCD 上的显示时的 SPD，可以看出，不同显示器显示相同颜色，所得的 SPD 却非常不同，因此用户提供的用于表示特定颜色的 RGB 值，只有知道显示器特性才有意义。

给定 SPD 的（xλ，yλ，zλ）表示形式，并选择一组特定的 SPD 来定义显示器的红色，绿色和蓝色，我们可以将（xλ，yλ，zλ）转换为对应显示器的 RGB 系数。给定显示器的光谱响应曲线 R(λ), G(λ), B(λ)，对于特定的显示器 RGB 系数可以通过下面方法来计算：

  inline void XYZToRGB(const Float xyz[3], Float rgb[3])
  {
      rgb[0] = 3.240479f * xyz[0] - 1.537150f * xyz[1] - 0.498535f * xyz[2];
      rgb[1] = -0.969256f * xyz[0] + 1.875991f * xyz[1] + 0.041556f * xyz[2];
      rgb[2] = 0.055648f * xyz[0] - 0.204043f * xyz[1] + 1.057311f * xyz[2];
  }

将 RGB 或 XYZ 值转化为 SPD 比较难，无数多个 SPD 有相同的 xλ，yλ，zλ系数。我们希望转化函数可以满足如下条件：

1. 当 RGB 系数相同时，得到的 SPD 应该是一个常数
   
2. 通常，希望计算得到的 SPD 是平滑的。大多数真实世界的物体都具有平滑的 SPD。不平滑的光谱主要来源于光源，特别是荧光灯。
   

通过将显示器的 R(λ), G(λ), and B(λ) SPDs 加权求和来构造 SPD 不是一个很好的解决方案，因为显示器的 R(λ), G(λ), and B(λ) SPDs 通常是不规律、不平滑的，加权求和得到的 SPD 也不会平滑。
Smits(1999)发现了一种将 RGBs 转化为 SPDs 的方法。分别为 red、green、blue 计算各自的 SPDs，并保证这些 SPDs 光滑，然后用给定的 RGB 系数对这些 SPDs 进行加权求和，从而得到平滑的 SPD，再将 SPD 转化为 RGB 可以得到一个和原始 RGB 系数接近的系数.可以对这种基本方式做下面两个额外的改进：

1. 不要使用 red\green\blue 的 SPDs 的和来表示常量光谱，直接通过常量光谱表示常量 SPDs。
   
2. 混合颜色是由两种元颜色混合而来，混合颜色最好通过预计算的光滑的 SPD 来表示，而不是通过元颜色的和表示
   

  SampledSpectrum SampledSpectrum::FromRGB(const Float rgb[3], SpectrumType type)
  {
      SampledSpectrum r;
      if (type == SpectrumType::Reflectance) {
          // Convert reflectance spectrum to RGB
          if (rgb[0] <= rgb[1] && rgb[0] <= rgb[2]) {
              // Compute reflectance _SampledSpectrum_ with _rgb[0]_ as minimum
              r += rgb[0] * rgbRefl2SpectWhite;
              if (rgb[1] <= rgb[2]) {
                  r += (rgb[1] - rgb[0]) * rgbRefl2SpectCyan;
                  r += (rgb[2] - rgb[1]) * rgbRefl2SpectBlue;
              } else {
                  r += (rgb[2] - rgb[0]) * rgbRefl2SpectCyan;
                  r += (rgb[1] - rgb[2]) * rgbRefl2SpectGreen;
              }
          }
      }
      // ...
  }

• LCD 显示器工作原理到底是怎样的？ https://www.zhihu.com/question/39452323
  
• 显示器如何显示准确的颜色？ https://zhuanlan.zhihu.com/p/43467096
  
• An RGB to Spectrum Conversion for Reflectances http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.40.9608&rep=rep1&type=pdf
  

下图为计算点 p 处的 irradiance(辐照度) n 为点 p 处表面的 normal，Li(p,w)为 radiance(入射辐射度)函数，cos(θ)为 dA⊥项。θ 为 w 和 n 的夹角。辐照度通常只计算 normal 指向的半球范围。

辐射量积分中的各种 cosine 项经常会分散其要表示的内容。使用投影立体角而不是立体角来度量积分对象覆盖的区域可以避免此问题。
物体的投影立体角是通过先将物体投影到单位球上，然后再将物体投影到垂直与 normal 的 disc 上得到的。如下图所示：

在半球范围内，从某个物体发出的总的通量可以通过下面公式计算，其中 A 为该物体的表面积：

将立体角积分转化为球面坐标系积分经常很方便。
下图将一个向量转化为球坐标系形式：

下图将立体角转化为球面坐标系形式：

则点 p 处的 irradiance(辐照度)可以通过如下方式计算：

将在方向上的积分转化为在区域(面积)上的积分可以简化计算。有一个正方形发出恒定的 radiance(辐射率)，计算点 p 处的 irradiance(辐照度)。在各个方向上积分来求该值是不太直接的，因为给定一个方向需要判断该方向上正方形是否可见。在四边形的区域上积分来计算辐照度则更容易。
微分区域和微分立体角有如下关系，其中 θ 为正方形微分表面 dA 的法线和指向 p 点的向量的夹角，r为 p 和正方形微分表面 dA 之间的距离。

BRDF(bidirectional reflectance distribution function) 形式化描述了表面的反射。其忽略了光在次表面传播的效果。
BTDF(bidirectional transmittance distribution function) 描述了穿过表面的光的分布。BTDF 不满足互换性(reciprocity).
BSDF(bidirectional scattering distribution function) 描述了表面反射和透射，其同时考虑了 BRDF 和 BTDF.

下图为 BRDF BTDF BSDF 的定义:

• Solid Angle
  
• brdf 能量守恒条件证明
  

Tips: brdf 的范围为 0 到无穷大。例如：完美镜面反射的 brdf 就是一个 δ 函数，其在完美反射方向的值为无穷大。

• https://www.gamedev.net/forums/topic/257668-brdf-range----greater-than-01-/
  
• http://resources.mpi-inf.mpg.de/departments/d4/teaching/ws200708/cg/slides/CG07-Brdf+Texture.pdf
  

BSSRDF(bidirectional scattering surface reflectance distribution function) 描述了材质的散射中展现了明显的次表面光传播。

下图为 BSSRDF 的定义()：

// 计算透视投影矩阵
Transform Perspective(Float fov, Float n, Float f)
{
    // Perform projective divide for perspective projection
    Matrix4x4 persp(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, f / (f - n), -f * n / (f - n),
                    0, 0, 1, 0);

    // Scale canonical perspective view to specified field of view
    Float invTanAng = 1 / std::tan(Radians(fov) / 2);
    return Scale(invTanAng, invTanAng, 1) * Transform(persp);
}

// 透视摄像机的 CameraToScreen变换矩阵为 Perspective(fov, 1e-2f, 1000.f)

上面透视投影变换的求解分为两步：

1. 将点 p 投影到近平面得到点 p'。如下图所示：
   
   
2. 利用用户指定的 FOV，将投影后的视景体内的坐标转化到[-1,1]范围内的坐标。
   对于正方形图片来说，x和 y 都在[-1,1]范围内。
   对于矩形图片来说，比较窄的方向会被映射到[-1, 1]范围，比较宽的方向会同比例映射到更大的范围。
   

真实的相机有镜头系统，其将通过有限尺寸光圈的光聚焦到胶片平面上。相机设计者面临一个取舍：光圈越大，就会有越多的光进入胶片，需要曝光的时间也就越短。然而，镜头只会聚焦在单个平面上（焦点平面），离焦点平面越远的物体会越模糊。光圈越大，产生的模糊越明显。

下面是成像系统相关的几个概念：
焦点（像方交点）：在物理学上指平行光线经透镜折射或曲面镜反射后的会聚点。
物方交点：
焦距：平行光从透镜的光心到光聚集之焦点的距离，即透镜中心到焦点的距离。
焦平面：与成像系统的光轴垂直、且包含成像系统焦点的平面。
镜屏距：透镜中心到屏幕（底片）的距离。
物距：物体到透镜中心的距离。
像距：给定物距和焦距情况下，像到透镜中心的距离。
合焦距离：屏幕上的像最清晰时，屏幕据透镜中心的距离。物距给定，合焦距离即是像距。

• https://www.zhihu.com/question/20086562/answer/15877139
  
• https://zhuanlan.zhihu.com/p/21672067
  
• https://baike.baidu.com/item/%E5%87%B8%E9%80%8F%E9%95%9C
  
• 工程光学（六）——几何光学（进阶）https://zhuanlan.zhihu.com/p/48980327
  

下图为凸透镜原理：

对于目前我们使用的简单摄像机模型，我们可以应用光学上经典的近似（薄镜头近似），使用计算机图形学投影模型来模拟有限光圈。薄镜头近似模拟拥有单个球面镜头的光学系统，镜头的厚薄和曲面镜头的直径有很小关系。
在薄镜头近似中，平行的入射光线经过镜头后会聚焦在一点，这点被称为焦点。焦点到镜头的距离被称为焦距 f。如果胶片被放置在镜头后焦点的位置上，则无线远的物体会被聚焦到胶片上的一个点。如下图所示：

下面链接中可以进行动态演示：
https://ophysics.com/l12.html
./PhysicallyBasedRendering/2020_06_12_ConcaveAndConvexLenses.ggb

对于场景中距离镜头为 z 的一点 p，镜头的焦距为 f，p 点被镜头聚焦的点为 p'，其到镜头的距离为 z'(像距)，由成像公式可的如下关系：

• 成像公式证明 https://baike.baidu.com/item/%E6%88%90%E5%83%8F%E5%85%AC%E5%BC%8F
  

上图中的 z'就是 z 对应的像距，z和 z'是一一对应的，如果点没有放置在 z 聚焦位置上，点在胶片上的投影为一个圆盘(disk)而不是一个点。这个圆盘的边界被称为 circle of confusion。圆盘的尺寸由光圈的直径，焦距以及物体到镜头的距离决定。
在实践中，物体不需要正好在焦点平面上才投影出清晰的像，只要 circle of confusion 大致小于胶片上每个像素的大小，物像就是清晰的。物体聚焦的距离范围被称为镜头的 depth of field(景深).
成像公式可用于计算 circle of confusion 的直径。如下图所示：

焦距 f 为 50mm 的镜头，25mm 光圈，清晰像的物距 Zf 为 1m，随着物距 z 的变化，circle of confusion 直径的变化如下图：

从上图可以看出，Z小于 Zf 时，Dc 变化比较快，而 Z 大于 Zf 时，Dc 变化比较慢。

在射线追踪中模拟一个薄镜头非常简单直接。设点 O 为胶片上一点，在镜头上选一个点 A，然后计算经过镜头中心点 Ao 的射线和焦距平面的交点 Z，则 AZ 就是胶片上 O 对应的射线。如下图所示：

通常为了模拟平滑的景深需要为每个像素点追踪很多个射线。
下面代码为镜头半径大于 0 时，镜头发出射线的重新计算代码。

*ray = Ray(Point3f(0, 0, 0), Normalize(Vector3f(pCamera)));
if (lensRadius > 0) {
    // Sample point on lens
    Point2f pLens = lensRadius * ConcentricSampleDisk(sample.pLens);

    // Compute point on plane of focus
    // focalDistance 为胶片到聚焦平面的距离
    Float ft = focalDistance / ray->d.z;
    Point3f pFocus = (*ray)(ft);  // (*ray)(ft) = ray.o + ft * ray.d

    // Update ray for effect of lens
    ray->o = Point3f(pLens.x, pLens.y, 0);
    ray->d = Normalize(pFocus - ray->o);
}

傅里叶变换是傅里叶分析的基础之一，傅里叶变换是函数在频率域上的表示(通常函数是在空间域上表示)。下图表示了两个函数在空间域和频率域上的图像。

下图为傅里叶级数和傅里叶变换的公式：

大多数函数可以被分解为多个带相位的正弦函数的权重和。在频率空间表示一个函数，便于深入了解该函数的特性(正弦函数的频率分布对应了原始函数的频率分布)。使用这种形式，可以使用傅里叶分析透彻理解，采样引入的错误，以及重建(reconstruction)的过程，以及如何减少感知到的错误效果。

频率域表示和空间域表示互逆，下图为一些重要的函数，以及其在频率域的表示。

• 正弦型函数 https://baike.baidu.com/item/%E6%AD%A3%E5%BC%A6%E5%9E%8B%E5%87%BD%E6%95%B0
  
• 傅里叶分析之掐死教程 https://zhuanlan.zhihu.com/p/19763358
  
• 如何通俗地理解傅立叶变换？ https://www.matongxue.com/madocs/473.html
  
• 如何理解傅立叶级数公式？ https://www.matongxue.com/madocs/619.html
  
• 从傅立叶级数到傅立叶变换 https://www.matongxue.com/madocs/712.html
  
• 如何理解傅里叶变换公式？https://www.zhihu.com/question/19714540/answer/1119070975
  
• 傅里叶系列（一）傅里叶级数的推导 https://zhuanlan.zhihu.com/p/41455378
  
• 傅里叶系列（二）傅里叶变换的推导 https://zhuanlan.zhihu.com/p/41875010
  
• Explanation Fourier http://www.jezzamon.com/fourier/
  
• Fourier 变换的性质 https://zhuanlan.zhihu.com/p/148484210
  
• 如何理解离散傅里叶变换及 Z 变换 https://zhuanlan.zhihu.com/p/45114376
  

使用频率空间分析，我们可以形式化地研究采样的性质。采样过程需要我们选择一组等间距的采样位置，并且计算在这些位置上函数的值。将采样过程形式化表示就是和 shah(狄拉克梳状函数)相乘，下图为该函数定义以及对应的采样过程图示：

这些采样值可以用于定义一个重建函数~f(x)，只需要选择一个重建过滤函数 r(x)，然后和采样值进行卷积运算就可以了。为了重建，卷积对重建滤波在采样点中心的实例进行了加权求和缩放。
下图展示了重建函数的数学描述，以及三角滤波重建的一个例子：

上面讲述的获得重建函数的过程显得不必要的复杂，直观地讲，可以通过在样本之间以某种方式插值来得到重建函数。上面形式化的表述，可以方便地在重建过程中应用 Fourier analysis。

通过在频域中分析样本函数，我们可以更深入地认识采样过程。特别是，我们可以确定在怎样的条件下，可以使用样本数据精确重建原始函数。对于此处的讨论，我们假定原始函数 f(x)的带宽为有限的，也就是说存在一个频率 w0，f(x)的频率都小于 w0。带宽有限的函数，其频域表示函数为紧支撑的，也就是说，|w| > w0 时，F(w) = 0。

Fourier Transform 的一个重要性质是两个函数乘积的傅里叶变换为，这两个函数各自傅里叶变换后再卷积；而两个函数卷积的傅里叶变换为，这两个函数各自傅里叶变换后再相乘。通过该性质可得出，对原函数的采样，在时域空间是原函数(时域表示)和 Shah 函数的乘积，在频域空间是原函数（频域表示）和另一个 Shah 函数的卷积。周期为 T 的 Shah 函数，其频域表示(傅里叶变换)为周期为 1/T 的 Shah 函数。这说明时域空间上样本间隔越远，频域空间上样本间隔越近。

频域空间的样本表示是通过 F(w)和一个新的 shah 函数的卷积得到的，将一个函数和 delta 函数卷积得到的是该函数的一个复制，所以，将函数和 shah 函数卷积的结果为无数个原函数的复制，这些复制之间的间隔为新的 shah 函数的周期(1/T)。
得到无穷多个原始函数分布的复制后，我们只要保留中心在原点的副本，忽略掉其他副本就可以得到原始函数了。为了达到这个目标，我们需要使用一个 box 函数。下图为该过程(重建过程)的图示：

由此可得，我们只需要对原函数进行一系列固定间隔的采样，就可以精确得到原函数在频域空间的表示，除了需要假定原函数为带宽有限的外，不需要知道其他任何信息。
在时域空间应用相同的过程将会精确重建原函数。因为 Box 函数的傅里叶逆变换为 sinc 函数，所以时域空间中理想的重建函数为：

不幸的是，sinc 函数定义域是[-Infinity, Infinity]，必须使用所有的样本值 f(i)才能计算出时域空间特定的~f(x)值。实践中，通常选择有限边界的滤波函数，尽管这会导致无法精确重建原始函数。
在图形学中，通常使用 Box 函数作为滤波函数，其可以有效地在 x 附近的区域对所有样本值进行平均。从 Box 滤波函数在频率域的行为就可以看出，这是一个非常差的选择，这种技术试图通过乘一个 sinc 函数来隔离函数分布的中心副本，这不仅对于选择中心副本来说效果很差，而且还包含了其他无限多个副本的高频部分。

• 如何通俗易懂地解释卷积？ https://www.zhihu.com/question/22298352
  
• 图示卷积 https://www.desmos.com/calculator/ea96vohtuq
  
• 图示卷积 ./PhysicallyBasedRendering/2020_06_21_convolution.ggb
  
• 图示卷积 https://www.geogebra.org/m/KMrAFBxN
  
• 卷积神经网络的卷积核 和 数学中的卷积 https://www.zhihu.com/question/52237725
  
• 《傅里叶光学（九）》 信号采样与重建 https://zhuanlan.zhihu.com/p/72079283
  
• 麻省理工学院公开课：信号与系统：模拟与数字信号处理 https://open.163.com/newview/movie/courseintro?newurl=%2Fspecial%2Fopencourse%2Fsignals.html
  
• 实变函数中的支撑（support）是什么意思？ https://www.zhihu.com/question/329984717
  

除了 sinc 函数无限边界的问题外，理想采样和重建方式的另一个严重的问题是假定信号的带宽是有限的。对于那些带宽无限的信号，或者相对信号自身的频率来说采样率不足够高的信号，前面描述的重建过程会得到一个和原函数不同的函数。
成功重建的关键是，将采样分布乘一个适当宽度的 Box 函数。信号的副本分布是由空格分割的，所以完美的重建是有可能的。如果对原函数的采样使用比较低的频率。这意味这时域空间中样本之间的间隔增大，频域空间中样本之间的间隔缩小，也就是 F(w)的副本之间靠的更近。如果这些副本靠的太近，他们就会重叠。下图展示了采样率太低的情况：

因为这些副本会叠加在一起，最终的分布将看起来就不像是原始函数的多个副本了。当新的分布被乘一个 box 函数，得到的分布看起来像，但不再等于原始 F(w):原始信号的高频细节漏到了重建信号分布的低频区域了。这些新的低频 artifacts 被称为 aliases（走样）(因为高频伪装成了低频)，并且结果信号被称为 aliased（走样了）。下图展示了这种走样现象：

解决频谱重叠问题的一个可行方案是简单增加采样率，直到频谱的副本分开足够远而没有重叠，这样就消除了走样。采样理论可以精确告诉我们需要多大的采样率。理论表示，当采样率大于两倍的信号最大频率，就可以完美重建原始信号。这个最小的采样率被称为 Nyquist freqqency.

对于带宽无限的信号，无法通过足够高的采样率来重建。没有带宽限制的信号其频谱是无限的，所以无论如何分离这些副本，他们都始终叠在一起。不幸的是，图形学中只有少数几个有趣的函数带宽是有限的。特别是，任何不连续的函数都不是带宽有限的，因此，我们无法完美地采样和重建它们。这是有道理的，因为不连续的函数不连续的地方总是落在两个样本之间，而样本无法提供不连续的信息。因此，需要应用不同的方法来抵消这种错误。

如果不小心地采样和重建，大量的 artifacts 将在最终的图片上出现。有时候区分采样和重建导致的 artifacts 很有用，当我们希望精确区分的时候，我们称采样导致的 artifacts 为 prealiasing，而重建导致的 artifacts 为 postaliasing。任何尝试纠正这些错误的方法被统称为 antialiasing.

运用这些理论于 2D 情况下采样和重建渲染场景的图片是非常直接的。我们可以把图片当作 2D 图片位置到辐射率值 L 的函数。
好消息是，使用我们的射线追踪，我们可以计算该函数在任意一点的值。坏消息是，无法在采样之前预过滤 f 函数来移除高频部分。因此，本章的采样器会使用 Nonuniform Sampling 和 Adaptive Sampling 策略。
将场景函数的定义一般化为更高维度的函数是非常有用的，让其也依赖时间 t 以及其采样的镜头的位置(u,v)。因为所有的从摄像机发射的射线都基于 5 个量，改变其中任意一个都会得到不同的射线，也就可能得到不同的值。对于特定的图片位置，该点的辐射率值通常会随时间而变换(场景中有运动的物体)，也会随镜头位置而变化（如果摄像机包含一个带有限光圈的镜头）。
更一般地，因为 14 到 16 章的很多积分器使用了统计技术来估算沿着给定射线的辐射量，当重复给定同一个射线时，他们可能返回一个不同的辐射率值。如果我们将来扩展了场景辐射率函数来包含积分器使用的样本值(例如，在区域光源上选择点来计算照明)我们将需要一个更高维的函数。
高效采样所有这些维度是生成高质量图片重要的一部分。

几何体是渲染图片中最常见的导致 aliasing 的原因。当投影到图片平面时，对象的边缘引入了一个阶跃函数（图片函数突然从一个值跳跃为另一个值）。阶跃函数不仅有无限频率内容，而且完美重建滤波函数应用于走样的样本时会导致 artifacts：ringing artifacts 在重建函数中出现，这种效果被称为 Gibbs phenomenon（吉布斯现象）。面对 aliasing 时，选择一个有效的重建滤波函数需要结合科学、艺术以及个人品味。下图展示了吉布斯现象：

场景中非常小的物体也会导致几何 aliasing。如果几何体足够小其落在了样本之间，其会出乎意料的消失，在过一个动画过几帧后又会出来。

另一个 aliasing 的来源为物体的贴图和材质。没有被正确过滤的贴图映射或者闪耀表面上小的亮点都可能导致渲染的 aliasing。如果采样率没有足够高，无法足够采样这些特性，结果就会 aliasing。物体投射的尖锐的阴影会引入另一种阶跃函数。尽管可以通过几何边缘在图片平面上确定阶跃函数的位置，从阴影边界来检测阶跃函数是非常难的。

我们无法去掉所有导致 aliasing 的源由，所以我们必须发展技术缓和其在最终图片上的影响。

在本章剩余部分，需要铭记两个关于像素的观念。
第一个：组成图片的像素是图片函数的很多点的样本，这些点样本在图片平面上为离散的点。注意一个像素关联的不是区域，将像素当作小的正方形区域是错误的模型，会导致一系列错误。
第二个：最终图片上的像素被定义在像素网格上的离散整数坐标(x,y)上，但是本章的 Sampler 生成的图片样本在连续的浮点坐标(x,y)上。一种自然的在这两种不同的域上映射的方式是将连续的浮点坐标舍入为最近的离散坐标。这种方式很有吸引力，因为其会将刚好等于离散坐标值的连续的坐标映射为离散坐标。然而，结果是给定一组间隔为[x0,x1]的离散坐标，一组覆盖这个范围的连续的坐标的范围为[x0-1/2,x1+1/2).这样，为一个给定离散像素范围生成的连续样本坐标有 1/2 的偏移。
按照下面方式转化离散坐标和连续坐标，离散的范围[x0,x1]会被转化为连续的范围[x0,x1+1)，这样更自然，并且使得代码更简单：

傅里叶分析给我们提供了评价采样模式质量的方法，但是其只能让我们可以量化在带宽有限的情况下均匀增加更多采样所带来的改进。给定带宽无限(图片中物体边缘的地方)以及大于 2 维采样向量(Monte Carlo 光照传播算法)的情况，只使用傅里叶分析是不够的。

给定一个渲染器和一个候选的放置样本的算法，一种验证该算法有效性的方法是使用其采样模式渲染一张图片，然后比较该图片和另一张使用了更多数量样本渲染得到的图片。我们将使用这种方式来比较采样算法。

傅里叶分析之外，数学家发明了一种被称为 Discrepancy（差异）的概念，其可用于评估一个 n 维采样模式的质量。分布比较好的模式其具备低差异值。

差异的基础理念为，可以通过下面方法来评估一组在 n 维空间[0,1)^n 上的样本点的质量，计算每个区域中样本点的数量，比较每个区域的体积和其中的样本点数量。下图在二维空间展示了这种理念：

Discrepancy 的数学定义如下。关于 Discrepancy，Star Discrepancy，以及均匀分布的一维样本模式的 Star Discrepancy 请参考 知乎-PBRT-E7.2-采样接口 https://zhuanlan.zhihu.com/p/73943687 文章的描述。

只使用差异来衡量采样模式，并不是一个很好的标准。一些低差异的采样模式展现出样本的聚集。7.7 节中的 Sobol 采样就有该问题。直观地看，靠的太近的样本就是没有充分利用样本资源：一个样本越靠近另一个，那么该样本就越可能没有提供新的信息。因此计算任意两个样本的最小距离也是一个有用的衡量标准，距离越大采样模式的质量越好。

有很多算法可以生成 Poisson Disk 采样模式，按照样本点距离的标准，这种模式评分很高。人眼的视杆细胞和视锥细胞就是按这种模式分布的，这进一步说明这种方法对于成像来说是一个好的模式。实践中发现 Poisson disk 模式对于 2D 图片来说非常好，但是对于高维采样来说，low discrepancy 模式要更有效。

• 低差异序列（一）- 常见序列的定义及性质 https://zhuanlan.zhihu.com/p/20197323
  
• Star Discrepancy https://mathworld.wolfram.com/StarDiscrepancy.html
  
• PBRT-E7.2-采样接口 https://zhuanlan.zhihu.com/p/73943687
  

Sampler 的实现存储了一系列当前样本的状态，例如：哪个像素被采样，样本向量的维度，等等。因此单个 Sampler 对象被多个线程并行使用时，自然会不安全。Sampler 类实现了 Clone 函数，每个线程都会克隆初始化的 Sampler 对象，每个 Sampler 对象都携带一个随机种子，这样不同的线程就会得到不同的随机数序列。在多个 image tiles 之间重用相同的随机数序列会导致 image artifacts。

Sampler 实现同时处理了获取样本向量各个分量的请求，其为这些值的存储申请空间。
std::vector<int> samples1DArraySizes; 存储请求的样本数组的大小。
std::vector<int> samples2DArraySizes;
std::vector<vector<Float>> sampleArray1D; 存储整个像素的样本数组。
std::vector<vector<Point2f>> sampleArray2D;
size_t array1DOffset; 记录下一个样本对应的样本向量数组的偏移。
size_t array2DOffset;

Request1DArray(int n); n 为请求访问的样本数量
Request2DArray(int n);

Get1DArray(int n); 获取当前样本对应的各个维度数据(样本向量)。
Get2DArray(int n);

Halton 和 Hammersley 序列是两种紧密相关的低差异点集。它们都基于一种被称为 radical inverse 的构造方法。这种方法的基础是任何一个正整数都可以被表示为 b 进制的数字序列。
Halton 和 Hammersley 的定义可以参考原书。

inline uint32_t ReverseBits32(uint32_t n)
{
    n = (n << 16) | (n >> 16);
    n = ((n & 0x00ff00ff) << 8) | ((n & 0xff00ff00) >> 8);
    n = ((n & 0x0f0f0f0f) << 4) | ((n & 0xf0f0f0f0) >> 4);
    n = ((n & 0x33333333) << 2) | ((n & 0xcccccccc) >> 2);
    n = ((n & 0x55555555) << 1) | ((n & 0xaaaaaaaa) >> 1);
    return n;
}

inline uint64_t ReverseBits64(uint64_t n)
{
    uint64_t n0 = ReverseBits32((uint32_t)n);
    uint64_t n1 = ReverseBits32((uint32_t)(n >> 32));
    return (n0 << 32) | n1;
}

Hammersley 和 Halton 序列有一个缺点，就是随着基数 b 的增加，样本值会展现出极其规律的模式。这个问题可以通过扰乱（scrambled）Halton 和 Hammersley 序列来解决，在计算基数取反时对数字进行一次排列。

• 低差异序列（一）- 常见序列的定义及性质 https://zhuanlan.zhihu.com/p/20197323
  

crop window 指定了实际被存储和写出的像素的边界。Crop window 有利于调试，并且将大图片划分为多个小图使其可以在多个电脑上被渲染。Crop window 为 NDC 空间的坐标值。

  struct Pixel
  {
    Pixel() { xyz[0] = xyz[1] = xyz[2] = filterWeightSum = 0; }
    // 运行时对像素有贡献的光谱权重和被表示为XYZ颜色并且被存储在Pixel对象的xyz成员变量中
    Float xyz[3];
    // filterWeightSum存储了样本对像素贡献的过滤器权重和
    Float filterWeightSum;
    // splatXYZ存储了无权重的样本和
    AtomicFloat splatXYZ[3];
    // 该成员没有实际用途，只是确保Pixel对象为32字节而不是28字节。
    // 这可以让Pixel对象在cache中对齐，避免缓存失败
    Float pad;
  };

  // 存储Film对应的像素数据
  std::unique_ptr<Pixel[]> Film::pixels;
  // filter 预计算值表
  static int filterTableWidth = 16;
  Float filterTable[filterTableWidth*filterTableWidth];

选择 XYZ 颜色而不是 RGB 是为了强调 XYZ 是独立于显示器的颜色表示方式，而 RGB 需要假定一组特定的显示器响应曲线。

Film 预计算了一个过滤器值表，从而可以避免耗时的虚函数 Filter::Evaluate()调用。没有用每个样本的精确位置来计算过滤器值会引入一些错误，但是实践上这种错误并不明显。
Film 的责任是确定整数像素值的范围，Sampler 的责任是为这些像素值生成样本。Film::GetSampleBounds 返回被采样的区域，因为像素重建过滤器通常跨越多个像素，因此 Sampler 生成的图片样本必须稍为超出被输出的像素的范围。这样即使是边界处的像素也会有相同的样本密度。这对于通过 crop windows 一片一片渲染图片来说也很重要，其避免了子图片边界处的 artifacts。

Bounds2i Film::GetSampleBounds() const
{
    Bounds2f floatBounds(Floor(Point2f(croppedPixelBounds.pMin) +
                               Vector2f(0.5f, 0.5f) - filter->radius),
                         Ceil(Point2f(croppedPixelBounds.pMax) -
                              Vector2f(0.5f, 0.5f) + filter->radius));
    return (Bounds2i)floatBounds;
}

有三种方式将样本的贡献提供给胶片。第一种方式是由 Sampler 在图片的 tiles 上生成的样本驱动的。尽管这种最简单直接的方法允许渲染器提供一个胶片像素位置以及一个直接指向胶片像素的射线携带的光谱，但是在多线程下无法为这种方式提供一个高效的实现，因为这种情况下多个线程可能会同时更新图片的同一部分。

Film 定义了一个接口，多个线程可以指定其生成的样本在某个像素范围。给定样本范围，GetFileTile()返回指向 FileTile 对象的指针，FileTile 中存储了在对应区域中，样本对该区域像素的贡献。FileTile 和其存储的数据对于调用者来说是独自占有的，因此线程可以提供样本值给 FilmTile 而不需要担心其他线程争夺。当线程在 tile 上的工作完成后，其会将完成的 tile 传回给 Film，Film 会将 tile 合并到最终的图片中。

  // pFilm为胶片上样本的位置
  // L为该样本对应的辐射率
  // sampleWeight 为该样本的权重
  void FilmTile::AddSample(const Point2f &pFilm, const Spectrum &L, Float sampleWeight = 1.)
  {
      // 计算当前样本影响的像素区域
      // Compute sample's raster bounds
      Point2f pFilmDiscrete = pFilm - Vector2f(0.5f, 0.5f);
      Point2i p0 = (Point2i)Ceil(pFilmDiscrete - filterRadius);
      Point2i p1 = (Point2i)Floor(pFilmDiscrete + filterRadius) + Point2i(1, 1);
      p0 = Max(p0, pixelBounds.pMin);
      p1 = Min(p1, pixelBounds.pMax);

      // Loop over filter support and add sample to pixel arrays

      // 预计算对于每个像素该样本对应的过滤器权重
      // Precompute $x$ and $y$ filter table offsets
      int *ifx = ALLOCA(int, p1.x - p0.x);
      for (int x = p0.x; x < p1.x; ++x)
      {
          // 查找表索引计算方法 distance / Radius * filterTableWidth
          // filterTableWidth = filterTableHeight = filterTableSize
          Float fx = std::abs((x - pFilmDiscrete.x) * invFilterRadius.x * filterTableSize);
          ifx[x - p0.x] = std::min((int)std::floor(fx), filterTableSize - 1);
      }
      int *ify = ALLOCA(int, p1.y - p0.y);
      for (int y = p0.y; y < p1.y; ++y)
      {
          Float fy = std::abs((y - pFilmDiscrete.y) * invFilterRadius.y * filterTableSize);
          ify[y - p0.y] = std::min((int)std::floor(fy), filterTableSize - 1);
      }

      for (int y = p0.y; y < p1.y; ++y)
      {
          for (int x = p0.x; x < p1.x; ++x)
          {
              // Evaluate filter value at $(x,y)$ pixel
              int offset = ify[y - p0.y] * filterTableSize + ifx[x - p0.x];
              Float filterWeight = filterTable[offset];

              // 计算过滤后的样本对于像素的贡献
              // Update pixel values with filtered sample contribution
              FilmTilePixel &pixel = GetPixel(Point2i(x, y));
              pixel.contribSum += L * sampleWeight * filterWeight;
              pixel.filterWeightSum += filterWeight;
          }
      }
  }

给定显示设备的响应特性，可以将像素值从设备独立的 XYZ 值转化为依赖于设备的 RGB 值。这种转化对应于另一种光谱基的改变，新的基是由显示设备的光谱响应曲线决定的。此处，从 XYZ 到设备 RGB 的转化使用的权重是基于 sRGB 的。sRGB 是标准的颜色空间，其被所有显示设备和打印设备所支持。

随着输出的 RGB 值被初始化，最终的值由，计算像素过滤方程得到即 contribSum/filterWeightSum。这个转变可能会导致 RGB 值的某些分量为负值；这些是色域外的颜色，其无法被选择的显示基色表示。有多种方案可以解决该问题，例如，将值截断到 0，或者偏移所有分量使其落在色域内，或者基于所有的像素执行一个全局的优化。重建的像素也可能会为负值，这是由于重建过滤函数中存在负波瓣。此处，将颜色值截断为 0 来解决同时这两种情况导致的负值。

考虑 4 维 BRDF 或 BTDF 聚合行为是很有用的。定义一个在多对方向上的函数，然后将其降低为只在一个方向上的 2D 函数来描述整体的散射行为。

Hemispherical-directional reflectance (半球-方向 反射)为一个 2D 函数，对于半球空间上为常量照明的情况下，其可以给出在给定方向上的总的反射。等价地，其可以给出由于从给定方向上的入射光导致的在半球空间上的整体反射。其定义如下：

推导：

// rho 就是希腊字母的 \rho 读“肉” 表示 reflectance ^_^
BxDF::rho(const Vector3f &wo, int nSamples, const Point2f *samples)

Hemispherical-hemispherical reflectance（半球-半球 反射）当入射光在各个方向上都相同时，其可以给出被表面反射的入射光部分。其定义如下：

推导：

BxDF::rho(int nSamples, const Point2f *samples1, const Point2f *samples2)

Tips: 此处的 \(\rho\) 描述的是某种光照情况下整体的散射行为。

• \(\rho_{hd}\) 表示入射光为特定方向，反射光在半球空间均匀分布的情况下，此时的反射率。或者入射光在半球空间均匀分布，出射光为特定方向
  
• Reflectance https://en.wikipedia.org/wiki/Reflectance
  
• pbr 中什么时候调用 rho？
  pbr 中由 BSDF::rho 依次调用其包含的 BxDF::rho，但是源代码中并没有调用 BSDF::rho 的地方。
  

除了反射和透射方向，计算入射光的反射和透射量也很重要。菲涅尔方程描述了表面反射的光的量。菲涅尔方程为 Maxwell 方程在光滑表面下求解出来的。

不同材质分类：

1. dielectrics 电介质。这种材质不会导电。
   他们的折射率为实数，通常其折射率范围为 1 到 3。
   这种材质会透射一部分照明。典型的电介质有玻璃，石油，水以及空气。
   
2. 导体。这种材质中，价电子可以在其原子晶格内自由移动，从而使电流从一个地方流到另一个地方。当导体受到电磁辐射(如可见光)时，其行为会非常不同于电介质：
   材质会表现为不透明。会反射回明显的一部分光。另一部分光会进入导体内部，并被很快吸收（通常在材质表面 0.1 微米内就被完全吸收了，因此只有从非常薄的金属透过的光才能被可感知到）。因此在 pbrt 中忽略导体的透射。
   和电介质不同，导体的折射率为复数。n' = n + ik
   
3. 半导体。硅、锗等属于半导体。pbrt 中不考虑这种材质。
   

导体和电介质都遵循相同的菲涅尔方程。尽管如此，我们更倾向于为电介质创建一个特殊的方程，从而从这种特定的简单形式中获得便利。
下图为电介质对应的菲涅尔方程：

下图为部分电介质的折射率：

当光从一种介质 A 到达另一种折射率比较低的介质 B，在入射角接近 90 度时，将没有光进入介质 B。发生这种全反射的最小入射角被称为临界角。当入射角大于临界角时，发生全反射，此时不需要计算菲涅尔方程。

电介质和导体边界处的菲涅尔方程为：

FrDielectric(Float cosThetaI, Float etaI, float etaT) 函数用来为电介质材质计算 Fresnel 反射公式。
FrConductor(Float cosThetaI, cost Spectrum &etai, const Spectrum &etat, const Spectrum &k) 函数用来为导体材质计算 Fresnel 反射公式。
方便起见，实现了 Fresnel FresnelConductors FresnelDielectrics 三个类来用于计算菲涅尔反射的系数：

  class Fresnel {
  public:
    // Fresnel Interface
    virtual ~Fresnel();
    virtual Spectrum Evaluate(Float cosI) const = 0;
    virtual std::string ToString() const = 0;
  };

  class FresnelConductor : public Fresnel {
  public:
    Spectrum Evaluate(Float cosThetaI) const
    {
      return FrConductor(std::abs(cosThetaI), etaI, etaT, k);
    }
    FresnelConductor(const Spectrum &etaI, const Spectrum &etaT,
                     const Spectrum &k)
      : etaI(etaI), etaT(etaT), k(k) {}
    std::string ToString() const;
  private:
    Spectrum etaI, etaT, k;
  };

  class FresnelDielectric : public Fresnel {
  public:
    Spectrum Evaluate(Float cosThetaI) const
    {
      return FrDielectric(cosThetaI, etaI, etaT);
    }
    FresnelDielectric(Float etaI, Float etaT) : etaI(etaI), etaT(etaT) {}
    std::string ToString() const;

  private:
    Float etaI, etaT;
  };

下图展示了 Gold 的折射率实数部分和虚数部分，这两个值都和波长相关：

• 菲涅耳方程 https://zhuanlan.zhihu.com/p/31534769
  
• 什么是偏振光？http://www.coozhi.com/shenghuojiaju/shenghuochangshi/110199.html
  

可以通过 Schlick's approximation 来近似计算菲涅尔方程：
\(F = F0 + (1-F0)(1-cos(I))^5\)
\(F0 = ((n1-n2)/(n1+n2))^2\)

I 为入射角
F0 为入射角为 0 时，菲涅尔方程的值
入射角很大时，Fresnel 现象会很明显，此时 Fresnel 值几乎为 1。上面的 Schlick 近似可以看作是对 I=0 （F=F0）和 I=90 (F=1)时的插值，只不过插值方式不是线性的。

入射光垂直于表面的情况下（入射角度为 0 时），金属对不同波长的光的反射量不同，因此入射光为白色时金属的反射光会有颜色；电介质对不同波长的光的反射量相同，因此反射光只会发生亮度变化。

入射角度为 0 时，电介质的反射量很小，金属的反射量比较大。随着入射角增大，在靠近掠射角的地方，电介质的反射量迅速增大，入射角为 90 度时，全部反射。

• https://en.wikipedia.org/wiki/Schlick%27s_approximation
  

Tips:
完美镜面反射下，fr(双向反射分布函数)和 Fr(菲涅尔方程)并不相等！双向反射分布函数是 dLo 和 dEi 的比值，其需要考虑所有入射光方向综合的情况。菲涅尔方程是 dLo 和 dLi 的比值，其不需要考虑各个方向综合的情况。

// SpecularReflection BxDF 为 反射+完美镜面
SpecularReflection(const Spectrum &R, Fresnel *fresnel)
: BxDF(BxDFType(BSDF_REFLECTION | BSDF_SPECULAR)),
  R(R),
  fresnel(fresnel) {}

// 对于任意给定的入射方向和出射方向，没有反射
Spectrum SpecularReflection::f(const Vector3f &wo, const Vector3f &wi) const
{
    return Spectrum(0.f);
}

// 对于给定出射方向，返回期望的入射光方向，以及反射的光谱
Spectrum SpecularReflection::Sample_f(const Vector3f &wo, Vector3f *wi,
                                      const Point2f &sample, Float *pdf,
                                      BxDFType *sampledType) const {
    // Compute perfect specular reflection direction
    *wi = Vector3f(-wo.x, -wo.y, wo.z);
    *pdf = 1;
    return fresnel->Evaluate(CosTheta(*wi)) * R / AbsCosTheta(*wi);
}

镜面透射向量通过折射定律来获得：

• 折射定律 https://baike.baidu.com/item/%E6%8A%98%E5%B0%84%E5%AE%9A%E5%BE%8B
  

// SpecularTransmission 的BxDF类型为 透射+完美镜面
SpecularTransmission(const Spectrum &T, Float etaA, Float etaB, TransportMode mode)
: BxDF(BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR)),
  T(T),
  etaA(etaA),
  etaB(etaB),
  fresnel(etaA, etaB),
  mode(mode) {}

  // FresnelSpecular 的BxDF类型为 反射+透射+完美镜面
  FresnelSpecular(const Spectrum &R, const Spectrum &T, Float etaA,
                      Float etaB, TransportMode mode)
  : BxDF(BxDFType(BSDF_REFLECTION | BSDF_TRANSMISSION | BSDF_SPECULAR)),
    R(R),
    T(T),
    etaA(etaA),
    etaB(etaB),
    fresnel(etaA, etaB),
    mode(mode) {}

  // FresnelSpecular f 返回0，因为是完美镜面，所以指定入射方向和出射方向对应的辐射率为0
  Spectrum FresnelSpecular::f(const Vector3f &wo, const Vector3f &wi) const {
      return Spectrum(0.f);
  }

FresnelSpecular 详细实现请参考下面部分：
FresnelSpecular Sample_f 以及 Pdf

下图展示了 Oren-Nayar 模型

微表面分布函数定义在和 BSDF 一样的坐标系下。

    // 微表面分布
    class MicrofacetDistribution
    {
    public:
        // MicrofacetDistribution Public Methods
        virtual ~MicrofacetDistribution();
        // 法线分布
        virtual Float D(const Vector3f &wh) const = 0;
        virtual Float Lambda(const Vector3f &w) const = 0;
        Float G1(const Vector3f &w) const {
          return 1 / (1 + Lambda(w));
        }
        // 几何分布
        Float G(const Vector3f &wo, const Vector3f &wi) const {
          return 1 / (1 + Lambda(wo) + Lambda(wi));
        }
        virtual Vector3f Sample_wh(const Vector3f &wo, const Point2f &u) const = 0;
        // wo出射方向，wh法线方向对应的概率密度
        Float Pdf(const Vector3f &wo, const Vector3f &wh) const;
        virtual std::string ToString() const = 0;

    protected:
        // MicrofacetDistribution Protected Methods
        MicrofacetDistribution(bool sampleVisibleArea)
          : sampleVisibleArea(sampleVisibleArea) {}

        // MicrofacetDistribution Protected Data
        const bool sampleVisibleArea;
    };

  Float MicrofacetDistribution::Pdf(const Vector3f &wo, const Vector3f &wh) const
  {
      if (sampleVisibleArea)
          return D(wh) * G1(wo) * AbsDot(wo, wh) / AbsCosTheta(wo);
      else
          return D(wh) * AbsCosTheta(wh);
  }

为了保证物理上的合理性，微表面函数必须要归一化的（即遵守一定的约束条件）。直观地讲，如果我们考虑在法线方向上，微表面上的入射射线，每个射线和微表面只相交一次。形式化的表示就是，给定一个微表面微分区域 dA,微表面在该区域的投影必须等于 dA.如下图所示：

下面链接中详细描述了法线分布需要遵守的约束条件：

• https://www.cs.cornell.edu/~srm/publications/EGSR07-btdf.pdf
  
• Microfacet Distribution Function 微表面分布函数 D
  

下图展示了两种法线分布函数： 入射角比较小时(入射角和微表面法线夹角比较小时)，GGX 比 Beckmann 衰减要快，表现在画面上，GGX 的高光亮点比 Beckmann 要聚集，入射角达到一定数值后，GGX 的衰减比 Beckmann 要慢，表现在画面上，GGX 的高光范围要比 Beckmann 要大。

下面文件为 Beckmann 和 GGX 法线分布函数的可视化：
./PhysicallyBasedRendering/00_03_normal-distribution.ggb

只用微表面的法线分布来描述微表面的性质是不够的。从给定的方向观察或照射时，因为一些微表面背对着指定方向，所以这些微表面是不可见的，另一些微表面由于其他微表面的遮挡也可能是不可见的，这种效果可以通过 Smith 的 masking-shadowing 函数 G1(w, wh)来描述，其给出对于指定的法线 wh，从 w 方向观察或照射，可见的微表面。(0<=G1(w, wh)<=1)。通常来说，微表面可见的概率独立于他们的朝向 wh。

下图展示了 G1 的归一化约束：

注意：正对 w 方向区域并不都是可见的。

给定 D(wh),假设微表面上相邻点的高度无关联，可得如下公式：

G(ωo, ωi) 给出在ωo 和ωi 方向上都可见的微表面。定义 G 需要额外的假设。我们假设两个方向上可见的可能性互相独立，则得到 G 的方程：G(ωo, ωi) = G1(ωo) G1(ωi)
实践中，两个方向上微表面可见性并不是独立的。考虑特殊情况 ωo = ωi ;此时 G(ωo, ωi) = G1(ωo) = G1(ωi); ωo 和ωi 夹角越小相关性越大。假定给定的点越高微表面可见性越大，在此假设下可推出如下 G 的方程：

• Shadowing-Masking Function 几何项 G
  
• 法线分布函数相关总结 https://zhuanlan.zhihu.com/p/69380665
  

Torrance-Sparrow Model 描述的是粗糙表面的反射。
Torrance-Sparrow Model 推导参考 No description for this link

\(dw_h dw_o\) 关系的第二种推导方式：

\(dw_h dw_o\) 关系的第三种推导方式：

• Importance Sampling PDFs (VNDF, Spherical Caps) https://zhuanlan.zhihu.com/p/682281086
  

粗糙表面的透射可以参考下面论文，其中也包含了 brdf 的推导。

• Microfacet Models for Refraction through Rough Surfaces http://www.cs.cornell.edu/~srm/publications/EGSR07-btdf.pdf
  
• Microfacet models for reflection and refraction https://www.cs.cornell.edu/courses/cs6630/2012sp/slides/05ufacet.pdf
  

  MicrofacetReflection(const Spectrum &R, MicrofacetDistribution *distribution, Fresnel *fresnel)
          : BxDF(BxDFType(BSDF_REFLECTION | BSDF_GLOSSY)),
            R(R),
            distribution(distribution),
            fresnel(fresnel) {}
  Spectrum MicrofacetReflection::f(const Vector3f &wo, const Vector3f &wi) const
  {
      Float cosThetaO = AbsCosTheta(wo), cosThetaI = AbsCosTheta(wi);
      Vector3f wh = wi + wo;
      // Handle degenerate cases for microfacet reflection
      if (cosThetaI == 0 || cosThetaO == 0) return Spectrum(0.);
      if (wh.x == 0 && wh.y == 0 && wh.z == 0) return Spectrum(0.);
      wh = Normalize(wh);
      Spectrum F = fresnel->Evaluate(Dot(wi, wh));
      return R * distribution->D(wh) * distribution->G(wo, wi) * F / (4 * cosThetaI * cosThetaO);
  }

MicrofacetTransmission(const Spectrum &T, MicrofacetDistribution *distribution, Float etaA, Float etaB, TransportMode mode)
          : BxDF(BxDFType(BSDF_TRANSMISSION | BSDF_GLOSSY)),
            T(T),
            distribution(distribution),
            etaA(etaA),
            etaB(etaB),
            fresnel(etaA, etaB),
            mode(mode) {}

Spectrum MicrofacetTransmission::f(const Vector3f &wo,
                                   const Vector3f &wi) const {
    if (SameHemisphere(wo, wi)) return 0;  // transmission only

    Float cosThetaO = CosTheta(wo);
    Float cosThetaI = CosTheta(wi);
    if (cosThetaI == 0 || cosThetaO == 0) return Spectrum(0);

    // Compute $\wh$ from $\wo$ and $\wi$ for microfacet transmission
    Float eta = CosTheta(wo) > 0 ? (etaB / etaA) : (etaA / etaB);
    Vector3f wh = Normalize(wo + wi * eta);
    if (wh.z < 0) wh = -wh;

    Spectrum F = fresnel.Evaluate(Dot(wo, wh));

    Float sqrtDenom = Dot(wo, wh) + eta * Dot(wi, wh);
    Float factor = (mode == TransportMode::Radiance) ? (1 / eta) : 1;

    return (Spectrum(1.f) - F) * T * std::abs(distribution->D(wh) * distribution->G(wo, wi) * eta * eta * AbsDot(wi, wh) * AbsDot(wo, wh) * factor * factor / (cosThetaI * cosThetaO * sqrtDenom * sqrtDenom));
}

已知函数 f 在 x0,x1,…,xk 处的一组函数值，和导函数值。对于每一段[x(i),x(i+1)]，我们可以使用一个三次曲线进行近似 pi(x) = 3x^3 + bx^2 + cx + d;
选择 pi(x)函数时，其需要满足以下条件：
pi(xi) = f(xi)
pi(xi+1)=f(xi+1)
pi'(xi) = f'(xi)
pi'(xi+1) = f'(xi+1)
为了简化讨论，我们只考虑[x0,x1]，并且假设[x0,x1]=[0,1]
求解 a,b,c,d 系数可得：
a = f'(x0) + f'(x1) + 2f(x0) - 2f(x1)
b = 3f(x1) - 3f(x0) - 2f'(x0) - f'(x1)
c = f'(x0)
d = f(x0)
带入后可得:
p(x) = (2x^3 - 3x^2 + 1)f(x0) + (-2x^3+3x^2)f(x1) + (x^3-2x^2+x)f'(x0) + (x^3-x^2)f'(x1)
这种插值方式很方便，但是限制依然太强，因为通常我们无法知道导函数的信息。因此我们通过相邻的函数值 f(xi-1)和 f(xi+1)来计算近似的 f'(xi),这样可以计算处最终的 p(x)如下：

• catmull-rom 样条(spline)公式推导 https://zhuanlan.zhihu.com/p/111708587
  

对于每个和几何体相交的射线，Integrator 在计算射线携带的辐射率时都会创建一个或多个 BSDF 对象，每个 BSDF 对象中又包含多个由 Materials 在交点处创建的 BxDFs。如果使用 new delete 为这些 BSDF 和 BxDFs 对象分配内存，效率会很差。使用 MemoryArena 来解决该问题。

BSDF *b = new BSDF;
BxDF *lam = new LambertianReflection(Spectrum(0.5f));

BSDF *b = ARENA_ALLOC(arena, BSDF);
BxDF *lam = ARENA_ALLOC(arena, LambertianReflection)(Spectrum(0.5f));

// 使用arena申请的内存，使用new执行类对象的初始化
#define ARENA_ALLOC(arena, Type) new (arena.Alloc(sizeof(Type))) Type

Matte Material 用于描述完全的漫反射材质

  // Kd 为漫反射光谱值
  // sigma 为粗糙度值(roughness)
  // bumpMap 为法线贴图
  MatteMaterial(const std::shared_ptr<Texture<Spectrum>> &Kd,
                const std::shared_ptr<Texture<Float>> &sigma,
                const std::shared_ptr<Texture<Float>> &bumpMap)
    : Kd(Kd), sigma(sigma), bumpMap(bumpMap) {}

  void MatteMaterial::ComputeScatteringFunctions(SurfaceInteraction *si,
                                                 MemoryArena &arena,
                                                 TransportMode mode,
                                                 bool allowMultipleLobes) const
  {
    // Perform bump mapping with _bumpMap_, if present
    if (bumpMap) Bump(bumpMap, si);

    // Evaluate textures for _MatteMaterial_ material and allocate BRDF
    si->bsdf = ARENA_ALLOC(arena, BSDF)(*si);
    Spectrum r = Kd->Evaluate(*si).Clamp();
    Float sig = Clamp(sigma->Evaluate(*si), 0, 90);
    if (!r.IsBlack())
    {
        // 粗糙度为0 则bsdf为LambertianReflection, 否则为OrenNayar
        if (sig == 0)
          si->bsdf->Add(ARENA_ALLOC(arena, LambertianReflection)(r));
        else
          si->bsdf->Add(ARENA_ALLOC(arena, OrenNayar)(r, sig));
    }
  }

Plastic 可用于模拟反射效果中混合了漫反射和有光泽散射的材质。其实就是典型的电介质材质。

  // PlasticMaterial Method Definitions
  void PlasticMaterial::ComputeScatteringFunctions(SurfaceInteraction *si, MemoryArena &arena, TransportMode mode, bool allowMultipleLobes) const
  {
      // Perform bump mapping with _bumpMap_, if present
      if (bumpMap) Bump(bumpMap, si);
      si->bsdf = ARENA_ALLOC(arena, BSDF)(*si);
      // Initialize diffuse component of plastic material
      Spectrum kd = Kd->Evaluate(*si).Clamp();
      if (!kd.IsBlack())
      {
          // 叠加一个Lambertian反射，kd为 diffuse color
          // LambertianReflection 的Pdf为默认的BxDF::Pdf 此处没有使用kd*(1-ks)来实现粗略能量守恒
          si->bsdf->Add(ARENA_ALLOC(arena, LambertianReflection)(kd));
      }

      // Initialize specular component of plastic material
      Spectrum ks = Ks->Evaluate(*si).Clamp();
      if (!ks.IsBlack()) {
          // 使用FresnelDielectric来计算菲涅尔项
          // 塑料的折射率为 1.5f 空气的折射率为 1.0f (FrDielectric函数内部根据入射角来判断哪个是入射介质，哪个是折射介质)
          Fresnel *fresnel = ARENA_ALLOC(arena, FresnelDielectric)(1.5f, 1.f);
          // Create microfacet distribution _distrib_ for plastic material
          Float rough = roughness->Evaluate(*si);
          if (remapRoughness)
              rough = TrowbridgeReitzDistribution::RoughnessToAlpha(rough);
          MicrofacetDistribution *distrib = ARENA_ALLOC(arena, TrowbridgeReitzDistribution)(rough, rough);
          // 叠加一个微表面镜面反射，ks为反射光光谱分布 specular color
          BxDF *spec = ARENA_ALLOC(arena, MicrofacetReflection)(ks, distrib, fresnel);
          si->bsdf->Add(spec);
      }
  }

  PlasticMaterial *CreatePlasticMaterial(const TextureParams &mp)
  {
      // 默认diffuse color 为0.25,0.25,0.25
      std::shared_ptr<Texture<Spectrum>> Kd = mp.GetSpectrumTexture("Kd", Spectrum(0.25f));
      // 默认specular color 为0.25,0.25,0.25
      std::shared_ptr<Texture<Spectrum>> Ks = mp.GetSpectrumTexture("Ks", Spectrum(0.25f));
      std::shared_ptr<Texture<Float>> roughness = mp.GetFloatTexture("roughness", .1f);
      std::shared_ptr<Texture<Float>> bumpMap = mp.GetFloatTextureOrNull("bumpmap");
      bool remapRoughness = mp.FindBool("remaproughness", true);
      return new PlasticMaterial(Kd, Ks, roughness, bumpMap, remapRoughness);
  }

Mix Material 实现按照权重值将两个 Material 合在一起。

void MixMaterial::ComputeScatteringFunctions(SurfaceInteraction *si,
                                             MemoryArena &arena,
                                             TransportMode mode,
                                             bool allowMultipleLobes) const {
    // scale 存储了混合权重
    // Compute weights and original _BxDF_s for mix material
    Spectrum s1 = scale->Evaluate(*si).Clamp();
    Spectrum s2 = (Spectrum(1.f) - s1).Clamp();
    SurfaceInteraction si2 = *si;
    m1->ComputeScatteringFunctions(si, arena, mode, allowMultipleLobes);   // 计算材质 1
    m2->ComputeScatteringFunctions(&si2, arena, mode, allowMultipleLobes); // 计算材质 2

    // Initialize _si->bsdf_ with weighted mixture of _BxDF_s
    int n1 = si->bsdf->NumComponents(), n2 = si2.bsdf->NumComponents();
    // bxdfs
    for (int i = 0; i < n1; ++i)
        si->bsdf->bxdfs[i] = ARENA_ALLOC(arena, ScaledBxDF)(si->bsdf->bxdfs[i], s1);
    for (int i = 0; i < n2; ++i)
        si->bsdf->Add(ARENA_ALLOC(arena, ScaledBxDF)(si2.bsdf->bxdfs[i], s2));
}

FourierMaterial 支持将测量或人工生成的 BSDF 表格数据转化为基于方向的 BSDF。

考虑任意的贴图函数，其为位置的函数，定义在场景中物体表面上。用 f(x,y)表示渲染图片的点(x,y)到物体表面的点的映射，则贴图函数可以表示为 T(f(x,y))。
考虑贴图函数 T(s,t)被应用于垂直于 z 轴的正方形，该正方形 4 个顶点世界坐标分别为(0,0,0) (1,0,0) (1,1,0) (0,1,0).如果摄像机为正交投影，并放在恰当位置，使得该正方形铺满渲染图片，此时正方形上的点 p 对应的贴图坐标为 s = px t = py;则贴图坐标(s,t)与屏幕像素坐标(x,y)的关系为：s=x/x_r t=y/y_r (渲染图片的分辨率为 x_r,y_r)。给定一个像素的一个样本间隔，贴图参数空间上的样本间隔为(1/x_r,1/y_r),贴图函数必须移除高于当前可表示采样率的细节。

像素坐标和贴图坐标的关系，以及像素采样率和贴图采样率之间的关系是确定贴图允许的最大频率内容的关键。给定一个三角形上的一点，其对应的贴图坐标为(s,t),摄像机为透视投影，可以通过解析的方法得到该点对应的渲染图片上样本点的差值。这是图形处理器进行贴图映射时实现抗锯齿的基础。
对于复杂场景、透视投影的情况，精确地找到渲染图片位置和贴图坐标值的关系是非常困难的。幸运的是，对于贴图抗锯齿，我们不需要对于任意的渲染图片上的点(x,y)求解对应物体表面上的点 f(x,y).我们只需要知道改变像素采样点位置导致的贴图样本位置的变化。像素图片上(x,y)变为(x',y')对应物体上的点从 f(x,y)变为 f(x',y')，通过下面公式计算得到 f(x',y')：

如果这些偏导数对于对应的 x'-x 和 y'-y 变化不大，上面的近似就是合理的。更重要的是，这些偏导数的值给出了，当像素在 x,y 方向上偏移 1 个单位时，贴图样本位置的变化，也就是贴图的采样率。
我们可以使用 offset 射线来估算，
下图展示了计算 dudx dvdx dudy dvdy (这些值就是贴图采样率)的原理：

上面有 3 个方程 2 个未知数(∂u 和∂v)，其中有一个方程是退化的。计算∂p/∂u 和 ∂p/∂v 的叉乘，判断所得向量的最大分量，选另外两个分量(这两个分量上的 uv 变化最大)对应的方程进行求解。

为了计算不带 aliasing 的 T(f(x,y)) ，我们需要先将该函数和 sinc 过滤器卷积，将超过 Nyquist 限制的频率移除掉，这样就得到基带有限的函数。再将基带有限的函数和屏幕上的像素过滤器函数进行卷积，从而得到贴图函数：

在实践中，对于上面的过程可以做很多简化，例如，使用一个 BoxFilter 代替 sinc 用于移除高频率内容，并且将第二步完全忽略，这相当于将 BoxFilter 当作像素过滤器。这样整个 antialiasing 工作就完全可以在贴图空间进行了，极大简化了贴图函数的实现。
理想的 sinc 过滤器可以让低于 Nyquist 限制的频率分量通过，但是移除高于 Nyquist 限制的频率分量。因此，如果我们知道贴图函数的内容频率，我们可以将高频项用其平均值代替，这样就高效地完成了 sinc 预过滤器的工作。
supersampling 技术：贴图函数在主点旁边的多个位置上被求解和过滤，这样就对应增加了贴图空间的采样率。如果 Box 过滤器被应用于这些样本值，这等价于对函数值进行平均。如果贴图函数的求解很复杂，则该方法会很昂贵。但是，其比增加渲染图片采样率要高效很多，因为其不会导致追踪多个穿越场景的射线。

为了计算在表面交点处，反射或透射的微分射线，我们需要对射线 r 进行近似得到 r'，r'和表面有新的交点，则对 r'进行射线追踪。如下图：

r'和表面的交点可以通过对 r 和表面的交点进行偏移得到，偏移量为∂p/∂x 和 ∂p/∂y。
r'的方向可以通过下面公式计算得到：


对应的代码实现在 SamplerIntegrator::SpecularReflect 函数中。

SphereMapping 的原理是将一个球放在物体中心，物体上每个点都沿着该点到球心的方向投影到球表面，使用球上投影点的 uv 作为物体上对应点的 uv。

将物体投影到圆柱体上。

我们可以使用两个不平行的向量 Vs 和 Vt 以及两个偏移量 Ds 和 Dt 定义一个参数化平面。物体上每一点到原点的向量和 Vs、Vt 分别点积再分别进行 Ds、Dt 偏移，从而得到对应的贴图坐标。
假如选择 Vs=(1,0,0) Vt=(0,1,0) Ds=Dt=0，这相当于将 z=0 的平面作为投影平面，此时 s=Px t=Py.

对世界坐标的点进行线性变换(通常取世界空间到物体空间的变换)，直接拿变换后的世界坐标作为 3D 贴图坐标。

ScaleTexture 取两个贴图，返回这两个贴图值的乘积。
ScaleTexture 忽略抗锯齿，将抗锯齿工作留给它的两个子贴图。尽管两个带宽有限的函数的积也是带宽有限的，但是乘积的最高频率可能比其中任何一个子贴图的频率都高。因此，即使两个贴图都进行了抗锯齿操作，它们的乘积依然可能会有锯齿。
幸好该类通常用于拿一个常数去缩放另一张贴图，此时被缩放贴图做抗锯齿就够了。

MixTexture 为更一般化的 ScaleTexture。其取三个贴图作为输入，其中两个可以为任意类型第三个的返回值必须为浮点值。浮点贴图为前两张贴图线性插值的权重。

BilerpTexture 提供在 4 个常数之间进行线性插值。这 4 个常数定义在(s,t)参数空间的(0,0) (1,0) (0,1) (1,1)四个位置处。

在一个场景中，每个图片贴图可能被重复使用多次，因此 pbrt 维护了一张图片表，用于记录已经被加载的图片，这样即使图片被多个 ImageTexture 使用，其也只需要加载一次。ImageTexture 构造函数通过调用 static ImageTexture::GetTexture() 函数得到所需贴图的一个 MipMap 表示。
下图为 sRGB gamma 曲线函数，其为一个分段函数:

将图片数据从文件中读出来后，先转化为 Tmemory 类型的数据，然后再执行缩放和 gamma 矫正，最后将得到的数据存到 mipmap 中。
Tips: scale 参数并不是改变图片分辨率，而是直接对贴图图元值进行缩放。

贴图函数包含的频率过高，贴图采样率不足以表示，会导致 aliasing。对贴图函数进行求值之前，需要通过预过滤操作，移除任何高于 Nyquist 限制的频率。
ImageTexture 的贴图图元是对贴图函数进行固定频率采样得到的一组样本。对贴图进行采样时，需要过滤的区域由中心点的贴图坐标(s, t)以及估算到的邻接图片样本对应的贴图坐标的偏移共同决定。这些偏移值就是估算的贴图采样率，我们必须移除频率高于两倍邻接图片样本的内容，来满足 Nyquist 标准。

贴图采样和重建过程与第七章的图片采样过程有关键的不同。对贴图采样来说，获取样本值是不昂贵的(只需要对数组进行查询就可以了，图片采样则需要追踪额外数量的射线并计算辐射率)。贴图函数完全由一组样本定义，并且贴图函数在样本之间的行为是确定的。这些差异，使得我们可以在采样之前移除贴图的细节，从而消除 aliasing。

MIPMap 实现了两种方法用于过滤宽度可变的高效贴图过滤。trilinear interpolation 和 elliptically weighted averaging(比第一个方法更慢更复杂，但是效果非常好)。
为了限制需要访问的图元数量，这两种过滤方法都使用了 image pyramid 增加原图片的低分辨率预过滤版本来加速运算。原始图片的图元在金字塔底层，一直到最上层，只有一个图元，其为原始图片的所有图元的平均值。这个图片集合需要多于原始图片 1/3 的内存来存储。

pbrt 中实现 image pyramid 时，会判断原始图片的分辨率是否为 2 的幂次，如果不是，会将图片放大到最近的 2 的幂次。放大图片涉及到了第 7 章中讲的采样和重建的理论。我们有一个图片函数，其由按照某个采样率进行采样得到的一组样本表示。我们使用原始样本重建一个连续图片函数，然后对新的图片函数按照一组新的样本位置(新的采样率)进行采样。因为新的采样增加了采样频率，因此不需要担心采样过疏导致高频分量出现。下图展示了，这种放大图片的 1D 图示：

MIPMap 使用了可分离的重建过滤器，可分离的过滤器可以写为两个 1D 过滤器的积：f(x,y) = f(x)f(y). 使用可分离过滤器的一个优点是，当我们对图片进行重新采样使其分辨率从(s,t)变为(s',t')时，可以使用两个 1D 步骤来实现对图片的重新采样，这样我们需要访问的图元数量为过滤器宽度的线性函数，而不是二次函数。
重新构建原始图片函数，然后对其在新的贴图图元位置进行采样。这等价于将重建过滤器核放在新的图元位置，然后加权附近的贴图图元。这样，每个新的图元就是原始图片中一小部分贴图图元的加权平均。

考虑原始分辨率为 3，那么新的分辨率为 4。MIPMap::resampleWeights 生成重新采样权重的原理如下图：

为了计算过滤后的值，每一次贴图查询通常需要使用 8-20 个贴图图元，因此精心选择贴图图元的内存分布来减少访问时的 cache miss 可以极大提高性能。
MIPMap 使用了一个简单的 BoxFilter 对上一层级的 4 个图元进行平均得到当前层级的图元的值。使用 Lanczos Filter 可以得到一个更好的结果。

各项同性的过滤算法的主要缺点是，当以倾斜的角度看贴图时会显得有些模糊。因为倾斜角度看贴图时，沿不同的轴采样率需要不同，而各项同性的过滤算法，各个方向采样率相同。
以很宽的过滤器对很多图元进行过滤会非常低效。因此 MIPMap 会选择某个层级的贴图使得过滤区域只包含 4 个图元。其原理如下图：

双线性过滤的原理：

计算出的 MipMap 层级 level 为 float 值，对 floor(level) 和 floor(level)+1 层贴图执行三角滤波，然后对滤波后的结果进行线性插值就得到最终的贴图值。

if (level < 0)
    return triangle(0, st);
else if (level >= Levels() - 1)
    return Texel(Levels() - 1, 0, 0);
else {
    int iLevel = std::floor(level);
    Float delta = level - iLevel;
    return Lerp(delta, triangle(iLevel, st), triangle(iLevel + 1, st));
}

椭圆加权平均（EWA）算法将椭圆拟合到纹理空间中的两个轴, 这两个轴由纹理坐标差指定，然后使用高斯滤波函数对纹理进行滤波。不像上一节的三角滤波器，该滤波器支持任意朝向的贴图区域，其在不同方向上由不同的过滤边界。这种类型的滤波器被称为各项异性的。

椭圆的离心率表示了椭圆的圆扁程度。离心率越小越接近于圆，越大则椭圆就越扁。EWA 是按照 minor axis 选择 MIPMap 层级的，当离心率大时，意味着需要过滤的图元数量会比较多。将 minor axis 的长度放大来限制椭圆的离心率。
下图为所有在该椭圆内的点的方程：

上面的隐式方程的一个比较好的特性是在特定图元处的值为 从椭圆中心到图元距离 和 从椭圆中心到椭圆边界距离 的平方比。这个值可用于对预计算的高斯过滤器函数值表进行索引。其原理如下图：

关于该节的详细内容可以参考下面 知乎文章 EWA 的讲解。

• 椭圆 https://baike.baidu.com/item/%E6%A4%AD%E5%9C%86
  
• EWA https://zhuanlan.zhihu.com/p/105167411
  

Checkerboard2DTexture 抗锯齿实现：
当整个过滤区域都在同一个 check 中时，不需要执行抗锯齿。通过计算过滤区域的 BoundingBox，如果 BoundingBox 的边界在同一个 check 内，则说明过滤区域在同一个 check 内。此处使用了轴对齐的 BoundingBox，这简化了实现，但增加了过滤值的模糊度。

当整个过滤区域不在同一个 check 中时。首先计算过滤区域覆盖的两种 check 的部分为多少。这等价于计算一个 2D 分段函数在过滤区域上的均值。得到均值后，我们就可以在这两种 check 之间进行混合了。

我们可以基于 3D 贴图坐标定义一个立体的 CheckerBoard。

Perlin noise 在所有的整型位置(x,y,z)处的值都为 0，其变化性来源于每个格子的梯度向量。
梯度向量对所求点的影响通过下图方法来计算：

指定整型格子的梯度向量是通过索引一个预计算的表(NoisePerm)得到的。NoisePerm 表中的值的低 4 位用于确定使用 16 个梯度向量中的哪个。
NoisePerm 的生成方法为，先使用 0 到 NoisePermSize-1 对 NoisePerm 前 NoisePermSize 项进行填充，然后对这些数据进行随机排列，最后用这些数据填充后 NoisePermSize 项。

下面文件为 3 维空间的梯度向量：
./PhysicallyBasedRendering/perlin-3d-grad-vector.ggb

• https://zh.wikipedia.org/wiki/Perlin%E5%99%AA%E5%A3%B0
  
• https://zhuanlan.zhihu.com/p/22337544
  

polka dot texture 将贴图空间分为矩形 cells。每个 cell 有 50%几率包含一个点，这样点就会随机放置在不同的 cell 中。

noise function 是基带有限的函数，这意味着可以通过缩放定义域来改变其频率内容。例如，Noise(2*p)的频率为 Noise(p)频率的二倍。
可以通过下面的公式，使用频谱合成来计算模式从而丰富 noise 的细节：

典型地，对于上面的公式，缩放因子 si、权重因子 wi 选择几何级数(等比数列)，si=2s(i-1)，权重因子 wi = w(i-1)/2。这使得高频率项对整体的影响相对要小。这里附加的项被称为一个八度 octave noise，因为其频率为前一个的两倍。当这个模式应用于 Perlin noise 时，结果被称为分形布朗运动 fractional Brownian motion.

Fractional Brownian motion 对于程序化生成贴图是非常有用的构建要素，其比简单的 noise 有更复杂的变化，同时易于计算并且依然有很好的频率。
fBm 的实现中使用了 clamping to antialias 的技术。其原理为，当我们计算一个基于很多项和的值时，已知每一项的频率，对于频率高于 Nyquist 限制的项，应该将该项的平均值用于求和。Noise()的平均值为 0。
Noise()内容的最大频率粗略为 1。每个子项表示 2 倍频率的内容。假设 noise 空间的采样频率为 s，由下面公式可以计算出最多可以有多少项：
第 1 项的频率为 \(s/2^{n-1} = 1\)
第 2 项的频率为 \(s/2^{n-2}\)
…
第 n-1 项的频率为 \(s/2\)
第 n 项的频率为 \(s\)

和 fBm() 类似的一个函数为 Turbulence()。绝对值会导致一阶不连续性。其图像如下：

fBm 和 Turbulence 函数可用于实现 Bump 贴图和 Wrinkeled 贴图。

fBm 可用于实现波纹。

noise function 的另一种用途为扰乱贴图坐标。

吸收是通过介质的吸收横截面来描述的，其为光在介质中传播单位距离被吸收的概率。尽管吸收截面通常为 p(位置)的函数，但通常会随 p(位置)和 w(光方向)变化。其也会随光的频率而变化。

上图中有一个自然底数 e 的积分，其指数部分的积分表示吸收率的积分。其意义的详细解释可以参考《An Intuitive Guide To Exponential Functions & e》。

• https://baike.baidu.com/item/%E5%90%B8%E6%94%B6%E6%88%AA%E9%9D%A2
  
• https://zh.wikipedia.org/wiki/%E5%90%B8%E6%94%B6%E6%88%AA%E9%9D%A2
  
• An Intuitive Guide To Exponential Functions & e https://betterexplained.com/articles/an-intuitive-guide-to-exponential-functions-e/
  

上面微分方程求解参考如下链接：

• 体渲染数学原理 https://zhuanlan.zhihu.com/p/56710440
  

光线经过介质时，由于化学、热或原子能等过程，将能量转化为可见光，这会增强光线。

光碰到粒子后被散射到不同方向，从而减少该方向上光的辐射量，这种现象被称为 out-scattering。
由于吸收和外散射导致辐射的减少被统称为衰减或消光.
散射系数和消光系数的比值被称为 albedo 反照率，其描述了散射相对于吸收的概率。
1/消光系数被称为平均自由程，其描述了光线在介质中传播时，碰到介质粒子之前所经过的平均距离。

在均匀介质中(homogeneous medium)，衰减系数为常数。

• An Intuitive Guide To Exponential Functions & e https://betterexplained.com/articles/an-intuitive-guide-to-exponential-functions-e/
  
• 用 C 语言画光（七）：比尔-朗伯定律 https://zhuanlan.zhihu.com/p/31901449
  

上面微分方程求解参考如下链接：

• 体渲染数学原理 https://zhuanlan.zhihu.com/p/56710440
  

由于将其他方向的光散射到当前方向上，因此内散射会增加光的辐射量。
假定粒子之间的间隔至少为他们半径的几倍，此时描述特定点的散射时，就可以忽略掉粒子之间的交互。在这种假设下，phase function 可以描述散射辐射在某点的角度分布，其为 BSDF 的体积模拟。

• https://baike.baidu.com/item/%E7%9B%B8%E5%87%BD%E6%95%B0
  

上面微分方程求解参考如下链接：

• 体渲染数学原理 https://zhuanlan.zhihu.com/p/56710440
  

MediumInteraction 表示在散射介质中，某点的相互作用。

  class MediumInteraction : public Interaction
  {
    public:
      // MediumInteraction Public Methods
      MediumInteraction() : phase(nullptr) {}
      MediumInteraction(const Point3f &p, const Vector3f &wo, Float time,
                        const Medium *medium, const PhaseFunction *phase)
          : Interaction(p, wo, time, medium), phase(phase) {}
      bool IsValid() const { return phase != nullptr; }

      // MediumInteraction Public Data
      const PhaseFunction *phase;
  };

齐次介质是最简单的介质。在这种介质区域内，光传播单位距离被吸收的概率σ_a为常量值，光传播单位距离被散射的概率σ_s也为常量值。

// HomogeneousMedium Declarations
class HomogeneousMedium : public Medium {
  public:
    // HomogeneousMedium Public Methods
    HomogeneousMedium(const Spectrum &sigma_a, const Spectrum &sigma_s, Float g)
        : sigma_a(sigma_a), // 吸收概率
          sigma_s(sigma_s), // 散射概率
          sigma_t(sigma_s + sigma_a), // 消光系数
          g(g) {}
    Spectrum Tr(const Ray &ray, Sampler &sampler) const;
    Spectrum Sample(const Ray &ray, Sampler &sampler, MemoryArena &arena,
                    MediumInteraction *mi) const;

  private:
    // HomogeneousMedium Private Data
    const Spectrum sigma_a, sigma_s, sigma_t;
    const Float g;
};

GridDensityMedium 类将介质的密度存储为一个 3D 格子。通过对样本进行插值来得到特定点的介质密度。

对于高反照率介质(high-albedo media),散射辐射分布通常是各项同性的，并且 Fresnel 透射比是定义最终方向分布的最重要的因素。
对于低反照率介质(low-albedo media), 方向变化是有意义的，下面这种近似会少一些精确性。


Tips: Sp 项即 profile 项，其概念和皮肤渲染中的 diffuse profile 应该一样。

解耦 S 的空间和方向参数明显减低了 S 的维度，但是没有从根本上解决难题，对于任意形状的表面 Sp 项很难求解，此处引入近似的方法，用 Sr 来近似 Sp，Sr 的具体计算请参考下一节 Tabulated BSSRDF。

• https://baike.baidu.com/item/%E5%8F%8D%E7%85%A7%E7%8E%87
  
• https://baike.baidu.com/item/%E7%9F%A9
  
• 一阶矩 二阶矩 三阶矩 https://www.cnblogs.com/mld-code-life/p/11172796.html
  

TabulatedBSSRDF 类提供了访问表格表示的 BSSRDF，其可以处理很多散射 profile，包括从真实世界测量得到的 BSSRDF。
需要特别注意的是，当所有的 BSSRDF 的材质属性固定时，radial profile Sr 只是一个 1D 函数。更一般地，其依赖于四个额外的参数：折射率η，散射各向异性 g，反照率ρ，以及消光系数σt。完整的 Sr 函数为 Sr(η, g, ρ, σt , r)，对于离散化来说，这个函数维度过高了。我们需要移除或者固定一些参数。
考虑唯一的有物理单位的参数σt, 该参数用于度量单位距离上散射和吸收交互的概率。其效果只是控制 BSRDF profile 的空间缩放。为了减少维度，我们将σt 固定为 1.
运行时，当查询给定σt 和半径 r 的 Sr 时，我们使用如下方程进行计算：

实践中，我们也会将η, g 固定下来，这意味着这些参数对于物体来说不能贴图化。这些简化，使我们得到一个非常易于管理的 2D 函数，其只需要在ρ(Tips:反照率 ρ=σs/σt), optical radius 两个参数上进行离散化。

注意 TabulatedBSSRDF::rho 成员变量给出了单次散射后能量的减少。这和材质整体的反照率(albedo)不同，材质整体的反照率将所有阶的散射都计算在内。将单次散射反照率记为ρ,将材质整体的反照率记为ρeff。

ρeff 为ρ的非线性单调递增函数。

// 用于计算交点处的散射属性
SubsurfaceMaterial::ComputeScatteringFunctions(SurfaceInteraction *si, MemoryArena &arena, TransportMode mode, bool allowMultipleLobes) const

直接设置吸收和散射系数来达到希望的表现效果是很难的。这些参数都不是线性的，而且也不直观。KdSubsurfaceMaterial 允许用户按照表面漫反射的术语来设置次表面散射属性，以及平均自由程(mean free path)1/σt。需要注意的是，表面属性的变化并不对应于介质中属性的变化。

黑体是一个完美的发射者：其完美地将能量转化为电磁辐射。黑体发射有封闭形式的表达式，其为温度和波长的函数。
黑体吸收所有入射的能量，不会反射出任何能量。

普朗克黑体辐射定律描述，在任意温度 T 下，从一个黑体中发射出的电磁辐射的辐射率与频率彼此之间的关系。

斯特藩-玻尔兹曼定律（Stefan–Boltzmann law）给出了空间中 p 点黑体的出射辐射。
维恩位移定律（Wien’s displacement law）给出了给定温度下黑体辐射的波长。

基尔霍夫定律(Kirchoff ’s law) 描述了非黑体的发射行为，其描述了物体的发射率与吸收比之间的关系。

黑体发射分布提供了有用的度量标准用于描述非黑体发射特征，这种度量标准概念称为色温。如果一个发射体的 SPD 和同温度的黑体的分布相同，我们就称该发射体有对应的色温。

• 黑体 https://baike.baidu.com/item/%E9%BB%91%E4%BD%93
  
• 黑体辐射 https://baike.baidu.com/item/%E9%BB%91%E4%BD%93%E8%BE%90%E5%B0%84
  
• 普朗克黑体辐射定律 https://zh.wikipedia.org/zh-hans/%E6%99%AE%E6%9C%97%E5%85%8B%E9%BB%91%E4%BD%93%E8%BE%90%E5%B0%84%E5%AE%9A%E5%BE%8B
  
• 斯特藩-玻尔兹曼定律 https://baike.baidu.com/item/%E6%96%AF%E7%89%B9%E8%97%A9-%E7%8E%BB%E5%B0%94%E5%85%B9%E6%9B%BC%E5%AE%9A%E5%BE%8B
  
• 维恩位移定律 https://baike.baidu.com/item/%E7%BB%B4%E6%81%A9%E4%BD%8D%E7%A7%BB%E5%AE%9A%E5%BE%8B
  
• 基尔霍夫定律 https://baike.baidu.com/item/%E5%9F%BA%E5%B0%94%E9%9C%8D%E5%A4%AB%E5%AE%9A%E5%BE%8B/15860277
  
• 色温 https://baike.baidu.com/item/%E8%89%B2%E6%B8%A9
  

https://baike.baidu.com/item/CIF%E6%A0%87%E5%87%86%E7%85%A7%E6%98%8E%E4%BD%93

  PointLight(const Transform &LightToWorld,
             const MediumInterface &mediumInterface, const Spectrum &I)
    : Light((int)LightFlags::DeltaPosition, LightToWorld, mediumInterface),
    pLight(LightToWorld(Point3f(0, 0, 0))),
    I(I) {}

严格来说使用辐射率来描述从点光源到达某个点的光是不正确的。辐射强度才是合适的单位。然而，我们滥用了术语，使用 Sample_Li 来计算从各种不同类型的光源到达某点的照明，将辐射强度除以距离的平方来转换单位。

• 辐射率 https://zh.wikipedia.org/wiki/%E8%BE%90%E5%B0%84%E7%8E%87
  
• 辐射强度 https://baike.baidu.com/item/%E8%BE%90%E5%B0%84%E5%BC%BA%E5%BA%A6
  

在渲染中，连续型随机变量比离散型随机变量更常用，其取值范围为连续区域。例如，实数、单位球上的方向，场景中物体的表面。

经典均匀分布随机变量为一个特别重要的随机变量，我们将其记为ξ，其等概率地在[0,1)区域上取值。
之所以重要有以下两个原因：

1. 在软件中易于生成符合该分布的随机变量
   
2. 可以通过先生成经典均匀分布的样本，然后再执行合适的变换来得到任意其他分布的样本。
   

概率密度函数(probability density function)描述了随机变量取特定值的概率。其为随机变量累计分布函数的导数。p(x) = dP(x)/dx.

./PhysicallyBasedRendering/00_13_uniform_pdf_cdf.ggb

• https://baike.baidu.com/item/%E5%9D%87%E5%8C%80%E5%88%86%E5%B8%83/954451
  

函数 f 的期望值 Ep[f(x)]定义为该函数在某个分布上的平均值。蒙特卡洛积分需要计算任意积分函数的期望值。

cos(x)在 0-π定义域内的期望值如下，其中 p 为均匀分布：

函数的方差为函数和该函数期望偏差平方的期望，其形式化定义如下：

inversion method 使用一个或多个均匀随机变量，并将其映射为目标分布的随机变量。下面使用离散的例子来解释其过程：

对于一些函数，可能无法对其积分获得 CDF，或者无法求得其 CDF 的反函数。此时可以使用 Rejection Method 来按照某分布函数来生成样本。Rejection Method 的本质是投飞镖方式。假设我们需要从 f(x)分布抽取样本，我们有 PDF p(x)满足 f(x)<cp(x)，而我们知道如何从 p(x)分布抽取样本。那么 Rejection Method 操作如下：

  // loop forever
  while(true)
    // 从 p 分布生成样本 X
    X = gen_sample(p);
    // 如果(X, ξ c p(X)) 在 f(X)下面，则接受生成的X样本，否则拒绝X样本
    if (ξ < f(X)/(c p(X))) then
      return X

该方法的效率依赖于 cp(x)和 f(x)靠近的程度。

Point2f RejectionSampleDisk(RNG &rng) {
    Point2f p;
    do {
        p.x = 1 - 2 * rng.UniformFloat();
        p.y = 1 - 2 * rng.UniformFloat();
    } while (p.x * p.x + p.y * p.y > 1);
    return p;
}

Metropolis 算法从函数 f 生成一组样本 Xi，函数 f 定义在任意维度状态空间 \(\Omega\) 上，并且返回实数值。选中第一个样本 X0 ∈ \(\Omega\) 后，在生成后续样本 Xi 时，对 X(i-1)进行随机突变得到一个建议样本 X'。若该突变被接受，则 Xi = X'，否则 Xi = X(i-1)。当选择从一个状态转换为另一个状态时需要满足一些条件，最终 Xi 的分布会达到一个均衡分布；这个分布为静态分布。在极限情况下，这组样本 Xi ∈ \(\Omega\) 的分布恰好符合 f(x) 的概率密度函数:

为了生成分布正确的样本，必须生成合适的突变，并且按照一些约束接受或拒绝突变。假设我们有一个突变方法，其可以将一个给定状态 X 转化为一个新的状态 X'。我们必须可以计算一个转换函数 T(X→X')，其可以给出从 X 突变为 X'的概率密度。
给定一个转变函数，可以定义一个 acceptance probability a(X→X') 其给出接受从 X 到 X' 突变的概率。如果分布已经均衡，则满足如下条件，这个属性被称为 detailed balance.

f 和 T 是设定好的，从上面公式可以得出 a 的定义。特别地，为了最大化达到均衡的速率，a的定义应该为：

接受概率的推导可以参考下面文章描述

• 随机采样方法整理与讲解（MCMC、Gibbs Sampling 等）https://www.cnblogs.com/xbinworld/p/4266146.html
  

Metropolis sampling 算法的伪代码如下：

  // 选定第一个样本
  X = X0;
  // 生成后续样本
  for(int i=1; i<n; i++)
  {
      // 对当前样本进行突变
      X' = mutate(X);
      // 计算接受当前样本的概率
      a = accept(X, X');
      // 比较转换概率和接受概率
      if (random() < a)
      {
          X = X'
      }
      record(X);
  }

因为 Metropolis 算法自然地避免了 \(\Omega\) 空间上部分区域，这些区域上 f(x)的值相对比较小，只有很少样本在该区域被累加。为了获取 f(x)在这些区域的更多信息，可以使用数学期望技术来增强 Metropolis 算法。增强后的 Metropolis 算法伪代码如下：

  // 选定第一个样本
  X = X0;
  // 生成后续样本
  for(int i=1; i<n; i++)
  {
      // 对当前样本进行突变
      X' = mutate(X);
      // 计算接受当前样本的概率
      a = accept(X, X');
      record(X, 1-a);
      record(X', a);
      // 比较转换概率和接受概率
      if (random() < a)
      {
        X = X'
      }
  }

通常希望突变有合适的大的变化而不是比较小的变化。这样可以快速的探索整个状态空间，小的突变会导致局限在状态空间的一个小的区域内。然而，当 f(X)的值在当前样本 X 处相对大时，很多建议的突变会被拒绝(从上节的接受概率的公式可以看出，此时的接受概率会变低)。我们希望避免很多样本是相同的，这样同样是为了更好地探索状态空间。在 f 的值相对大的地方，小的突变可能是建议的样本 X'，其可以得到更高的接受度。
因此，一个有效的突变方式是对当前样本 X 进行随机扰动。如果样本 X 表示一个实数向量(x0,x1,…)，则可以对该向量的某些维度或所有维度进行扰动。可以利用如下算法进行扰动：
xi' = (xi ± sξ) mod 1; 该方法是均匀对称的，因此不需要计算转换密度。
另一个相关的突变方式是，直接忽略当前的样本，使用均匀随机值生成一个新的样本：xi = ξ.完全随机地生成一个新的样本可以保证我们不会卡在状态空间的某个部分。一般地，能够到达状态空间的所有状态是必要的，这个属性被称为可遍历性。特殊地，为了保证遍历性，需要满足 T(X->X') > 0。
另一种方式是，使用利用概率密度函数进行突变，该 PDF 和被采样函数的部分是匹配的。如果我们有一个 PDF p(x),其和 f 的某个分量相同。此时的转换函数为：T (X→X') = p(X')。也就是说，当前的状态 X 对于转换密度的计算是无关的，转换到 X'的密度值依赖于 X'，和当前的状态无关。

如何计算初始的 X0 样本？如果使用一个不符合 f 分布的样本会导致 start-up bias.
一个通用的解决方案是，以一个任意的状态，执行 Meropolis 采样算法多次，忽略其生成的样本，假设前面的计算已经得到了一个合适的 X，然后继续执行该算法。该方式有两个缺点：1 消耗比较高，其需要忽略一些样本。2. 需要忽略多少样本只能靠猜测。
另一种方式是，如果另一种采样方法可用，我们可以使用任意密度函数 X0~p(x)来得到 X0.我们从状态 X0 开始 Markov chain，但是我们对所有生成的样本分配如下权重：w=f(X0)/p(X0).该方法可以以可预见的方式完全消除 start-up bias。
如果对于我们选取的 X0, f(X0)=0.此时，所有的样本的权重都为 0.为了避免该错误，我们可以采样 N 个样本，Y1, …,Yn,定义每个样本的权重为 wi=f(Yi)/p(Yi). 然后，从 Yi 中选择 X0，计算 Yi 的平均权重 w，将其作为 Metropolis 算法生成样本的权重。

多维的情况如下：

行列式是线性变换的伸缩因子。行列式=1，图形面积不变；行列式>1，图形被放大；0<行列式<1，图形被缩小；行列式<0，图形镜像（改变了基的左右手法则）。
因此，使用雅可比行列式缩放概率密度，可以表示坐标变换引起的概率密度变换。

• 行列式的本质是什么？ https://www.zhihu.com/question/36966326
  

以两维极坐标为例，假如我们从某个密度函数 p(r, θ)抽取样本,其对应的 p(x,y) 密度函数是什么？

雅可比矩阵表示的是从(x,y)坐标系到(r, θ)坐标系的变换。
雅可比行列式表示从(x,y)坐标系变换为(r, θ)坐标系后，空间的缩放系数。
(x, y) -> (r, θ) scaleFactor = r
p(x,y)r = p(r, θ)

Tips:
雅可比行列式为负数并不表示空间缩小了，而是说空间镜像了。0<行列式绝对值<1 才表示空间缩小了。

：
下面为三维球坐标系的情况：
方法 1：

方法 2：

• 球面立体角 https://baike.baidu.com/item/%E7%AB%8B%E4%BD%93%E8%A7%92
  

散射方程中 BSDF 和入射辐射率的乘积会再乘一个 cosine 项，因此生成的方向靠近半球顶端的可能性更大一些会很有用。
球面立体角 w 的密度函数 请参考链接 No description for this link , 下面为包含 cosine 项的球面立体角密度函数的推导：

下面使用了 Malley 方法生成 cosine 权重的立体角分布：

下面是我自己的推导：

对于 SpotLight 和基于 Sphere 的区域灯光来说，在一个 cone 内进行均匀采样是很有用的。

球面立体角 w 的密度函数 请参考链接 No description for this link。

• https://mathworld.wolfram.com/Cone.html
  

下面使用面积为 1/2 的等腰直角三角形推导对三角形采样的方法，其使用了重心坐标，通用于一般的三角形。

从离散的 2D 分布采样适用于生成贴图或环境贴图的样本。

再次考虑只有直接光照的反射问题，忽略像素过滤后，该问题可以写成如下形式：

通常会生成 N 个样本来计算积分，每个样本由一个像素位置(x,y)一个朝向光源的方向 w 组成。如果场景中由很多个光源，或者有一个区域光源用于产生软阴影，此时可能需要很多样本才能达到可接受的效果。不幸的是每个样本需要追踪两个射线，一个射线为像素(x,y)到看见的点，另一个阴影射线从看见的点到光源。
假设 N=100,则有 200 个射线要追踪：100 个摄像机射线，100 个阴影射线。100 个摄像机射线太多了，其对结果影响不大。使用 Splitting 方法，我们只需要取 5 个像素样本但对于每个像素样本取 20 个光照样本。

分层采样将整个积分区域 A 划分为 n 个不重叠的区域 A1,A2,…,An.每个区域被称为一个层。 从 A 中抽取样本时，我们按照每个层中的密度函数 pi，从每个 Ai 中抽取对应的 ni 个样本。
超级采样一个像素，就是分层采样的一个简单例子。将一个像素周围区域划分为 kxk 个格子，从每个格子中抽取一个样本。这比随机抽取 k*k 个样本要好，因为完全随机抽取的样本可能会聚集在一起。
选择一个不分层的样本等价于按照离散概率分布随机选择一个层 I，该离散概率分布由 vi 定义，然后在 Ai 层随机选择一个样本 X。这种情况下，X在 I 条件下进行选择，由条件概率可得其方差。

分层采样的主要缺点为随着维度提高，需要的样本数量变多。
一种方案是只对一些维度进行独立的分层，然后随机地将不同维度的样本关联起来，如：Stratified Sampling 部分描述。
另一种方案是使用 Latin Hypercube sampling，不论维度多少，其可生成任意数量样本。不幸的是，其无法像分层采样那样可以有效减小方差。但是其比均匀随机采样效果要好。

低差异采样技术被称为 Quasi Monte Carlo.拟蒙特卡洛技术的关键是其用精心设计的确定性算法生成低差异点集，用这些点集代替标准蒙特卡洛方法中的伪随机数。
拟蒙特卡洛比标准蒙特卡洛有渐进更快的收敛速度。使用标准蒙特卡洛的大多数技术都可以用拟蒙特卡洛来代替，有一些则不行(例如：rejection sampling)。图形学中，对于不连续的积分渐进收敛速度通常是无法接受的，这种情况下，拟蒙特卡洛要比标准蒙特卡洛表现要好。

• https://baike.baidu.com/item/%E6%8B%9F%E8%92%99%E7%89%B9%E5%8D%A1%E7%BD%97%E6%96%B9%E6%B3%95
  

pbrt 中应用分层采样或 low-discrepancy 采样时，会在[0,1)^2 空间上生成一系列样本，然后再利用 13.5 和 13.6 节中的转换方法将这些样本转换到需要的空间。转换过程中需要保持样本的分层属性，也就是说靠近的样本应该映射为物体表面上靠近的位置，分离开的样本应该映射为分离开的位置。如果无法保持该性质，则分层的好处就丢失了。

Multiple Importance Sampling 的基本理念为，估算积分时，我们可以从多个样本分布抽取样本，多个样本中至少有一个会和积分的的形状相匹配，即使我们不知道匹配的样本是哪个。MIS 提供了一种方法给样本分配权重，其可以减少方差峰值。只考虑不常见的特殊情况的特化采样程序是被鼓励的，因为当这些情况发生时，其可以减少方差，而且只是用先对少的开销。
权重函数考虑生成一个样本 Xi 或 Yj 所有不同的方式，而不只是特定的实际用到的一个。平衡启发式是权重函数的一个好的选择。其可以为样本分配权重减少方差。

  virtual Spectrum BxDF::Sample_f(const Vector3f &wo, Vector3f *wi, const Point2f &sample, Float *pdf, BxDFType *sampledType = nullptr) const;
  virtual Float BxDF::Pdf(const Vector3f &wo, const Vector3f &wi) const;

考虑间接光照时，在 BxDF::Sample_f 方法中，首先生成入射方向样本 wi，然后调用 BxDF::f 方法返回 wi 入射方向，wo 出射方向下对应的辐射率。例如：整个过程为 WhittedIntegrator::Li -> SamplerIntegrator::SpecularReflect -> BSDF::Sample_f -> BxDF::Sample_f -> BSDF::f。 Tips: 考虑光在表面的一次反射，其实就是包含了一次间接光照。

BxDF::Sample_f 方法按照分布选择一个方向，该分布和对应的散射方程相似。该方法用于在完美镜面表面处查找反射和透射射线。
BxDF::Sample_f 会取[0,1)^2 范围内的两个样本值，这两个样本值可以通过分层采样或低差异采样技术来生成。这两个样本值可用于基于 inversion method 的采样算法。
BxDF::Sample_f 方法返回给定方向的 BSDF 值，以及采样的方向 wi 和 wi 方向对应的概率密度函数 p(wi)的值 pdf。pdf 值为立体角，wi 和 wo 在标准反射坐标系下。

BxDF::Pdf() 方法返回给定入射方向和出射方向对应的概率密度 pdf。该方法对于 multiple importance sampling 非常有用。

基于微表面的反射模型会基于一个微表面分布 D(wh)，每一个微表面都展现出完美的镜面反射或透射。D(wh)函数对 Torrance-Sparrow BSDF 函数的形状起决定性作用，因此按照 D(wh)函数的分布进行采样会非常高效。使用这种方式，首先从微表面分布 D(wh)生成一个微表面的朝向样本，然后入射方向通过镜面反射或透射公式来获得。

// 实践中，采样可见微表面区域比采样整个分布要更高效
// 下面变量用于指定是否只采样可见微表面区域
  const bool MicrofacetDistribution::sampleVisibleArea;

Vector3f Sample_wh(const Vector3f &wo, const Point2f &u) const;

下面为 BeckmannDistribution::Sample_wh 对于各项同性情况下，在整个分布上进行采样的原理：

下图为从倾斜方向看可见微表面分布和整体微表面分布的区别：

下面为 从 wh 的分布概率密度函数推导 wi 的分布概率密度函数：

FresnelBlend 类为 diffuse 和 glossy 的混合。采样该 BRDF 的一种直接的方式为同时对 cosine-weighted 分布和 microfacet 分布进行采样。

前面章节使用 Dirac delta 分布来定义 specular 反射的 BRDF 和 specular 透射的 BTDF。对于当前的采样框架其非常适用，Dirac delta 分布在 x!=0 时，其值为 0，其只有一个样本，因此为该分布生成样本非常简单。但是，确定其 PDF 则不是很简单。严格来说，delta 分布不是一个真正的函数，其必须定义为另一个函数的极限(如 面积为 1 的 box 函数，其宽度无限接近于 0 对应函数，此时δ(0)趋近于无穷大)，当然，计算 PDF 时，返回一个无穷大或很大的值，在渲染时，不会得到正确的结果。

基于上面的分析，Sample_f()中分子分母中，因此返回的 pdf 为 1，而 Pdf()返回的 pdf 为 0，表示其他方向的概率密度为 0.

FresnelSpecular 类同时封装了镜面反射和镜面透射，其通过电介质 Fresnel 项来调节镜面反射和镜面透射的量。例如，在掠射角度上，反射比较高，此时使用 MonteCarlo 采样时应该更大可能返回反射方向而不是透射方向，这种方式改进了 MonteCarlo 的效果。

FresnelSpecular::Sample_f() 当样本为 Specular Reflection 时，pdf=F，当样本为 Specular Transmission 时，pdf=1-F。F 为 Fresnel 反射量，计算 F 值的函数为 FrDielectric函数 FrConductor函数。

  FresnelSpecular(const Spectrum &R, const Spectrum &T, Float etaA,
                  Float etaB, TransportMode mode)
    : BxDF(BxDFType(BSDF_REFLECTION | BSDF_TRANSMISSION | BSDF_SPECULAR)),
    R(R),
    T(T),
    etaA(etaA),
    etaB(etaB),
    fresnel(etaA, etaB),
    mode(mode) {}


  Spectrum FresnelSpecular::Sample_f(const Vector3f &wo, Vector3f *wi,
                                     const Point2f &u, Float *pdf,
                                     BxDFType *sampledType) const {
    Float F = FrDielectric(CosTheta(wo), etaA, etaB);
    // 对反射进行采样
    if (u[0] < F) {
      // Compute specular reflection for _FresnelSpecular_

      // Compute perfect specular reflection direction
      *wi = Vector3f(-wo.x, -wo.y, wo.z);
      if (sampledType)
        *sampledType = BxDFType(BSDF_SPECULAR | BSDF_REFLECTION);
      *pdf = F;
      return F * R / AbsCosTheta(*wi);
    }
    else // 对透射进行采样
    {
      // Compute specular transmission for _FresnelSpecular_

      // Figure out which $\eta$ is incident and which is transmitted
      bool entering = CosTheta(wo) > 0;
      Float etaI = entering ? etaA : etaB;
      Float etaT = entering ? etaB : etaA;

      // Compute ray direction for specular transmission
      if (!Refract(wo, Faceforward(Normal3f(0, 0, 1), wo), etaI / etaT, wi))
        return 0;
      Spectrum ft = T * (1 - F);

      // Account for non-symmetry with transmission to different medium
      if (mode == TransportMode::Radiance)
        ft *= (etaI * etaI) / (etaT * etaT);
      if (sampledType)
        *sampledType = BxDFType(BSDF_SPECULAR | BSDF_TRANSMISSION);
      // 注意：透射的pdf = 1-F，这保证了能量守恒
      *pdf = 1 - F;
      return ft / AbsCosTheta(*wi);
    }
  }

hemispherical-directional reflectance半球方向反射
hemispherical-hemispherical reflectance 半球半球反射

BSDF 包含了所有的 BxDF。BxDF* BSDF::bxdfs[MaxBxDFs] 成员变量存储所有 bxdf 的引用。

BSDF::Pdf

和完美的镜面反射和镜面透射一样，被定义为 delta 分布的光源也适用于当前采样框架，但是需要注意返回的辐射率和 PDF 中所隐含的 delta 分布，大多数情况下，这些 delta 分布会被约掉，但是 Multiple Importance Sampling 必须注意这种情况。
PointLight 是由 delta 分布来描述的，其只会从单个方向照明一个点。由于是 delta 分布，因此 PointLight::Pdf_Li()返回 0.

pbrt 中，area lights 都有一个 Shape 对象来指定区域光的形状，为了从这样的区域光采样入射照明，需要在区域光形状上生成样本。因此，pbrt 为 Shape 类增加了采样方法。
有两个采样方法，第一个方法如下：

// 该方法除了初始化采样点的位置p，法线n以外，其还需要处理浮点值的舍入错误
virtual Shape::Interaction Sample(const Point2f &u) const = 0;

// 形状上的采样几乎都是均匀的，因此其pdf为面积的倒数
virtual Float Shape::Pdf(const Interaction &) const { return 1 / Area(); }

第二个采样方法如下，该方法对于光照来说特别有用，因为调用者可以传入被照亮的点，其允许区域光的形状实现保证只采样对于被照亮点可见的部分。

virtual Interaction Shape::Sample(const Interaction &ref, const Point2f &u) const {return Sample(u);}
// 默认实现中，将定义在区域上的密度转化为定义在立体角上的密度
// 给定一个交点和wi，若从交点发出wi方向的射线和区域光形状不相交，则pdf为0，若相交其计算公式如下图：
virtual Float Shape::Pdf(const Interaction &ref, const Vector3f &wi) const;