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