Cycles: port curve-ray intersection from Embree for use in Cycles GPU

This keeps render results compatible for combined CPU + GPU rendering.
Peformance and quality primitives is quite different than before. There
are now two options:

* Rounded Ribbon: render hair as flat ribbon with (fake) rounded normals, for
  fast rendering. Hair curves are subdivided with a fixed number of user
  specified subdivisions.

  This gives relatively good results, especially when used with the Principled
  Hair BSDF and hair viewed from a typical distance. There are artifacts when
  viewed closed up, though this was also the case with all previous primitives
  (but different ones).

* 3D Curve: render hair as 3D curve, for accurate results when viewing hair
  close up. This automatically subdivides the curve until it is smooth.

  This gives higher quality than any of the previous primitives, but does come
  at a performance cost and is somewhat slower than our previous Thick curves.

The main problem here is performance. For CPU and OpenCL rendering performance
seems usually quite close or better for similar quality results.

However for CUDA and Optix, performance of 3D curve intersection is problematic,
with e.g. 1.45x longer render time in Koro (though there is no equivalent quality
and rounded ribbons seem fine for that scene). Any help or ideas to optimize this
are welcome.

Ref T73778

Depends on D8012

Maniphest Tasks: T73778

Differential Revision: https://developer.blender.org/D8013

13 files changed, 758 insertions, 410 deletions

diff --git a/intern/cycles/blender/addon/properties.py b/intern/cycles/blender/addon/properties.py
index f0f7d24002f..ceaeda6e798 100644
--- a/intern/cycles/blender/addon/properties.py
+++ b/intern/cycles/blender/addon/properties.py
@@ -73,8 +73,8 @@ enum_panorama_types = (
 )
 
 enum_curve_shape = (
-    ('RIBBONS', "Ribbons", "Ignore thickness of each hair"),
-    ('THICK', "Thick", "Use thickness of hair when rendering"),
+    ('RIBBONS', "Rounded Ribbons", "Render hair as flat ribbon with rounded normals, for fast rendering"),
+    ('THICK', "3D Curves", "Render hair as 3D curve, for accurate results when viewing hair close up"),
 )
 
 enum_tile_order = (
@@ -1223,7 +1223,7 @@ class CyclesCurveRenderSettings(bpy.types.PropertyGroup):
         name="Shape",
         description="Form of hair",
         items=enum_curve_shape,
-        default='THICK',
+        default='RIBBONS',
     )
     use_curves: BoolProperty(
         name="Use Cycles Hair Rendering",
@@ -1234,7 +1234,7 @@ class CyclesCurveRenderSettings(bpy.types.PropertyGroup):
         name="Subdivisions",
         description="Number of subdivisions used in Cardinal curve intersection (power of 2)",
         min=0, max=24,
-        default=4,
+        default=2,
     )
 
     @classmethod
diff --git a/intern/cycles/blender/addon/ui.py b/intern/cycles/blender/addon/ui.py
index 0859a8a82b0..daab26d8ec0 100644
--- a/intern/cycles/blender/addon/ui.py
+++ b/intern/cycles/blender/addon/ui.py
@@ -407,7 +407,6 @@ class CYCLES_RENDER_PT_hair(CyclesButtonsPanel, Panel):
         col = layout.column()
         col.prop(ccscene, "shape", text="Shape")
         if ccscene.shape == 'RIBBONS':
-            # TODO: use for embree
             col.prop(ccscene, "subdivisions", text="Curve subdivisions")
 
 
