在光线追踪中优化BVH遍历

Question

在我的CPU光线追踪器(实际上是路径追踪器)中，大部分CPU时间都花费在BVH遍历函数中。根据我的性能分析器，75%的光线追踪时间都花费在这个函数及其调用的函数上，而35%的时间则花费在函数本身上。其他40%的时间则花费在它调用的不同相交测试上。

基本上，代码通过所有与之相交的边界框和三角形进行深度优先搜索遍历。它使用静态分配的堆栈上的数组来保存要探索的节点(BVHSTACKSIZE设置为32，大多数情况下不需要该空间)，以便不会动态分配内存。然而，我觉得35%的时间在这里花费太高了。我已经花了一段时间来优化代码，目前已经是我能做到的最快速度，但它仍然是我的程序中最大的瓶颈。

有没有人有更好的优化建议？我已经有一个不错的BVH构造算法了，所以我认为使用不同的BVH也不会有任何加速效果。有没有人有关于如何在Mac上进行逐行性能分析的最佳方法？

供参考，在一个示例场景中，这段代码根据相交数量需要花费1微秒到40微秒不等的时间，而while循环则运行1到约400次迭代(也取决于相交数量)。

谢谢！

bool BVHAccel::intersect(Ray& ray) const {
  bool hit = false;

  BVHNode* to_intersect[BVHSTACKSIZE];
  int head = 0;
  to_intersect[head++] = root;

  while (head != 0) {
    assert(head < BVHSTACKSIZE);
    BVHNode* cur = to_intersect[--head];

    if (cur->bb.intersect(ray)) { // Does not modify the ray
      if (cur->isLeaf()) {
        for (const auto& primitive : cur->primitives) {
          hit |= primitive->intersect(ray); // Modifies the ray!
        }
      } else {
        to_intersect[head++] = cur->r;
        to_intersect[head++] = cur->l;
      }
    }
  }

  return hit;
}

bool BBox::intersect(const Ray& r) const {
  double txmin = (min.x - r.o.x) * r.inv_d.x;
  double txmax = (max.x - r.o.x) * r.inv_d.x;
  double tymin = (min.y - r.o.y) * r.inv_d.y;
  double tymax = (max.y - r.o.y) * r.inv_d.y;
  double tzmin = (min.z - r.o.z) * r.inv_d.z;
  double tzmax = (max.z - r.o.z) * r.inv_d.z;

  ascending(txmin, txmax);
  ascending(tymin, tymax);
  ascending(tzmin, tzmax);

  double t0 = std::max(txmin, std::max(tymin, tzmin));
  double t1 = std::min(txmax, std::min(tymax, tzmax));

  if (t1 < t0 || t0 > r.max_t || t1 < r.min_t) {
    return false;
  }

  return true;
}

void ascending(double& a, double& b) {
  if (a > b) {
    std::swap(a, b);
  }
}

- Leo Adberg

---

你能否发布一份完整的可运行代码？我可以大致了解你的代码是做什么的，但如果我们有一份可以复制粘贴并运行的代码，那将会更容易调试。 - Nelewout

我会稍微扩展一下这个例子，但不幸的是我不能发布一个自包含的例子，因为我没有编写周围代码的大部分并且不能发布它。 - Leo Adberg

但是可以将BVH视为二叉搜索树，它以深度优先的方式遍历。该树大约有10层，每个遍历的树节点都需要进行BBox交集计算，而叶节点则需要进行特殊的交集计算（例如三角形、球体）。绝大部分（约70%）的时间都花费在DFS遍历和BBox交集计算上，其中DFS遍历所需的时间最长。 - Leo Adberg

我会关注这个问题，并明天再考虑一下。感谢您对交集方法的澄清！ - Nelewout

没问题，谢谢你的帮助！ - Leo Adberg

你的原始类有边界球吗？在进行边界矩形检查之前，你应该能够快速地通过球来检查射线。 - 1201ProgramAlarm

Answer 1

你的代码似乎至少有一个问题。 复制 primitive 可能会是一项昂贵的操作。

bool BVHAccel::intersect(Ray ray) const {
  bool hit = false;

  BVHNode* to_intersect[BVHSTACKSIZE];
  int head = 0;
  to_intersect[head++] = root;

  while (head != 0) {
    assert(head < BVHSTACKSIZE);
    BVHNode* cur = to_intersect[--head];

    if (cur->bb.intersect(ray)) { // Does not modify the ray
      if (cur->isLeaf()) {
        for (const auto& primitive : cur->primitives) { // this code made a copy of primitives on every call!
          hit |= primitive->intersect(ray); // Modifies the ray!
        }
      } else {
        to_intersect[head++] = cur->r;
        to_intersect[head++] = cur->l;
      }
    }
  }

  return hit;
}

为什么需要修改光线的副本？

编辑1：我们可以假设BVHNode长这样吗？

constexpr auto BVHSTACKSIZE = 32;

struct Primitive;

struct BVHNode {
    std::vector<Primitive> primitives;
    AABB        bb;   
    BVHNode*    r = nullptr;
    BVHNode*    l = nullptr;

    bool isLeaf() const { return r == nullptr && l == nullptr; }
};

