Cycles/Eevee: Implement disk and ellipse shapes for area lamps

The implementation is pretty straightforward.

In Cycles, sampling the shapes is currently done w.r.t. area instead of solid angle.

There is a paper on solid angle sampling for disks [1], but the described algorithm is based on
simply sampling the enclosing square and rejecting samples outside of the disk, which is not exactly
great for Cycles' RNG (we'd need to setup a LCG for the repeated sampling) and for GPU divergence.

Even worse, the algorithm is only defined for disks. For ellipses, the basic idea still works, but a
way to analytically calculate the solid angle is required. This is technically possible [2], but the
calculation is extremely complex and still requires a lookup table for the Heuman Lambda function.

Therefore, I've decided to not implement that for now, we could still look into it later on.

In Eevee, the code uses the existing ltc_evaluate_disk to implement the lighting calculations.

[1]: "Solid Angle Sampling of Disk and Cylinder Lights"
[2]: "Analytical solution for the solid angle subtended at any point by an ellipse via a point source radiation vector potential"

Reviewers: sergey, brecht, fclem

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

6 files changed, 158 insertions, 84 deletions

diff --git a/intern/cycles/blender/addon/ui.py b/intern/cycles/blender/addon/ui.py
index eeb059f79f1..ff8ca0adc18 100644
--- a/intern/cycles/blender/addon/ui.py
+++ b/intern/cycles/blender/addon/ui.py
@@ -909,9 +909,9 @@ class CYCLES_LAMP_PT_lamp(CyclesButtonsPanel, Panel):
             col.prop(lamp, "shape", text="")
             sub = col.column(align=True)
 
-            if lamp.shape == 'SQUARE':
+            if lamp.shape in {'SQUARE', 'DISK'}:
                 sub.prop(lamp, "size")
-            elif lamp.shape == 'RECTANGLE':
+            elif lamp.shape in {'RECTANGLE', 'ELLIPSE'}:
                 sub.prop(lamp, "size", text="Size X")
                 sub.prop(lamp, "size_y", text="Size Y")
 
diff --git a/intern/cycles/blender/blender_object.cpp b/intern/cycles/blender/blender_object.cpp
index d0df8e1800f..090682c8e14 100644
--- a/intern/cycles/blender/blender_object.cpp
+++ b/intern/cycles/blender/blender_object.cpp
@@ -162,10 +162,24 @@ void BlenderSync::sync_light(BL::Object& b_parent,
 			light->axisu = transform_get_column(&tfm, 0);
 			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();
-			else
-				light->sizev = light->sizeu;
+			switch(b_area_lamp.shape()) {
+				case BL::AreaLamp::shape_SQUARE:
+					light->sizev = light->sizeu;
+					light->round = false;
+					break;
+				case BL::AreaLamp::shape_RECTANGLE:
+					light->sizev = b_area_lamp.size_y();
+					light->round = false;
+					break;
+				case BL::AreaLamp::shape_DISK:
+					light->sizev = light->sizeu;
+					light->round = true;
+					break;
+				case BL::AreaLamp::shape_ELLIPSE:
+					light->sizev = b_area_lamp.size_y();
+					light->round = true;
+					break;
+			}
 			light->type = LIGHT_AREA;
 			break;
 		}
