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