快速SSE射线 - 四个三角形交叉点

Question

我目前正在研究Path Tracer，并且正在寻找优化我的射线三角形交点的方法。我目前使用下面的SSE4实施穆勒 - Trumbore算法：

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个三角形上执行相交测试。但如何做到这一点？我不认为如何以比我的顺序方式更有效的方式来实现它。

Here我的渲​​染器有关的小代码。

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

你基本上有4个条件中的方法：

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中的压缩/扩展指令帮助了很多！）

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

成绩

结果是真棒。对于一个10万的三角场景我有84％的加速 !!。对于一个非常小的场景（20个三角形），我有13％的表现损失。但没关系，因为这些并不常见。