Rayon SSE rapide - intersection de 4 triangles

Question

Je travaille actuellement sur un Path Tracer, et je cherche des façons d'optimiser mes intersections rayon-triangle. J'utilise actuellement la mise en œuvre de l'algorithme suivant SSE4 Moller-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; 
} 

(Si quelqu'un voit un moyen d'optimiser me le faire savoir). Ensuite, j'ai lu qu'il est possible d'effectuer des tests d'intersection sur 4 triangles simultanément en utilisant SIMD intructions. Mais comment faire ça? Je ne vois pas comment il serait possible de l'implanter d'une manière qui soit plus efficace que mon mode séquentiel.

Here le petit code concerné de mon moteur de rendu.

Answer 1

Il est possible de faire jusqu'à 16 triangles avec AVX512, 8 avec AVX2 et 4 avec SSE. L'astuce cependant, est de s'assurer que les données sont au format SOA. L'autre astuce consiste à ne pas "retourner faux" à tout moment (juste filtrer les résultats à la fin). Donc, vous entrez triangle ressemblerait à quelque chose comme:

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

Et votre rayon ressemblerait à ceci:

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

Le code mathématique commence alors à regarder bien plus agréable (attention les _mm_dp_ps n'est pas la fonction la plus rapide jamais écrit - et sachez également que l'accès à l'implémentation interne des types __m128/__ m256/__ m512 n'est pas portable).

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

Vous avez essentiellement 4 conditions dans la méthode:

if (a > negativeEpsilon && a < positiveEpsilon) 

if (u < 0.0f) 

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

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

Si l'une de ces conditions sont vraies, alors il n'y a pas d'intersection. Cela exige essentiellement un peu refactoring (dans le code pseudo)

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

Donc finalement, si une intersection a eu lieu, le résultat sera t. Si une intersection N'A PAS se, alors nous devons mettre le résultat à quelque chose de mal (un nombre négatif est peut-être un bon choix dans ce cas!)

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

Alors que le code final peut sembler un peu méchant, il finira par être beaucoup plus rapide que votre approche. Le diable est dans les détails cependant ....

Un problème majeur avec cette approche est que vous avez le choix entre tester 1 rayon contre 8 triangles, ou tester 8 rayons contre 1 triangle. Pour les rayons principalement ce n'est probablement pas une grosse affaire. Pour les rayons secondaires qui ont l'habitude de se disperser dans des directions différentes), les choses peuvent commencer à devenir un peu gênantes. Il y a de bonnes chances que la plus grande partie du code de suivi de rayon se termine par un modèle de: test -> sort -> batch -> test -> sort -> lot

Si vous ne suivez pas ce modèle, vous êtes à peu près jamais tirer le meilleur parti des unités vectorielles. (Heureusement, les instructions de compression/expansion dans AVX512 aider beaucoup à ce sujet!)

Answer 2

J'ai fini avec le code de travail suivant

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

RÉSULTATS

Les résultats sont impressionnants. Pour une scène en triangle de 100k j'ai une accélération de 84% !! Pour une très petite scène (20 triangles) j'ai une perte de performance de 13%. Mais ça va parce que ceux-ci ne sont pas habituels.