8888000f0e65740d5af83f7b9ef31f4f826024a0
[blender.git] / intern / cycles / kernel / geom / geom_curve.h
1 /*
2  * Licensed under the Apache License, Version 2.0 (the "License");
3  * you may not use this file except in compliance with the License.
4  * You may obtain a copy of the License at
5  *
6  * http://www.apache.org/licenses/LICENSE-2.0
7  *
8  * Unless required by applicable law or agreed to in writing, software
9  * distributed under the License is distributed on an "AS IS" BASIS,
10  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11  * See the License for the specific language governing permissions and
12  * limitations under the License.
13  */
14
15 CCL_NAMESPACE_BEGIN
16
17 /* Curve Primitive
18  *
19  * Curve primitive for rendering hair and fur. These can be render as flat ribbons
20  * or curves with actual thickness. The curve can also be rendered as line segments
21  * rather than curves for better performance */
22
23 #ifdef __HAIR__
24
25 #if defined(__KERNEL_CUDA__) && (__CUDA_ARCH__ < 300)
26 #  define ccl_device_curveintersect ccl_device
27 #else
28 #  define ccl_device_curveintersect ccl_device_forceinline
29 #endif
30
31 /* Reading attributes on various curve elements */
32
33 ccl_device float curve_attribute_float(KernelGlobals *kg, const ShaderData *sd, const AttributeDescriptor desc, float *dx, float *dy)
34 {
35         if(desc.element == ATTR_ELEMENT_CURVE) {
36 #ifdef __RAY_DIFFERENTIALS__
37                 if(dx) *dx = 0.0f;
38                 if(dy) *dy = 0.0f;
39 #endif
40
41                 return kernel_tex_fetch(__attributes_float, desc.offset + sd->prim);
42         }
43         else if(desc.element == ATTR_ELEMENT_CURVE_KEY || desc.element == ATTR_ELEMENT_CURVE_KEY_MOTION) {
44                 float4 curvedata = kernel_tex_fetch(__curves, sd->prim);
45                 int k0 = __float_as_int(curvedata.x) + PRIMITIVE_UNPACK_SEGMENT(sd->type);
46                 int k1 = k0 + 1;
47
48                 float f0 = kernel_tex_fetch(__attributes_float, desc.offset + k0);
49                 float f1 = kernel_tex_fetch(__attributes_float, desc.offset + k1);
50
51 #ifdef __RAY_DIFFERENTIALS__
52                 if(dx) *dx = sd->du.dx*(f1 - f0);
53                 if(dy) *dy = 0.0f;
54 #endif
55
56                 return (1.0f - sd->u)*f0 + sd->u*f1;
57         }
58         else {
59 #ifdef __RAY_DIFFERENTIALS__
60                 if(dx) *dx = 0.0f;
61                 if(dy) *dy = 0.0f;
62 #endif
63
64                 return 0.0f;
65         }
66 }
67
68 ccl_device float3 curve_attribute_float3(KernelGlobals *kg, const ShaderData *sd, const AttributeDescriptor desc, float3 *dx, float3 *dy)
69 {
70         if(desc.element == ATTR_ELEMENT_CURVE) {
71                 /* idea: we can't derive any useful differentials here, but for tiled
72                  * mipmap image caching it would be useful to avoid reading the highest
73                  * detail level always. maybe a derivative based on the hair density
74                  * could be computed somehow? */
75 #ifdef __RAY_DIFFERENTIALS__
76                 if(dx) *dx = make_float3(0.0f, 0.0f, 0.0f);
77                 if(dy) *dy = make_float3(0.0f, 0.0f, 0.0f);
78 #endif
79
80                 return float4_to_float3(kernel_tex_fetch(__attributes_float3, desc.offset + sd->prim));
81         }
82         else if(desc.element == ATTR_ELEMENT_CURVE_KEY || desc.element == ATTR_ELEMENT_CURVE_KEY_MOTION) {
83                 float4 curvedata = kernel_tex_fetch(__curves, sd->prim);
84                 int k0 = __float_as_int(curvedata.x) + PRIMITIVE_UNPACK_SEGMENT(sd->type);
85                 int k1 = k0 + 1;
86
87                 float3 f0 = float4_to_float3(kernel_tex_fetch(__attributes_float3, desc.offset + k0));
88                 float3 f1 = float4_to_float3(kernel_tex_fetch(__attributes_float3, desc.offset + k1));
89
90 #ifdef __RAY_DIFFERENTIALS__
91                 if(dx) *dx = sd->du.dx*(f1 - f0);
92                 if(dy) *dy = make_float3(0.0f, 0.0f, 0.0f);
93 #endif
94
95                 return (1.0f - sd->u)*f0 + sd->u*f1;
96         }
97         else {
98 #ifdef __RAY_DIFFERENTIALS__
99                 if(dx) *dx = make_float3(0.0f, 0.0f, 0.0f);
100                 if(dy) *dy = make_float3(0.0f, 0.0f, 0.0f);
101 #endif
102
103                 return make_float3(0.0f, 0.0f, 0.0f);
104         }
105 }
106
107 /* Curve thickness */
108
109 ccl_device float curve_thickness(KernelGlobals *kg, ShaderData *sd)
110 {
111         float r = 0.0f;
112
113         if(sd->type & PRIMITIVE_ALL_CURVE) {
114                 float4 curvedata = kernel_tex_fetch(__curves, sd->prim);
115                 int k0 = __float_as_int(curvedata.x) + PRIMITIVE_UNPACK_SEGMENT(sd->type);
116                 int k1 = k0 + 1;
117
118                 float4 P_curve[2];
119
120                 if(sd->type & PRIMITIVE_CURVE) {
121                         P_curve[0]= kernel_tex_fetch(__curve_keys, k0);
122                         P_curve[1]= kernel_tex_fetch(__curve_keys, k1);
123                 }
124                 else {
125                         motion_curve_keys(kg, sd->object, sd->prim, sd->time, k0, k1, P_curve);
126                 }
127
128                 r = (P_curve[1].w - P_curve[0].w) * sd->u + P_curve[0].w;
129         }
130
131         return r*2.0f;
132 }
133
134 /* Curve location for motion pass, linear interpolation between keys and
135  * ignoring radius because we do the same for the motion keys */
136
137 ccl_device float3 curve_motion_center_location(KernelGlobals *kg, ShaderData *sd)
138 {
139         float4 curvedata = kernel_tex_fetch(__curves, sd->prim);
140         int k0 = __float_as_int(curvedata.x) + PRIMITIVE_UNPACK_SEGMENT(sd->type);
141         int k1 = k0 + 1;
142
143         float4 P_curve[2];
144
145         P_curve[0]= kernel_tex_fetch(__curve_keys, k0);
146         P_curve[1]= kernel_tex_fetch(__curve_keys, k1);
147
148         return float4_to_float3(P_curve[1]) * sd->u + float4_to_float3(P_curve[0]) * (1.0f - sd->u);
149 }
150
151 /* Curve tangent normal */
152
153 ccl_device float3 curve_tangent_normal(KernelGlobals *kg, ShaderData *sd)
154 {       
155         float3 tgN = make_float3(0.0f,0.0f,0.0f);
156
157         if(sd->type & PRIMITIVE_ALL_CURVE) {
158
159                 tgN = -(-sd->I - sd->dPdu * (dot(sd->dPdu,-sd->I) / len_squared(sd->dPdu)));
160                 tgN = normalize(tgN);
161
162                 /* need to find suitable scaled gd for corrected normal */
163 #if 0
164                 tgN = normalize(tgN - gd * sd->dPdu);
165 #endif
166         }
167
168         return tgN;
169 }
170
171 /* Curve bounds utility function */
172
173 ccl_device_inline void curvebounds(float *lower, float *upper, float *extremta, float *extrema, float *extremtb, float *extremb, float p0, float p1, float p2, float p3)
174 {
175         float halfdiscroot = (p2 * p2 - 3 * p3 * p1);
176         float ta = -1.0f;
177         float tb = -1.0f;
178
179         *extremta = -1.0f;
180         *extremtb = -1.0f;
181         *upper = p0;
182         *lower = (p0 + p1) + (p2 + p3);
183         *extrema = *upper;
184         *extremb = *lower;
185
186         if(*lower >= *upper) {
187                 *upper = *lower;
188                 *lower = p0;
189         }
190
191         if(halfdiscroot >= 0) {
192                 float inv3p3 = (1.0f/3.0f)/p3;
193                 halfdiscroot = sqrtf(halfdiscroot);
194                 ta = (-p2 - halfdiscroot) * inv3p3;
195                 tb = (-p2 + halfdiscroot) * inv3p3;
196         }
197
198         float t2;
199         float t3;
200
201         if(ta > 0.0f && ta < 1.0f) {
202                 t2 = ta * ta;
203                 t3 = t2 * ta;
204                 *extremta = ta;
205                 *extrema = p3 * t3 + p2 * t2 + p1 * ta + p0;
206
207                 *upper = fmaxf(*extrema, *upper);
208                 *lower = fminf(*extrema, *lower);
209         }
210
211         if(tb > 0.0f && tb < 1.0f) {
212                 t2 = tb * tb;
213                 t3 = t2 * tb;
214                 *extremtb = tb;
215                 *extremb = p3 * t3 + p2 * t2 + p1 * tb + p0;
216
217                 *upper = fmaxf(*extremb, *upper);
218                 *lower = fminf(*extremb, *lower);
219         }
220 }
221
222 #ifdef __KERNEL_SSE2__
223 ccl_device_inline ssef transform_point_T3(const ssef t[3], const ssef &a)
224 {
225         return madd(shuffle<0>(a), t[0], madd(shuffle<1>(a), t[1], shuffle<2>(a) * t[2]));
226 }
227 #endif
228
229 #ifdef __KERNEL_SSE2__
230 /* Pass P and dir by reference to aligned vector */
231 ccl_device_curveintersect bool bvh_cardinal_curve_intersect(KernelGlobals *kg, Intersection *isect,
232         const float3 &P, const float3 &dir, uint visibility, int object, int curveAddr, float time, int type, uint *lcg_state, float difl, float extmax)
233 #else
234 ccl_device_curveintersect bool bvh_cardinal_curve_intersect(KernelGlobals *kg, Intersection *isect,
235         float3 P, float3 dir, uint visibility, int object, int curveAddr, float time,int type, uint *lcg_state, float difl, float extmax)
236 #endif
237 {
238         const bool is_curve_primitive = (type & PRIMITIVE_CURVE);
239
240         if(!is_curve_primitive && kernel_data.bvh.use_bvh_steps) {
241                 const float2 prim_time = kernel_tex_fetch(__prim_time, curveAddr);
242                 if(time < prim_time.x || time > prim_time.y) {
243                         return false;
244                 }
245         }
246
247         int segment = PRIMITIVE_UNPACK_SEGMENT(type);
248         float epsilon = 0.0f;
249         float r_st, r_en;
250
251         int depth = kernel_data.curve.subdivisions;
252         int flags = kernel_data.curve.curveflags;
253         int prim = kernel_tex_fetch(__prim_index, curveAddr);
254
255 #ifdef __KERNEL_SSE2__
256         ssef vdir = load4f(dir);
257         ssef vcurve_coef[4];
258         const float3 *curve_coef = (float3 *)vcurve_coef;
259         
260         {
261                 ssef dtmp = vdir * vdir;
262                 ssef d_ss = mm_sqrt(dtmp + shuffle<2>(dtmp));
263                 ssef rd_ss = load1f_first(1.0f) / d_ss;
264
265                 ssei v00vec = load4i((ssei *)&kg->__curves.data[prim]);
266                 int2 &v00 = (int2 &)v00vec;
267
268                 int k0 = v00.x + segment;
269                 int k1 = k0 + 1;
270                 int ka = max(k0 - 1, v00.x);
271                 int kb = min(k1 + 1, v00.x + v00.y - 1);
272
273 #if defined(__KERNEL_AVX2__) && defined(__KERNEL_SSE__) && (!defined(_MSC_VER) || _MSC_VER > 1800)
274                 avxf P_curve_0_1, P_curve_2_3;
275                 if(is_curve_primitive) {
276                         P_curve_0_1 = _mm256_loadu2_m128(&kg->__curve_keys.data[k0].x, &kg->__curve_keys.data[ka].x);
277                         P_curve_2_3 = _mm256_loadu2_m128(&kg->__curve_keys.data[kb].x, &kg->__curve_keys.data[k1].x);
278                 }
279                 else {
280                         int fobject = (object == OBJECT_NONE) ? kernel_tex_fetch(__prim_object, curveAddr) : object;
281                         motion_cardinal_curve_keys_avx(kg, fobject, prim, time, ka, k0, k1, kb, &P_curve_0_1,&P_curve_2_3);
282                 }
283 #else  /* __KERNEL_AVX2__ */
284                 ssef P_curve[4];
285
286                 if(is_curve_primitive) {
287                         P_curve[0] = load4f(&kg->__curve_keys.data[ka].x);
288                         P_curve[1] = load4f(&kg->__curve_keys.data[k0].x);
289                         P_curve[2] = load4f(&kg->__curve_keys.data[k1].x);
290                         P_curve[3] = load4f(&kg->__curve_keys.data[kb].x);
291                 }
292                 else {
293                         int fobject = (object == OBJECT_NONE)? kernel_tex_fetch(__prim_object, curveAddr): object;
294                         motion_cardinal_curve_keys(kg, fobject, prim, time, ka, k0, k1, kb, (float4*)&P_curve);
295                 }
296 #endif  /* __KERNEL_AVX2__ */
297
298                 ssef rd_sgn = set_sign_bit<0, 1, 1, 1>(shuffle<0>(rd_ss));
299                 ssef mul_zxxy = shuffle<2, 0, 0, 1>(vdir) * rd_sgn;
300                 ssef mul_yz = shuffle<1, 2, 1, 2>(vdir) * mul_zxxy;
301                 ssef mul_shuf = shuffle<0, 1, 2, 3>(mul_zxxy, mul_yz);
302                 ssef vdir0 = vdir & cast(ssei(0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0));
303
304                 ssef htfm0 = shuffle<0, 2, 0, 3>(mul_shuf, vdir0);
305                 ssef htfm1 = shuffle<1, 0, 1, 3>(load1f_first(extract<0>(d_ss)), vdir0);
306                 ssef htfm2 = shuffle<1, 3, 2, 3>(mul_shuf, vdir0);
307
308 #if defined(__KERNEL_AVX2__) && defined(__KERNEL_SSE__) && (!defined(_MSC_VER) || _MSC_VER > 1800)
309                 const avxf vPP = _mm256_broadcast_ps(&P.m128);
310                 const avxf htfm00 = avxf(htfm0.m128, htfm0.m128);
311                 const avxf htfm11 = avxf(htfm1.m128, htfm1.m128);
312                 const avxf htfm22 = avxf(htfm2.m128, htfm2.m128);
313
314                 const avxf p01 = madd(shuffle<0>(P_curve_0_1 - vPP),
315                                       htfm00,
316                                       madd(shuffle<1>(P_curve_0_1 - vPP),
317                                            htfm11,
318                                            shuffle<2>(P_curve_0_1 - vPP) * htfm22));
319                 const avxf p23 = madd(shuffle<0>(P_curve_2_3 - vPP),
320                                       htfm00,
321                                       madd(shuffle<1>(P_curve_2_3 - vPP),
322                                            htfm11,
323                                            shuffle<2>(P_curve_2_3 - vPP)*htfm22));
324
325                 const ssef p0 = _mm256_castps256_ps128(p01);
326                 const ssef p1 = _mm256_extractf128_ps(p01, 1);
327                 const ssef p2 = _mm256_castps256_ps128(p23);
328                 const ssef p3 = _mm256_extractf128_ps(p23, 1);
329
330                 const ssef P_curve_1 = _mm256_extractf128_ps(P_curve_0_1, 1);
331                 r_st = ((float4 &)P_curve_1).w;
332                 const ssef P_curve_2 = _mm256_castps256_ps128(P_curve_2_3);
333                 r_en = ((float4 &)P_curve_2).w;
334 #else  /* __KERNEL_AVX2__ */
335                 ssef htfm[] = { htfm0, htfm1, htfm2 };
336                 ssef vP = load4f(P);
337                 ssef p0 = transform_point_T3(htfm, P_curve[0] - vP);
338                 ssef p1 = transform_point_T3(htfm, P_curve[1] - vP);
339                 ssef p2 = transform_point_T3(htfm, P_curve[2] - vP);
340                 ssef p3 = transform_point_T3(htfm, P_curve[3] - vP);
341
342                 r_st = ((float4 &)P_curve[1]).w;
343                 r_en = ((float4 &)P_curve[2]).w;
344 #endif  /* __KERNEL_AVX2__ */
345
346                 float fc = 0.71f;
347                 ssef vfc = ssef(fc);
348                 ssef vfcxp3 = vfc * p3;
349
350                 vcurve_coef[0] = p1;
351                 vcurve_coef[1] = vfc * (p2 - p0);
352                 vcurve_coef[2] = madd(ssef(fc * 2.0f), p0, madd(ssef(fc - 3.0f), p1, msub(ssef(3.0f - 2.0f * fc), p2, vfcxp3)));
353                 vcurve_coef[3] = msub(ssef(fc - 2.0f), p2 - p1, msub(vfc, p0, vfcxp3));
354
355         }
356 #else
357         float3 curve_coef[4];
358
359         /* curve Intersection check */
360         /* obtain curve parameters */
361         {
362                 /* ray transform created - this should be created at beginning of intersection loop */
363                 Transform htfm;
364                 float d = sqrtf(dir.x * dir.x + dir.z * dir.z);
365                 htfm = make_transform(
366                         dir.z / d, 0, -dir.x /d, 0,
367                         -dir.x * dir.y /d, d, -dir.y * dir.z /d, 0,
368                         dir.x, dir.y, dir.z, 0,
369                         0, 0, 0, 1);
370
371                 float4 v00 = kernel_tex_fetch(__curves, prim);
372
373                 int k0 = __float_as_int(v00.x) + segment;
374                 int k1 = k0 + 1;
375
376                 int ka = max(k0 - 1,__float_as_int(v00.x));
377                 int kb = min(k1 + 1,__float_as_int(v00.x) + __float_as_int(v00.y) - 1);
378
379                 float4 P_curve[4];
380
381                 if(is_curve_primitive) {
382                         P_curve[0] = kernel_tex_fetch(__curve_keys, ka);
383                         P_curve[1] = kernel_tex_fetch(__curve_keys, k0);
384                         P_curve[2] = kernel_tex_fetch(__curve_keys, k1);
385                         P_curve[3] = kernel_tex_fetch(__curve_keys, kb);
386                 }
387                 else {
388                         int fobject = (object == OBJECT_NONE)? kernel_tex_fetch(__prim_object, curveAddr): object;
389                         motion_cardinal_curve_keys(kg, fobject, prim, time, ka, k0, k1, kb, P_curve);
390                 }
391
392                 float3 p0 = transform_point(&htfm, float4_to_float3(P_curve[0]) - P);
393                 float3 p1 = transform_point(&htfm, float4_to_float3(P_curve[1]) - P);
394                 float3 p2 = transform_point(&htfm, float4_to_float3(P_curve[2]) - P);
395                 float3 p3 = transform_point(&htfm, float4_to_float3(P_curve[3]) - P);
396
397                 float fc = 0.71f;
398                 curve_coef[0] = p1;
399                 curve_coef[1] = -fc*p0 + fc*p2;
400                 curve_coef[2] = 2.0f * fc * p0 + (fc - 3.0f) * p1 + (3.0f - 2.0f * fc) * p2 - fc * p3;
401                 curve_coef[3] = -fc * p0 + (2.0f - fc) * p1 + (fc - 2.0f) * p2 + fc * p3;
402                 r_st = P_curve[1].w;
403                 r_en = P_curve[2].w;
404         }
405 #endif
406
407         float r_curr = max(r_st, r_en);
408
409         if((flags & CURVE_KN_RIBBONS) || !(flags & CURVE_KN_BACKFACING))
410                 epsilon = 2 * r_curr;
411
412         /* find bounds - this is slow for cubic curves */
413         float upper, lower;
414
415         float zextrem[4];
416         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);
417         if(lower - r_curr > isect->t || upper + r_curr < epsilon)
418                 return false;
419
420         /* minimum width extension */
421         float mw_extension = min(difl * fabsf(upper), extmax);
422         float r_ext = mw_extension + r_curr;
423
424         float xextrem[4];
425         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);
426         if(lower > r_ext || upper < -r_ext)
427                 return false;
428
429         float yextrem[4];
430         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);
431         if(lower > r_ext || upper < -r_ext)
432                 return false;
433
434         /* setup recurrent loop */
435         int level = 1 << depth;
436         int tree = 0;
437         float resol = 1.0f / (float)level;
438         bool hit = false;
439
440         /* begin loop */
441         while(!(tree >> (depth))) {
442                 const float i_st = tree * resol;
443                 const float i_en = i_st + (level * resol);
444
445 #ifdef __KERNEL_SSE2__
446                 ssef vi_st = ssef(i_st), vi_en = ssef(i_en);
447                 ssef vp_st = madd(madd(madd(vcurve_coef[3], vi_st, vcurve_coef[2]), vi_st, vcurve_coef[1]), vi_st, vcurve_coef[0]);
448                 ssef vp_en = madd(madd(madd(vcurve_coef[3], vi_en, vcurve_coef[2]), vi_en, vcurve_coef[1]), vi_en, vcurve_coef[0]);
449
450                 ssef vbmin = min(vp_st, vp_en);
451                 ssef vbmax = max(vp_st, vp_en);
452
453                 float3 &bmin = (float3 &)vbmin, &bmax = (float3 &)vbmax;
454                 float &bminx = bmin.x, &bminy = bmin.y, &bminz = bmin.z;
455                 float &bmaxx = bmax.x, &bmaxy = bmax.y, &bmaxz = bmax.z;
456                 float3 &p_st = (float3 &)vp_st, &p_en = (float3 &)vp_en;
457 #else
458                 float3 p_st = ((curve_coef[3] * i_st + curve_coef[2]) * i_st + curve_coef[1]) * i_st + curve_coef[0];
459                 float3 p_en = ((curve_coef[3] * i_en + curve_coef[2]) * i_en + curve_coef[1]) * i_en + curve_coef[0];
460                 
461                 float bminx = min(p_st.x, p_en.x);
462                 float bmaxx = max(p_st.x, p_en.x);
463                 float bminy = min(p_st.y, p_en.y);
464                 float bmaxy = max(p_st.y, p_en.y);
465                 float bminz = min(p_st.z, p_en.z);
466                 float bmaxz = max(p_st.z, p_en.z);
467 #endif
468
469                 if(xextrem[0] >= i_st && xextrem[0] <= i_en) {
470                         bminx = min(bminx,xextrem[1]);
471                         bmaxx = max(bmaxx,xextrem[1]);
472                 }
473                 if(xextrem[2] >= i_st && xextrem[2] <= i_en) {
474                         bminx = min(bminx,xextrem[3]);
475                         bmaxx = max(bmaxx,xextrem[3]);
476                 }
477                 if(yextrem[0] >= i_st && yextrem[0] <= i_en) {
478                         bminy = min(bminy,yextrem[1]);
479                         bmaxy = max(bmaxy,yextrem[1]);
480                 }
481                 if(yextrem[2] >= i_st && yextrem[2] <= i_en) {
482                         bminy = min(bminy,yextrem[3]);
483                         bmaxy = max(bmaxy,yextrem[3]);
484                 }
485                 if(zextrem[0] >= i_st && zextrem[0] <= i_en) {
486                         bminz = min(bminz,zextrem[1]);
487                         bmaxz = max(bmaxz,zextrem[1]);
488                 }
489                 if(zextrem[2] >= i_st && zextrem[2] <= i_en) {
490                         bminz = min(bminz,zextrem[3]);
491                         bmaxz = max(bmaxz,zextrem[3]);
492                 }
493
494                 float r1 = r_st + (r_en - r_st) * i_st;
495                 float r2 = r_st + (r_en - r_st) * i_en;
496                 r_curr = max(r1, r2);
497
498                 mw_extension = min(difl * fabsf(bmaxz), extmax);
499                 float r_ext = mw_extension + r_curr;
500                 float coverage = 1.0f;
501
502                 if(bminz - r_curr > isect->t || bmaxz + r_curr < epsilon || bminx > r_ext|| bmaxx < -r_ext|| bminy > r_ext|| bmaxy < -r_ext) {
503                         /* the bounding box does not overlap the square centered at O */
504                         tree += level;
505                         level = tree & -tree;
506                 }
507                 else if(level == 1) {
508
509                         /* the maximum recursion depth is reached.
510                          * check if dP0.(Q-P0)>=0 and dPn.(Pn-Q)>=0.
511                          * dP* is reversed if necessary.*/
512                         float t = isect->t;
513                         float u = 0.0f;
514                         float gd = 0.0f;
515
516                         if(flags & CURVE_KN_RIBBONS) {
517                                 float3 tg = (p_en - p_st);
518 #ifdef __KERNEL_SSE__
519                                 const float3 tg_sq = tg * tg;
520                                 float w = tg_sq.x + tg_sq.y;
521 #else
522                                 float w = tg.x * tg.x + tg.y * tg.y;
523 #endif
524                                 if(w == 0) {
525                                         tree++;
526                                         level = tree & -tree;
527                                         continue;
528                                 }
529 #ifdef __KERNEL_SSE__
530                                 const float3 p_sttg = p_st * tg;
531                                 w = -(p_sttg.x + p_sttg.y) / w;
532 #else
533                                 w = -(p_st.x * tg.x + p_st.y * tg.y) / w;
534 #endif
535                                 w = saturate(w);
536
537                                 /* compute u on the curve segment */
538                                 u = i_st * (1 - w) + i_en * w;
539                                 r_curr = r_st + (r_en - r_st) * u;
540                                 /* compare x-y distances */
541                                 float3 p_curr = ((curve_coef[3] * u + curve_coef[2]) * u + curve_coef[1]) * u + curve_coef[0];
542
543                                 float3 dp_st = (3 * curve_coef[3] * i_st + 2 * curve_coef[2]) * i_st + curve_coef[1];
544                                 if(dot(tg, dp_st)< 0)
545                                         dp_st *= -1;
546                                 if(dot(dp_st, -p_st) + p_curr.z * dp_st.z < 0) {
547                                         tree++;
548                                         level = tree & -tree;
549                                         continue;
550                                 }
551                                 float3 dp_en = (3 * curve_coef[3] * i_en + 2 * curve_coef[2]) * i_en + curve_coef[1];
552                                 if(dot(tg, dp_en) < 0)
553                                         dp_en *= -1;
554                                 if(dot(dp_en, p_en) - p_curr.z * dp_en.z < 0) {
555                                         tree++;
556                                         level = tree & -tree;
557                                         continue;
558                                 }
559
560                                 /* compute coverage */
561                                 float r_ext = r_curr;
562                                 coverage = 1.0f;
563                                 if(difl != 0.0f) {
564                                         mw_extension = min(difl * fabsf(bmaxz), extmax);
565                                         r_ext = mw_extension + r_curr;
566 #ifdef __KERNEL_SSE__
567                                         const float3 p_curr_sq = p_curr * p_curr;
568                                         const float3 dxxx = _mm_sqrt_ss(_mm_hadd_ps(p_curr_sq.m128, p_curr_sq.m128));
569                                         float d = dxxx.x;
570 #else
571                                         float d = sqrtf(p_curr.x * p_curr.x + p_curr.y * p_curr.y);
572 #endif
573                                         float d0 = d - r_curr;
574                                         float d1 = d + r_curr;
575                                         float inv_mw_extension = 1.0f/mw_extension;
576                                         if(d0 >= 0)
577                                                 coverage = (min(d1 * inv_mw_extension, 1.0f) - min(d0 * inv_mw_extension, 1.0f)) * 0.5f;
578                                         else // inside
579                                                 coverage = (min(d1 * inv_mw_extension, 1.0f) + min(-d0 * inv_mw_extension, 1.0f)) * 0.5f;
580                                 }
581                                 
582                                 if(p_curr.x * p_curr.x + p_curr.y * p_curr.y >= r_ext * r_ext || p_curr.z <= epsilon || isect->t < p_curr.z) {
583                                         tree++;
584                                         level = tree & -tree;
585                                         continue;
586                                 }
587
588                                 t = p_curr.z;
589
590                                 /* stochastic fade from minimum width */
591                                 if(difl != 0.0f && lcg_state) {
592                                         if(coverage != 1.0f && (lcg_step_float(lcg_state) > coverage))
593                                                 return hit;
594                                 }
595                         }
596                         else {
597                                 float l = len(p_en - p_st);
598                                 /* minimum width extension */
599                                 float or1 = r1;
600                                 float or2 = r2;
601
602                                 if(difl != 0.0f) {
603                                         mw_extension = min(len(p_st - P) * difl, extmax);
604                                         or1 = r1 < mw_extension ? mw_extension : r1;
605                                         mw_extension = min(len(p_en - P) * difl, extmax);
606                                         or2 = r2 < mw_extension ? mw_extension : r2;
607                                 }
608                                 /* --- */
609                                 float invl = 1.0f/l;
610                                 float3 tg = (p_en - p_st) * invl;
611                                 gd = (or2 - or1) * invl;
612                                 float difz = -dot(p_st,tg);
613                                 float cyla = 1.0f - (tg.z * tg.z * (1 + gd*gd));
614                                 float invcyla = 1.0f/cyla;
615                                 float halfb = (-p_st.z - tg.z*(difz + gd*(difz*gd + or1)));
616                                 float tcentre = -halfb*invcyla;
617                                 float zcentre = difz + (tg.z * tcentre);
618                                 float3 tdif = - p_st;
619                                 tdif.z += tcentre;
620                                 float tdifz = dot(tdif,tg);
621                                 float tb = 2*(tdif.z - tg.z*(tdifz + gd*(tdifz*gd + or1)));
622                                 float tc = dot(tdif,tdif) - tdifz * tdifz * (1 + gd*gd) - or1*or1 - 2*or1*tdifz*gd;
623                                 float td = tb*tb - 4*cyla*tc;
624                                 if(td < 0.0f) {
625                                         tree++;
626                                         level = tree & -tree;
627                                         continue;
628                                 }
629                                 
630                                 float rootd = sqrtf(td);
631                                 float correction = (-tb - rootd) * 0.5f * invcyla;
632                                 t = tcentre + correction;
633
634                                 float3 dp_st = (3 * curve_coef[3] * i_st + 2 * curve_coef[2]) * i_st + curve_coef[1];
635                                 if(dot(tg, dp_st)< 0)
636                                         dp_st *= -1;
637                                 float3 dp_en = (3 * curve_coef[3] * i_en + 2 * curve_coef[2]) * i_en + curve_coef[1];
638                                 if(dot(tg, dp_en) < 0)
639                                         dp_en *= -1;
640
641                                 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)) {
642                                         correction = (-tb + rootd) * 0.5f * invcyla;
643                                         t = tcentre + correction;
644                                 }                       
645
646                                 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) {
647                                         tree++;
648                                         level = tree & -tree;
649                                         continue;
650                                 }
651
652                                 float w = (zcentre + (tg.z * correction)) * invl;
653                                 w = saturate(w);
654                                 /* compute u on the curve segment */
655                                 u = i_st * (1 - w) + i_en * w;
656
657                                 /* stochastic fade from minimum width */
658                                 if(difl != 0.0f && lcg_state) {
659                                         r_curr = r1 + (r2 - r1) * w;
660                                         r_ext = or1 + (or2 - or1) * w;
661                                         coverage = r_curr/r_ext;
662
663                                         if(coverage != 1.0f && (lcg_step_float(lcg_state) > coverage))
664                                                 return hit;
665                                 }
666                         }
667                         /* we found a new intersection */
668
669 #ifdef __VISIBILITY_FLAG__
670                         /* visibility flag test. we do it here under the assumption
671                          * that most triangles are culled by node flags */
672                         if(kernel_tex_fetch(__prim_visibility, curveAddr) & visibility)
673 #endif
674                         {
675                                 /* record intersection */
676                                 isect->t = t;
677                                 isect->u = u;
678                                 isect->v = gd;
679                                 isect->prim = curveAddr;
680                                 isect->object = object;
681                                 isect->type = type;
682                                 hit = true;
683                         }
684                         
685                         tree++;
686                         level = tree & -tree;
687                 }
688                 else {
689                         /* split the curve into two curves and process */
690                         level = level >> 1;
691                 }
692         }
693
694         return hit;
695 }
696
697 ccl_device_curveintersect bool bvh_curve_intersect(KernelGlobals *kg, Intersection *isect,
698         float3 P, float3 direction, uint visibility, int object, int curveAddr, float time, int type, uint *lcg_state, float difl, float extmax)
699 {
700         /* define few macros to minimize code duplication for SSE */
701 #ifndef __KERNEL_SSE2__
702 #  define len3_squared(x) len_squared(x)
703 #  define len3(x) len(x)
704 #  define dot3(x, y) dot(x, y)
705 #endif
706
707         const bool is_curve_primitive = (type & PRIMITIVE_CURVE);
708
709         if(!is_curve_primitive && kernel_data.bvh.use_bvh_steps) {
710                 const float2 prim_time = kernel_tex_fetch(__prim_time, curveAddr);
711                 if(time < prim_time.x || time > prim_time.y) {
712                         return false;
713                 }
714         }
715
716         int segment = PRIMITIVE_UNPACK_SEGMENT(type);
717         /* curve Intersection check */
718         int flags = kernel_data.curve.curveflags;
719
720         int prim = kernel_tex_fetch(__prim_index, curveAddr);
721         float4 v00 = kernel_tex_fetch(__curves, prim);
722
723         int cnum = __float_as_int(v00.x);
724         int k0 = cnum + segment;
725         int k1 = k0 + 1;
726
727 #ifndef __KERNEL_SSE2__
728         float4 P_curve[2];
729
730         if(is_curve_primitive) {
731                 P_curve[0] = kernel_tex_fetch(__curve_keys, k0);
732                 P_curve[1] = kernel_tex_fetch(__curve_keys, k1);
733         }
734         else {
735                 int fobject = (object == OBJECT_NONE)? kernel_tex_fetch(__prim_object, curveAddr): object;
736                 motion_curve_keys(kg, fobject, prim, time, k0, k1, P_curve);
737         }
738
739         float or1 = P_curve[0].w;
740         float or2 = P_curve[1].w;
741         float3 p1 = float4_to_float3(P_curve[0]);
742         float3 p2 = float4_to_float3(P_curve[1]);
743
744         /* minimum width extension */
745         float r1 = or1;
746         float r2 = or2;
747         float3 dif = P - p1;
748         float3 dif_second = P - p2;
749         if(difl != 0.0f) {
750                 float pixelsize = min(len3(dif) * difl, extmax);
751                 r1 = or1 < pixelsize ? pixelsize : or1;
752                 pixelsize = min(len3(dif_second) * difl, extmax);
753                 r2 = or2 < pixelsize ? pixelsize : or2;
754         }
755         /* --- */
756
757         float3 p21_diff = p2 - p1;
758         float3 sphere_dif1 = (dif + dif_second) * 0.5f;
759         float3 dir = direction;
760         float sphere_b_tmp = dot3(dir, sphere_dif1);
761         float3 sphere_dif2 = sphere_dif1 - sphere_b_tmp * dir;
762 #else
763         ssef P_curve[2];
764         
765         if(is_curve_primitive) {
766                 P_curve[0] = load4f(&kg->__curve_keys.data[k0].x);
767                 P_curve[1] = load4f(&kg->__curve_keys.data[k1].x);
768         }
769         else {
770                 int fobject = (object == OBJECT_NONE)? kernel_tex_fetch(__prim_object, curveAddr): object;
771                 motion_curve_keys(kg, fobject, prim, time, k0, k1, (float4*)&P_curve);
772         }
773
774         const ssef or12 = shuffle<3, 3, 3, 3>(P_curve[0], P_curve[1]);
775
776         ssef r12 = or12;
777         const ssef vP = load4f(P);
778         const ssef dif = vP - P_curve[0];
779         const ssef dif_second = vP - P_curve[1];
780         if(difl != 0.0f) {
781                 const ssef len1_sq = len3_squared_splat(dif);
782                 const ssef len2_sq = len3_squared_splat(dif_second);
783                 const ssef len12 = mm_sqrt(shuffle<0, 0, 0, 0>(len1_sq, len2_sq));
784                 const ssef pixelsize12 = min(len12 * difl, ssef(extmax));
785                 r12 = max(or12, pixelsize12);
786         }
787         float or1 = extract<0>(or12), or2 = extract<0>(shuffle<2>(or12));
788         float r1 = extract<0>(r12), r2 = extract<0>(shuffle<2>(r12));
789
790         const ssef p21_diff = P_curve[1] - P_curve[0];
791         const ssef sphere_dif1 = (dif + dif_second) * 0.5f;
792         const ssef dir = load4f(direction);
793         const ssef sphere_b_tmp = dot3_splat(dir, sphere_dif1);
794         const ssef sphere_dif2 = nmadd(sphere_b_tmp, dir, sphere_dif1);
795 #endif
796
797         float mr = max(r1, r2);
798         float l = len3(p21_diff);
799         float invl = 1.0f / l;
800         float sp_r = mr + 0.5f * l;
801
802         float sphere_b = dot3(dir, sphere_dif2);
803         float sdisc = sphere_b * sphere_b - len3_squared(sphere_dif2) + sp_r * sp_r;
804
805         if(sdisc < 0.0f)
806                 return false;
807
808         /* obtain parameters and test midpoint distance for suitable modes */
809 #ifndef __KERNEL_SSE2__
810         float3 tg = p21_diff * invl;
811 #else
812         const ssef tg = p21_diff * invl;
813 #endif
814         float gd = (r2 - r1) * invl;
815
816         float dirz = dot3(dir, tg);
817         float difz = dot3(dif, tg);
818
819         float a = 1.0f - (dirz*dirz*(1 + gd*gd));
820
821         float halfb = dot3(dir, dif) - dirz*(difz + gd*(difz*gd + r1));
822
823         float tcentre = -halfb/a;
824         float zcentre = difz + (dirz * tcentre);
825
826         if((tcentre > isect->t) && !(flags & CURVE_KN_ACCURATE))
827                 return false;
828         if((zcentre < 0 || zcentre > l) && !(flags & CURVE_KN_ACCURATE) && !(flags & CURVE_KN_INTERSECTCORRECTION))
829                 return false;
830
831         /* test minimum separation */
832 #ifndef __KERNEL_SSE2__
833         float3 cprod = cross(tg, dir);
834         float cprod2sq = len3_squared(cross(tg, dif));
835 #else
836         const ssef cprod = cross(tg, dir);
837         float cprod2sq = len3_squared(cross_zxy(tg, dif));
838 #endif
839         float cprodsq = len3_squared(cprod);
840         float distscaled = dot3(cprod, dif);
841
842         if(cprodsq == 0)
843                 distscaled = cprod2sq;
844         else
845                 distscaled = (distscaled*distscaled)/cprodsq;
846
847         if(distscaled > mr*mr)
848                 return false;
849
850         /* calculate true intersection */
851 #ifndef __KERNEL_SSE2__
852         float3 tdif = dif + tcentre * dir;
853 #else
854         const ssef tdif = madd(ssef(tcentre), dir, dif);
855 #endif
856         float tdifz = dot3(tdif, tg);
857         float tdifma = tdifz*gd + r1;
858         float tb = 2*(dot3(dir, tdif) - dirz*(tdifz + gd*tdifma));
859         float tc = dot3(tdif, tdif) - tdifz*tdifz - tdifma*tdifma;
860         float td = tb*tb - 4*a*tc;
861
862         if(td < 0.0f)
863                 return false;
864
865         float rootd = 0.0f;
866         float correction = 0.0f;
867         if(flags & CURVE_KN_ACCURATE) {
868                 rootd = sqrtf(td);
869                 correction = ((-tb - rootd)/(2*a));
870         }
871
872         float t = tcentre + correction;
873
874         if(t < isect->t) {
875
876                 if(flags & CURVE_KN_INTERSECTCORRECTION) {
877                         rootd = sqrtf(td);
878                         correction = ((-tb - rootd)/(2*a));
879                         t = tcentre + correction;
880                 }
881
882                 float z = zcentre + (dirz * correction);
883                 // bool backface = false;
884
885                 if(flags & CURVE_KN_BACKFACING && (t < 0.0f || z < 0 || z > l)) {
886                         // backface = true;
887                         correction = ((-tb + rootd)/(2*a));
888                         t = tcentre + correction;
889                         z = zcentre + (dirz * correction);
890                 }
891
892                 /* stochastic fade from minimum width */
893                 float adjradius = or1 + z * (or2 - or1) * invl;
894                 adjradius = adjradius / (r1 + z * gd);
895                 if(lcg_state && adjradius != 1.0f) {
896                         if(lcg_step_float(lcg_state) > adjradius)
897                                 return false;
898                 }
899                 /* --- */
900
901                 if(t > 0.0f && t < isect->t && z >= 0 && z <= l) {
902
903                         if(flags & CURVE_KN_ENCLOSEFILTER) {
904                                 float enc_ratio = 1.01f;
905                                 if((difz > -r1 * enc_ratio) && (dot3(dif_second, tg) < r2 * enc_ratio)) {
906                                         float a2 = 1.0f - (dirz*dirz*(1 + gd*gd*enc_ratio*enc_ratio));
907                                         float c2 = dot3(dif, dif) - difz * difz * (1 + gd*gd*enc_ratio*enc_ratio) - r1*r1*enc_ratio*enc_ratio - 2*r1*difz*gd*enc_ratio;
908                                         if(a2*c2 < 0.0f)
909                                                 return false;
910                                 }
911                         }
912
913 #ifdef __VISIBILITY_FLAG__
914                         /* visibility flag test. we do it here under the assumption
915                          * that most triangles are culled by node flags */
916                         if(kernel_tex_fetch(__prim_visibility, curveAddr) & visibility)
917 #endif
918                         {
919                                 /* record intersection */
920                                 isect->t = t;
921                                 isect->u = z*invl;
922                                 isect->v = gd;
923                                 isect->prim = curveAddr;
924                                 isect->object = object;
925                                 isect->type = type;
926
927                                 return true;
928                         }
929                 }
930         }
931
932         return false;
933
934 #ifndef __KERNEL_SSE2__
935 #  undef len3_squared
936 #  undef len3
937 #  undef dot3
938 #endif
939 }
940
941 ccl_device_inline float3 curvetangent(float t, float3 p0, float3 p1, float3 p2, float3 p3)
942 {
943         float fc = 0.71f;
944         float data[4];
945         float t2 = t * t;
946         data[0] = -3.0f * fc          * t2  + 4.0f * fc * t                  - fc;
947         data[1] =  3.0f * (2.0f - fc) * t2  + 2.0f * (fc - 3.0f) * t;
948         data[2] =  3.0f * (fc - 2.0f) * t2  + 2.0f * (3.0f - 2.0f * fc) * t  + fc;
949         data[3] =  3.0f * fc          * t2  - 2.0f * fc * t;
950         return data[0] * p0 + data[1] * p1 + data[2] * p2 + data[3] * p3;
951 }
952
953 ccl_device_inline float3 curvepoint(float t, float3 p0, float3 p1, float3 p2, float3 p3)
954 {
955         float data[4];
956         float fc = 0.71f;
957         float t2 = t * t;
958         float t3 = t2 * t;
959         data[0] = -fc          * t3  + 2.0f * fc          * t2 - fc * t;
960         data[1] =  (2.0f - fc) * t3  + (fc - 3.0f)        * t2 + 1.0f;
961         data[2] =  (fc - 2.0f) * t3  + (3.0f - 2.0f * fc) * t2 + fc * t;
962         data[3] =  fc          * t3  - fc * t2;
963         return data[0] * p0 + data[1] * p1 + data[2] * p2 + data[3] * p3;
964 }
965
966 ccl_device_inline float3 bvh_curve_refine(KernelGlobals *kg, ShaderData *sd, const Intersection *isect, const Ray *ray)
967 {
968         int flag = kernel_data.curve.curveflags;
969         float t = isect->t;
970         float3 P = ray->P;
971         float3 D = ray->D;
972
973         if(isect->object != OBJECT_NONE) {
974 #ifdef __OBJECT_MOTION__
975                 Transform tfm = sd->ob_itfm;
976 #else
977                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
978 #endif
979
980                 P = transform_point(&tfm, P);
981                 D = transform_direction(&tfm, D*t);
982                 D = normalize_len(D, &t);
983         }
984
985         int prim = kernel_tex_fetch(__prim_index, isect->prim);
986         float4 v00 = kernel_tex_fetch(__curves, prim);
987
988         int k0 = __float_as_int(v00.x) + PRIMITIVE_UNPACK_SEGMENT(sd->type);
989         int k1 = k0 + 1;
990
991         float3 tg;
992
993         if(flag & CURVE_KN_INTERPOLATE) {
994                 int ka = max(k0 - 1,__float_as_int(v00.x));
995                 int kb = min(k1 + 1,__float_as_int(v00.x) + __float_as_int(v00.y) - 1);
996
997                 float4 P_curve[4];
998
999                 if(sd->type & PRIMITIVE_CURVE) {
1000                         P_curve[0] = kernel_tex_fetch(__curve_keys, ka);
1001                         P_curve[1] = kernel_tex_fetch(__curve_keys, k0);
1002                         P_curve[2] = kernel_tex_fetch(__curve_keys, k1);
1003                         P_curve[3] = kernel_tex_fetch(__curve_keys, kb);
1004                 }
1005                 else {
1006                         motion_cardinal_curve_keys(kg, sd->object, sd->prim, sd->time, ka, k0, k1, kb, P_curve);
1007                 }
1008
1009                 float3 p[4];
1010                 p[0] = float4_to_float3(P_curve[0]);
1011                 p[1] = float4_to_float3(P_curve[1]);
1012                 p[2] = float4_to_float3(P_curve[2]);
1013                 p[3] = float4_to_float3(P_curve[3]);
1014
1015                 P = P + D*t;
1016
1017 #ifdef __UV__
1018                 sd->u = isect->u;
1019                 sd->v = 0.0f;
1020 #endif
1021
1022                 tg = normalize(curvetangent(isect->u, p[0], p[1], p[2], p[3]));
1023
1024                 if(kernel_data.curve.curveflags & CURVE_KN_RIBBONS) {
1025                         sd->Ng = normalize(-(D - tg * (dot(tg, D))));
1026                 }
1027                 else {
1028                         /* direction from inside to surface of curve */
1029                         float3 p_curr = curvepoint(isect->u, p[0], p[1], p[2], p[3]);   
1030                         sd->Ng = normalize(P - p_curr);
1031
1032                         /* adjustment for changing radius */
1033                         float gd = isect->v;
1034
1035                         if(gd != 0.0f) {
1036                                 sd->Ng = sd->Ng - gd * tg;
1037                                 sd->Ng = normalize(sd->Ng);
1038                         }
1039                 }
1040
1041                 /* todo: sometimes the normal is still so that this is detected as
1042                  * backfacing even if cull backfaces is enabled */
1043
1044                 sd->N = sd->Ng;
1045         }
1046         else {
1047                 float4 P_curve[2];
1048
1049                 if(sd->type & PRIMITIVE_CURVE) {
1050                         P_curve[0]= kernel_tex_fetch(__curve_keys, k0);
1051                         P_curve[1]= kernel_tex_fetch(__curve_keys, k1);
1052                 }
1053                 else {
1054                         motion_curve_keys(kg, sd->object, sd->prim, sd->time, k0, k1, P_curve);
1055                 }
1056
1057                 float l = 1.0f;
1058                 tg = normalize_len(float4_to_float3(P_curve[1] - P_curve[0]), &l);
1059                 
1060                 P = P + D*t;
1061
1062                 float3 dif = P - float4_to_float3(P_curve[0]);
1063
1064 #ifdef __UV__
1065                 sd->u = dot(dif,tg)/l;
1066                 sd->v = 0.0f;
1067 #endif
1068
1069                 if(flag & CURVE_KN_TRUETANGENTGNORMAL) {
1070                         sd->Ng = -(D - tg * dot(tg, D));
1071                         sd->Ng = normalize(sd->Ng);
1072                 }
1073                 else {
1074                         float gd = isect->v;
1075
1076                         /* direction from inside to surface of curve */
1077                         sd->Ng = (dif - tg * sd->u * l) / (P_curve[0].w + sd->u * l * gd);
1078
1079                         /* adjustment for changing radius */
1080                         if(gd != 0.0f) {
1081                                 sd->Ng = sd->Ng - gd * tg;
1082                                 sd->Ng = normalize(sd->Ng);
1083                         }
1084                 }
1085
1086                 sd->N = sd->Ng;
1087         }
1088
1089 #ifdef __DPDU__
1090         /* dPdu/dPdv */
1091         sd->dPdu = tg;
1092         sd->dPdv = cross(tg, sd->Ng);
1093 #endif
1094
1095         if(isect->object != OBJECT_NONE) {
1096 #ifdef __OBJECT_MOTION__
1097                 Transform tfm = sd->ob_tfm;
1098 #else
1099                 Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
1100 #endif
1101
1102                 P = transform_point(&tfm, P);
1103         }
1104
1105         return P;
1106 }
1107
1108 #endif
1109
1110 CCL_NAMESPACE_END
1111