diff --git a/intern/cycles/kernel/kernel_light.h b/intern/cycles/kernel/kernel_light.h
index efab69ee37d..cd879e27e72 100644
--- a/intern/cycles/kernel/kernel_light.h
+++ b/intern/cycles/kernel/kernel_light.h
@@ -44,7 +44,7 @@ typedef struct LightSample {
  *
  * Note: light_p is modified when sample_coord is true.
  */
-ccl_device_inline float area_light_sample(float3 P,
+ccl_device_inline float rect_light_sample(float3 P,
                                           float3 *light_p,
                                           float3 axisu, float3 axisv,
                                           float randu, float randv,
@@ -125,6 +125,60 @@ ccl_device_inline float area_light_sample(float3 P,
 		return 0.0f;
 }
 
+ccl_device_inline float3 ellipse_sample(float3 ru, float3 rv, float randu, float randv)
+{
+	to_unit_disk(&randu, &randv);
+	return ru*randu + rv*randv;
+}
+
+ccl_device float3 disk_light_sample(float3 v, float randu, float randv)
+{
+	float3 ru, rv;
+
+	make_orthonormals(v, &ru, &rv);
+
+	return ellipse_sample(ru, rv, randu, randv);
+}
+
+ccl_device float3 distant_light_sample(float3 D, float radius, float randu, float randv)
+{
+	return normalize(D + disk_light_sample(D, randu, randv)*radius);
+}
+
+ccl_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;
+}
+
+ccl_device float spot_light_attenuation(float3 dir, float spot_angle, float spot_smooth, LightSample *ls)
+{
+	float3 I = ls->Ng;
+
+	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;
+}
+
+ccl_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;
+}
+
 /* Background Light */
 
 #ifdef __BACKGROUND_MIS__
@@ -295,11 +349,19 @@ ccl_device_inline float background_portal_pdf(KernelGlobals *kg,
 		const ccl_global KernelLight *klight = &kernel_tex_fetch(__lights, portal);
 		float3 axisu = make_float3(klight->area.axisu[0], klight->area.axisu[1], klight->area.axisu[2]);
 		float3 axisv = make_float3(klight->area.axisv[0], klight->area.axisv[1], klight->area.axisv[2]);
+		bool is_round = (klight->area.invarea < 0.0f);
 
-		if(!ray_quad_intersect(P, direction, 1e-4f, FLT_MAX, lightpos, axisu, axisv, dir, NULL, NULL, NULL, NULL))
+		if(!ray_quad_intersect(P, direction, 1e-4f, FLT_MAX, lightpos, axisu, axisv, dir, NULL, NULL, NULL, NULL, is_round))
 			continue;
 
-		portal_pdf += area_light_sample(P, &lightpos, axisu, axisv, 0.0f, 0.0f, false);
+		if(is_round) {
+			float t;
+			float3 D = normalize_len(lightpos - P, &t);
+			portal_pdf += fabsf(klight->area.invarea) * lamp_light_pdf(kg, dir, -D, t);
+		}
+		else {
+			portal_pdf += rect_light_sample(P, &lightpos, axisu, axisv, 0.0f, 0.0f, false);
+		}
 	}
 
 	if(ignore_portal >= 0) {
@@ -349,15 +411,26 @@ ccl_device float3 background_portal_sample(KernelGlobals *kg,
 			const ccl_global KernelLight *klight = &kernel_tex_fetch(__lights, portal);
 			float3 axisu = make_float3(klight->area.axisu[0], klight->area.axisu[1], klight->area.axisu[2]);
 			float3 axisv = make_float3(klight->area.axisv[0], klight->area.axisv[1], klight->area.axisv[2]);
-
-			*pdf = area_light_sample(P, &lightpos,
-			                         axisu, axisv,
-			                         randu, randv,
-			                         true);
+			bool is_round = (klight->area.invarea < 0.0f);
+
+			float3 D;
+			if(is_round) {
+				lightpos += ellipse_sample(axisu*0.5f, axisv*0.5f, randu, randv);
+				float t;
+				D = normalize_len(lightpos - P, &t);
+				*pdf = fabsf(klight->area.invarea) * lamp_light_pdf(kg, dir, -D, t);
+			}
+			else {
+				*pdf = rect_light_sample(P, &lightpos,
+				                         axisu, axisv,
+				                         randu, randv,
+				                         true);
+				D = normalize(lightpos - P);
+			}
 
 			*pdf /= num_possible;
 			*sampled_portal = p;
-			return normalize(lightpos - P);
+			return D;
 		}
 
 		portal--;
@@ -458,55 +531,6 @@ ccl_device float background_light_pdf(KernelGlobals *kg, float3 P, float3 direct
 
 /* Regular Light */
 
-ccl_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;
-}
-
-ccl_device float3 distant_light_sample(float3 D, float radius, float randu, float randv)
-{
-	return normalize(D + disk_light_sample(D, randu, randv)*radius);
-}
-
-ccl_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;
-}
-
-ccl_device float spot_light_attenuation(float3 dir, float spot_angle, float spot_smooth, LightSample *ls)
-{
-	float3 I = ls->Ng;
-
-	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;
-}
-
-ccl_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;
-}
-
 ccl_device_inline bool lamp_light_sample(KernelGlobals *kg,
                                          int lamp,
                                          float randu, float randv,
@@ -601,26 +625,39 @@ ccl_device_inline bool lamp_light_sample(KernelGlobals *kg,
 			float3 D = make_float3(klight->area.dir[0],
 			                       klight->area.dir[1],
 			                       klight->area.dir[2]);
+			float invarea = fabsf(klight->area.invarea);
+			bool is_round = (klight->area.invarea < 0.0f);
 
 			if(dot(ls->P - P, D) > 0.0f) {
 				return false;
 			}
 
-			float3 inplane = ls->P;
-			ls->pdf = area_light_sample(P, &ls->P,
-			                          axisu, axisv,
-			                          randu, randv,
-			                          true);
+			float3 inplane;
+
+			if(is_round) {
+				inplane = ellipse_sample(axisu*0.5f, axisv*0.5f, randu, randv);
+				ls->P += inplane;
+				ls->pdf = invarea;
+			}
+			else {
+				inplane = ls->P;
+				ls->pdf = rect_light_sample(P, &ls->P,
+				                            axisu, axisv,
+				                            randu, randv,
+				                            true);
+				inplane = ls->P - inplane;
+			}
 
-			inplane = ls->P - inplane;
 			ls->u = dot(inplane, axisu) * (1.0f / dot(axisu, axisu)) + 0.5f;
 			ls->v = dot(inplane, axisv) * (1.0f / dot(axisv, axisv)) + 0.5f;
 
 			ls->Ng = D;
 			ls->D = normalize_len(ls->P - P, &ls->t);
 
-			float invarea = klight->area.invarea;
 			ls->eval_fac = 0.25f*invarea;
+			if(is_round) {
+				ls->pdf *= lamp_light_pdf(kg, D, -ls->D, ls->t);
+			}
 		}
 	}
 
