Cycles: Tweak to previous commit, experimental sm_52 works on Linux but not OSX
[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 /* 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_CUDA__)
53 #  if (defined(i386) || defined(_M_IX86))
54 #    if __CUDA_ARCH__ > 500
55 ccl_device_noinline
56 #    else  /* __CUDA_ARCH__ > 500 */
57 ccl_device_inline
58 #    endif  /* __CUDA_ARCH__ > 500 */
59 #  else  /* (defined(i386) || defined(_M_IX86)) */
60 #    if defined(__KERNEL_EXPERIMENTAL__) && (__CUDA_ARCH__ > 500)
61 ccl_device_noinline
62 #    else
63 ccl_device_inline
64 #    endif
65 #  endif  /* (defined(i386) || defined(_M_IX86)) */
66 #elif defined(__KERNEL_OPENCL_APPLE__)
67 ccl_device_noinline
68 #else  /* defined(__KERNEL_OPENCL_APPLE__) */
69 ccl_device_inline
70 #endif  /* defined(__KERNEL_OPENCL_APPLE__) */
71 void triangle_intersect_precalc(float3 dir,
72                                 IsectPrecalc *isect_precalc)
73 {
74         /* Calculate dimension where the ray direction is maximal. */
75         int kz = util_max_axis(make_float3(fabsf(dir.x),
76                                            fabsf(dir.y),
77                                            fabsf(dir.z)));
78         int kx = kz + 1; if(kx == 3) kx = 0;
79         int ky = kx + 1; if(ky == 3) ky = 0;
80
81         /* Swap kx and ky dimensions to preserve winding direction of triangles. */
82         if(IDX(dir, kz) < 0.0f) {
83                 int tmp = kx;
84                 kx = ky;
85                 ky = tmp;
86         }
87
88         /* Calculate the shear constants. */
89         float inv_dir_z = 1.0f / IDX(dir, kz);
90         isect_precalc->Sx = IDX(dir, kx) * inv_dir_z;
91         isect_precalc->Sy = IDX(dir, ky) * inv_dir_z;
92         isect_precalc->Sz = inv_dir_z;
93
94         /* Store the dimensions. */
95         isect_precalc->kx = kx;
96         isect_precalc->ky = ky;
97         isect_precalc->kz = kz;
98 }
99
100 /* TODO(sergey): Make it general utility function. */
101 ccl_device_inline float xor_signmast(float x, int y)
102 {
103         return __int_as_float(__float_as_int(x) ^ y);
104 }
105
106 ccl_device_inline bool triangle_intersect(KernelGlobals *kg,
107                                           const IsectPrecalc *isect_precalc,
108                                           Intersection *isect,
109                                           float3 P,
110                                           uint visibility,
111                                           int object,
112                                           int triAddr)
113 {
114         const int kx = isect_precalc->kx;
115         const int ky = isect_precalc->ky;
116         const int kz = isect_precalc->kz;
117         const float Sx = isect_precalc->Sx;
118         const float Sy = isect_precalc->Sy;
119         const float Sz = isect_precalc->Sz;
120
121         /* Calculate vertices relative to ray origin. */
122         const float4 tri_a = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+0),
123                      tri_b = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+1),
124                      tri_c = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+2);
125         const float3 A = make_float3(tri_a.x - P.x, tri_a.y - P.y, tri_a.z - P.z);
126         const float3 B = make_float3(tri_b.x - P.x, tri_b.y - P.y, tri_b.z - P.z);
127         const float3 C = make_float3(tri_c.x - P.x, tri_c.y - P.y, tri_c.z - P.z);
128
129         const float A_kx = IDX(A, kx), A_ky = IDX(A, ky), A_kz = IDX(A, kz);
130         const float B_kx = IDX(B, kx), B_ky = IDX(B, ky), B_kz = IDX(B, kz);
131         const float C_kx = IDX(C, kx), C_ky = IDX(C, ky), C_kz = IDX(C, kz);
132
133         /* Perform shear and scale of vertices. */
134         const float Ax = A_kx - Sx * A_kz;
135         const float Ay = A_ky - Sy * A_kz;
136         const float Bx = B_kx - Sx * B_kz;
137         const float By = B_ky - Sy * B_kz;
138         const float Cx = C_kx - Sx * C_kz;
139         const float Cy = C_ky - Sy * C_kz;
140
141         /* Calculate scaled barycentric coordinates. */
142         float U = Cx * By - Cy * Bx;
143         int sign_mask = (__float_as_int(U) & 0x80000000);
144         float V = Ax * Cy - Ay * Cx;
145         if(sign_mask != (__float_as_int(V) & 0x80000000)) {
146                 return false;
147         }
148         float W = Bx * Ay - By * Ax;
149         if(sign_mask != (__float_as_int(W) & 0x80000000)) {
150                 return false;
151         }
152
153         /* Calculate determinant. */
154         float det = U + V + W;
155         if(UNLIKELY(det == 0.0f)) {
156                 return false;
157         }
158
159         /* Calculate scaled z−coordinates of vertices and use them to calculate
160          * the hit distance.
161          */
162         const float T = (U * A_kz + V * B_kz + W * C_kz) * Sz;
163         const float sign_T = xor_signmast(T, sign_mask);
164         if((sign_T < 0.0f) ||
165            (sign_T > isect->t * xor_signmast(det, sign_mask)))
166         {
167                 return false;
168         }
169
170 #ifdef __VISIBILITY_FLAG__
171         /* visibility flag test. we do it here under the assumption
172          * that most triangles are culled by node flags */
173         if(kernel_tex_fetch(__prim_visibility, triAddr) & visibility)
174 #endif
175         {
176                 /* Normalize U, V, W, and T. */
177                 const float inv_det = 1.0f / det;
178                 isect->prim = triAddr;
179                 isect->object = object;
180                 isect->type = PRIMITIVE_TRIANGLE;
181                 isect->u = U * inv_det;
182                 isect->v = V * inv_det;
183                 isect->t = T * inv_det;
184                 return true;
185         }
186         return false;
187 }
188
189 /* Special ray intersection routines for subsurface scattering. In that case we
190  * only want to intersect with primitives in the same object, and if case of
191  * multiple hits we pick a single random primitive as the intersection point.
192  */
193
194 #ifdef __SUBSURFACE__
195 ccl_device_inline void triangle_intersect_subsurface(
196         KernelGlobals *kg,
197         const IsectPrecalc *isect_precalc,
198         Intersection *isect_array,
199         float3 P,
200         int object,
201         int triAddr,
202         float tmax,
203         uint *num_hits,
204         uint *lcg_state,
205         int max_hits)
206 {
207         const int kx = isect_precalc->kx;
208         const int ky = isect_precalc->ky;
209         const int kz = isect_precalc->kz;
210         const float Sx = isect_precalc->Sx;
211         const float Sy = isect_precalc->Sy;
212         const float Sz = isect_precalc->Sz;
213
214         /* Calculate vertices relative to ray origin. */
215         const float4 tri_a = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+0),
216                      tri_b = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+1),
217                      tri_c = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+2);
218         const float3 A = make_float3(tri_a.x - P.x, tri_a.y - P.y, tri_a.z - P.z);
219         const float3 B = make_float3(tri_b.x - P.x, tri_b.y - P.y, tri_b.z - P.z);
220         const float3 C = make_float3(tri_c.x - P.x, tri_c.y - P.y, tri_c.z - P.z);
221
222         const float A_kx = IDX(A, kx), A_ky = IDX(A, ky), A_kz = IDX(A, kz);
223         const float B_kx = IDX(B, kx), B_ky = IDX(B, ky), B_kz = IDX(B, kz);
224         const float C_kx = IDX(C, kx), C_ky = IDX(C, ky), C_kz = IDX(C, kz);
225
226         /* Perform shear and scale of vertices. */
227         const float Ax = A_kx - Sx * A_kz;
228         const float Ay = A_ky - Sy * A_kz;
229         const float Bx = B_kx - Sx * B_kz;
230         const float By = B_ky - Sy * B_kz;
231         const float Cx = C_kx - Sx * C_kz;
232         const float Cy = C_ky - Sy * C_kz;
233
234         /* Calculate scaled barycentric coordinates. */
235         float U = Cx * By - Cy * Bx;
236         int sign_mask = (__float_as_int(U) & 0x80000000);
237         float V = Ax * Cy - Ay * Cx;
238         if(sign_mask != (__float_as_int(V) & 0x80000000)) {
239                 return;
240         }
241         float W = Bx * Ay - By * Ax;
242         if(sign_mask != (__float_as_int(W) & 0x80000000)) {
243                 return;
244         }
245
246         /* Calculate determinant. */
247         float det = U + V + W;
248         if(UNLIKELY(det == 0.0f)) {
249                 return;
250         }
251
252         /* Calculate scaled z−coordinates of vertices and use them to calculate
253          * the hit distance.
254          */
255         const float T = (U * A_kz + V * B_kz + W * C_kz) * Sz;
256         const float sign_T = xor_signmast(T, sign_mask);
257         if((sign_T < 0.0f) ||
258            (sign_T > tmax * xor_signmast(det, sign_mask)))
259         {
260                 return;
261         }
262
263         /* Normalize U, V, W, and T. */
264         const float inv_det = 1.0f / det;
265
266         (*num_hits)++;
267         int hit;
268
269         if(*num_hits <= max_hits) {
270                 hit = *num_hits - 1;
271         }
272         else {
273                 /* reservoir sampling: if we are at the maximum number of
274                  * hits, randomly replace element or skip it */
275                 hit = lcg_step_uint(lcg_state) % *num_hits;
276
277                 if(hit >= max_hits)
278                         return;
279         }
280
281         /* record intersection */
282         Intersection *isect = &isect_array[hit];
283         isect->prim = triAddr;
284         isect->object = object;
285         isect->type = PRIMITIVE_TRIANGLE;
286         isect->u = U * inv_det;
287         isect->v = V * inv_det;
288         isect->t = T * inv_det;
289 }
290 #endif
291
292 /* Refine triangle intersection to more precise hit point. For rays that travel
293  * far the precision is often not so good, this reintersects the primitive from
294  * a closer distance. */
295
296 /* Reintersections uses the paper:
297  *
298  * Tomas Moeller
299  * Fast, minimum storage ray/triangle intersection
300  * http://www.cs.virginia.edu/~gfx/Courses/2003/ImageSynthesis/papers/Acceleration/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf
301  */
302
303 ccl_device_inline float3 triangle_refine(KernelGlobals *kg,
304                                          ShaderData *sd,
305                                          const Intersection *isect,
306                                          const Ray *ray)
307 {
308         float3 P = ray->P;
309         float3 D = ray->D;
310         float t = isect->t;
311
312 #ifdef __INTERSECTION_REFINE__
313         if(isect->object != OBJECT_NONE) {
314                 if(UNLIKELY(t == 0.0f)) {
315                         return P;
316                 }
317 #ifdef __OBJECT_MOTION__
318                 Transform tfm = ccl_fetch(sd, ob_itfm);
319 #else
320                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
321 #endif
322
323                 P = transform_point(&tfm, P);
324                 D = transform_direction(&tfm, D*t);
325                 D = normalize_len(D, &t);
326         }
327
328         P = P + D*t;
329
330         const float4 tri_a = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+0),
331                      tri_b = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+1),
332                      tri_c = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+2);
333         float3 edge1 = make_float3(tri_a.x - tri_c.x, tri_a.y - tri_c.y, tri_a.z - tri_c.z);
334         float3 edge2 = make_float3(tri_b.x - tri_c.x, tri_b.y - tri_c.y, tri_b.z - tri_c.z);
335         float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
336         float3 qvec = cross(tvec, edge1);
337         float3 pvec = cross(D, edge2);
338         float rt = dot(edge2, qvec) / dot(edge1, pvec);
339
340         P = P + D*rt;
341
342         if(isect->object != OBJECT_NONE) {
343 #ifdef __OBJECT_MOTION__
344                 Transform tfm = ccl_fetch(sd, ob_tfm);
345 #else
346                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
347 #endif
348
349                 P = transform_point(&tfm, P);
350         }
351
352         return P;
353 #else
354         return P + D*t;
355 #endif
356 }
357
358 /* Same as above, except that isect->t is assumed to be in object space for
359  * instancing.
360  */
361 ccl_device_inline float3 triangle_refine_subsurface(KernelGlobals *kg,
362                                                     ShaderData *sd,
363                                                     const Intersection *isect,
364                                                     const Ray *ray)
365 {
366         float3 P = ray->P;
367         float3 D = ray->D;
368         float t = isect->t;
369
370 #ifdef __INTERSECTION_REFINE__
371         if(isect->object != OBJECT_NONE) {
372 #ifdef __OBJECT_MOTION__
373                 Transform tfm = ccl_fetch(sd, ob_itfm);
374 #else
375                 Transform tfm = object_fetch_transform(kg,
376                                                        isect->object,
377                                                        OBJECT_INVERSE_TRANSFORM);
378 #endif
379
380                 P = transform_point(&tfm, P);
381                 D = transform_direction(&tfm, D);
382                 D = normalize(D);
383         }
384
385         P = P + D*t;
386
387         const float4 tri_a = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+0),
388                      tri_b = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+1),
389                      tri_c = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+2);
390         float3 edge1 = make_float3(tri_a.x - tri_c.x, tri_a.y - tri_c.y, tri_a.z - tri_c.z);
391         float3 edge2 = make_float3(tri_b.x - tri_c.x, tri_b.y - tri_c.y, tri_b.z - tri_c.z);
392         float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
393         float3 qvec = cross(tvec, edge1);
394         float3 pvec = cross(D, edge2);
395         float rt = dot(edge2, qvec) / dot(edge1, pvec);
396
397         P = P + D*rt;
398
399         if(isect->object != OBJECT_NONE) {
400 #ifdef __OBJECT_MOTION__
401                 Transform tfm = ccl_fetch(sd, ob_tfm);
402 #else
403                 Transform tfm = object_fetch_transform(kg,
404                                                        isect->object,
405                                                        OBJECT_TRANSFORM);
406 #endif
407
408                 P = transform_point(&tfm, P);
409         }
410
411         return P;
412 #else
413         return P + D*t;
414 #endif
415 }
416
417 #undef IDX
418
419 CCL_NAMESPACE_END