- Remo

---

你在两个方面都是完全正确的，结果我在试图使我的代码在SO上更易读时搞砸了（实际上它使用现有的迭代器来迭代原语）。然而，在真正的代码中它不会复制任何内容。感谢你发现了这一点！ - Leo Adberg

Answer 2

我看到你可以进行三个改进。

第一个大问题（很难）是你的代码中有许多条件分支，这肯定会减慢CPU速度，因为它无法预测代码路径（这在编译时也是正确的）。例如，我看到你首先相交，然后测试节点是否为叶子节点，然后与所有主体相交。能否先测试是否为叶子节点，然后再进行适当的相交？这将稍微减少分支。

其次，你的BVH的内存布局是什么？能否优化它以使遍历更加友好。你可以尝试查看遍历过程中发生的缓存未命中数，这将是你的内存是否具有适当布局的良好指标。虽然不是直接相关的，但了解你的平台和底层硬件是很好的。我建议阅读this。

最后，这是你在性能方面将产生最大影响的地方，使用SSE/AVX！通过对你的交叉代码进行一些重构，你可以同时相交四个边界框，从而在应用程序中获得良好的提升。你可以看看embree（英特尔跟踪器）的数学库做了什么。

此外，我刚刚发现你正在使用double。有什么原因吗？我们的路径追踪器根本不使用double，因为在渲染过程中没有真正需要那种精度的情况。

希望它有所帮助！

编辑：如果你想尝试一下，我已经制作了一个bbox交叉的sse版本。它部分基于我们的代码，但我不太确定它是否会起作用，你应该进行基准测试和测试！

#include <xmmintrin.h>
#include <emmintrin.h>
#include <smmintrin.h>

#include <cmath>
#include <limits>

constexpr float pos_inf = std::numeric_limits<float>::max();
constexpr float neg_inf = std::numeric_limits<float>::min();

size_t __bsf(size_t v)
{
  size_t r = 0; asm ("bsf %1,%0" : "=r"(r) : "r"(v));
  return r;
}

__m128 mini(const __m128 a, const __m128 b)
{
  return _mm_castsi128_ps(_mm_min_epi32(_mm_castps_si128(a),_mm_castps_si128(b)));
}

__m128 maxi(const __m128 a, const __m128 b)
{
  return _mm_castsi128_ps(_mm_max_epi32(_mm_castps_si128(a),_mm_castps_si128(b)));
}

__m128 abs(const __m128 a)
{
  return _mm_andnot_ps(_mm_set1_ps(-0.0f), a);
}

__m128 select(const __m128 mask, const __m128 t, const __m128 f)
{ 
  return _mm_blendv_ps(f, t, mask); 
}

template<size_t i0, size_t i1, size_t i2, size_t i3>
__m128 shuffle(const __m128 b)
{
  return _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(b), _MM_SHUFFLE(i3, i2, i1, i0)));
}

__m128 min(const __m128 a, const __m128 b) { return _mm_min_ps(a, b); }

__m128 max(const __m128 a, const __m128 b) { return _mm_max_ps(a, b); }

__m128 vreduce_min(const __m128 v)
{ 
  __m128 h = min(shuffle<1,0,3,2>(v),v);
  return min(shuffle<2,3,0,1>(h),h);
}

__m128 vreduce_max(const __m128 v)
{ 
  __m128 h = max(shuffle<1,0,3,2>(v),v);
  return max(shuffle<2,3,0,1>(h),h);
}

size_t select_min(__m128 valid, __m128 v)
{
  const __m128 a = select(valid, v, _mm_set_ps1(pos_inf));
  return __bsf(_mm_movemask_ps(_mm_and_ps(valid, (a == vreduce_min(a)))));
}

size_t select_max(const __m128 valid, const __m128 v)
{ 
  const __m128 a = select(valid, v, _mm_set_ps1(neg_inf));
  return __bsf(_mm_movemask_ps(_mm_and_ps(valid, (a == vreduce_max(a)))));
}

struct Ray
{
  vec3 o, inv_d;

  float min_t, max_t;
};

struct BBox
{
  vec3 min, max;

  bool intersect(const Ray& r) const;
};

bool BBox::intersect(const Ray& r) const
{
  const __m128 lowerSlab = _mm_mul_ps(_mm_sub_ps(max.m128, r.o.m128), r.inv_d.m128);
  const __m128 upperSlab = _mm_mul_ps(_mm_sub_ps(min.m128, r.o.m128), r.inv_d.m128);

  __m128 tmin = mini(lowerSlab, upperSlab);
  __m128 tmax = maxi(lowerSlab, upperSlab);

  reinterpret_cast<float*>(&tmin)[3] = r.min_t;
  reinterpret_cast<float*>(&tmax)[3] = r.max_t;

  const __m128 maskmin = _mm_castsi128_ps(_mm_cmpeq_epi32(tmin, tmin));
  const __m128 maskmax = _mm_castsi128_ps(_mm_cmpeq_epi32(tmax, tmax));

  const float tNear = abs(tmin[select_max(maskmin, tmin)]); // select the max non NaN value and ensure the result is positive using abs
  const float tFar  =     tmax[select_min(maskmax, tmax)]; // select the min non NaN value
  return tNear <= tFar;
}

