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