Cycles: add bevel shader, for raytrace based rounded edges.

The algorithm averages normals from nearby surfaces. It uses the same
sampling strategy as BSSRDFs, casting rays along the normal and two
orthogonal axes, and combining the samples with MIS.

The main concern here is that we are introducing raytracing inside
shader evaluation, which could be quite bad for GPU performance and
stack memory usage. In practice it doesn't seem so bad though.

Note that using this feature can easily slow down renders 20%, and
that if you care about performance then it's better to use a bevel
modifier. Mainly this is useful for baking, and for cases where the
mesh topology makes it difficult for the bevel modifier to work well.

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

16 files changed, 438 insertions, 59 deletions

diff --git a/intern/cycles/blender/blender_shader.cpp b/intern/cycles/blender/blender_shader.cpp
index cd6c9f319db..a90e0d3ddc4 100644
--- a/intern/cycles/blender/blender_shader.cpp
+++ b/intern/cycles/blender/blender_shader.cpp
@@ -866,6 +866,12 @@ static ShaderNode *add_node(Scene *scene,
 			        transform_inverse(get_transform(b_ob.matrix_world()));
 		}
 	}
+	else if(b_node.is_a(&RNA_ShaderNodeBevel)) {
+		BL::ShaderNodeBevel b_bevel_node(b_node);
+		BevelNode *bevel = new BevelNode();
+		bevel->samples = b_bevel_node.samples();
+		node = bevel;
+	}
 
 	if(node) {
 		node->name = b_node.name();
diff --git a/intern/cycles/kernel/CMakeLists.txt b/intern/cycles/kernel/CMakeLists.txt
index 9d52cef1f2c..de056ce97f0 100644
--- a/intern/cycles/kernel/CMakeLists.txt
+++ b/intern/cycles/kernel/CMakeLists.txt
@@ -160,6 +160,7 @@ set(SRC_CLOSURE_HEADERS
 set(SRC_SVM_HEADERS
 	svm/svm.h
 	svm/svm_attribute.h
+	svm/svm_bevel.h
 	svm/svm_blackbody.h
 	svm/svm_bump.h
 	svm/svm_camera.h
diff --git a/intern/cycles/kernel/closure/bsdf.h b/intern/cycles/kernel/closure/bsdf.h
index 86a00d2124d..0a3f9ff23fe 100644
--- a/intern/cycles/kernel/closure/bsdf.h
+++ b/intern/cycles/kernel/closure/bsdf.h
@@ -29,9 +29,7 @@
 #include "kernel/closure/bsdf_hair.h"
 #include "kernel/closure/bsdf_principled_diffuse.h"
 #include "kernel/closure/bsdf_principled_sheen.h"
-#ifdef __SUBSURFACE__
-#  include "kernel/closure/bssrdf.h"
-#endif
+#include "kernel/closure/bssrdf.h"
 #ifdef __VOLUME__
 #  include "kernel/closure/volume.h"
 #endif
diff --git a/intern/cycles/kernel/closure/bssrdf.h b/intern/cycles/kernel/closure/bssrdf.h
index 267aeea6e86..6791c0b83cc 100644
--- a/intern/cycles/kernel/closure/bssrdf.h
+++ b/intern/cycles/kernel/closure/bssrdf.h
@@ -39,12 +39,11 @@ typedef ccl_addr_space struct Bssrdf {
 /* paper suggests 1/12.46 which is much too small, suspect it's *12.46 */
 #define GAUSS_TRUNCATE 12.46f
 
-ccl_device float bssrdf_gaussian_eval(const ShaderClosure *sc, float r)
+ccl_device float bssrdf_gaussian_eval(const float radius, float r)
 {
 	/* integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) from 0 to Rm
 	 * = 1 - exp(-Rm*Rm/(2*v)) */
-	const Bssrdf *bssrdf = (const Bssrdf*)sc;
-	const float v = bssrdf->radius*bssrdf->radius*(0.25f*0.25f);
+	const float v = radius*radius*(0.25f*0.25f);
 	const float Rm = sqrtf(v*GAUSS_TRUNCATE);
 
 	if(r >= Rm)
@@ -53,20 +52,19 @@ ccl_device float bssrdf_gaussian_eval(const ShaderClosure *sc, float r)
 	return expf(-r*r/(2.0f*v))/(2.0f*M_PI_F*v);
 }
 
-ccl_device float bssrdf_gaussian_pdf(const ShaderClosure *sc, float r)
+ccl_device float bssrdf_gaussian_pdf(const float radius, float r)
 {
 	/* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
 	const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
 
-	return bssrdf_gaussian_eval(sc, r) * (1.0f/(area_truncated));
+	return bssrdf_gaussian_eval(radius, r) * (1.0f/(area_truncated));
 }
 
-ccl_device void bssrdf_gaussian_sample(const ShaderClosure *sc, float xi, float *r, float *h)
+ccl_device void bssrdf_gaussian_sample(const float radius, float xi, float *r, float *h)
 {
 	/* xi = integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) = -exp(-r^2/(2*v))
 	 * r = sqrt(-2*v*logf(xi)) */
-	const Bssrdf *bssrdf = (const Bssrdf*)sc;
-	const float v = bssrdf->radius*bssrdf->radius*(0.25f*0.25f);
+	const float v = radius*radius*(0.25f*0.25f);
 	const float Rm = sqrtf(v*GAUSS_TRUNCATE);
 
 	/* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
@@ -87,13 +85,10 @@ ccl_device void bssrdf_gaussian_sample(const ShaderClosure *sc, float xi, float
  * far as I can tell has no closed form solution. So we get an iterative solution
  * instead with newton-raphson. */
 
-ccl_device float bssrdf_cubic_eval(const ShaderClosure *sc, float r)
+ccl_device float bssrdf_cubic_eval(const float radius, const float sharpness, float r)
 {
-	const Bssrdf *bssrdf = (const Bssrdf*)sc;
-	const float sharpness = bssrdf->sharpness;
-
 	if(sharpness == 0.0f) {
-		const float Rm = bssrdf->radius;
+		const float Rm = radius;
 
 		if(r >= Rm)
 			return 0.0f;
@@ -107,7 +102,7 @@ ccl_device float bssrdf_cubic_eval(const ShaderClosure *sc, float r)
 
 	}
 	else {
-		float Rm = bssrdf->radius*(1.0f + sharpness);
+		float Rm = radius*(1.0f + sharpness);
 
 		if(r >= Rm)
 			return 0.0f;
@@ -135,9 +130,9 @@ ccl_device float bssrdf_cubic_eval(const ShaderClosure *sc, float r)
 	}
 }
 
-ccl_device float bssrdf_cubic_pdf(const ShaderClosure *sc, float r)
+ccl_device float bssrdf_cubic_pdf(const float radius, const float sharpness, float r)
 {
-	return bssrdf_cubic_eval(sc, r);
+	return bssrdf_cubic_eval(radius, sharpness, r);
 }
 
 /* solve 10x^2 - 20x^3 + 15x^4 - 4x^5 - xi == 0 */
@@ -168,11 +163,9 @@ ccl_device_forceinline float bssrdf_cubic_quintic_root_find(float xi)
 	return x;
 }
 
-ccl_device void bssrdf_cubic_sample(const ShaderClosure *sc, float xi, float *r, float *h)
+ccl_device void bssrdf_cubic_sample(const float radius, const float sharpness, float xi, float *r, float *h)
 {
-	const Bssrdf *bssrdf = (const Bssrdf*)sc;
-	const float sharpness = bssrdf->sharpness;
-	float Rm = bssrdf->radius;
+	float Rm = radius;
 	float r_ = bssrdf_cubic_quintic_root_find(xi);
 
 	if(sharpness != 0.0f) {
@@ -224,10 +217,8 @@ ccl_device void bssrdf_burley_setup(Bssrdf *bssrdf)
 	bssrdf->d = d;
 }
 
-ccl_device float bssrdf_burley_eval(const ShaderClosure *sc, float r)
+ccl_device float bssrdf_burley_eval(const float d, float r)
 {
-	const Bssrdf *bssrdf = (const Bssrdf*)sc;
-	const float d = bssrdf->d;
 	const float Rm = BURLEY_TRUNCATE * d;
 
 	if(r >= Rm)
@@ -246,9 +237,9 @@ ccl_device float bssrdf_burley_eval(const ShaderClosure *sc, float r)
 	return (exp_r_d + exp_r_3_d) / (4.0f*d);
 }
 
-ccl_device float bssrdf_burley_pdf(const ShaderClosure *sc, float r)
+ccl_device float bssrdf_burley_pdf(const float d, float r)
 {
-	return bssrdf_burley_eval(sc, r) * (1.0f/BURLEY_TRUNCATE_CDF);
+	return bssrdf_burley_eval(d, r) * (1.0f/BURLEY_TRUNCATE_CDF);
 }
 
 /* Find the radius for desired CDF value.
@@ -291,13 +282,11 @@ ccl_device_forceinline float bssrdf_burley_root_find(float xi)
 	return r;
 }
 
-ccl_device void bssrdf_burley_sample(const ShaderClosure *sc,
+ccl_device void bssrdf_burley_sample(const float d,
                                      float xi,
                                      float *r,
                                      float *h)
 {
-	const Bssrdf *bssrdf = (const Bssrdf*)sc;
-	const float d = bssrdf->d;
 	const float Rm = BURLEY_TRUNCATE * d;
 	const float r_ = bssrdf_burley_root_find(xi * BURLEY_TRUNCATE_CDF) * d;
 
@@ -311,29 +300,26 @@ ccl_device void bssrdf_burley_sample(const ShaderClosure *sc,
  *
  * Samples distributed over disk with no falloff, for reference. */
 
-ccl_device float bssrdf_none_eval(const ShaderClosure *sc, float r)
+ccl_device float bssrdf_none_eval(const float radius, float r)
 {
-	const Bssrdf *bssrdf = (const Bssrdf*)sc;
-	const float Rm = bssrdf->radius;
+	const float Rm = radius;
 	return (r < Rm)? 1.0f: 0.0f;
 }
 
-ccl_device float bssrdf_none_pdf(const ShaderClosure *sc, float r)
+ccl_device float bssrdf_none_pdf(const float radius, float r)
 {
 	/* integrate (2*pi*r)/(pi*Rm*Rm) from 0 to Rm = 1 */
-	const Bssrdf *bssrdf = (const Bssrdf*)sc;
-	const float Rm = bssrdf->radius;
+	const float Rm = radius;
 	const float area = (M_PI_F*Rm*Rm);
 
-	return bssrdf_none_eval(sc, r) / area;
+	return bssrdf_none_eval(radius, r) / area;
 }
 
-ccl_device void bssrdf_none_sample(const ShaderClosure *sc, float xi, float *r, float *h)
+ccl_device void bssrdf_none_sample(const float radius, float xi, float *r, float *h)
 {
 	/* xi = integrate (2*pi*r)/(pi*Rm*Rm) = r^2/Rm^2
 	 * r = sqrt(xi)*Rm */
-	const Bssrdf *bssrdf = (const Bssrdf*)sc;
-	const float Rm = bssrdf->radius;
+	const float Rm = radius;
 	const float r_ = sqrtf(xi)*Rm;
 
 	*r = r_;
@@ -406,22 +392,26 @@ ccl_device int bssrdf_setup(Bssrdf *bssrdf, ClosureType type)
 
 ccl_device void bssrdf_sample(const ShaderClosure *sc, float xi, float *r, float *h)
 {
+	const Bssrdf *bssrdf = (const Bssrdf*)sc;
+
 	if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
-		bssrdf_cubic_sample(sc, xi, r, h);
+		bssrdf_cubic_sample(bssrdf->radius, bssrdf->sharpness, xi, r, h);
 	else if(sc->type == CLOSURE_BSSRDF_GAUSSIAN_ID)
-		bssrdf_gaussian_sample(sc, xi, r, h);
+		bssrdf_gaussian_sample(bssrdf->radius, xi, r, h);
 	else /*if(sc->type == CLOSURE_BSSRDF_BURLEY_ID || sc->type == CLOSURE_BSSRDF_PRINCIPLED_ID)*/
-		bssrdf_burley_sample(sc, xi, r, h);
+		bssrdf_burley_sample(bssrdf->d, xi, r, h);
 }
 
 ccl_device_forceinline float bssrdf_pdf(const ShaderClosure *sc, float r)
 {
+	const Bssrdf *bssrdf = (const Bssrdf*)sc;
+
 	if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
-		return bssrdf_cubic_pdf(sc, r);
+		return bssrdf_cubic_pdf(bssrdf->radius, bssrdf->sharpness, r);
 	else if(sc->type == CLOSURE_BSSRDF_GAUSSIAN_ID)
-		return bssrdf_gaussian_pdf(sc, r);
+		return bssrdf_gaussian_pdf(bssrdf->radius, r);
 	else /*if(sc->type == CLOSURE_BSSRDF_BURLEY_ID || sc->type == CLOSURE_BSSRDF_PRINCIPLED_ID)*/
-		return bssrdf_burley_pdf(sc, r);
+		return bssrdf_burley_pdf(bssrdf->d, r);
 }
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/geom/geom_motion_triangle.h b/intern/cycles/kernel/geom/geom_motion_triangle.h
index cd28b75c22c..7ac6807e749 100644
--- a/intern/cycles/kernel/geom/geom_motion_triangle.h
+++ b/intern/cycles/kernel/geom/geom_motion_triangle.h
@@ -117,4 +117,39 @@ ccl_device_inline void motion_triangle_vertices(KernelGlobals *kg, int object, i
 	verts[2] = (1.0f - t)*verts[2] + t*next_verts[2];
 }
 
+ccl_device_inline float3 motion_triangle_smooth_normal(KernelGlobals *kg, float3 Ng, int object, int prim, float u, float v, float time)
+{
+	/* get motion info */
+	int numsteps, numverts;
+	object_motion_info(kg, object, &numsteps, &numverts, NULL);
+
+	/* figure out which steps we need to fetch and their interpolation factor */
+	int maxstep = numsteps*2;
+	int step = min((int)(time*maxstep), maxstep-1);
+	float t = time*maxstep - step;
+
+	/* find attribute */
+	AttributeElement elem;
+	int offset = find_attribute_motion(kg, object, ATTR_STD_MOTION_VERTEX_NORMAL, &elem);
+	kernel_assert(offset != ATTR_STD_NOT_FOUND);
+
+	/* fetch normals */
+	float3 normals[3], next_normals[3];
+	uint4 tri_vindex = kernel_tex_fetch(__tri_vindex, prim);
+
+	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 w = 1.0f - u - v;
+	float3 N = safe_normalize(u*normals[0] + v*normals[1] + w*normals[2]);
+
+	return is_zero(N)? Ng: N;
+}
+
 CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/kernel_bake.h b/intern/cycles/kernel/kernel_bake.h
index 5cb32545c6e..73cddeb27f7 100644
--- a/intern/cycles/kernel/kernel_bake.h
+++ b/intern/cycles/kernel/kernel_bake.h
@@ -316,6 +316,13 @@ ccl_device void kernel_bake_evaluate(KernelGlobals *kg, ccl_global uint4 *input,
 	sd.dv.dx = dvdx;
 	sd.dv.dy = dvdy;
 
+	/* set RNG state for shaders that use sampling */
+	state.rng_hash = rng_hash;
+	state.rng_offset = 0;
+	state.sample = sample;
+	state.num_samples = num_samples;
+	state.min_ray_pdf = FLT_MAX;
+
 	/* light passes if we need more than color */
 	if(pass_filter & ~BAKE_FILTER_COLOR)
 		compute_light_pass(kg, &sd, &L, rng_hash, pass_filter, sample);
diff --git a/intern/cycles/kernel/kernel_types.h b/intern/cycles/kernel/kernel_types.h
index 73cf63d42be..fc3e7b3da98 100644
--- a/intern/cycles/kernel/kernel_types.h
+++ b/intern/cycles/kernel/kernel_types.h
@@ -306,6 +306,9 @@ enum PathTraceDimension {
 	PRNG_PHASE_CHANNEL = 6,
 	PRNG_SCATTER_DISTANCE = 7,
 	PRNG_BOUNCE_NUM = 8,
+
+	PRNG_BEVEL_U = 6, /* reuse volume dimension, correlation won't harm */
+	PRNG_BEVEL_V = 7,
 };
 
 enum SamplingPattern {
diff --git a/intern/cycles/kernel/osl/osl_services.cpp b/intern/cycles/kernel/osl/osl_services.cpp
index 8ae004031e1..3789073f344 100644
--- a/intern/cycles/kernel/osl/osl_services.cpp
+++ b/intern/cycles/kernel/osl/osl_services.cpp
@@ -117,6 +117,7 @@ ustring OSLRenderServices::u_I("I");
 ustring OSLRenderServices::u_u("u");
 ustring OSLRenderServices::u_v("v");
 ustring OSLRenderServices::u_empty;
+ustring OSLRenderServices::u_at_bevel("@bevel");
 
 OSLRenderServices::OSLRenderServices()
 {
@@ -958,20 +959,36 @@ bool OSLRenderServices::texture(ustring filename,
 		return true;
 	}
 #endif
-	bool status;
+	bool status = false;
 
 	if(filename.length() && filename[0] == '@') {
-		int slot = atoi(filename.c_str() + 1);
-		float4 rgba = kernel_tex_image_interp(kg, slot, s, 1.0f - t);
-
-		result[0] = rgba[0];
-		if(nchannels > 1)
-			result[1] = rgba[1];
-		if(nchannels > 2)
-			result[2] = rgba[2];
-		if(nchannels > 3)
-			result[3] = rgba[3];
-		status = true;
+		if(filename == u_at_bevel) {
+			/* Bevel shader hack. */
+			if(nchannels >= 3) {
+				PathState *state = sd->osl_path_state;
+				int num_samples = (int)s;
+				float radius = t;
+				float3 N = svm_bevel(kg, sd, state, radius, num_samples);
+				result[0] = N.x;
+				result[1] = N.y;
+				result[2] = N.z;
+				status = true;
+			}
+		}
+		else {
+			/* Packed texture. */
+			int slot = atoi(filename.c_str() + 1);
+			float4 rgba = kernel_tex_image_interp(kg, slot, s, 1.0f - t);
+
+			result[0] = rgba[0];
+			if(nchannels > 1)
+				result[1] = rgba[1];
+			if(nchannels > 2)
+				result[2] = rgba[2];
+			if(nchannels > 3)
+				result[3] = rgba[3];
+			status = true;
+		}
 	}
 	else {
 		if(texture_handle != NULL) {
diff --git a/intern/cycles/kernel/osl/osl_services.h b/intern/cycles/kernel/osl/osl_services.h
index ec34ca77115..642529655bf 100644
--- a/intern/cycles/kernel/osl/osl_services.h
+++ b/intern/cycles/kernel/osl/osl_services.h
@@ -179,6 +179,7 @@ public:
 	static ustring u_u;
 	static ustring u_v;
 	static ustring u_empty;
+	static ustring u_at_bevel;
 
 private:
 	KernelGlobals *kernel_globals;
diff --git a/intern/cycles/kernel/shaders/CMakeLists.txt b/intern/cycles/kernel/shaders/CMakeLists.txt
index 1a8ed4c884a..83cfdbcf8ba 100644
--- a/intern/cycles/kernel/shaders/CMakeLists.txt
+++ b/intern/cycles/kernel/shaders/CMakeLists.txt
@@ -7,6 +7,7 @@ set(SRC_OSL
 	node_anisotropic_bsdf.osl
 	node_attribute.osl
 	node_background.osl
+	node_bevel.osl
 	node_brick_texture.osl
 	node_brightness.osl
 	node_bump.osl
diff --git a/intern/cycles/kernel/shaders/node_bevel.osl b/intern/cycles/kernel/shaders/node_bevel.osl
new file mode 100644
index 00000000000..a5b185b6b4c
--- /dev/null
+++ b/intern/cycles/kernel/shaders/node_bevel.osl
@@ -0,0 +1,31 @@
+/*
+ * Copyright 2011-2013 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.
+ */
+
+#include "stdosl.h"
+
+shader node_bevel(
+	int samples = 4,
+	float Radius = 0.05,
+	normal NormalIn = N,
+	output normal NormalOut = N)
+{
+	/* Abuse texture call with special @bevel token. */
+	NormalOut = (normal)(color)texture("@bevel", samples, Radius);
+
+	/* Preserve input normal. */
+	NormalOut = normalize(N + (NormalOut - NormalIn));
+}
+
diff --git a/intern/cycles/kernel/svm/svm.h b/intern/cycles/kernel/svm/svm.h
index f0b3adcdad5..9ff02c1586b 100644
--- a/intern/cycles/kernel/svm/svm.h
+++ b/intern/cycles/kernel/svm/svm.h
@@ -183,6 +183,10 @@ CCL_NAMESPACE_END
 #include "kernel/svm/svm_voxel.h"
 #include "kernel/svm/svm_bump.h"
 
+#ifdef __SHADER_RAYTRACE__
+#  include "kernel/svm/svm_bevel.h"
+#endif
+
 CCL_NAMESPACE_BEGIN
 
 #define NODES_GROUP(group) ((group) <= __NODES_MAX_GROUP__)
@@ -464,6 +468,11 @@ ccl_device_noinline void svm_eval_nodes(KernelGlobals *kg, ShaderData *sd, ccl_a
 				svm_node_tex_voxel(kg, sd, stack, node, &offset);
 				break;
 #  endif  /* NODES_FEATURE(NODE_FEATURE_VOLUME) */
+#  ifdef __SHADER_RAYTRACE__
+			case NODE_BEVEL:
+				svm_node_bevel(kg, sd, state, stack, node);
+				break;
+#  endif  /* __SHADER_RAYTRACE__ */
 #endif  /* NODES_GROUP(NODE_GROUP_LEVEL_3) */
 			case NODE_END:
 				return;
diff --git a/intern/cycles/kernel/svm/svm_bevel.h b/intern/cycles/kernel/svm/svm_bevel.h
new file mode 100644
index 00000000000..65afe1f74ec
--- /dev/null
+++ b/intern/cycles/kernel/svm/svm_bevel.h
@@ -0,0 +1,227 @@
+/*
+ * Copyright 2011-2013 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.
+ */
+
+CCL_NAMESPACE_BEGIN
+
+/* Bevel shader averaging normals from nearby surfaces.
+ *
+ * Sampling strategy from: BSSRDF Importance Sampling, SIGGRAPH 2013
+ * http://library.imageworks.com/pdfs/imageworks-library-BSSRDF-sampling.pdf
+ */
+
+ccl_device_noinline float3 svm_bevel(
+	KernelGlobals *kg,
+	ShaderData *sd,
+	ccl_addr_space PathState *state,
+	float radius,
+	int num_samples)
+{
+	/* Early out if no sampling needed. */
+	if(radius <= 0.0f || num_samples < 1 || sd->object == OBJECT_NONE) {
+		return sd->N;
+	}
+
+	/* Don't bevel for blurry indirect rays. */
+	if(state->min_ray_pdf < 8.0f) {
+		return sd->N;
+	}
+
+	/* Setup for multi intersection. */
+	LocalIntersection isect;
+	uint lcg_state = lcg_state_init_addrspace(state, 0x64c6a40e);
+
+	/* Sample normals from surrounding points on surface. */
+	float3 sum_N = make_float3(0.0f, 0.0f, 0.0f);
+
+	for(int sample = 0; sample < num_samples; sample++) {
+		float disk_u, disk_v;
+		path_branched_rng_2D(kg, state->rng_hash, state, sample, num_samples,
+		                     PRNG_BEVEL_U, &disk_u, &disk_v);
+
+		/* Pick random axis in local frame and point on disk. */
+		float3 disk_N, disk_T, disk_B;
+		float pick_pdf_N, pick_pdf_T, pick_pdf_B;
+
+		disk_N = sd->Ng;
+		make_orthonormals(disk_N, &disk_T, &disk_B);
+
+		float axisu = disk_u;
+
+		if(axisu < 0.5f) {
+			pick_pdf_N = 0.5f;
+			pick_pdf_T = 0.25f;
+			pick_pdf_B = 0.25f;
+			disk_u *= 2.0f;
+		}
+		else if(axisu < 0.75f) {
+			float3 tmp = disk_N;
+			disk_N = disk_T;
+			disk_T = tmp;
+			pick_pdf_N = 0.25f;
+			pick_pdf_T = 0.5f;
+			pick_pdf_B = 0.25f;
+			disk_u = (disk_u - 0.5f)*4.0f;
+		}
+		else {
+			float3 tmp = disk_N;
+			disk_N = disk_B;
+			disk_B = tmp;
+			pick_pdf_N = 0.25f;
+			pick_pdf_T = 0.25f;
+			pick_pdf_B = 0.5f;
+			disk_u = (disk_u - 0.75f)*4.0f;
+		}
+
+		/* Sample point on disk. */
+		float phi = M_2PI_F * disk_u;
+		float disk_r = disk_v;
+		float disk_height;
+
+		/* Perhaps find something better than Cubic BSSRDF, but happens to work well. */
+		bssrdf_cubic_sample(radius, 0.0f, disk_r, &disk_r, &disk_height);
+
+		float3 disk_P = (disk_r*cosf(phi)) * disk_T + (disk_r*sinf(phi)) * disk_B;
+
+		/* Create ray. */
+		Ray *ray = &isect.ray;
+		ray->P = sd->P + disk_N*disk_height + disk_P;
+		ray->D = -disk_N;
+		ray->t = 2.0f*disk_height;
+		ray->dP = sd->dP;
+		ray->dD = differential3_zero();
+		ray->time = sd->time;
+
+		/* Intersect with the same object. if multiple intersections are found it
+		 * will use at most LOCAL_MAX_HITS hits, a random subset of all hits. */
+		scene_intersect_local(kg,
+		                      *ray,
+		                      &isect,
+		                      sd->object,
+		                      &lcg_state,
+		                      LOCAL_MAX_HITS);
+
+		int num_eval_hits = min(isect.num_hits, LOCAL_MAX_HITS);
+
+		for(int hit = 0; hit < num_eval_hits; hit++) {
+			/* Quickly retrieve P and Ng without setting up ShaderData. */
+			float3 hit_P;
+			if(sd->type & PRIMITIVE_TRIANGLE) {
+				hit_P = triangle_refine_local(kg,
+				                              sd,
+				                              &isect.hits[hit],
+				                              ray);
+			}
+#ifdef __OBJECT_MOTION__
+			else  if(sd->type & PRIMITIVE_MOTION_TRIANGLE) {
+				float3 verts[3];
+				motion_triangle_vertices(
+					kg,
+					sd->object,
+					kernel_tex_fetch(__prim_index, isect.hits[hit].prim),
+					sd->time,
+					verts);
+				hit_P = motion_triangle_refine_local(kg,
+				                                     sd,
+				                                     &isect.hits[hit],
+				                                     ray,
+				                                     verts);
+			}
+#endif  /* __OBJECT_MOTION__ */
+
+			float3 hit_Ng = isect.Ng[hit];
+
+			/* Compute smooth normal. */
+			float3 N = hit_Ng;
+			int prim = kernel_tex_fetch(__prim_index, isect.hits[hit].prim);
+			int shader = kernel_tex_fetch(__tri_shader, prim);
+
+			if (shader & SHADER_SMOOTH_NORMAL) {
+				float u = isect.hits[hit].u;
+				float v = isect.hits[hit].v;
+
+				if (sd->type & PRIMITIVE_TRIANGLE) {
+					N = triangle_smooth_normal(kg, N, prim, u, v);
+				}
+#ifdef __OBJECT_MOTION__
+				else if(sd->type & PRIMITIVE_MOTION_TRIANGLE) {
+					N = motion_triangle_smooth_normal(kg, N, sd->object, prim, u, v, sd->time);
+				}
+#endif  /* __OBJECT_MOTION__ */
+			}
+
+			/* Transform normals to world space. */
+			if(isect.hits[hit].object != OBJECT_NONE) {
+				object_normal_transform(kg, sd, &N);
+				object_normal_transform(kg, sd, &hit_Ng);
+			}
+
+			/* Probability densities for local frame axes. */
+			float pdf_N = pick_pdf_N * fabsf(dot(disk_N, hit_Ng));
+			float pdf_T = pick_pdf_T * fabsf(dot(disk_T, hit_Ng));
+			float pdf_B = pick_pdf_B * fabsf(dot(disk_B, hit_Ng));
+
+			/* Multiple importance sample between 3 axes, power heuristic
+			 * found to be slightly better than balance heuristic. */
+			float mis_weight = power_heuristic_3(pdf_N, pdf_T, pdf_B);
+
+			/* Real distance to sampled point. */
+			float r = len(hit_P - sd->P);
+
+			/* Compute weight. */
+			float w = mis_weight / pdf_N;
+			if(isect.num_hits > LOCAL_MAX_HITS) {
+				w *= isect.num_hits/(float)LOCAL_MAX_HITS;
+			}
+
+			float pdf = bssrdf_cubic_pdf(radius, 0.0f, r);
+			float disk_pdf = bssrdf_cubic_pdf(radius, 0.0f, disk_r);
+
+			w *= pdf / disk_pdf;
+
+			/* Sum normal and weight. */
+			sum_N += w * N;
+		}
+	}
+
+	/* Normalize. */
+	float3 N = safe_normalize(sum_N);
+	return is_zero(N) ? sd->N : N;
+}
+
+ccl_device void svm_node_bevel(
+	KernelGlobals *kg,
+	ShaderData *sd,
+	ccl_addr_space PathState *state,
+	float *stack,
+	uint4 node)
+{
+	uint num_samples, radius_offset, normal_offset, out_offset;
+	decode_node_uchar4(node.y, &num_samples, &radius_offset, &normal_offset, &out_offset);
+
+	float radius = stack_load_float(stack, radius_offset);
+	float3 bevel_N = svm_bevel(kg, sd, state, radius, num_samples);
+
+	if(stack_valid(normal_offset)) {
+		/* Preserve input normal. */
+		float3 ref_N = stack_load_float3(stack, normal_offset);
+		bevel_N = normalize(sd->N + (bevel_N - ref_N));
+	}
+
+	stack_store_float3(stack, out_offset, bevel_N);
+}
+
+CCL_NAMESPACE_END
+
diff --git a/intern/cycles/kernel/svm/svm_types.h b/intern/cycles/kernel/svm/svm_types.h
index e81c9a211b0..f08ec76c055 100644
--- a/intern/cycles/kernel/svm/svm_types.h
+++ b/intern/cycles/kernel/svm/svm_types.h
@@ -132,6 +132,7 @@ typedef enum ShaderNodeType {
 	NODE_TEX_VOXEL,
 	NODE_ENTER_BUMP_EVAL,
 	NODE_LEAVE_BUMP_EVAL,
+	NODE_BEVEL,
 } ShaderNodeType;
 
 typedef enum NodeAttributeType {
diff --git a/intern/cycles/render/nodes.cpp b/intern/cycles/render/nodes.cpp
index 2b682756c6a..4452fadaf02 100644
--- a/intern/cycles/render/nodes.cpp
+++ b/intern/cycles/render/nodes.cpp
@@ -5604,4 +5604,44 @@ void TangentNode::compile(OSLCompiler& compiler)
 	compiler.add(this, "node_tangent"); 
 }
 
+/* Bevel */
+
+NODE_DEFINE(BevelNode)
+{
+	NodeType* type = NodeType::add("bevel", create, NodeType::SHADER);
+
+	SOCKET_INT(samples, "Samples", 4);
+
+	SOCKET_IN_FLOAT(radius, "Radius", 0.05f);
+	SOCKET_IN_NORMAL(normal, "Normal", make_float3(0.0f, 0.0f, 0.0f), SocketType::LINK_NORMAL);
+
+	SOCKET_OUT_NORMAL(bevel, "Normal");
+
+	return type;
+}
+
+BevelNode::BevelNode()
+: ShaderNode(node_type)
+{
+}
+
+void BevelNode::compile(SVMCompiler& compiler)
+{
+	ShaderInput *radius_in = input("Radius");
+	ShaderInput *normal_in = input("Normal");
+	ShaderOutput *normal_out = output("Normal");
+
+	compiler.add_node(NODE_BEVEL,
+		compiler.encode_uchar4(samples,
+		                       compiler.stack_assign(radius_in),
+		                       compiler.stack_assign_if_linked(normal_in),
+		                       compiler.stack_assign(normal_out)));
+}
+
+void BevelNode::compile(OSLCompiler& compiler)
+{
+	compiler.parameter(this, "samples");
+	compiler.add(this, "node_bevel");
+}
+
 CCL_NAMESPACE_END
diff --git a/intern/cycles/render/nodes.h b/intern/cycles/render/nodes.h
index 3b8a8bcd4b5..64d7522a23a 100644
--- a/intern/cycles/render/nodes.h
+++ b/intern/cycles/render/nodes.h
@@ -1015,6 +1015,18 @@ public:
 	float3 normal_osl;
 };
 
+class BevelNode : public ShaderNode {
+public:
+	SHADER_NODE_CLASS(BevelNode)
+	bool has_spatial_varying() { return true; }
+	virtual int get_group() { return NODE_GROUP_LEVEL_3; }
+	virtual bool has_raytrace() { return true; }
+
+	float radius;
+	float3 normal;
+	int samples;
+};
+
 CCL_NAMESPACE_END
 
 #endif /* __NODES_H__ */