@@ -731,7 +768,8 @@ ccl_device bool lamp_light_eval(KernelGlobals *kg, int lamp, float3 P, float3 D,
 	}
 	else if(type == LIGHT_AREA) {
 		/* area light */
-		float invarea = klight->area.invarea;
+		float invarea = fabsf(klight->area.invarea);
+		bool is_round = (klight->area.invarea < 0.0f);
 		if(invarea == 0.0f)
 			return false;
 
@@ -754,14 +792,20 @@ ccl_device bool lamp_light_eval(KernelGlobals *kg, int lamp, float3 P, float3 D,
 		if(!ray_quad_intersect(P, D, 0.0f, t, light_P,
 		                       axisu, axisv, Ng,
 		                       &ls->P, &ls->t,
-		                       &ls->u, &ls->v))
+		                       &ls->u, &ls->v,
+		                       is_round))
 		{
 			return false;
 		}
 
 		ls->D = D;
 		ls->Ng = Ng;
-		ls->pdf = area_light_sample(P, &light_P, axisu, axisv, 0, 0, false);
+		if(is_round) {
+			ls->pdf = invarea * lamp_light_pdf(kg, Ng, -D, ls->t);
+		}
+		else {
+			ls->pdf = rect_light_sample(P, &light_P, axisu, axisv, 0, 0, false);
+		}
 		ls->eval_fac = 0.25f*invarea;
 	}
 	else {
diff --git a/intern/cycles/render/light.cpp b/intern/cycles/render/light.cpp
index 8dec7e4ea64..01a8b50cb21 100644
--- a/intern/cycles/render/light.cpp
+++ b/intern/cycles/render/light.cpp
@@ -117,6 +117,7 @@ NODE_DEFINE(Light)
 	SOCKET_FLOAT(sizeu, "Size U", 1.0f);
 	SOCKET_VECTOR(axisv, "Axis V", make_float3(0.0f, 0.0f, 0.0f));
 	SOCKET_FLOAT(sizev, "Size V", 1.0f);
+	SOCKET_BOOLEAN(round, "Round", false);
 
 	SOCKET_INT(map_resolution, "Map Resolution", 512);
 
@@ -730,12 +731,15 @@ void LightManager::device_update_points(Device *,
 			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: 1.0f;
+			if(light->round) {
+				area *= -M_PI_4_F;
+			}
+			float invarea = (area != 0.0f)? 1.0f/area: 1.0f;
 			float3 dir = light->dir;
 			
 			dir = safe_normalize(dir);
 
-			if(light->use_mis && area > 0.0f)
+			if(light->use_mis && area != 0.0f)
 				shader_id |= SHADER_USE_MIS;
 
 			klights[light_index].co[0] = co.x;
@@ -803,7 +807,10 @@ void LightManager::device_update_points(Device *,
 		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: 1.0f;
+		if(light->round) {
+			area *= -M_PI_4_F;
+		}
+		float invarea = (area != 0.0f)? 1.0f/area: 1.0f;
 		float3 dir = light->dir;
 
 		dir = safe_normalize(dir);
diff --git a/intern/cycles/render/light.h b/intern/cycles/render/light.h
index 97b7b971c73..539a33f4e41 100644
--- a/intern/cycles/render/light.h
+++ b/intern/cycles/render/light.h
@@ -49,6 +49,7 @@ public:
 	float sizeu;
 	float3 axisv;
 	float sizev;
+	bool round;
 
 	Transform tfm;
 
diff --git a/intern/cycles/util/util_math_intersect.h b/intern/cycles/util/util_math_intersect.h
index 61ddcc38f50..5f15487df57 100644
--- a/intern/cycles/util/util_math_intersect.h
+++ b/intern/cycles/util/util_math_intersect.h
@@ -186,12 +186,17 @@ ccl_device_forceinline bool ray_triangle_intersect(
 #undef dot3
 }
 
+/* Tests for an intersection between a ray and a quad defined by
+ * its midpoint, normal and sides.
+ * If ellipse is true, hits outside the ellipse that's enclosed by the
+ * quad are rejected.
+ */
 ccl_device bool ray_quad_intersect(float3 ray_P, float3 ray_D,
                                    float ray_mint, float ray_maxt,
                                    float3 quad_P,
                                    float3 quad_u, float3 quad_v, float3 quad_n,
                                    float3 *isect_P, float *isect_t,
-                                   float *isect_u, float *isect_v)
+                                   float *isect_u, float *isect_v, bool ellipse)
 {
 	/* Perform intersection test. */
 	float t = -(dot(ray_P, quad_n) - dot(quad_P, quad_n)) / dot(ray_D, quad_n);
@@ -200,20 +205,23 @@ ccl_device bool ray_quad_intersect(float3 ray_P, float3 ray_D,
 	}
 	const float3 hit = ray_P + t*ray_D;
 	const float3 inplane = hit - quad_P;
-	const float u = dot(inplane, quad_u) / dot(quad_u, quad_u) + 0.5f;
-	if(u < 0.0f || u > 1.0f) {
+	const float u = dot(inplane, quad_u) / dot(quad_u, quad_u);
+	if(u < -0.5f || u > 0.5f) {
+		return false;
+	}
+	const float v = dot(inplane, quad_v) / dot(quad_v, quad_v);
+	if(v < -0.5f || v > 0.5f) {
 		return false;
 	}
-	const float v = dot(inplane, quad_v) / dot(quad_v, quad_v) + 0.5f;
-	if(v < 0.0f || v > 1.0f) {
+	if(ellipse && (u*u + v*v > 0.25f)) {
 		return false;
 	}
 	/* Store the result. */
 	/* TODO(sergey): Check whether we can avoid some checks here. */
 	if(isect_P != NULL) *isect_P = hit;
 	if(isect_t != NULL) *isect_t = t;
-	if(isect_u != NULL) *isect_u = u;
-	if(isect_v != NULL) *isect_v = v;
+	if(isect_u != NULL) *isect_u = u + 0.5f;
+	if(isect_v != NULL) *isect_v = v + 0.5f;
 	return true;
 }