Cycles: multiple importance sampling for lamps, which helps reduce noise for

big lamps and sharp glossy reflections. This was already supported for mesh
lights and the background, so lamps should do it too.

This is not for free and it's a bit slower than I hoped even though there is
no extra BVH ray intersection. I'll try to optimize it more later.

* Area lights look a bit different now, they had the wrong shape before.
* Also fixes a sampling issue in the non-progressive integrator.
* Only enabled for the CPU, will test on the GPU later.
* An option to disable this will be added for situations where it does not help.

Same time comparison before/after:
http://www.pasteall.org/pic/show.php?id=43313
http://www.pasteall.org/pic/show.php?id=43314

intern/cycles/blender/blender_object.cpp

@@ -127,8 +127,8 @@ void BlenderSync::sync_light(BL::Object b_parent, int persistent_id[OBJECT_PERSI
 		case BL::Lamp::type_AREA: {
 			BL::AreaLamp b_area_lamp(b_lamp);
 			light->size = 1.0f;
-			light->axisu = make_float3(tfm.x.x, tfm.y.x, tfm.z.x);
+			light->axisu = transform_get_column(&tfm, 0);
-			light->axisv = make_float3(tfm.x.y, tfm.y.y, tfm.z.y);
+			light->axisv = transform_get_column(&tfm, 1);
 			light->sizeu = b_area_lamp.size();
 			if(b_area_lamp.shape() == BL::AreaLamp::shape_RECTANGLE)
 				light->sizev = b_area_lamp.size_y();
@@ -140,8 +140,8 @@ void BlenderSync::sync_light(BL::Object b_parent, int persistent_id[OBJECT_PERSI
 	}
 
 	/* location and (inverted!) direction */
-	light->co = make_float3(tfm.x.w, tfm.y.w, tfm.z.w);
+	light->co = transform_get_column(&tfm, 3);
-	light->dir = -make_float3(tfm.x.z, tfm.y.z, tfm.z.z);
+	light->dir = -transform_get_column(&tfm, 2);
 
 	/* shader */
 	vector<uint> used_shaders;

intern/cycles/kernel/kernel_camera.h

@@ -199,7 +199,6 @@ __device void camera_sample_panorama(KernelGlobals *kg, float raster_x, float ra
 
 	Pcamera = transform_perspective(&rastertocamera, make_float3(raster_x, raster_y + 1.0f, 0.0f));
 	ray->dD.dy = normalize(transform_direction(&cameratoworld, panorama_to_direction(kg, Pcamera.x, Pcamera.y))) - ray->D;
-
 #endif
 }
 

intern/cycles/kernel/kernel_emission.h

@@ -66,6 +66,8 @@ __device float3 direct_emissive_eval(KernelGlobals *kg, float rando,
 		else
 			eval = make_float3(0.0f, 0.0f, 0.0f);
 	}
+	
+	eval *= ls->eval_fac;
 
 	shader_release(kg, &sd);
 
@@ -74,29 +76,29 @@ __device float3 direct_emissive_eval(KernelGlobals *kg, float rando,
 
 __device bool direct_emission(KernelGlobals *kg, ShaderData *sd, int lindex,
 	float randt, float rando, float randu, float randv, Ray *ray, BsdfEval *eval,
-	bool *is_lamp)
+	int *lamp)
 {
 	LightSample ls;
 
-	float pdf = -1.0f;
-
 #ifdef __NON_PROGRESSIVE__
 	if(lindex != -1) {
 		/* sample position on a specified light */
-		light_select(kg, lindex, randu, randv, sd->P, &ls, &pdf);
+		light_select(kg, lindex, randu, randv, sd->P, &ls);
 	}
 	else
 #endif
 	{
 		/* sample a light and position on int */
-		light_sample(kg, randt, randu, randv, sd->time, sd->P, &ls, &pdf);
+		light_sample(kg, randt, randu, randv, sd->time, sd->P, &ls);
 	}
 
-	/* compute pdf */
+	/* return lamp index for MIS */
-	if(pdf < 0.0f)
+	if(ls.use_mis)
-		pdf = light_sample_pdf(kg, &ls, -ls.D, ls.t);
+		*lamp = ls.lamp;
+	else
+		*lamp= ~0;
 
-	if(pdf == 0.0f)
+	if(ls.pdf == 0.0f)
 		return false;
 
 	/* evaluate closure */
@@ -112,13 +114,13 @@ __device bool direct_emission(KernelGlobals *kg, ShaderData *sd, int lindex,
 
 	shader_bsdf_eval(kg, sd, ls.D, eval, &bsdf_pdf);
 
-	if(ls.prim != ~0 || ls.type == LIGHT_BACKGROUND) {
+	if(ls.use_mis) {
 		/* multiple importance sampling */
-		float mis_weight = power_heuristic(pdf, bsdf_pdf);
+		float mis_weight = power_heuristic(ls.pdf, bsdf_pdf);
 		light_eval *= mis_weight;
 	}
 	
-	bsdf_eval_mul(eval, light_eval*(ls.eval_fac/pdf));
+	bsdf_eval_mul(eval, light_eval/ls.pdf);
 
 	if(bsdf_eval_is_zero(eval))
 		return false;
@@ -144,14 +146,12 @@ __device bool direct_emission(KernelGlobals *kg, ShaderData *sd, int lindex,
 		ray->t = 0.0f;
 	}
 
-	*is_lamp = (ls.prim == ~0);
-
 	return true;
 }
 
-/* Indirect Emission */
+/* Indirect Primitive Emission */
 
-__device float3 indirect_emission(KernelGlobals *kg, ShaderData *sd, float t, int path_flag, float bsdf_pdf)
+__device float3 indirect_primitive_emission(KernelGlobals *kg, ShaderData *sd, float t, int path_flag, float bsdf_pdf)
 {
 	/* evaluate emissive closure */
 	float3 L = shader_emissive_eval(kg, sd);
@@ -172,6 +172,35 @@ __device float3 indirect_emission(KernelGlobals *kg, ShaderData *sd, float t, in
 	return L;
 }
 
+/* Indirect Lamp Emission */
+
+__device bool indirect_lamp_emission(KernelGlobals *kg, Ray *ray, int path_flag, float bsdf_pdf, float randt, float3 *emission)
+{
+	LightSample ls;
+	int lamp = lamp_light_eval_sample(kg, randt);
+
+	if(lamp == ~0)
+		return false;
+
+	if(!lamp_light_eval(kg, lamp, ray->P, ray->D, ray->t, &ls))
+		return false;
+	
+	/* todo: missing texture coordinates */
+	float u = 0.0f;
+	float v = 0.0f;
+	float3 L = direct_emissive_eval(kg, 0.0f, &ls, u, v, -ray->D, ls.t, ray->time);
+
+	if(!(path_flag & PATH_RAY_MIS_SKIP)) {
+		/* multiple importance sampling, get regular light pdf,
+		 * and compute weight with respect to BSDF pdf */
+		float mis_weight = power_heuristic(bsdf_pdf, ls.pdf);
+		L *= mis_weight;
+	}
+
+	*emission = L;
+	return true;
+}
+
 /* Indirect Background */
 
 __device float3 indirect_background(KernelGlobals *kg, Ray *ray, int path_flag, float bsdf_pdf)

intern/cycles/kernel/kernel_light.h

@@ -18,49 +18,27 @@
 
 CCL_NAMESPACE_BEGIN
 
+/* Light Sample result */
+
 typedef struct LightSample {
-	float3 P;
+	float3 P;			/* position on light, or direction for distant light */
-	float3 D;
+	float3 Ng;			/* normal on light */
-	float3 Ng;
+	float3 D;			/* direction from shading point to light */
-	float t;
+	float t;			/* distance to light (FLT_MAX for distant light) */
-	float eval_fac;
+	float pdf;			/* light sampling probability density function */
-	int object;
+	float eval_fac;		/* intensity multiplier */
-	int prim;
+	int object;			/* object id for triangle/curve lights */
-	int shader;
+	int prim;			/* primitive id for triangle/curve ligths */
-	LightType type;
+	int shader;			/* shader id */
+	int lamp;			/* lamp id */
+	int use_mis;		/* for lamps with size zero */
+	LightType type;		/* type of light */
 } LightSample;
 
-/* Regular Light */
+/* Background Light */
-
-__device float3 disk_light_sample(float3 v, float randu, float randv)
-{
-	float3 ru, rv;
-
-	make_orthonormals(v, &ru, &rv);
-	to_unit_disk(&randu, &randv);
-
-	return ru*randu + rv*randv;
-}
-
-__device float3 distant_light_sample(float3 D, float size, float randu, float randv)
-{
-	return normalize(D + disk_light_sample(D, randu, randv)*size);
-}
-
-__device float3 sphere_light_sample(float3 P, float3 center, float size, float randu, float randv)
-{
-	return disk_light_sample(normalize(P - center), randu, randv)*size;
-}
-
-__device float3 area_light_sample(float3 axisu, float3 axisv, float randu, float randv)
-{
-	randu = randu - 0.5f;
-	randv = randv - 0.5f;
-
-	return axisu*randu + axisv*randv;
-}
 
 #ifdef __BACKGROUND_MIS__
+
 __device float3 background_light_sample(KernelGlobals *kg, float randu, float randv, float *pdf)
 {
 	/* for the following, the CDF values are actually a pair of floats, with the
@@ -169,33 +147,108 @@ __device float background_light_pdf(KernelGlobals *kg, float3 direction)
 }
 #endif
 
-__device void regular_light_sample(KernelGlobals *kg, int point,
+/* Regular Light */
-	float randu, float randv, float3 P, LightSample *ls, float *pdf)
+
+__device float3 disk_light_sample(float3 v, float randu, float randv)
 {
-	float4 data0 = kernel_tex_fetch(__light_data, point*LIGHT_SIZE + 0);
+	float3 ru, rv;
-	float4 data1 = kernel_tex_fetch(__light_data, point*LIGHT_SIZE + 1);
+
+	make_orthonormals(v, &ru, &rv);
+	to_unit_disk(&randu, &randv);
+
+	return ru*randu + rv*randv;
+}
+
+__device float3 distant_light_sample(float3 D, float radius, float randu, float randv)
+{
+	return normalize(D + disk_light_sample(D, randu, randv)*radius);
+}
+
+__device float3 sphere_light_sample(float3 P, float3 center, float radius, float randu, float randv)
+{
+	return disk_light_sample(normalize(P - center), randu, randv)*radius;
+}
+
+__device float3 area_light_sample(float3 axisu, float3 axisv, float randu, float randv)
+{
+	randu = randu - 0.5f;
+	randv = randv - 0.5f;
+
+	return axisu*randu + axisv*randv;
+}
+
+__device float spot_light_attenuation(float4 data1, float4 data2, LightSample *ls)
+{
+	float3 dir = make_float3(data2.y, data2.z, data2.w);
+	float3 I = ls->Ng;
+
+	float spot_angle = data1.w;
+	float spot_smooth = data2.x;
+
+	float attenuation = dot(dir, I);
+
+	if(attenuation <= spot_angle) {
+		attenuation = 0.0f;
+	}
+	else {
+		float t = attenuation - spot_angle;
+
+		if(t < spot_smooth && spot_smooth != 0.0f)
+			attenuation *= smoothstepf(t/spot_smooth);
+	}
+
+	return attenuation;
+}
+
+__device float lamp_light_pdf(KernelGlobals *kg, const float3 Ng, const float3 I, float t)
+{
+	float cos_pi = dot(Ng, I);
+
+	if(cos_pi <= 0.0f)
+		return 0.0f;
+	
+	return t*t/cos_pi;
+}
+
+__device void lamp_light_sample(KernelGlobals *kg, int lamp,
+	float randu, float randv, float3 P, LightSample *ls)
+{
+	float4 data0 = kernel_tex_fetch(__light_data, lamp*LIGHT_SIZE + 0);
+	float4 data1 = kernel_tex_fetch(__light_data, lamp*LIGHT_SIZE + 1);
 
 	LightType type = (LightType)__float_as_int(data0.x);
 	ls->type = type;
+#ifdef __LAMP_MIS__
+	ls->use_mis = true;
+#else
+	ls->use_mis = false;
+#endif
 
 	if(type == LIGHT_DISTANT) {
 		/* distant light */
-		float3 D = make_float3(data0.y, data0.z, data0.w);
+		float3 lightD = make_float3(data0.y, data0.z, data0.w);
-		float size = data1.y;
+		float3 D = lightD;
+		float radius = data1.y;
+		float invarea = data1.w;
 
-		if(size > 0.0f)
+		if(radius > 0.0f)
-			D = distant_light_sample(D, size, randu, randv);
+			D = distant_light_sample(D, radius, randu, randv);
+		else
+			ls->use_mis = false;
 
 		ls->P = D;
 		ls->Ng = D;
 		ls->D = -D;
 		ls->t = FLT_MAX;
-		ls->eval_fac = 1.0f;
+
+		float costheta = dot(lightD, D);
+		ls->pdf = invarea/(costheta*costheta*costheta);
+		ls->eval_fac = ls->pdf;
 	}
 #ifdef __BACKGROUND_MIS__
 	else if(type == LIGHT_BACKGROUND) {
 		/* infinite area light (e.g. light dome or env light) */
-		float3 D = background_light_sample(kg, randu, randv, pdf);
+		float3 D = background_light_sample(kg, randu, randv, &ls->pdf);
 
 		ls->P = D;
 		ls->Ng = D;
@@ -207,85 +260,208 @@ __device void regular_light_sample(KernelGlobals *kg, int point,
 	else {
 		ls->P = make_float3(data0.y, data0.z, data0.w);
 
-		if(type == LIGHT_POINT) {
+		if(type == LIGHT_POINT || type == LIGHT_SPOT) {
-			float size = data1.y;
+			float radius = data1.y;
 
-			/* sphere light */
+			if(radius > 0.0f)
-			if(size > 0.0f)
+				/* sphere light */
-				ls->P += sphere_light_sample(P, ls->P, size, randu, randv);
+				ls->P += sphere_light_sample(P, ls->P, radius, randu, randv);
+			else
+				ls->use_mis = false;
 
-			ls->Ng = normalize(P - ls->P);
+			ls->D = normalize_len(ls->P - P, &ls->t);
-			ls->eval_fac = 0.25f*M_1_PI_F;
+			ls->Ng = -ls->D;
-		}
-		else if(type == LIGHT_SPOT) {
-			float4 data2 = kernel_tex_fetch(__light_data, point*LIGHT_SIZE + 2);
-			float size = data1.y;
 
-			/* spot light */
+			float invarea = data1.z;
-			if(size > 0.0f)
+			ls->eval_fac = (0.25f*M_1_PI_F)*invarea;
-				ls->P += sphere_light_sample(P, ls->P, size, randu, randv);
+			ls->pdf = invarea;
 
-			float3 dir = make_float3(data1.z, data1.w, data2.x);
+			if(type == LIGHT_SPOT) {
-			float3 I = normalize(P - ls->P);
+				/* spot light attentuation */
-
+				float4 data2 = kernel_tex_fetch(__light_data, lamp*LIGHT_SIZE + 2);
-			float spot_angle = data2.y;
+				ls->eval_fac *= spot_light_attenuation(data1, data2, ls);
-			float spot_smooth = data2.z;
-
-			float eval_fac = dot(dir, I);
-
-			if(eval_fac <= spot_angle) {
-				eval_fac = 0.0f;
 			}
-			else {
-				float t = eval_fac - spot_angle;
-
-				if(t < spot_smooth && spot_smooth != 0.0f)
-					eval_fac *= smoothstepf(t/spot_smooth);
-			}
-
-			ls->Ng = I;
-			ls->eval_fac = eval_fac*0.25f*M_1_PI_F;
 		}
 		else {
 			/* area light */
-			float4 data2 = kernel_tex_fetch(__light_data, point*LIGHT_SIZE + 2);
+			float4 data2 = kernel_tex_fetch(__light_data, lamp*LIGHT_SIZE + 2);
-			float4 data3 = kernel_tex_fetch(__light_data, point*LIGHT_SIZE + 3);
+			float4 data3 = kernel_tex_fetch(__light_data, lamp*LIGHT_SIZE + 3);
 
-			float3 axisu = make_float3(data1.y, data1.z, data2.w);
+			float3 axisu = make_float3(data1.y, data1.z, data1.w);
 			float3 axisv = make_float3(data2.y, data2.z, data2.w);
 			float3 D = make_float3(data3.y, data3.z, data3.w);
 
 			ls->P += area_light_sample(axisu, axisv, randu, randv);
 			ls->Ng = D;
-			ls->eval_fac = 0.25f;
+			ls->D = normalize_len(ls->P - P, &ls->t);
-		}
 
-		ls->t = 0.0f;
+			float invarea = data2.x;
+
+			if(invarea == 0.0f) {
+				ls->use_mis = false;
+				invarea = 1.0f;
+			}
+
+			ls->pdf = invarea;
+			ls->eval_fac = 0.25f*ls->pdf;
+		}
 	}
 
 	ls->shader = __float_as_int(data1.x);
 	ls->object = ~0;
 	ls->prim = ~0;
+	ls->lamp = lamp;
+
+	/* compute pdf */
+	if(ls->t != FLT_MAX)
+		ls->pdf *= lamp_light_pdf(kg, ls->Ng, -ls->D, ls->t);
+	
+	/* this is a bit weak, but we don't want this as part of the pdf for
+	 * multiple importance sampling */
+	ls->eval_fac *= kernel_data.integrator.inv_pdf_lights;
 }
 
-__device float regular_light_pdf(KernelGlobals *kg,
+__device bool lamp_light_eval(KernelGlobals *kg, int lamp, float3 P, float3 D, float t, LightSample *ls)
-	const float3 Ng, const float3 I, float t)
 {
-	float pdf = kernel_data.integrator.pdf_lights;
+	float4 data0 = kernel_tex_fetch(__light_data, lamp*LIGHT_SIZE + 0);
+	float4 data1 = kernel_tex_fetch(__light_data, lamp*LIGHT_SIZE + 1);
 
-	if(t == FLT_MAX)
+	LightType type = (LightType)__float_as_int(data0.x);
-		return pdf;
+	ls->type = type;
+	ls->shader = __float_as_int(data1.x);
+	ls->object = ~0;
+	ls->prim = ~0;
+	ls->lamp = lamp;
+	ls->use_mis = false; /* flag not used for eval */
 
-	float cos_pi = dot(Ng, I);
+	if(type == LIGHT_DISTANT) {
+		/* distant light */
+		float radius = data1.y;
 
-	if(cos_pi <= 0.0f)
+		if(radius == 0.0f)
-		return 0.0f;
+			return false;
+		if(t != FLT_MAX)
+			return false;
 
-	return t*t*pdf/cos_pi;
+		/* a distant light is infinitely far away, but equivalent to a disk
+		 * shaped light exactly 1 unit away from the current shading point.
+		 *
+		 *     radius              t^2/cos(theta)
+		 *  <---------->           t = sqrt(1^2 + tan(theta)^2)
+		 *       tan(th)           area = radius*radius*pi
+		 *       <----->
+		 *        \    |           (1 + tan(theta)^2)/cos(theta)
+		 *         \   |           (1 + tan(acos(cos(theta)))^2)/cos(theta)
+		 *       t  \th| 1         simplifies to
+		 *           \-|           1/(cos(theta)^3)
+		 *            \|           magic!
+		 *             P
+		 */
+
+		float3 lightD = make_float3(data0.y, data0.z, data0.w);
+		float costheta = dot(-lightD, D);
+		float cosangle = data1.z;
+
+		if(costheta < cosangle)
+			return false;
+
+		ls->P = -D;
+		ls->Ng = -D;
+		ls->D = D;
+		ls->t = FLT_MAX;
+
+		float invarea = data1.w;
+		ls->pdf = invarea/(costheta*costheta*costheta);
+		ls->eval_fac = ls->pdf;
+	}
+	else if(type == LIGHT_POINT || type == LIGHT_SPOT) {
+		float3 lightP = make_float3(data0.y, data0.z, data0.w);
+		float radius = data1.y;
+
+		/* sphere light */
+		if(radius == 0.0f)
+			return false;
+
+		if(!ray_aligned_disk_intersect(P, D, t,
+			lightP, radius, &ls->P, &ls->t))
+			return false;
+
+		ls->Ng = -D;
+		ls->D = D;
+
+		float invarea = data1.z;
+		ls->eval_fac = (0.25f*M_1_PI_F)*invarea;
+		ls->pdf = invarea;
+
+		if(type == LIGHT_SPOT) {
+			/* spot light attentuation */
+			float4 data2 = kernel_tex_fetch(__light_data, lamp*LIGHT_SIZE + 2);
+			ls->eval_fac *= spot_light_attenuation(data1, data2, ls);
+
+			if(ls->eval_fac == 0.0f)
+				return false;
+		}
+	}
+	else if(type == LIGHT_AREA) {
+		/* area light */
+		float4 data2 = kernel_tex_fetch(__light_data, lamp*LIGHT_SIZE + 2);
+		float4 data3 = kernel_tex_fetch(__light_data, lamp*LIGHT_SIZE + 3);
+
+		float invarea = data2.x;
+		if(invarea == 0.0f)
+			return false;
+
+		float3 axisu = make_float3(data1.y, data1.z, data1.w);
+		float3 axisv = make_float3(data2.y, data2.z, data2.w);
+		float3 Ng = make_float3(data3.y, data3.z, data3.w);
+
+		/* one sided */
+		if(dot(D, Ng) >= 0.0f)
+			return false;
+
+		ls->P = make_float3(data0.y, data0.z, data0.w);
+
+		if(!ray_quad_intersect(P, D, t,
+			ls->P, axisu, axisv, &ls->P, &ls->t))
+			return false;
+
+		ls->D = D;
+		ls->Ng = Ng;
+		ls->pdf = invarea;
+		ls->eval_fac = 0.25f*ls->pdf;
+	}
+	else
+		return false;
+
+	/* compute pdf */
+	if(ls->t != FLT_MAX)
+		ls->pdf *= lamp_light_pdf(kg, ls->Ng, -ls->D, ls->t);
+	ls->eval_fac *= kernel_data.integrator.inv_pdf_lights;
+
+	return true;
 }
 
 /* Triangle Light */
 
+__device void object_transform_light_sample(KernelGlobals *kg, LightSample *ls, int object, float time)
+{
+#ifdef __INSTANCING__
+	/* instance transform */
+	if(object >= 0) {
+#ifdef __OBJECT_MOTION__
+		Transform itfm;
+		Transform tfm = object_fetch_transform_motion_test(kg, object, time, &itfm);
+#else
+		Transform tfm = object_fetch_transform(kg, object, OBJECT_TRANSFORM);
+		Transform itfm = object_fetch_transform(kg, object, OBJECT_INVERSE_TRANSFORM);
+#endif
+
+		ls->P = transform_point(&tfm, ls->P);
+		ls->Ng = normalize(transform_direction(&tfm, ls->Ng));
+	}
+#endif
+}
+
 __device void triangle_light_sample(KernelGlobals *kg, int prim, int object,
 	float randu, float randv, float time, LightSample *ls)
 {
@@ -294,40 +470,30 @@ __device void triangle_light_sample(KernelGlobals *kg, int prim, int object,
 	ls->Ng = triangle_normal_MT(kg, prim, &ls->shader);
 	ls->object = object;
 	ls->prim = prim;
+	ls->lamp = ~0;
+	ls->use_mis = true;
 	ls->t = 0.0f;
 	ls->type = LIGHT_AREA;
 	ls->eval_fac = 1.0f;
 
-#ifdef __INSTANCING__
+	object_transform_light_sample(kg, ls, object, time);
-	/* instance transform */
-	if(ls->object >= 0) {
-#ifdef __OBJECT_MOTION__
-		Transform itfm;
-		Transform tfm = object_fetch_transform_motion_test(kg, object, time, &itfm);
-#else
-		Transform tfm = object_fetch_transform(kg, ls->object, OBJECT_TRANSFORM);
-		Transform itfm = object_fetch_transform(kg, ls->object, OBJECT_INVERSE_TRANSFORM);
-#endif
-
-		ls->P = transform_point(&tfm, ls->P);
-		ls->Ng = normalize(transform_direction_transposed(&itfm, ls->Ng));
-	}
-#endif
 }
 
 __device float triangle_light_pdf(KernelGlobals *kg,
 	const float3 Ng, const float3 I, float t)
 {
+	float pdf = kernel_data.integrator.pdf_triangles;
 	float cos_pi = fabsf(dot(Ng, I));
 
 	if(cos_pi == 0.0f)
 		return 0.0f;
 	
-	return (t*t*kernel_data.integrator.pdf_triangles)/cos_pi;
+	return t*t*pdf/cos_pi;
 }
 
+/* Curve Light */
+
 #ifdef __HAIR__
-/* Strand Light */
 
 __device void curve_segment_light_sample(KernelGlobals *kg, int prim, int object,
 	int segment, float randu, float randv, float time, LightSample *ls)
@@ -358,27 +524,16 @@ __device void curve_segment_light_sample(KernelGlobals *kg, int prim, int object
 	ls->P = randu * l * tg + (gd * l + r1) * ls->Ng;
 	ls->object = object;
 	ls->prim = prim;
+	ls->lamp = ~0;
+	ls->use_mis = true;
 	ls->t = 0.0f;
 	ls->type = LIGHT_STRAND;
 	ls->eval_fac = 1.0f;
 	ls->shader = __float_as_int(v00.z);
 
-#ifdef __INSTANCING__
+	object_transform_light_sample(kg, ls, object, time);
-	/* instance transform */
-	if(ls->object >= 0) {
-#ifdef __OBJECT_MOTION__
-		Transform itfm;
-		Transform tfm = object_fetch_transform_motion_test(kg, object, time, &itfm);
-#else
-		Transform tfm = object_fetch_transform(kg, ls->object, OBJECT_TRANSFORM);
-		Transform itfm = object_fetch_transform(kg, ls->object, OBJECT_INVERSE_TRANSFORM);
-#endif
-
-		ls->P = transform_point(&tfm, ls->P);
-		ls->Ng = normalize(transform_direction(&tfm, ls->Ng));
-	}
-#endif
 }
+
 #endif
 
 /* Light Distribution */
@@ -412,7 +567,7 @@ __device int light_distribution_sample(KernelGlobals *kg, float randt)
 
 /* Generic Light */
 
-__device void light_sample(KernelGlobals *kg, float randt, float randu, float randv, float time, float3 P, LightSample *ls, float *pdf)
+__device void light_sample(KernelGlobals *kg, float randt, float randu, float randv, float time, float3 P, LightSample *ls)
 {
 	/* sample index */
 	int index = light_distribution_sample(kg, randt);
@@ -420,12 +575,12 @@ __device void light_sample(KernelGlobals *kg, float randt, float randu, float ra
 	/* fetch light data */
 	float4 l = kernel_tex_fetch(__light_distribution, index);
 	int prim = __float_as_int(l.y);
-#ifdef __HAIR__
-	int segment = __float_as_int(l.z);
-#endif
 
 	if(prim >= 0) {
 		int object = __float_as_int(l.w);
+#ifdef __HAIR__
+		int segment = __float_as_int(l.z);
+#endif
 
 #ifdef __HAIR__
 		if (segment != ~0)
@@ -433,27 +588,15 @@ __device void light_sample(KernelGlobals *kg, float randt, float randu, float ra
 		else
 #endif
 			triangle_light_sample(kg, prim, object, randu, randv, time, ls);
+
+		/* compute incoming direction, distance and pdf */
+		ls->D = normalize_len(ls->P - P, &ls->t);
+		ls->pdf = triangle_light_pdf(kg, ls->Ng, -ls->D, ls->t);
 	}
 	else {
-		int point = -prim-1;
+		int lamp = -prim-1;
-		regular_light_sample(kg, point, randu, randv, P, ls, pdf);
+		lamp_light_sample(kg, lamp, randu, randv, P, ls);
 	}
-
-	/* compute incoming direction and distance */
-	if(ls->t != FLT_MAX)
-		ls->D = normalize_len(ls->P - P, &ls->t);
-}
-
-__device float light_sample_pdf(KernelGlobals *kg, LightSample *ls, float3 I, float t)
-{
-	float pdf;
-
-	if(ls->prim != ~0)
-		pdf = triangle_light_pdf(kg, ls->Ng, I, t);
-	else
-		pdf = regular_light_pdf(kg, ls->Ng, I, t);
-	
-	return pdf;
 }
 
 __device int light_select_num_samples(KernelGlobals *kg, int index)
@@ -462,18 +605,26 @@ __device int light_select_num_samples(KernelGlobals *kg, int index)
 	return __float_as_int(data3.x);
 }
 
-__device void light_select(KernelGlobals *kg, int index, float randu, float randv, float3 P, LightSample *ls, float *pdf)
+__device void light_select(KernelGlobals *kg, int index, float randu, float randv, float3 P, LightSample *ls)
 {
-	regular_light_sample(kg, index, randu, randv, P, ls, pdf);
+	lamp_light_sample(kg, index, randu, randv, P, ls);
-
-	/* compute incoming direction and distance */
-	if(ls->t != FLT_MAX)
-		ls->D = normalize_len(ls->P - P, &ls->t);
 }
 
-__device float light_select_pdf(KernelGlobals *kg, LightSample *ls, float3 I, float t)
+__device int lamp_light_eval_sample(KernelGlobals *kg, float randt)
 {
-	return regular_light_pdf(kg, ls->Ng, I, t);
+	/* sample index */
+	int index = light_distribution_sample(kg, randt);
+
+	/* fetch light data */
+	float4 l = kernel_tex_fetch(__light_distribution, index);
+	int prim = __float_as_int(l.y);
+
+	if(prim < 0) {
+		int lamp = -prim-1;
+		return lamp;
+	}
+	else
+		return ~0;
 }
 
 CCL_NAMESPACE_END

intern/cycles/kernel/kernel_path.h

@@ -238,6 +238,7 @@ __device float4 kernel_path_progressive(KernelGlobals *kg, RNG *rng, int sample,
 
 	float min_ray_pdf = FLT_MAX;
 	float ray_pdf = 0.0f;
+	float ray_t = 0.0f;
 	PathState state;
 	int rng_offset = PRNG_BASE_NUM;
 
@@ -248,8 +249,29 @@ __device float4 kernel_path_progressive(KernelGlobals *kg, RNG *rng, int sample,
 		/* intersect scene */
 		Intersection isect;
 		uint visibility = path_state_ray_visibility(kg, &state);
+		bool hit = scene_intersect(kg, &ray, visibility, &isect);
 
-		if(!scene_intersect(kg, &ray, visibility, &isect)) {
+#ifdef __LAMP_MIS__
+		if(kernel_data.integrator.pdf_lights > 0.0f && !(state.flag & PATH_RAY_CAMERA)) {
+			/* ray starting from previous non-transparent bounce */
+			Ray light_ray;
+
+			light_ray.P = ray.P - ray_t*ray.D;
+			ray_t += isect.t;
+			light_ray.D = ray.D;
+			light_ray.t = ray_t;
+			light_ray.time = ray.time;
+
+			/* intersect with lamp */
+			float light_t = path_rng(kg, rng, sample, rng_offset + PRNG_LIGHT);
+			float3 emission;
+
+			if(indirect_lamp_emission(kg, &light_ray, state.flag, ray_pdf, light_t, &emission))
+				path_radiance_accum_emission(&L, throughput, emission, state.bounce);
+		}
+#endif
+
+		if(!hit) {
 			/* eval background shader if nothing hit */
 			if(kernel_data.background.transparent && (state.flag & PATH_RAY_CAMERA)) {
 				L_transparent += average(throughput);
@@ -313,7 +335,8 @@ __device float4 kernel_path_progressive(KernelGlobals *kg, RNG *rng, int sample,
 #ifdef __EMISSION__
 		/* emission */
 		if(sd.flag & SD_EMISSION) {
-			float3 emission = indirect_emission(kg, &sd, isect.t, state.flag, ray_pdf);
+			/* todo: is isect.t wrong here for transparent surfaces? */
+			float3 emission = indirect_primitive_emission(kg, &sd, isect.t, state.flag, ray_pdf);
 			path_radiance_accum_emission(&L, throughput, emission, state.bounce);
 		}
 #endif
@@ -374,18 +397,19 @@ __device float4 kernel_path_progressive(KernelGlobals *kg, RNG *rng, int sample,
 
 				Ray light_ray;
 				BsdfEval L_light;
-				bool is_lamp;
+				int lamp;
 
 #ifdef __OBJECT_MOTION__
 				light_ray.time = sd.time;
 #endif
 
-				if(direct_emission(kg, &sd, -1, light_t, light_o, light_u, light_v, &light_ray, &L_light, &is_lamp)) {
+				if(direct_emission(kg, &sd, -1, light_t, light_o, light_u, light_v, &light_ray, &L_light, &lamp)) {
 					/* trace shadow ray */
 					float3 shadow;
 
 					if(!shadow_blocked(kg, &state, &light_ray, &shadow)) {
 						/* accumulate */
+						bool is_lamp = (lamp != ~0);
 						path_radiance_accum_light(&L, throughput, &L_light, shadow, state.bounce, is_lamp);
 					}
 				}
@@ -422,6 +446,7 @@ __device float4 kernel_path_progressive(KernelGlobals *kg, RNG *rng, int sample,
 		/* set labels */
 		if(!(label & LABEL_TRANSPARENT)) {
 			ray_pdf = bsdf_pdf;
+			ray_t = 0.0f;
 			min_ray_pdf = fminf(bsdf_pdf, min_ray_pdf);
 		}
 
@@ -459,13 +484,36 @@ __device float4 kernel_path_progressive(KernelGlobals *kg, RNG *rng, int sample,
 __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray ray, __global float *buffer,
 	float3 throughput, float min_ray_pdf, float ray_pdf, PathState state, int rng_offset, PathRadiance *L)
 {
+	float ray_t = 0.0f;
+
 	/* path iteration */
 	for(;; rng_offset += PRNG_BOUNCE_NUM) {
 		/* intersect scene */
 		Intersection isect;
 		uint visibility = path_state_ray_visibility(kg, &state);
+		bool hit = scene_intersect(kg, &ray, visibility, &isect);
 
-		if(!scene_intersect(kg, &ray, visibility, &isect)) {
+#ifdef __LAMP_MIS__
+		if(kernel_data.integrator.pdf_lights > 0.0f && !(state.flag & PATH_RAY_CAMERA)) {
+			/* ray starting from previous non-transparent bounce */
+			Ray light_ray;
+
+			light_ray.P = ray.P - ray_t*ray.D;
+			ray_t += isect.t;
+			light_ray.D = ray.D;
+			light_ray.t = ray_t;
+			light_ray.time = ray.time;
+
+			/* intersect with lamp */
+			float light_t = path_rng(kg, rng, sample, rng_offset + PRNG_LIGHT);
+			float3 emission;
+
+			if(indirect_lamp_emission(kg, &light_ray, state.flag, ray_pdf, light_t, &emission))
+				path_radiance_accum_emission(L, throughput, emission, state.bounce);
+		}
+#endif
+
+		if(!hit) {
 #ifdef __BACKGROUND__
 			/* sample background shader */
 			float3 L_background = indirect_background(kg, &ray, state.flag, ray_pdf);
@@ -496,7 +544,7 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 #ifdef __EMISSION__
 		/* emission */
 		if(sd.flag & SD_EMISSION) {
-			float3 emission = indirect_emission(kg, &sd, isect.t, state.flag, ray_pdf);
+			float3 emission = indirect_primitive_emission(kg, &sd, isect.t, state.flag, ray_pdf);
 			path_radiance_accum_emission(L, throughput, emission, state.bounce);
 		}
 #endif
@@ -557,19 +605,20 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 
 				Ray light_ray;
 				BsdfEval L_light;
-				bool is_lamp;
+				int lamp;
 
 #ifdef __OBJECT_MOTION__
 				light_ray.time = sd.time;
 #endif
 
 				/* sample random light */
-				if(direct_emission(kg, &sd, -1, light_t, light_o, light_u, light_v, &light_ray, &L_light, &is_lamp)) {
+				if(direct_emission(kg, &sd, -1, light_t, light_o, light_u, light_v, &light_ray, &L_light, &lamp)) {
 					/* trace shadow ray */
 					float3 shadow;
 
 					if(!shadow_blocked(kg, &state, &light_ray, &shadow)) {
 						/* accumulate */
+						bool is_lamp = (lamp != ~0);
 						path_radiance_accum_light(L, throughput, &L_light, shadow, state.bounce, is_lamp);
 					}
 				}
@@ -606,6 +655,7 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 		/* set labels */
 		if(!(label & LABEL_TRANSPARENT)) {
 			ray_pdf = bsdf_pdf;
+			ray_t = 0.0f;
 			min_ray_pdf = fminf(bsdf_pdf, min_ray_pdf);
 		}
 
@@ -697,7 +747,7 @@ __device float4 kernel_path_non_progressive(KernelGlobals *kg, RNG *rng, int sam
 #ifdef __EMISSION__
 		/* emission */
 		if(sd.flag & SD_EMISSION) {
-			float3 emission = indirect_emission(kg, &sd, isect.t, state.flag, ray_pdf);
+			float3 emission = indirect_primitive_emission(kg, &sd, isect.t, state.flag, ray_pdf);
 			path_radiance_accum_emission(&L, throughput, emission, state.bounce);
 		}
 #endif
@@ -760,7 +810,7 @@ __device float4 kernel_path_non_progressive(KernelGlobals *kg, RNG *rng, int sam
 		if(sd.flag & SD_BSDF_HAS_EVAL) {
 			Ray light_ray;
 			BsdfEval L_light;
-			bool is_lamp;
+			int lamp;
 
 #ifdef __OBJECT_MOTION__
 			light_ray.time = sd.time;
@@ -778,12 +828,13 @@ __device float4 kernel_path_non_progressive(KernelGlobals *kg, RNG *rng, int sam
 					float light_u = path_rng(kg, rng, sample*num_samples + j, rng_offset + PRNG_LIGHT_U);
 					float light_v = path_rng(kg, rng, sample*num_samples + j, rng_offset + PRNG_LIGHT_V);
 
-					if(direct_emission(kg, &sd, i, 0.0f, 0.0f, light_u, light_v, &light_ray, &L_light, &is_lamp)) {
+					if(direct_emission(kg, &sd, i, 0.0f, 0.0f, light_u, light_v, &light_ray, &L_light, &lamp)) {
 						/* trace shadow ray */
 						float3 shadow;
 
 						if(!shadow_blocked(kg, &state, &light_ray, &shadow)) {
 							/* accumulate */
+							bool is_lamp = (lamp != ~0);
 							path_radiance_accum_light(&L, throughput*num_samples_inv, &L_light, shadow, state.bounce, is_lamp);
 						}
 					}
@@ -807,12 +858,13 @@ __device float4 kernel_path_non_progressive(KernelGlobals *kg, RNG *rng, int sam
 					if(kernel_data.integrator.num_all_lights)
 						light_t = 0.5f*light_t;
 
-					if(direct_emission(kg, &sd, -1, light_t, 0.0f, light_u, light_v, &light_ray, &L_light, &is_lamp)) {
+					if(direct_emission(kg, &sd, -1, light_t, 0.0f, light_u, light_v, &light_ray, &L_light, &lamp)) {
 						/* trace shadow ray */
 						float3 shadow;
 
 						if(!shadow_blocked(kg, &state, &light_ray, &shadow)) {
 							/* accumulate */
+							bool is_lamp = (lamp != ~0);
 							path_radiance_accum_light(&L, throughput*num_samples_inv, &L_light, shadow, state.bounce, is_lamp);
 						}
 					}
@@ -885,7 +937,7 @@ __device float4 kernel_path_non_progressive(KernelGlobals *kg, RNG *rng, int sam
 				bsdf_ray.time = sd.time;
 #endif
 
-				kernel_path_indirect(kg, rng, sample*num_samples, bsdf_ray, buffer,
+				kernel_path_indirect(kg, rng, sample*num_samples + j, bsdf_ray, buffer,
 					tp*num_samples_inv, min_ray_pdf, bsdf_pdf, ps, rng_offset+PRNG_BOUNCE_NUM, &L);
 			}
 		}

intern/cycles/kernel/kernel_types.h

@@ -48,6 +48,7 @@ CCL_NAMESPACE_BEGIN
 #endif
 #define __NON_PROGRESSIVE__
 #define __HAIR__
+#define __LAMP_MIS__
 #endif
 
 #ifdef __KERNEL_CUDA__
@@ -384,7 +385,7 @@ typedef enum AttributeStandard {
 
 /* Closure data */
 
-#define MAX_CLOSURE 8
+#define MAX_CLOSURE 16
 
 typedef struct ShaderClosure {
 	ClosureType type;
@@ -636,6 +637,7 @@ typedef struct KernelIntegrator {
 	int num_all_lights;
 	float pdf_triangles;
 	float pdf_lights;
+	float inv_pdf_lights;
 	int pdf_background_res;
 
 	/* bounces */
@@ -671,7 +673,7 @@ typedef struct KernelIntegrator {
 	int transmission_samples;
 	int ao_samples;
 	int mesh_light_samples;
-	int pad1, pad2;
+	int pad1;
 } KernelIntegrator;
 
 typedef struct KernelBVH {

intern/cycles/kernel/osl/osl_services.cpp

@@ -605,7 +605,11 @@ bool OSLRenderServices::get_object_standard_attribute(KernelGlobals *kg, ShaderD
 		return set_attribute_int(3, type, derivatives, val);
 	}
 	else if ((name == u_geom_trianglevertices || name == u_geom_polyvertices)
+#ifdef __HAIR__
 		     && sd->segment == ~0) {
+#else
+		) {
+#endif
 		float3 P[3];
 		triangle_vertices(kg, sd->prim, P);
 
@@ -675,7 +679,11 @@ bool OSLRenderServices::get_attribute(void *renderstate, bool derivatives, ustri
 	else {
 		object = sd->object;
 		prim = sd->prim;
+#ifdef __HAIR__
 		segment = sd->segment;
+#else
+		segment = ~0;
+#endif
 
 		if (object == ~0)
 			return get_background_attribute(kg, sd, name, type, derivatives, val);

intern/cycles/kernel/osl/osl_shader.cpp

@@ -457,7 +457,11 @@ float3 OSLShader::volume_eval_phase(const ShaderClosure *sc, const float3 omega_
 int OSLShader::find_attribute(KernelGlobals *kg, const ShaderData *sd, uint id, AttributeElement *elem)
 {
 	/* for OSL, a hash map is used to lookup the attribute by name. */
-	int object = sd->object*ATTR_PRIM_TYPES + (sd->segment != ~0);
+	int object = sd->object*ATTR_PRIM_TYPES;
+#ifdef __HAIR__
+	if(sd->segment != ~0) object += ATTR_PRIM_CURVE;
+#endif
+
 	OSLGlobals::AttributeMap &attr_map = kg->osl->attribute_map[object];
 	ustring stdname(std::string("geom:") + std::string(Attribute::standard_name((AttributeStandard)id)));
 	OSLGlobals::AttributeMap::const_iterator it = attr_map.find(stdname);

intern/cycles/render/light.cpp

@@ -323,6 +323,7 @@ void LightManager::device_update_distribution(Device *device, DeviceScene *dscen
 		/* precompute pdfs */
 		kintegrator->pdf_triangles = 0.0f;
 		kintegrator->pdf_lights = 0.0f;
+		kintegrator->inv_pdf_lights = 0.0f;
 
 		/* sample one, with 0.5 probability of light or triangle */
 		kintegrator->num_all_lights = num_lights;
@@ -337,6 +338,8 @@ void LightManager::device_update_distribution(Device *device, DeviceScene *dscen
 			kintegrator->pdf_lights = 1.0f/num_lights;
 			if(trianglearea > 0.0f)
 				kintegrator->pdf_lights *= 0.5f;
+
+			kintegrator->inv_pdf_lights = 1.0f/kintegrator->pdf_lights;
 		}
 
 		/* CDF */
@@ -349,6 +352,7 @@ void LightManager::device_update_distribution(Device *device, DeviceScene *dscen
 		kintegrator->num_all_lights = 0;
 		kintegrator->pdf_triangles = 0.0f;
 		kintegrator->pdf_lights = 0.0f;
+		kintegrator->inv_pdf_lights = 0.0f;
 	}
 }
 
@@ -475,16 +479,25 @@ void LightManager::device_update_points(Device *device, DeviceScene *dscene, Sce
 		if(light->type == LIGHT_POINT) {
 			shader_id &= ~SHADER_AREA_LIGHT;
 
+			float radius = light->size;
+			float invarea = (radius > 0.0f)? 1.0f/(M_PI_F*radius*radius): 1.0f;
+
 			light_data[i*LIGHT_SIZE + 0] = make_float4(__int_as_float(light->type), co.x, co.y, co.z);
-			light_data[i*LIGHT_SIZE + 1] = make_float4(__int_as_float(shader_id), light->size, 0.0f, 0.0f);
+			light_data[i*LIGHT_SIZE + 1] = make_float4(__int_as_float(shader_id), radius, invarea, 0.0f);
 			light_data[i*LIGHT_SIZE + 2] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
 			light_data[i*LIGHT_SIZE + 3] = make_float4(samples, 0.0f, 0.0f, 0.0f);
 		}
 		else if(light->type == LIGHT_DISTANT) {
 			shader_id &= ~SHADER_AREA_LIGHT;
 
+			float radius = light->size;
+			float angle = atanf(radius);
+			float cosangle = cosf(angle);
+			float area = M_PI_F*radius*radius;
+			float invarea = (area > 0.0f)? 1.0f/area: 1.0f;
+
 			light_data[i*LIGHT_SIZE + 0] = make_float4(__int_as_float(light->type), dir.x, dir.y, dir.z);
-			light_data[i*LIGHT_SIZE + 1] = make_float4(__int_as_float(shader_id), light->size, 0.0f, 0.0f);
+			light_data[i*LIGHT_SIZE + 1] = make_float4(__int_as_float(shader_id), radius, cosangle, invarea);
 			light_data[i*LIGHT_SIZE + 2] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
 			light_data[i*LIGHT_SIZE + 3] = make_float4(samples, 0.0f, 0.0f, 0.0f);
 		}
@@ -499,21 +512,25 @@ void LightManager::device_update_points(Device *device, DeviceScene *dscene, Sce
 		else if(light->type == LIGHT_AREA) {
 			float3 axisu = light->axisu*(light->sizeu*light->size);
 			float3 axisv = light->axisv*(light->sizev*light->size);
+			float area = len(axisu)*len(axisv);
+			float invarea = (area > 0.0f)? 1.0f/area: 0.0f;
 
 			light_data[i*LIGHT_SIZE + 0] = make_float4(__int_as_float(light->type), co.x, co.y, co.z);
 			light_data[i*LIGHT_SIZE + 1] = make_float4(__int_as_float(shader_id), axisu.x, axisu.y, axisu.z);
-			light_data[i*LIGHT_SIZE + 2] = make_float4(0.0f, axisv.x, axisv.y, axisv.z);
+			light_data[i*LIGHT_SIZE + 2] = make_float4(invarea, axisv.x, axisv.y, axisv.z);
 			light_data[i*LIGHT_SIZE + 3] = make_float4(samples, dir.x, dir.y, dir.z);
 		}
 		else if(light->type == LIGHT_SPOT) {
 			shader_id &= ~SHADER_AREA_LIGHT;
 
+			float radius = light->size;
+			float invarea = (radius > 0.0f)? 1.0f/(M_PI_F*radius*radius): 1.0f;
 			float spot_angle = cosf(light->spot_angle*0.5f);
 			float spot_smooth = (1.0f - spot_angle)*light->spot_smooth;
 
 			light_data[i*LIGHT_SIZE + 0] = make_float4(__int_as_float(light->type), co.x, co.y, co.z);
-			light_data[i*LIGHT_SIZE + 1] = make_float4(__int_as_float(shader_id), light->size, dir.x, dir.y);
+			light_data[i*LIGHT_SIZE + 1] = make_float4(__int_as_float(shader_id), radius, invarea, spot_angle);
-			light_data[i*LIGHT_SIZE + 2] = make_float4(dir.z, spot_angle, spot_smooth, 0.0f);
+			light_data[i*LIGHT_SIZE + 2] = make_float4(spot_smooth, dir.x, dir.y, dir.z);
 			light_data[i*LIGHT_SIZE + 3] = make_float4(samples, 0.0f, 0.0f, 0.0f);
 		}
 	}

intern/cycles/util/util_math.h

@@ -1161,6 +1161,130 @@ __device float safe_divide(float a, float b)
 	return result;
 }
 
+/* Ray Intersection */
+
+__device bool ray_sphere_intersect(
+	float3 ray_P, float3 ray_D, float ray_t,
+	float3 sphere_P, float sphere_radius,
+	float3 *isect_P, float *isect_t)
+{
+	float3 d = sphere_P - ray_P;
+	float radiussq = sphere_radius*sphere_radius;
+	float tsq = dot(d, d);
+
+	if(tsq > radiussq) { /* ray origin outside sphere */
+		float tp = dot(d, ray_D);
+
+		if(tp < 0.0f) /* dir points away from sphere */
+			return false;
+
+		float dsq = tsq - tp*tp; /* pythagoras */
+
+		if(dsq > radiussq) /* closest point on ray outside sphere */
+			return false;
+
+		float t = tp - sqrtf(radiussq - dsq); /* pythagoras */
+
+		if(t < ray_t) {
+			*isect_t = t;
+			*isect_P = ray_P + ray_D*t;
+			return true;
+		}
+	}
+
+	return false;
+}
+
+__device bool ray_aligned_disk_intersect(
+	float3 ray_P, float3 ray_D, float ray_t,
+	float3 disk_P, float disk_radius,
+	float3 *isect_P, float *isect_t)
+{
+	/* aligned disk normal */
+	float disk_t;
+	float3 disk_N = normalize_len(ray_P - disk_P, &disk_t);
+	float div = dot(ray_D, disk_N);
+
+	if(div == 0.0f)
+		return false;
+
+	/* compute t to intersection point */
+	float t = -disk_t/div;
+	if(t < 0.0f || t > ray_t)
+		return false;
+	
+	/* test if within radius */
+	float3 P = ray_P + ray_D*t;
+	if(len_squared(P - disk_P) > disk_radius*disk_radius)
+		return false;
+
+	*isect_P = P;
+	*isect_t = t;
+
+	return true;
+}
+
+__device bool ray_triangle_intersect(
+	float3 ray_P, float3 ray_D, float ray_t,
+	float3 v0, float3 v1, float3 v2,
+	float3 *isect_P, float *isect_t)
+{
+	/* Calculate intersection */
+	float3 e1 = v1 - v0;
+	float3 e2 = v2 - v0;
+	float3 s1 = cross(ray_D, e2);
+
+	const float divisor = dot(s1, e1);
+	if(divisor == 0.0f)
+		return false;
+
+	const float invdivisor = 1.0f/divisor;
+
+	/* compute first barycentric coordinate */
+	const float3 d = ray_P - v0;
+	const float u = dot(d, s1)*invdivisor;
+	if(u < 0.0f)
+		return false;
+
+	/* Compute second barycentric coordinate */
+	const float3 s2 = cross(d, e1);
+	const float v = dot(ray_D, s2)*invdivisor;
+	if(v < 0.0f)
+		return false;
+
+	const float b0 = 1.0f - u - v;
+	if(b0 < 0.0f)
+		return false;
+
+	/* compute t to intersection point */
+	const float t = dot(e2, s2)*invdivisor;
+	if(t < 0.0f || t > ray_t)
+		return false;
+
+	*isect_t = t;
+	*isect_P = ray_P + ray_D*t;
+
+	return true;
+}
+
+__device bool ray_quad_intersect(
+	float3 ray_P, float3 ray_D, float ray_t,
+	float3 quad_P, float3 quad_u, float3 quad_v,
+	float3 *isect_P, float *isect_t)
+{
+	float3 v0 = quad_P - quad_u*0.5f - quad_v*0.5f;
+	float3 v1 = quad_P + quad_u*0.5f - quad_v*0.5f;
+	float3 v2 = quad_P + quad_u*0.5f + quad_v*0.5f;
+	float3 v3 = quad_P - quad_u*0.5f + quad_v*0.5f;
+
+	if(ray_triangle_intersect(ray_P, ray_D, ray_t, v0, v1, v2, isect_P, isect_t))
+		return true;
+	else if(ray_triangle_intersect(ray_P, ray_D, ray_t, v0, v2, v3, isect_P, isect_t))
+		return true;
+	
+	return false;
+}
+
 CCL_NAMESPACE_END
 
 #endif /* __UTIL_MATH_H__ */