diff options
Diffstat (limited to 'intern/cycles/kernel/geom/geom_curve_intersect.h')
| -rw-r--r-- | intern/cycles/kernel/geom/geom_curve_intersect.h | 1354 |
1 files changed, 637 insertions, 717 deletions
diff --git a/intern/cycles/kernel/geom/geom_curve_intersect.h b/intern/cycles/kernel/geom/geom_curve_intersect.h index 7a770470150..c04dbee52cc 100644 --- a/intern/cycles/kernel/geom/geom_curve_intersect.h +++ b/intern/cycles/kernel/geom/geom_curve_intersect.h @@ -1,4 +1,7 @@ /* + * Copyright 2009-2020 Intel Corporation. Adapted from Embree with + * with modifications. + * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at @@ -14,501 +17,620 @@ CCL_NAMESPACE_BEGIN -/* Curve primitive intersection functions. */ +/* Curve primitive intersection functions. + * + * The code here was adapted from curve_intersector_sweep.h in Embree, to get + * an exact match between Embree CPU ray-tracing and our GPU ray-tracing. */ + +#define CURVE_NUM_BEZIER_SUBDIVISIONS 3 +#define CURVE_NUM_BEZIER_SUBDIVISIONS_UNSTABLE (CURVE_NUM_BEZIER_SUBDIVISIONS + 1) +#define CURVE_NUM_BEZIER_STEPS 2 +#define CURVE_NUM_JACOBIAN_ITERATIONS 5 #ifdef __HAIR__ -# ifdef __KERNEL_SSE2__ -ccl_device_inline ssef transform_point_T3(const ssef t[3], const ssef &a) +/* Catmull-rom curve evaluation. */ + +ccl_device_inline float4 catmull_rom_basis_eval(const float4 curve[4], float u) { - return madd(shuffle<0>(a), t[0], madd(shuffle<1>(a), t[1], shuffle<2>(a) * t[2])); + const float t = u; + const float s = 1.0f - u; + const float n0 = -t * s * s; + const float n1 = 2.0f + t * t * (3.0f * t - 5.0f); + const float n2 = 2.0f + s * s * (3.0f * s - 5.0f); + const float n3 = -s * t * t; + return 0.5f * (curve[0] * n0 + curve[1] * n1 + curve[2] * n2 + curve[3] * n3); } -# endif -/* On CPU pass P and dir by reference to aligned vector. */ -ccl_device_forceinline bool cardinal_curve_intersect(KernelGlobals *kg, - Intersection *isect, - const float3 ccl_ref P, - const float3 ccl_ref dir, - uint visibility, - int object, - int curveAddr, - float time, - int type) +ccl_device_inline float4 catmull_rom_basis_derivative(const float4 curve[4], float u) { - const bool is_curve_primitive = (type & PRIMITIVE_CURVE); + const float t = u; + const float s = 1.0f - u; + const float n0 = -s * s + 2.0f * s * t; + const float n1 = 2.0f * t * (3.0f * t - 5.0f) + 3.0f * t * t; + const float n2 = 2.0f * s * (3.0f * t + 2.0f) - 3.0f * s * s; + const float n3 = -2.0f * s * t + t * t; + return 0.5f * (curve[0] * n0 + curve[1] * n1 + curve[2] * n2 + curve[3] * n3); +} -# ifndef __KERNEL_OPTIX__ /* see OptiX motion flag OPTIX_MOTION_FLAG_[START|END]_VANISH */ - if (!is_curve_primitive && kernel_data.bvh.use_bvh_steps) { - const float2 prim_time = kernel_tex_fetch(__prim_time, curveAddr); - if (time < prim_time.x || time > prim_time.y) { - return false; - } - } -# endif +ccl_device_inline float4 catmull_rom_basis_derivative2(const float4 curve[4], float u) +{ - int segment = PRIMITIVE_UNPACK_SEGMENT(type); - float epsilon = 0.0f; - float r_st, r_en; + const float t = u; + const float n0 = -3.0f * t + 2.0f; + const float n1 = 9.0f * t - 5.0f; + const float n2 = -9.0f * t + 4.0f; + const float n3 = 3.0f * t - 1.0f; + return (curve[0] * n0 + curve[1] * n1 + curve[2] * n2 + curve[3] * n3); +} - int depth = kernel_data.curve.subdivisions; - int flags = kernel_data.curve.curveflags; - int prim = kernel_tex_fetch(__prim_index, curveAddr); +/* Thick Curve */ -# ifdef __KERNEL_SSE2__ - ssef vdir = load4f(dir); - ssef vcurve_coef[4]; - const float3 *curve_coef = (float3 *)vcurve_coef; +ccl_device_inline float3 dnormalize(const float3 p, const float3 dp) +{ + const float pp = dot(p, p); + const float pdp = dot(p, dp); + return (pp * dp - pdp * p) / (pp * sqrtf(pp)); +} - { - ssef dtmp = vdir * vdir; - ssef d_ss = mm_sqrt(dtmp + shuffle<2>(dtmp)); - ssef rd_ss = load1f_first(1.0f) / d_ss; - - ssei v00vec = load4i((ssei *)&kg->__curves.data[prim]); - int2 &v00 = (int2 &)v00vec; - - int k0 = v00.x + segment; - int k1 = k0 + 1; - int ka = max(k0 - 1, v00.x); - int kb = min(k1 + 1, v00.x + v00.y - 1); - -# if defined(__KERNEL_AVX2__) && defined(__KERNEL_SSE__) && \ - (!defined(_MSC_VER) || _MSC_VER > 1800) - avxf P_curve_0_1, P_curve_2_3; - if (is_curve_primitive) { - P_curve_0_1 = _mm256_loadu2_m128(&kg->__curve_keys.data[k0].x, &kg->__curve_keys.data[ka].x); - P_curve_2_3 = _mm256_loadu2_m128(&kg->__curve_keys.data[kb].x, &kg->__curve_keys.data[k1].x); - } - else { - int fobject = (object == OBJECT_NONE) ? kernel_tex_fetch(__prim_object, curveAddr) : object; - motion_cardinal_curve_keys_avx( - kg, fobject, prim, time, ka, k0, k1, kb, &P_curve_0_1, &P_curve_2_3); - } -# else /* __KERNEL_AVX2__ */ - ssef P_curve[4]; - - if (is_curve_primitive) { - P_curve[0] = load4f(&kg->__curve_keys.data[ka].x); - P_curve[1] = load4f(&kg->__curve_keys.data[k0].x); - P_curve[2] = load4f(&kg->__curve_keys.data[k1].x); - P_curve[3] = load4f(&kg->__curve_keys.data[kb].x); +ccl_device_inline float sqr_point_to_line_distance(const float3 PmQ0, const float3 Q1mQ0) +{ + const float3 N = cross(PmQ0, Q1mQ0); + const float3 D = Q1mQ0; + return dot(N, N) / dot(D, D); +} + +ccl_device_inline bool cylinder_intersect(const float3 cylinder_start, + const float3 cylinder_end, + const float cylinder_radius, + const float3 ray_dir, + float2 *t_o, + float *u0_o, + float3 *Ng0_o, + float *u1_o, + float3 *Ng1_o) +{ + /* Calculate quadratic equation to solve. */ + const float rl = 1.0f / len(cylinder_end - cylinder_start); + const float3 P0 = cylinder_start, dP = (cylinder_end - cylinder_start) * rl; + const float3 O = -P0, dO = ray_dir; + + const float dOdO = dot(dO, dO); + const float OdO = dot(dO, O); + const float OO = dot(O, O); + const float dOz = dot(dP, dO); + const float Oz = dot(dP, O); + + const float A = dOdO - sqr(dOz); + const float B = 2.0f * (OdO - dOz * Oz); + const float C = OO - sqr(Oz) - sqr(cylinder_radius); + + /* We miss the cylinder if determinant is smaller than zero. */ + const float D = B * B - 4.0f * A * C; + if (!(D >= 0.0f)) { + *t_o = make_float2(FLT_MAX, -FLT_MAX); + return false; + } + + /* Special case for rays that are parallel to the cylinder. */ + const float eps = 16.0f * FLT_EPSILON * max(fabsf(dOdO), fabsf(sqr(dOz))); + if (fabsf(A) < eps) { + if (C <= 0.0f) { + *t_o = make_float2(-FLT_MAX, FLT_MAX); + return true; } else { - int fobject = (object == OBJECT_NONE) ? kernel_tex_fetch(__prim_object, curveAddr) : object; - motion_cardinal_curve_keys(kg, fobject, prim, time, ka, k0, k1, kb, (float4 *)&P_curve); + *t_o = make_float2(-FLT_MAX, FLT_MAX); + return false; } -# endif /* __KERNEL_AVX2__ */ - - ssef rd_sgn = set_sign_bit<0, 1, 1, 1>(shuffle<0>(rd_ss)); - ssef mul_zxxy = shuffle<2, 0, 0, 1>(vdir) * rd_sgn; - ssef mul_yz = shuffle<1, 2, 1, 2>(vdir) * mul_zxxy; - ssef mul_shuf = shuffle<0, 1, 2, 3>(mul_zxxy, mul_yz); - ssef vdir0 = vdir & cast(ssei(0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0)); - - ssef htfm0 = shuffle<0, 2, 0, 3>(mul_shuf, vdir0); - ssef htfm1 = shuffle<1, 0, 1, 3>(load1f_first(extract<0>(d_ss)), vdir0); - ssef htfm2 = shuffle<1, 3, 2, 3>(mul_shuf, vdir0); - -# if defined(__KERNEL_AVX2__) && defined(__KERNEL_SSE__) && \ - (!defined(_MSC_VER) || _MSC_VER > 1800) - const avxf vPP = _mm256_broadcast_ps(&P.m128); - const avxf htfm00 = avxf(htfm0.m128, htfm0.m128); - const avxf htfm11 = avxf(htfm1.m128, htfm1.m128); - const avxf htfm22 = avxf(htfm2.m128, htfm2.m128); - - const avxf p01 = madd( - shuffle<0>(P_curve_0_1 - vPP), - htfm00, - madd(shuffle<1>(P_curve_0_1 - vPP), htfm11, shuffle<2>(P_curve_0_1 - vPP) * htfm22)); - const avxf p23 = madd( - shuffle<0>(P_curve_2_3 - vPP), - htfm00, - madd(shuffle<1>(P_curve_2_3 - vPP), htfm11, shuffle<2>(P_curve_2_3 - vPP) * htfm22)); - - const ssef p0 = _mm256_castps256_ps128(p01); - const ssef p1 = _mm256_extractf128_ps(p01, 1); - const ssef p2 = _mm256_castps256_ps128(p23); - const ssef p3 = _mm256_extractf128_ps(p23, 1); - - const ssef P_curve_1 = _mm256_extractf128_ps(P_curve_0_1, 1); - r_st = ((float4 &)P_curve_1).w; - const ssef P_curve_2 = _mm256_castps256_ps128(P_curve_2_3); - r_en = ((float4 &)P_curve_2).w; -# else /* __KERNEL_AVX2__ */ - ssef htfm[] = {htfm0, htfm1, htfm2}; - ssef vP = load4f(P); - ssef p0 = transform_point_T3(htfm, P_curve[0] - vP); - ssef p1 = transform_point_T3(htfm, P_curve[1] - vP); - ssef p2 = transform_point_T3(htfm, P_curve[2] - vP); - ssef p3 = transform_point_T3(htfm, P_curve[3] - vP); - - r_st = ((float4 &)P_curve[1]).w; - r_en = ((float4 &)P_curve[2]).w; -# endif /* __KERNEL_AVX2__ */ - - float fc = 0.71f; - ssef vfc = ssef(fc); - ssef vfcxp3 = vfc * p3; - - vcurve_coef[0] = p1; - vcurve_coef[1] = vfc * (p2 - p0); - vcurve_coef[2] = madd( - ssef(fc * 2.0f), p0, madd(ssef(fc - 3.0f), p1, msub(ssef(3.0f - 2.0f * fc), p2, vfcxp3))); - vcurve_coef[3] = msub(ssef(fc - 2.0f), p2 - p1, msub(vfc, p0, vfcxp3)); } -# else - float3 curve_coef[4]; - /* curve Intersection check */ - /* obtain curve parameters */ + /* Standard case for rays that are not parallel to the cylinder. */ + const float Q = sqrtf(D); + const float rcp_2A = 1.0f / (2.0f * A); + const float t0 = (-B - Q) * rcp_2A; + const float t1 = (-B + Q) * rcp_2A; + + /* Calculates u and Ng for near hit. */ { - /* ray transform created - this should be created at beginning of intersection loop */ - Transform htfm; - float d = sqrtf(dir.x * dir.x + dir.z * dir.z); - htfm = make_transform(dir.z / d, - 0, - -dir.x / d, - 0, - -dir.x * dir.y / d, - d, - -dir.y * dir.z / d, - 0, - dir.x, - dir.y, - dir.z, - 0); - - float4 v00 = kernel_tex_fetch(__curves, prim); - - int k0 = __float_as_int(v00.x) + segment; - int k1 = k0 + 1; - - int ka = max(k0 - 1, __float_as_int(v00.x)); - int kb = min(k1 + 1, __float_as_int(v00.x) + __float_as_int(v00.y) - 1); - - float4 P_curve[4]; - - if (is_curve_primitive) { - P_curve[0] = kernel_tex_fetch(__curve_keys, ka); - P_curve[1] = kernel_tex_fetch(__curve_keys, k0); - P_curve[2] = kernel_tex_fetch(__curve_keys, k1); - P_curve[3] = kernel_tex_fetch(__curve_keys, kb); - } - else { - int fobject = (object == OBJECT_NONE) ? kernel_tex_fetch(__prim_object, curveAddr) : object; - motion_cardinal_curve_keys(kg, fobject, prim, time, ka, k0, k1, kb, P_curve); - } + *u0_o = (t0 * dOz + Oz) * rl; + const float3 Pr = t0 * ray_dir; + const float3 Pl = (*u0_o) * (cylinder_end - cylinder_start) + cylinder_start; + *Ng0_o = Pr - Pl; + } - float3 p0 = transform_point(&htfm, float4_to_float3(P_curve[0]) - P); - float3 p1 = transform_point(&htfm, float4_to_float3(P_curve[1]) - P); - float3 p2 = transform_point(&htfm, float4_to_float3(P_curve[2]) - P); - float3 p3 = transform_point(&htfm, float4_to_float3(P_curve[3]) - P); - - float fc = 0.71f; - curve_coef[0] = p1; - curve_coef[1] = -fc * p0 + fc * p2; - curve_coef[2] = 2.0f * fc * p0 + (fc - 3.0f) * p1 + (3.0f - 2.0f * fc) * p2 - fc * p3; - curve_coef[3] = -fc * p0 + (2.0f - fc) * p1 + (fc - 2.0f) * p2 + fc * p3; - r_st = P_curve[1].w; - r_en = P_curve[2].w; + /* Calculates u and Ng for far hit. */ + { + *u1_o = (t1 * dOz + Oz) * rl; + const float3 Pr = t1 * ray_dir; + const float3 Pl = (*u1_o) * (cylinder_end - cylinder_start) + cylinder_start; + *Ng1_o = Pr - Pl; } -# endif - float r_curr = max(r_st, r_en); - - if ((flags & CURVE_KN_RIBBONS) || !(flags & CURVE_KN_BACKFACING)) - epsilon = 2 * r_curr; - - /* find bounds - this is slow for cubic curves */ - float upper, lower; - - float zextrem[4]; - 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); - if (lower - r_curr > isect->t || upper + r_curr < epsilon) - return false; + *t_o = make_float2(t0, t1); - /* minimum width extension */ - float xextrem[4]; - 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); - if (lower > r_curr || upper < -r_curr) - return false; + return true; +} - float yextrem[4]; - 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); - if (lower > r_curr || upper < -r_curr) - return false; +ccl_device_inline float2 half_plane_intersect(const float3 P, const float3 N, const float3 ray_dir) +{ + const float3 O = -P; + const float3 D = ray_dir; + const float ON = dot(O, N); + const float DN = dot(D, N); + const float min_rcp_input = 1e-18f; + const bool eps = fabsf(DN) < min_rcp_input; + const float t = -ON / DN; + const float lower = (eps || DN < 0.0f) ? -FLT_MAX : t; + const float upper = (eps || DN > 0.0f) ? FLT_MAX : t; + return make_float2(lower, upper); +} - /* setup recurrent loop */ - int level = 1 << depth; - int tree = 0; - float resol = 1.0f / (float)level; - bool hit = false; - - /* begin loop */ - while (!(tree >> (depth))) { - const float i_st = tree * resol; - const float i_en = i_st + (level * resol); - -# ifdef __KERNEL_SSE2__ - ssef vi_st = ssef(i_st), vi_en = ssef(i_en); - ssef vp_st = madd(madd(madd(vcurve_coef[3], vi_st, vcurve_coef[2]), vi_st, vcurve_coef[1]), - vi_st, - vcurve_coef[0]); - ssef vp_en = madd(madd(madd(vcurve_coef[3], vi_en, vcurve_coef[2]), vi_en, vcurve_coef[1]), - vi_en, - vcurve_coef[0]); - - ssef vbmin = min(vp_st, vp_en); - ssef vbmax = max(vp_st, vp_en); - - float3 &bmin = (float3 &)vbmin, &bmax = (float3 &)vbmax; - float &bminx = bmin.x, &bminy = bmin.y, &bminz = bmin.z; - float &bmaxx = bmax.x, &bmaxy = bmax.y, &bmaxz = bmax.z; - float3 &p_st = (float3 &)vp_st, &p_en = (float3 &)vp_en; -# else - float3 p_st = ((curve_coef[3] * i_st + curve_coef[2]) * i_st + curve_coef[1]) * i_st + - curve_coef[0]; - float3 p_en = ((curve_coef[3] * i_en + curve_coef[2]) * i_en + curve_coef[1]) * i_en + - curve_coef[0]; - - float bminx = min(p_st.x, p_en.x); - float bmaxx = max(p_st.x, p_en.x); - float bminy = min(p_st.y, p_en.y); - float bmaxy = max(p_st.y, p_en.y); - float bminz = min(p_st.z, p_en.z); - float bmaxz = max(p_st.z, p_en.z); -# endif +ccl_device bool curve_intersect_iterative(const float3 ray_dir, + const float dt, + const float4 curve[4], + float u, + float t, + const bool use_backfacing, + Intersection *isect) +{ + const float length_ray_dir = len(ray_dir); + + /* Error of curve evaluations is proportional to largest coordinate. */ + const float4 box_min = min(min(curve[0], curve[1]), min(curve[2], curve[3])); + const float4 box_max = max(min(curve[0], curve[1]), max(curve[2], curve[3])); + const float4 box_abs = max(fabs(box_min), fabs(box_max)); + const float P_err = 16.0f * FLT_EPSILON * + max(box_abs.x, max(box_abs.y, max(box_abs.z, box_abs.w))); + const float radius_max = box_max.w; + + for (int i = 0; i < CURVE_NUM_JACOBIAN_ITERATIONS; i++) { + const float3 Q = ray_dir * t; + const float3 dQdt = ray_dir; + const float Q_err = 16.0f * FLT_EPSILON * length_ray_dir * t; + + const float4 P4 = catmull_rom_basis_eval(curve, u); + const float4 dPdu4 = catmull_rom_basis_derivative(curve, u); + + const float3 P = float4_to_float3(P4); + const float3 dPdu = float4_to_float3(dPdu4); + const float radius = P4.w; + const float dradiusdu = dPdu4.w; + + const float3 ddPdu = float4_to_float3(catmull_rom_basis_derivative2(curve, u)); + + const float3 R = Q - P; + const float len_R = len(R); + const float R_err = max(Q_err, P_err); + const float3 dRdu = -dPdu; + const float3 dRdt = dQdt; + + const float3 T = normalize(dPdu); + const float3 dTdu = dnormalize(dPdu, ddPdu); + const float cos_err = P_err / len(dPdu); + + const float f = dot(R, T); + const float f_err = len_R * P_err + R_err + cos_err * (1.0f + len_R); + const float dfdu = dot(dRdu, T) + dot(R, dTdu); + const float dfdt = dot(dRdt, T); + + const float K = dot(R, R) - sqr(f); + const float dKdu = (dot(R, dRdu) - f * dfdu); + const float dKdt = (dot(R, dRdt) - f * dfdt); + const float rsqrt_K = inversesqrtf(K); + + const float g = sqrtf(K) - radius; + const float g_err = R_err + f_err + 16.0f * FLT_EPSILON * radius_max; + const float dgdu = dKdu * rsqrt_K - dradiusdu; + const float dgdt = dKdt * rsqrt_K; + + const float invdet = 1.0f / (dfdu * dgdt - dgdu * dfdt); + u -= (dgdt * f - dfdt * g) * invdet; + t -= (-dgdu * f + dfdu * g) * invdet; + + if (fabsf(f) < f_err && fabsf(g) < g_err) { + t += dt; + if (!(0.0f <= t && t <= isect->t)) { + return false; /* Rejects NaNs */ + } + if (!(u >= 0.0f && u <= 1.0f)) { + return false; /* Rejects NaNs */ + } - if (xextrem[0] >= i_st && xextrem[0] <= i_en) { - bminx = min(bminx, xextrem[1]); - bmaxx = max(bmaxx, xextrem[1]); - } - if (xextrem[2] >= i_st && xextrem[2] <= i_en) { - bminx = min(bminx, xextrem[3]); - bmaxx = max(bmaxx, xextrem[3]); - } - if (yextrem[0] >= i_st && yextrem[0] <= i_en) { - bminy = min(bminy, yextrem[1]); - bmaxy = max(bmaxy, yextrem[1]); - } - if (yextrem[2] >= i_st && yextrem[2] <= i_en) { - bminy = min(bminy, yextrem[3]); - bmaxy = max(bmaxy, yextrem[3]); - } - if (zextrem[0] >= i_st && zextrem[0] <= i_en) { - bminz = min(bminz, zextrem[1]); - bmaxz = max(bmaxz, zextrem[1]); - } - if (zextrem[2] >= i_st && zextrem[2] <= i_en) { - bminz = min(bminz, zextrem[3]); - bmaxz = max(bmaxz, zextrem[3]); - } + /* Backface culling. */ + const float3 R = normalize(Q - P); + const float3 U = dradiusdu * R + dPdu; + const float3 V = cross(dPdu, R); + const float3 Ng = cross(V, U); + if (!use_backfacing && dot(ray_dir, Ng) > 0.0f) { + return false; + } - float r1 = r_st + (r_en - r_st) * i_st; - float r2 = r_st + (r_en - r_st) * i_en; - r_curr = max(r1, r2); + /* Record intersection. */ + isect->t = t; + isect->u = u; + isect->v = 0.0f; - if (bminz - r_curr > isect->t || bmaxz + r_curr < epsilon || bminx > r_curr || - bmaxx < -r_curr || bminy > r_curr || bmaxy < -r_curr) { - /* the bounding box does not overlap the square centered at O */ - tree += level; - level = tree & -tree; + return true; } - else if (level == 1) { - - /* the maximum recursion depth is reached. - * check if dP0.(Q-P0)>=0 and dPn.(Pn-Q)>=0. - * dP* is reversed if necessary.*/ - float t = isect->t; - float u = 0.0f; - float gd = 0.0f; - - if (flags & CURVE_KN_RIBBONS) { - float3 tg = (p_en - p_st); -# ifdef __KERNEL_SSE__ - const float3 tg_sq = tg * tg; - float w = tg_sq.x + tg_sq.y; -# else - float w = tg.x * tg.x + tg.y * tg.y; -# endif - if (w == 0) { - tree++; - level = tree & -tree; - continue; - } -# ifdef __KERNEL_SSE__ - const float3 p_sttg = p_st * tg; - w = -(p_sttg.x + p_sttg.y) / w; + } + return false; +} + +ccl_device bool curve_intersect_recursive(const float3 ray_orig, + const float3 ray_dir, + float4 curve[4], + Intersection *isect) +{ + /* Move ray closer to make intersection stable. */ + const float3 center = float4_to_float3(0.25f * (curve[0] + curve[1] + curve[2] + curve[3])); + const float dt = dot(center - ray_orig, ray_dir) / dot(ray_dir, ray_dir); + const float3 ref = ray_orig + ray_dir * dt; + const float4 ref4 = make_float4(ref.x, ref.y, ref.z, 0.0f); + curve[0] -= ref4; + curve[1] -= ref4; + curve[2] -= ref4; + curve[3] -= ref4; + + const bool use_backfacing = false; + const float step_size = 1.0f / (float)(CURVE_NUM_BEZIER_STEPS); + + int depth = 0; + + /* todo: optimize stack for GPU somehow? Possibly some bitflags are enough, and + * u0/u1 can be derived from the depth. */ + struct { + float u0, u1; + int i; + } stack[CURVE_NUM_BEZIER_SUBDIVISIONS_UNSTABLE]; + + bool found = false; + + float u0 = 0.0f; + float u1 = 1.0f; + int i = 0; + + while (1) { + for (; i < CURVE_NUM_BEZIER_STEPS; i++) { + const float step = i * step_size; + + /* Subdivide curve. */ + const float dscale = (u1 - u0) * (1.0f / 3.0f) * step_size; + const float vu0 = mix(u0, u1, step); + const float vu1 = mix(u0, u1, step + step_size); + + const float4 P0 = catmull_rom_basis_eval(curve, vu0); + const float4 dP0du = dscale * catmull_rom_basis_derivative(curve, vu0); + const float4 P3 = catmull_rom_basis_eval(curve, vu1); + const float4 dP3du = dscale * catmull_rom_basis_derivative(curve, vu1); + + const float4 P1 = P0 + dP0du; + const float4 P2 = P3 - dP3du; + + /* Calculate bounding cylinders. */ + const float rr1 = sqr_point_to_line_distance(float4_to_float3(dP0du), + float4_to_float3(P3 - P0)); + const float rr2 = sqr_point_to_line_distance(float4_to_float3(dP3du), + float4_to_float3(P3 - P0)); + const float maxr12 = sqrtf(max(rr1, rr2)); + const float one_plus_ulp = 1.0f + 2.0f * FLT_EPSILON; + const float one_minus_ulp = 1.0f - 2.0f * FLT_EPSILON; + float r_outer = max(max(P0.w, P1.w), max(P2.w, P3.w)) + maxr12; + float r_inner = min(min(P0.w, P1.w), min(P2.w, P3.w)) - maxr12; + r_outer = one_plus_ulp * r_outer; + r_inner = max(0.0f, one_minus_ulp * r_inner); + bool valid = true; + + /* Intersect with outer cylinder. */ + float2 tc_outer; + float u_outer0, u_outer1; + float3 Ng_outer0, Ng_outer1; + valid = cylinder_intersect(float4_to_float3(P0), + float4_to_float3(P3), + r_outer, + ray_dir, + &tc_outer, + &u_outer0, + &Ng_outer0, + &u_outer1, + &Ng_outer1); + if (!valid) { + continue; + } + + /* Intersect with cap-planes. */ + float2 tp = make_float2(-dt, isect->t - dt); + tp = make_float2(max(tp.x, tc_outer.x), min(tp.y, tc_outer.y)); + const float2 h0 = half_plane_intersect( + float4_to_float3(P0), float4_to_float3(dP0du), ray_dir); + tp = make_float2(max(tp.x, h0.x), min(tp.y, h0.y)); + const float2 h1 = half_plane_intersect( + float4_to_float3(P3), -float4_to_float3(dP3du), ray_dir); + tp = make_float2(max(tp.x, h1.x), min(tp.y, h1.y)); + valid = tp.x <= tp.y; + if (!valid) { + continue; + } + + /* Clamp and correct u parameter. */ + u_outer0 = clamp(u_outer0, 0.0f, 1.0f); + u_outer1 = clamp(u_outer1, 0.0f, 1.0f); + u_outer0 = mix(u0, u1, (step + u_outer0) * (1.0f / (float)(CURVE_NUM_BEZIER_STEPS + 1))); + u_outer1 = mix(u0, u1, (step + u_outer1) * (1.0f / (float)(CURVE_NUM_BEZIER_STEPS + 1))); + + /* Intersect with inner cylinder. */ + float2 tc_inner; + float u_inner0, u_inner1; + float3 Ng_inner0, Ng_inner1; + const bool valid_inner = cylinder_intersect(float4_to_float3(P0), + float4_to_float3(P3), + r_inner, + ray_dir, + &tc_inner, + &u_inner0, + &Ng_inner0, + &u_inner1, + &Ng_inner1); + + /* At the unstable area we subdivide deeper. */ +# if 0 + const bool unstable0 = (!valid_inner) | + (fabsf(dot(normalize(ray_dir), normalize(Ng_inner0))) < 0.3f); + const bool unstable1 = (!valid_inner) | + (fabsf(dot(normalize(ray_dir), normalize(Ng_inner1))) < 0.3f); # else - w = -(p_st.x * tg.x + p_st.y * tg.y) / w; -# endif - w = saturate(w); - - /* compute u on the curve segment */ - u = i_st * (1 - w) + i_en * w; - r_curr = r_st + (r_en - r_st) * u; - /* compare x-y distances */ - float3 p_curr = ((curve_coef[3] * u + curve_coef[2]) * u + curve_coef[1]) * u + - curve_coef[0]; - - float3 dp_st = (3 * curve_coef[3] * i_st + 2 * curve_coef[2]) * i_st + curve_coef[1]; - if (dot(tg, dp_st) < 0) - dp_st *= -1; - if (dot(dp_st, -p_st) + p_curr.z * dp_st.z < 0) { - tree++; - level = tree & -tree; - continue; - } - float3 dp_en = (3 * curve_coef[3] * i_en + 2 * curve_coef[2]) * i_en + curve_coef[1]; - if (dot(tg, dp_en) < 0) - dp_en *= -1; - if (dot(dp_en, p_en) - p_curr.z * dp_en.z < 0) { - tree++; - level = tree & -tree; - continue; - } + /* On the GPU appears to be a little faster if always enabled. */ + (void)valid_inner; - if (p_curr.x * p_curr.x + p_curr.y * p_curr.y >= r_curr * r_curr || p_curr.z <= epsilon || - isect->t < p_curr.z) { - tree++; - level = tree & -tree; - continue; - } + const bool unstable0 = true; + const bool unstable1 = true; +# endif - t = p_curr.z; + /* Subtract the inner interval from the current hit interval. */ + float2 tp0 = make_float2(tp.x, min(tp.y, tc_inner.x)); + float2 tp1 = make_float2(max(tp.x, tc_inner.y), tp.y); + bool valid0 = valid && (tp0.x <= tp0.y); + bool valid1 = valid && (tp1.x <= tp1.y); + if (!(valid0 || valid1)) { + continue; } - else { - float l = len(p_en - p_st); - float invl = 1.0f / l; - float3 tg = (p_en - p_st) * invl; - gd = (r2 - r1) * invl; - float difz = -dot(p_st, tg); - float cyla = 1.0f - (tg.z * tg.z * (1 + gd * gd)); - float invcyla = 1.0f / cyla; - float halfb = (-p_st.z - tg.z * (difz + gd * (difz * gd + r1))); - float tcentre = -halfb * invcyla; - float zcentre = difz + (tg.z * tcentre); - float3 tdif = -p_st; - tdif.z += tcentre; - float tdifz = dot(tdif, tg); - float tb = 2 * (tdif.z - tg.z * (tdifz + gd * (tdifz * gd + r1))); - float tc = dot(tdif, tdif) - tdifz * tdifz * (1 + gd * gd) - r1 * r1 - 2 * r1 * tdifz * gd; - float td = tb * tb - 4 * cyla * tc; - if (td < 0.0f) { - tree++; - level = tree & -tree; - continue; - } - float rootd = sqrtf(td); - float correction = (-tb - rootd) * 0.5f * invcyla; - t = tcentre + correction; - - float3 dp_st = (3 * curve_coef[3] * i_st + 2 * curve_coef[2]) * i_st + curve_coef[1]; - if (dot(tg, dp_st) < 0) - dp_st *= -1; - float3 dp_en = (3 * curve_coef[3] * i_en + 2 * curve_coef[2]) * i_en + curve_coef[1]; - if (dot(tg, dp_en) < 0) - dp_en *= -1; - - 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)) { - correction = (-tb + rootd) * 0.5f * invcyla; - t = tcentre + correction; + /* Process one or two hits. */ + bool recurse = false; + if (valid0) { + const int termDepth = unstable0 ? CURVE_NUM_BEZIER_SUBDIVISIONS_UNSTABLE : + CURVE_NUM_BEZIER_SUBDIVISIONS; + if (depth >= termDepth) { + found |= curve_intersect_iterative( + ray_dir, dt, curve, u_outer0, tp0.x, use_backfacing, isect); } - - 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) { - tree++; - level = tree & -tree; - continue; + else { + recurse = true; } + } - float w = (zcentre + (tg.z * correction)) * invl; - w = saturate(w); - /* compute u on the curve segment */ - u = i_st * (1 - w) + i_en * w; + if (valid1 && (tp1.x + dt <= isect->t)) { + const int termDepth = unstable1 ? CURVE_NUM_BEZIER_SUBDIVISIONS_UNSTABLE : + CURVE_NUM_BEZIER_SUBDIVISIONS; + if (depth >= termDepth) { + found |= curve_intersect_iterative( + ray_dir, dt, curve, u_outer1, tp1.y, use_backfacing, isect); + } + else { + recurse = true; + } } - /* we found a new intersection */ -# ifdef __VISIBILITY_FLAG__ - /* visibility flag test. we do it here under the assumption - * that most triangles are culled by node flags */ - if (kernel_tex_fetch(__prim_visibility, curveAddr) & visibility) -# endif - { - /* record intersection */ - isect->t = t; - isect->u = u; - isect->v = gd; - isect->prim = curveAddr; - isect->object = object; - isect->type = type; - hit = true; + if (recurse) { + stack[depth].u0 = u0; + stack[depth].u1 = u1; + stack[depth].i = i + 1; + depth++; + + u0 = vu0; + u1 = vu1; + i = -1; } + } - tree++; - level = tree & -tree; + if (depth > 0) { + depth--; + u0 = stack[depth].u0; + u1 = stack[depth].u1; + i = stack[depth].i; } else { - /* split the curve into two curves and process */ - level = level >> 1; + break; } } - return hit; + return found; +} + +/* Ribbons */ + +ccl_device_inline bool cylinder_culling_test(const float2 p1, const float2 p2, const float r) +{ + /* Performs culling against a cylinder. */ + const float2 dp = p2 - p1; + const float num = dp.x * p1.y - dp.y * p1.x; + const float den2 = dot(p2 - p1, p2 - p1); + return num * num <= r * r * den2; +} + +/*! Intersects a ray with a quad with backface culling + * enabled. The quad v0,v1,v2,v3 is split into two triangles + * v0,v1,v3 and v2,v3,v1. The edge v1,v2 decides which of the two + * triangles gets intersected. */ +ccl_device_inline bool ribbon_intersect_quad(const float ray_tfar, + const float3 quad_v0, + const float3 quad_v1, + const float3 quad_v2, + const float3 quad_v3, + float *u_o, + float *v_o, + float *t_o) +{ + /* Calculate vertices relative to ray origin? */ + const float3 O = make_float3(0.0f, 0.0f, 0.0f); + const float3 D = make_float3(0.0f, 0.0f, 1.0f); + const float3 va = quad_v0 - O; + const float3 vb = quad_v1 - O; + const float3 vc = quad_v2 - O; + const float3 vd = quad_v3 - O; + + const float3 edb = vb - vd; + const float WW = dot(cross(vd, edb), D); + const float3 v0 = (WW <= 0.0f) ? va : vc; + const float3 v1 = (WW <= 0.0f) ? vb : vd; + const float3 v2 = (WW <= 0.0f) ? vd : vb; + + /* Calculate edges? */ + const float3 e0 = v2 - v0; + const float3 e1 = v0 - v1; + + /* perform edge tests */ + const float U = dot(cross(v0, e0), D); + const float V = dot(cross(v1, e1), D); + if (!(max(U, V) <= 0.0f)) { + return false; + } + + /* Calculate geometry normal and denominator? */ + const float3 Ng = cross(e1, e0); + const float den = dot(Ng, D); + const float rcpDen = 1.0f / den; + + /* Perform depth test? */ + const float t = rcpDen * dot(v0, Ng); + if (!(0.0f <= t && t <= ray_tfar)) { + return false; + } + + /* Avoid division by 0? */ + if (!(den != 0.0f)) { + return false; + } + + /* Update hit information? */ + *t_o = t; + *u_o = U * rcpDen; + *v_o = V * rcpDen; + *u_o = (WW <= 0.0f) ? *u_o : 1.0f - *u_o; + *v_o = (WW <= 0.0f) ? *v_o : 1.0f - *v_o; + return true; +} + +ccl_device_inline void ribbon_ray_space(const float3 ray_dir, float3 ray_space[3]) +{ + const float3 dx0 = make_float3(0, ray_dir.z, -ray_dir.y); + const float3 dx1 = make_float3(-ray_dir.z, 0, ray_dir.x); + ray_space[0] = normalize(dot(dx0, dx0) > dot(dx1, dx1) ? dx0 : dx1); + ray_space[1] = normalize(cross(ray_dir, ray_space[0])); + ray_space[2] = ray_dir; +} + +ccl_device_inline float4 ribbon_to_ray_space(const float3 ray_space[3], + const float3 ray_org, + const float4 P4) +{ + float3 P = float4_to_float3(P4) - ray_org; + return make_float4(dot(ray_space[0], P), dot(ray_space[1], P), dot(ray_space[2], P), P4.w); +} + +ccl_device_inline bool ribbon_intersect(const float3 ray_org, + const float3 ray_dir, + const float ray_tfar, + const int N, + float4 curve[4], + Intersection *isect) +{ + /* Transform control points into ray space. */ + float3 ray_space[3]; + ribbon_ray_space(ray_dir, ray_space); + + curve[0] = ribbon_to_ray_space(ray_space, ray_org, curve[0]); + curve[1] = ribbon_to_ray_space(ray_space, ray_org, curve[1]); + curve[2] = ribbon_to_ray_space(ray_space, ray_org, curve[2]); + curve[3] = ribbon_to_ray_space(ray_space, ray_org, curve[3]); + + const float4 mx = max(max(fabs(curve[0]), fabs(curve[1])), max(fabs(curve[2]), fabs(curve[3]))); + const float eps = 4.0f * FLT_EPSILON * max(max(mx.x, mx.y), max(mx.z, mx.w)); + const float step_size = 1.0f / (float)N; + + /* Evaluate first point and radius scaled normal direction. */ + float4 p0 = catmull_rom_basis_eval(curve, 0.0f); + float3 dp0dt = float4_to_float3(catmull_rom_basis_derivative(curve, 0.0f)); + if (max3(fabs(dp0dt)) < eps) { + const float4 p1 = catmull_rom_basis_eval(curve, step_size); + dp0dt = float4_to_float3(p1 - p0); + } + float3 wn0 = normalize(make_float3(dp0dt.y, -dp0dt.x, 0.0f)) * p0.w; + + /* Evaluate the bezier curve. */ + for (int i = 0; i < N; i++) { + const float u = i * step_size; + const float4 p1 = catmull_rom_basis_eval(curve, u + step_size); + bool valid = cylinder_culling_test( + make_float2(p0.x, p0.y), make_float2(p1.x, p1.y), max(p0.w, p1.w)); + if (!valid) { + continue; + } + + /* Evaluate next point. */ + float3 dp1dt = float4_to_float3(catmull_rom_basis_derivative(curve, u + step_size)); + dp1dt = (max3(fabs(dp1dt)) < eps) ? float4_to_float3(p1 - p0) : dp1dt; + const float3 wn1 = normalize(make_float3(dp1dt.y, -dp1dt.x, 0.0f)) * p1.w; + + /* Construct quad coordinates. */ + const float3 lp0 = float4_to_float3(p0) + wn0; + const float3 lp1 = float4_to_float3(p1) + wn1; + const float3 up0 = float4_to_float3(p0) - wn0; + const float3 up1 = float4_to_float3(p1) - wn1; + + /* Intersect quad. */ + float vu, vv, vt; + bool valid0 = ribbon_intersect_quad(isect->t, lp0, lp1, up1, up0, &vu, &vv, &vt); + + if (valid0) { + /* ignore self intersections */ + const float avoidance_factor = 2.0f; + if (avoidance_factor != 0.0f) { + float r = mix(p0.w, p1.w, vu); + valid0 = vt > avoidance_factor * r; + } + + if (valid0) { + vv = 2.0f * vv - 1.0f; + + /* Record intersection. */ + isect->t = vt; + isect->u = u + vu * step_size; + isect->v = vv; + return true; + } + } + + p0 = p1; + wn0 = wn1; + } + return false; } ccl_device_forceinline bool curve_intersect(KernelGlobals *kg, Intersection *isect, - float3 P, - float3 direction, + const float3 P, + const float3 dir, uint visibility, int object, int curveAddr, float time, int type) { - /* define few macros to minimize code duplication for SSE */ -# ifndef __KERNEL_SSE2__ -# define len3_squared(x) len_squared(x) -# define len3(x) len(x) -# define dot3(x, y) dot(x, y) -# endif - - const bool is_curve_primitive = (type & PRIMITIVE_CURVE); + const bool is_motion = (type & PRIMITIVE_ALL_MOTION); -# ifndef __KERNEL_OPTIX__ /* see OptiX motion flag OPTIX_MOTION_FLAG_[START|END]_VANISH */ - if (!is_curve_primitive && kernel_data.bvh.use_bvh_steps) { +# ifndef __KERNEL_OPTIX__ /* See OptiX motion flag OPTIX_MOTION_FLAG_[START|END]_VANISH */ + if (is_motion && kernel_data.bvh.use_bvh_steps) { const float2 prim_time = kernel_tex_fetch(__prim_time, curveAddr); if (time < prim_time.x || time > prim_time.y) { return false; @@ -517,210 +639,63 @@ ccl_device_forceinline bool curve_intersect(KernelGlobals *kg, # endif int segment = PRIMITIVE_UNPACK_SEGMENT(type); - /* curve Intersection check */ - int flags = kernel_data.curve.curveflags; - int prim = kernel_tex_fetch(__prim_index, curveAddr); + float4 v00 = kernel_tex_fetch(__curves, prim); - int cnum = __float_as_int(v00.x); - int k0 = cnum + segment; + int k0 = __float_as_int(v00.x) + segment; int k1 = k0 + 1; -# ifndef __KERNEL_SSE2__ - float4 P_curve[2]; + int ka = max(k0 - 1, __float_as_int(v00.x)); + int kb = min(k1 + 1, __float_as_int(v00.x) + __float_as_int(v00.y) - 1); - if (is_curve_primitive) { - P_curve[0] = kernel_tex_fetch(__curve_keys, k0); - P_curve[1] = kernel_tex_fetch(__curve_keys, k1); + float4 curve[4]; + if (!is_motion) { + curve[0] = kernel_tex_fetch(__curve_keys, ka); + curve[1] = kernel_tex_fetch(__curve_keys, k0); + curve[2] = kernel_tex_fetch(__curve_keys, k1); + curve[3] = kernel_tex_fetch(__curve_keys, kb); } else { int fobject = (object == OBJECT_NONE) ? kernel_tex_fetch(__prim_object, curveAddr) : object; - motion_curve_keys(kg, fobject, prim, time, k0, k1, P_curve); - } - - float r1 = P_curve[0].w; - float r2 = P_curve[1].w; - float3 p1 = float4_to_float3(P_curve[0]); - float3 p2 = float4_to_float3(P_curve[1]); - - /* minimum width extension */ - float3 dif = P - p1; - float3 dif_second = P - p2; - - float3 p21_diff = p2 - p1; - float3 sphere_dif1 = (dif + dif_second) * 0.5f; - float3 dir = direction; - float sphere_b_tmp = dot3(dir, sphere_dif1); - float3 sphere_dif2 = sphere_dif1 - sphere_b_tmp * dir; -# else - ssef P_curve[2]; - - if (is_curve_primitive) { - P_curve[0] = load4f(&kg->__curve_keys.data[k0].x); - P_curve[1] = load4f(&kg->__curve_keys.data[k1].x); + motion_curve_keys(kg, fobject, prim, time, ka, k0, k1, kb, curve); } - else { - int fobject = (object == OBJECT_NONE) ? kernel_tex_fetch(__prim_object, curveAddr) : object; - motion_curve_keys(kg, fobject, prim, time, k0, k1, (float4 *)&P_curve); - } - - ssef r12 = shuffle<3, 3, 3, 3>(P_curve[0], P_curve[1]); - const ssef vP = load4f(P); - const ssef dif = vP - P_curve[0]; - const ssef dif_second = vP - P_curve[1]; - float r1 = extract<0>(r12), r2 = extract<0>(shuffle<2>(r12)); - - const ssef p21_diff = P_curve[1] - P_curve[0]; - const ssef sphere_dif1 = (dif + dif_second) * 0.5f; - const ssef dir = load4f(direction); - const ssef sphere_b_tmp = dot3_splat(dir, sphere_dif1); - const ssef sphere_dif2 = nmadd(sphere_b_tmp, dir, sphere_dif1); -# endif - - float mr = max(r1, r2); - float l = len3(p21_diff); - float invl = 1.0f / l; - float sp_r = mr + 0.5f * l; - float sphere_b = dot3(dir, sphere_dif2); - float sdisc = sphere_b * sphere_b - len3_squared(sphere_dif2) + sp_r * sp_r; - - if (sdisc < 0.0f) - return false; - - /* obtain parameters and test midpoint distance for suitable modes */ -# ifndef __KERNEL_SSE2__ - float3 tg = p21_diff * invl; -# else - const ssef tg = p21_diff * invl; -# endif - float gd = (r2 - r1) * invl; - - float dirz = dot3(dir, tg); - float difz = dot3(dif, tg); - - float a = 1.0f - (dirz * dirz * (1 + gd * gd)); - - float halfb = dot3(dir, dif) - dirz * (difz + gd * (difz * gd + r1)); - - float tcentre = -halfb / a; - float zcentre = difz + (dirz * tcentre); - - if ((tcentre > isect->t) && !(flags & CURVE_KN_ACCURATE)) - return false; - if ((zcentre < 0 || zcentre > l) && !(flags & CURVE_KN_ACCURATE) && - !(flags & CURVE_KN_INTERSECTCORRECTION)) +# ifdef __VISIBILITY_FLAG__ + if (!(kernel_tex_fetch(__prim_visibility, curveAddr) & visibility)) { return false; - - /* test minimum separation */ -# ifndef __KERNEL_SSE2__ - float3 cprod = cross(tg, dir); - float cprod2sq = len3_squared(cross(tg, dif)); -# else - const ssef cprod = cross(tg, dir); - float cprod2sq = len3_squared(cross_zxy(tg, dif)); + } # endif - float cprodsq = len3_squared(cprod); - float distscaled = dot3(cprod, dif); - - if (cprodsq == 0) - distscaled = cprod2sq; - else - distscaled = (distscaled * distscaled) / cprodsq; - - if (distscaled > mr * mr) - return false; - /* calculate true intersection */ -# ifndef __KERNEL_SSE2__ - float3 tdif = dif + tcentre * dir; -# else - const ssef tdif = madd(ssef(tcentre), dir, dif); -# endif - float tdifz = dot3(tdif, tg); - float tdifma = tdifz * gd + r1; - float tb = 2 * (dot3(dir, tdif) - dirz * (tdifz + gd * tdifma)); - float tc = dot3(tdif, tdif) - tdifz * tdifz - tdifma * tdifma; - float td = tb * tb - 4 * a * tc; + if (type & (PRIMITIVE_CURVE_RIBBON | PRIMITIVE_MOTION_CURVE_RIBBON)) { + /* todo: adaptive number of subdivisions could help performance here. */ + const int subdivisions = kernel_data.bvh.curve_subdivisions; + if (ribbon_intersect(P, dir, isect->t, subdivisions, curve, isect)) { + isect->prim = curveAddr; + isect->object = object; + isect->type = type; + return true; + } - if (td < 0.0f) return false; - - float rootd = 0.0f; - float correction = 0.0f; - if (flags & CURVE_KN_ACCURATE) { - rootd = sqrtf(td); - correction = ((-tb - rootd) / (2 * a)); } - - float t = tcentre + correction; - - if (t < isect->t) { - - if (flags & CURVE_KN_INTERSECTCORRECTION) { - rootd = sqrtf(td); - correction = ((-tb - rootd) / (2 * a)); - t = tcentre + correction; - } - - float z = zcentre + (dirz * correction); - // bool backface = false; - - if (flags & CURVE_KN_BACKFACING && (t < 0.0f || z < 0 || z > l)) { - // backface = true; - correction = ((-tb + rootd) / (2 * a)); - t = tcentre + correction; - z = zcentre + (dirz * correction); + else { + if (curve_intersect_recursive(P, dir, curve, isect)) { + isect->prim = curveAddr; + isect->object = object; + isect->type = type; + return true; } - if (t > 0.0f && t < isect->t && z >= 0 && z <= l) { - - if (flags & CURVE_KN_ENCLOSEFILTER) { - float enc_ratio = 1.01f; - if ((difz > -r1 * enc_ratio) && (dot3(dif_second, tg) < r2 * enc_ratio)) { - float a2 = 1.0f - (dirz * dirz * (1 + gd * gd * enc_ratio * enc_ratio)); - 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; - if (a2 * c2 < 0.0f) - return false; - } - } - -# ifdef __VISIBILITY_FLAG__ - /* visibility flag test. we do it here under the assumption - * that most triangles are culled by node flags */ - if (kernel_tex_fetch(__prim_visibility, curveAddr) & visibility) -# endif - { - /* record intersection */ - isect->t = t; - isect->u = z * invl; - isect->v = gd; - isect->prim = curveAddr; - isect->object = object; - isect->type = type; - - return true; - } - } + return false; } - - return false; - -# ifndef __KERNEL_SSE2__ -# undef len3_squared -# undef len3 -# undef dot3 -# endif } -ccl_device_inline float3 curve_refine(KernelGlobals *kg, - ShaderData *sd, - const Intersection *isect, - const Ray *ray) +ccl_device_inline void curve_shader_setup(KernelGlobals *kg, + ShaderData *sd, + const Intersection *isect, + const Ray *ray) { - int flag = kernel_data.curve.curveflags; float t = isect->t; float3 P = ray->P; float3 D = ray->D; @@ -743,118 +718,60 @@ ccl_device_inline float3 curve_refine(KernelGlobals *kg, int k0 = __float_as_int(v00.x) + PRIMITIVE_UNPACK_SEGMENT(sd->type); int k1 = k0 + 1; - float3 tg; + int ka = max(k0 - 1, __float_as_int(v00.x)); + int kb = min(k1 + 1, __float_as_int(v00.x) + __float_as_int(v00.y) - 1); - if (flag & CURVE_KN_INTERPOLATE) { - int ka = max(k0 - 1, __float_as_int(v00.x)); - int kb = min(k1 + 1, __float_as_int(v00.x) + __float_as_int(v00.y) - 1); + float4 P_curve[4]; - float4 P_curve[4]; + if (!(sd->type & PRIMITIVE_ALL_MOTION)) { + P_curve[0] = kernel_tex_fetch(__curve_keys, ka); + P_curve[1] = kernel_tex_fetch(__curve_keys, k0); + P_curve[2] = kernel_tex_fetch(__curve_keys, k1); + P_curve[3] = kernel_tex_fetch(__curve_keys, kb); + } + else { + motion_curve_keys(kg, sd->object, sd->prim, sd->time, ka, k0, k1, kb, P_curve); + } - if (sd->type & PRIMITIVE_CURVE) { - P_curve[0] = kernel_tex_fetch(__curve_keys, ka); - P_curve[1] = kernel_tex_fetch(__curve_keys, k0); - P_curve[2] = kernel_tex_fetch(__curve_keys, k1); - P_curve[3] = kernel_tex_fetch(__curve_keys, kb); - } - else { - motion_cardinal_curve_keys(kg, sd->object, sd->prim, sd->time, ka, k0, k1, kb, P_curve); - } + sd->u = isect->u; + sd->v = isect->v; - float3 p[4]; - p[0] = float4_to_float3(P_curve[0]); - p[1] = float4_to_float3(P_curve[1]); - p[2] = float4_to_float3(P_curve[2]); - p[3] = float4_to_float3(P_curve[3]); + P = P + D * t; - P = P + D * t; + const float4 dPdu4 = catmull_rom_basis_derivative(P_curve, isect->u); + const float3 dPdu = float4_to_float3(dPdu4); -# ifdef __UV__ - sd->u = isect->u; - sd->v = 0.0f; -# endif + if (sd->type & (PRIMITIVE_CURVE_RIBBON | PRIMITIVE_MOTION_CURVE_RIBBON)) { + /* Rounded smooth normals for ribbons, to approximate thick curve shape. */ + const float3 tangent = normalize(dPdu); + const float3 bitangent = normalize(cross(tangent, -D)); + const float sine = isect->v; + const float cosine = safe_sqrtf(1.0f - sine * sine); - tg = normalize(curvetangent(isect->u, p[0], p[1], p[2], p[3])); + sd->N = normalize(sine * bitangent - cosine * normalize(cross(tangent, bitangent))); + sd->Ng = -D; - if (kernel_data.curve.curveflags & CURVE_KN_RIBBONS) { - sd->Ng = normalize(-(D - tg * (dot(tg, D)))); - } - else { -# ifdef __EMBREE__ - if (kernel_data.bvh.scene) { - sd->Ng = normalize(isect->Ng); - } - else +# if 0 + /* This approximates the position and geometric normal of a thick curve too, + * but gives too many issues with wrong self intersections. */ + const float dPdu_radius = dPdu4.w; + sd->Ng = sd->N; + P += sd->N * dPdu_radius; # endif - { - /* direction from inside to surface of curve */ - float3 p_curr = curvepoint(isect->u, p[0], p[1], p[2], p[3]); - sd->Ng = normalize(P - p_curr); - - /* adjustment for changing radius */ - float gd = isect->v; - - if (gd != 0.0f) { - sd->Ng = sd->Ng - gd * tg; - sd->Ng = normalize(sd->Ng); - } - } - } - - /* todo: sometimes the normal is still so that this is detected as - * backfacing even if cull backfaces is enabled */ - - sd->N = sd->Ng; } else { - float4 P_curve[2]; - - if (sd->type & PRIMITIVE_CURVE) { - P_curve[0] = kernel_tex_fetch(__curve_keys, k0); - P_curve[1] = kernel_tex_fetch(__curve_keys, k1); - } - else { - motion_curve_keys(kg, sd->object, sd->prim, sd->time, k0, k1, P_curve); - } - - float l = 1.0f; - tg = normalize_len(float4_to_float3(P_curve[1] - P_curve[0]), &l); - - P = P + D * t; - - float3 dif = P - float4_to_float3(P_curve[0]); - -# ifdef __UV__ - sd->u = dot(dif, tg) / l; - sd->v = 0.0f; -# endif - - if (flag & CURVE_KN_TRUETANGENTGNORMAL) { - sd->Ng = -(D - tg * dot(tg, D)); - sd->Ng = normalize(sd->Ng); - } - else { - float gd = isect->v; - - /* direction from inside to surface of curve */ - float denom = fmaxf(P_curve[0].w + sd->u * l * gd, 1e-8f); - sd->Ng = (dif - tg * sd->u * l) / denom; - - /* adjustment for changing radius */ - if (gd != 0.0f) { - sd->Ng = sd->Ng - gd * tg; - } - - sd->Ng = normalize(sd->Ng); - } - + /* Thick curves, compute normal using direction from inside the curve. + * This could be optimized by recording the normal in the intersection, + * however for Optix this would go beyond the size of the payload. */ + const float3 P_inside = float4_to_float3(catmull_rom_basis_eval(P_curve, isect->u)); + sd->Ng = normalize(P - P_inside); sd->N = sd->Ng; } # ifdef __DPDU__ /* dPdu/dPdv */ - sd->dPdu = tg; - sd->dPdv = cross(tg, sd->Ng); + sd->dPdu = dPdu; + sd->dPdv = cross(dPdu, sd->Ng); # endif if (isect->object != OBJECT_NONE) { @@ -867,7 +784,10 @@ ccl_device_inline float3 curve_refine(KernelGlobals *kg, P = transform_point(&tfm, P); } - return P; + sd->P = P; + + float4 curvedata = kernel_tex_fetch(__curves, sd->prim); + sd->shader = __float_as_int(curvedata.z); } #endif |