快速SSE射线与4个三角形的相交

Question

我目前正在开发一款路径追踪器，并寻求优化射线三角形相交的方法。我目前使用的是Moller-Trumbore算法的sse4实现：

bool Ray::intersectTriangle(const Triangle tri, float& result) const
{
    __m128 q = _mm_cross_ps(m_directionM128, tri.e2);

    float a = _mm_dp_ps(tri.e1, q, dotProductMask).m128_f32[0];

    if (a > negativeEpsilon && a < positiveEpsilon)
        return false;

    float f = 1.0f / a;

    __m128 s = _mm_sub_ps(m_originM128, tri.v0);
    float u = f * _mm_dp_ps(s, q, dotProductMask).m128_f32[0];

    if (u < 0.0f)
        return false;

    __m128 r = _mm_cross_ps(s, tri.e1);
    float v = f * _mm_dp_ps(m_directionM128, r, dotProductMask).m128_f32[0];

    if (v < 0.0f || (u + v > 1.0f))
        return false;

    float t = f * _mm_dp_ps(tri.e2, r, dotProductMask).m128_f32[0];
    if (t < 0.0f || t > m_length)
        return false;

    result = t;

    return true;
}

（如果有人看到了优化的方法，请告诉我。）然后，我读到可以使用SIMD指令同时执行4个三角形的交集测试。但是如何做到这一点呢？我不知道如何实现它，以使其比我的顺序方式更加高效。

在这里，这是我渲染器相关的小代码。（点击链接）

- Cloyz

---

你看过Embree了吗？https://embree.github.io/ - Roman Kolesnikov

谢谢！我会看一下的。 - Cloyz

Answer 1

使用AVX512可以处理最多16个三角形，AVX2可以处理8个，SSE可以处理4个。但是，关键在于确保数据采用SOA格式。另一个技巧是不要在任何时候“返回false”（只需在最后过滤结果即可）。因此，您的三角形输入应该类似于：

  struct Tri {
    __m256 e1[3];
    __m256 e2[3];
    __m256 v0[3];
  };

你的光线会像这样：

  struct Ray {
    __m256 dir[3];
    __m256 pos[3];
  };

数学代码看起来更加美观（请注意，_mm_dp_ps不是有史以来最快的函数 - 还要注意访问__m128/__m256/__m512类型的内部实现不具可移植性）。

  #define or8f _mm256_or_ps
  #define mul _mm256_mul_ps
  #define fmsub _mm256_fmsub_ps
  #define fmadd _mm256_fmadd_ps

  void cross(__m256 result[3], const __m256 a[3], const __m256 b[3])
  {
    result[0] = fmsub(a[1], b[2], mul(b[1], a[2]));
    result[1] = fmsub(a[2], b[0], mul(b[2], a[0]));
    result[2] = fmsub(a[0], b[1], mul(b[0], a[1]));
  }
  __m256 dot(const __m256 a[3], const __m256 b[3])
  {
    return fmadd(a[2], b[2], fmadd(a[1], b[1], mul(a[0], b[0])));
  }

该方法基本上有四个条件：

if (a > negativeEpsilon && a < positiveEpsilon)

if (u < 0.0f)

if (v < 0.0f || (u + v > 1.0f))

if (t < 0.0f || t > m_length)

如果这些条件中有任何一个为真，则没有交集。这基本上需要进行一些重构（伪代码）。

__m256 condition0 = (a > negativeEpsilon && a < positiveEpsilon);

__m256 condition1 = (u < 0.0f)

__m256 condition2 = (v < 0.0f || (u + v > 1.0f))

__m256 condition3 = (t < 0.0f || t > m_length)

// combine all conditions that can cause failure.
__m256 failed = or8f(or8f(condition0, condition1), or8f(condition2, condition3));

最终，如果发生了交集，结果将是t。如果没有交集，则需要将结果设置为错误的值（在这种情况下，负数可能是一个不错的选择！）

// if(failed) return -1;
// else return t;
return _mm256_blendv_ps(t, _mm256_set1_ps(-1.0f), failed);

虽然最终代码可能看起来有点丑陋，但它最终会比您的方法快得多。但魔鬼在于细节......

这种方法的一个主要问题是，您可以选择测试1个光线与8个三角形，或测试8个光线与1个三角形。对于主要光线，这可能并不是什么大问题。对于次要光线（它们有在不同方向上散射的习惯），事情可能开始变得有点烦人。很有可能大部分光线跟踪代码最终都会遵循以下模式：测试->排序->批处理->测试->排序->批处理。

如果您不遵循该模式，那么您几乎永远无法充分利用向量单元。（幸运的是，AVX512中的压缩/展开指令在很大程度上有所帮助！）

- robthebloke

