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