diff --git a/intern/cycles/device/device.h b/intern/cycles/device/device.h
index 67828103394..69f22aeb35c 100644
--- a/intern/cycles/device/device.h
+++ b/intern/cycles/device/device.h
@@ -132,6 +132,7 @@ class DeviceRequestedFeatures {
 
   /* BVH/sampling kernel features. */
   bool use_hair;
+  bool use_hair_thick;
   bool use_object_motion;
   bool use_camera_motion;
 
@@ -178,6 +179,7 @@ class DeviceRequestedFeatures {
     max_nodes_group = 0;
     nodes_features = 0;
     use_hair = false;
+    use_hair_thick = false;
     use_object_motion = false;
     use_camera_motion = false;
     use_baking = false;
@@ -200,6 +202,7 @@ class DeviceRequestedFeatures {
              max_nodes_group == requested_features.max_nodes_group &&
              nodes_features == requested_features.nodes_features &&
              use_hair == requested_features.use_hair &&
+             use_hair_thick == requested_features.use_hair_thick &&
              use_object_motion == requested_features.use_object_motion &&
              use_camera_motion == requested_features.use_camera_motion &&
              use_baking == requested_features.use_baking &&
diff --git a/intern/cycles/device/device_optix.cpp b/intern/cycles/device/device_optix.cpp
index bc2aeb0ae90..7aab2c96db4 100644
--- a/intern/cycles/device/device_optix.cpp
+++ b/intern/cycles/device/device_optix.cpp
@@ -428,11 +428,20 @@ class OptiXDevice : public CUDADevice {
     group_descs[PG_HITS].hitgroup.entryFunctionNameAH = "__anyhit__kernel_optix_shadow_all_hit";
 
     if (requested_features.use_hair) {
-      // Add curve intersection programs
       group_descs[PG_HITD].hitgroup.moduleIS = optix_module;
-      group_descs[PG_HITD].hitgroup.entryFunctionNameIS = "__intersection__curve";
       group_descs[PG_HITS].hitgroup.moduleIS = optix_module;
-      group_descs[PG_HITS].hitgroup.entryFunctionNameIS = "__intersection__curve";
+
+      // Add curve intersection programs
+      if (requested_features.use_hair_thick) {
+        // Slower programs for thick hair since that also slows down ribbons.
+        // Ideally this should not be needed.
+        group_descs[PG_HITD].hitgroup.entryFunctionNameIS = "__intersection__curve_all";
+        group_descs[PG_HITS].hitgroup.entryFunctionNameIS = "__intersection__curve_all";
+      }
+      else {
+        group_descs[PG_HITD].hitgroup.entryFunctionNameIS = "__intersection__curve_ribbon";
+        group_descs[PG_HITS].hitgroup.entryFunctionNameIS = "__intersection__curve_ribbon";
+      }
     }
 
     if (requested_features.use_subsurface || requested_features.use_shader_raytrace) {
diff --git a/intern/cycles/kernel/closure/bsdf.h b/intern/cycles/kernel/closure/bsdf.h
index 4cc61e8ee71..6070fd983f5 100644
--- a/intern/cycles/kernel/closure/bsdf.h
+++ b/intern/cycles/kernel/closure/bsdf.h
@@ -119,13 +119,16 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
                                   differential3 *domega_in,
                                   float *pdf)
 {
+  /* For curves use the smooth normal, particularly for ribbons the geometric
+   * normal gives too much darkening otherwise. */
   int label;
+  const float3 Ng = (sd->type & PRIMITIVE_ALL_CURVE) ? sc->N : sd->Ng;
 
   switch (sc->type) {
     case CLOSURE_BSDF_DIFFUSE_ID:
     case CLOSURE_BSDF_BSSRDF_ID:
       label = bsdf_diffuse_sample(sc,
-                                  sd->Ng,
+                                  Ng,
                                   sd->I,
                                   sd->dI.dx,
                                   sd->dI.dy,
@@ -140,7 +143,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
 #ifdef __SVM__
     case CLOSURE_BSDF_OREN_NAYAR_ID:
       label = bsdf_oren_nayar_sample(sc,
-                                     sd->Ng,
+                                     Ng,
                                      sd->I,
                                      sd->dI.dx,
                                      sd->dI.dy,
@@ -155,7 +158,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
 #  ifdef __OSL__
     case CLOSURE_BSDF_PHONG_RAMP_ID:
       label = bsdf_phong_ramp_sample(sc,
-                                     sd->Ng,
+                                     Ng,
                                      sd->I,
                                      sd->dI.dx,
                                      sd->dI.dy,
@@ -169,7 +172,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_DIFFUSE_RAMP_ID:
       label = bsdf_diffuse_ramp_sample(sc,
-                                       sd->Ng,
+                                       Ng,
                                        sd->I,
                                        sd->dI.dx,
                                        sd->dI.dy,
@@ -184,7 +187,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
 #  endif
     case CLOSURE_BSDF_TRANSLUCENT_ID:
       label = bsdf_translucent_sample(sc,
-                                      sd->Ng,
+                                      Ng,
                                       sd->I,
                                       sd->dI.dx,
                                       sd->dI.dy,
@@ -198,7 +201,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_REFLECTION_ID:
       label = bsdf_reflection_sample(sc,
-                                     sd->Ng,
+                                     Ng,
                                      sd->I,
                                      sd->dI.dx,
                                      sd->dI.dy,
@@ -212,7 +215,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_REFRACTION_ID:
       label = bsdf_refraction_sample(sc,
-                                     sd->Ng,
+                                     Ng,
                                      sd->I,
                                      sd->dI.dx,
                                      sd->dI.dy,
@@ -226,7 +229,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_TRANSPARENT_ID:
       label = bsdf_transparent_sample(sc,
-                                      sd->Ng,
+                                      Ng,
                                       sd->I,
                                       sd->dI.dx,
                                       sd->dI.dy,
@@ -244,7 +247,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
     case CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID:
       label = bsdf_microfacet_ggx_sample(kg,
                                          sc,
-                                         sd->Ng,
+                                         Ng,
                                          sd->I,
                                          sd->dI.dx,
                                          sd->dI.dy,
@@ -260,7 +263,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
     case CLOSURE_BSDF_MICROFACET_MULTI_GGX_FRESNEL_ID:
       label = bsdf_microfacet_multi_ggx_sample(kg,
                                                sc,
-                                               sd->Ng,
+                                               Ng,
                                                sd->I,
                                                sd->dI.dx,
                                                sd->dI.dy,
@@ -277,7 +280,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
     case CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_FRESNEL_ID:
       label = bsdf_microfacet_multi_ggx_glass_sample(kg,
                                                      sc,
-                                                     sd->Ng,
+                                                     Ng,
                                                      sd->I,
                                                      sd->dI.dx,
                                                      sd->dI.dy,
@@ -294,7 +297,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
     case CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID:
       label = bsdf_microfacet_beckmann_sample(kg,
                                               sc,
-                                              sd->Ng,
+                                              Ng,
                                               sd->I,
                                               sd->dI.dx,
                                               sd->dI.dy,
@@ -308,7 +311,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_ASHIKHMIN_SHIRLEY_ID:
       label = bsdf_ashikhmin_shirley_sample(sc,
-                                            sd->Ng,
+                                            Ng,
                                             sd->I,
                                             sd->dI.dx,
                                             sd->dI.dy,
@@ -322,7 +325,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_ASHIKHMIN_VELVET_ID:
       label = bsdf_ashikhmin_velvet_sample(sc,
-                                           sd->Ng,
+                                           Ng,
                                            sd->I,
                                            sd->dI.dx,
                                            sd->dI.dy,
@@ -336,7 +339,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_DIFFUSE_TOON_ID:
       label = bsdf_diffuse_toon_sample(sc,
-                                       sd->Ng,
+                                       Ng,
                                        sd->I,
                                        sd->dI.dx,
                                        sd->dI.dy,
@@ -350,7 +353,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_GLOSSY_TOON_ID:
       label = bsdf_glossy_toon_sample(sc,
-                                      sd->Ng,
+                                      Ng,
                                       sd->I,
                                       sd->dI.dx,
                                       sd->dI.dy,
@@ -364,7 +367,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_HAIR_REFLECTION_ID:
       label = bsdf_hair_reflection_sample(sc,
-                                          sd->Ng,
+                                          Ng,
                                           sd->I,
                                           sd->dI.dx,
                                           sd->dI.dy,
@@ -378,7 +381,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_HAIR_TRANSMISSION_ID:
       label = bsdf_hair_transmission_sample(sc,
-                                            sd->Ng,
+                                            Ng,
                                             sd->I,
                                             sd->dI.dx,
                                             sd->dI.dy,
@@ -398,7 +401,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
     case CLOSURE_BSDF_PRINCIPLED_DIFFUSE_ID:
     case CLOSURE_BSDF_BSSRDF_PRINCIPLED_ID:
       label = bsdf_principled_diffuse_sample(sc,
-                                             sd->Ng,
+                                             Ng,
                                              sd->I,
                                              sd->dI.dx,
                                              sd->dI.dy,
@@ -412,7 +415,7 @@ ccl_device_inline int bsdf_sample(KernelGlobals *kg,
       break;
     case CLOSURE_BSDF_PRINCIPLED_SHEEN_ID:
       label = bsdf_principled_sheen_sample(sc,
-                                           sd->Ng,
+                                           Ng,
                                            sd->I,
                                            sd->dI.dx,
                                            sd->dI.dy,
@@ -485,9 +488,12 @@ ccl_device_inline
               const float3 omega_in,
               float *pdf)
 {
+  /* For curves use the smooth normal, particularly for ribbons the geometric
+   * normal gives too much darkening otherwise. */
+  const float3 Ng = (sd->type & PRIMITIVE_ALL_CURVE) ? sd->N : sd->Ng;
   float3 eval;
 
-  if (dot(sd->Ng, omega_in) >= 0.0f) {
+  if (dot(Ng, omega_in) >= 0.0f) {
     switch (sc->type) {
       case CLOSURE_BSDF_DIFFUSE_ID:
       case CLOSURE_BSDF_BSSRDF_ID:
diff --git a/intern/cycles/kernel/closure/bsdf_hair_principled.h b/intern/cycles/kernel/closure/bsdf_hair_principled.h
index f78bbeb5d9d..18934d86cc8 100644
--- a/intern/cycles/kernel/closure/bsdf_hair_principled.h
+++ b/intern/cycles/kernel/closure/bsdf_hair_principled.h
@@ -206,9 +206,6 @@ ccl_device int bsdf_principled_hair_setup(ShaderData *sd, PrincipledHairBSDF *bs
   float3 X = safe_normalize(sd->dPdu);
   float3 Y = safe_normalize(cross(X, sd->I));
   float3 Z = safe_normalize(cross(X, Y));
-  /* TODO: the solution below works where sd->Ng is the normal
-   * pointing from the center of the curve to the shading point.
-   * It doesn't work for triangles, see https://developer.blender.org/T43625 */
 
   /* h -1..0..1 means the rays goes from grazing the hair, to hitting it at
    * the center, to grazing the other edge. This is the sine of the angle
diff --git a/intern/cycles/kernel/geom/geom_curve.h b/intern/cycles/kernel/geom/geom_curve.h
index 928cad58452..d2ac2d60435 100644
--- a/intern/cycles/kernel/geom/geom_curve.h
+++ b/intern/cycles/kernel/geom/geom_curve.h
@@ -23,33 +23,6 @@ CCL_NAMESPACE_BEGIN
 
 #ifdef __HAIR__
 
-/* Interpolation of curve geometry */
-
-ccl_device_inline float3 curvetangent(float t, float3 p0, float3 p1, float3 p2, float3 p3)
-{
-  float fc = 0.71f;
-  float data[4];
-  float t2 = t * t;
-  data[0] = -3.0f * fc * t2 + 4.0f * fc * t - fc;
-  data[1] = 3.0f * (2.0f - fc) * t2 + 2.0f * (fc - 3.0f) * t;
-  data[2] = 3.0f * (fc - 2.0f) * t2 + 2.0f * (3.0f - 2.0f * fc) * t + fc;
-  data[3] = 3.0f * fc * t2 - 2.0f * fc * t;
-  return data[0] * p0 + data[1] * p1 + data[2] * p2 + data[3] * p3;
-}
-
-ccl_device_inline float3 curvepoint(float t, float3 p0, float3 p1, float3 p2, float3 p3)
-{
-  float data[4];
-  float fc = 0.71f;
-  float t2 = t * t;
-  float t3 = t2 * t;
-  data[0] = -fc * t3 + 2.0f * fc * t2 - fc * t;
-  data[1] = (2.0f - fc) * t3 + (fc - 3.0f) * t2 + 1.0f;
-  data[2] = (fc - 2.0f) * t3 + (3.0f - 2.0f * fc) * t2 + fc * t;
-  data[3] = fc * t3 - fc * t2;
-  return data[0] * p0 + data[1] * p1 + data[2] * p2 + data[3] * p3;
-}
-
 /* Reading attributes on various curve elements */
 
 ccl_device float curve_attribute_float(
@@ -243,7 +216,7 @@ ccl_device float curve_thickness(KernelGlobals *kg, ShaderData *sd)
       P_curve[1] = kernel_tex_fetch(__curve_keys, k1);
     }
     else {
-      motion_curve_keys(kg, sd->object, sd->prim, sd->time, k0, k1, P_curve);
+      motion_curve_keys_linear(kg, sd->object, sd->prim, sd->time, k0, k1, P_curve);
     }
 
     r = (P_curve[1].w - P_curve[0].w) * sd->u + P_curve[0].w;
diff --git a/intern/cycles/kernel/geom/geom_curve_intersect.h b/intern/cycles/kernel/geom/geom_curve_intersect.h
index 87ed0bf201f..c4a614ab676 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,337 +17,685 @@
 
 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 betwee 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__
 
-/* On CPU pass P and dir by reference to aligned vector. */
-ccl_device_forceinline bool 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)
+/* Catmull-rom curve evaluation. */
+
+ccl_device_inline float4 catmull_rom_basis_eval(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 = -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);
+}
 
-#  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) {
+ccl_device_inline float4 catmull_rom_basis_derivative(const float4 curve[4], float u)
+{
+  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);
+}
+
+ccl_device_inline float4 catmull_rom_basis_derivative2(const float4 curve[4], float u)
+{
+
+  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);
+}
+
+/* Thick Curve */
+
+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));
+}
+
+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 {
+      *t_o = make_float2(-FLT_MAX, FLT_MAX);
       return false;
     }
   }
+
+  /* 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. */
+  {
+    *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;
+  }
+
+  /* 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;
+  }
+
+  *t_o = make_float2(t0, t1);
+
+  return true;
+}
+
+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);
+}
+
+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 propertional 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 */
+      }
+
+      /* 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;
+      }
+
+      /* Record intersection. */
+      isect->t = t;
+      isect->u = u;
+      isect->v = 0.0f;
+
+      return true;
+    }
+  }
+  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
+      /* On the GPU appears to be a little faster if always enabled. */
+      (void)valid_inner;
+
+      const bool unstable0 = true;
+      const bool unstable1 = true;
 #  endif
 
-  int segment = PRIMITIVE_UNPACK_SEGMENT(type);
-  float epsilon = 0.0f;
-  float r_st, r_en;
+      /* 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;
+      }
 
-  int depth = kernel_data.curve.subdivisions;
-  int flags = kernel_data.curve.curveflags;
-  int prim = kernel_tex_fetch(__prim_index, curveAddr);
+      /* 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);
+        }
+        else {
+          recurse = true;
+        }
+      }
 
-  float3 curve_coef[4];
+      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;
+        }
+      }
 
-  /* curve Intersection check */
-  /* obtain curve parameters */
-  {
-    /* 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);
+      if (recurse) {
+        stack[depth].u0 = u0;
+        stack[depth].u1 = u1;
+        stack[depth].i = i + 1;
+        depth++;
+
+        u0 = vu0;
+        u1 = vu1;
+        i = -1;
+      }
+    }
+
+    if (depth > 0) {
+      depth--;
+      u0 = stack[depth].u0;
+      u1 = stack[depth].u1;
+      i = stack[depth].i;
     }
     else {
-      int fobject = (object == OBJECT_NONE) ? kernel_tex_fetch(__prim_object, curveAddr) : object;
-      motion_curve_keys(kg, fobject, prim, time, ka, k0, k1, kb, P_curve);
+      break;
     }
-
-    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;
   }
 
-  float r_curr = max(r_st, r_en);
-
-  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 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;
 
-  /* 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)
+  /* Perform depth test? */
+  const float t = rcpDen * dot(v0, Ng);
+  if (!(0.0f <= t && t <= ray_tfar)) {
     return false;
+  }
 
-  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)
+  /* Avoid division by 0? */
+  if (!(den != 0.0f)) {
     return false;
+  }
 
-  /* 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);
-
-    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);
-
-    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]);
+  /* 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;
     }
-    if (zextrem[2] >= i_st && zextrem[2] <= i_en) {
-      bminz = min(bminz, zextrem[3]);
-      bmaxz = max(bmaxz, zextrem[3]);
+
+    /* 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;
+      }
     }
 
-    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);
+    p0 = p1;
+    wn0 = wn1;
+  }
+  return false;
+}
 
-    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;
+ccl_device_forceinline bool curve_intersect(KernelGlobals *kg,
+                                            Intersection *isect,
+                                            const float3 P,
+                                            const float3 dir,
+                                            uint visibility,
+                                            int object,
+                                            int curveAddr,
+                                            float time,
+                                            int type)
+{
+  const bool is_curve_primitive = (type & PRIMITIVE_CURVE);
+
+#  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;
     }
-    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);
-        float w = tg.x * tg.x + tg.y * tg.y;
-        if (w == 0) {
-          tree++;
-          level = tree & -tree;
-          continue;
-        }
-        w = -(p_st.x * tg.x + p_st.y * tg.y) / w;
-        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;
-        }
+  }
+#  endif
 
-        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;
-        }
+  int segment = PRIMITIVE_UNPACK_SEGMENT(type);
+  int prim = kernel_tex_fetch(__prim_index, curveAddr);
 
-        t = p_curr.z;
-      }
-      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;
-        }
+  float4 v00 = kernel_tex_fetch(__curves, prim);
 
-        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 (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;
-        }
+  int k0 = __float_as_int(v00.x) + segment;
+  int k1 = k0 + 1;
 
-        float w = (zcentre + (tg.z * correction)) * invl;
-        w = saturate(w);
-        /* compute u on the curve segment */
-        u = i_st * (1 - w) + i_en * w;
-      }
-      /* we found a new intersection */
+  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 curve[4];
+  if (is_curve_primitive) {
+    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, ka, k0, k1, kb, curve);
+  }
 
 #  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)
+  if (!(kernel_tex_fetch(__prim_visibility, curveAddr) & visibility)) {
+    return false;
+  }
 #  endif
-      {
-        /* record intersection */
-        isect->t = t;
-        isect->u = u;
-        isect->v = gd;
-        isect->prim = curveAddr;
-        isect->object = object;
-        isect->type = type;
-        hit = true;
-      }
 
-      tree++;
-      level = tree & -tree;
-    }
-    else {
-      /* split the curve into two curves and process */
-      level = level >> 1;
+  const bool use_ribbon = (kernel_data.curve.curveflags & CURVE_KN_RIBBONS) != 0;
+  if (use_ribbon) {
+    /* todo: adaptive number of subdivisions could help performance here. */
+    const int subdivisions = kernel_data.curve.subdivisions;
+    if (ribbon_intersect(P, dir, isect->t, subdivisions, curve, isect)) {
+      isect->prim = curveAddr;
+      isect->object = object;
+      isect->type = type;
+      return true;
     }
+
+    return false;
   }
+  else {
+    if (curve_intersect_recursive(P, dir, curve, isect)) {
+      isect->prim = curveAddr;
+      isect->object = object;
+      isect->type = type;
+      return true;
+    }
 
-  return hit;
+    return false;
+  }
 }
 
-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)
 {
   float t = isect->t;
   float3 P = ray->P;
@@ -383,53 +734,45 @@ ccl_device_inline float3 curve_refine(KernelGlobals *kg,
     motion_curve_keys(kg, sd->object, sd->prim, sd->time, ka, k0, k1, kb, P_curve);
   }
 
-  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]);
+  sd->u = isect->u;
+  sd->v = isect->v;
 
   P = P + D * t;
 
-  sd->u = isect->u;
-  sd->v = 0.0f;
-
-  float3 tg = normalize(curvetangent(isect->u, p[0], p[1], p[2], p[3]));
+  const float4 dPdu4 = catmull_rom_basis_derivative(P_curve, isect->u);
+  const float3 dPdu = float4_to_float3(dPdu4);
 
   if (kernel_data.curve.curveflags & CURVE_KN_RIBBONS) {
-    sd->Ng = normalize(-(D - tg * (dot(tg, D))));
+    /* 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);
+
+    sd->N = normalize(sine * bitangent - cosine * normalize(cross(tangent, bitangent)));
+    sd->Ng = -D;
+
+#  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
   }
   else {
-#  ifdef __EMBREE__
-    if (kernel_data.bvh.scene) {
-      sd->Ng = normalize(isect->Ng);
-    }
-    else
-#  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);
-      }
-    }
+    /* 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;
   }
 
-  /* todo: sometimes the normal is still so that this is detected as
-   * backfacing even if cull backfaces is enabled */
-
-  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) {
@@ -442,7 +785,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
diff --git a/intern/cycles/kernel/geom/geom_motion_curve.h b/intern/cycles/kernel/geom/geom_motion_curve.h
index 0e7a05eaac2..0f66f4af755 100644
--- a/intern/cycles/kernel/geom/geom_motion_curve.h
+++ b/intern/cycles/kernel/geom/geom_motion_curve.h
@@ -50,14 +50,14 @@ ccl_device_inline int find_attribute_curve_motion(KernelGlobals *kg,
   return (attr_map.y == ATTR_ELEMENT_NONE) ? (int)ATTR_STD_NOT_FOUND : (int)attr_map.z;
 }
 
-ccl_device_inline void motion_curve_keys_for_step(KernelGlobals *kg,
-                                                  int offset,
-                                                  int numkeys,
-                                                  int numsteps,
-                                                  int step,
-                                                  int k0,
-                                                  int k1,
-                                                  float4 keys[2])
+ccl_device_inline void motion_curve_keys_for_step_linear(KernelGlobals *kg,
+                                                         int offset,
+                                                         int numkeys,
+                                                         int numsteps,
+                                                         int step,
+                                                         int k0,
+                                                         int k1,
+                                                         float4 keys[2])
 {
   if (step == numsteps) {
     /* center step: regular key location */
@@ -77,7 +77,7 @@ ccl_device_inline void motion_curve_keys_for_step(KernelGlobals *kg,
 }
 
 /* return 2 curve key locations */
-ccl_device_inline void motion_curve_keys(
+ccl_device_inline void motion_curve_keys_linear(
     KernelGlobals *kg, int object, int prim, float time, int k0, int k1, float4 keys[2])
 {
   /* get motion info */
@@ -97,8 +97,8 @@ ccl_device_inline void motion_curve_keys(
   /* fetch key coordinates */
   float4 next_keys[2];
 
-  motion_curve_keys_for_step(kg, offset, numkeys, numsteps, step, k0, k1, keys);
-  motion_curve_keys_for_step(kg, offset, numkeys, numsteps, step + 1, k0, k1, next_keys);
+  motion_curve_keys_for_step_linear(kg, offset, numkeys, numsteps, step, k0, k1, keys);
+  motion_curve_keys_for_step_linear(kg, offset, numkeys, numsteps, step + 1, k0, k1, next_keys);
 
   /* interpolate between steps */
   keys[0] = (1.0f - t) * keys[0] + t * next_keys[0];
diff --git a/intern/cycles/kernel/kernel_shader.h b/intern/cycles/kernel/kernel_shader.h
index 535253bc762..3d9f787f267 100644
--- a/intern/cycles/kernel/kernel_shader.h
+++ b/intern/cycles/kernel/kernel_shader.h
@@ -86,10 +86,7 @@ ccl_device_noinline
 #ifdef __HAIR__
   if (sd->type & PRIMITIVE_ALL_CURVE) {
     /* curve */
-    float4 curvedata = kernel_tex_fetch(__curves, sd->prim);
-
-    sd->shader = __float_as_int(curvedata.z);
-    sd->P = curve_refine(kg, sd, isect, ray);
+    curve_shader_setup(kg, sd, isect, ray);
   }
   else
 #endif
diff --git a/intern/cycles/kernel/kernels/optix/kernel_optix.cu b/intern/cycles/kernel/kernels/optix/kernel_optix.cu
index 59fece73b3c..c730d952ed4 100644
--- a/intern/cycles/kernel/kernels/optix/kernel_optix.cu
+++ b/intern/cycles/kernel/kernels/optix/kernel_optix.cu
@@ -256,11 +256,9 @@ extern "C" __global__ void __closesthit__kernel_optix_hit()
 }
 
 #ifdef __HAIR__
-extern "C" __global__ void __intersection__curve()
+ccl_device_inline void optix_intersection_curve(const uint prim, const uint type)
 {
-  const uint prim = optixGetPrimitiveIndex();
   const uint object = get_object_id<true>();
-  const uint type = kernel_tex_fetch(__prim_type, prim);
   const uint visibility = optixGetPayload_4();
 
   float3 P = optixGetObjectRayOrigin();
@@ -288,6 +286,24 @@ extern "C" __global__ void __intersection__curve()
                             __float_as_int(isect.u),   // Attribute_0
                             __float_as_int(isect.v));  // Attribute_1
   }
+
+}
+
+extern "C" __global__ void __intersection__curve_ribbon()
+{
+  const uint prim = optixGetPrimitiveIndex();
+  const uint type = kernel_tex_fetch(__prim_type, prim);
+
+  if (type & (PRIMITIVE_CURVE_RIBBON | PRIMITIVE_MOTION_CURVE_RIBBON)) {
+    optix_intersection_curve(prim, type);
+  }
+}
+
+extern "C" __global__ void __intersection__curve_all()
+{
+  const uint prim = optixGetPrimitiveIndex();
+  const uint type = kernel_tex_fetch(__prim_type, prim);
+  optix_intersection_curve(prim, type);
 }
 #endif
 
diff --git a/intern/cycles/render/curves.cpp b/intern/cycles/render/curves.cpp
index 0d2be71eed8..7c9bcaa2549 100644
--- a/intern/cycles/render/curves.cpp
+++ b/intern/cycles/render/curves.cpp
@@ -36,13 +36,12 @@ void curvebounds(float *lower, float *upper, float3 *p, int dim)
   float *p2 = &p[2].x;
   float *p3 = &p[3].x;
 
-  float fc = 0.71f;
+  /* Catmull-Rom weights. */
   float curve_coef[4];
   curve_coef[0] = p1[dim];
-  curve_coef[1] = -fc * p0[dim] + fc * p2[dim];
-  curve_coef[2] = 2.0f * fc * p0[dim] + (fc - 3.0f) * p1[dim] + (3.0f - 2.0f * fc) * p2[dim] -
-                  fc * p3[dim];
-  curve_coef[3] = -fc * p0[dim] + (2.0f - fc) * p1[dim] + (fc - 2.0f) * p2[dim] + fc * p3[dim];
+  curve_coef[1] = 0.5f * (-p0[dim] + p2[dim]);
+  curve_coef[2] = 0.5f * (2 * p0[dim] - 5 * p1[dim] + 4 * p2[dim] - p3[dim]);
+  curve_coef[3] = 0.5f * (-p0[dim] + 3 * p1[dim] - 3 * p2[dim] + p3[dim]);
 
   float discroot = curve_coef[2] * curve_coef[2] - 3 * curve_coef[3] * curve_coef[1];
   float ta = -1.0f;
@@ -115,7 +114,8 @@ void CurveSystemManager::device_update(Device *device,
       kcurve->curveflags |= CURVE_KN_RIBBONS;
     }
 
-    kcurve->subdivisions = subdivisions;
+    /* Matching the tesselation rate limit in Embree. */
+    kcurve->subdivisions = clamp(1 << subdivisions, 1, 16);
   }
 
   if (progress.get_cancel())
diff --git a/intern/cycles/render/session.cpp b/intern/cycles/render/session.cpp
index f5bfebbaf78..6caa686847e 100644
--- a/intern/cycles/render/session.cpp
+++ b/intern/cycles/render/session.cpp
@@ -21,6 +21,7 @@
 #include "render/bake.h"
 #include "render/buffers.h"
 #include "render/camera.h"
+#include "render/curves.h"
 #include "render/graph.h"
 #include "render/integrator.h"
 #include "render/light.h"
@@ -773,6 +774,7 @@ DeviceRequestedFeatures Session::get_requested_device_features()
    */
   bool use_motion = scene->need_motion() == Scene::MotionType::MOTION_BLUR;
   requested_features.use_hair = false;
+  requested_features.use_hair_thick = (scene->curve_system_manager->curve_shape == CURVE_THICK);
   requested_features.use_object_motion = false;
   requested_features.use_camera_motion = use_motion && scene->camera->use_motion();
   foreach (Object *object, scene->objects) {