优化光线追踪器中的 BVH 遍历

Question

在我的 CPU 光线追踪器（嗯，路径追踪器）中，大部分 CPU 时间花费在 BVH 遍历函数中。根据我的分析器，75% 的光线追踪时间花在这个函数和它调用的函数上，而 35% 的时间花在函数本身上。另外 40% 在它调用的不同交叉点测试中。

基本上，代码对所有与之相交的边界框和三角形进行 DFS 遍历。它使用堆栈上的静态分配数组来保存要探索的节点（BVHSTACKSIZE 设置为 32，它从不需要大部分空间），因此不会动态分配内存。然而，对我来说，35% 的时间都花在这里似乎很疯狂。我花了一段时间优化代码，目前它是我能做到的最快速度，但它仍然是我程序中最大的单一瓶颈。

有没有人有更多优化这个的提示？我已经有了一个不错的 BVH 构造算法，所以我认为使用不同的 BVH 不会得到任何加速。有没有人有关于如何在 Mac 上最好地进行逐行分析的提示？

作为参考，示例场景中的此代码所需时间从

谢谢！

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);
  }
}

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; }
};

Answer 2

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

我相信你感兴趣的是你的射线上的第一击。您不会从沿射线的多次命中中进行任何累积，对吗？ （好像原语半透明或什么的）。是吗？

所以 - 如果上述情况属实，我发现您的遍历顺序存在问题。如果我正确阅读了您的代码，那么您将无条件地先遍历cur->l，然后再遍历cur->r。如果你的光线首先击中左边的盒子，这很好。但是 - 它可能来自场景的另一侧 - 然后你执行的交叉点比需要的多，基本上是沿着你的光线来回遍历。

相反，你想沿着你的射线从前到后遍历，希望你早点击中一些东西，这样你就可以跳过大部分进一步的遍历测试。

幸运的是，到边界框的距离已经在交集函数中计算出来，您可以轻松检查哪个首先相交。它只需要返回一个到交叉点的float 距离而不是bool（例如+infinity 失败）。您需要在将内容推送到您的 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

一旦您实施了上述操作，您可能会注意到还有一项改进需要改进。如果你在 push 和 pop 之间点击了一些原语并且你的 ray.max_t 被减少了，那么“我们已经知道我们想要进入这个节点”的条件不一定是真的。考虑到这一点，您可能需要在堆栈中存储一对 (t_enter, node)。然后，当你弹出时，你用你当前的ray.max_t 仔细检查t_enter。这可以为您节省最多 h 次遍历检查，其中 h 是树的高度。

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

Answer 3

我看到您可以进行三项改进。

第一个大问题（很难）是您的代码中有许多条件分支，这肯定会减慢您的 CPU，因为它无法很好地预测代码路径（编译时也是如此）。例如，我看到你先相交，然后测试节点是否是叶子，然后与所有 prims 相交。你能先测试它是否是一片叶子，然后做正确的交叉点吗？这会稍微减少分支。

其次，你的 BVH 的内存布局是什么？您能否对其进行优化以使其对您的遍历友好。您可以尝试查看遍历期间发生的缓存未命中次数，这将很好地表明您的内存是否具有正确的布局。虽然没有直接关系，但现在你的平台和底层硬件很好。我推荐阅读this。

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

另外，我刚刚看到您使用的是double。这有什么原因吗？我们的 pathtracer 根本不使用 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;
}