diff options
author | Sergey Sharybin <sergey.vfx@gmail.com> | 2017-03-27 18:06:37 +0300 |
---|---|---|
committer | Sergey Sharybin <sergey.vfx@gmail.com> | 2017-03-28 18:26:47 +0300 |
commit | 6ea54fe9ffe2b2514990fdf3489ca53d05ce449a (patch) | |
tree | 3026ce049c43bf8e0d56fbad57e3cbb844444fc4 /intern/cycles/util | |
parent | 69aa6577b3dfea5d8a6d915fad7fb7650d8b6278 (diff) |
Cycles: Switch to reformulated Pluecker ray/triangle intersection
The intention of this commit it to address issues mentioned in the
reports T43865,T50164 and T50452.
The code is based on Embree code with some extra vectorization
to speed up single ray to single triangle intersection.
Unfortunately, such a fix is not coming for free. There is some
slowdown for AVX2 processors, mainly due to different vectorization
code, which caused different number of instructions to be executed
and different instructions-per-cycle counters. But on another hand
this commit makes pre-AVX2 platforms such as AVX and SSE4.1 a bit
faster. The prerformance goes as following:
2.78c AVX2 2.78c AVX Patch AVX2 Patch AVX
BMW 05:21.09 06:05.34 05:32.97 (+3.5%) 05:34.97 (-8.5%)
Classroom 16:55.36 18:24.51 17:10.41 (+1.4%) 17:15.87 (-6.3%)
Fishy Cat 08:08.49 08:36.26 08:09.19 (+0.2%) 08:12.25 (-4.7%
Koro 11:22.54 11:45.24 11:13.25 (-1.5%) 11:43.81 (-0.3%)
Barcelone 14:18.32 16:09.46 14:15.20 (-0.4%) 14:25.15 (-10.8%)
On GPU the performance is about 1.5-2% slower in my tests on GTX1080
but afraid we can't do much as a part of this chaneg here and
consider it a price to pay for more proper intersection check.
Made in collaboration with Maxym Dmytrychenko, big thanks to him!
Reviewers: brecht, juicyfruit, lukasstockner97, dingto
Differential Revision: https://developer.blender.org/D1574
Diffstat (limited to 'intern/cycles/util')
-rw-r--r-- | intern/cycles/util/util_math_intersect.h | 264 |
1 files changed, 77 insertions, 187 deletions
diff --git a/intern/cycles/util/util_math_intersect.h b/intern/cycles/util/util_math_intersect.h index 9e0587e1afb..f0f505c311e 100644 --- a/intern/cycles/util/util_math_intersect.h +++ b/intern/cycles/util/util_math_intersect.h @@ -79,216 +79,106 @@ ccl_device bool ray_aligned_disk_intersect( return true; } -/* Optimized watertight ray-triangle intersection. - * - * Sven Woop - * Watertight Ray/Triangle Intersection - * - * http://jcgt.org/published/0002/01/05/paper.pdf - */ - -/* Precalculated data for the ray->tri intersection. */ -typedef struct TriangleIsectPrecalc { - /* Maximal dimension kz, and orthogonal dimensions. */ - int kx, ky, kz; - - /* Shear constants. */ - float Sx, Sy, Sz; -} TriangleIsectPrecalc; - -/* Workaround stupidness of CUDA/OpenCL which doesn't allow to access indexed - * component of float3 value. - */ -#ifdef __KERNEL_GPU__ -# define IDX(vec, idx) \ - ((idx == 0) ? ((vec).x) : ( (idx == 1) ? ((vec).y) : ((vec).z) )) -#else -# define IDX(vec, idx) ((vec)[idx]) -#endif - -#if (defined(__KERNEL_OPENCL_APPLE__)) || \ - (defined(__KERNEL_CUDA__) && (defined(i386) || defined(_M_IX86))) -ccl_device_noinline -#else -ccl_device_inline -#endif -void ray_triangle_intersect_precalc(float3 dir, - TriangleIsectPrecalc *isect_precalc) -{ - /* Calculate dimension where the ray direction is maximal. */ -#ifndef __KERNEL_SSE__ - int kz = util_max_axis(make_float3(fabsf(dir.x), - fabsf(dir.y), - fabsf(dir.z))); - int kx = kz + 1; if(kx == 3) kx = 0; - int ky = kx + 1; if(ky == 3) ky = 0; -#else - int kx, ky, kz; - /* Avoiding mispredicted branch on direction. */ - kz = util_max_axis(fabs(dir)); - static const char inc_xaxis[] = {1, 2, 0, 55}; - static const char inc_yaxis[] = {2, 0, 1, 55}; - kx = inc_xaxis[kz]; - ky = inc_yaxis[kz]; -#endif - - float dir_kz = IDX(dir, kz); - - /* Swap kx and ky dimensions to preserve winding direction of triangles. */ - if(dir_kz < 0.0f) { - int tmp = kx; - kx = ky; - ky = tmp; - } - - /* Calculate the shear constants. */ - float inv_dir_z = 1.0f / dir_kz; - isect_precalc->Sx = IDX(dir, kx) * inv_dir_z; - isect_precalc->Sy = IDX(dir, ky) * inv_dir_z; - isect_precalc->Sz = inv_dir_z; - - /* Store the dimensions. */ - isect_precalc->kx = kx; - isect_precalc->ky = ky; - isect_precalc->kz = kz; -} - ccl_device_forceinline bool ray_triangle_intersect( - const TriangleIsectPrecalc *isect_precalc, - float3 ray_P, float ray_t, -#if defined(__KERNEL_AVX2__) && defined(__KERNEL_SSE__) + float3 ray_P, float3 ray_dir, float ray_t, +#if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__) const ssef *ssef_verts, #else const float3 tri_a, const float3 tri_b, const float3 tri_c, #endif float *isect_u, float *isect_v, float *isect_t) { - const int kx = isect_precalc->kx; - const int ky = isect_precalc->ky; - const int kz = isect_precalc->kz; - const float Sx = isect_precalc->Sx; - const float Sy = isect_precalc->Sy; - const float Sz = isect_precalc->Sz; - -#if defined(__KERNEL_AVX2__) && defined(__KERNEL_SSE__) - const avxf avxf_P(ray_P.m128, ray_P.m128); - const avxf tri_ab(_mm256_loadu_ps((float *)(ssef_verts))); - const avxf tri_bc(_mm256_loadu_ps((float *)(ssef_verts + 1))); - - const avxf AB = tri_ab - avxf_P; - const avxf BC = tri_bc - avxf_P; - - const __m256i permute_mask = _mm256_set_epi32(0x3, kz, ky, kx, 0x3, kz, ky, kx); - - const avxf AB_k = shuffle(AB, permute_mask); - const avxf BC_k = shuffle(BC, permute_mask); - - /* Akz, Akz, Bkz, Bkz, Bkz, Bkz, Ckz, Ckz */ - const avxf ABBC_kz = shuffle<2>(AB_k, BC_k); - - /* Akx, Aky, Bkx, Bky, Bkx,Bky, Ckx, Cky */ - const avxf ABBC_kxy = shuffle<0,1,0,1>(AB_k, BC_k); - - const avxf Sxy(Sy, Sx, Sy, Sx); - - /* Ax, Ay, Bx, By, Bx, By, Cx, Cy */ - const avxf ABBC_xy = nmadd(ABBC_kz, Sxy, ABBC_kxy); - - float ABBC_kz_array[8]; - _mm256_storeu_ps((float*)&ABBC_kz_array, ABBC_kz); - - const float A_kz = ABBC_kz_array[0]; - const float B_kz = ABBC_kz_array[2]; - const float C_kz = ABBC_kz_array[6]; - - /* By, Bx, Cy, Cx, By, Bx, Ay, Ax */ - const avxf BCBA_yx = permute<3,2,7,6,3,2,1,0>(ABBC_xy); - - const avxf neg_mask(0,0,0,0,0x80000000, 0x80000000, 0x80000000, 0x80000000); - - /* W U V - * (AxBy-AyBx) (BxCy-ByCx) XX XX (BxBy-ByBx) (CxAy-CyAx) XX XX - */ - const avxf WUxxxxVxx_neg = _mm256_hsub_ps(ABBC_xy * BCBA_yx, neg_mask /* Dont care */); - - const avxf WUVWnegWUVW = permute<0,1,5,0,0,1,5,0>(WUxxxxVxx_neg) ^ neg_mask; - - /* Calculate scaled barycentric coordinates. */ - float WUVW_array[4]; - _mm_storeu_ps((float*)&WUVW_array, _mm256_castps256_ps128 (WUVWnegWUVW)); - - const float W = WUVW_array[0]; - const float U = WUVW_array[1]; - const float V = WUVW_array[2]; - - const int WUVW_mask = 0x7 & _mm256_movemask_ps(WUVWnegWUVW); - const int WUVW_zero = 0x7 & _mm256_movemask_ps(_mm256_cmp_ps(WUVWnegWUVW, - _mm256_setzero_ps(), 0)); - - if(!((WUVW_mask == 7) || (WUVW_mask == 0)) && ((WUVW_mask | WUVW_zero) != 7)) { - return false; - } +#if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__) + typedef ssef float3; + const float3 tri_a(ssef_verts[0]); + const float3 tri_b(ssef_verts[1]); + const float3 tri_c(ssef_verts[2]); + const float3 P(ray_P); + const float3 dir(ray_dir); #else - /* Calculate vertices relative to ray origin. */ - const float3 A = make_float3(tri_a.x - ray_P.x, tri_a.y - ray_P.y, tri_a.z - ray_P.z); - const float3 B = make_float3(tri_b.x - ray_P.x, tri_b.y - ray_P.y, tri_b.z - ray_P.z); - const float3 C = make_float3(tri_c.x - ray_P.x, tri_c.y - ray_P.y, tri_c.z - ray_P.z); - - const float A_kx = IDX(A, kx), A_ky = IDX(A, ky), A_kz = IDX(A, kz); - const float B_kx = IDX(B, kx), B_ky = IDX(B, ky), B_kz = IDX(B, kz); - const float C_kx = IDX(C, kx), C_ky = IDX(C, ky), C_kz = IDX(C, kz); - - /* Perform shear and scale of vertices. */ - const float Ax = A_kx - Sx * A_kz; - const float Ay = A_ky - Sy * A_kz; - const float Bx = B_kx - Sx * B_kz; - const float By = B_ky - Sy * B_kz; - const float Cx = C_kx - Sx * C_kz; - const float Cy = C_ky - Sy * C_kz; +# define dot3(a, b) dot(a, b) + const float3 P = ray_P; + const float3 dir = ray_dir; +#endif - /* Calculate scaled barycentric coordinates. */ - float U = Cx * By - Cy * Bx; - float V = Ax * Cy - Ay * Cx; - float W = Bx * Ay - By * Ax; - if((U < 0.0f || V < 0.0f || W < 0.0f) && - (U > 0.0f || V > 0.0f || W > 0.0f)) - { + /* Calculate vertices relative to ray origin. */ + const float3 v0 = tri_c - P; + const float3 v1 = tri_a - P; + const float3 v2 = tri_b - P; + + /* Calculate triangle edges. */ + const float3 e0 = v2 - v0; + const float3 e1 = v0 - v1; + const float3 e2 = v1 - v2; + + /* Perform edge tests. */ +#ifdef __KERNEL_SSE2__ + const float3 crossU = cross(v2 + v0, e0); + const float3 crossV = cross(v0 + v1, e1); + const float3 crossW = cross(v1 + v2, e2); +# ifndef __KERNEL_SSE__ + const ssef crossX(crossU.x, crossV.x, crossW.x, crossW.x); + const ssef crossY(crossU.y, crossV.y, crossW.y, crossW.y); + const ssef crossZ(crossU.z, crossV.z, crossW.z, crossW.z); +# else + ssef crossX(crossU); + ssef crossY(crossV); + ssef crossZ(crossW); + ssef zero = _mm_setzero_ps(); + _MM_TRANSPOSE4_PS(crossX, crossY, crossZ, zero); +# endif + const ssef dirX(ray_dir.x); + const ssef dirY(ray_dir.y); + const ssef dirZ(ray_dir.z); + /*const*/ ssef UVWW = crossX*dirX + crossY*dirY + crossZ*dirZ; + const float minUVW = reduce_min(UVWW); + const float maxUVW = reduce_max(UVWW); +#else /* __KERNEL_SSE2__ */ + const float U = dot(cross(v2 + v0, e0), ray_dir); + const float V = dot(cross(v0 + v1, e1), ray_dir); + const float W = dot(cross(v1 + v2, e2), ray_dir); + const float minUVW = min(U, min(V, W)); + const float maxUVW = max(U, max(V, W)); +#endif /* __KERNEL_SSE2__ */ + + if(minUVW < 0.0f && maxUVW > 0.0f) { return false; } -#endif - /* Calculate determinant. */ - float det = U + V + W; - if(UNLIKELY(det == 0.0f)) { + /* Calculate geometry normal and denominator. */ + const float3 Ng1 = cross(e1, e0); + //const Vec3vfM Ng1 = stable_triangle_normal(e2,e1,e0); + const float3 Ng = Ng1 + Ng1; + const float den = dot3(Ng, dir); + /* Avoid division by 0. */ + if(UNLIKELY(den == 0.0f)) { return false; } - /* Calculate scaled z-coordinates of vertices and use them to calculate - * the hit distance. - */ - const float T = (U * A_kz + V * B_kz + W * C_kz) * Sz; - const int sign_det = (__float_as_int(det) & 0x80000000); - const float sign_T = xor_signmask(T, sign_det); + /* Perform depth test. */ + const float T = dot3(v0, Ng); + const int sign_den = (__float_as_int(den) & 0x80000000); + const float sign_T = xor_signmask(T, sign_den); if((sign_T < 0.0f) || - (sign_T > ray_t * xor_signmask(det, sign_det))) + (sign_T > ray_t * xor_signmask(den, sign_den))) { return false; } - /* Workaround precision error on CUDA. */ -#ifdef __KERNEL_CUDA__ - if(A == B && B == C) { - return false; - } + const float inv_den = 1.0f / den; +#ifdef __KERNEL_SSE2__ + UVWW *= inv_den; + _mm_store_ss(isect_u, UVWW); + _mm_store_ss(isect_v, shuffle<1,1,3,3>(UVWW)); +#else + *isect_u = U * inv_den; + *isect_v = V * inv_den; #endif - const float inv_det = 1.0f / det; - *isect_u = U * inv_det; - *isect_v = V * inv_det; - *isect_t = T * inv_det; + *isect_t = T * inv_den; return true; -} -#undef IDX +#undef dot3 +} ccl_device bool ray_quad_intersect(float3 ray_P, float3 ray_D, float ray_mint, float ray_maxt, |