intern/cycles/kernel/kernel_bvh.h

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