diff options
author | Sergey Sharybin <sergey.vfx@gmail.com> | 2014-10-03 16:27:08 +0400 |
---|---|---|
committer | Sergey Sharybin <sergey.vfx@gmail.com> | 2014-10-03 16:28:44 +0400 |
commit | 1e4d99368bba53edb6f903036b49edda92f27e31 (patch) | |
tree | 033e43b12b9419c9ef0848e6cd67817064b74c02 /intern/cycles/kernel | |
parent | 72f557d34e21a5fe341de0a88797334810a3c66c (diff) |
Cycles: Use more accurate implementation of erf() and erfinv()
This functions are orders of magnitude more accurate than the old ones,
and they're around the same complexity to compute.
Diffstat (limited to 'intern/cycles/kernel')
-rw-r--r-- | intern/cycles/kernel/closure/bsdf_microfacet.h | 145 |
1 files changed, 22 insertions, 123 deletions
diff --git a/intern/cycles/kernel/closure/bsdf_microfacet.h b/intern/cycles/kernel/closure/bsdf_microfacet.h index a0c59e6cbc0..c6aed1205e0 100644 --- a/intern/cycles/kernel/closure/bsdf_microfacet.h +++ b/intern/cycles/kernel/closure/bsdf_microfacet.h @@ -35,139 +35,38 @@ CCL_NAMESPACE_BEGIN -/* Approximate erf and erfinv implementations +/* Approximate erf and erfinv implementations. + * Implementation comes stright from the wikipedia: * - * Adapted from code (C) Copyright John Maddock 2006. - * Use, modification and distribution are subject to the - * Boost Software License, Version 1.0. (See accompanying file - * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) */ - -ccl_device float approx_erff_impl(float z) -{ - float result; - - if(z < 0.5f) { - if(z < 1e-10f) { - if(z == 0) { - result = 0; - } - else { - float c = 0.0033791670f; - result = z * 1.125f + z * c; - } - } - else { - float Y = 1.044948577f; - - float zz = z * z; - float num = (((-0.007727583f * zz) + -0.050999073f)*zz + -0.338165134f)*zz + 0.083430589f; - float denom = (((0.000370900f * zz) + 0.008585719f)*zz + 0.087522260f)*zz + 0.455004033f; - result = z * (Y + num / denom); - } - } - else if(z < 2.5f) { - if(z < 1.5f) { - float Y = 0.4059357643f; - float fz = z - 0.5f; - - float num = (((0.088890036f * fz) + 0.191003695f)*fz + 0.178114665f)*fz + -0.098090592f; - float denom = (((0.123850974f * fz) + 0.578052804f)*fz + 1.426280048f)*fz + 1.847590709f; - - result = Y + num / denom; - result *= expf(-z * z) / z; - } - else { - float Y = 0.506728172f; - float fz = z - 1.5f; - float num = (((0.017567943f * fz) + 0.043948189f)*fz + 0.038654037f)*fz + -0.024350047f; - float denom = (((0.325732924f * fz) + 0.982403709f)*fz + 1.539914949f)*fz + 1; - - result = Y + num / denom; - result *= expf(-z * z) / z; - } - - result = 1 - result; - } - else { - result = 1; - } - - return result; -} + * http://en.wikipedia.org/wiki/Error_function + * + * Some constants are baked into the code. + */ -ccl_device float approx_erff(float z) +ccl_device_inline float approx_erff(float x) { float s = 1.0f; - - if(z < 0.0f) { + if(x < 0.0f) { s = -1.0f; - z = -z; - } - - return s * approx_erff_impl(z); -} - -ccl_device float approx_erfinvf_impl(float p, float q) -{ - float result = 0; - - if(p <= 0.5f) { - float Y = 0.089131474f; - float g = p * (p + 10); - float num = (((-0.012692614f * p) + 0.033480662f)*p + -0.008368748f)*p + -0.000508781f; - float denom = (((1.562215583f * p) + -1.565745582f)*p + -0.970005043f)*p + 1.0f; - float r = num / denom; - result = g * Y + g * r; + x = -x; } - else if(q >= 0.25f) { - float Y = 2.249481201f; - float g = sqrtf(-2 * logf(q)); - float xs = q - 0.25f; - float num = (((17.644729840f * xs) + 8.370503283f)*xs + 0.105264680f)*xs + -0.202433508f; - float denom = (((-28.660818049f * xs) + 3.971343795f)*xs + 6.242641248f)*xs + 1.0f; - float r = num / denom; - result = g / (Y + r); + /* Such a clamp doesn't give much distortion to the output value + * and gives quite a few of the speedup. + */ + if(x > 3.0f) { + return s; } - else { - float x = sqrtf(-logf(q)); - - if(x < 3) { - float Y = 0.807220458f; - float xs = x - 1.125f; - float num = (((0.387079738f * xs) + 0.117030156f)*xs + -0.163794047f)*xs + -0.131102781f; - float denom = (((4.778465929f * xs) + 5.381683457f)*xs + 3.466254072f)*xs + 1.0f; - float R = num / denom; - result = Y * x + R * x; - } - else { - float Y = 0.939955711f; - float xs = x - 3; - float num = (((0.009508047f * xs) + 0.018557330f)*xs + -0.002224265f)*xs + -0.035035378f; - float denom = (((0.220091105f * xs) + 0.762059164f)*xs + 1.365334981f)*xs + 1.0f; - float R = num / denom; - result = Y * x + R * x; - } - } - - return result; + float t = 1.0f / (1.0f + 0.47047f*x); + return s * (1.0f - + t*(0.3480242f + t*(-0.0958798f + t*0.7478556f)) * expf(-x*x)); } -ccl_device float approx_erfinvf(float z) +ccl_device_inline float approx_erfinvf(float x) { - float p, q, s; - - if(z < 0) { - p = -z; - q = 1 - p; - s = -1; - } - else { - p = z; - q = 1 - z; - s = 1; - } - - return s * approx_erfinvf_impl(p, q); + float ln1_x2 = logf(1.0f - x*x); + float term = 4.546884979448f + ln1_x2 * 0.5f; + return copysignf(1.0f, x) * + sqrtf(sqrtf(term*term - ln1_x2 * 7.142230224076f) - term); } /* Beckmann and GGX microfacet importance sampling from: |