[论文阅读] Embree: A Kernel Framework for Efficient CPU Ray Tracing

Embree 结构

Embree 的结构：计算硬件 – Embree Common Infrastructure – Embree Kernel Layer – Embree API – 渲染器。

Embree 的关键组件为：

BVH 构建核心（BVH Construction Kernels）；
光线追踪核心（Traversal and Intersection Kernels）：
- 单光线三角面求交（Single-Ray Triangle Intersection）；
- 单光线 BVH 求交（Single-Ray BVH Traversal）；
- ······

光线追踪核心

光线追踪核心（Traversal and Intersection Kernels）。

Embree 支持的 BVH 格式与基元类型、矢量化方法和指令集架构（ISA）很多，也可以只编译面剔除（face culling）、光线遮罩（ray masks）部分。

单光线三角面求交

单光线三角面求交（Single-Ray Triangle Intersection）。

对于不同处理器架构有不同的并行实现：

4-wide SSE 使用 triangle1n、triangle4n 类型。
8-wide AVX / AVX2 使用 triangle4n、triangle8n 类型。
16-wide IMCI 上可使用 triangle16n 类型······

会通过 __AVX__、__AVX512F__ 启用不同的基类，有对应的实现。

单光线 BVH 求交

单光线 BVH 求交（Single-Ray BVH Traversal）。

对于四叉BVH（Quad-BVH）结构，支持向量宽度为 4 的处理器架构可以实现高效的遍历。Embree 没有对 AVX / AVX2 做映射到 8 宽 SIMD 的优化，只是改进了运算细节（融合加法和乘法）。

光线包 BVH 遍历

光线包 BVH 遍历（Packet Traversal and Intersection）。

Embree 使用经典光线包求交，而不是单程序多数据（SPMD）的方法。在每个 BVH 的节点上通过单指令多数据（SIMD）地分摊标量计算，Embree 在这一步对于不同 ISA 可以生成最佳指令组合。

混合 BVH 遍历

混合 BVH 遍历（Hybrid Traversal and Intersection）。

在混合相干（coherent）光线与非相干光线时，如果支持光线包和单光线的切换，可以加速求解。但针对光线包优化的 BVH 的内存存储顺序未必对于单光线是最佳的（反之亦然）。Embree 设定了阈值用于切换两种方式的使用。

//kernels/bvh/bvh_intersector_hybrid.h#L33
static const size_t switchThresholdIncoherent =
    (K == 4)    ? 3
    : (K == 8)  ? ((N == 4) ? 5 : 7)
    : (K == 16) ? 14 : // 14 seems to work best for KNL due to better ordered chunk traversal
    0;

//kernels/bvh/bvh_intersector_hybrid.cpp#L679
/* switch to single ray traversal */
#if (!defined(__WIN32__) || defined(__X86_64__)) && ((defined(__aarch64__)) || defined(__SSE4_2__))
#if FORCE_SINGLE_MODE == 0
if (single)
#endif
{
  size_t bits = movemask(active);
#if FORCE_SINGLE_MODE == 0
  if (unlikely(popcnt(bits) <= switchThreshold))
#endif
  {
    for (; bits != 0;) {
      const size_t i = bscf(bits);
      if (occluded1(This, bvh, cur, i, pre, ray, tray, context)) set(terminated, i);
    }
    if (all(terminated)) break;
    tray.tfar = select(terminated, vfloat
<K>(neg_inf), tray.tfar);
    continue;
  }
}
#endif

BVH 构建核心

BVH 构建核心（BVH Construction Kernels）。

桶式表面积启发式 BVH 构建

桶式表面积启发式 BVH 构建（Binned SAH BVH Construction）。

Embree 构建了针对优化遍历成本的 BVH，使用表面积启发式（Surface Area Heuristic, SAH）和空间分割的对象分区的技术，在构建时对节点进行估价，产生最低预估成本的分区。

表面积启发式（Surface Area Heuristic, SAH）

