Cycles: Implement AVX2 version of triangle_intersect
[blender-staging.git] / intern / cycles / kernel / geom / geom_triangle_intersect.h
1 /*
2  * Copyright 2014, Blender Foundation.
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16
17 /* Triangle/Ray intersections.
18  *
19  * For BVH ray intersection we use a precomputed triangle storage to accelerate
20  * intersection at the cost of more memory usage.
21  */
22
23 CCL_NAMESPACE_BEGIN
24
25 /* Workaround stupidness of CUDA/OpenCL which doesn't allow to access indexed
26  * component of float3 value.
27  */
28 #ifndef __KERNEL_CPU__
29 #  define IDX(vec, idx) \
30     ((idx == 0) ? ((vec).x) : ( (idx == 1) ? ((vec).y) : ((vec).z) ))
31 #else
32 #  define IDX(vec, idx) ((vec)[idx])
33 #endif
34
35 /* Ray-Triangle intersection for BVH traversal
36  *
37  * Sven Woop
38  * Watertight Ray/Triangle Intersection
39  *
40  * http://jcgt.org/published/0002/01/05/paper.pdf
41  */
42
43 /* Precalculated data for the ray->tri intersection. */
44 typedef struct IsectPrecalc {
45         /* Maximal dimension kz, and orthogonal dimensions. */
46         int kx, ky, kz;
47
48         /* Shear constants. */
49         float Sx, Sy, Sz;
50 } IsectPrecalc;
51
52 #if (defined(__KERNEL_OPENCL_APPLE__)) || \
53     (defined(__KERNEL_CUDA__) && (defined(i386) || defined(_M_IX86)))
54 ccl_device_noinline
55 #else
56 ccl_device_inline
57 #endif
58 void triangle_intersect_precalc(float3 dir,
59                                 IsectPrecalc *isect_precalc)
60 {
61         /* Calculate dimension where the ray direction is maximal. */
62         int kz = util_max_axis(make_float3(fabsf(dir.x),
63                                            fabsf(dir.y),
64                                            fabsf(dir.z)));
65         int kx = kz + 1; if(kx == 3) kx = 0;
66         int ky = kx + 1; if(ky == 3) ky = 0;
67
68         /* Swap kx and ky dimensions to preserve winding direction of triangles. */
69         if(IDX(dir, kz) < 0.0f) {
70                 int tmp = kx;
71                 kx = ky;
72                 ky = tmp;
73         }
74
75         /* Calculate the shear constants. */
76         float inv_dir_z = 1.0f / IDX(dir, kz);
77         isect_precalc->Sx = IDX(dir, kx) * inv_dir_z;
78         isect_precalc->Sy = IDX(dir, ky) * inv_dir_z;
79         isect_precalc->Sz = inv_dir_z;
80
81         /* Store the dimensions. */
82         isect_precalc->kx = kx;
83         isect_precalc->ky = ky;
84         isect_precalc->kz = kz;
85 }
86
87 /* TODO(sergey): Make it general utility function. */
88 ccl_device_inline float xor_signmask(float x, int y)
89 {
90         return __int_as_float(__float_as_int(x) ^ y);
91 }
92
93 ccl_device_inline bool triangle_intersect(KernelGlobals *kg,
94                                           const IsectPrecalc *isect_precalc,
95                                           Intersection *isect,
96                                           float3 P,
97                                           uint visibility,
98                                           int object,
99                                           int triAddr)
100 {
101         const int kx = isect_precalc->kx;
102         const int ky = isect_precalc->ky;
103         const int kz = isect_precalc->kz;
104         const float Sx = isect_precalc->Sx;
105         const float Sy = isect_precalc->Sy;
106         const float Sz = isect_precalc->Sz;
107
108         /* Calculate vertices relative to ray origin. */
109         const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, triAddr);
110
111 #if defined(__KERNEL_AVX2__)
112         const avxf avxf_P(P.m128, P.m128);
113
114         const avxf tri_ab = kernel_tex_fetch_avxf(__prim_tri_verts, tri_vindex + 0);
115         const avxf tri_bc = kernel_tex_fetch_avxf(__prim_tri_verts, tri_vindex + 1);
116
117         const avxf AB = tri_ab - avxf_P;
118         const avxf BC = tri_bc - avxf_P;
119
120         const __m256i permuteMask = _mm256_set_epi32(0x3, kz, ky, kx, 0x3, kz, ky, kx);
121
122         const avxf AB_k = shuffle(AB, permuteMask);
123         const avxf BC_k = shuffle(BC, permuteMask);
124
125         /* Akz, Akz, Bkz, Bkz, Bkz, Bkz, Ckz, Ckz */
126         const avxf ABBC_kz = shuffle<2>(AB_k, BC_k);
127
128         /* Akx, Aky, Bkx, Bky, Bkx,Bky, Ckx, Cky */
129         const avxf ABBC_kxy = shuffle<0,1,0,1>(AB_k, BC_k);
130
131         const avxf Sxy(Sy, Sx, Sy, Sx);
132
133         /* Ax, Ay, Bx, By, Bx, By, Cx, Cy */
134         const avxf ABBC_xy = nmadd(ABBC_kz, Sxy, ABBC_kxy);
135
136         float ABBC_kz_array[8];
137         _mm256_storeu_ps((float*)&ABBC_kz_array, ABBC_kz);
138
139         const float A_kz = ABBC_kz_array[0];
140         const float B_kz = ABBC_kz_array[2];
141         const float C_kz = ABBC_kz_array[6];
142
143         /* By, Bx, Cy, Cx, By, Bx, Ay, Ax */
144         const avxf BCBA_yx = permute<3,2,7,6,3,2,1,0>(ABBC_xy);
145
146         const avxf negMask(0,0,0,0,0x80000000, 0x80000000, 0x80000000, 0x80000000);
147
148         /* W           U                             V
149          * (AxBy-AyBx) (BxCy-ByCx) XX XX (BxBy-ByBx) (CxAy-CyAx) XX XX
150          */
151         const avxf WUxxxxVxx_neg = _mm256_hsub_ps(ABBC_xy * BCBA_yx, negMask /* Dont care */);
152
153         const avxf WUVWnegWUVW = permute<0,1,5,0,0,1,5,0>(WUxxxxVxx_neg) ^ negMask;
154
155         /* Calculate scaled barycentric coordinates. */
156         float WUVW_array[4];
157         _mm_storeu_ps((float*)&WUVW_array, _mm256_castps256_ps128 (WUVWnegWUVW));
158
159         const float W = WUVW_array[0];
160         const float U = WUVW_array[1];
161         const float V = WUVW_array[2];
162
163         const int WUVW_mask = 0x7 & _mm256_movemask_ps(WUVWnegWUVW);
164         const int WUVW_zero = 0x7 & _mm256_movemask_ps(_mm256_cmp_ps(WUVWnegWUVW,
165                                                        _mm256_setzero_ps(), 0));
166
167         if(!((WUVW_mask == 7) || (WUVW_mask == 0)) && ((WUVW_mask | WUVW_zero) != 7)) {
168                 return false;
169         }
170 #else
171         const float4 tri_a = kernel_tex_fetch(__prim_tri_verts, tri_vindex+0),
172                      tri_b = kernel_tex_fetch(__prim_tri_verts, tri_vindex+1),
173                      tri_c = kernel_tex_fetch(__prim_tri_verts, tri_vindex+2);
174         const float3 A = make_float3(tri_a.x - P.x, tri_a.y - P.y, tri_a.z - P.z);
175         const float3 B = make_float3(tri_b.x - P.x, tri_b.y - P.y, tri_b.z - P.z);
176         const float3 C = make_float3(tri_c.x - P.x, tri_c.y - P.y, tri_c.z - P.z);
177
178         const float A_kx = IDX(A, kx), A_ky = IDX(A, ky), A_kz = IDX(A, kz);
179         const float B_kx = IDX(B, kx), B_ky = IDX(B, ky), B_kz = IDX(B, kz);
180         const float C_kx = IDX(C, kx), C_ky = IDX(C, ky), C_kz = IDX(C, kz);
181
182         /* Perform shear and scale of vertices. */
183         const float Ax = A_kx - Sx * A_kz;
184         const float Ay = A_ky - Sy * A_kz;
185         const float Bx = B_kx - Sx * B_kz;
186         const float By = B_ky - Sy * B_kz;
187         const float Cx = C_kx - Sx * C_kz;
188         const float Cy = C_ky - Sy * C_kz;
189
190         /* Calculate scaled barycentric coordinates. */
191         float U = Cx * By - Cy * Bx;
192         float V = Ax * Cy - Ay * Cx;
193         float W = Bx * Ay - By * Ax;
194         if((U < 0.0f || V < 0.0f || W < 0.0f) &&
195            (U > 0.0f || V > 0.0f || W > 0.0f))
196         {
197                 return false;
198         }
199 #endif
200
201         /* Calculate determinant. */
202         float det = U + V + W;
203         if(UNLIKELY(det == 0.0f)) {
204                 return false;
205         }
206
207         /* Calculate scaled z-coordinates of vertices and use them to calculate
208          * the hit distance.
209          */
210         const float T = (U * A_kz + V * B_kz + W * C_kz) * Sz;
211         const int sign_det = (__float_as_int(det) & 0x80000000);
212         const float sign_T = xor_signmask(T, sign_det);
213         if((sign_T < 0.0f) ||
214            (sign_T > isect->t * xor_signmask(det, sign_det)))
215         {
216                 return false;
217         }
218
219 #ifdef __VISIBILITY_FLAG__
220         /* visibility flag test. we do it here under the assumption
221          * that most triangles are culled by node flags */
222         if(kernel_tex_fetch(__prim_visibility, triAddr) & visibility)
223 #endif
224         {
225 #ifdef __KERNEL_CUDA__
226                 if(A == B && B == C) {
227                         return false;
228                 }
229 #endif
230                 /* Normalize U, V, W, and T. */
231                 const float inv_det = 1.0f / det;
232                 isect->prim = triAddr;
233                 isect->object = object;
234                 isect->type = PRIMITIVE_TRIANGLE;
235                 isect->u = U * inv_det;
236                 isect->v = V * inv_det;
237                 isect->t = T * inv_det;
238                 return true;
239         }
240         return false;
241 }
242
243 /* Special ray intersection routines for subsurface scattering. In that case we
244  * only want to intersect with primitives in the same object, and if case of
245  * multiple hits we pick a single random primitive as the intersection point.
246  */
247
248 #ifdef __SUBSURFACE__
249 ccl_device_inline void triangle_intersect_subsurface(
250         KernelGlobals *kg,
251         const IsectPrecalc *isect_precalc,
252         SubsurfaceIntersection *ss_isect,
253         float3 P,
254         int object,
255         int triAddr,
256         float tmax,
257         uint *lcg_state,
258         int max_hits)
259 {
260         const int kx = isect_precalc->kx;
261         const int ky = isect_precalc->ky;
262         const int kz = isect_precalc->kz;
263         const float Sx = isect_precalc->Sx;
264         const float Sy = isect_precalc->Sy;
265         const float Sz = isect_precalc->Sz;
266
267         /* Calculate vertices relative to ray origin. */
268         const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, triAddr);
269         const float4 tri_a = kernel_tex_fetch(__prim_tri_verts, tri_vindex+0),
270                      tri_b = kernel_tex_fetch(__prim_tri_verts, tri_vindex+1),
271                      tri_c = kernel_tex_fetch(__prim_tri_verts, tri_vindex+2);
272         const float3 A = make_float3(tri_a.x - P.x, tri_a.y - P.y, tri_a.z - P.z);
273         const float3 B = make_float3(tri_b.x - P.x, tri_b.y - P.y, tri_b.z - P.z);
274         const float3 C = make_float3(tri_c.x - P.x, tri_c.y - P.y, tri_c.z - P.z);
275
276         const float A_kx = IDX(A, kx), A_ky = IDX(A, ky), A_kz = IDX(A, kz);
277         const float B_kx = IDX(B, kx), B_ky = IDX(B, ky), B_kz = IDX(B, kz);
278         const float C_kx = IDX(C, kx), C_ky = IDX(C, ky), C_kz = IDX(C, kz);
279
280         /* Perform shear and scale of vertices. */
281         const float Ax = A_kx - Sx * A_kz;
282         const float Ay = A_ky - Sy * A_kz;
283         const float Bx = B_kx - Sx * B_kz;
284         const float By = B_ky - Sy * B_kz;
285         const float Cx = C_kx - Sx * C_kz;
286         const float Cy = C_ky - Sy * C_kz;
287
288         /* Calculate scaled barycentric coordinates. */
289         float U = Cx * By - Cy * Bx;
290         float V = Ax * Cy - Ay * Cx;
291         float W = Bx * Ay - By * Ax;
292
293         if((U < 0.0f || V < 0.0f || W < 0.0f) &&
294            (U > 0.0f || V > 0.0f || W > 0.0f))
295         {
296                 return;
297         }
298
299         /* Calculate determinant. */
300         float det = U + V + W;
301         if(UNLIKELY(det == 0.0f)) {
302                 return;
303         }
304
305         /* Calculate scaled z−coordinates of vertices and use them to calculate
306          * the hit distance.
307          */
308         const int sign_det = (__float_as_int(det) & 0x80000000);
309         const float T = (U * A_kz + V * B_kz + W * C_kz) * Sz;
310         const float sign_T = xor_signmask(T, sign_det);
311         if((sign_T < 0.0f) ||
312            (sign_T > tmax * xor_signmask(det, sign_det)))
313         {
314                 return;
315         }
316
317         /* Normalize U, V, W, and T. */
318         const float inv_det = 1.0f / det;
319
320         const float t = T * inv_det;
321         for(int i = min(max_hits, ss_isect->num_hits) - 1; i >= 0; --i) {
322                 if(ss_isect->hits[i].t == t) {
323                         return;
324                 }
325         }
326
327         ss_isect->num_hits++;
328         int hit;
329
330         if(ss_isect->num_hits <= max_hits) {
331                 hit = ss_isect->num_hits - 1;
332         }
333         else {
334                 /* reservoir sampling: if we are at the maximum number of
335                  * hits, randomly replace element or skip it */
336                 hit = lcg_step_uint(lcg_state) % ss_isect->num_hits;
337
338                 if(hit >= max_hits)
339                         return;
340         }
341
342         /* record intersection */
343         Intersection *isect = &ss_isect->hits[hit];
344         isect->prim = triAddr;
345         isect->object = object;
346         isect->type = PRIMITIVE_TRIANGLE;
347         isect->u = U * inv_det;
348         isect->v = V * inv_det;
349         isect->t = t;
350
351         /* Record geometric normal. */
352         /* TODO(sergey): Use float4_to_float3() on just an edges. */
353         const float3 v0 = float4_to_float3(tri_a);
354         const float3 v1 = float4_to_float3(tri_b);
355         const float3 v2 = float4_to_float3(tri_c);
356         ss_isect->Ng[hit] = normalize(cross(v1 - v0, v2 - v0));
357 }
358 #endif
359
360 /* Refine triangle intersection to more precise hit point. For rays that travel
361  * far the precision is often not so good, this reintersects the primitive from
362  * a closer distance. */
363
364 /* Reintersections uses the paper:
365  *
366  * Tomas Moeller
367  * Fast, minimum storage ray/triangle intersection
368  * http://www.cs.virginia.edu/~gfx/Courses/2003/ImageSynthesis/papers/Acceleration/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf
369  */
370
371 ccl_device_inline float3 triangle_refine(KernelGlobals *kg,
372                                          ShaderData *sd,
373                                          const Intersection *isect,
374                                          const Ray *ray)
375 {
376         float3 P = ray->P;
377         float3 D = ray->D;
378         float t = isect->t;
379
380 #ifdef __INTERSECTION_REFINE__
381         if(isect->object != OBJECT_NONE) {
382                 if(UNLIKELY(t == 0.0f)) {
383                         return P;
384                 }
385 #  ifdef __OBJECT_MOTION__
386                 Transform tfm = ccl_fetch(sd, ob_itfm);
387 #  else
388                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
389 #  endif
390
391                 P = transform_point(&tfm, P);
392                 D = transform_direction(&tfm, D*t);
393                 D = normalize_len(D, &t);
394         }
395
396         P = P + D*t;
397
398         const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, isect->prim);
399         const float4 tri_a = kernel_tex_fetch(__prim_tri_verts, tri_vindex+0),
400                      tri_b = kernel_tex_fetch(__prim_tri_verts, tri_vindex+1),
401                      tri_c = kernel_tex_fetch(__prim_tri_verts, tri_vindex+2);
402         float3 edge1 = make_float3(tri_a.x - tri_c.x, tri_a.y - tri_c.y, tri_a.z - tri_c.z);
403         float3 edge2 = make_float3(tri_b.x - tri_c.x, tri_b.y - tri_c.y, tri_b.z - tri_c.z);
404         float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
405         float3 qvec = cross(tvec, edge1);
406         float3 pvec = cross(D, edge2);
407         float det = dot(edge1, pvec);
408         if(det != 0.0f) {
409                 /* If determinant is zero it means ray lies in the plane of
410                  * the triangle. It is possible in theory due to watertight
411                  * nature of triangle intersection. For such cases we simply
412                  * don't refine intersection hoping it'll go all fine.
413                  */
414                 float rt = dot(edge2, qvec) / det;
415                 P = P + D*rt;
416         }
417
418         if(isect->object != OBJECT_NONE) {
419 #  ifdef __OBJECT_MOTION__
420                 Transform tfm = ccl_fetch(sd, ob_tfm);
421 #  else
422                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
423 #  endif
424
425                 P = transform_point(&tfm, P);
426         }
427
428         return P;
429 #else
430         return P + D*t;
431 #endif
432 }
433
434 /* Same as above, except that isect->t is assumed to be in object space for
435  * instancing.
436  */
437 ccl_device_inline float3 triangle_refine_subsurface(KernelGlobals *kg,
438                                                     ShaderData *sd,
439                                                     const Intersection *isect,
440                                                     const Ray *ray)
441 {
442         float3 P = ray->P;
443         float3 D = ray->D;
444         float t = isect->t;
445
446         if(isect->object != OBJECT_NONE) {
447 #ifdef __OBJECT_MOTION__
448                 Transform tfm = ccl_fetch(sd, ob_itfm);
449 #else
450                 Transform tfm = object_fetch_transform(kg,
451                                                        isect->object,
452                                                        OBJECT_INVERSE_TRANSFORM);
453 #endif
454
455                 P = transform_point(&tfm, P);
456                 D = transform_direction(&tfm, D);
457                 D = normalize(D);
458         }
459
460         P = P + D*t;
461
462 #ifdef __INTERSECTION_REFINE__
463         const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, isect->prim);
464         const float4 tri_a = kernel_tex_fetch(__prim_tri_verts, tri_vindex+0),
465                      tri_b = kernel_tex_fetch(__prim_tri_verts, tri_vindex+1),
466                      tri_c = kernel_tex_fetch(__prim_tri_verts, tri_vindex+2);
467         float3 edge1 = make_float3(tri_a.x - tri_c.x, tri_a.y - tri_c.y, tri_a.z - tri_c.z);
468         float3 edge2 = make_float3(tri_b.x - tri_c.x, tri_b.y - tri_c.y, tri_b.z - tri_c.z);
469         float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
470         float3 qvec = cross(tvec, edge1);
471         float3 pvec = cross(D, edge2);
472         float det = dot(edge1, pvec);
473         if(det != 0.0f) {
474                 /* If determinant is zero it means ray lies in the plane of
475                  * the triangle. It is possible in theory due to watertight
476                  * nature of triangle intersection. For such cases we simply
477                  * don't refine intersection hoping it'll go all fine.
478                  */
479                 float rt = dot(edge2, qvec) / det;
480                 P = P + D*rt;
481         }
482 #endif  /* __INTERSECTION_REFINE__ */
483
484         if(isect->object != OBJECT_NONE) {
485 #ifdef __OBJECT_MOTION__
486                 Transform tfm = ccl_fetch(sd, ob_tfm);
487 #else
488                 Transform tfm = object_fetch_transform(kg,
489                                                        isect->object,
490                                                        OBJECT_TRANSFORM);
491 #endif
492
493                 P = transform_point(&tfm, P);
494         }
495
496         return P;
497 }
498
499 #undef IDX
500
501 CCL_NAMESPACE_END