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