368ea25d94103543a04304a0eb40e28da6417c88
[blender.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 ccl_device_inline bool triangle_intersect(KernelGlobals *kg,
26                                           Intersection *isect,
27                                           float3 P,
28                                           float3 dir,
29                                           uint visibility,
30                                           int object,
31                                           int prim_addr)
32 {
33         const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, prim_addr);
34 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
35         const ssef *ssef_verts = (ssef*)&kg->__prim_tri_verts.data[tri_vindex];
36 #else
37         const float4 tri_a = kernel_tex_fetch(__prim_tri_verts, tri_vindex+0),
38                      tri_b = kernel_tex_fetch(__prim_tri_verts, tri_vindex+1),
39                      tri_c = kernel_tex_fetch(__prim_tri_verts, tri_vindex+2);
40 #endif
41         float t, u, v;
42         if(ray_triangle_intersect(P,
43                                   dir,
44                                   isect->t,
45 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
46                                   ssef_verts,
47 #else
48                                   float4_to_float3(tri_a),
49                                   float4_to_float3(tri_b),
50                                   float4_to_float3(tri_c),
51 #endif
52                                   &u, &v, &t))
53         {
54 #ifdef __VISIBILITY_FLAG__
55                 /* Visibility flag test. we do it here under the assumption
56                  * that most triangles are culled by node flags.
57                  */
58                 if(kernel_tex_fetch(__prim_visibility, prim_addr) & visibility)
59 #endif
60                 {
61                         isect->prim = prim_addr;
62                         isect->object = object;
63                         isect->type = PRIMITIVE_TRIANGLE;
64                         isect->u = u;
65                         isect->v = v;
66                         isect->t = t;
67                         return true;
68                 }
69         }
70         return false;
71 }
72
73 #ifdef __KERNEL_AVX2__
74 #define cross256(A,B, C,D) _mm256_fmsub_ps(A,B, _mm256_mul_ps(C,D))
75 ccl_device_inline int ray_triangle_intersect8(
76             KernelGlobals *kg,
77             float3 ray_P,
78             float3 ray_dir,
79             Intersection **isect,
80             uint visibility,
81             int object,
82             __m256 *triA,
83             __m256 *triB,
84             __m256 *triC,
85             int prim_addr,
86             int prim_num,
87             uint *num_hits,
88             uint max_hits,
89             int *num_hits_in_instance,
90             float isec_t)
91 {
92
93         const unsigned char prim_num_mask = (1 << prim_num) - 1;
94
95         const __m256i zero256 = _mm256_setzero_si256();
96
97         const __m256 Px256 = _mm256_set1_ps(ray_P.x);
98         const __m256 Py256 = _mm256_set1_ps(ray_P.y);
99         const __m256 Pz256 = _mm256_set1_ps(ray_P.z);
100
101         const __m256 dirx256 = _mm256_set1_ps(ray_dir.x);
102         const __m256 diry256 = _mm256_set1_ps(ray_dir.y);
103         const __m256 dirz256 = _mm256_set1_ps(ray_dir.z);
104
105         /* Calculate vertices relative to ray origin. */
106         /*      const float3 v0 = tri_c - P;
107         const float3 v1 = tri_a - P;
108         const float3 v2 = tri_b - P; */
109
110         __m256 v0_x_256 = _mm256_sub_ps(triC[0], Px256);
111         __m256 v0_y_256 = _mm256_sub_ps(triC[1], Py256);
112         __m256 v0_z_256 = _mm256_sub_ps(triC[2], Pz256);
113
114         __m256 v1_x_256 = _mm256_sub_ps(triA[0], Px256);
115         __m256 v1_y_256 = _mm256_sub_ps(triA[1], Py256);
116         __m256 v1_z_256 = _mm256_sub_ps(triA[2], Pz256);
117
118         __m256 v2_x_256 = _mm256_sub_ps(triB[0], Px256);
119         __m256 v2_y_256 = _mm256_sub_ps(triB[1], Py256);
120         __m256 v2_z_256 = _mm256_sub_ps(triB[2], Pz256);
121
122         __m256 v0_v1_x_256 = _mm256_add_ps(v0_x_256, v1_x_256);
123         __m256 v0_v1_y_256 = _mm256_add_ps(v0_y_256, v1_y_256);
124         __m256 v0_v1_z_256 = _mm256_add_ps(v0_z_256, v1_z_256);
125
126         __m256 v0_v2_x_256 = _mm256_add_ps(v0_x_256, v2_x_256);
127         __m256 v0_v2_y_256 = _mm256_add_ps(v0_y_256, v2_y_256);
128         __m256 v0_v2_z_256 = _mm256_add_ps(v0_z_256, v2_z_256);
129
130         __m256 v1_v2_x_256 = _mm256_add_ps(v1_x_256, v2_x_256);
131         __m256 v1_v2_y_256 = _mm256_add_ps(v1_y_256, v2_y_256);
132         __m256 v1_v2_z_256 = _mm256_add_ps(v1_z_256, v2_z_256);
133
134         /* Calculate triangle edges.
135         const float3 e0 = v2 - v0;
136         const float3 e1 = v0 - v1;
137         const float3 e2 = v1 - v2;*/
138
139         __m256 e0_x_256 = _mm256_sub_ps(v2_x_256, v0_x_256);
140         __m256 e0_y_256 = _mm256_sub_ps(v2_y_256, v0_y_256);
141         __m256 e0_z_256 = _mm256_sub_ps(v2_z_256, v0_z_256);
142
143         __m256 e1_x_256 = _mm256_sub_ps(v0_x_256, v1_x_256);
144         __m256 e1_y_256 = _mm256_sub_ps(v0_y_256, v1_y_256);
145         __m256 e1_z_256 = _mm256_sub_ps(v0_z_256, v1_z_256);
146
147         __m256 e2_x_256 = _mm256_sub_ps(v1_x_256, v2_x_256);
148         __m256 e2_y_256 = _mm256_sub_ps(v1_y_256, v2_y_256);
149         __m256 e2_z_256 = _mm256_sub_ps(v1_z_256, v2_z_256);
150
151         /* Perform edge tests.
152         const float U = dot(cross(v2 + v0, e0), ray_dir);
153         const float V = dot(cross(v0 + v1, e1), ray_dir);
154         const float W = dot(cross(v1 + v2, e2), ray_dir);*/
155
156         //cross (AyBz - AzBy, AzBx -AxBz,  AxBy - AyBx)
157         __m256 U_x_256 = cross256(v0_v2_y_256, e0_z_256, v0_v2_z_256, e0_y_256);
158         __m256 U_y_256 = cross256(v0_v2_z_256, e0_x_256, v0_v2_x_256, e0_z_256);
159         __m256 U_z_256 = cross256(v0_v2_x_256, e0_y_256, v0_v2_y_256, e0_x_256);
160         //vertical dot
161         __m256 U_256 = _mm256_mul_ps(U_x_256, dirx256);
162         U_256 = _mm256_fmadd_ps(U_y_256, diry256, U_256); //_mm256_add_ps(U_256, _mm256_mul_ps(U_y_256, diry256));
163         U_256 = _mm256_fmadd_ps(U_z_256, dirz256, U_256); //_mm256_add_ps(U_256, _mm256_mul_ps(U_z_256, dirz256));
164
165         __m256 V_x_256 = cross256(v0_v1_y_256, e1_z_256, v0_v1_z_256, e1_y_256);
166         __m256 V_y_256 = cross256(v0_v1_z_256, e1_x_256, v0_v1_x_256, e1_z_256);
167         __m256 V_z_256 = cross256(v0_v1_x_256, e1_y_256, v0_v1_y_256, e1_x_256);
168         //vertical dot
169         __m256 V_256 = _mm256_mul_ps(V_x_256, dirx256);
170         V_256 = _mm256_fmadd_ps(V_y_256, diry256, V_256);// _mm256_add_ps(V_256, _mm256_mul_ps(V_y_256, diry256));
171         V_256 = _mm256_fmadd_ps(V_z_256, dirz256, V_256);// _mm256_add_ps(V_256, _mm256_mul_ps(V_z_256, dirz256));
172
173         __m256 W_x_256 = cross256(v1_v2_y_256, e2_z_256, v1_v2_z_256, e2_y_256);
174         __m256 W_y_256 = cross256(v1_v2_z_256, e2_x_256, v1_v2_x_256, e2_z_256);
175         __m256 W_z_256 = cross256(v1_v2_x_256, e2_y_256, v1_v2_y_256, e2_x_256);
176         //vertical dot
177         __m256 W_256 = _mm256_mul_ps(W_x_256, dirx256);
178         W_256 = _mm256_fmadd_ps(W_y_256, diry256,W_256);//_mm256_add_ps(W_256, _mm256_mul_ps(W_y_256, diry256));
179         W_256 = _mm256_fmadd_ps(W_z_256, dirz256,W_256);//_mm256_add_ps(W_256, _mm256_mul_ps(W_z_256, dirz256));
180
181         //const float minUVW = min(U, min(V, W));
182         //const float maxUVW = max(U, max(V, W));
183 #if 0
184         __m256 minUVW_256 = _mm256_min_ps(U_256, _mm256_min_ps(V_256, W_256));
185         __m256 maxUVW_256 = _mm256_max_ps(U_256, _mm256_max_ps(V_256, W_256));
186
187         //if(minUVW < 0.0f && maxUVW > 0.0f)
188         __m256i mask_minmaxUVW_256 = _mm256_and_si256(
189                 _mm256_cmpgt_epi32(zero256, _mm256_castps_si256(minUVW_256)),
190                 //_mm256_castps_si256(minUVW_256),
191                 _mm256_cmpgt_epi32(_mm256_castps_si256(maxUVW_256), zero256));
192 #else
193         __m256i U_256_1 = _mm256_srli_epi32(_mm256_castps_si256(U_256), 31);
194         __m256i V_256_1 = _mm256_srli_epi32(_mm256_castps_si256(V_256), 31);
195         __m256i W_256_1 = _mm256_srli_epi32(_mm256_castps_si256(W_256), 31);
196         __m256i UVW_256_1 = _mm256_add_epi32(_mm256_add_epi32(U_256_1, V_256_1), W_256_1);
197
198         const __m256i one256 = _mm256_set1_epi32(1);
199         const __m256i two256 = _mm256_set1_epi32(2);
200
201         __m256i mask_minmaxUVW_256 = _mm256_or_si256(
202                 _mm256_cmpeq_epi32(one256, UVW_256_1),
203                 _mm256_cmpeq_epi32(two256, UVW_256_1) );
204 #endif
205
206         unsigned char mask_minmaxUVW_pos = _mm256_movemask_ps(_mm256_castsi256_ps(mask_minmaxUVW_256));
207         if((mask_minmaxUVW_pos & prim_num_mask) == prim_num_mask) { //all bits set
208                 return false;
209         }
210
211         /* Calculate geometry normal and denominator. */
212         //                      const float3 Ng1 = cross(e1, e0);
213         //const Vec3vfM Ng1 = stable_triangle_normal(e2,e1,e0);
214
215         __m256 Ng1_x_256 = cross256(e1_y_256, e0_z_256, e1_z_256, e0_y_256);
216         __m256 Ng1_y_256 = cross256(e1_z_256, e0_x_256, e1_x_256, e0_z_256);
217         __m256 Ng1_z_256 = cross256(e1_x_256, e0_y_256, e1_y_256, e0_x_256);
218
219         //const float3 Ng = Ng1 + Ng1;
220         Ng1_x_256 = _mm256_add_ps(Ng1_x_256, Ng1_x_256);
221         Ng1_y_256 = _mm256_add_ps(Ng1_y_256, Ng1_y_256);
222         Ng1_z_256 = _mm256_add_ps(Ng1_z_256, Ng1_z_256);
223
224         //const float den = dot3(Ng, dir);
225         //vertical dot
226         __m256 den_256 = _mm256_mul_ps(Ng1_x_256, dirx256);
227         den_256 = _mm256_fmadd_ps(Ng1_y_256, diry256,den_256);//_mm256_add_ps(den_256, _mm256_mul_ps(Ng1_y_256, diry256));
228         den_256 = _mm256_fmadd_ps(Ng1_z_256, dirz256,den_256);//_mm256_add_ps(den_256, _mm256_mul_ps(Ng1_z_256, dirz256));
229
230         // __m256i maskden256 = _mm256_cmpeq_epi32(_mm256_castps_si256(den_256), zero256);
231
232         /* Perform depth test. */
233         //const float T = dot3(v0, Ng);
234         __m256 T_256 = _mm256_mul_ps(Ng1_x_256, v0_x_256);
235         T_256 = _mm256_fmadd_ps(Ng1_y_256, v0_y_256,T_256);//_mm256_add_ps(T_256, _mm256_mul_ps(Ng1_y_256, v0_y_256));
236         T_256 = _mm256_fmadd_ps(Ng1_z_256, v0_z_256,T_256);//_mm256_add_ps(T_256, _mm256_mul_ps(Ng1_z_256, v0_z_256));
237
238         //const int sign_den = (__float_as_int(den) & 0x80000000);
239         const __m256i c0x80000000 = _mm256_set1_epi32(0x80000000);
240         __m256i sign_den_256 = _mm256_and_si256(_mm256_castps_si256(den_256), c0x80000000);
241
242         //const float sign_T = xor_signmask(T, sign_den);
243         __m256 sign_T_256 = _mm256_castsi256_ps(_mm256_xor_si256(_mm256_castps_si256(T_256), sign_den_256));
244
245         /*if((sign_T < 0.0f) || mask_minmaxUVW_pos {    return false;}  */
246         unsigned char mask_sign_T = _mm256_movemask_ps(sign_T_256);
247         if(((mask_minmaxUVW_pos | mask_sign_T) & prim_num_mask) == prim_num_mask) {
248                 return false;
249         } /**/
250
251         __m256 xor_signmask_256 = _mm256_castsi256_ps(_mm256_xor_si256(_mm256_castps_si256(den_256), sign_den_256));
252
253
254         ccl_align(32) float den8[8], U8[8], V8[8], T8[8], sign_T8[8], xor_signmask8[8];
255         ccl_align(32) unsigned int mask_minmaxUVW8[8];
256
257         if(visibility == PATH_RAY_SHADOW_OPAQUE){
258                         __m256i mask_final_256 = _mm256_cmpeq_epi32(mask_minmaxUVW_256, zero256);//~mask_minmaxUVW_256
259
260                         __m256i maskden256 = _mm256_cmpeq_epi32(_mm256_castps_si256(den_256), zero256);
261
262                         __m256i mask0 = _mm256_cmpgt_epi32(zero256, _mm256_castps_si256(sign_T_256));
263                         __m256 rayt_256 = _mm256_set1_ps((*isect)->t);
264
265                         __m256i mask1 = _mm256_cmpgt_epi32(_mm256_castps_si256(sign_T_256),
266                                 _mm256_castps_si256(
267                                         _mm256_mul_ps(_mm256_castsi256_ps(_mm256_xor_si256(_mm256_castps_si256(den_256), sign_den_256)), rayt_256)
268                                 )
269                         );
270                         /*      __m256i mask1 = _mm256_castps_si256(_mm256_cmp_ps(sign_T_256,
271                         _mm256_mul_ps(_mm256_castsi256_ps(_mm256_xor_si256(_mm256_castps_si256(den_256), sign_den_256)), rayt_256),
272                         _CMP_GT_OS
273                         ) );*/
274
275                         mask0 = _mm256_or_si256(mask1, mask0);
276                         //unsigned char mask = _mm256_movemask_ps(_mm256_castsi256_ps(mask0));
277                         //unsigned char maskden = _mm256_movemask_ps(_mm256_castsi256_ps(maskden256));
278                         //unsigned char mask_final = ((~mask) & (~maskden) & (~mask_minmaxUVW_pos));
279                         mask_final_256 = _mm256_andnot_si256(mask0, mask_final_256); //(~mask_minmaxUVW_pos) &(~mask)
280                         mask_final_256 = _mm256_andnot_si256(maskden256, mask_final_256); //(~mask_minmaxUVW_pos) &(~mask) & (~maskden)
281
282                         unsigned char mask_final = _mm256_movemask_ps(_mm256_castsi256_ps(mask_final_256));
283                         if((mask_final & prim_num_mask) == 0) { //all bits NOT set
284                                 return false;
285                         }               /**/
286
287                         unsigned long i = 0;
288 #if defined(_MSC_VER)
289                         unsigned char res = _BitScanForward(&i, (unsigned long)mask_final);
290 #else
291             i = __builtin_ffs(mask_final)-1;
292 #endif
293
294                         den_256 = _mm256_rcp_ps(den_256); //inv_den
295                         U_256 = _mm256_mul_ps(U_256, den_256); //*inv_den
296                         V_256 = _mm256_mul_ps(V_256, den_256); //*inv_den
297                         T_256 = _mm256_mul_ps(T_256, den_256); //*inv_den
298
299                         _mm256_store_ps(U8, U_256);
300                         _mm256_store_ps(V8, V_256);
301                         _mm256_store_ps(T8, T_256);
302
303
304                         //here we assume (kernel_tex_fetch(__prim_visibility, (prim_addr +i)) & visibility) is always true
305
306                         (*isect)->u = U8[i];
307                         (*isect)->v = V8[i];
308                         (*isect)->t = T8[i];
309
310                         (*isect)->prim = (prim_addr + i);
311                         (*isect)->object = object;
312                         (*isect)->type = PRIMITIVE_TRIANGLE;
313
314                         return true;
315                 }
316         else {
317                         _mm256_store_ps(den8, den_256);
318                         _mm256_store_ps(U8, U_256);
319                         _mm256_store_ps(V8, V_256);
320                         _mm256_store_ps(T8, T_256);
321
322                         _mm256_store_ps(sign_T8, sign_T_256);
323                         _mm256_store_ps(xor_signmask8, xor_signmask_256);
324                         _mm256_store_si256((__m256i*)mask_minmaxUVW8, mask_minmaxUVW_256);
325
326                         int ret = false;
327
328                         if(visibility == PATH_RAY_SHADOW) {
329                                 for(int i = 0; i < prim_num; i++) {
330                                         if(!mask_minmaxUVW8[i]) {
331 #ifdef __VISIBILITY_FLAG__
332                                                 if(kernel_tex_fetch(__prim_visibility, (prim_addr + i)) & visibility)
333 #endif
334                                                 {
335                                                         if((sign_T8[i] >= 0.0f) &&
336                                                            (sign_T8[i] <= (*isect)->t * xor_signmask8[i]))
337                                                         {
338                                                                 if(den8[i]) {
339                                                                         const float inv_den = 1.0f / den8[i];
340
341                                                                         (*isect)->u = U8[i] * inv_den;
342                                                                         (*isect)->v = V8[i] * inv_den;
343                                                                         (*isect)->t = T8[i] * inv_den;
344
345                                                                         (*isect)->prim = (prim_addr + i);
346                                                                         (*isect)->object = object;
347                                                                         (*isect)->type = PRIMITIVE_TRIANGLE;
348
349                                                                         int prim = kernel_tex_fetch(__prim_index, (*isect)->prim);
350                                                                         int shader = 0;
351
352 #ifdef __HAIR__
353                                                                         if(kernel_tex_fetch(__prim_type, (*isect)->prim) & PRIMITIVE_ALL_TRIANGLE)
354 #endif
355                                                                         {
356                                                                                 shader = kernel_tex_fetch(__tri_shader, prim);
357                                                                         }
358 #ifdef __HAIR__
359                                                                         else {
360                                                                                 float4 str = kernel_tex_fetch(__curves, prim);
361                                                                                 shader = __float_as_int(str.z);
362                                                                         }
363 #endif
364                                                                         int flag = kernel_tex_fetch(__shaders, (shader & SHADER_MASK)).flags;
365
366                                                                         /* if no transparent shadows, all light is blocked */
367                                                                         if(!(flag & SD_HAS_TRANSPARENT_SHADOW)) {
368                                                                                 return 2;
369                                                                         }
370                                                                         /* if maximum number of hits reached, block all light */
371                                                                         else if(*num_hits == max_hits) {
372                                                                                 return 2;
373                                                                         }
374                                                                         /* move on to next entry in intersections array */
375                                                                         ret = true;
376
377                                                                         (*isect)++;
378                                                                         (*num_hits)++;
379
380                                                                         (*num_hits_in_instance)++;
381
382                                                                         (*isect)->t = isec_t;
383
384                                                                 } //den
385                                                         } //if sign
386                                                 } //vis
387                                         }//if mask
388                                 } //for
389                 }
390                 else { //default case
391                         for(int i = 0; i < prim_num; i++) {
392                                 if(!mask_minmaxUVW8[i]) {
393 #ifdef __VISIBILITY_FLAG__
394                                         if(kernel_tex_fetch(__prim_visibility, (prim_addr + i)) & visibility)
395 #endif
396                                         {
397                                                 if((sign_T8[i] >= 0.0f) &&
398                                                    (sign_T8[i] <= (*isect)->t * xor_signmask8[i]))
399                                                 {
400                                                         if(den8[i]) {
401                                                                 const float inv_den = 1.0f / den8[i];
402
403                                                                 (*isect)->u = U8[i] * inv_den;
404                                                                 (*isect)->v = V8[i] * inv_den;
405                                                                 (*isect)->t = T8[i] * inv_den;
406
407                                                                 (*isect)->prim = (prim_addr + i);
408                                                                 (*isect)->object = object;
409                                                                 (*isect)->type = PRIMITIVE_TRIANGLE;
410
411                                                                 ret = true;
412                                                         } //den
413                                                 } //if sign
414                                         } //vis
415                                 }//if mask
416                         } //for
417                 } //default
418         return ret;
419 }// else PATH_RAY_SHADOW_OPAQUE
420
421 }
422
423 ccl_device_inline int triangle_intersect8(
424         KernelGlobals *kg,
425         Intersection **isect,
426         float3 P,
427         float3 dir,
428         uint visibility,
429         int object,
430         int prim_addr,
431         int prim_num,
432         uint *num_hits,
433         uint max_hits,
434         int *num_hits_in_instance,
435         float isec_t)
436  {
437         __m128 tri_a[8], tri_b[8], tri_c[8];
438         __m256  tritmp[12], tri[12];
439         __m256  triA[3], triB[3], triC[3];
440
441         int i, r;
442
443         uint tri_vindex = kernel_tex_fetch(__prim_tri_index, prim_addr);
444         for(i = 0; i < prim_num; i++) {
445                 tri_a[i] = *(__m128*)&kg->__prim_tri_verts.data[tri_vindex++];
446                 tri_b[i] = *(__m128*)&kg->__prim_tri_verts.data[tri_vindex++];
447                 tri_c[i] = *(__m128*)&kg->__prim_tri_verts.data[tri_vindex++];
448          }
449         //create 9 or  12 placeholders
450         tri[0] = _mm256_castps128_ps256(tri_a[0]);    //_mm256_zextps128_ps256
451         tri[1] = _mm256_castps128_ps256(tri_b[0]);//_mm256_zextps128_ps256
452         tri[2] = _mm256_castps128_ps256(tri_c[0]);//_mm256_zextps128_ps256
453
454         tri[3] = _mm256_castps128_ps256(tri_a[1]);    //_mm256_zextps128_ps256
455         tri[4] = _mm256_castps128_ps256(tri_b[1]);//_mm256_zextps128_ps256
456         tri[5] = _mm256_castps128_ps256(tri_c[1]);//_mm256_zextps128_ps256
457
458         tri[6] = _mm256_castps128_ps256(tri_a[2]);    //_mm256_zextps128_ps256
459         tri[7] = _mm256_castps128_ps256(tri_b[2]);//_mm256_zextps128_ps256
460         tri[8] = _mm256_castps128_ps256(tri_c[2]);//_mm256_zextps128_ps256
461
462         if(prim_num > 3) {
463                 tri[9] =  _mm256_castps128_ps256(tri_a[3]);    //_mm256_zextps128_ps256
464                 tri[10] = _mm256_castps128_ps256(tri_b[3]);//_mm256_zextps128_ps256
465                 tri[11] = _mm256_castps128_ps256(tri_c[3]);//_mm256_zextps128_ps256
466         }
467
468         for(i = 4, r = 0; i < prim_num; i ++, r += 3) {
469                 tri[r] =     _mm256_insertf128_ps(tri[r] , tri_a[i], 1);
470                 tri[r + 1] = _mm256_insertf128_ps(tri[r + 1], tri_b[i], 1);
471                 tri[r + 2] = _mm256_insertf128_ps(tri[r + 2], tri_c[i], 1);
472          }
473
474         //------------------------------------------------
475         //0!  Xa0 Ya0 Za0 1 Xa4 Ya4 Za4  1
476         //1!  Xb0 Yb0 Zb0 1 Xb4 Yb4 Zb4 1
477         //2!  Xc0 Yc0 Zc0 1 Xc4 Yc4 Zc4 1
478
479         //3!  Xa1 Ya1 Za1 1 Xa5 Ya5 Za5 1
480         //4!  Xb1 Yb1 Zb1 1 Xb5 Yb5 Zb5  1
481         //5!  Xc1 Yc1 Zc1 1 Xc5 Yc5 Zc5 1
482
483         //6!  Xa2 Ya2 Za2 1 Xa6 Ya6 Za6 1
484         //7!  Xb2 Yb2 Zb2 1 Xb6 Yb6 Zb6  1
485         //8!  Xc2 Yc2 Zc2 1 Xc6 Yc6 Zc6 1
486
487         //9!  Xa3 Ya3 Za3 1 Xa7 Ya7 Za7  1
488         //10! Xb3 Yb3 Zb3 1 Xb7 Yb7 Zb7  1
489         //11! Xc3 Yc3 Zc3 1 Xc7 Yc7 Zc7  1
490
491         //"transpose"
492         tritmp[0] = _mm256_unpacklo_ps(tri[0], tri[3]);   //0!  Xa0 Xa1 Ya0 Ya1 Xa4 Xa5 Ya4 Ya5
493         tritmp[1] = _mm256_unpackhi_ps(tri[0], tri[3]);   //1!  Za0 Za1 1   1   Za4 Za5  1   1
494
495         tritmp[2] = _mm256_unpacklo_ps(tri[6], tri[9]);   //2!  Xa2 Xa3 Ya2 Ya3 Xa6 Xa7 Ya6 Ya7
496         tritmp[3] = _mm256_unpackhi_ps(tri[6], tri[9]);   //3!  Za2 Za3  1   1  Za6 Za7  1   1
497
498         tritmp[4] = _mm256_unpacklo_ps(tri[1], tri[4]);   //4!  Xb0 Xb1 Yb0 Yb1 Xb4 Xb5 Yb4 Yb5
499         tritmp[5] = _mm256_unpackhi_ps(tri[1], tri[4]);   //5!  Zb0 Zb1  1  1   Zb4 Zb5  1   1
500
501         tritmp[6] = _mm256_unpacklo_ps(tri[7], tri[10]);  //6!  Xb2 Xb3 Yb2 Yb3 Xb6 Xb7 Yb6 Yb7
502         tritmp[7] = _mm256_unpackhi_ps(tri[7], tri[10]);  //7!  Zb2 Zb3  1    1 Zb6 Zb7  1   1
503
504         tritmp[8] = _mm256_unpacklo_ps(tri[2], tri[5]);   //8!  Xc0 Xc1 Yc0 Yc1 Xc4 Xc5 Yc4 Yc5
505         tritmp[9] = _mm256_unpackhi_ps(tri[2], tri[5]);   //9!  Zc0 Zc1  1   1  Zc4 Zc5  1   1
506
507         tritmp[10] = _mm256_unpacklo_ps(tri[8], tri[11]); //10! Xc2 Xc3 Yc2 Yc3 Xc6 Xc7 Yc6 Yc7
508         tritmp[11] = _mm256_unpackhi_ps(tri[8], tri[11]); //11! Zc2 Zc3  1   1  Zc6 Zc7  1   1
509
510                                 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
511         triA[0] = _mm256_castpd_ps(_mm256_unpacklo_pd(_mm256_castps_pd(tritmp[0]), _mm256_castps_pd(tritmp[2])));       //  Xa0 Xa1 Xa2 Xa3 Xa4 Xa5 Xa6 Xa7
512         triA[1] = _mm256_castpd_ps(_mm256_unpackhi_pd(_mm256_castps_pd(tritmp[0]), _mm256_castps_pd(tritmp[2])));   //  Ya0 Ya1 Ya2 Ya3 Ya4 Ya5 Ya6 Ya7
513         triA[2] = _mm256_castpd_ps(_mm256_unpacklo_pd(_mm256_castps_pd(tritmp[1]), _mm256_castps_pd(tritmp[3])));   //  Za0 Za1 Za2 Za3 Za4 Za5 Za6 Za7
514
515         triB[0] = _mm256_castpd_ps(_mm256_unpacklo_pd(_mm256_castps_pd(tritmp[4]), _mm256_castps_pd(tritmp[6])));       //  Xb0 Xb1  Xb2 Xb3 Xb4 Xb5 Xb5 Xb7
516         triB[1] = _mm256_castpd_ps(_mm256_unpackhi_pd(_mm256_castps_pd(tritmp[4]), _mm256_castps_pd(tritmp[6])));   //  Yb0 Yb1  Yb2 Yb3 Yb4 Yb5 Yb5 Yb7
517         triB[2] = _mm256_castpd_ps(_mm256_unpacklo_pd(_mm256_castps_pd(tritmp[5]), _mm256_castps_pd(tritmp[7]))); //    Zb0 Zb1  Zb2 Zb3 Zb4 Zb5 Zb5 Zb7
518
519         triC[0] = _mm256_castpd_ps(_mm256_unpacklo_pd(_mm256_castps_pd(tritmp[8]), _mm256_castps_pd(tritmp[10])));     //Xc0 Xc1 Xc2 Xc3 Xc4 Xc5 Xc6 Xc7
520         triC[1] = _mm256_castpd_ps(_mm256_unpackhi_pd(_mm256_castps_pd(tritmp[8]), _mm256_castps_pd(tritmp[10])));     //Yc0 Yc1 Yc2 Yc3 Yc4 Yc5 Yc6 Yc7
521         triC[2] = _mm256_castpd_ps(_mm256_unpacklo_pd(_mm256_castps_pd(tritmp[9]), _mm256_castps_pd(tritmp[11])));     //Zc0 Zc1 Zc2 Zc3 Zc4 Zc5 Zc6 Zc7
522
523                           /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
524
525         int result = ray_triangle_intersect8(kg, P,
526                                              dir,
527                                              isect,
528                                              visibility, object,
529                                              triA,
530                                              triB,
531                                              triC,
532                                              prim_addr,
533                                              prim_num,
534                                              num_hits,
535                                              max_hits,
536                                              num_hits_in_instance,
537                                              isec_t);
538         return result;
539 }
540
541 #endif  /* __KERNEL_AVX2__ */
542
543 /* Special ray intersection routines for subsurface scattering. In that case we
544  * only want to intersect with primitives in the same object, and if case of
545  * multiple hits we pick a single random primitive as the intersection point.
546  * Returns whether traversal should be stopped.
547  */
548
549 #ifdef __BVH_LOCAL__
550 ccl_device_inline bool triangle_intersect_local(
551         KernelGlobals *kg,
552         LocalIntersection *local_isect,
553         float3 P,
554         float3 dir,
555         int object,
556         int local_object,
557         int prim_addr,
558         float tmax,
559         uint *lcg_state,
560         int max_hits)
561 {
562         /* Only intersect with matching object, for instanced objects we
563          * already know we are only intersecting the right object. */
564         if(object == OBJECT_NONE) {
565                 if(kernel_tex_fetch(__prim_object, prim_addr) != local_object) {
566                         return false;
567                 }
568         }
569
570         const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, prim_addr);
571 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
572         const ssef *ssef_verts = (ssef*)&kg->__prim_tri_verts.data[tri_vindex];
573 #else
574         const float3 tri_a = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex+0)),
575                      tri_b = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex+1)),
576                      tri_c = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex+2));
577 #endif
578         float t, u, v;
579         if(!ray_triangle_intersect(P,
580                                    dir,
581                                    tmax,
582 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
583                                    ssef_verts,
584 #else
585                                    tri_a, tri_b, tri_c,
586 #endif
587                                    &u, &v, &t))
588         {
589                 return false;
590         }
591
592         /* If no actual hit information is requested, just return here. */
593         if(max_hits == 0) {
594                 return true;
595         }
596
597         int hit;
598         if(lcg_state) {
599                 /* Record up to max_hits intersections. */
600                 for(int i = min(max_hits, local_isect->num_hits) - 1; i >= 0; --i) {
601                         if(local_isect->hits[i].t == t) {
602                                 return false;
603                         }
604                 }
605
606                 local_isect->num_hits++;
607
608                 if(local_isect->num_hits <= max_hits) {
609                         hit = local_isect->num_hits - 1;
610                 }
611                 else {
612                         /* reservoir sampling: if we are at the maximum number of
613                          * hits, randomly replace element or skip it */
614                         hit = lcg_step_uint(lcg_state) % local_isect->num_hits;
615
616                         if(hit >= max_hits)
617                                 return false;
618                 }
619         }
620         else {
621                 /* Record closest intersection only. */
622                 if(local_isect->num_hits && t > local_isect->hits[0].t) {
623                         return false;
624                 }
625
626                 hit = 0;
627                 local_isect->num_hits = 1;
628         }
629
630         /* Record intersection. */
631         Intersection *isect = &local_isect->hits[hit];
632         isect->prim = prim_addr;
633         isect->object = object;
634         isect->type = PRIMITIVE_TRIANGLE;
635         isect->u = u;
636         isect->v = v;
637         isect->t = t;
638
639         /* Record geometric normal. */
640 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
641         const float3 tri_a = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex+0)),
642                      tri_b = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex+1)),
643                      tri_c = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex+2));
644 #endif
645         local_isect->Ng[hit] = normalize(cross(tri_b - tri_a, tri_c - tri_a));
646
647         return false;
648 }
649 #endif  /* __BVH_LOCAL__ */
650
651 /* Refine triangle intersection to more precise hit point. For rays that travel
652  * far the precision is often not so good, this reintersects the primitive from
653  * a closer distance. */
654
655 /* Reintersections uses the paper:
656  *
657  * Tomas Moeller
658  * Fast, minimum storage ray/triangle intersection
659  * http://www.cs.virginia.edu/~gfx/Courses/2003/ImageSynthesis/papers/Acceleration/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf
660  */
661
662 ccl_device_inline float3 triangle_refine(KernelGlobals *kg,
663                                          ShaderData *sd,
664                                          const Intersection *isect,
665                                          const Ray *ray)
666 {
667         float3 P = ray->P;
668         float3 D = ray->D;
669         float t = isect->t;
670
671 #ifdef __INTERSECTION_REFINE__
672         if(isect->object != OBJECT_NONE) {
673                 if(UNLIKELY(t == 0.0f)) {
674                         return P;
675                 }
676 #  ifdef __OBJECT_MOTION__
677                 Transform tfm = sd->ob_itfm;
678 #  else
679                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
680 #  endif
681
682                 P = transform_point(&tfm, P);
683                 D = transform_direction(&tfm, D*t);
684                 D = normalize_len(D, &t);
685         }
686
687         P = P + D*t;
688
689         const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, isect->prim);
690         const float4 tri_a = kernel_tex_fetch(__prim_tri_verts, tri_vindex+0),
691                      tri_b = kernel_tex_fetch(__prim_tri_verts, tri_vindex+1),
692                      tri_c = kernel_tex_fetch(__prim_tri_verts, tri_vindex+2);
693         float3 edge1 = make_float3(tri_a.x - tri_c.x, tri_a.y - tri_c.y, tri_a.z - tri_c.z);
694         float3 edge2 = make_float3(tri_b.x - tri_c.x, tri_b.y - tri_c.y, tri_b.z - tri_c.z);
695         float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
696         float3 qvec = cross(tvec, edge1);
697         float3 pvec = cross(D, edge2);
698         float det = dot(edge1, pvec);
699         if(det != 0.0f) {
700                 /* If determinant is zero it means ray lies in the plane of
701                  * the triangle. It is possible in theory due to watertight
702                  * nature of triangle intersection. For such cases we simply
703                  * don't refine intersection hoping it'll go all fine.
704                  */
705                 float rt = dot(edge2, qvec) / det;
706                 P = P + D*rt;
707         }
708
709         if(isect->object != OBJECT_NONE) {
710 #  ifdef __OBJECT_MOTION__
711                 Transform tfm = sd->ob_tfm;
712 #  else
713                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
714 #  endif
715
716                 P = transform_point(&tfm, P);
717         }
718
719         return P;
720 #else
721         return P + D*t;
722 #endif
723 }
724
725 /* Same as above, except that isect->t is assumed to be in object space for
726  * instancing.
727  */
728 ccl_device_inline float3 triangle_refine_local(KernelGlobals *kg,
729                                                ShaderData *sd,
730                                                const Intersection *isect,
731                                                const Ray *ray)
732 {
733         float3 P = ray->P;
734         float3 D = ray->D;
735         float t = isect->t;
736
737         if(isect->object != OBJECT_NONE) {
738 #ifdef __OBJECT_MOTION__
739                 Transform tfm = sd->ob_itfm;
740 #else
741                 Transform tfm = object_fetch_transform(kg,
742                                                        isect->object,
743                                                        OBJECT_INVERSE_TRANSFORM);
744 #endif
745
746                 P = transform_point(&tfm, P);
747                 D = transform_direction(&tfm, D);
748                 D = normalize(D);
749         }
750
751         P = P + D*t;
752
753 #ifdef __INTERSECTION_REFINE__
754         const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, isect->prim);
755         const float4 tri_a = kernel_tex_fetch(__prim_tri_verts, tri_vindex+0),
756                      tri_b = kernel_tex_fetch(__prim_tri_verts, tri_vindex+1),
757                      tri_c = kernel_tex_fetch(__prim_tri_verts, tri_vindex+2);
758         float3 edge1 = make_float3(tri_a.x - tri_c.x, tri_a.y - tri_c.y, tri_a.z - tri_c.z);
759         float3 edge2 = make_float3(tri_b.x - tri_c.x, tri_b.y - tri_c.y, tri_b.z - tri_c.z);
760         float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
761         float3 qvec = cross(tvec, edge1);
762         float3 pvec = cross(D, edge2);
763         float det = dot(edge1, pvec);
764         if(det != 0.0f) {
765                 /* If determinant is zero it means ray lies in the plane of
766                  * the triangle. It is possible in theory due to watertight
767                  * nature of triangle intersection. For such cases we simply
768                  * don't refine intersection hoping it'll go all fine.
769                  */
770                 float rt = dot(edge2, qvec) / det;
771                 P = P + D*rt;
772         }
773 #endif  /* __INTERSECTION_REFINE__ */
774
775         if(isect->object != OBJECT_NONE) {
776 #ifdef __OBJECT_MOTION__
777                 Transform tfm = sd->ob_tfm;
778 #else
779                 Transform tfm = object_fetch_transform(kg,
780                                                        isect->object,
781                                                        OBJECT_TRANSFORM);
782 #endif
783
784                 P = transform_point(&tfm, P);
785         }
786
787         return P;
788 }
789
790 CCL_NAMESPACE_END