Cycles: Fix rays parallel to the surface in the triangle refine and MultiGGX code

In the triangle intersection refinement code, rays that are parallel to the triangle caused a divide by zero.
These rays might initially hit the triangle due to the watertight intersection test, but are very rare - therefore, just skipping the refinement for them works fine.

Also, a few remaining issues in the MultiGGX code are fixed that were caused by rays parallel to the surface (which happened more often there due to smooth shading).

2 files changed, 28 insertions, 11 deletions

diff --git a/intern/cycles/kernel/closure/bsdf_microfacet_multi.h b/intern/cycles/kernel/closure/bsdf_microfacet_multi.h
index deb3d24bf3d..acb50ce6faa 100644
--- a/intern/cycles/kernel/closure/bsdf_microfacet_multi.h
+++ b/intern/cycles/kernel/closure/bsdf_microfacet_multi.h
@@ -42,7 +42,7 @@ ccl_device_inline float D_ggx_aniso(const float3 wm, const float2 alpha)
 /* Sample slope distribution (based on page 14 of the supplemental implementation). */
 ccl_device_inline float2 mf_sampleP22_11(const float cosI, const float2 randU)
 {
-	if(cosI > 0.9999f) {
+	if(cosI > 0.9999f || cosI < 1e-6f) {
 		const float r = sqrtf(randU.x / (1.0f - randU.x));
 		const float phi = M_2PI_F * randU.y;
 		return make_float2(r*cosf(phi), r*sinf(phi));
@@ -117,7 +117,7 @@ ccl_device_inline float3 mf_eval_phase_glossy(const float3 w, const float lambda
 	if(dotW_WH < 0.0f)
 		return make_float3(0.0f, 0.0f, 0.0f);
 
-	float phase = max(0.0f, dotW_WH) * 0.25f / (pArea * dotW_WH);
+	float phase = max(0.0f, dotW_WH) * 0.25f / max(pArea * dotW_WH, 1e-7f);
 	if(alpha.x == alpha.y)
 		phase *= D_ggx(wh, alpha.x);
 	else
@@ -200,9 +200,9 @@ ccl_device_inline float mf_lambda(const float3 w, const float2 alpha)
 	if(w.z > 0.9999f)
 		return 0.0f;
 	else if(w.z < -0.9999f)
-		return -1.0f;
+		return -0.9999f;
 
-	const float inv_wz2 = 1.0f / (w.z*w.z);
+	const float inv_wz2 = 1.0f / max(w.z*w.z, 1e-7f);
 	const float2 wa = make_float2(w.x, w.y)*alpha;
 	float v = sqrtf(1.0f + dot(wa, wa) * inv_wz2);
 	if(w.z <= 0.0f)
@@ -271,7 +271,10 @@ ccl_device_inline float mf_ggx_albedo(float r)
 
 ccl_device_inline float mf_ggx_pdf(const float3 wi, const float3 wo, const float alpha)
 {
-	return 0.25f * D_ggx(normalize(wi+wo), alpha) / ((1.0f + mf_lambda(wi, make_float2(alpha, alpha))) * wi.z) + (1.0f - mf_ggx_albedo(alpha)) * wo.z;
+	float D = D_ggx(normalize(wi+wo), alpha);
+	float lambda = mf_lambda(wi, make_float2(alpha, alpha));
+	float albedo = mf_ggx_albedo(alpha);
+	return 0.25f * D / max((1.0f + lambda) * wi.z, 1e-7f) + (1.0f - albedo) * wo.z;
 }
 
 ccl_device_inline float mf_ggx_aniso_pdf(const float3 wi, const float3 wo, const float2 alpha)
diff --git a/intern/cycles/kernel/geom/geom_triangle_intersect.h b/intern/cycles/kernel/geom/geom_triangle_intersect.h
index caa6c9d9a5b..eb0decc800b 100644
--- a/intern/cycles/kernel/geom/geom_triangle_intersect.h
+++ b/intern/cycles/kernel/geom/geom_triangle_intersect.h
@@ -342,9 +342,16 @@ ccl_device_inline float3 triangle_refine(KernelGlobals *kg,
 	float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
 	float3 qvec = cross(tvec, edge1);
 	float3 pvec = cross(D, edge2);
-	float rt = dot(edge2, qvec) / dot(edge1, pvec);
-
-	P = P + D*rt;
+	float det = dot(edge1, pvec);
+	if(det != 0.0f) {
+		/* If determinant is zero it means ray lies in the plane of
+		 * the triangle. It is possible in theory due to watertight
+		 * nature of triangle intersection. For suc hcases we simply
+		 * don't refine intersection hoping it'll go all fine.
+		 */
+		float rt = dot(edge2, qvec) / det;
+		P = P + D*rt;
+	}
 
 	if(isect->object != OBJECT_NONE) {
 #  ifdef __OBJECT_MOTION__
@@ -400,9 +407,16 @@ ccl_device_inline float3 triangle_refine_subsurface(KernelGlobals *kg,
 	float3 tvec = make_float3(P.x - tri_c.x, P.y - tri_c.y, P.z - tri_c.z);
 	float3 qvec = cross(tvec, edge1);
 	float3 pvec = cross(D, edge2);
-	float rt = dot(edge2, qvec) / dot(edge1, pvec);
-
-	P = P + D*rt;
+	float det = dot(edge1, pvec);
+	if(det != 0.0f) {
+		/* If determinant is zero it means ray lies in the plane of
+		 * the triangle. It is possible in theory due to watertight
+		 * nature of triangle intersection. For such cases we simply
+		 * don't refine intersection hoping it'll go all fine.
+		 */
+		float rt = dot(edge2, qvec) / det;
+		P = P + D*rt;
+	}
 #endif  /* __INTERSECTION_REFINE__ */
 
 	if(isect->object != OBJECT_NONE) {