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