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

1 files changed, 124 insertions, 0 deletions

diff --git a/intern/cycles/util/util_math.h b/intern/cycles/util/util_math.h
index 8932c85db07..c37fa1a4dc6 100644
--- a/intern/cycles/util/util_math.h
+++ b/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__ */