Merge branch 'master' into blender2.8

Conflicts:
	source/blender/editors/space_action/action_draw.c

34 files changed, 1162 insertions, 404 deletions

diff --git a/intern/cycles/CMakeLists.txt b/intern/cycles/CMakeLists.txt
index 97854a88e84..79c1c3e3e82 100644
--- a/intern/cycles/CMakeLists.txt
+++ b/intern/cycles/CMakeLists.txt
@@ -74,7 +74,6 @@ elseif(CMAKE_COMPILER_IS_GNUCC)
 	if(CXX_HAS_AVX2)
 		set(CYCLES_AVX2_KERNEL_FLAGS "-ffast-math -msse -msse2 -msse3 -mssse3 -msse4.1 -mavx -mavx2 -mfma -mlzcnt -mbmi -mbmi2 -mf16c -mfpmath=sse")
 	endif()
-	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ffast-math")
 elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
 	check_cxx_compiler_flag(-msse CXX_HAS_SSE)
 	check_cxx_compiler_flag(-mavx CXX_HAS_AVX)
@@ -90,7 +89,6 @@ elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
 	if(CXX_HAS_AVX2)
 		set(CYCLES_AVX2_KERNEL_FLAGS "-ffast-math -msse -msse2 -msse3 -mssse3 -msse4.1 -mavx -mavx2 -mfma -mlzcnt -mbmi -mbmi2 -mf16c")
 	endif()
-	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ffast-math")
 endif()
 
 if(CXX_HAS_SSE)
diff --git a/intern/cycles/blender/addon/__init__.py b/intern/cycles/blender/addon/__init__.py
index 29388317873..1fc3758ad4d 100644
--- a/intern/cycles/blender/addon/__init__.py
+++ b/intern/cycles/blender/addon/__init__.py
@@ -28,6 +28,20 @@ bl_info = {
     "support": 'OFFICIAL',
     "category": "Render"}
 
+# Support 'reload' case.
+if "bpy" in locals():
+    import importlib
+    if "engine" in locals():
+        importlib.reload(engine)
+    if "version_update" in locals():
+        importlib.reload(version_update)
+    if "ui" in locals():
+        importlib.reload(ui)
+    if "properties" in locals():
+        importlib.reload(properties)
+    if "presets" in locals():
+        importlib.reload(presets)
+
 import bpy
 
 from . import (
diff --git a/intern/cycles/blender/addon/properties.py b/intern/cycles/blender/addon/properties.py
index 3616b13e751..802b9b76c5d 100644
--- a/intern/cycles/blender/addon/properties.py
+++ b/intern/cycles/blender/addon/properties.py
@@ -528,6 +528,12 @@ class CyclesRenderSettings(bpy.types.PropertyGroup):
                 description="Use special type BVH optimized for hair (uses more ram but renders faster)",
                 default=True,
                 )
+        cls.debug_bvh_time_steps = IntProperty(
+                name="BVH Time Steps",
+                description="Split BVH primitives by this number of time steps to speed up render time in cost of memory",
+                default=0,
+                min=0, max=16,
+                )
         cls.tile_order = EnumProperty(
                 name="Tile Order",
                 description="Tile order for rendering",
diff --git a/intern/cycles/blender/addon/ui.py b/intern/cycles/blender/addon/ui.py
index a573fa1ce22..7a59698e15e 100644
--- a/intern/cycles/blender/addon/ui.py
+++ b/intern/cycles/blender/addon/ui.py
@@ -432,6 +432,10 @@ class CyclesRender_PT_performance(CyclesButtonsPanel, Panel):
         col.prop(cscene, "debug_use_spatial_splits")
         col.prop(cscene, "debug_use_hair_bvh")
 
+        row = col.row()
+        row.active = not cscene.debug_use_spatial_splits
+        row.prop(cscene, "debug_bvh_time_steps")
+
 
 class CyclesRender_PT_layer_options(CyclesButtonsPanel, Panel):
     bl_label = "Layer"
@@ -787,10 +791,13 @@ class CyclesObject_PT_cycles_settings(CyclesButtonsPanel, Panel):
         col = layout.column()
         col.label(text="Performance:")
         row = col.row()
-        row.active = scene.render.use_simplify and cscene.use_camera_cull
-        row.prop(cob, "use_camera_cull")
-        row.active = scene.render.use_simplify and cscene.use_distance_cull
-        row.prop(cob, "use_distance_cull")
+        sub = row.row()
+        sub.active = scene.render.use_simplify and cscene.use_camera_cull
+        sub.prop(cob, "use_camera_cull")
+
+        sub = row.row()
+        sub.active = scene.render.use_simplify and cscene.use_distance_cull
+        sub.prop(cob, "use_distance_cull")
 
 
 class CYCLES_OT_use_shading_nodes(Operator):
diff --git a/intern/cycles/blender/blender_sync.cpp b/intern/cycles/blender/blender_sync.cpp
index 38b2ce19e8a..f8f2303ec76 100644
--- a/intern/cycles/blender/blender_sync.cpp
+++ b/intern/cycles/blender/blender_sync.cpp
@@ -498,6 +498,7 @@ SceneParams BlenderSync::get_scene_params(BL::Scene& b_scene,
 
 	params.use_bvh_spatial_split = RNA_boolean_get(&cscene, "debug_use_spatial_splits");
 	params.use_bvh_unaligned_nodes = RNA_boolean_get(&cscene, "debug_use_hair_bvh");
+	params.num_bvh_time_steps = RNA_int_get(&cscene, "debug_bvh_time_steps");
 
 	if(background && params.shadingsystem != SHADINGSYSTEM_OSL)
 		params.persistent_data = r.use_persistent_data();
diff --git a/intern/cycles/bvh/bvh.cpp b/intern/cycles/bvh/bvh.cpp
index 4851de5b481..874a4246d1d 100644
--- a/intern/cycles/bvh/bvh.cpp
+++ b/intern/cycles/bvh/bvh.cpp
@@ -845,6 +845,8 @@ void QBVH::pack_aligned_inner(const BVHStackEntry& e,
 	                  bounds,
 	                  child,
 	                  e.node->m_visibility,
+	                  e.node->m_time_from,
+	                  e.node->m_time_to,
 	                  num);
 }
 
@@ -852,12 +854,17 @@ void QBVH::pack_aligned_node(int idx,
                              const BoundBox *bounds,
                              const int *child,
                              const uint visibility,
+                             const float time_from,
+                             const float time_to,
                              const int num)
 {
 	float4 data[BVH_QNODE_SIZE];
 	memset(data, 0, sizeof(data));
 
 	data[0].x = __uint_as_float(visibility & ~PATH_RAY_NODE_UNALIGNED);
+	data[0].y = time_from;
+	data[0].z = time_to;
+
 	for(int i = 0; i < num; i++) {
 		float3 bb_min = bounds[i].min;
 		float3 bb_max = bounds[i].max;
@@ -908,6 +915,8 @@ void QBVH::pack_unaligned_inner(const BVHStackEntry& e,
 	                    bounds,
 	                    child,
 	                    e.node->m_visibility,
+	                    e.node->m_time_from,
+	                    e.node->m_time_to,
 	                    num);
 }
 
@@ -916,12 +925,16 @@ void QBVH::pack_unaligned_node(int idx,
                                const BoundBox *bounds,
                                const int *child,
                                const uint visibility,
+                               const float time_from,
+                               const float time_to,
                                const int num)
 {
 	float4 data[BVH_UNALIGNED_QNODE_SIZE];
 	memset(data, 0, sizeof(data));
 
 	data[0].x = __uint_as_float(visibility | PATH_RAY_NODE_UNALIGNED);
+	data[0].y = time_from;
+	data[0].z = time_to;
 
 	for(int i = 0; i < num; i++) {
 		Transform space = BVHUnaligned::compute_node_transform(
@@ -1207,6 +1220,8 @@ void QBVH::refit_node(int idx, bool leaf, BoundBox& bbox, uint& visibility)
 			                    child_bbox,
 			                    &c[0],
 			                    visibility,
+			                    0.0f,
+			                    1.0f,
 			                    4);
 		}
 		else {
@@ -1214,6 +1229,8 @@ void QBVH::refit_node(int idx, bool leaf, BoundBox& bbox, uint& visibility)
 			                  child_bbox,
 			                  &c[0],
 			                  visibility,
+			                  0.0f,
+			                  1.0f,
 			                  4);
 		}
 	}
diff --git a/intern/cycles/bvh/bvh.h b/intern/cycles/bvh/bvh.h
index f8fb3b568ca..35f4d305883 100644
--- a/intern/cycles/bvh/bvh.h
+++ b/intern/cycles/bvh/bvh.h
@@ -175,6 +175,8 @@ protected:
 	                       const BoundBox *bounds,
 	                       const int *child,
 	                       const uint visibility,
+	                       const float time_from,
+	                       const float time_to,
 	                       const int num);
 
 	void pack_unaligned_inner(const BVHStackEntry& e,
@@ -185,6 +187,8 @@ protected:
 	                         const BoundBox *bounds,
 	                         const int *child,
 	                         const uint visibility,
+	                         const float time_from,
+	                         const float time_to,
 	                         const int num);
 
 	/* refit */
diff --git a/intern/cycles/bvh/bvh_build.cpp b/intern/cycles/bvh/bvh_build.cpp
index 8cf1495b33d..1ce3a754460 100644
--- a/intern/cycles/bvh/bvh_build.cpp
+++ b/intern/cycles/bvh/bvh_build.cpp
@@ -112,84 +112,237 @@ BVHBuild::~BVHBuild()
 
 /* Adding References */
 
