Cysles: Avoid having ShaderData on the stack
[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_signmask(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         float V = Ax * Cy - Ay * Cx;
144         float W = Bx * Ay - By * Ax;
145         if((U < 0.0f || V < 0.0f || W < 0.0f) &&
146            (U > 0.0f || V > 0.0f || W > 0.0f))
147         {
148                 return false;
149         }
150
151         /* Calculate determinant. */
152         float det = U + V + W;
153         if(UNLIKELY(det == 0.0f)) {
154                 return false;
155         }
156
157         /* Calculate scaled z-coordinates of vertices and use them to calculate
158          * the hit distance.
159          */
160         const float T = (U * A_kz + V * B_kz + W * C_kz) * Sz;
161         const int sign_det = (__float_as_int(det) & 0x80000000);
162         const float sign_T = xor_signmask(T, sign_det);
163         if((sign_T < 0.0f) ||
164            (sign_T > isect->t * xor_signmask(det, sign_det)))
165         {
166                 return false;
167         }
168
169 #ifdef __VISIBILITY_FLAG__
170         /* visibility flag test. we do it here under the assumption
171          * that most triangles are culled by node flags */
172         if(kernel_tex_fetch(__prim_visibility, triAddr) & visibility)
173 #endif
174         {
175 #ifdef __KERNEL_GPU__
176                 float4 a = tri_b - tri_a, b = tri_c - tri_a;
177                 if(len_squared(make_float3(a.y*b.z - a.z*b.y,
178                                            a.z*b.x - a.x*b.z,
179                                            a.x*b.y - a.y*b.x)) == 0.0f)
180                 {
181                         return false;
182                 }
183 #endif
184
185                 /* Normalize U, V, W, and T. */
186                 const float inv_det = 1.0f / det;
187                 isect->prim = triAddr;
188                 isect->object = object;
189                 isect->type = PRIMITIVE_TRIANGLE;
190                 isect->u = U * inv_det;
191                 isect->v = V * inv_det;
192                 isect->t = T * inv_det;
193                 return true;
194         }
195         return false;
196 }
197
198 /* Special ray intersection routines for subsurface scattering. In that case we
199  * only want to intersect with primitives in the same object, and if case of
200  * multiple hits we pick a single random primitive as the intersection point.
201  */
202
203 #ifdef __SUBSURFACE__
204 ccl_device_inline void triangle_intersect_subsurface(
205         KernelGlobals *kg,
206         const IsectPrecalc *isect_precalc,
207         SubsurfaceIntersection *ss_isect,
208         float3 P,
209         int object,
210         int triAddr,
211         float tmax,
212         uint *lcg_state,
213         int max_hits)
214 {
215         const int kx = isect_precalc->kx;
216         const int ky = isect_precalc->ky;
217         const int kz = isect_precalc->kz;
218         const float Sx = isect_precalc->Sx;
219         const float Sy = isect_precalc->Sy;
220         const float Sz = isect_precalc->Sz;
221
222         /* Calculate vertices relative to ray origin. */
223         const float4 tri_a = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+0),
224                      tri_b = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+1),
225                      tri_c = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+2);
226         const float3 A = make_float3(tri_a.x - P.x, tri_a.y - P.y, tri_a.z - P.z);
227         const float3 B = make_float3(tri_b.x - P.x, tri_b.y - P.y, tri_b.z - P.z);
228         const float3 C = make_float3(tri_c.x - P.x, tri_c.y - P.y, tri_c.z - P.z);
229
230         const float A_kx = IDX(A, kx), A_ky = IDX(A, ky), A_kz = IDX(A, kz);
231         const float B_kx = IDX(B, kx), B_ky = IDX(B, ky), B_kz = IDX(B, kz);
232         const float C_kx = IDX(C, kx), C_ky = IDX(C, ky), C_kz = IDX(C, kz);
233
234         /* Perform shear and scale of vertices. */
235         const float Ax = A_kx - Sx * A_kz;
236         const float Ay = A_ky - Sy * A_kz;
237         const float Bx = B_kx - Sx * B_kz;
238         const float By = B_ky - Sy * B_kz;
239         const float Cx = C_kx - Sx * C_kz;
240         const float Cy = C_ky - Sy * C_kz;
241
242         /* Calculate scaled barycentric coordinates. */
243         float U = Cx * By - Cy * Bx;
244         float V = Ax * Cy - Ay * Cx;
245         float W = Bx * Ay - By * Ax;
246
247         if((U < 0.0f || V < 0.0f || W < 0.0f) &&
248            (U > 0.0f || V > 0.0f || W > 0.0f))
249         {
250                 return;
251         }
252
253         /* Calculate determinant. */
254         float det = U + V + W;
255         if(UNLIKELY(det == 0.0f)) {
256                 return;
257         }
258
259         /* Calculate scaled z−coordinates of vertices and use them to calculate
260          * the hit distance.
261          */
262         const int sign_det = (__float_as_int(det) & 0x80000000);
263         const float T = (U * A_kz + V * B_kz + W * C_kz) * Sz;
264         const float sign_T = xor_signmask(T, sign_det);
265         if((sign_T < 0.0f) ||
266            (sign_T > tmax * xor_signmask(det, sign_det)))
267         {
268                 return;
269         }
270
271         /* Normalize U, V, W, and T. */
272         const float inv_det = 1.0f / det;
273
274         ss_isect->num_hits++;
275         int hit;
276
277         if(ss_isect->num_hits <= max_hits) {
278                 hit = ss_isect->num_hits - 1;
279         }
280         else {
281                 /* reservoir sampling: if we are at the maximum number of
282                  * hits, randomly replace element or skip it */
283                 hit = lcg_step_uint(lcg_state) % ss_isect->num_hits;
284
285                 if(hit >= max_hits)
286                         return;
287         }
288
289         /* record intersection */
290         Intersection *isect = &ss_isect->hits[hit];
291         isect->prim = triAddr;
292         isect->object = object;
293         isect->type = PRIMITIVE_TRIANGLE;
294         isect->u = U * inv_det;
295         isect->v = V * inv_det;
296         isect->t = T * inv_det;
297
298         /* Record geometric normal. */
299         /* TODO(sergey): Use float4_to_float3() on just an edges. */
300         const float3 v0 = float4_to_float3(tri_a);
301         const float3 v1 = float4_to_float3(tri_b);
302         const float3 v2 = float4_to_float3(tri_c);
303         ss_isect->Ng[hit] = normalize(cross(v1 - v0, v2 - v0));
304 }
305 #endif
306
307 /* Refine triangle intersection to more precise hit point. For rays that travel
308  * far the precision is often not so good, this reintersects the primitive from
309  * a closer distance. */
310
311 /* Reintersections uses the paper:
312  *
313  * Tomas Moeller
314  * Fast, minimum storage ray/triangle intersection
315  * http://www.cs.virginia.edu/~gfx/Courses/2003/ImageSynthesis/papers/Acceleration/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf
316  */
317
318 ccl_device_inline float3 triangle_refine(KernelGlobals *kg,
319                                          ShaderData *sd,
320                                          const Intersection *isect,
321                                          const Ray *ray)
322 {
323         float3 P = ray->P;
324         float3 D = ray->D;
325         float t = isect->t;
326
327 #ifdef __INTERSECTION_REFINE__
328         if(isect->object != OBJECT_NONE) {
329                 if(UNLIKELY(t == 0.0f)) {
330                         return P;
331                 }
332 #ifdef __OBJECT_MOTION__
333                 Transform tfm = ccl_fetch(sd, ob_itfm);
334 #else
335                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
336 #endif
337
338                 P = transform_point(&tfm, P);
339                 D = transform_direction(&tfm, D*t);
340                 D = normalize_len(D, &t);
341         }
342
343         P = P + D*t;
344
345         const float4 tri_a = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+0),
346                      tri_b = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+1),
347                      tri_c = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+2);
348         float3 edge1 = make_float3(tri_a.x - tri_c.x, tri_a.y - tri_c.y, tri_a.z - tri_c.z);
349         float3 edge2 = make_float3(tri_b.x - tri_c.x, tri_b.y - tri_c.y, tri_b.z - tri_c.z);
350         float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
351         float3 qvec = cross(tvec, edge1);
352         float3 pvec = cross(D, edge2);
353         float rt = dot(edge2, qvec) / dot(edge1, pvec);
354
355         P = P + D*rt;
356
357         if(isect->object != OBJECT_NONE) {
358 #ifdef __OBJECT_MOTION__
359                 Transform tfm = ccl_fetch(sd, ob_tfm);
360 #else
361                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
362 #endif
363
364                 P = transform_point(&tfm, P);
365         }
366
367         return P;
368 #else
369         return P + D*t;
370 #endif
371 }
372
373 /* Same as above, except that isect->t is assumed to be in object space for
374  * instancing.
375  */
376 ccl_device_inline float3 triangle_refine_subsurface(KernelGlobals *kg,
377                                                     ShaderData *sd,
378                                                     const Intersection *isect,
379                                                     const Ray *ray)
380 {
381         float3 P = ray->P;
382         float3 D = ray->D;
383         float t = isect->t;
384
385 #ifdef __INTERSECTION_REFINE__
386         if(isect->object != OBJECT_NONE) {
387 #ifdef __OBJECT_MOTION__
388                 Transform tfm = ccl_fetch(sd, ob_itfm);
389 #else
390                 Transform tfm = object_fetch_transform(kg,
391                                                        isect->object,
392                                                        OBJECT_INVERSE_TRANSFORM);
393 #endif
394
395                 P = transform_point(&tfm, P);
396                 D = transform_direction(&tfm, D);
397                 D = normalize(D);
398         }
399
400         P = P + D*t;
401
402         const float4 tri_a = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+0),
403                      tri_b = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+1),
404                      tri_c = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+2);
405         float3 edge1 = make_float3(tri_a.x - tri_c.x, tri_a.y - tri_c.y, tri_a.z - tri_c.z);
406         float3 edge2 = make_float3(tri_b.x - tri_c.x, tri_b.y - tri_c.y, tri_b.z - tri_c.z);
407         float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
408         float3 qvec = cross(tvec, edge1);
409         float3 pvec = cross(D, edge2);
410         float rt = dot(edge2, qvec) / dot(edge1, pvec);
411
412         P = P + D*rt;
413
414         if(isect->object != OBJECT_NONE) {
415 #ifdef __OBJECT_MOTION__
416                 Transform tfm = ccl_fetch(sd, ob_tfm);
417 #else
418                 Transform tfm = object_fetch_transform(kg,
419                                                        isect->object,
420                                                        OBJECT_TRANSFORM);
421 #endif
422
423                 P = transform_point(&tfm, P);
424         }
425
426         return P;
427 #else
428         return P + D*t;
429 #endif
430 }
431
432 #undef IDX
433
434 CCL_NAMESPACE_END