add BLI_strcpy_rlen, replace strcat, which was used in misleading way.
[blender.git] / intern / cycles / kernel / kernel_bvh.h
1 /*
2  * Adapted from code Copyright 2009-2010 NVIDIA Corporation
3  * Modifications Copyright 2011, Blender Foundation.
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  */
17
18 CCL_NAMESPACE_BEGIN
19
20 /*
21  * "Persistent while-while kernel" used in:
22  *
23  * "Understanding the Efficiency of Ray Traversal on GPUs",
24  * Timo Aila and Samuli Laine,
25  * Proc. High-Performance Graphics 2009
26  */
27
28 /* bottom-most stack entry, indicating the end of traversal */
29
30 #define ENTRYPOINT_SENTINEL 0x76543210
31 /* 64 object BVH + 64 mesh BVH + 64 object node splitting */
32 #define BVH_STACK_SIZE 192
33 #define BVH_NODE_SIZE 4
34 #define TRI_NODE_SIZE 3
35
36 /* silly workaround for float extended precision that happens when compiling
37  * without sse support on x86, it results in different results for float ops
38  * that you would otherwise expect to compare correctly */
39 #if !defined(__i386__) || defined(__SSE__)
40 #define NO_EXTENDED_PRECISION
41 #else
42 #define NO_EXTENDED_PRECISION volatile
43 #endif
44
45 __device_inline float3 bvh_inverse_direction(float3 dir)
46 {
47         /* avoid divide by zero (ooeps = exp2f(-80.0f)) */
48         float ooeps = 0.00000000000000000000000082718061255302767487140869206996285356581211090087890625f;
49         float3 idir;
50
51         idir.x = 1.0f/((fabsf(dir.x) > ooeps)? dir.x: copysignf(ooeps, dir.x));
52         idir.y = 1.0f/((fabsf(dir.y) > ooeps)? dir.y: copysignf(ooeps, dir.y));
53         idir.z = 1.0f/((fabsf(dir.z) > ooeps)? dir.z: copysignf(ooeps, dir.z));
54
55         return idir;
56 }
57
58 __device_inline void bvh_instance_push(KernelGlobals *kg, int object, const Ray *ray, float3 *P, float3 *idir, float *t, const float tmax)
59 {
60         Transform tfm = object_fetch_transform(kg, object, OBJECT_INVERSE_TRANSFORM);
61
62         *P = transform_point(&tfm, ray->P);
63
64         float3 dir = transform_direction(&tfm, ray->D);
65
66         float len;
67         dir = normalize_len(dir, &len);
68
69         *idir = bvh_inverse_direction(dir);
70
71         if(*t != FLT_MAX)
72                 *t *= len;
73 }
74
75 __device_inline void bvh_instance_pop(KernelGlobals *kg, int object, const Ray *ray, float3 *P, float3 *idir, float *t, const float tmax)
76 {
77         if(*t != FLT_MAX) {
78                 Transform tfm = object_fetch_transform(kg, object, OBJECT_TRANSFORM);
79                 *t *= len(transform_direction(&tfm, 1.0f/(*idir)));
80         }
81
82         *P = ray->P;
83         *idir = bvh_inverse_direction(ray->D);
84 }
85
86 #ifdef __OBJECT_MOTION__
87 __device_inline void bvh_instance_motion_push(KernelGlobals *kg, int object, const Ray *ray, float3 *P, float3 *idir, float *t, Transform *tfm, const float tmax)
88 {
89         Transform itfm;
90         *tfm = object_fetch_transform_motion_test(kg, object, ray->time, &itfm);
91
92         *P = transform_point(&itfm, ray->P);
93
94         float3 dir = transform_direction(&itfm, ray->D);
95
96         float len;
97         dir = normalize_len(dir, &len);
98
99         *idir = bvh_inverse_direction(dir);
100
101         if(*t != FLT_MAX)
102                 *t *= len;
103 }
104
105 __device_inline void bvh_instance_motion_pop(KernelGlobals *kg, int object, const Ray *ray, float3 *P, float3 *idir, float *t, Transform *tfm, const float tmax)
106 {
107         if(*t != FLT_MAX)
108                 *t *= len(transform_direction(tfm, 1.0f/(*idir)));
109
110         *P = ray->P;
111         *idir = bvh_inverse_direction(ray->D);
112 }
113 #endif
114
115 /* intersect two bounding boxes */
116 #ifdef __HAIR__
117 __device_inline void bvh_node_intersect(KernelGlobals *kg,
118         bool *traverseChild0, bool *traverseChild1,
119         bool *closestChild1, int *nodeAddr0, int *nodeAddr1,
120         float3 P, float3 idir, float t, uint visibility, int nodeAddr, float difl, float extmax)
121 {
122 #else
123 __device_inline void bvh_node_intersect(KernelGlobals *kg,
124         bool *traverseChild0, bool *traverseChild1,
125         bool *closestChild1, int *nodeAddr0, int *nodeAddr1,
126         float3 P, float3 idir, float t, uint visibility, int nodeAddr)
127 {
128 #endif
129
130         /* fetch node data */
131         float4 n0xy = kernel_tex_fetch(__bvh_nodes, nodeAddr*BVH_NODE_SIZE+0);
132         float4 n1xy = kernel_tex_fetch(__bvh_nodes, nodeAddr*BVH_NODE_SIZE+1);
133         float4 nz = kernel_tex_fetch(__bvh_nodes, nodeAddr*BVH_NODE_SIZE+2);
134         float4 cnodes = kernel_tex_fetch(__bvh_nodes, nodeAddr*BVH_NODE_SIZE+3);
135
136         /* intersect ray against child nodes */
137         float3 ood = P * idir;
138         NO_EXTENDED_PRECISION float c0lox = n0xy.x * idir.x - ood.x;
139         NO_EXTENDED_PRECISION float c0hix = n0xy.y * idir.x - ood.x;
140         NO_EXTENDED_PRECISION float c0loy = n0xy.z * idir.y - ood.y;
141         NO_EXTENDED_PRECISION float c0hiy = n0xy.w * idir.y - ood.y;
142         NO_EXTENDED_PRECISION float c0loz = nz.x * idir.z - ood.z;
143         NO_EXTENDED_PRECISION float c0hiz = nz.y * idir.z - ood.z;
144         NO_EXTENDED_PRECISION float c0min = max4(min(c0lox, c0hix), min(c0loy, c0hiy), min(c0loz, c0hiz), 0.0f);
145         NO_EXTENDED_PRECISION float c0max = min4(max(c0lox, c0hix), max(c0loy, c0hiy), max(c0loz, c0hiz), t);
146
147         NO_EXTENDED_PRECISION float c1loz = nz.z * idir.z - ood.z;
148         NO_EXTENDED_PRECISION float c1hiz = nz.w * idir.z - ood.z;
149         NO_EXTENDED_PRECISION float c1lox = n1xy.x * idir.x - ood.x;
150         NO_EXTENDED_PRECISION float c1hix = n1xy.y * idir.x - ood.x;
151         NO_EXTENDED_PRECISION float c1loy = n1xy.z * idir.y - ood.y;
152         NO_EXTENDED_PRECISION float c1hiy = n1xy.w * idir.y - ood.y;
153         NO_EXTENDED_PRECISION float c1min = max4(min(c1lox, c1hix), min(c1loy, c1hiy), min(c1loz, c1hiz), 0.0f);
154         NO_EXTENDED_PRECISION float c1max = min4(max(c1lox, c1hix), max(c1loy, c1hiy), max(c1loz, c1hiz), t);
155
156 #ifdef __HAIR__
157         if(difl != 0.0f) {
158                 float hdiff = 1.0f + difl;
159                 float ldiff = 1.0f - difl;
160                 if(__float_as_int(cnodes.z) & PATH_RAY_CURVE) {
161                         c0min = max(ldiff * c0min, c0min - extmax);
162                         c0max = min(hdiff * c0max, c0max + extmax);
163                 }
164                 if(__float_as_int(cnodes.w) & PATH_RAY_CURVE) {
165                         c1min = max(ldiff * c1min, c1min - extmax);
166                         c1max = min(hdiff * c1max, c1max + extmax);
167                 }
168         }
169 #endif
170
171         /* decide which nodes to traverse next */
172 #ifdef __VISIBILITY_FLAG__
173         /* this visibility test gives a 5% performance hit, how to solve? */
174         *traverseChild0 = (c0max >= c0min) && (__float_as_uint(cnodes.z) & visibility);
175         *traverseChild1 = (c1max >= c1min) && (__float_as_uint(cnodes.w) & visibility);
176 #else
177         *traverseChild0 = (c0max >= c0min);
178         *traverseChild1 = (c1max >= c1min);
179 #endif
180
181         *nodeAddr0 = __float_as_int(cnodes.x);
182         *nodeAddr1 = __float_as_int(cnodes.y);
183
184         *closestChild1 = (c1min < c0min);
185 }
186
187 /* Sven Woop's algorithm */
188 __device_inline void bvh_triangle_intersect(KernelGlobals *kg, Intersection *isect,
189         float3 P, float3 idir, uint visibility, int object, int triAddr)
190 {
191         /* compute and check intersection t-value */
192         float4 v00 = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+0);
193         float4 v11 = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+1);
194         float3 dir = 1.0f/idir;
195
196         float Oz = v00.w - P.x*v00.x - P.y*v00.y - P.z*v00.z;
197         float invDz = 1.0f/(dir.x*v00.x + dir.y*v00.y + dir.z*v00.z);
198         float t = Oz * invDz;
199
200         if(t > 0.0f && t < isect->t) {
201                 /* compute and check barycentric u */
202                 float Ox = v11.w + P.x*v11.x + P.y*v11.y + P.z*v11.z;
203                 float Dx = dir.x*v11.x + dir.y*v11.y + dir.z*v11.z;
204                 float u = Ox + t*Dx;
205
206                 if(u >= 0.0f) {
207                         /* compute and check barycentric v */
208                         float4 v22 = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+2);
209                         float Oy = v22.w + P.x*v22.x + P.y*v22.y + P.z*v22.z;
210                         float Dy = dir.x*v22.x + dir.y*v22.y + dir.z*v22.z;
211                         float v = Oy + t*Dy;
212
213                         if(v >= 0.0f && u + v <= 1.0f) {
214 #ifdef __VISIBILITY_FLAG__
215                                 /* visibility flag test. we do it here under the assumption
216                                  * that most triangles are culled by node flags */
217                                 if(kernel_tex_fetch(__prim_visibility, triAddr) & visibility)
218 #endif
219                                 {
220                                         /* record intersection */
221                                         isect->prim = triAddr;
222                                         isect->object = object;
223                                         isect->u = u;
224                                         isect->v = v;
225                                         isect->t = t;
226                                 }
227                         }
228                 }
229         }
230 }
231
232 #ifdef __HAIR__
233 __device_inline void curvebounds(float *lower, float *upper, float *extremta, float *extrema, float *extremtb, float *extremb, float p0, float p1, float p2, float p3)
234 {
235         float halfdiscroot = (p2 * p2 - 3 * p3 * p1);
236         float ta = -1.0f;
237         float tb = -1.0f;
238         *extremta = -1.0f;
239         *extremtb = -1.0f;
240         *upper = p0;
241         *lower = p0 + p1 + p2 + p3;
242         *extrema = *upper;
243         *extremb = *lower;
244         if(*lower >= *upper) {
245                 *upper = *lower;
246                 *lower = p0;
247         }
248
249         if(halfdiscroot >= 0) {
250                 halfdiscroot = sqrt(halfdiscroot);
251                 ta = (-p2 - halfdiscroot) / (3 * p3);
252                 tb = (-p2 + halfdiscroot) / (3 * p3);
253         }
254
255         float t2;
256         float t3;
257         if(ta > 0.0f && ta < 1.0f) {
258                 t2 = ta * ta;
259                 t3 = t2 * ta;
260                 *extremta = ta;
261                 *extrema = p3 * t3 + p2 * t2 + p1 * ta + p0;
262                 if(*extrema > *upper) {
263                         *upper = *extrema;
264                 }
265                 if(*extrema < *lower) {
266                         *lower = *extrema;
267                 }
268         }
269         if(tb > 0.0f && tb < 1.0f) {
270                 t2 = tb * tb;
271                 t3 = t2 * tb;
272                 *extremtb = tb;
273                 *extremb = p3 * t3 + p2 * t2 + p1 * tb + p0;
274                 if(*extremb >= *upper) {
275                         *upper = *extremb;
276                 }
277                 if(*extremb <= *lower) {
278                         *lower = *extremb;
279                 }
280         }
281 }
282
283 __device_inline void bvh_cardinal_curve_intersect(KernelGlobals *kg, Intersection *isect,
284         float3 P, float3 idir, uint visibility, int object, int curveAddr, int segment, uint *lcg_state, float difl, float extmax)
285 {
286         float epsilon = 0.0f;
287         int depth = kernel_data.curve_kernel_data.subdivisions;
288
289         /* curve Intersection check */
290         float3 dir = 1.0f/idir;
291         
292         int flags = kernel_data.curve_kernel_data.curveflags;
293
294         int prim = kernel_tex_fetch(__prim_index, curveAddr);
295
296         float3 curve_coef[4];
297         float r_st,r_en;
298
299         /*obtain curve parameters*/
300         {
301                 /*ray transform created - this should be created at beginning of intersection loop*/
302                 Transform htfm;
303                 float d = sqrtf(dir.x * dir.x + dir.z * dir.z);
304                 htfm = make_transform(
305                         dir.z / d, 0, -dir.x /d, 0,
306                         -dir.x * dir.y /d, d, -dir.y * dir.z /d, 0,
307                         dir.x, dir.y, dir.z, 0,
308                         0, 0, 0, 1);
309
310                 float4 v00 = kernel_tex_fetch(__curves, prim);
311
312                 int k0 = __float_as_int(v00.x) + segment;
313                 int k1 = k0 + 1;
314
315                 int ka = max(k0 - 1,__float_as_int(v00.x));
316                 int kb = min(k1 + 1,__float_as_int(v00.x) + __float_as_int(v00.y) - 1);
317
318                 float4 P0 = kernel_tex_fetch(__curve_keys, ka);
319                 float4 P1 = kernel_tex_fetch(__curve_keys, k0);
320                 float4 P2 = kernel_tex_fetch(__curve_keys, k1);
321                 float4 P3 = kernel_tex_fetch(__curve_keys, kb);
322
323                 float3 p0 = transform_point(&htfm, float4_to_float3(P0) - P);
324                 float3 p1 = transform_point(&htfm, float4_to_float3(P1) - P);
325                 float3 p2 = transform_point(&htfm, float4_to_float3(P2) - P);
326                 float3 p3 = transform_point(&htfm, float4_to_float3(P3) - P);
327
328                 float fc = 0.71f;
329                 curve_coef[0] = p1;
330                 curve_coef[1] = -fc*p0 + fc*p2;
331                 curve_coef[2] = 2.0f * fc * p0 + (fc - 3.0f) * p1 + (3.0f - 2.0f * fc) * p2 - fc * p3;
332                 curve_coef[3] = -fc * p0 + (2.0f - fc) * p1 + (fc - 2.0f) * p2 + fc * p3;
333                 r_st = P1.w;
334                 r_en = P2.w;
335         }
336         
337
338         float r_curr = max(r_st, r_en);
339
340         if((flags & CURVE_KN_RIBBONS) || !(flags & CURVE_KN_BACKFACING))
341                 epsilon = 2 * r_curr;
342
343         /*find bounds - this is slow for cubic curves*/
344         float upper,lower;
345
346         float zextrem[4];
347         curvebounds(&lower, &upper, &zextrem[0], &zextrem[1], &zextrem[2], &zextrem[3], curve_coef[0].z, curve_coef[1].z, curve_coef[2].z, curve_coef[3].z);
348         if(lower - r_curr > isect->t || upper + r_curr < epsilon)
349                 return;
350
351         /*minimum width extension*/
352         float mw_extension = min(difl * fabsf(upper), extmax);
353         float r_ext = mw_extension + r_curr;
354
355         float xextrem[4];
356         curvebounds(&lower, &upper, &xextrem[0], &xextrem[1], &xextrem[2], &xextrem[3], curve_coef[0].x, curve_coef[1].x, curve_coef[2].x, curve_coef[3].x);
357         if(lower > r_ext || upper < -r_ext)
358                 return;
359
360         float yextrem[4];
361         curvebounds(&lower, &upper, &yextrem[0], &yextrem[1], &yextrem[2], &yextrem[3], curve_coef[0].y, curve_coef[1].y, curve_coef[2].y, curve_coef[3].y);
362         if(lower > r_ext || upper < -r_ext)
363                 return;
364
365         /*setup recurrent loop*/
366         int level = 1 << depth;
367         int tree = 0;
368         float resol = 1.0f / (float)level;
369
370         /*begin loop*/
371         while(!(tree >> (depth))) {
372                 float i_st = tree * resol;
373                 float i_en = i_st + (level * resol);
374                 float3 p_st = ((curve_coef[3] * i_st + curve_coef[2]) * i_st + curve_coef[1]) * i_st + curve_coef[0];
375                 float3 p_en = ((curve_coef[3] * i_en + curve_coef[2]) * i_en + curve_coef[1]) * i_en + curve_coef[0];
376                 
377                 float bminx = min(p_st.x, p_en.x);
378                 float bmaxx = max(p_st.x, p_en.x);
379                 float bminy = min(p_st.y, p_en.y);
380                 float bmaxy = max(p_st.y, p_en.y);
381                 float bminz = min(p_st.z, p_en.z);
382                 float bmaxz = max(p_st.z, p_en.z);
383
384                 if(xextrem[0] >= i_st && xextrem[0] <= i_en) {
385                         bminx = min(bminx,xextrem[1]);
386                         bmaxx = max(bmaxx,xextrem[1]);
387                 }
388                 if(xextrem[2] >= i_st && xextrem[2] <= i_en) {
389                         bminx = min(bminx,xextrem[3]);
390                         bmaxx = max(bmaxx,xextrem[3]);
391                 }
392                 if(yextrem[0] >= i_st && yextrem[0] <= i_en) {
393                         bminy = min(bminy,yextrem[1]);
394                         bmaxy = max(bmaxy,yextrem[1]);
395                 }
396                 if(yextrem[2] >= i_st && yextrem[2] <= i_en) {
397                         bminy = min(bminy,yextrem[3]);
398                         bmaxy = max(bmaxy,yextrem[3]);
399                 }
400                 if(zextrem[0] >= i_st && zextrem[0] <= i_en) {
401                         bminz = min(bminz,zextrem[1]);
402                         bmaxz = max(bmaxz,zextrem[1]);
403                 }
404                 if(zextrem[2] >= i_st && zextrem[2] <= i_en) {
405                         bminz = min(bminz,zextrem[3]);
406                         bmaxz = max(bmaxz,zextrem[3]);
407                 }
408
409                 float r1 = r_st + (r_en - r_st) * i_st;
410                 float r2 = r_st + (r_en - r_st) * i_en;
411                 r_curr = max(r1, r2);
412
413                 mw_extension = min(difl * fabsf(bmaxz), extmax);
414                 float r_ext = mw_extension + r_curr;
415                 float coverage = 1.0f;
416
417                 if (bminz - r_curr > isect->t || bmaxz + r_curr < epsilon || bminx > r_ext|| bmaxx < -r_ext|| bminy > r_ext|| bmaxy < -r_ext) {
418                         /* the bounding box does not overlap the square centered at O.*/
419                         tree += level;
420                         level = tree & -tree;
421                 }
422                 else if (level == 1) {
423
424                         /* the maximum recursion depth is reached.
425                         * check if dP0.(Q-P0)>=0 and dPn.(Pn-Q)>=0.
426                         * dP* is reversed if necessary.*/
427                         float t = isect->t;
428                         float u = 0.0f;
429                         if(flags & CURVE_KN_RIBBONS) {
430                                 float3 tg = (p_en - p_st);
431                                 float w = tg.x * tg.x + tg.y * tg.y;
432                                 if (w == 0) {
433                                         tree++;
434                                         level = tree & -tree;
435                                         continue;
436                                 }
437                                 w = -(p_st.x * tg.x + p_st.y * tg.y) / w;
438                                 w = clamp((float)w, 0.0f, 1.0f);
439
440                                 /* compute u on the curve segment.*/
441                                 u = i_st * (1 - w) + i_en * w;
442                                 r_curr = r_st + (r_en - r_st) * u;
443                                 /* compare x-y distances.*/
444                                 float3 p_curr = ((curve_coef[3] * u + curve_coef[2]) * u + curve_coef[1]) * u + curve_coef[0];
445
446                                 float3 dp_st = (3 * curve_coef[3] * i_st + 2 * curve_coef[2]) * i_st + curve_coef[1];
447                                 if (dot(tg, dp_st)< 0)
448                                         dp_st *= -1;
449                                 if (dot(dp_st, -p_st) + p_curr.z * dp_st.z < 0) {
450                                         tree++;
451                                         level = tree & -tree;
452                                         continue;
453                                 }
454                                 float3 dp_en = (3 * curve_coef[3] * i_en + 2 * curve_coef[2]) * i_en + curve_coef[1];
455                                 if (dot(tg, dp_en) < 0)
456                                         dp_en *= -1;
457                                 if (dot(dp_en, p_en) - p_curr.z * dp_en.z < 0) {
458                                         tree++;
459                                         level = tree & -tree;
460                                         continue;
461                                 }
462
463                                 /* compute coverage */
464                                 float r_ext = r_curr;
465                                 coverage = 1.0f;
466                                 if(difl != 0.0f) {
467                                         mw_extension = min(difl * fabsf(bmaxz), extmax);
468                                         r_ext = mw_extension + r_curr;
469                                         float d = sqrtf(p_curr.x * p_curr.x + p_curr.y * p_curr.y);
470                                         float d0 = d - r_curr;
471                                         float d1 = d + r_curr;
472                                         if (d0 >= 0)
473                                                 coverage = (min(d1 / mw_extension, 1.0f) - min(d0 / mw_extension, 1.0f)) * 0.5f;
474                                         else // inside
475                                                 coverage = (min(d1 / mw_extension, 1.0f) + min(-d0 / mw_extension, 1.0f)) * 0.5f;
476                                 }
477                                 
478                                 if (p_curr.x * p_curr.x + p_curr.y * p_curr.y >= r_ext * r_ext || p_curr.z <= epsilon) {
479                                         tree++;
480                                         level = tree & -tree;
481                                         continue;
482                                 }
483                                 /* compare z distances.*/
484                                 if (isect->t < p_curr.z) {
485                                         tree++;
486                                         level = tree & -tree;
487                                         continue;
488                                 }
489                                 t = p_curr.z;
490                         }
491                         else {
492                                 float l = len(p_en - p_st);
493                                 /*minimum width extension*/
494                                 float or1 = r1;
495                                 float or2 = r2;
496                                 if(difl != 0.0f) {
497                                         mw_extension = min(len(p_st - P) * difl, extmax);
498                                         or1 = r1 < mw_extension ? mw_extension : r1;
499                                         mw_extension = min(len(p_en - P) * difl, extmax);
500                                         or2 = r2 < mw_extension ? mw_extension : r2;
501                                 }
502                                 /* --- */
503                                 float3 tg = (p_en - p_st) / l;
504                                 float gd = (or2 - or1) / l;
505                                 float difz = -dot(p_st,tg);
506                                 float cyla = 1.0f - (tg.z * tg.z * (1 + gd*gd));
507                                 float halfb = (-p_st.z - tg.z*(difz + gd*(difz*gd + or1)));
508                                 float tcentre = -halfb/cyla;
509                                 float zcentre = difz + (tg.z * tcentre);
510                                 float3 tdif = - p_st;
511                                 tdif.z += tcentre;
512                                 float tdifz = dot(tdif,tg);
513                                 float tb = 2*(tdif.z - tg.z*(tdifz + gd*(tdifz*gd + or1)));
514                                 float tc = dot(tdif,tdif) - tdifz * tdifz * (1 + gd*gd) - or1*or1 - 2*or1*tdifz*gd;
515                                 float td = tb*tb - 4*cyla*tc;
516                                 if (td < 0.0f){
517                                         tree++;
518                                         level = tree & -tree;
519                                         continue;
520                                 }
521                                 
522                                 float rootd = sqrtf(td);
523                                 float correction = ((-tb - rootd)/(2*cyla));
524                                 t = tcentre + correction;
525                                 float w = (zcentre + (tg.z * correction))/l;
526
527                                 float3 dp_st = (3 * curve_coef[3] * i_st + 2 * curve_coef[2]) * i_st + curve_coef[1];
528                                 if (dot(tg, dp_st)< 0)
529                                         dp_st *= -1;
530                                 float3 dp_en = (3 * curve_coef[3] * i_en + 2 * curve_coef[2]) * i_en + curve_coef[1];
531                                 if (dot(tg, dp_en) < 0)
532                                         dp_en *= -1;
533
534
535                                 if(flags & CURVE_KN_BACKFACING && (dot(dp_st, -p_st) + t * dp_st.z < 0 || dot(dp_en, p_en) - t * dp_en.z < 0 || isect->t < t || t <= 0.0f)) {
536                                         correction = ((-tb + rootd)/(2*cyla));
537                                         t = tcentre + correction;
538                                         w = (zcentre + (tg.z * correction))/l;
539                                 }                       
540
541                                 if (dot(dp_st, -p_st) + t * dp_st.z < 0 || dot(dp_en, p_en) - t * dp_en.z < 0 || isect->t < t || t <= 0.0f) {
542                                         tree++;
543                                         level = tree & -tree;
544                                         continue;
545                                 }
546
547                                 w = clamp((float)w, 0.0f, 1.0f);
548                                 /* compute u on the curve segment.*/
549                                 u = i_st * (1 - w) + i_en * w;
550                                 r_curr = r1 + (r2 - r1) * w;
551                                 r_ext = or1 + (or2 - or1) * w;
552                                 coverage = r_curr/r_ext;
553
554                         }
555                         /* we found a new intersection.*/
556
557                         /*stochastic fade from minimum width*/
558                         if(lcg_state && coverage != 1.0f) {
559                                 if(lcg_step(lcg_state) > coverage)
560                                         return;
561                         }
562
563 #ifdef __VISIBILITY_FLAG__
564                         /* visibility flag test. we do it here under the assumption
565                          * that most triangles are culled by node flags */
566                         if(kernel_tex_fetch(__prim_visibility, curveAddr) & visibility)
567 #endif
568                         {
569                                 /* record intersection */
570                                 isect->prim = curveAddr;
571                                 isect->segment = segment;
572                                 isect->object = object;
573                                 isect->u = u;
574                                 isect->v = 0.0f;
575                                 /*isect->v = 1.0f - coverage; */
576                                 isect->t = t;
577                         }
578                         
579                         tree++;
580                         level = tree & -tree;
581                 }
582                 else {
583                         /* split the curve into two curves and process */
584                         level = level >> 1;
585                 }
586         }
587 }
588
589 __device_inline void bvh_curve_intersect(KernelGlobals *kg, Intersection *isect,
590         float3 P, float3 idir, uint visibility, int object, int curveAddr, int segment, uint *lcg_state, float difl, float extmax)
591 {
592         /* curve Intersection check */
593         int flags = kernel_data.curve_kernel_data.curveflags;
594
595         int prim = kernel_tex_fetch(__prim_index, curveAddr);
596         float4 v00 = kernel_tex_fetch(__curves, prim);
597
598         int cnum = __float_as_int(v00.x);
599         int k0 = cnum + segment;
600         int k1 = k0 + 1;
601
602         float4 P1 = kernel_tex_fetch(__curve_keys, k0);
603         float4 P2 = kernel_tex_fetch(__curve_keys, k1);
604
605         float or1 = P1.w;
606         float or2 = P2.w;
607         float3 p1 = float4_to_float3(P1);
608         float3 p2 = float4_to_float3(P2);
609
610         /*minimum width extension*/
611         float r1 = or1;
612         float r2 = or2;
613         if(difl != 0.0f) {
614                 float pixelsize = min(len(p1 - P) * difl, extmax);
615                 r1 = or1 < pixelsize ? pixelsize : or1;
616                 pixelsize = min(len(p2 - P) * difl, extmax);
617                 r2 = or2 < pixelsize ? pixelsize : or2;
618         }
619         /* --- */
620
621         float mr = max(r1,r2);
622         float3 dif = P - p1;
623         float3 dir = 1.0f/idir;
624         float l = len(p2 - p1);
625
626         float sp_r = mr + 0.5f * l;
627         float3 sphere_dif = P - ((p1 + p2) * 0.5f);
628         float sphere_b = dot(dir,sphere_dif);
629         sphere_dif = sphere_dif - sphere_b * dir;
630         sphere_b = dot(dir,sphere_dif);
631         float sdisc = sphere_b * sphere_b - len_squared(sphere_dif) + sp_r * sp_r;
632         if(sdisc < 0.0f)
633                 return;
634
635         /* obtain parameters and test midpoint distance for suitable modes*/
636         float3 tg = (p2 - p1) / l;
637         float gd = (r2 - r1) / l;
638         float dirz = dot(dir,tg);
639         float difz = dot(dif,tg);
640
641         float a = 1.0f - (dirz*dirz*(1 + gd*gd));
642         float halfb = dot(dir,dif) - dirz*(difz + gd*(difz*gd + r1));
643
644         float tcentre = -halfb/a;
645         float zcentre = difz + (dirz * tcentre);
646
647         if((tcentre > isect->t) && !(flags & CURVE_KN_ACCURATE))
648                 return;
649         if((zcentre < 0 || zcentre > l) && !(flags & CURVE_KN_ACCURATE) && !(flags & CURVE_KN_INTERSECTCORRECTION))
650                 return;
651
652         /* test minimum separation*/
653         float3 cprod = cross(tg, dir);
654         float3 cprod2 = cross(tg, dif);
655         float cprodsq = len_squared(cprod);
656         float cprod2sq = len_squared(cprod2);
657         float distscaled = dot(cprod,dif);
658
659         if(cprodsq == 0)
660                 distscaled = cprod2sq;
661         else
662                 distscaled = (distscaled*distscaled)/cprodsq;
663
664         if(distscaled > mr*mr)
665                 return;
666
667         /* calculate true intersection*/
668         float3 tdif = P - p1 + tcentre * dir;
669         float tdifz = dot(tdif,tg);
670         float tb = 2*(dot(dir,tdif) - dirz*(tdifz + gd*(tdifz*gd + r1)));
671         float tc = dot(tdif,tdif) - tdifz * tdifz * (1 + gd*gd) - r1*r1 - 2*r1*tdifz*gd;
672         float td = tb*tb - 4*a*tc;
673
674         if (td < 0.0f)
675                 return;
676
677         float rootd = 0.0f;
678         float correction = 0.0f;
679         if(flags & CURVE_KN_ACCURATE) {
680                 rootd = sqrtf(td);
681                 correction = ((-tb - rootd)/(2*a));
682         }
683
684         float t = tcentre + correction;
685
686         if(t < isect->t) {
687
688                 if(flags & CURVE_KN_INTERSECTCORRECTION) {
689                         rootd = sqrtf(td);
690                         correction = ((-tb - rootd)/(2*a));
691                         t = tcentre + correction;
692                 }
693
694                 float z = zcentre + (dirz * correction);
695                 bool backface = false;
696
697                 if(flags & CURVE_KN_BACKFACING && (t < 0.0f || z < 0 || z > l)) {
698                         backface = true;
699                         correction = ((-tb + rootd)/(2*a));
700                         t = tcentre + correction;
701                         z = zcentre + (dirz * correction);
702                 }
703
704                 /*stochastic fade from minimum width*/
705                 float adjradius = or1 + z * (or2 - or1) / l;
706                 adjradius = adjradius / (r1 + z * gd);
707                 if(lcg_state && adjradius != 1.0f) {
708                         if(lcg_step(lcg_state) > adjradius)
709                                 return;
710                 }
711                 /* --- */
712
713                 if(t > 0.0f && t < isect->t && z >= 0 && z <= l) {
714
715                         if (flags & CURVE_KN_ENCLOSEFILTER) {
716
717                                 float enc_ratio = kernel_data.curve_kernel_data.encasing_ratio;
718                                 if((dot(P - p1, tg) > -r1 * enc_ratio) && (dot(P - p2, tg) < r2 * enc_ratio)) {
719                                         float a2 = 1.0f - (dirz*dirz*(1 + gd*gd*enc_ratio*enc_ratio));
720                                         float c2 = dot(dif,dif) - difz * difz * (1 + gd*gd*enc_ratio*enc_ratio) - r1*r1*enc_ratio*enc_ratio - 2*r1*difz*gd*enc_ratio;
721                                         if(a2*c2 < 0.0f)
722                                                 return;
723                                 }
724                         }
725
726 #ifdef __VISIBILITY_FLAG__
727                         /* visibility flag test. we do it here under the assumption
728                          * that most triangles are culled by node flags */
729                         if(kernel_tex_fetch(__prim_visibility, curveAddr) & visibility)
730 #endif
731                         {
732                                 /* record intersection */
733                                 isect->prim = curveAddr;
734                                 isect->segment = segment;
735                                 isect->object = object;
736                                 isect->u = z/l;
737                                 isect->v = td/(4*a*a);
738                                 /*isect->v = 1.0f - adjradius;*/
739                                 isect->t = t;
740
741                                 if(backface) 
742                                         isect->u = -isect->u;
743                         }
744                 }
745         }
746 }
747 #endif
748
749 #ifdef __SUBSURFACE__
750 /* Special ray intersection routines for subsurface scattering. In that case we
751  * only want to intersect with primitives in the same object, and if case of
752  * multiple hits we pick a single random primitive as the intersection point. */
753
754 __device_inline void bvh_triangle_intersect_subsurface(KernelGlobals *kg, Intersection *isect,
755         float3 P, float3 idir, int object, int triAddr, float tmax, int *num_hits, float subsurface_random)
756 {
757         /* compute and check intersection t-value */
758         float4 v00 = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+0);
759         float4 v11 = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+1);
760         float3 dir = 1.0f/idir;
761
762         float Oz = v00.w - P.x*v00.x - P.y*v00.y - P.z*v00.z;
763         float invDz = 1.0f/(dir.x*v00.x + dir.y*v00.y + dir.z*v00.z);
764         float t = Oz * invDz;
765
766         if(t > 0.0f && t < tmax) {
767                 /* compute and check barycentric u */
768                 float Ox = v11.w + P.x*v11.x + P.y*v11.y + P.z*v11.z;
769                 float Dx = dir.x*v11.x + dir.y*v11.y + dir.z*v11.z;
770                 float u = Ox + t*Dx;
771
772                 if(u >= 0.0f) {
773                         /* compute and check barycentric v */
774                         float4 v22 = kernel_tex_fetch(__tri_woop, triAddr*TRI_NODE_SIZE+2);
775                         float Oy = v22.w + P.x*v22.x + P.y*v22.y + P.z*v22.z;
776                         float Dy = dir.x*v22.x + dir.y*v22.y + dir.z*v22.z;
777                         float v = Oy + t*Dy;
778
779                         if(v >= 0.0f && u + v <= 1.0f) {
780                                 (*num_hits)++;
781
782                                 if(subsurface_random * (*num_hits) <= 1.0f) {
783                                         /* record intersection */
784                                         isect->prim = triAddr;
785                                         isect->object = object;
786                                         isect->u = u;
787                                         isect->v = v;
788                                         isect->t = t;
789                                 }
790                         }
791                 }
792         }
793 }
794 #endif
795
796 /* BVH intersection function variations */
797
798 #define BVH_INSTANCING                  1
799 #define BVH_MOTION                              2
800 #define BVH_HAIR                                4
801 #define BVH_HAIR_MINIMUM_WIDTH  8
802 #define BVH_SUBSURFACE                  16
803
804 #define BVH_FUNCTION_NAME bvh_intersect
805 #define BVH_FUNCTION_FEATURES 0
806 #include "kernel_bvh_traversal.h"
807
808 #if defined(__INSTANCING__)
809 #define BVH_FUNCTION_NAME bvh_intersect_instancing
810 #define BVH_FUNCTION_FEATURES BVH_INSTANCING
811 #include "kernel_bvh_traversal.h"
812 #endif
813
814 #if defined(__HAIR__)
815 #define BVH_FUNCTION_NAME bvh_intersect_hair
816 #define BVH_FUNCTION_FEATURES BVH_INSTANCING|BVH_HAIR|BVH_HAIR_MINIMUM_WIDTH
817 #include "kernel_bvh_traversal.h"
818 #endif
819
820 #if defined(__OBJECT_MOTION__)
821 #define BVH_FUNCTION_NAME bvh_intersect_motion
822 #define BVH_FUNCTION_FEATURES BVH_INSTANCING|BVH_MOTION
823 #include "kernel_bvh_traversal.h"
824 #endif
825
826 #if defined(__HAIR__) && defined(__OBJECT_MOTION__)
827 #define BVH_FUNCTION_NAME bvh_intersect_hair_motion
828 #define BVH_FUNCTION_FEATURES BVH_INSTANCING|BVH_HAIR|BVH_HAIR_MINIMUM_WIDTH|BVH_MOTION
829 #include "kernel_bvh_traversal.h"
830 #endif
831
832 #if defined(__SUBSURFACE__)
833 #define BVH_FUNCTION_NAME bvh_intersect_subsurface
834 #define BVH_FUNCTION_FEATURES BVH_SUBSURFACE
835 #include "kernel_bvh_traversal.h"
836 #endif
837
838 #if defined(__SUBSURFACE__) && defined(__INSTANCING__)
839 #define BVH_FUNCTION_NAME bvh_intersect_subsurface_instancing
840 #define BVH_FUNCTION_FEATURES BVH_INSTANCING|BVH_SUBSURFACE
841 #include "kernel_bvh_traversal.h"
842 #endif
843
844 #if defined(__SUBSURFACE__) && defined(__HAIR__)
845 #define BVH_FUNCTION_NAME bvh_intersect_subsurface_hair
846 #define BVH_FUNCTION_FEATURES BVH_INSTANCING|BVH_SUBSURFACE|BVH_HAIR|BVH_HAIR_MINIMUM_WIDTH
847 #include "kernel_bvh_traversal.h"
848 #endif
849
850 #if defined(__SUBSURFACE__) && defined(__OBJECT_MOTION__)
851 #define BVH_FUNCTION_NAME bvh_intersect_subsurface_motion
852 #define BVH_FUNCTION_FEATURES BVH_INSTANCING|BVH_SUBSURFACE|BVH_MOTION
853 #include "kernel_bvh_traversal.h"
854 #endif
855
856 #if defined(__SUBSURFACE__) && defined(__HAIR__) && defined(__OBJECT_MOTION__)
857 #define BVH_FUNCTION_NAME bvh_intersect_subsurface_hair_motion
858 #define BVH_FUNCTION_FEATURES BVH_INSTANCING|BVH_SUBSURFACE|BVH_HAIR|BVH_HAIR_MINIMUM_WIDTH|BVH_MOTION
859 #include "kernel_bvh_traversal.h"
860 #endif
861
862
863 #ifdef __HAIR__ 
864 __device_inline bool scene_intersect(KernelGlobals *kg, const Ray *ray, const uint visibility, Intersection *isect, uint *lcg_state, float difl, float extmax)
865 #else
866 __device_inline bool scene_intersect(KernelGlobals *kg, const Ray *ray, const uint visibility, Intersection *isect)
867 #endif
868 {
869 #ifdef __OBJECT_MOTION__
870         if(kernel_data.bvh.have_motion) {
871 #ifdef __HAIR__
872                 if(kernel_data.bvh.have_curves)
873                         return bvh_intersect_hair_motion(kg, ray, isect, visibility, lcg_state, difl, extmax);
874 #endif /* __HAIR__ */
875
876                 return bvh_intersect_motion(kg, ray, isect, visibility);
877         }
878 #endif /* __OBJECT_MOTION__ */
879
880 #ifdef __HAIR__ 
881         if(kernel_data.bvh.have_curves)
882                 return bvh_intersect_hair(kg, ray, isect, visibility, lcg_state, difl, extmax);
883 #endif /* __HAIR__ */
884
885 #ifdef __KERNEL_CPU__
886
887 #ifdef __INSTANCING__
888         if(kernel_data.bvh.have_instancing)
889                 return bvh_intersect_instancing(kg, ray, isect, visibility);
890 #endif /* __INSTANCING__ */
891
892         return bvh_intersect(kg, ray, isect, visibility);
893 #else /* __KERNEL_CPU__ */
894
895 #ifdef __INSTANCING__
896         return bvh_intersect_instancing(kg, ray, isect, visibility);
897 #else
898         return bvh_intersect(kg, ray, isect, visibility);
899 #endif /* __INSTANCING__ */
900
901 #endif /* __KERNEL_CPU__ */
902 }
903
904 #ifdef __SUBSURFACE__
905 __device_inline int scene_intersect_subsurface(KernelGlobals *kg, const Ray *ray, Intersection *isect, int subsurface_object, float subsurface_random)
906 {
907 #ifdef __OBJECT_MOTION__
908         if(kernel_data.bvh.have_motion) {
909 #ifdef __HAIR__
910                 if(kernel_data.bvh.have_curves)
911                         return bvh_intersect_subsurface_hair_motion(kg, ray, isect, subsurface_object, subsurface_random);
912 #endif /* __HAIR__ */
913
914                 return bvh_intersect_subsurface_motion(kg, ray, isect, subsurface_object, subsurface_random);
915         }
916 #endif /* __OBJECT_MOTION__ */
917
918 #ifdef __HAIR__ 
919         if(kernel_data.bvh.have_curves)
920                 return bvh_intersect_subsurface_hair(kg, ray, isect, subsurface_object, subsurface_random);
921 #endif /* __HAIR__ */
922
923 #ifdef __KERNEL_CPU__
924
925 #ifdef __INSTANCING__
926         if(kernel_data.bvh.have_instancing)
927                 return bvh_intersect_subsurface_instancing(kg, ray, isect, subsurface_object, subsurface_random);
928 #endif /* __INSTANCING__ */
929
930         return bvh_intersect_subsurface(kg, ray, isect, subsurface_object, subsurface_random);
931 #else /* __KERNEL_CPU__ */
932
933 #ifdef __INSTANCING__
934         return bvh_intersect_subsurface_instancing(kg, ray, isect, subsurface_object, subsurface_random);
935 #else
936         return bvh_intersect_subsurface(kg, ray, isect, subsurface_object, subsurface_random);
937 #endif /* __INSTANCING__ */
938
939 #endif /* __KERNEL_CPU__ */
940 }
941 #endif
942
943 /* Ray offset to avoid self intersection */
944
945 __device_inline float3 ray_offset(float3 P, float3 Ng)
946 {
947 #ifdef __INTERSECTION_REFINE__
948         const float epsilon_f = 1e-5f;
949         /* ideally this should match epsilon_f, but instancing/mblur
950          * precision makes it problematic */
951         const float epsilon_test = 1.0f;
952         const int epsilon_i = 32;
953
954         float3 res;
955
956         /* x component */
957         if(fabsf(P.x) < epsilon_test) {
958                 res.x = P.x + Ng.x*epsilon_f;
959         }
960         else {
961                 uint ix = __float_as_uint(P.x);
962                 ix += ((ix ^ __float_as_uint(Ng.x)) >> 31)? -epsilon_i: epsilon_i;
963                 res.x = __uint_as_float(ix);
964         }
965
966         /* y component */
967         if(fabsf(P.y) < epsilon_test) {
968                 res.y = P.y + Ng.y*epsilon_f;
969         }
970         else {
971                 uint iy = __float_as_uint(P.y);
972                 iy += ((iy ^ __float_as_uint(Ng.y)) >> 31)? -epsilon_i: epsilon_i;
973                 res.y = __uint_as_float(iy);
974         }
975
976         /* z component */
977         if(fabsf(P.z) < epsilon_test) {
978                 res.z = P.z + Ng.z*epsilon_f;
979         }
980         else {
981                 uint iz = __float_as_uint(P.z);
982                 iz += ((iz ^ __float_as_uint(Ng.z)) >> 31)? -epsilon_i: epsilon_i;
983                 res.z = __uint_as_float(iz);
984         }
985
986         return res;
987 #else
988         const float epsilon_f = 1e-4f;
989         return P + epsilon_f*Ng;
990 #endif
991 }
992
993 /* Refine triangle intersection to more precise hit point. For rays that travel
994  * far the precision is often not so good, this reintersects the primitive from
995  * a closer distance. */
996
997 __device_inline float3 bvh_triangle_refine(KernelGlobals *kg, ShaderData *sd, const Intersection *isect, const Ray *ray)
998 {
999         float3 P = ray->P;
1000         float3 D = ray->D;
1001         float t = isect->t;
1002
1003 #ifdef __INTERSECTION_REFINE__
1004         if(isect->object != ~0) {
1005 #ifdef __OBJECT_MOTION__
1006                 Transform tfm = sd->ob_itfm;
1007 #else
1008                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
1009 #endif
1010
1011                 P = transform_point(&tfm, P);
1012                 D = transform_direction(&tfm, D*t);
1013                 D = normalize_len(D, &t);
1014         }
1015
1016         P = P + D*t;
1017
1018         float4 v00 = kernel_tex_fetch(__tri_woop, isect->prim*TRI_NODE_SIZE+0);
1019         float Oz = v00.w - P.x*v00.x - P.y*v00.y - P.z*v00.z;
1020         float invDz = 1.0f/(D.x*v00.x + D.y*v00.y + D.z*v00.z);
1021         float rt = Oz * invDz;
1022
1023         P = P + D*rt;
1024
1025         if(isect->object != ~0) {
1026 #ifdef __OBJECT_MOTION__
1027                 Transform tfm = sd->ob_tfm;
1028 #else
1029                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
1030 #endif
1031
1032                 P = transform_point(&tfm, P);
1033         }
1034
1035         return P;
1036 #else
1037         return P + D*t;
1038 #endif
1039 }
1040
1041 #ifdef __HAIR__
1042
1043 __device_inline float3 curvetangent(float t, float3 p0, float3 p1, float3 p2, float3 p3)
1044 {
1045         float fc = 0.71f;
1046         float data[4];
1047         float t2 = t * t;
1048     data[0] = -3.0f * fc          * t2  + 4.0f * fc * t                  - fc;
1049     data[1] =  3.0f * (2.0f - fc) * t2  + 2.0f * (fc - 3.0f) * t;
1050     data[2] =  3.0f * (fc - 2.0f) * t2  + 2.0f * (3.0f - 2.0f * fc) * t  + fc;
1051     data[3] =  3.0f * fc          * t2  - 2.0f * fc * t;
1052         return data[0] * p0 + data[1] * p1 + data[2] * p2 + data[3] * p3;
1053 }
1054
1055 __device_inline float3 curvepoint(float t, float3 p0, float3 p1, float3 p2, float3 p3)
1056 {
1057         float data[4];
1058         float fc = 0.71f;
1059         float t2 = t * t;
1060         float t3 = t2 * t;
1061         data[0] = -fc          * t3  + 2.0f * fc          * t2 - fc * t;
1062         data[1] =  (2.0f - fc) * t3  + (fc - 3.0f)        * t2 + 1.0f;
1063         data[2] =  (fc - 2.0f) * t3  + (3.0f - 2.0f * fc) * t2 + fc * t;
1064         data[3] =  fc          * t3  - fc * t2;
1065         return data[0] * p0 + data[1] * p1 + data[2] * p2 + data[3] * p3;
1066 }
1067
1068 __device_inline float3 bvh_curve_refine(KernelGlobals *kg, ShaderData *sd, const Intersection *isect, const Ray *ray, float t)
1069 {
1070         int flag = kernel_data.curve_kernel_data.curveflags;
1071         float3 P = ray->P;
1072         float3 D = ray->D;
1073
1074         if(isect->object != ~0) {
1075 #ifdef __OBJECT_MOTION__
1076                 Transform tfm = sd->ob_itfm;
1077 #else
1078                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
1079 #endif
1080
1081                 P = transform_point(&tfm, P);
1082                 D = transform_direction(&tfm, D*t);
1083                 D = normalize_len(D, &t);
1084         }
1085
1086         int prim = kernel_tex_fetch(__prim_index, isect->prim);
1087         float4 v00 = kernel_tex_fetch(__curves, prim);
1088
1089         int k0 = __float_as_int(v00.x) + isect->segment;
1090         int k1 = k0 + 1;
1091
1092         float4 P1 = kernel_tex_fetch(__curve_keys, k0);
1093         float4 P2 = kernel_tex_fetch(__curve_keys, k1);
1094         float l = 1.0f;
1095         float3 tg = normalize_len(float4_to_float3(P2 - P1),&l);
1096         float r1 = P1.w;
1097         float r2 = P2.w;
1098         float gd = ((r2 - r1)/l);
1099         
1100         P = P + D*t;
1101
1102         if(flag & CURVE_KN_INTERPOLATE) {
1103                 int ka = max(k0 - 1,__float_as_int(v00.x));
1104                 int kb = min(k1 + 1,__float_as_int(v00.x) + __float_as_int(v00.y) - 1);
1105
1106                 float4 P0 = kernel_tex_fetch(__curve_keys, ka);
1107                 float4 P3 = kernel_tex_fetch(__curve_keys, kb);
1108
1109                 float3 p[4];
1110                 p[0] = float4_to_float3(P0);
1111                 p[1] = float4_to_float3(P1);
1112                 p[2] = float4_to_float3(P2);
1113                 p[3] = float4_to_float3(P3);
1114
1115                 tg = normalize(curvetangent(isect->u,p[0],p[1],p[2],p[3]));
1116                 float3 p_curr = curvepoint(isect->u,p[0],p[1],p[2],p[3]);
1117
1118 #ifdef __UV__
1119                 sd->u = isect->u;
1120                 sd->v = 0.0f;
1121 #endif
1122
1123                 if(kernel_data.curve_kernel_data.curveflags & CURVE_KN_RIBBONS)
1124                         sd->Ng = normalize(-(D - tg * (dot(tg,D))));
1125                 else {
1126                         sd->Ng = normalize(P - p_curr);
1127                         sd->Ng = sd->Ng - gd * tg;
1128                         sd->Ng = normalize(sd->Ng);
1129                 }
1130                 sd->N = sd->Ng;
1131         }
1132         else {
1133                 float3 dif = P - float4_to_float3(P1);
1134
1135 #ifdef __UV__
1136                 sd->u = dot(dif,tg)/l;
1137                 sd->v = 0.0f;
1138 #endif
1139
1140                 if (flag & CURVE_KN_TRUETANGENTGNORMAL) {
1141                         sd->Ng = -(D - tg * (dot(tg,D) * kernel_data.curve_kernel_data.normalmix));
1142                         sd->Ng = normalize(sd->Ng);
1143                         if (flag & CURVE_KN_NORMALCORRECTION) {
1144                                 sd->Ng = sd->Ng - gd * tg;
1145                                 sd->Ng = normalize(sd->Ng);
1146                         }
1147                 }
1148                 else {
1149                         sd->Ng = (dif - tg * sd->u * l) / (P1.w + sd->u * l * gd);
1150                         if (gd != 0.0f) {
1151                                 sd->Ng = sd->Ng - gd * tg ;
1152                                 sd->Ng = normalize(sd->Ng);
1153                         }
1154                 }
1155
1156                 sd->N = sd->Ng;
1157
1158                 if (flag & CURVE_KN_TANGENTGNORMAL && !(flag & CURVE_KN_TRUETANGENTGNORMAL)) {
1159                         sd->N = -(D - tg * (dot(tg,D) * kernel_data.curve_kernel_data.normalmix));
1160                         sd->N = normalize(sd->N);
1161                         if (flag & CURVE_KN_NORMALCORRECTION) {
1162                                 sd->N = sd->N - gd * tg;
1163                                 sd->N = normalize(sd->N);
1164                         }
1165                 }
1166                 if (!(flag & CURVE_KN_TANGENTGNORMAL) && flag & CURVE_KN_TRUETANGENTGNORMAL) {
1167                         sd->N = (dif - tg * sd->u * l) / (P1.w + sd->u * l * gd);
1168                         if (gd != 0.0f) {
1169                                 sd->N = sd->N - gd * tg ;
1170                                 sd->N = normalize(sd->N);
1171                         }
1172                 }
1173         }
1174
1175 #ifdef __DPDU__
1176         /* dPdu/dPdv */
1177         sd->dPdu = tg;
1178         sd->dPdv = cross(tg,sd->Ng);
1179 #endif
1180
1181         /*add fading parameter for minimum pixel width with transparency bsdf*/
1182         /*sd->curve_transparency = isect->v;*/
1183         /*sd->curve_radius = sd->u * gd * l + r1;*/
1184
1185         if(isect->object != ~0) {
1186 #ifdef __OBJECT_MOTION__
1187                 Transform tfm = sd->ob_tfm;
1188 #else
1189                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
1190 #endif
1191
1192                 P = transform_point(&tfm, P);
1193         }
1194
1195         return P;
1196 }
1197 #endif
1198
1199 CCL_NAMESPACE_END
1200