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 
triangle_intersect(KernelGlobals * kg,Intersection * isect,float3 P,float3 dir,uint visibility,int object,int prim_addr)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 /* Special ray intersection routines for subsurface scattering. In that case we
75  * only want to intersect with primitives in the same object, and if case of
76  * multiple hits we pick a single random primitive as the intersection point.
77  * Returns whether traversal should be stopped.
78  */
79 
80 #ifdef __BVH_LOCAL__
triangle_intersect_local(KernelGlobals * kg,LocalIntersection * local_isect,float3 P,float3 dir,int object,int local_object,int prim_addr,float tmax,uint * lcg_state,int max_hits)81 ccl_device_inline bool triangle_intersect_local(KernelGlobals *kg,
82                                                 LocalIntersection *local_isect,
83                                                 float3 P,
84                                                 float3 dir,
85                                                 int object,
86                                                 int local_object,
87                                                 int prim_addr,
88                                                 float tmax,
89                                                 uint *lcg_state,
90                                                 int max_hits)
91 {
92   /* Only intersect with matching object, for instanced objects we
93    * already know we are only intersecting the right object. */
94   if (object == OBJECT_NONE) {
95     if (kernel_tex_fetch(__prim_object, prim_addr) != local_object) {
96       return false;
97     }
98   }
99 
100   const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, prim_addr);
101 #  if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
102   const ssef *ssef_verts = (ssef *)&kg->__prim_tri_verts.data[tri_vindex];
103 #  else
104   const float3 tri_a = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex + 0)),
105                tri_b = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex + 1)),
106                tri_c = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex + 2));
107 #  endif
108   float t, u, v;
109   if (!ray_triangle_intersect(P,
110                               dir,
111                               tmax,
112 #  if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
113                               ssef_verts,
114 #  else
115                               tri_a,
116                               tri_b,
117                               tri_c,
118 #  endif
119                               &u,
120                               &v,
121                               &t)) {
122     return false;
123   }
124 
125   /* If no actual hit information is requested, just return here. */
126   if (max_hits == 0) {
127     return true;
128   }
129 
130   int hit;
131   if (lcg_state) {
132     /* Record up to max_hits intersections. */
133     for (int i = min(max_hits, local_isect->num_hits) - 1; i >= 0; --i) {
134       if (local_isect->hits[i].t == t) {
135         return false;
136       }
137     }
138 
139     local_isect->num_hits++;
140 
141     if (local_isect->num_hits <= max_hits) {
142       hit = local_isect->num_hits - 1;
143     }
144     else {
145       /* reservoir sampling: if we are at the maximum number of
146        * hits, randomly replace element or skip it */
147       hit = lcg_step_uint(lcg_state) % local_isect->num_hits;
148 
149       if (hit >= max_hits)
150         return false;
151     }
152   }
153   else {
154     /* Record closest intersection only. */
155     if (local_isect->num_hits && t > local_isect->hits[0].t) {
156       return false;
157     }
158 
159     hit = 0;
160     local_isect->num_hits = 1;
161   }
162 
163   /* Record intersection. */
164   Intersection *isect = &local_isect->hits[hit];
165   isect->prim = prim_addr;
166   isect->object = object;
167   isect->type = PRIMITIVE_TRIANGLE;
168   isect->u = u;
169   isect->v = v;
170   isect->t = t;
171 
172   /* Record geometric normal. */
173 #  if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
174   const float3 tri_a = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex + 0)),
175                tri_b = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex + 1)),
176                tri_c = float4_to_float3(kernel_tex_fetch(__prim_tri_verts, tri_vindex + 2));
177 #  endif
178   local_isect->Ng[hit] = normalize(cross(tri_b - tri_a, tri_c - tri_a));
179 
180   return false;
181 }
182 #endif /* __BVH_LOCAL__ */
183 
184 /* Refine triangle intersection to more precise hit point. For rays that travel
185  * far the precision is often not so good, this reintersects the primitive from
186  * a closer distance. */
187 
188 /* Reintersections uses the paper:
189  *
190  * Tomas Moeller
191  * Fast, minimum storage ray/triangle intersection
192  * http://www.cs.virginia.edu/~gfx/Courses/2003/ImageSynthesis/papers/Acceleration/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf
193  */
194 
triangle_refine(KernelGlobals * kg,ShaderData * sd,const Intersection * isect,const Ray * ray)195 ccl_device_inline float3 triangle_refine(KernelGlobals *kg,
196                                          ShaderData *sd,
197                                          const Intersection *isect,
198                                          const Ray *ray)
199 {
200   float3 P = ray->P;
201   float3 D = ray->D;
202   float t = isect->t;
203 
204 #ifdef __INTERSECTION_REFINE__
205   if (isect->object != OBJECT_NONE) {
206     if (UNLIKELY(t == 0.0f)) {
207       return P;
208     }
209 #  ifdef __OBJECT_MOTION__
210     Transform tfm = sd->ob_itfm;
211 #  else
212     Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
213 #  endif
214 
215     P = transform_point(&tfm, P);
216     D = transform_direction(&tfm, D * t);
217     D = normalize_len(D, &t);
218   }
219 
220   P = P + D * t;
221 
222   const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, isect->prim);
223   const float4 tri_a = kernel_tex_fetch(__prim_tri_verts, tri_vindex + 0),
224                tri_b = kernel_tex_fetch(__prim_tri_verts, tri_vindex + 1),
225                tri_c = kernel_tex_fetch(__prim_tri_verts, tri_vindex + 2);
226   float3 edge1 = make_float3(tri_a.x - tri_c.x, tri_a.y - tri_c.y, tri_a.z - tri_c.z);
227   float3 edge2 = make_float3(tri_b.x - tri_c.x, tri_b.y - tri_c.y, tri_b.z - tri_c.z);
228   float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
229   float3 qvec = cross(tvec, edge1);
230   float3 pvec = cross(D, edge2);
231   float det = dot(edge1, pvec);
232   if (det != 0.0f) {
233     /* If determinant is zero it means ray lies in the plane of
234      * the triangle. It is possible in theory due to watertight
235      * nature of triangle intersection. For such cases we simply
236      * don't refine intersection hoping it'll go all fine.
237      */
238     float rt = dot(edge2, qvec) / det;
239     P = P + D * rt;
240   }
241 
242   if (isect->object != OBJECT_NONE) {
243 #  ifdef __OBJECT_MOTION__
244     Transform tfm = sd->ob_tfm;
245 #  else
246     Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
247 #  endif
248 
249     P = transform_point(&tfm, P);
250   }
251 
252   return P;
253 #else
254   return P + D * t;
255 #endif
256 }
257 
258 /* Same as above, except that isect->t is assumed to be in object space for
259  * instancing.
260  */
triangle_refine_local(KernelGlobals * kg,ShaderData * sd,const Intersection * isect,const Ray * ray)261 ccl_device_inline float3 triangle_refine_local(KernelGlobals *kg,
262                                                ShaderData *sd,
263                                                const Intersection *isect,
264                                                const Ray *ray)
265 {
266 #ifdef __KERNEL_OPTIX__
267   /* isect->t is always in world space with OptiX. */
268   return triangle_refine(kg, sd, isect, ray);
269 #else
270   float3 P = ray->P;
271   float3 D = ray->D;
272   float t = isect->t;
273 
274   if (isect->object != OBJECT_NONE) {
275 #  ifdef __OBJECT_MOTION__
276     Transform tfm = sd->ob_itfm;
277 #  else
278     Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
279 #  endif
280 
281     P = transform_point(&tfm, P);
282     D = transform_direction(&tfm, D);
283     D = normalize(D);
284   }
285 
286   P = P + D * t;
287 
288 #  ifdef __INTERSECTION_REFINE__
289   const uint tri_vindex = kernel_tex_fetch(__prim_tri_index, isect->prim);
290   const float4 tri_a = kernel_tex_fetch(__prim_tri_verts, tri_vindex + 0),
291                tri_b = kernel_tex_fetch(__prim_tri_verts, tri_vindex + 1),
292                tri_c = kernel_tex_fetch(__prim_tri_verts, tri_vindex + 2);
293   float3 edge1 = make_float3(tri_a.x - tri_c.x, tri_a.y - tri_c.y, tri_a.z - tri_c.z);
294   float3 edge2 = make_float3(tri_b.x - tri_c.x, tri_b.y - tri_c.y, tri_b.z - tri_c.z);
295   float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
296   float3 qvec = cross(tvec, edge1);
297   float3 pvec = cross(D, edge2);
298   float det = dot(edge1, pvec);
299   if (det != 0.0f) {
300     /* If determinant is zero it means ray lies in the plane of
301      * the triangle. It is possible in theory due to watertight
302      * nature of triangle intersection. For such cases we simply
303      * don't refine intersection hoping it'll go all fine.
304      */
305     float rt = dot(edge2, qvec) / det;
306     P = P + D * rt;
307   }
308 #  endif /* __INTERSECTION_REFINE__ */
309 
310   if (isect->object != OBJECT_NONE) {
311 #  ifdef __OBJECT_MOTION__
312     Transform tfm = sd->ob_tfm;
313 #  else
314     Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
315 #  endif
316 
317     P = transform_point(&tfm, P);
318   }
319 
320   return P;
321 #endif
322 }
323 
324 CCL_NAMESPACE_END
325