- Thomas Caissard

---

谢谢！这些是很好的建议。我真的不知道如何减少分支的数量，对于isLeaf()，我仍然必须进行交集，所以我不知道我该如何改变它（而且在大多数迭代中它都是false，所以它可能会被正确预测）。BVH的内存布局有点随机。目前每个节点都是单独malloc的，但我会尝试将它们存储在数组中。最后，我被告知双精度浮点数与单精度浮点数一样快（只是在内存带宽方面稍慢，这里没有太多），但我会尝试自己检查一下。 - Leo Adberg

将所有节点粘贴到一个连续的数组中，在一堆运行中平均提高了约2％的速度，谢谢！ - Leo Adberg

你应该确定将节点存储在数组中或使用某种内存池，这确实可以提高缓存行的性能。关于float和double速度的比较，这是一个棘手的问题，它真正取决于你正在做什么。如果你缺少很多L1缓存，那么double会比float更耗费资源，因为硬件需要加载更多的数据到寄存器中。如果你决定使用SSE，那么float是一个真正的优势，因为CPU可以对4个操作进行计算，而double只能计算两个。 - Thomas Caissard

这也让我想到：你的交点代码中是否使用了标准数学函数？如果是，你应该考虑将它们更改为不太健壮但更快速的版本（我主要考虑exp和sqrt）。 - Thomas Caissard

到目前为止，最常见的交集是BBox交集，它只使用*，-，+，<，>。其他交集函数（三角形和球体）确实使用了更复杂的函数（除法，sqrt），但它们仅占执行时间的7％，所以我还没有费心去优化它们。 - Leo Adberg

哇，感谢提供SSE代码！我还应该提到，在实施CygnusX1的算法建议后，实现您的建议同时测试多个边界框变得更容易了。我编写了一段代码，可以同时处理两个（左右子节点），然后clang自动向量化了它：https://pastebin.com/DKsQ4tFH - Leo Adberg

Answer 3

我认为在你开始考虑硬件并挤出每个字节和周期之前，你可以进行算法改进。

我认为您感兴趣的是光线的第一个命中点。您没有对沿射线多次命中进行累加，对吗？（好像基元是半透明的或其他什么）。是的吗？

因此-如果以上内容属实，我认为您的遍历顺序存在问题。如果我正确理解了您的代码，您无条件地首先遍历cur->l，然后遍历cur->r。如果您的光线首先击中左侧框，则这很好。但是-它可能来自场景的另一侧-那么您将执行比所需更多的交点，从而沿着光线从后到前地遍历。

相反，您希望沿着光线从前到后遍历，希望您尽早命中某些东西，使您能够跳过大部分进一步的遍历测试。

幸运的是，在您的交集函数中已经计算出了与边界框的距离，您可以轻松地检查哪个优先相交。它只需要返回一个float交点距离而不是bool（例如在失败时返回+无穷大）。并且您需要在将东西推入to_intersect堆栈之前调用它。

因此，以伪代码的形式，我会像这样遍历：

while stack not empty:
    cur = pop from top of stack;
    //we already know that we want to enter this node!
    if cur is leaf:
        intersect primitives
    else:
        t_left = intersect bbox of cur->l
        t_right = intersect bbox of cur->r
        if both intersected:
            if t_left < t_right:
                push cur->r, cur->l in that order (so that cur->l will be on top)
            else:
                push cur->l, cur->r in that order (so that cur->r will be on top)
        else if one intersected:
            push only that one
        else:
            push nothing

一旦您实施了上述内容，您可能会注意到还有一个改进需要进行。如果在推入和弹出之间击中某些原始物体，并且您的ray.max_t被减少，则条件“我们已经知道要进入此节点”不一定成立。为了解决这个问题，您可以在堆栈中存储一对(t_enter，node)。然后，当您弹出时，您需要使用当前的ray.max_t再次检查t_enter。这可以节省高度为h的遍历检查，其中h是您的树的高度。

可能会出现的陷阱 - 您可能已经知道这一点，但以防万一 - 仅因为您从前面到后面遍历并不意味着您可以在找到第一个命中后立即终止遍历。BVH节点可能会重叠。这就是为什么要比较t_enter和ray.max_t的原因。

- CygnusX1

---

这是一个我没有想到的好主意！我稍后会尝试一下并告诉你加速效果。 - Leo Adberg

刚刚完成测试，它提高了约15%的速度！谢谢，如果没有其他更好的方案，你将获得奖励。 - Leo Adberg

坦白说，我希望能超过15％，但我想这很大程度上取决于屏幕。 - CygnusX1

是的，这些场景比较密集，大多数边界框有一点重叠。 - Leo Adberg

此外，15% 仍然值得注意！我之前已经进行了一堆优化，每个优化都小于 15%。 - Leo Adberg

由于某种原因，函数中的对偶初始化占用了大量时间。由于我不需要任何初始值，所以我只是将其切换为两个堆栈（一个用于节点，一个用于距离），这样速度更快！ - Leo Adberg