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