假设：光线进入体积的概率与其表面积成正比。则对于划分 $V\to {L,R}$ ，记两边三角形数量为 $N_L, N_R$ ，相关体积为 $V_L,V_R$ ，体积 $V$ 的表面积为 $SA(V)$ ，可以量化预期遍历成本为： $\text{Cost}(V\to{L,R})=K_T+K_I\left(\frac{SA(V_L)}{SA(V)}N_L+\frac{SA(V_R)}{SA(V)}N_R\right)$ 其中 $K_T$ 表示遍历子节点的开销， $K_I$ 表示求三角形交点的成本。

对于 KD 树来说，分区垂直于特定坐标轴，可以证明只有 $6N$ 个合理的分割位置。而对于 BVH ，分割会有 $2^N-2$ 种方案，所以需要设计划分的方案，在 PBRT 的实现里，作者将分割维度的 span 分成 12 段，取期望开销最小的一种。

Embree 的 SAH BVH 具体实现在 /kernels/builders/bvh_builder_sah.h 中的 recurse 方法：

  const float leafSAH  = cfg.intCost*current.prims.leafSAH(cfg.logBlockSize);
  const float splitSAH = cfg.travCost*halfArea(current.prims.geomBounds)+cfg.intCost*split.splitSAH();

  /*! compute leaf and split cost */
  if (current.prims.size() <= cfg.minLeafSize || current.depth+MIN_LARGE_LEAF_LEVELS >= cfg.maxDepth || (current.prims.size() <= cfg.maxLeafSize && leafSAH <= splitSAH)) {
    heuristic.deterministic_order(current.prims);
    return createLargeLeaf(current,alloc);
  }
  // ... some codes here ...

  /*! split until node is full or SAH tells us to stop */
  while (numChildren < cfg.branchingFactor) {
    /*! find best child to split */
    float bestArea = neg_inf;
    ssize_t bestChild = -1;
    for (size_t i = 0; i < numChildren; i++) {
      /* ignore leaves as they cannot get split */
      if (children[i].prims.size() <= cfg.minLeafSize) continue;

      /* find child with largest surface area */
      if (halfArea(children[i].prims.geomBounds) > bestArea) {
        bestChild = i;
        bestArea = halfArea(children[i].prims.geomBounds);
      }
    }
    if (bestChild == -1) break;

    /* perform best found split */
    BuildRecord& brecord = children[bestChild];
    BuildRecord lrecord(current.depth + 1);
    BuildRecord rrecord(current.depth + 1);
    auto split = heuristic.find(brecord.prims, cfg.logBlockSize);
    heuristic.split(split, brecord.prims, lrecord.prims, rrecord.prims);
    children[bestChild] = lrecord;
    children[numChildren] = rrecord;
    numChildren++;
  }

桶式 BVH 构建（Binned BVH Building）

对质心边界框宽度最大的轴进行划分，均匀分为 $K$ 个相等的桶。根据三角面的重心分到对应组内，但三角形通常会超出桶的边界，这里对所有投影到的桶的数量 $n_i$ 与其边界框 $B_i$ 进行追踪，然后将 $n_i$ 应用到对应桶的初始边界，然后进行合并边界。同时为了防止恶化，加入超参数用来调整该层构建。

上面 auto split = heuristic.find(brecord.prims, cfg.logBlockSize); 这里的启发式寻找使用了桶式 BVH 构建，Embree 实现如下：

//kernels/builders/heuristic_binning.h#L339
__forceinline Split best(const BinMapping
<BINS>& mapping, const size_t blocks_shift) const {
  // ... some codes here ...
  vuint4 blocks_add = (1 << blocks_shift)-1;
  for (size_t i = 1; i < mapping.size(); i++, ii += 1) {
    count += counts(i - 1);
    bx.extend(bounds(i - 1, 0)); float Ax = expectedApproxHalfArea(bx);
    by.extend(bounds(i - 1, 1)); float Ay = expectedApproxHalfArea(by);
    bz.extend(bounds(i - 1, 2)); float Az = expectedApproxHalfArea(bz);
    const vfloat4 lArea = vfloat4(Ax, Ay, Az, Az);
    const vfloat4 rArea = rAreas[i];
    const vuint4 lCount =
        (count + blocks_add) >> (unsigned int)(blocks_shift);
    const vuint4 rCount =
        (rCounts[i] + blocks_add) >> (unsigned int)(blocks_shift);
    const vfloat4 sah = madd(lArea, vfloat4(lCount), rArea * vfloat4(rCount));

    vbestPos = select(sah < vbestSAH, ii, vbestPos);
    vbestSAH = select(sah < vbestSAH, sah, vbestSAH);
  }
  // ... some codes here ...
  return Split(bestSAH,bestDim,bestPos,mapping);
}

