From d9cc3ea2c6e8cfd17612c35a82050cb01d8070ec Mon Sep 17 00:00:00 2001 From: Lukas Stockner Date: Mon, 25 Jul 2016 16:02:25 +0200 Subject: 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). --- .../cycles/kernel/closure/bsdf_microfacet_multi.h | 13 ++++++----- .../cycles/kernel/geom/geom_triangle_intersect.h | 26 +++++++++++++++++----- 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) { -- cgit v1.2.3