Commit: 65b25df8018cdef8952d2e49ad66f55c6b0b3e2b Author: Lukas Stockner Date: Sun Oct 21 03:41:31 2018 +0200 Branches: master https://developer.blender.org/rB65b25df8018cdef8952d2e49ad66f55c6b0b3e2b
Cycles: Overhaul ensure_valid_reflection to fix issues with normal- and bumpmapping This function is supposed to prevent the black artifacts caused by strong normal- or bumpmapping, but failed in some cases. Now the code correctly handles all test files and previous issues I am aware of and also has extensive comments describing the algorithm and the math behind it. Basically, the main problem was that there can be multiple valid solutions that fulfil the reflection angle criterium, but I had assumed that only one would exist and therefore simply picked the first solution with a positive term in srqt(). Now, the code uses additional validity checks and a simple heuristic to pick the best valid solution. Additionally, the code messed up very shallow reflections even if the normal map strength was zero due to the constant limit for the outgoing ray angle, which caused shallow incoming rays to fail the initial test even when reflected directly on Ng. Now, the code accounts for this by reducing the threshold in the case of a shallow incoming ray, ensuring that at least N=Ng is always a valid solution. Reviewers: brecht Differential Revision: https://developer.blender.org/D3816 =================================================================== M intern/cycles/kernel/kernel_montecarlo.h M intern/cycles/kernel/shaders/stdosl.h =================================================================== diff --git a/intern/cycles/kernel/kernel_montecarlo.h b/intern/cycles/kernel/kernel_montecarlo.h index 9b96bb80c32..ce68aa16af8 100644 --- a/intern/cycles/kernel/kernel_montecarlo.h +++ b/intern/cycles/kernel/kernel_montecarlo.h @@ -187,7 +187,10 @@ ccl_device float2 regular_polygon_sample(float corners, float rotation, float u, ccl_device float3 ensure_valid_reflection(float3 Ng, float3 I, float3 N) { float3 R = 2*dot(N, I)*N - I; - if(dot(Ng, R) >= 0.05f) { + + /* Reflection rays may always be at least as shallow as the incoming ray. */ + float threshold = min(0.9f*dot(Ng, I), 0.01f); + if(dot(Ng, R) >= threshold) { return N; } @@ -195,22 +198,86 @@ ccl_device float3 ensure_valid_reflection(float3 Ng, float3 I, float3 N) * The X axis is found by normalizing the component of N that's orthogonal to Ng. * The Y axis isn't actually needed. */ - float3 X = normalize(N - dot(N, Ng)*Ng); - - /* Calculate N.z and N.x in the local coordinate system. */ - float Iz = dot(I, Ng); - float Ix2 = sqr(dot(I, X)), Iz2 = sqr(Iz); - float Ix2Iz2 = Ix2 + Iz2; - - float a = safe_sqrtf(Ix2*(Ix2Iz2 - sqr(0.05f))); - float b = Iz*0.05f + Ix2Iz2; - float c = (a + b > 0.0f)? (a + b) : (-a + b); + float NdotNg = dot(N, Ng); + float3 X = normalize(N - NdotNg*Ng); + + /* Calculate N.z and N.x in the local coordinate system. + * + * The goal of this computation is to find a N' that is rotated towards Ng just enough + * to lift R' above the threshold (here called t), therefore dot(R', Ng) = t. + * + * According to the standard reflection equation, this means that we want dot(2*dot(N', I)*N' - I, Ng) = t. + * + * Since the Z axis of our local coordinate system is Ng, dot(x, Ng) is just x.z, so we get 2*dot(N', I)*N'.z - I.z = t. + * + * The rotation is simple to express in the coordinate system we formed - since N lies in the X-Z-plane, we know that + * N' will also lie in the X-Z-plane, so N'.y = 0 and therefore dot(N', I) = N'.x*I.x + N'.z*I.z . + * + * Furthermore, we want N' to be normalized, so N'.x = sqrt(1 - N'.z^2). + * + * With these simplifications, we get the final equation 2*(sqrt(1 - N'.z^2)*I.x + N'.z*I.z)*N'.z - I.z = t. + * + * The only unknown here is N'.z, so we can solve for that. + * + * The equation has four solutions in general: + * + * N'.z = +-sqrt(0.5*(+-sqrt(I.x^2*(I.x^2 + I.z^2 - t^2)) + t*I.z + I.x^2 + I.z^2)/(I.x^2 + I.z^2)) + * We can simplify this expression a bit by grouping terms: + * + * a = I.x^2 + I.z^2 + * b = sqrt(I.x^2 * (a - t^2)) + * c = I.z*t + a + * N'.z = +-sqrt(0.5*(+-b + c)/a) + * + * Two solutions can immediately be discarded because they're negative so N' would lie in the lower hemisphere. + */ + float Ix = dot(I, X), Iz = dot(I, Ng); + float Ix2 = sqr(Ix), Iz2 = sqr(Iz); + float a = Ix2 + Iz2; + + float b = safe_sqrtf(Ix2*(a - sqr(threshold))); + float c = Iz*threshold + a; + + /* Evaluate both solutions. + * In many cases one can be immediately discarded (if N'.z would be imaginary or larger than one), so check for that first. + * If no option is viable (might happen in extreme cases like N being in the wrong hemisphere), give up and return Ng. */ + float fac = 0.5f/a; + float N1_z2 = fac*(b+c), N2_z2 = fac*(-b+c); + bool valid1 = (N1_z2 > 1e-5f) && (N1_z2 <= (1.0f + 1e-5f)); + bool valid2 = (N2_z2 > 1e-5f) && (N2_z2 <= (1.0f + 1e-5f)); + + float2 N_new; + if(valid1 && valid2) { + /* If both are possible, do the expensive reflection-based check. */ + float2 N1 = make_float2(safe_sqrtf(1.0f - N1_z2), safe_sqrtf(N1_z2)); + float2 N2 = make_float2(safe_sqrtf(1.0f - N2_z2), safe_sqrtf(N2_z2)); + + float R1 = 2*(N1.x*Ix + N1.y*Iz)*N1.y - Iz; + float R2 = 2*(N2.x*Ix + N2.y*Iz)*N2.y - Iz; + + valid1 = (R1 >= 1e-5f); + valid2 = (R2 >= 1e-5f); + if(valid1 && valid2) { + /* If both solutions are valid, return the one with the shallower reflection since it will be closer to the input + * (if the original reflection wasn't shallow, we would not be in this part of the function). */ + N_new = (R1 < R2)? N1 : N2; + } + else { + /* If only one reflection is valid (= positive), pick that one. */ + N_new = (R1 > R2)? N1 : N2; + } - float Nz = safe_sqrtf(0.5f * c * (1.0f / Ix2Iz2)); - float Nx = safe_sqrtf(1.0f - sqr(Nz)); + } + else if(valid1 || valid2) { + /* Only one solution passes the N'.z criterium, so pick that one. */ + float Nz2 = valid1? N1_z2 : N2_z2; + N_new = make_float2(safe_sqrtf(1.0f - Nz2), safe_sqrtf(Nz2)); + } + else { + return Ng; + } - /* Transform back into global coordinates. */ - return Nx*X + Nz*Ng; + return N_new.x*X + N_new.y*Ng; } CCL_NAMESPACE_END diff --git a/intern/cycles/kernel/shaders/stdosl.h b/intern/cycles/kernel/shaders/stdosl.h index 4a8378796ba..f1235500f2b 100644 --- a/intern/cycles/kernel/shaders/stdosl.h +++ b/intern/cycles/kernel/shaders/stdosl.h @@ -284,33 +284,63 @@ point rotate (point p, float angle, point a, point b) normal ensure_valid_reflection(normal Ng, vector I, normal N) { + /* The implementation here mirrors the one in kernel_montecarlo.h, + * check there for an explanation of the algorithm. */ + float sqr(float x) { return x*x; } vector R = 2*dot(N, I)*N - I; - if (dot(Ng, R) >= 0.05) { + + float threshold = min(0.9*dot(Ng, I), 0.01); + if(dot(Ng, R) >= threshold) { return N; } - /* Form coordinate system with Ng as the Z axis and N inside the X-Z-plane. - * The X axis is found by normalizing the component of N that's orthogonal to Ng. - * The Y axis isn't actually needed. - */ - vector X = normalize(N - dot(N, Ng)*Ng); + float NdotNg = dot(N, Ng); + vector X = normalize(N - NdotNg*Ng); - /* Calculate N.z and N.x in the local coordinate system. */ float Ix = dot(I, X), Iz = dot(I, Ng); - float Ix2 = sqr(dot(I, X)), Iz2 = sqr(dot(I, Ng)); - float Ix2Iz2 = Ix2 + Iz2; - - float a = sqrt(Ix2*(Ix2Iz2 - sqr(0.05))); - float b = Iz*0.05 + Ix2Iz2; - float c = (a + b > 0.0)? (a + b) : (-a + b); + float Ix2 = sqr(Ix), Iz2 = sqr(Iz); + float a = Ix2 + Iz2; + + float b = sqrt(Ix2*(a - sqr(threshold))); + float c = Iz*threshold + a; + + float fac = 0.5/a; + float N1_z2 = fac*(b+c), N2_z2 = fac*(-b+c); + int valid1 = (N1_z2 > 1e-5) && (N1_z2 <= (1.0 + 1e-5)); + int valid2 = (N2_z2 > 1e-5) && (N2_z2 <= (1.0 + 1e-5)); + + float N_new_x, N_new_z; + if(valid1 && valid2) { + float N1_x = sqrt(1.0 - N1_z2), N1_z = sqrt(N1_z2); + float N2_x = sqrt(1.0 - N2_z2), N2_z = sqrt(N2_z2); + + float R1 = 2*(N1_x*Ix + N1_z*Iz)*N1_z - Iz; + float R2 = 2*(N2_x*Ix + N2_z*Iz)*N2_z - Iz; + + valid1 = (R1 >= 1e-5); + valid2 = (R2 >= 1e-5); + if(valid1 && valid2) { + N_new_x = (R1 < R2)? N1_x : N2_x; + N_new_z = (R1 < R2)? N1_z : N2_z; + } + else { + N_new_x = (R1 > R2)? N1_x : N2_x; + N_new_z = (R1 > R2)? N1_z : N2_z; + } - float Nz = sqrt(0.5 * c * (1.0 / Ix2Iz2)); - float Nx = sqrt(1.0 - sqr(Nz)); + } + else if(valid1 || valid2) { + float Nz2 = valid1? N1_z2 : N2_z2; + N_new_x = sqrt(1.0 - Nz2); + N_new_z = sqrt(Nz2); + } + else { + return Ng; + } - /* Transform back into global coordinates. */ - return Nx*X + Nz*Ng; + return N_new_x*X + N_new_z*Ng; } _______________________________________________ Bf-blender-cvs mailing list Bf-blender-cvs@blender.org https://lists.blender.org/mailman/listinfo/bf-blender-cvs