对于 bin 的划分可以在 /kernels/builders/ 下 heuristic_ 开头的文件找到，这里不再赘述。

空间分割 BVH（Split Bounding Volume Hierarchy, SBVH）

将 KD 树和 BVH 更深度地结合，以引用的方式允许将一个物体分到两侧（将这个代价添加到成本函数中）。

多线程 BVH 构建

多线程 BVH 构建（Morton-Code BVH Construction）。

高维坐标映射（Morton Code）

对于高维整点坐标 $(x,y,z)$ ，写成二进制形式： $\left(\overline{x_{n} x_{n-1} \cdots x_{1} x_{0}}, \overline{y_{n} y_{n-1} \cdots y_{1} y_{0}}, \overline{z_{n} z_{n-1} \cdots z_{1} z_{0}}\right)$ ，然后交错排列，得到 $\overline{z_{n}y_{n}x_{n}z_{n-1}y_{n-1}x_{n-1}\cdots z_{1}y_{1}x_{1}z_{0}y_{0}x_{0}}$ ，再转换成无符号整数，就得到了这一坐标全序的 Morton code 表示。Morton code 有很好的局部性，以下图的二维 Morton code 为例：

(a) 根据最高位划分、(b) 根据次高位划分、(c) 根据一、二位划分、(d) 根据一、三位划分。

实际计算 Morton code 时，我们将包含所有物体的全局包围盒在各个象限上先归一化再 $1024$ 等分，将每个物体的包围盒中心坐标转换到 $[0,2^{10})$ 上去，然后就可以得到 $30$ 位的 Morton code，正好可以放进一个 uint32_t 里。

得到所有 Morton code 之后，我们对这些点进行一次基数排序，然后按照固定步长将所有物体划分为若干个组，就得到了彼此不交而组内点的距离也比较接近的很多个不同组，组内建 BVH 的操作在组间可以并行。

Embree 的 Morton code 存储在 64 位值中，具体的并行算法为：初始所有线程协作自同一节点划分三角形，当子节点的数量大于线程数量时，将节点划给线程，在三角形数少于 4 个时在此叶子节点处中止。

Embree 的 Morton 代码实现为：

//kernels/builders/bvh_builder_morton.h#L295
/* recalculate morton codes */
MortonCodeMapping mapping(centBounds);
parallel_for(current.begin(), current.end(), unsigned(1024), [&](const range
<unsigned>& r) { for (size_t i = r.begin(); i < r.end(); i++) { morton[i].code = mapping.code(calculateBounds(morton[i])); } });

/*! sort morton codes */
#if defined(TASKING_TBB)
tbb::parallel_sort(morton + current.begin(), morton + current.end());
#else
radixsort32(morton + current.begin(), current.size());
#endif

线程方面代码实现为（singleThreadThreshold 被初始化为 (size_t)1024）：

//kernels/builders/bvh_builder_morton.h#L413
if (current.size() > singleThreadThreshold) {
  /*! parallel_for is faster than spawning sub-tasks */
  parallel_for(size_t(0), numChildren, [&](const range
<size_t>& r) {
    for (size_t i = r.begin(); i < r.end(); i++) {
      bounds[i] = recurse(depth + 1, children[i], nullptr, true);
      _mm_mfence();  // to allow non-temporal stores during build
    }
  });
}

Embree 应用