-void BVHBuild::add_reference_mesh(BoundBox& root, BoundBox& center, Mesh *mesh, int i)
+void BVHBuild::add_reference_triangles(BoundBox& root, BoundBox& center, Mesh *mesh, int i)
 {
-	if(params.primitive_mask & PRIMITIVE_ALL_TRIANGLE) {
-		Attribute *attr_mP = NULL;
-
-		if(mesh->has_motion_blur())
-			attr_mP = mesh->attributes.find(ATTR_STD_MOTION_VERTEX_POSITION);
-
-		size_t num_triangles = mesh->num_triangles();
-		for(uint j = 0; j < num_triangles; j++) {
-			Mesh::Triangle t = mesh->get_triangle(j);
+	const Attribute *attr_mP = NULL;
+	if(mesh->has_motion_blur()) {
+		attr_mP = mesh->attributes.find(ATTR_STD_MOTION_VERTEX_POSITION);
+	}
+	const size_t num_triangles = mesh->num_triangles();
+	for(uint j = 0; j < num_triangles; j++) {
+		Mesh::Triangle t = mesh->get_triangle(j);
+		const float3 *verts = &mesh->verts[0];
+		if(attr_mP == NULL) {
 			BoundBox bounds = BoundBox::empty;
-			PrimitiveType type = PRIMITIVE_TRIANGLE;
-
-			t.bounds_grow(&mesh->verts[0], bounds);
-
-			/* motion triangles */
-			if(attr_mP) {
-				const size_t mesh_size = mesh->verts.size();
-				const size_t num_steps = mesh->motion_steps - 1;
-				const float3 *vert_steps = attr_mP->data_float3();
-
-				for(size_t step = 0; step < num_steps; step++) {
-					t.bounds_grow(vert_steps + step*mesh_size, bounds);
-				}
-
-				type = PRIMITIVE_MOTION_TRIANGLE;
+			t.bounds_grow(verts, bounds);
+			if(bounds.valid()) {
+				references.push_back(BVHReference(bounds,
+				                                  j,
+				                                  i,
+				                                  PRIMITIVE_TRIANGLE));
+				root.grow(bounds);
+				center.grow(bounds.center2());
+			}
+		}
+		else if(params.num_motion_triangle_steps == 0 || params.use_spatial_split) {
+			/* Motion triangles, simple case: single node for the whole
+			 * primitive. Lowest memory footprint and faster BVH build but
+			 * least optimal ray-tracing.
+			 */
+			/* TODO(sergey): Support motion steps for spatially split BVH. */
+			const size_t num_verts = mesh->verts.size();
+			const size_t num_steps = mesh->motion_steps;
+			const float3 *vert_steps = attr_mP->data_float3();
+			BoundBox bounds = BoundBox::empty;
+			t.bounds_grow(verts, bounds);
+			for(size_t step = 0; step < num_steps - 1; step++) {
+				t.bounds_grow(vert_steps + step*num_verts, bounds);
 			}
-
 			if(bounds.valid()) {
-				references.push_back(BVHReference(bounds, j, i, type));
+				references.push_back(
+				        BVHReference(bounds,
+				                     j,
+				                     i,
+				                     PRIMITIVE_MOTION_TRIANGLE));
 				root.grow(bounds);
 				center.grow(bounds.center2());
 			}
 		}
+		else {
+			/* Motion triangles, trace optimized case:  we split triangle
+			 * primitives into separate nodes for each of the time steps.
+			 * This way we minimize overlap of neighbor curve primitives.
+			 */
+			const int num_bvh_steps = params.num_motion_curve_steps * 2 + 1;
+			const float num_bvh_steps_inv_1 = 1.0f / (num_bvh_steps - 1);
+			const size_t num_verts = mesh->verts.size();
+			const size_t num_steps = mesh->motion_steps;
+			const float3 *vert_steps = attr_mP->data_float3();
+			/* Calculate bounding box of the previous time step.
+			 * Will be reused later to avoid duplicated work on
+			 * calculating BVH time step boundbox.
+			 */
+			float3 prev_verts[3];
+			t.motion_verts(verts,
+			               vert_steps,
+			               num_verts,
+			               num_steps,
+			               0.0f,
+			               prev_verts);
+			BoundBox prev_bounds = BoundBox::empty;
+			prev_bounds.grow(prev_verts[0]);
+			prev_bounds.grow(prev_verts[1]);
+			prev_bounds.grow(prev_verts[2]);
+			/* Create all primitive time steps, */
+			for(int bvh_step = 1; bvh_step < num_bvh_steps; ++bvh_step) {
+				const float curr_time = (float)(bvh_step) * num_bvh_steps_inv_1;
+				float3 curr_verts[3];
+				t.motion_verts(verts,
+				               vert_steps,
+				               num_verts,
+				               num_steps,
+				               curr_time,
+				               curr_verts);
+				BoundBox curr_bounds = BoundBox::empty;
+				curr_bounds.grow(curr_verts[0]);
+				curr_bounds.grow(curr_verts[1]);
+				curr_bounds.grow(curr_verts[2]);
+				BoundBox bounds = prev_bounds;
+				bounds.grow(curr_bounds);
+				if(bounds.valid()) {
+					const float prev_time = (float)(bvh_step - 1) * num_bvh_steps_inv_1;
+					references.push_back(
+					        BVHReference(bounds,
+					                     j,
+					                     i,
+					                     PRIMITIVE_MOTION_TRIANGLE,
+					                     prev_time,
+					                     curr_time));
+					root.grow(bounds);
+					center.grow(bounds.center2());
+				}
+				/* Current time boundbox becomes previous one for the
+				 * next time step.
+				 */
+				prev_bounds = curr_bounds;
+			}
+		}
 	}
+}
 
-	if(params.primitive_mask & PRIMITIVE_ALL_CURVE) {
-		Attribute *curve_attr_mP = NULL;
-
-		if(mesh->has_motion_blur())
-			curve_attr_mP = mesh->curve_attributes.find(ATTR_STD_MOTION_VERTEX_POSITION);
-
-		size_t num_curves = mesh->num_curves();
-		for(uint j = 0; j < num_curves; j++) {
-			const Mesh::Curve curve = mesh->get_curve(j);
-			PrimitiveType type = PRIMITIVE_CURVE;
-			const float *curve_radius = &mesh->curve_radius[0];
-
-			for(int k = 0; k < curve.num_keys - 1; k++) {
+void BVHBuild::add_reference_curves(BoundBox& root, BoundBox& center, Mesh *mesh, int i)
+{
+	const Attribute *curve_attr_mP = NULL;
+	if(mesh->has_motion_blur()) {
+		curve_attr_mP = mesh->curve_attributes.find(ATTR_STD_MOTION_VERTEX_POSITION);
+	}
+	const size_t num_curves = mesh->num_curves();
+	for(uint j = 0; j < num_curves; j++) {
+		const Mesh::Curve curve = mesh->get_curve(j);
+		const float *curve_radius = &mesh->curve_radius[0];
+		for(int k = 0; k < curve.num_keys - 1; k++) {
+			if(curve_attr_mP == NULL) {
+				/* Really simple logic for static hair. */
 				BoundBox bounds = BoundBox::empty;
 				curve.bounds_grow(k, &mesh->curve_keys[0], curve_radius, bounds);
-
-				/* motion curve */
-				if(curve_attr_mP) {
-					const size_t mesh_size = mesh->curve_keys.size();
-					const size_t num_steps = mesh->motion_steps - 1;
-					const float3 *key_steps = curve_attr_mP->data_float3();
-
-					for(size_t step = 0; step < num_steps; step++) {
-						curve.bounds_grow(k, key_steps + step*mesh_size, curve_radius, bounds);
-					}
-
-					type = PRIMITIVE_MOTION_CURVE;
-				}
-
 				if(bounds.valid()) {
-					int packed_type = PRIMITIVE_PACK_SEGMENT(type, k);
-
+					int packed_type = PRIMITIVE_PACK_SEGMENT(PRIMITIVE_CURVE, k);
 					references.push_back(BVHReference(bounds, j, i, packed_type));
 					root.grow(bounds);
 					center.grow(bounds.center2());
 				}
 			}
+			else if(params.num_motion_curve_steps == 0 || params.use_spatial_split) {
+				/* Simple case of motion curves: single node for the while
+				 * shutter time. Lowest memory usage but less optimal
+				 * rendering.
+				 */
+				/* TODO(sergey): Support motion steps for spatially split BVH. */
+				BoundBox bounds = BoundBox::empty;
+				curve.bounds_grow(k, &mesh->curve_keys[0], curve_radius, bounds);
+				const size_t num_keys = mesh->curve_keys.size();
+				const size_t num_steps = mesh->motion_steps;
+				const float3 *key_steps = curve_attr_mP->data_float3();
+				for(size_t step = 0; step < num_steps - 1; step++) {
+					curve.bounds_grow(k,
+					                  key_steps + step*num_keys,
+					                  curve_radius,
+					                  bounds);
+				}
+				if(bounds.valid()) {
+					int packed_type = PRIMITIVE_PACK_SEGMENT(PRIMITIVE_MOTION_CURVE, k);
+					references.push_back(BVHReference(bounds,
+					                                  j,
+					                                  i,
+					                                  packed_type));
+					root.grow(bounds);
+					center.grow(bounds.center2());
+				}
+			}
+			else {
+				/* Motion curves, trace optimized case:  we split curve keys
+				 * primitives into separate nodes for each of the time steps.
+				 * This way we minimize overlap of neighbor curve primitives.
+				 */
+				const int num_bvh_steps = params.num_motion_curve_steps * 2 + 1;
+				const float num_bvh_steps_inv_1 = 1.0f / (num_bvh_steps - 1);
+				const size_t num_steps = mesh->motion_steps;
+				const float3 *curve_keys = &mesh->curve_keys[0];
+				const float3 *key_steps = curve_attr_mP->data_float3();
+				const size_t num_keys = mesh->curve_keys.size();
+				/* Calculate bounding box of the previous time step.
+				 * Will be reused later to avoid duplicated work on
+				 * calculating BVH time step boundbox.
+				 */
+				float4 prev_keys[4];
+				curve.cardinal_motion_keys(curve_keys,
+				                           curve_radius,
+				                           key_steps,
+				                           num_keys,
+				                           num_steps,
+				                           0.0f,
+				                           k - 1, k, k + 1, k + 2,
+				                           prev_keys);
+				BoundBox prev_bounds = BoundBox::empty;
+				curve.bounds_grow(prev_keys, prev_bounds);
+				/* Create all primitive time steps, */
+				for(int bvh_step = 1; bvh_step < num_bvh_steps; ++bvh_step) {
+					const float curr_time = (float)(bvh_step) * num_bvh_steps_inv_1;
+					float4 curr_keys[4];
+					curve.cardinal_motion_keys(curve_keys,
+					                           curve_radius,
+					                           key_steps,
+					                           num_keys,
+					                           num_steps,
+					                           curr_time,
+					                           k - 1, k, k + 1, k + 2,
+					                           curr_keys);
+					BoundBox curr_bounds = BoundBox::empty;
+					curve.bounds_grow(curr_keys, curr_bounds);
+					BoundBox bounds = prev_bounds;
+					bounds.grow(curr_bounds);
+					if(bounds.valid()) {
+						const float prev_time = (float)(bvh_step - 1) * num_bvh_steps_inv_1;
+						int packed_type = PRIMITIVE_PACK_SEGMENT(PRIMITIVE_MOTION_CURVE, k);
+						references.push_back(BVHReference(bounds,
+						                                  j,
+						                                  i,
+						                                  packed_type,
+						                                  prev_time,
+						                                  curr_time));
+						root.grow(bounds);
+						center.grow(bounds.center2());
+					}
+					/* Current time boundbox becomes previous one for the
+					 * next time step.
+					 */
+					prev_bounds = curr_bounds;
+				}
+			}
 		}
 	}
 }
 
+void BVHBuild::add_reference_mesh(BoundBox& root, BoundBox& center, Mesh *mesh, int i)
+{
+	if(params.primitive_mask & PRIMITIVE_ALL_TRIANGLE) {
+		add_reference_triangles(root, center, mesh, i);
+	}
+	if(params.primitive_mask & PRIMITIVE_ALL_CURVE) {
+		add_reference_curves(root, center, mesh, i);
+	}
+}
+
 void BVHBuild::add_reference_object(BoundBox& root, BoundBox& center, Object *ob, int i)
 {
 	references.push_back(BVHReference(ob->bounds, -1, i, 0));
@@ -203,7 +356,7 @@ static size_t count_curve_segments(Mesh *mesh)
 
 	for(size_t i = 0; i < num_curves; i++)
 		num += mesh->get_curve(i).num_keys - 1;
-	
+
 	return num;
 }
 
@@ -347,6 +500,7 @@ BVHNode* BVHBuild::run()
 		else {
 			/*rotate(rootnode, 4, 5);*/
 			rootnode->update_visibility();
+			rootnode->update_time();
 		}
 		if(rootnode != NULL) {
 			VLOG(1) << "BVH build statistics:\n"
@@ -374,7 +528,7 @@ void BVHBuild::progress_update()
 {
 	if(time_dt() - progress_start_time < 0.25)
 		return;
-	
+
 	double progress_start = (double)progress_count/(double)progress_total;
 	double duplicates = (double)(progress_total - progress_original_total)/(double)progress_total;
 
@@ -382,7 +536,7 @@ void BVHBuild::progress_update()
 	                           progress_start * 100.0, duplicates * 100.0);
 
 	progress.set_substatus(msg);
-	progress_start_time = time_dt(); 
+	progress_start_time = time_dt();
 }
 
 void BVHBuild::thread_build_node(InnerNode *inner,
@@ -696,18 +850,24 @@ BVHNode *BVHBuild::create_object_leaf_nodes(const BVHReference *ref, int start,
 		prim_object[start] = ref->prim_object();
 
 		uint visibility = objects[ref->prim_object()]->visibility;
-		return new LeafNode(ref->bounds(), visibility, start, start+1);
+		BVHNode *leaf_node =  new LeafNode(ref->bounds(), visibility, start, start+1);
+		leaf_node->m_time_from = ref->time_from();
+		leaf_node->m_time_to = ref->time_to();
+		return leaf_node;
 	}
 	else {
 		int mid = num/2;
-		BVHNode *leaf0 = create_object_leaf_nodes(ref, start, mid); 
-		BVHNode *leaf1 = create_object_leaf_nodes(ref+mid, start+mid, num-mid); 
+		BVHNode *leaf0 = create_object_leaf_nodes(ref, start, mid);
+		BVHNode *leaf1 = create_object_leaf_nodes(ref+mid, start+mid, num-mid);
 
 		BoundBox bounds = BoundBox::empty;
 		bounds.grow(leaf0->m_bounds);
 		bounds.grow(leaf1->m_bounds);
 
-		return new InnerNode(bounds, leaf0, leaf1);
+		BVHNode *inner_node = new InnerNode(bounds, leaf0, leaf1);
+		inner_node->m_time_from = min(leaf0->m_time_from, leaf1->m_time_from);
+		inner_node->m_time_to = max(leaf0->m_time_to, leaf1->m_time_to);
+		return inner_node;
 	}
 }
 
@@ -811,6 +971,16 @@ BVHNode* BVHBuild::create_leaf_node(const BVHRange& range,
 			                                   visibility[i],
 			                                   start_index,
 			                                   start_index + num);
+			if(true) {
+				float time_from = 1.0f, time_to = 0.0f;
+				for(int j = 0; j < num; ++j) {
+					const BVHReference &ref = p_ref[i][j];
+					time_from = min(time_from, ref.time_from());
+					time_to = max(time_to, ref.time_to());
+				}
+				leaf_node->m_time_from = time_from;
+				leaf_node->m_time_to = time_to;
+			}
 			if(alignment_found) {
 				/* Need to recalculate leaf bounds with new alignment. */
 				leaf_node->m_bounds = BoundBox::empty;
@@ -958,7 +1128,7 @@ void BVHBuild::rotate(BVHNode *node, int max_depth)
 	/* nothing to rotate if we reached a leaf node. */
 	if(node->is_leaf() || max_depth < 0)
 		return;
-	
+
 	InnerNode *parent = (InnerNode*)node;
 
 	/* rotate all children first */
diff --git a/intern/cycles/bvh/bvh_build.h b/intern/cycles/bvh/bvh_build.h
index 64180349935..ee3cde66a2f 100644
--- a/intern/cycles/bvh/bvh_build.h
+++ b/intern/cycles/bvh/bvh_build.h
@@ -63,6 +63,8 @@ protected:
 	friend class BVHObjectBinning;
 
 	/* Adding references. */
+	void add_reference_triangles(BoundBox& root, BoundBox& center, Mesh *mesh, int i);
+	void add_reference_curves(BoundBox& root, BoundBox& center, Mesh *mesh, int i);
 	void add_reference_mesh(BoundBox& root, BoundBox& center, Mesh *mesh, int i);
 	void add_reference_object(BoundBox& root, BoundBox& center, Object *ob, int i);
 	void add_references(BVHRange& root);
diff --git a/intern/cycles/bvh/bvh_node.cpp b/intern/cycles/bvh/bvh_node.cpp
index f5cd699bdf4..67580e1bc7b 100644
--- a/intern/cycles/bvh/bvh_node.cpp
+++ b/intern/cycles/bvh/bvh_node.cpp
@@ -176,6 +176,19 @@ uint BVHNode::update_visibility()
 	return m_visibility;
 }
 
+void BVHNode::update_time()
+{
+	if(!is_leaf()) {
+		InnerNode *inner = (InnerNode*)this;
+		BVHNode *child0 = inner->children[0];
+		BVHNode *child1 = inner->children[1];
+		child0->update_time();
+		child1->update_time();
+		m_time_from = min(child0->m_time_from, child1->m_time_from);
+		m_time_to =  max(child0->m_time_to, child1->m_time_to);
+	}
+}
+
 /* Inner Node */
 
 void InnerNode::print(int depth) const
diff --git a/intern/cycles/bvh/bvh_node.h b/intern/cycles/bvh/bvh_node.h
index 2faa40ab657..090c426de56 100644
--- a/intern/cycles/bvh/bvh_node.h
+++ b/intern/cycles/bvh/bvh_node.h
@@ -47,7 +47,9 @@ class BVHNode
 {
 public:
 	BVHNode() : m_is_unaligned(false),
-	            m_aligned_space(NULL)
+	            m_aligned_space(NULL),
+	            m_time_from(0.0f),
+	            m_time_to(1.0f)
 	{
 	}
 
@@ -91,12 +93,15 @@ public:
 	void deleteSubtree();
 
 	uint update_visibility();
+	void update_time();
 
 	bool m_is_unaligned;
 
 	// TODO(sergey): Can be stored as 3x3 matrix, but better to have some
 	// utilities and type defines in util_transform first.
 	Transform *m_aligned_space;
+
+	float m_time_from, m_time_to;
 };
 
 class InnerNode : public BVHNode
diff --git a/intern/cycles/bvh/bvh_params.h b/intern/cycles/bvh/bvh_params.h
index 6d426475737..65f9da1c194 100644
--- a/intern/cycles/bvh/bvh_params.h
+++ b/intern/cycles/bvh/bvh_params.h
@@ -61,6 +61,17 @@ public:
 	 */
 	bool use_unaligned_nodes;
 
+	/* Split time range to this number of steps and create leaf node for each
+	 * of this time steps.
+	 *
+	 * Speeds up rendering of motion curve primitives in the cost of higher
+	 * memory usage.
+	 */
+	int num_motion_curve_steps;
+
+	/* Same as above, but for triangle primitives. */
+	int num_motion_triangle_steps;
+
 	/* fixed parameters */
 	enum {
 		MAX_DEPTH = 64,
@@ -91,6 +102,8 @@ public:
 		use_unaligned_nodes = false;
 
 		primitive_mask = PRIMITIVE_ALL;
+
+		num_motion_curve_steps = 0;
 	}
 
 	/* SAH costs */
@@ -117,8 +130,15 @@ class BVHReference
 public:
 	__forceinline BVHReference() {}
 
-	__forceinline BVHReference(const BoundBox& bounds_, int prim_index_, int prim_object_, int prim_type)
-	: rbounds(bounds_)
+	__forceinline BVHReference(const BoundBox& bounds_,
+	                           int prim_index_,
+	                           int prim_object_,
+	                           int prim_type,
+	                           float time_from = 0.0f,
+	                           float time_to = 1.0f)
+	        : rbounds(bounds_),
+	          time_from_(time_from),
+	          time_to_(time_to)
 	{
 		rbounds.min.w = __int_as_float(prim_index_);
 		rbounds.max.w = __int_as_float(prim_object_);
@@ -129,6 +149,9 @@ public:
 	__forceinline int prim_index() const { return __float_as_int(rbounds.min.w); }
 	__forceinline int prim_object() const { return __float_as_int(rbounds.max.w); }
 	__forceinline int prim_type() const { return type; }
+	__forceinline float time_from() const { return time_from_; }
+	__forceinline float time_to() const { return time_to_; }
+
 
 	BVHReference& operator=(const BVHReference &arg) {
 		if(&arg != this) {
@@ -137,9 +160,11 @@ public:
 		return *this;
 	}
 
+
 protected:
 	BoundBox rbounds;
 	uint type;
+	float time_from_, time_to_;
 };
 
 /* BVH Range
diff --git a/intern/cycles/kernel/CMakeLists.txt b/intern/cycles/kernel/CMakeLists.txt
index 56bcafbce38..29e0f44841e 100644
--- a/intern/cycles/kernel/CMakeLists.txt
+++ b/intern/cycles/kernel/CMakeLists.txt
@@ -164,6 +164,8 @@ set(SRC_GEOM_HEADERS
 	geom/geom_curve.h
 	geom/geom_motion_curve.h
 	geom/geom_motion_triangle.h
+	geom/geom_motion_triangle_intersect.h
+	geom/geom_motion_triangle_shader.h
 	geom/geom_object.h
 	geom/geom_patch.h
 	geom/geom_primitive.h
diff --git a/intern/cycles/kernel/bvh/qbvh_shadow_all.h b/intern/cycles/kernel/bvh/qbvh_shadow_all.h
index b2e99725626..607295f9ed5 100644
--- a/intern/cycles/kernel/bvh/qbvh_shadow_all.h
+++ b/intern/cycles/kernel/bvh/qbvh_shadow_all.h
@@ -106,14 +106,20 @@ ccl_device bool BVH_FUNCTION_FULL_NAME(QBVH)(KernelGlobals *kg,
 			while(node_addr >= 0 && node_addr != ENTRYPOINT_SENTINEL) {
 				float4 inodes = kernel_tex_fetch(__bvh_nodes, node_addr+0);
 
+				if(false
 #ifdef __VISIBILITY_FLAG__
-				if((__float_as_uint(inodes.x) & PATH_RAY_SHADOW) == 0) {
+				   || ((__float_as_uint(inodes.x) & PATH_RAY_SHADOW) == 0)
+#endif
+#if BVH_FEATURE(BVH_MOTION)
+				   || UNLIKELY(ray->time < inodes.y)
+				   || UNLIKELY(ray->time > inodes.z)
+#endif
+				) {
 					/* Pop. */
 					node_addr = traversal_stack[stack_ptr].addr;
 					--stack_ptr;
 					continue;
 				}
-#endif
 
 				ssef dist;
 				int child_mask = NODE_INTERSECT(kg,
diff --git a/intern/cycles/kernel/bvh/qbvh_traversal.h b/intern/cycles/kernel/bvh/qbvh_traversal.h
index 1d5643ca540..10ae7bee852 100644
--- a/intern/cycles/kernel/bvh/qbvh_traversal.h
+++ b/intern/cycles/kernel/bvh/qbvh_traversal.h
@@ -117,6 +117,10 @@ ccl_device bool BVH_FUNCTION_FULL_NAME(QBVH)(KernelGlobals *kg,
 				float4 inodes = kernel_tex_fetch(__bvh_nodes, node_addr+0);
 
 				if(UNLIKELY(node_dist > isect->t)
+#if BVH_FEATURE(BVH_MOTION)
+				   || UNLIKELY(ray->time < inodes.y)
+				   || UNLIKELY(ray->time > inodes.z)
+#endif
 #ifdef __VISIBILITY_FLAG__
 				   || (__float_as_uint(inodes.x) & visibility) == 0)
 #endif
diff --git a/intern/cycles/kernel/geom/geom.h b/intern/cycles/kernel/geom/geom.h
index 24ced934c8b..6838e26c242 100644
--- a/intern/cycles/kernel/geom/geom.h
+++ b/intern/cycles/kernel/geom/geom.h
@@ -23,6 +23,8 @@
 #include "geom_subd_triangle.h"
 #include "geom_triangle_intersect.h"
 #include "geom_motion_triangle.h"
+#include "geom_motion_triangle_intersect.h"
+#include "geom_motion_triangle_shader.h"
 #include "geom_motion_curve.h"
 #include "geom_curve.h"
 #include "geom_volume.h"
diff --git a/intern/cycles/kernel/geom/geom_motion_curve.h b/intern/cycles/kernel/geom/geom_motion_curve.h
index 80b33fad68b..dc1388b6643 100644
--- a/intern/cycles/kernel/geom/geom_motion_curve.h
+++ b/intern/cycles/kernel/geom/geom_motion_curve.h
@@ -50,12 +50,12 @@ ccl_device_inline int find_attribute_curve_motion(KernelGlobals *kg, int object,
 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])
 {
 	if(step == numsteps) {
-		/* center step: regular vertex location */
+		/* center step: regular key location */
 		keys[0] = kernel_tex_fetch(__curve_keys, k0);
 		keys[1] = kernel_tex_fetch(__curve_keys, k1);
 	}
 	else {
-		/* center step not stored in this array */
+		/* center step is not stored in this array */
 		if(step > numsteps)
 			step--;
 
@@ -97,14 +97,14 @@ ccl_device_inline void motion_curve_keys(KernelGlobals *kg, int object, int prim
 ccl_device_inline void motion_cardinal_curve_keys_for_step(KernelGlobals *kg, int offset, int numkeys, int numsteps, int step, int k0, int k1, int k2, int k3, float4 keys[4])
 {
 	if(step == numsteps) {
-		/* center step: regular vertex location */
+		/* center step: regular key location */
 		keys[0] = kernel_tex_fetch(__curve_keys, k0);
 		keys[1] = kernel_tex_fetch(__curve_keys, k1);
 		keys[2] = kernel_tex_fetch(__curve_keys, k2);
 		keys[3] = kernel_tex_fetch(__curve_keys, k3);
 	}
 	else {
-		/* center step not store in this array */
+		/* center step is not stored in this array */
 		if(step > numsteps)
 			step--;
 
diff --git a/intern/cycles/kernel/geom/geom_motion_triangle.h b/intern/cycles/kernel/geom/geom_motion_triangle.h
index 538c332c63a..4e84aa97776 100644
--- a/intern/cycles/kernel/geom/geom_motion_triangle.h
+++ b/intern/cycles/kernel/geom/geom_motion_triangle.h
@@ -76,7 +76,7 @@ ccl_device_inline void motion_triangle_normals_for_step(KernelGlobals *kg, uint4
 		normals[2] = float4_to_float3(kernel_tex_fetch(__tri_vnormal, tri_vindex.z));
 	}
 	else {
-		/* center step not stored in this array */
+		/* center step is not stored in this array */
 		if(step > numsteps)
 			step--;
 
@@ -117,312 +117,4 @@ ccl_device_inline void motion_triangle_vertices(KernelGlobals *kg, int object, i
 	verts[2] = (1.0f - t)*verts[2] + t*next_verts[2];
 }
 
-/* Refine triangle intersection to more precise hit point. For rays that travel
- * far the precision is often not so good, this reintersects the primitive from
- * a closer distance. */
-
-ccl_device_inline float3 motion_triangle_refine(KernelGlobals *kg, ShaderData *sd, const Intersection *isect, const Ray *ray, float3 verts[3])
-{
-	float3 P = ray->P;
-	float3 D = ray->D;
-	float t = isect->t;
-
-#ifdef __INTERSECTION_REFINE__
-	if(isect->object != OBJECT_NONE) {
-		if(UNLIKELY(t == 0.0f)) {
-			return P;
-		}
-#  ifdef __OBJECT_MOTION__
-		Transform tfm = ccl_fetch(sd, ob_itfm);
-#  else
-		Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
-#  endif
-
-		P = transform_point(&tfm, P);
-		D = transform_direction(&tfm, D*t);
-		D = normalize_len(D, &t);
-	}
-
-	P = P + D*t;
-
-	/* compute refined intersection distance */
-	const float3 e1 = verts[0] - verts[2];
-	const float3 e2 = verts[1] - verts[2];
-	const float3 s1 = cross(D, e2);
-
-	const float invdivisor = 1.0f/dot(s1, e1);
-	const float3 d = P - verts[2];
-	const float3 s2 = cross(d, e1);
-	float rt = dot(e2, s2)*invdivisor;
-
-	/* compute refined position */
-	P = P + D*rt;
-
-	if(isect->object != OBJECT_NONE) {
-#  ifdef __OBJECT_MOTION__
-		Transform tfm = ccl_fetch(sd, ob_tfm);
-#  else
-		Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
-#  endif
-
-		P = transform_point(&tfm, P);
-	}
-
-	return P;
-#else
-	return P + D*t;
-#endif
-}
-
-/* Same as above, except that isect->t is assumed to be in object space for instancing */
-
-#ifdef __SUBSURFACE__
-#  if defined(__KERNEL_CUDA__) && (defined(i386) || defined(_M_IX86))
-ccl_device_noinline
-#  else
-ccl_device_inline
-#  endif
-float3 motion_triangle_refine_subsurface(KernelGlobals *kg, ShaderData *sd, const Intersection *isect, const Ray *ray, float3 verts[3])
-{
-	float3 P = ray->P;
-	float3 D = ray->D;
-	float t = isect->t;
-
-#  ifdef __INTERSECTION_REFINE__
-	if(isect->object != OBJECT_NONE) {
-#    ifdef __OBJECT_MOTION__
-		Transform tfm = ccl_fetch(sd, ob_itfm);
-#    else
-		Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_INVERSE_TRANSFORM);
-#    endif
-
-		P = transform_point(&tfm, P);
-		D = transform_direction(&tfm, D);
-		D = normalize(D);
-	}
-
-	P = P + D*t;
-
-	/* compute refined intersection distance */
-	const float3 e1 = verts[0] - verts[2];
-	const float3 e2 = verts[1] - verts[2];
-	const float3 s1 = cross(D, e2);
-
-	const float invdivisor = 1.0f/dot(s1, e1);
-	const float3 d = P - verts[2];
-	const float3 s2 = cross(d, e1);
-	float rt = dot(e2, s2)*invdivisor;
-
-	P = P + D*rt;
-
-	if(isect->object != OBJECT_NONE) {
-#    ifdef __OBJECT_MOTION__
-		Transform tfm = ccl_fetch(sd, ob_tfm);
-#    else
-		Transform tfm = object_fetch_transform(kg, isect->object, OBJECT_TRANSFORM);
-#    endif
-
-		P = transform_point(&tfm, P);
-	}
-
-	return P;
-#  else
-	return P + D*t;
-#  endif
-}
-#endif
-
-/* Setup of motion triangle specific parts of ShaderData, moved into this one
- * function to more easily share computation of interpolated positions and
- * normals */
-
-/* return 3 triangle vertex normals */
-ccl_device_noinline void motion_triangle_shader_setup(KernelGlobals *kg, ShaderData *sd, const Intersection *isect, const Ray *ray, bool subsurface)
-{
-	/* get shader */
-	ccl_fetch(sd, shader) = kernel_tex_fetch(__tri_shader, ccl_fetch(sd, prim));
-
-	/* get motion info */
-	int numsteps, numverts;
-	object_motion_info(kg, ccl_fetch(sd, object), &numsteps, &numverts, NULL);
-
-	/* figure out which steps we need to fetch and their interpolation factor */
-	int maxstep = numsteps*2;
-	int step = min((int)(ccl_fetch(sd, time)*maxstep), maxstep-1);
-	float t = ccl_fetch(sd, time)*maxstep - step;
-
-	/* find attribute */
-	AttributeElement elem;
-	int offset = find_attribute_motion(kg, ccl_fetch(sd, object), ATTR_STD_MOTION_VERTEX_POSITION, &elem);
-	kernel_assert(offset != ATTR_STD_NOT_FOUND);
-
-	/* fetch vertex coordinates */
-	float3 verts[3], next_verts[3];
-	uint4 tri_vindex = kernel_tex_fetch(__tri_vindex, ccl_fetch(sd, prim));
-
-	motion_triangle_verts_for_step(kg, tri_vindex, offset, numverts, numsteps, step, verts);
-	motion_triangle_verts_for_step(kg, tri_vindex, offset, numverts, numsteps, step+1, next_verts);
-
-	/* interpolate between steps */
-	verts[0] = (1.0f - t)*verts[0] + t*next_verts[0];
-	verts[1] = (1.0f - t)*verts[1] + t*next_verts[1];
-	verts[2] = (1.0f - t)*verts[2] + t*next_verts[2];
-
-	/* compute refined position */
-#ifdef __SUBSURFACE__
-	if(!subsurface)
-#endif
-		ccl_fetch(sd, P) = motion_triangle_refine(kg, sd, isect, ray, verts);
-#ifdef __SUBSURFACE__
-	else
-		ccl_fetch(sd, P) = motion_triangle_refine_subsurface(kg, sd, isect, ray, verts);
-#endif
-
-	/* compute face normal */
-	float3 Ng;
-	if(ccl_fetch(sd, flag) & SD_NEGATIVE_SCALE_APPLIED)
-		Ng = normalize(cross(verts[2] - verts[0], verts[1] - verts[0]));
-	else
-		Ng = normalize(cross(verts[1] - verts[0], verts[2] - verts[0]));
-
-	ccl_fetch(sd, Ng) = Ng;
-	ccl_fetch(sd, N) = Ng;
-
-	/* compute derivatives of P w.r.t. uv */
-#ifdef __DPDU__
-	ccl_fetch(sd, dPdu) = (verts[0] - verts[2]);
-	ccl_fetch(sd, dPdv) = (verts[1] - verts[2]);
-#endif
-
-	/* compute smooth normal */
-	if(ccl_fetch(sd, shader) & SHADER_SMOOTH_NORMAL) {
-		/* find attribute */
-		AttributeElement elem;
-		int offset = find_attribute_motion(kg, ccl_fetch(sd, object), ATTR_STD_MOTION_VERTEX_NORMAL, &elem);
-		kernel_assert(offset != ATTR_STD_NOT_FOUND);
-
-		/* fetch vertex coordinates */
-		float3 normals[3], next_normals[3];
-		motion_triangle_normals_for_step(kg, tri_vindex, offset, numverts, numsteps, step, normals);
-		motion_triangle_normals_for_step(kg, tri_vindex, offset, numverts, numsteps, step+1, next_normals);
-
-		/* interpolate between steps */
-		normals[0] = (1.0f - t)*normals[0] + t*next_normals[0];
-		normals[1] = (1.0f - t)*normals[1] + t*next_normals[1];
-		normals[2] = (1.0f - t)*normals[2] + t*next_normals[2];
-
-		/* interpolate between vertices */
-		float u = ccl_fetch(sd, u);
-		float v = ccl_fetch(sd, v);
-		float w = 1.0f - u - v;
-		ccl_fetch(sd, N) = (u*normals[0] + v*normals[1] + w*normals[2]);
-	}
-}
-
-/* Ray intersection. We simply compute the vertex positions at the given ray
- * time and do a ray intersection with the resulting triangle */
-
-ccl_device_inline bool motion_triangle_intersect(KernelGlobals *kg, Intersection *isect,
-	float3 P, float3 dir, float time, uint visibility, int object, int prim_addr)
-{
-	/* primitive index for vertex location lookup */
-	int prim = kernel_tex_fetch(__prim_index, prim_addr);
-	int fobject = (object == OBJECT_NONE)? kernel_tex_fetch(__prim_object, prim_addr): object;
-
-	/* get vertex locations for intersection */
-	float3 verts[3];
-	motion_triangle_vertices(kg, fobject, prim, time, verts);
-
-	/* ray-triangle intersection, unoptimized */
-	float t, u, v;
-
-	if(ray_triangle_intersect_uv(P, dir, isect->t, verts[2], verts[0], verts[1], &u, &v, &t)) {
-#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, prim_addr) & visibility)
-#endif
-		{
-			isect->t = t;
-			isect->u = u;
-			isect->v = v;
-			isect->prim = prim_addr;
-			isect->object = object;
-			isect->type = PRIMITIVE_MOTION_TRIANGLE;
-		
-			return true;
-		}
-	}
-
-	return false;
-}
-
-/* Special ray intersection routines for subsurface scattering. In that case we
- * only want to intersect with primitives in the same object, and if case of
- * multiple hits we pick a single random primitive as the intersection point. */
-
-#ifdef __SUBSURFACE__
-ccl_device_inline void motion_triangle_intersect_subsurface(
-        KernelGlobals *kg,
-        SubsurfaceIntersection *ss_isect,
-        float3 P,
-        float3 dir,
-        float time,
-        int object,
-        int prim_addr,
-        float tmax,
-        uint *lcg_state,
-        int max_hits)
-{
-	/* primitive index for vertex location lookup */
-	int prim = kernel_tex_fetch(__prim_index, prim_addr);
-	int fobject = (object == OBJECT_NONE)? kernel_tex_fetch(__prim_object, prim_addr): object;
-
-	/* get vertex locations for intersection */
-	float3 verts[3];
-	motion_triangle_vertices(kg, fobject, prim, time, verts);
-
-	/* ray-triangle intersection, unoptimized */
-	float t, u, v;
-
-	if(ray_triangle_intersect_uv(P, dir, tmax, verts[2], verts[0], verts[1], &u, &v, &t)) {
-		for(int i = min(max_hits, ss_isect->num_hits) - 1; i >= 0; --i) {
-			if(ss_isect->hits[i].t == t) {
-				return;
-			}
-		}
-
-		ss_isect->num_hits++;
-
-		int hit;
-
-		if(ss_isect->num_hits <= max_hits) {
-			hit = ss_isect->num_hits - 1;
-		}
-		else {
-			/* reservoir sampling: if we are at the maximum number of
-			 * hits, randomly replace element or skip it */
-			hit = lcg_step_uint(lcg_state) % ss_isect->num_hits;
-
-			if(hit >= max_hits)
-				return;
-		}
-
-		/* record intersection */
-		Intersection *isect = &ss_isect->hits[hit];
-		isect->t = t;
-		isect->u = u;
-		isect->v = v;
-		isect->prim = prim_addr;
-		isect->object = object;
-		isect->type = PRIMITIVE_MOTION_TRIANGLE;
-
-		/* Record geometric normal. */
-		ss_isect->Ng[hit] = normalize(cross(verts[1] - verts[0],
-		                                    verts[2] - verts[0]));
-	}
-}
-#endif
-
 CCL_NAMESPACE_END
-
diff --git a/intern/cycles/kernel/geom/geom_motion_triangle_intersect.h b/intern/cycles/kernel/geom/geom_motion_triangle_intersect.h
new file mode 100644
index 00000000000..d57d74ea882
--- /dev/null
+++ b/intern/cycles/kernel/geom/geom_motion_triangle_intersect.h
@@ -0,0 +1,280 @@
+/*
+ * Copyright 2011-2016 Blender Foundation
+ *
+ * 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
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/* Motion Triangle Primitive
+ *
+ * These are stored as regular triangles, plus extra positions and normals at
+ * times other than the frame center. Computing the triangle vertex positions
+ * or normals at a given ray time is a matter of interpolation of the two steps
+ * between which the ray time lies.
+ *
+ * The extra positions and normals are stored as ATTR_STD_MOTION_VERTEX_POSITION
+ * and ATTR_STD_MOTION_VERTEX_NORMAL mesh attributes.
+ */
+
+CCL_NAMESPACE_BEGIN
+
+/* Refine triangle intersection to more precise hit point. For rays that travel
+ * far the precision is often not so good, this reintersects the primitive from
+ * a closer distance.
+ */
+
+ccl_device_inline float3 motion_triangle_refine(KernelGlobals *kg,
+                                                ShaderData *sd,
+                                                const Intersection *isect,
+                                                const Ray *ray,
+                                                float3 verts[3])
+{
+	float3 P = ray->P;
+	float3 D = ray->D;
+	float t = isect->t;
+
+#ifdef __INTERSECTION_REFINE__
+	if(isect->object != OBJECT_NONE) {
+		if(UNLIKELY(t == 0.0f)) {
+			return P;
+		}
+#  ifdef __OBJECT_MOTION__
+		Transform tfm = ccl_fetch(sd, ob_itfm);
+#  else
+		Transform tfm = object_fetch_transform(kg,
+		                                       isect->object,
+		                                       OBJECT_INVERSE_TRANSFORM);
+#  endif
+
+		P = transform_point(&tfm, P);
+		D = transform_direction(&tfm, D*t);
+		D = normalize_len(D, &t);
+	}
+
+	P = P + D*t;
+
+	/* Compute refined intersection distance. */
+	const float3 e1 = verts[0] - verts[2];
+	const float3 e2 = verts[1] - verts[2];
+	const float3 s1 = cross(D, e2);
+
+	const float invdivisor = 1.0f/dot(s1, e1);
+	const float3 d = P - verts[2];
+	const float3 s2 = cross(d, e1);
+	float rt = dot(e2, s2)*invdivisor;
+
+	/* Compute refined position. */
+	P = P + D*rt;
+
+	if(isect->object != OBJECT_NONE) {
+#  ifdef __OBJECT_MOTION__
+		Transform tfm = ccl_fetch(sd, ob_tfm);
+#  else
+		Transform tfm = object_fetch_transform(kg,
+		                                       isect->object,
+		                                       OBJECT_TRANSFORM);
+#  endif
+
+		P = transform_point(&tfm, P);
+	}
+
+	return P;
+#else
+	return P + D*t;
+#endif
+}
+
+/* Same as above, except that isect->t is assumed to be in object space
+ * for instancing.
+ */
+
+#ifdef __SUBSURFACE__
+#  if defined(__KERNEL_CUDA__) && (defined(i386) || defined(_M_IX86))
+ccl_device_noinline
+#  else
+ccl_device_inline
+#  endif
+float3 motion_triangle_refine_subsurface(KernelGlobals *kg,
+                                         ShaderData *sd,
+                                         const Intersection *isect,
+                                         const Ray *ray,
+                                         float3 verts[3])
+{
+	float3 P = ray->P;
+	float3 D = ray->D;
+	float t = isect->t;
+
+#  ifdef __INTERSECTION_REFINE__
+	if(isect->object != OBJECT_NONE) {
+#    ifdef __OBJECT_MOTION__
+		Transform tfm = ccl_fetch(sd, ob_itfm);
+#    else
+		Transform tfm = object_fetch_transform(kg,
+		                                       isect->object,
+		                                       OBJECT_INVERSE_TRANSFORM);
+#    endif
+
+		P = transform_point(&tfm, P);
+		D = transform_direction(&tfm, D);
+		D = normalize(D);
+	}
+
+	P = P + D*t;
+
+	/* compute refined intersection distance */
+	const float3 e1 = verts[0] - verts[2];
+	const float3 e2 = verts[1] - verts[2];
+	const float3 s1 = cross(D, e2);
+
+	const float invdivisor = 1.0f/dot(s1, e1);
+	const float3 d = P - verts[2];
+	const float3 s2 = cross(d, e1);
+	float rt = dot(e2, s2)*invdivisor;
+
+	P = P + D*rt;
+
+	if(isect->object != OBJECT_NONE) {
+#    ifdef __OBJECT_MOTION__
+		Transform tfm = ccl_fetch(sd, ob_tfm);
+#    else
+		Transform tfm = object_fetch_transform(kg,
+		                                       isect->object,
+		                                       OBJECT_TRANSFORM);
+#    endif
+
+		P = transform_point(&tfm, P);
+	}
+
+	return P;
+#  else  /* __INTERSECTION_REFINE__ */
+	return P + D*t;
+#  endif  /* __INTERSECTION_REFINE__ */
+}
+#endif  /* __SUBSURFACE__ */
+
+
+/* Ray intersection. We simply compute the vertex positions at the given ray
+ * time and do a ray intersection with the resulting triangle.
+ */
+
+ccl_device_inline bool motion_triangle_intersect(KernelGlobals *kg,
+                                                 Intersection *isect,
+                                                 float3 P,
+                                                 float3 dir,
+                                                 float time,
+                                                 uint visibility,
+                                                 int object,
+                                                 int prim_addr)
+{
+	/* Primitive index for vertex location lookup. */
+	int prim = kernel_tex_fetch(__prim_index, prim_addr);
+	int fobject = (object == OBJECT_NONE)
+	                  ? kernel_tex_fetch(__prim_object, prim_addr)
+	                  : object;
+	/* Get vertex locations for intersection. */
+	float3 verts[3];
+	motion_triangle_vertices(kg, fobject, prim, time, verts);
+	/* Ray-triangle intersection, unoptimized. */
+	float t, u, v;
+	if(ray_triangle_intersect_uv(P,
+	                             dir,
+	                             isect->t,
+	                             verts[2], verts[0], verts[1],
+	                             &u, &v, &t))
+	{
+#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, prim_addr) & visibility)
+#endif
+		{
+			isect->t = t;
+			isect->u = u;
+			isect->v = v;
+			isect->prim = prim_addr;
+			isect->object = object;
+			isect->type = PRIMITIVE_MOTION_TRIANGLE;
+			return true;
+		}
+	}
+	return false;
+}
+
+/* Special ray intersection routines for subsurface scattering. In that case we
+ * only want to intersect with primitives in the same object, and if case of
+ * multiple hits we pick a single random primitive as the intersection point.
+ */
+#ifdef __SUBSURFACE__
+ccl_device_inline void motion_triangle_intersect_subsurface(
+        KernelGlobals *kg,
+        SubsurfaceIntersection *ss_isect,
+        float3 P,
+        float3 dir,
+        float time,
+        int object,
+        int prim_addr,
+        float tmax,
+        uint *lcg_state,
+        int max_hits)
+{
+	/* Primitive index for vertex location lookup. */
+	int prim = kernel_tex_fetch(__prim_index, prim_addr);
+	int fobject = (object == OBJECT_NONE)
+	                  ? kernel_tex_fetch(__prim_object, prim_addr)
+	                  : object;
+	/* Get vertex locations for intersection. */
+	float3 verts[3];
+	motion_triangle_vertices(kg, fobject, prim, time, verts);
+	/* Ray-triangle intersection, unoptimized. */
+	float t, u, v;
+	if(ray_triangle_intersect_uv(P,
+	                             dir,
+	                             tmax,
+	                             verts[2], verts[0], verts[1],
+	                             &u, &v, &t))
+	{
+		for(int i = min(max_hits, ss_isect->num_hits) - 1; i >= 0; --i) {
+			if(ss_isect->hits[i].t == t) {
+				return;
+			}
+		}
+		ss_isect->num_hits++;
+		int hit;
+		if(ss_isect->num_hits <= max_hits) {
+			hit = ss_isect->num_hits - 1;
+		}
+		else {
+			/* Reservoir sampling: if we are at the maximum number of
+			 * hits, randomly replace element or skip it.
+			 */
+			hit = lcg_step_uint(lcg_state) % ss_isect->num_hits;
+
+			if(hit >= max_hits)
+				return;
+		}
+		/* Record intersection. */
+		Intersection *isect = &ss_isect->hits[hit];
+		isect->t = t;
+		isect->u = u;
+		isect->v = v;
+		isect->prim = prim_addr;
+		isect->object = object;
+		isect->type = PRIMITIVE_MOTION_TRIANGLE;
+		/* Record geometric normal. */
+		ss_isect->Ng[hit] = normalize(cross(verts[1] - verts[0],
+		                                    verts[2] - verts[0]));
+	}
+}
+#endif  /* __SUBSURFACE__ */
+
+CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/geom/geom_motion_triangle_shader.h b/intern/cycles/kernel/geom/geom_motion_triangle_shader.h
new file mode 100644
index 00000000000..c5dbc6a2f52
--- /dev/null
+++ b/intern/cycles/kernel/geom/geom_motion_triangle_shader.h
@@ -0,0 +1,123 @@
+/*
+ * Copyright 2011-2016 Blender Foundation
+ *
+ * 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
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/* Motion Triangle Primitive
+ *
+ * These are stored as regular triangles, plus extra positions and normals at
+ * times other than the frame center. Computing the triangle vertex positions
+ * or normals at a given ray time is a matter of interpolation of the two steps
+ * between which the ray time lies.
+ *
+ * The extra positions and normals are stored as ATTR_STD_MOTION_VERTEX_POSITION
+ * and ATTR_STD_MOTION_VERTEX_NORMAL mesh attributes.
+ */
+
+CCL_NAMESPACE_BEGIN
+
+/* Setup of motion triangle specific parts of ShaderData, moved into this one
+ * function to more easily share computation of interpolated positions and
+ * normals */
+
+/* return 3 triangle vertex normals */
+ccl_device_noinline void motion_triangle_shader_setup(KernelGlobals *kg,
+                                                      ShaderData *sd, const
+                                                      Intersection *isect,
+                                                      const Ray *ray,
+                                                      bool subsurface)
+{
+	/* Get shader. */
+	ccl_fetch(sd, shader) = kernel_tex_fetch(__tri_shader, ccl_fetch(sd, prim));
+	/* Get motion info. */
+	/* TODO(sergey): This logic is really similar to motion_triangle_vertices(),
+	 * can we de-duplicate something here?
+	 */
+	int numsteps, numverts;
+	object_motion_info(kg, ccl_fetch(sd, object), &numsteps, &numverts, NULL);
+	/* Figure out which steps we need to fetch and their interpolation factor. */
+	int maxstep = numsteps*2;
+	int step = min((int)(ccl_fetch(sd, time)*maxstep), maxstep-1);
+	float t = ccl_fetch(sd, time)*maxstep - step;
+	/* Find attribute. */
+	AttributeElement elem;
+	int offset = find_attribute_motion(kg, ccl_fetch(sd, object),
+	                                   ATTR_STD_MOTION_VERTEX_POSITION,
+	                                   &elem);
+	kernel_assert(offset != ATTR_STD_NOT_FOUND);
+	/* Fetch vertex coordinates. */
+	float3 verts[3], next_verts[3];
+	uint4 tri_vindex = kernel_tex_fetch(__tri_vindex, ccl_fetch(sd, prim));
+	motion_triangle_verts_for_step(kg, tri_vindex, offset, numverts, numsteps, step, verts);
+	motion_triangle_verts_for_step(kg, tri_vindex, offset, numverts, numsteps, step+1, next_verts);
+	/* Interpolate between steps. */
+	verts[0] = (1.0f - t)*verts[0] + t*next_verts[0];
+	verts[1] = (1.0f - t)*verts[1] + t*next_verts[1];
+	verts[2] = (1.0f - t)*verts[2] + t*next_verts[2];
+	/* Compute refined position. */
+#ifdef __SUBSURFACE__
+	if(subsurface) {
+		ccl_fetch(sd, P) = motion_triangle_refine_subsurface(kg,
+		                                                     sd,
+		                                                     isect,
+		                                                     ray,
+		                                                     verts);
+	}
+	else
+#endif  /*  __SUBSURFACE__*/
+	{
+		ccl_fetch(sd, P) = motion_triangle_refine(kg, sd, isect, ray, verts);
+	}
+	/* Compute face normal. */
+	float3 Ng;
+	if(ccl_fetch(sd, flag) & SD_NEGATIVE_SCALE_APPLIED) {
+		Ng = normalize(cross(verts[2] - verts[0], verts[1] - verts[0]));
+	}
+	else {
+		Ng = normalize(cross(verts[1] - verts[0], verts[2] - verts[0]));
+	}
+	ccl_fetch(sd, Ng) = Ng;
+	ccl_fetch(sd, N) = Ng;
+	/* Compute derivatives of P w.r.t. uv. */
+#ifdef __DPDU__
+	ccl_fetch(sd, dPdu) = (verts[0] - verts[2]);
+	ccl_fetch(sd, dPdv) = (verts[1] - verts[2]);
+#endif
+	/* Compute smooth normal. */
+	if(ccl_fetch(sd, shader) & SHADER_SMOOTH_NORMAL) {
+		/* Find attribute. */
+		AttributeElement elem;
+		int offset = find_attribute_motion(kg,
+		                                   ccl_fetch(sd, object),
+		                                   ATTR_STD_MOTION_VERTEX_NORMAL,
+		                                   &elem);
+		kernel_assert(offset != ATTR_STD_NOT_FOUND);
+		/* Fetch vertex coordinates. */
+		float3 normals[3], next_normals[3];
+		motion_triangle_normals_for_step(kg, tri_vindex, offset, numverts, numsteps, step, normals);
+		motion_triangle_normals_for_step(kg, tri_vindex, offset, numverts, numsteps, step+1, next_normals);
+		/* Interpolate between steps. */
+		normals[0] = (1.0f - t)*normals[0] + t*next_normals[0];
+		normals[1] = (1.0f - t)*normals[1] + t*next_normals[1];
+		normals[2] = (1.0f - t)*normals[2] + t*next_normals[2];
+		/* Interpolate between vertices. */
+		float u = ccl_fetch(sd, u);
+		float v = ccl_fetch(sd, v);
+		float w = 1.0f - u - v;
+		ccl_fetch(sd, N) = (u*normals[0] + v*normals[1] + w*normals[2]);
+	}
+}
+
+CCL_NAMESPACE_END
+
diff --git a/intern/cycles/kernel/osl/osl_services.cpp b/intern/cycles/kernel/osl/osl_services.cpp
index 26543862b80..eeccf9a9641 100644
--- a/intern/cycles/kernel/osl/osl_services.cpp
+++ b/intern/cycles/kernel/osl/osl_services.cpp
@@ -102,6 +102,8 @@ ustring OSLRenderServices::u_curve_tangent_normal("geom:curve_tangent_normal");
 #endif
 ustring OSLRenderServices::u_path_ray_length("path:ray_length");
 ustring OSLRenderServices::u_path_ray_depth("path:ray_depth");
+ustring OSLRenderServices::u_path_diffuse_depth("path:diffuse_depth");
+ustring OSLRenderServices::u_path_glossy_depth("path:glossy_depth");
 ustring OSLRenderServices::u_path_transparent_depth("path:transparent_depth");
 ustring OSLRenderServices::u_path_transmission_depth("path:transmission_depth");
 ustring OSLRenderServices::u_trace("trace");
@@ -759,6 +761,24 @@ bool OSLRenderServices::get_background_attribute(KernelGlobals *kg, ShaderData *
 		int f = state->bounce;
 		return set_attribute_int(f, type, derivatives, val);
 	}
+	else if(name == u_path_diffuse_depth) {
+		/* Diffuse Ray Depth */
+		PathState *state = sd->osl_path_state;
+		int f = state->diffuse_bounce;
+		return set_attribute_int(f, type, derivatives, val);
+	}
+	else if(name == u_path_glossy_depth) {
+		/* Glossy Ray Depth */
+		PathState *state = sd->osl_path_state;
+		int f = state->glossy_bounce;
+		return set_attribute_int(f, type, derivatives, val);
+	}
+	else if(name == u_path_transmission_depth) {
+		/* Transmission Ray Depth */
+		PathState *state = sd->osl_path_state;
+		int f = state->transmission_bounce;
+		return set_attribute_int(f, type, derivatives, val);
+	}
 	else if(name == u_path_transparent_depth) {
 		/* Transparent Ray Depth */
 		PathState *state = sd->osl_path_state;
diff --git a/intern/cycles/kernel/osl/osl_services.h b/intern/cycles/kernel/osl/osl_services.h
index 0f2e02c62b0..ec34ca77115 100644
--- a/intern/cycles/kernel/osl/osl_services.h
+++ b/intern/cycles/kernel/osl/osl_services.h
@@ -165,6 +165,8 @@ public:
 	static ustring u_curve_tangent_normal;
 	static ustring u_path_ray_length;
 	static ustring u_path_ray_depth;
+	static ustring u_path_diffuse_depth;
+	static ustring u_path_glossy_depth;
 	static ustring u_path_transparent_depth;
 	static ustring u_path_transmission_depth;
 	static ustring u_trace;
diff --git a/intern/cycles/kernel/shaders/node_light_path.osl b/intern/cycles/kernel/shaders/node_light_path.osl
index a021a40467d..64fe4c20132 100644
--- a/intern/cycles/kernel/shaders/node_light_path.osl
+++ b/intern/cycles/kernel/shaders/node_light_path.osl
@@ -27,6 +27,8 @@ shader node_light_path(
 	output float IsVolumeScatterRay = 0.0,
 	output float RayLength = 0.0,
 	output float RayDepth = 0.0,
+	output float DiffuseDepth = 0.0,
+	output float GlossyDepth = 0.0,
 	output float TransparentDepth = 0.0,
 	output float TransmissionDepth = 0.0)
 {
@@ -45,6 +47,14 @@ shader node_light_path(
 	getattribute("path:ray_depth", ray_depth);
 	RayDepth = (float)ray_depth;
 
+	int diffuse_depth;
+	getattribute("path:diffuse_depth", diffuse_depth);
+	DiffuseDepth = (float)diffuse_depth;
+
+	int glossy_depth;
+	getattribute("path:glossy_depth", glossy_depth);
+	GlossyDepth = (float)glossy_depth;
+
 	int transparent_depth;
 	getattribute("path:transparent_depth", transparent_depth);
 	TransparentDepth = (float)transparent_depth;
diff --git a/intern/cycles/kernel/svm/svm_light_path.h b/intern/cycles/kernel/svm/svm_light_path.h
index f35ea05048b..04f6f623f18 100644
--- a/intern/cycles/kernel/svm/svm_light_path.h
+++ b/intern/cycles/kernel/svm/svm_light_path.h
@@ -34,6 +34,8 @@ ccl_device void svm_node_light_path(ShaderData *sd, ccl_addr_space PathState *st
 		case NODE_LP_backfacing: info = (ccl_fetch(sd, flag) & SD_BACKFACING)? 1.0f: 0.0f; break;
 		case NODE_LP_ray_length: info = ccl_fetch(sd, ray_length); break;
 		case NODE_LP_ray_depth: info = (float)state->bounce; break;
+		case NODE_LP_ray_diffuse: info = (float)state->diffuse_bounce; break;
+		case NODE_LP_ray_glossy: info = (float)state->glossy_bounce; break;
 		case NODE_LP_ray_transparent: info = (float)state->transparent_bounce; break;
 		case NODE_LP_ray_transmission: info = (float)state->transmission_bounce; break;
 	}
diff --git a/intern/cycles/kernel/svm/svm_types.h b/intern/cycles/kernel/svm/svm_types.h
index 5adf7d34f7f..47209ddfbab 100644
--- a/intern/cycles/kernel/svm/svm_types.h
+++ b/intern/cycles/kernel/svm/svm_types.h
@@ -188,6 +188,8 @@ typedef enum NodeLightPath {
 	NODE_LP_backfacing,
 	NODE_LP_ray_length,
 	NODE_LP_ray_depth,
+	NODE_LP_ray_diffuse,
+	NODE_LP_ray_glossy,
 	NODE_LP_ray_transparent,
 	NODE_LP_ray_transmission,
 } NodeLightPath;
diff --git a/intern/cycles/render/mesh.cpp b/intern/cycles/render/mesh.cpp
index df4327d021a..c42b32919d4 100644
--- a/intern/cycles/render/mesh.cpp
+++ b/intern/cycles/render/mesh.cpp
@@ -49,6 +49,64 @@ void Mesh::Triangle::bounds_grow(const float3 *verts, BoundBox& bounds) const
 	bounds.grow(verts[v[2]]);
 }
 
+void Mesh::Triangle::motion_verts(const float3 *verts,
+                                  const float3 *vert_steps,
+                                  size_t num_verts,
+                                  size_t num_steps,
+                                  float time,
+                                  float3 r_verts[3]) const
+{
+	/* Figure out which steps we need to fetch and their interpolation factor. */
+	const size_t max_step = num_steps - 1;
+	const size_t step = min((int)(time * max_step), max_step - 1);
+	const float t = time*max_step - step;
+	/* Fetch vertex coordinates. */
+	float3 curr_verts[3];
+	float3 next_verts[3];
+	verts_for_step(verts,
+	               vert_steps,
+	               num_verts,
+	               num_steps,
+	               step,
+	               curr_verts);
+	verts_for_step(verts,
+	               vert_steps,
+	               num_verts,
+	               num_steps,
+	               step + 1,
+	               next_verts);
+	/* Interpolate between steps. */
+	r_verts[0] = (1.0f - t)*curr_verts[0] + t*next_verts[0];
+	r_verts[1] = (1.0f - t)*curr_verts[1] + t*next_verts[1];
+	r_verts[2] = (1.0f - t)*curr_verts[2] + t*next_verts[2];
+}
+
+void Mesh::Triangle::verts_for_step(const float3 *verts,
+                                    const float3 *vert_steps,
+                                    size_t num_verts,
+                                    size_t num_steps,
+                                    size_t step,
+                                    float3 r_verts[3]) const
+{
+	const size_t center_step = ((num_steps - 1) / 2);
+	if(step == center_step) {
+		/* Center step: regular vertex location. */
+		r_verts[0] = verts[v[0]];
+		r_verts[1] = verts[v[1]];
+		r_verts[2] = verts[v[2]];
+	}
+	else {
+		/* Center step not stored in the attribute array array. */
+		if(step > center_step) {
+			step--;
+		}
+		size_t offset = step * num_verts;
+		r_verts[0] = vert_steps[offset + v[0]];
+		r_verts[1] = vert_steps[offset + v[1]];
+		r_verts[2] = vert_steps[offset + v[2]];
+	}
+}
+
 /* Curve */
 
 void Mesh::Curve::bounds_grow(const int k, const float3 *curve_keys, const float *curve_radius, BoundBox& bounds) const
@@ -104,6 +162,205 @@ void Mesh::Curve::bounds_grow(const int k,
 	bounds.grow(upper, mr);
 }
 
+void Mesh::Curve::bounds_grow(float4 keys[4], BoundBox& bounds) const
+{
+	float3 P[4] = {
+		float4_to_float3(keys[0]),
+		float4_to_float3(keys[1]),
+		float4_to_float3(keys[2]),
+		float4_to_float3(keys[3]),
+	};
+
+	float3 lower;
+	float3 upper;
+
+	curvebounds(&lower.x, &upper.x, P, 0);
+	curvebounds(&lower.y, &upper.y, P, 1);
+	curvebounds(&lower.z, &upper.z, P, 2);
+
+	float mr = max(keys[1].w, keys[2].w);
+
+	bounds.grow(lower, mr);
+	bounds.grow(upper, mr);
+}
+
+void Mesh::Curve::motion_keys(const float3 *curve_keys,
+                              const float *curve_radius,
+                              const float3 *key_steps,
+                              size_t num_curve_keys,
+                              size_t num_steps,
+                              float time,
+                              size_t k0, size_t k1,
+                              float4 r_keys[2]) const
+{
+	/* Figure out which steps we need to fetch and their interpolation factor. */
+	const size_t max_step = num_steps - 1;
+	const size_t step = min((int)(time * max_step), max_step - 1);
+	const float t = time*max_step - step;
+	/* Fetch vertex coordinates. */
+	float4 curr_keys[2];
+	float4 next_keys[2];
+	keys_for_step(curve_keys,
+	              curve_radius,
+	              key_steps,
+	              num_curve_keys,
+	              num_steps,
+	              step,
+	              k0, k1,
+	              curr_keys);
+	keys_for_step(curve_keys,
+	              curve_radius,
+	              key_steps,
+	              num_curve_keys,
+	              num_steps,
+	              step + 1,
+	              k0, k1,
+	              next_keys);
+	/* Interpolate between steps. */
+	r_keys[0] = (1.0f - t)*curr_keys[0] + t*next_keys[0];
+	r_keys[1] = (1.0f - t)*curr_keys[1] + t*next_keys[1];
+}
+
+void Mesh::Curve::cardinal_motion_keys(const float3 *curve_keys,
+                                       const float *curve_radius,
+                                       const float3 *key_steps,
+                                       size_t num_curve_keys,
+                                       size_t num_steps,
+                                       float time,
+                                       size_t k0, size_t k1,
+                                       size_t k2, size_t k3,
+                                       float4 r_keys[4]) const
+{
+	/* Figure out which steps we need to fetch and their interpolation factor. */
+	const size_t max_step = num_steps - 1;
+	const size_t step = min((int)(time * max_step), max_step - 1);
+	const float t = time*max_step - step;
+	/* Fetch vertex coordinates. */
+	float4 curr_keys[4];
+	float4 next_keys[4];
+	cardinal_keys_for_step(curve_keys,
+	                       curve_radius,
+	                       key_steps,
+	                       num_curve_keys,
+	                       num_steps,
+	                       step,
+	                       k0, k1, k2, k3,
+	                       curr_keys);
+	cardinal_keys_for_step(curve_keys,
+	                       curve_radius,
+	                       key_steps,
+	                       num_curve_keys,
+	                       num_steps,
+	                       step + 1,
+	                       k0, k1, k2, k3,
+	                       next_keys);
+	/* Interpolate between steps. */
+	r_keys[0] = (1.0f - t)*curr_keys[0] + t*next_keys[0];
+	r_keys[1] = (1.0f - t)*curr_keys[1] + t*next_keys[1];
+	r_keys[2] = (1.0f - t)*curr_keys[2] + t*next_keys[2];
+	r_keys[3] = (1.0f - t)*curr_keys[3] + t*next_keys[3];
+}
+
+void Mesh::Curve::keys_for_step(const float3 *curve_keys,
+                                const float *curve_radius,
+                                const float3 *key_steps,
+                                size_t num_curve_keys,
+                                size_t num_steps,
+                                size_t step,
+                                size_t k0, size_t k1,
+                                float4 r_keys[2]) const
+{
+	k0 = max(k0, 0);
+	k1 = min(k1, num_keys - 1);
+	const size_t center_step = ((num_steps - 1) / 2);
+	if(step == center_step) {
+		/* Center step: regular key location. */
+		/* TODO(sergey): Consider adding make_float4(float3, float)
+		 * function.
+		 */
+		r_keys[0] = make_float4(curve_keys[first_key + k0].x,
+		                        curve_keys[first_key + k0].y,
+		                        curve_keys[first_key + k0].z,
+		                        curve_radius[first_key + k0]);
+		r_keys[1] = make_float4(curve_keys[first_key + k1].x,
+		                        curve_keys[first_key + k1].y,
+		                        curve_keys[first_key + k1].z,
+		                        curve_radius[first_key + k1]);
+	}
+	else {
+		/* Center step is not stored in this array. */
+		if(step > center_step) {
+			step--;
+		}
+		const size_t offset = first_key + step * num_curve_keys;
+		r_keys[0] = make_float4(key_steps[offset + k0].x,
+		                        key_steps[offset + k0].y,
+		                        key_steps[offset + k0].z,
+		                        curve_radius[first_key + k0]);
+		r_keys[1] = make_float4(key_steps[offset + k1].x,
+		                        key_steps[offset + k1].y,
+		                        key_steps[offset + k1].z,
+		                        curve_radius[first_key + k1]);
+	}
+}
+
+void Mesh::Curve::cardinal_keys_for_step(const float3 *curve_keys,
+                                         const float *curve_radius,
+                                         const float3 *key_steps,
+                                         size_t num_curve_keys,
+                                         size_t num_steps,
+                                         size_t step,
+                                         size_t k0, size_t k1,
+                                         size_t k2, size_t k3,
+                                         float4 r_keys[4]) const
+{
+	k0 = max(k0, 0);
+	k3 = min(k3, num_keys - 1);
+	const size_t center_step = ((num_steps - 1) / 2);
+	if(step == center_step) {
+		/* Center step: regular key location. */
+		r_keys[0] = make_float4(curve_keys[first_key + k0].x,
+		                        curve_keys[first_key + k0].y,
+		                        curve_keys[first_key + k0].z,
+		                        curve_radius[first_key + k0]);
+		r_keys[1] = make_float4(curve_keys[first_key + k1].x,
+		                        curve_keys[first_key + k1].y,
+		                        curve_keys[first_key + k1].z,
+		                        curve_radius[first_key + k1]);
+		r_keys[2] = make_float4(curve_keys[first_key + k2].x,
+		                        curve_keys[first_key + k2].y,
+		                        curve_keys[first_key + k2].z,
+		                        curve_radius[first_key + k2]);
+		r_keys[3] = make_float4(curve_keys[first_key + k3].x,
+		                        curve_keys[first_key + k3].y,
+		                        curve_keys[first_key + k3].z,
+		                        curve_radius[first_key + k3]);
+	}
+	else {
+		/* Center step is not stored in this array. */
+		if(step > center_step) {
+			step--;
+		}
+		const size_t offset = first_key + step * num_curve_keys;
+		r_keys[0] = make_float4(key_steps[offset + k0].x,
+		                        key_steps[offset + k0].y,
+		                        key_steps[offset + k0].z,
+		                        curve_radius[first_key + k0]);
+		r_keys[1] = make_float4(key_steps[offset + k1].x,
+		                        key_steps[offset + k1].y,
+		                        key_steps[offset + k1].z,
+		                        curve_radius[first_key + k1]);
+		r_keys[2] = make_float4(key_steps[offset + k2].x,
+		                        key_steps[offset + k2].y,
+		                        key_steps[offset + k2].z,
+		                        curve_radius[first_key + k2]);
+		r_keys[3] = make_float4(key_steps[offset + k3].x,
+		                        key_steps[offset + k3].y,
+		                        key_steps[offset + k3].z,
+		                        curve_radius[first_key + k3]);
+	}
+}
+
 /* SubdFace */
 
 float3 Mesh::SubdFace::normal(const Mesh *mesh) const
@@ -394,7 +651,7 @@ void Mesh::compute_bounds()
 		if(use_motion_blur && attr) {
 			size_t steps_size = verts.size() * (motion_steps - 1);
 			float3 *vert_steps = attr->data_float3();
-	
+
 			for(size_t i = 0; i < steps_size; i++)
 				bnds.grow(vert_steps[i]);
 		}
@@ -403,7 +660,7 @@ void Mesh::compute_bounds()
 		if(use_motion_blur && curve_attr) {
 			size_t steps_size = curve_keys.size() * (motion_steps - 1);
 			float3 *key_steps = curve_attr->data_float3();
-	
+
 			for(size_t i = 0; i < steps_size; i++)
 				bnds.grow(key_steps[i]);
 		}
@@ -417,11 +674,11 @@ void Mesh::compute_bounds()
 
 			for(size_t i = 0; i < curve_keys_size; i++)
 				bnds.grow_safe(curve_keys[i], curve_radius[i]);
-			
+
 			if(use_motion_blur && attr) {
 				size_t steps_size = verts.size() * (motion_steps - 1);
 				float3 *vert_steps = attr->data_float3();
-		
+
 				for(size_t i = 0; i < steps_size; i++)
 					bnds.grow_safe(vert_steps[i]);
 			}
@@ -429,7 +686,7 @@ void Mesh::compute_bounds()
 			if(use_motion_blur && curve_attr) {
 				size_t steps_size = curve_keys.size() * (motion_steps - 1);
 				float3 *key_steps = curve_attr->data_float3();
-		
+
 				for(size_t i = 0; i < steps_size; i++)
 					bnds.grow_safe(key_steps[i]);
 			}
@@ -464,7 +721,7 @@ void Mesh::add_face_normals()
 	/* don't compute if already there */
 	if(attributes.find(ATTR_STD_FACE_NORMAL))
 		return;
-	
+
 	/* get attributes */
 	Attribute *attr_fN = attributes.add(ATTR_STD_FACE_NORMAL);
 	float3 *fN = attr_fN->data_float3();
@@ -796,6 +1053,8 @@ void Mesh::compute_bvh(DeviceScene *dscene,
 			bparams.use_qbvh = params->use_qbvh;
 			bparams.use_unaligned_nodes = dscene->data.bvh.have_curves &&
 			                              params->use_bvh_unaligned_nodes;
+			bparams.num_motion_triangle_steps = params->num_bvh_time_steps;
+			bparams.num_motion_curve_steps = params->num_bvh_time_steps;
 
 			delete bvh;
 			bvh = BVH::create(bparams, objects);
@@ -1002,7 +1261,7 @@ void MeshManager::update_svm_attributes(Device *device, DeviceScene *dscene, Sce
 
 	if(attr_map_stride == 0)
 		return;
-	
+
 	/* create attribute map */
 	uint4 *attr_map = dscene->attributes_map.resize(attr_map_stride*scene->objects.size());
 	memset(attr_map, 0, dscene->attributes_map.size()*sizeof(uint));
@@ -1564,6 +1823,8 @@ void MeshManager::device_update_bvh(Device *device, DeviceScene *dscene, Scene *
 	bparams.use_spatial_split = scene->params.use_bvh_spatial_split;
 	bparams.use_unaligned_nodes = dscene->data.bvh.have_curves &&
 	                              scene->params.use_bvh_unaligned_nodes;
+	bparams.num_motion_triangle_steps = scene->params.num_bvh_time_steps;
+	bparams.num_motion_curve_steps = scene->params.num_bvh_time_steps;
 
 	delete bvh;
 	bvh = BVH::create(bparams, scene->objects);
@@ -1946,14 +2207,14 @@ bool Mesh::need_attribute(Scene *scene, AttributeStandard std)
 {
 	if(std == ATTR_STD_NONE)
 		return false;
-	
+
 	if(scene->need_global_attribute(std))
 		return true;
 
 	foreach(Shader *shader, used_shaders)
 		if(shader->attributes.find(std))
 			return true;
-	
+
 	return false;
 }
 
@@ -1965,9 +2226,8 @@ bool Mesh::need_attribute(Scene * /*scene*/, ustring name)
 	foreach(Shader *shader, used_shaders)
 		if(shader->attributes.find(name))
 			return true;
-	
+
 	return false;
 }
 
 CCL_NAMESPACE_END
-
diff --git a/intern/cycles/render/mesh.h b/intern/cycles/render/mesh.h
index c0310f45840..5f33e30eac2 100644
--- a/intern/cycles/render/mesh.h
+++ b/intern/cycles/render/mesh.h
@@ -31,6 +31,7 @@
 
 CCL_NAMESPACE_BEGIN
 
+class Attribute;
 class BVH;
 class Device;
 class DeviceScene;
@@ -54,11 +55,27 @@ public:
 		int v[3];
 
 		void bounds_grow(const float3 *verts, BoundBox& bounds) const;
+
+		void motion_verts(const float3 *verts,
+		                  const float3 *vert_steps,
+		                  size_t num_verts,
+		                  size_t num_steps,
+		                  float time,
+		                  float3 r_verts[3]) const;
+
+		void verts_for_step(const float3 *verts,
+		                    const float3 *vert_steps,
+		                    size_t num_verts,
+		                    size_t num_steps,
+		                    size_t step,
+		                    float3 r_verts[3]) const;
 	};
 
 	Triangle get_triangle(size_t i) const
 	{
-		Triangle tri = {{triangles[i*3 + 0], triangles[i*3 + 1], triangles[i*3 + 2]}};
+		Triangle tri = {{triangles[i*3 + 0],
+		                 triangles[i*3 + 1],
+		                 triangles[i*3 + 2]}};
 		return tri;
 	}
 
@@ -78,11 +95,48 @@ public:
 		                 const float3 *curve_keys,
 		                 const float *curve_radius,
 		                 BoundBox& bounds) const;
+		void bounds_grow(float4 keys[4], BoundBox& bounds) const;
 		void bounds_grow(const int k,
 		                 const float3 *curve_keys,
 		                 const float *curve_radius,
 		                 const Transform& aligned_space,
 		                 BoundBox& bounds) const;
+
+		void motion_keys(const float3 *curve_keys,
+		                 const float *curve_radius,
+		                 const float3 *key_steps,
+		                 size_t num_curve_keys,
+		                 size_t num_steps,
+		                 float time,
+		                 size_t k0, size_t k1,
+		                 float4 r_keys[2]) const;
+		void cardinal_motion_keys(const float3 *curve_keys,
+		                          const float *curve_radius,
+		                          const float3 *key_steps,
+		                          size_t num_curve_keys,
+		                          size_t num_steps,
+		                          float time,
+		                          size_t k0, size_t k1,
+		                          size_t k2, size_t k3,
+		                          float4 r_keys[4]) const;
+
+		void keys_for_step(const float3 *curve_keys,
+		                   const float *curve_radius,
+		                   const float3 *key_steps,
+		                   size_t num_curve_keys,
+		                   size_t num_steps,
+		                   size_t step,
+		                   size_t k0, size_t k1,
+		                   float4 r_keys[2]) const;
+		void cardinal_keys_for_step(const float3 *curve_keys,
+		                            const float *curve_radius,
+		                            const float3 *key_steps,
+		                            size_t num_curve_keys,
+		                            size_t num_steps,
+		                            size_t step,
+		                            size_t k0, size_t k1,
+		                            size_t k2, size_t k3,
+		                            float4 r_keys[4]) const;
 	};
 
 	Curve get_curve(size_t i) const
diff --git a/intern/cycles/render/mesh_subdivision.cpp b/intern/cycles/render/mesh_subdivision.cpp
index 913c3c74b42..57c76a9f1c8 100644
--- a/intern/cycles/render/mesh_subdivision.cpp
+++ b/intern/cycles/render/mesh_subdivision.cpp
@@ -92,7 +92,7 @@ namespace Far {
 
 			if(vert_edges.size() == 2) {
 				float sharpness = refiner.getLevel(0).getEdgeSharpness(vert_edges[0]);
-				sharpness = min(sharpness, refiner.getLevel(0).getEdgeSharpness(vert_edges[1]));
+				sharpness = ccl::min(sharpness, refiner.getLevel(0).getEdgeSharpness(vert_edges[1]));
 
 				setBaseVertexSharpness(refiner, i, sharpness);
 			}
diff --git a/intern/cycles/render/nodes.cpp b/intern/cycles/render/nodes.cpp
index c7f37a13fba..1e4a9fd300c 100644
--- a/intern/cycles/render/nodes.cpp
+++ b/intern/cycles/render/nodes.cpp
@@ -3027,6 +3027,8 @@ NODE_DEFINE(LightPathNode)
 	SOCKET_OUT_FLOAT(is_volume_scatter_ray, "Is Volume Scatter Ray");
 	SOCKET_OUT_FLOAT(ray_length, "Ray Length");
 	SOCKET_OUT_FLOAT(ray_depth, "Ray Depth");
+	SOCKET_OUT_FLOAT(diffuse_depth, "Diffuse Depth");
+	SOCKET_OUT_FLOAT(glossy_depth, "Glossy Depth");
 	SOCKET_OUT_FLOAT(transparent_depth, "Transparent Depth");
 	SOCKET_OUT_FLOAT(transmission_depth, "Transmission Depth");
 
@@ -3093,6 +3095,16 @@ void LightPathNode::compile(SVMCompiler& compiler)
 		compiler.add_node(NODE_LIGHT_PATH, NODE_LP_ray_depth, compiler.stack_assign(out));
 	}
 
+	out = output("Diffuse Depth");
+	if(!out->links.empty()) {
+		compiler.add_node(NODE_LIGHT_PATH, NODE_LP_ray_diffuse, compiler.stack_assign(out));
+	}
+
+	out = output("Glossy Depth");
+	if(!out->links.empty()) {
+		compiler.add_node(NODE_LIGHT_PATH, NODE_LP_ray_glossy, compiler.stack_assign(out));
+	}
+
 	out = output("Transparent Depth");
 	if(!out->links.empty()) {
 		compiler.add_node(NODE_LIGHT_PATH, NODE_LP_ray_transparent, compiler.stack_assign(out));
diff --git a/intern/cycles/render/object.cpp b/intern/cycles/render/object.cpp
index 8b8b988b969..c592b62829e 100644
--- a/intern/cycles/render/object.cpp
+++ b/intern/cycles/render/object.cpp
@@ -166,7 +166,7 @@ void Object::apply_transform(bool apply_to_motion)
 		float3 c0 = transform_get_column(&tfm, 0);
 		float3 c1 = transform_get_column(&tfm, 1);
 		float3 c2 = transform_get_column(&tfm, 2);
-		float scalar = pow(fabsf(dot(cross(c0, c1), c2)), 1.0f/3.0f);
+		float scalar = powf(fabsf(dot(cross(c0, c1), c2)), 1.0f/3.0f);
 
 		/* apply transform to curve keys */
 		for(size_t i = 0; i < mesh->curve_keys.size(); i++) {
diff --git a/intern/cycles/render/scene.h b/intern/cycles/render/scene.h
index df9363cc768..8768682043f 100644
--- a/intern/cycles/render/scene.h
+++ b/intern/cycles/render/scene.h
@@ -143,6 +143,7 @@ public:
 	} bvh_type;
 	bool use_bvh_spatial_split;
 	bool use_bvh_unaligned_nodes;
+	int num_bvh_time_steps;
 	bool use_qbvh;
 	bool persistent_data;
 	int texture_limit;
@@ -153,6 +154,7 @@ public:
 		bvh_type = BVH_DYNAMIC;
 		use_bvh_spatial_split = false;
 		use_bvh_unaligned_nodes = true;
+		num_bvh_time_steps = 0;
 		use_qbvh = false;
 		persistent_data = false;
 		texture_limit = 0;
@@ -163,6 +165,7 @@ public:
 		&& bvh_type == params.bvh_type
 		&& use_bvh_spatial_split == params.use_bvh_spatial_split
 		&& use_bvh_unaligned_nodes == params.use_bvh_unaligned_nodes
+		&& num_bvh_time_steps == params.num_bvh_time_steps
 		&& use_qbvh == params.use_qbvh
 		&& persistent_data == params.persistent_data
 		&& texture_limit == params.texture_limit); }
diff --git a/intern/cycles/render/session.cpp b/intern/cycles/render/session.cpp
index 73caf93ea00..7c01934cfd8 100644
--- a/intern/cycles/render/session.cpp
+++ b/intern/cycles/render/session.cpp
@@ -458,6 +458,8 @@ void Session::release_tile(RenderTile& rtile)
 {
 	thread_scoped_lock tile_lock(tile_mutex);
 
+	progress.add_finished_tile();
+
 	if(write_render_tile_cb) {
 		if(params.progressive_refine == false) {
 			/* todo: optimize this by making it thread safe and removing lock */
diff --git a/intern/cycles/util/util_boundbox.h b/intern/cycles/util/util_boundbox.h
index 599222da9c5..dfe4977aef3 100644
--- a/intern/cycles/util/util_boundbox.h
+++ b/intern/cycles/util/util_boundbox.h
@@ -25,8 +25,6 @@
 #include "util_transform.h"
 #include "util_types.h"
 
-using namespace std;
-
 CCL_NAMESPACE_BEGIN
 
 /* 3D BoundBox */
diff --git a/intern/cycles/util/util_math.h b/intern/cycles/util/util_math.h
index 6cb68b53d16..2b81c8c498a 100644
--- a/intern/cycles/util/util_math.h
+++ b/intern/cycles/util/util_math.h
@@ -22,6 +22,11 @@
  * Basic math functions on scalar and vector types. This header is used by
  * both the kernel code when compiled as C++, and other C++ non-kernel code. */
 
+#ifndef __KERNEL_GPU__
+#  include <cmath>
+#endif
+
+
 #ifndef __KERNEL_OPENCL__
 
 #include <float.h>
@@ -97,6 +102,9 @@ ccl_device_inline float fminf(float a, float b)
 
 #ifndef __KERNEL_GPU__
 
+using std::isfinite;
+using std::isnan;
+
 ccl_device_inline int abs(int x)
 {
 	return (x > 0)? x: -x;
@@ -1233,6 +1241,20 @@ ccl_device_inline float __uint_as_float(uint i)
 	return u.f;
 }
 
+/* Versions of functions which are safe for fast math. */
+ccl_device_inline bool isnan_safe(float f)
+{
+	unsigned int x = __float_as_uint(f);
+	return (x << 1) > 0xff000000u;
+}
+
+ccl_device_inline bool isfinite_safe(float f)
+{
+	/* By IEEE 754 rule, 2*Inf equals Inf */
+	unsigned int x = __float_as_uint(f);
+	return (f == f) && (x == 0 || (f != 2.0f*f));
+}
+
 /* Interpolation */
 
 template<class A, class B> A lerp(const A& a, const A& b, const B& t)