---

非常感谢！我同意对一个三角形发射 8 条光线会是一项改进，因为我还可以从主要光线的迷你锥体剔除中受益。对不同路径的光线结果进行批处理和重新组织似乎相当具有挑战性 :) - Cloyz

对于最后一行，我认为前两个操作数应该交换位置。因此：_mm256_blendv_ps(t，_mm256_set1_ps（-1.0f），failed)； - Cloyz

Answer 2

我最终得到了以下可用代码。

struct PackedTriangles
{
    __m256 e1[3];
    __m256 e2[3];
    __m256 v0[3];
    __m256 inactiveMask; // Required. We cant always have 8 triangles per packet.
};

struct PackedIntersectionResult
{
    float t = Math::infinity<float>();
    int idx;
};

struct PackedRay
{
    __m256 m_origin[3];
    __m256 m_direction[3];
    __m256 m_length;

    bool intersect(const PackedTriangles& packedTris, PackedIntersectionResult& result) const;
};

#define or8f _mm256_or_ps
#define mul _mm256_mul_ps
#define fmsub _mm256_fmsub_ps
#define fmadd _mm256_fmadd_ps
#define cmp _mm256_cmp_ps
#define div _mm256_div_ps

void avx_multi_cross(__m256 result[3], const __m256 a[3], const __m256 b[3])
{
    result[0] = fmsub(a[1], b[2], mul(b[1], a[2]));
    result[1] = fmsub(a[2], b[0], mul(b[2], a[0]));
    result[2] = fmsub(a[0], b[1], mul(b[0], a[1]));
}

__m256 avx_multi_dot(const __m256 a[3], const __m256 b[3])
{
    return fmadd(a[2], b[2], fmadd(a[1], b[1], mul(a[0], b[0])));
}

void avx_multi_sub(__m256 result[3], const __m256 a[3], const __m256 b[3])
{
    result[0] = _mm256_sub_ps(a[0], b[0]);
    result[1] = _mm256_sub_ps(a[1], b[1]);
    result[2] = _mm256_sub_ps(a[2], b[2]);
}

const __m256 oneM256 = _mm256_set1_ps(1.0f);
const __m256 minusOneM256 = _mm256_set1_ps(-1.0f);
const __m256 positiveEpsilonM256 = _mm256_set1_ps(1e-6f);
const __m256 negativeEpsilonM256 = _mm256_set1_ps(-1e-6f);
const __m256 zeroM256 = _mm256_set1_ps(0.0f);

bool PackedRay::intersect(const PackedTriangles& packedTris, PackedIntersectionResult& result) const
{
    __m256 q[3];
    avx_multi_cross(q, m_direction, packedTris.e2);

    __m256 a = avx_multi_dot(packedTris.e1, q);

    __m256 f = div(oneM256, a);

    __m256 s[3];
    avx_multi_sub(s, m_origin, packedTris.v0);

    __m256 u = mul(f, avx_multi_dot(s, q));

    __m256 r[3];
    avx_multi_cross(r, s, packedTris.e1);

    __m256 v = mul(f, avx_multi_dot(m_direction, r));

    __m256 t = mul(f, avx_multi_dot(packedTris.e2, r));

    // Failure conditions
    __m256 failed = _mm256_and_ps(
        cmp(a, negativeEpsilonM256, _CMP_GT_OQ),
        cmp(a, positiveEpsilonM256, _CMP_LT_OQ)
    );

    failed = or8f(failed, cmp(u, zeroM256, _CMP_LT_OQ));
    failed = or8f(failed, cmp(v, zeroM256, _CMP_LT_OQ));
    failed = or8f(failed, cmp(_mm256_add_ps(u, v), oneM256, _CMP_GT_OQ));
    failed = or8f(failed, cmp(t, zeroM256, _CMP_LT_OQ));
    failed = or8f(failed, cmp(t, m_length, _CMP_GT_OQ));
    failed = or8f(failed, packedTris.inactiveMask);

    __m256 tResults = _mm256_blendv_ps(t, minusOneM256, failed);

    int mask = _mm256_movemask_ps(tResults);
    if (mask != 0xFF)
    {
        // There is at least one intersection
        result.idx = -1;

        float* ptr = (float*)&tResults;
        for (int i = 0; i < 8; ++i)
        {
            if (ptr[i] >= 0.0f && ptr[i] < result.t)
            {
                result.t = ptr[i];
                result.idx = i;
            }
        }

        return result.idx != -1;
    }

    return false;
}

结果

这些结果真的很棒。对于一个包含100,000个三角形的场景，我获得了84%的加速!! 对于非常小的场景（只有20个三角形），性能下降了13%。但这没关系，因为这种情况并不常见。

- Cloyz