Unreal Engine 5 中的 Lumen 可以实现全动态全局光照（Fully Dynamic Global Illumination），中间使用辐照度算法（Radiosity）来生成间接光（Indirect Lighting），理论上不是光线追踪而是光线投射（Ray Casting）。

传统的辐照度算法方法需要将场景划分为 Patch，而 Lumen 已经拥有了粗粒度的全局距离场（Global DF）以及粗粒度的体素光照（Voxel Lighting），因此可以直接从表面缓存（Surface Cache）上射出光线，与全局距离场进行光线追踪求交，交点采样上一帧的体素光照后转换为二阶球谐，最后再根据 Normal 计算漫反射（Diffuse Transfer）球谐系数，计算出最终的间接辐照度（Indirect Radiance）。

这个间接辐照度也保存在辐照度缓存中，称为间接光，并与直接光合并计算得到最终的渲染光照，而下一帧的体素光照又来自于这一帧的辐照度缓存，因此后续所有的 Lumen 光照流程自然具有了间接光照。

表面缓存是算法的核心部分，Lumen 创造了 Card 的概念，用于定义面在表面缓存中存储的信息。

Lumen 会从多个角度捕获每个网格体的材质属性。这些捕获位置（这里叫做：卡（Cards））是针对每个网格体脱机生成的。

这里的”多个角度”就是指 6 个轴对齐方向（上、下、左、右、前、后），”捕获”指的是在 6 个轴对齐方向上以正投影方式光栅化 Mesh，获取 Mesh 对应的材质中的各个 Material Attributes（Albedo、Opacity、Normal、Emissive）后存储到表面缓存上。如果网格体具有复杂结构还会生成多个 Card。

在访问表面缓存时也需要通过 Card 将 Hit Point 进行物理地址转换到表面缓存图集空间中再执行采样。因此 Card 并非是单一的位置数据，其核心数据是一个局部轴对齐的边界框（Bounding Box）。

在初始化阶段还是需要光线追踪来进行处理表面缓存，也就是计算 Card，需要对体素化的 Mesh 按 Cube 轴对齐方向执行光线步进生成 Surfel，记录 Surfel 的信息，最后根据这些 Surfel 计算每方向 Card 的边界，具体流程（MeshCardRepresentationUtilities.cpp 下的 GenerateSurfelsForDirection 方法）为：

遍历当前轴对齐方向的 XY 平面（注意每个方向的 Basis 不同），沿着 Z 方向在每个 Cell 中发出 32 条射线进行光线求交，如果相交则增加一个面元采样点（Surfel Sample），记录交点、Normal、交点所在的 Cell Z 和前一个交点的 Cell Z 记为 MinRayZ，直到超出边界结束。
将面元采样点列表按照 Z 升序排序，如果 Z 相同则按 MinRayZ 降序排序、统计每个 Cell 的面元采样点数量以及相邻 Cell（Z 方向）之间的面元采样点数量差。
遍历每个 Cell 判断其对应的 Surfel 是否有效，如果有效则增加真正的 Surfel，记录 Cell 坐标以及 MinRayZ 的中位数。有效性判断是通过光线追踪判断是否在 Geometry 内部，不在内部则认为有效。

参考文献

Embree: A Kernel Framework for Efficient CPU Ray Tracing: 提出 Embree 的论文。
embree/embree: Embree ray tracing kernels repository：Embree 代码。
Shallow bounding volume hierarchies for accelerated ray tracing: 验证 Quad-BVH 在 4 宽架构下的优化。
On fast Construction of SAH-based Bounding Volume Hierarchies: 将 SAH 与 Bin-Build 应用到 BVH。
Spatial Splits in Bounding Volume Hierarchies: SBVH，空间分割BVH。
渲染拾遗：BVH: 对 SAH 与 Morton Code 有形象地解释。
Lumen: Real-time Global Illumination in Unreal Engine 5 – SIGGRAPH 2022：Lumen 原论文。
Lumen Technical Details in Unreal Engine：Lumen 的技术细节。
游戏引擎随笔 0x30：UE5 Lumen 源码解析（二）Surface Cache 篇：详细描述 Card 的功能与实现。