blender
/
blender
136
406
Fork 589
Code Issues 6.2k Pull Requests 757 Projects 19 Wiki Activity

Cycles: Use faster and exact GGX VNDF sampling algorithm 

Based on "Sampling the GGX Distribution of Visible Normals" by Eric Heitz
(https://jcgt.org/published/0007/04/01/).

Also, this removes the lambdaI computation from the Beckmann sampling code and
just recomputes it below. We already need to recompute for two other cases
(GGX and clearcoat), so this makes the code more consistent.

In terms of performance, I don't expect a notable impact since the earlier
computation also was non-trivial, and while it probably was slightly more
accurate, I'd argue that being consistent between evaluation and sampling is
more important than absolute numerical accuracy anyways.

Differential Revision: https://developer.blender.org/D17100
This commit is contained in:
Lukas Stockner 2023-01-24 17:34:24 +01:00
parent fdcb55b285
commit e308b891c8
2 changed files with 132 additions and 175 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

302
intern/cycles/kernel/closure/bsdf_microfacet.h
View File

@ -37,189 +37,100 @@ typedef struct MicrofacetBsdf {
static_assert(sizeof(ShaderClosure) >= sizeof(MicrofacetBsdf), "MicrofacetBsdf is too large!");
/* Beckmann and GGX microfacet importance sampling. */
ccl_device_inline void microfacet_beckmann_sample_slopes(KernelGlobals kg,
const float cos_theta_i,
const float sin_theta_i,
float randu,
float randv,
ccl_private float *slope_x,
ccl_private float *slope_y,
ccl_private float *lambda_i)
{
/* Special case (normal incidence). */
if (cos_theta_i >= 0.99999f) {
const float r = sqrtf(-logf(randu));
const float phi = M_2PI_F * randv;
*slope_x = r * cosf(phi);
*slope_y = r * sinf(phi);
*lambda_i = 0.0f;
return;
}
/* Precomputations. */
const float tan_theta_i = sin_theta_i / cos_theta_i;
const float inv_a = tan_theta_i;
const float cot_theta_i = 1.0f / tan_theta_i;
const float erf_a = fast_erff(cot_theta_i);
const float exp_a2 = expf(-cot_theta_i * cot_theta_i);
const float SQRT_PI_INV = 0.56418958354f;
const float Lambda = 0.5f * (erf_a - 1.0f) + (0.5f * SQRT_PI_INV) * (exp_a2 * inv_a);
*lambda_i = Lambda;
/* Based on paper from Wenzel Jakob
* An Improved Visible Normal Sampling Routine for the Beckmann Distribution
*
* http://www.mitsuba-renderer.org/~wenzel/files/visnormal.pdf
*
* Reformulation from OpenShadingLanguage which avoids using inverse
* trigonometric functions.
*/
/* Sample slope X.
*
* Compute a coarse approximation using the approximation:
* exp(-ierf(x)^2) ~= 1 - x * x
* solve y = 1 + b + K * (1 - b * b)
*/
const float K = tan_theta_i * SQRT_PI_INV;
const float y_approx = randu * (1.0f + erf_a + K * (1 - erf_a * erf_a));
const float y_exact = randu * (1.0f + erf_a + K * exp_a2);
float b = K > 0 ? (0.5f - sqrtf(K * (K - y_approx + 1.0f) + 0.25f)) / K : y_approx - 1.0f;
float inv_erf = fast_ierff(b);
float2 begin = make_float2(-1.0f, -y_exact);
float2 end = make_float2(erf_a, 1.0f + erf_a + K * exp_a2 - y_exact);
float2 current = make_float2(b, 1.0f + b + K * expf(-sqr(inv_erf)) - y_exact);
/* Find root in a monotonic interval using newton method, under given precision and maximal
* iterations. Falls back to bisection if newton step produces results outside of the valid
* interval.*/
const float precision = 1e-6f;
const int max_iter = 3;
int iter = 0;
while (fabsf(current.y) > precision && iter++ < max_iter) {
if (signf(begin.y) == signf(current.y)) {
begin.x = current.x;
begin.y = current.y;
}
else {
end.x = current.x;
}
const float newton_x = current.x - current.y / (1.0f - inv_erf * tan_theta_i);
current.x = (newton_x >= begin.x && newton_x <= end.x) ? newton_x : 0.5f * (begin.x + end.x);
inv_erf = fast_ierff(current.x);
current.y = 1.0f + current.x + K * expf(-sqr(inv_erf)) - y_exact;
}
*slope_x = inv_erf;
*slope_y = fast_ierff(2.0f * randv - 1.0f);
}
/* GGX microfacet importance sampling from:
*
/* Beckmann VNDF importance sampling algorithm from:
* Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals.
* E. Heitz and E. d'Eon, EGSR 2014
*/
* Eric Heitz and Eugene d'Eon, EGSR 2014.
* https://hal.inria.fr/hal-00996995v2/document */
ccl_device_inline void microfacet_ggx_sample_slopes(const float cos_theta_i,
const float sin_theta_i,
float randu,
float randv,
ccl_private float *slope_x,
ccl_private float *slope_y,
ccl_private float *lambda_i)
{
/* Special case (normal incidence). */
if (cos_theta_i >= 0.99999f) {
const float r = sqrtf(randu / (1.0f - randu));
const float phi = M_2PI_F * randv;
*slope_x = r * cosf(phi);
*slope_y = r * sinf(phi);
*lambda_i = 0.0f;
return;
}
/* Precomputations. */
const float tan_theta_i = sin_theta_i / cos_theta_i;
const float G1_inv = 0.5f * (1.0f + safe_sqrtf(1.0f + tan_theta_i * tan_theta_i));
*lambda_i = G1_inv - 1.0f;
/* Sample slope_x. */
const float A = 2.0f * randu * G1_inv - 1.0f;
const float AA = A * A;
const float tmp = 1.0f / (AA - 1.0f);
const float B = tan_theta_i;
const float BB = B * B;
const float D = safe_sqrtf(BB * (tmp * tmp) - (AA - BB) * tmp);
const float slope_x_1 = B * tmp - D;
const float slope_x_2 = B * tmp + D;
*slope_x = (A < 0.0f || slope_x_2 * tan_theta_i > 1.0f) ? slope_x_1 : slope_x_2;
/* Sample slope_y. */
float S;
if (randv > 0.5f) {
S = 1.0f;
randv = 2.0f * (randv - 0.5f);
}
else {
S = -1.0f;
randv = 2.0f * (0.5f - randv);
}
const float z = (randv * (randv * (randv * 0.27385f - 0.73369f) + 0.46341f)) /
(randv * (randv * (randv * 0.093073f + 0.309420f) - 1.000000f) + 0.597999f);
*slope_y = S * z * safe_sqrtf(1.0f + (*slope_x) * (*slope_x));
}
template<MicrofacetType m_type>
ccl_device_forceinline float3 microfacet_sample_stretched(KernelGlobals kg,
const float3 wi,
const float alpha_x,
const float alpha_y,
const float randu,
const float randv,
ccl_private float *lambda_i)
ccl_device_forceinline float3 microfacet_beckmann_sample_vndf(KernelGlobals kg,
const float3 wi,
const float alpha_x,
const float alpha_y,
const float randu,
const float randv)
{
/* 1. stretch wi */
float3 wi_ = make_float3(alpha_x * wi.x, alpha_y * wi.y, wi.z);
wi_ = normalize(wi_);
/* Compute polar coordinates of wi_. */
float costheta_ = 1.0f;
float sintheta_ = 0.0f;
float cosphi_ = 1.0f;
float sinphi_ = 0.0f;
if (wi_.z < 0.99999f) {
costheta_ = wi_.z;
sintheta_ = sin_from_cos(costheta_);
float invlen = 1.0f / sintheta_;
cosphi_ = wi_.x * invlen;
sinphi_ = wi_.y * invlen;
}
/* 2. sample P22_{wi}(x_slope, y_slope, 1, 1) */
float slope_x, slope_y;
float cos_phi_i = 1.0f;
float sin_phi_i = 0.0f;
if (m_type == MicrofacetType::BECKMANN) {
microfacet_beckmann_sample_slopes(
kg, costheta_, sintheta_, randu, randv, &slope_x, &slope_y, lambda_i);
if (wi_.z >= 0.99999f) {
/* Special case (normal incidence). */
const float r = sqrtf(-logf(randu));
const float phi = M_2PI_F * randv;
slope_x = r * cosf(phi);
slope_y = r * sinf(phi);
}
else {
microfacet_ggx_sample_slopes(costheta_, sintheta_, randu, randv, &slope_x, &slope_y, lambda_i);
/* Precomputations. */
const float cos_theta_i = wi_.z;
const float sin_theta_i = sin_from_cos(cos_theta_i);
const float tan_theta_i = sin_theta_i / cos_theta_i;
const float cot_theta_i = 1.0f / tan_theta_i;
const float erf_a = fast_erff(cot_theta_i);
const float exp_a2 = expf(-cot_theta_i * cot_theta_i);
const float SQRT_PI_INV = 0.56418958354f;
float invlen = 1.0f / sin_theta_i;
cos_phi_i = wi_.x * invlen;
sin_phi_i = wi_.y * invlen;
/* Based on paper from Wenzel Jakob
* An Improved Visible Normal Sampling Routine for the Beckmann Distribution
*
* http://www.mitsuba-renderer.org/~wenzel/files/visnormal.pdf
*
* Reformulation from OpenShadingLanguage which avoids using inverse
* trigonometric functions.
*/
/* Sample slope X.
*
* Compute a coarse approximation using the approximation:
* exp(-ierf(x)^2) ~= 1 - x * x
* solve y = 1 + b + K * (1 - b * b)
*/
const float K = tan_theta_i * SQRT_PI_INV;
const float y_approx = randu * (1.0f + erf_a + K * (1 - erf_a * erf_a));
const float y_exact = randu * (1.0f + erf_a + K * exp_a2);
float b = K > 0 ? (0.5f - sqrtf(K * (K - y_approx + 1.0f) + 0.25f)) / K : y_approx - 1.0f;
float inv_erf = fast_ierff(b);
float2 begin = make_float2(-1.0f, -y_exact);
float2 end = make_float2(erf_a, 1.0f + erf_a + K * exp_a2 - y_exact);
float2 current = make_float2(b, 1.0f + b + K * expf(-sqr(inv_erf)) - y_exact);
/* Find root in a monotonic interval using newton method, under given precision and maximal
* iterations. Falls back to bisection if newton step produces results outside of the valid
* interval.*/
const float precision = 1e-6f;
const int max_iter = 3;
int iter = 0;
while (fabsf(current.y) > precision && iter++ < max_iter) {
if (signf(begin.y) == signf(current.y)) {
begin.x = current.x;
begin.y = current.y;
}
else {
end.x = current.x;
}
const float newton_x = current.x - current.y / (1.0f - inv_erf * tan_theta_i);
current.x = (newton_x >= begin.x && newton_x <= end.x) ? newton_x : 0.5f * (begin.x + end.x);
inv_erf = fast_ierff(current.x);
current.y = 1.0f + current.x + K * expf(-sqr(inv_erf)) - y_exact;
}
slope_x = inv_erf;
slope_y = fast_ierff(2.0f * randv - 1.0f);
}
/* 3. rotate */
float tmp = cosphi_ * slope_x - sinphi_ * slope_y;
slope_y = sinphi_ * slope_x + cosphi_ * slope_y;
float tmp = cos_phi_i * slope_x - sin_phi_i * slope_y;
slope_y = sin_phi_i * slope_x + cos_phi_i * slope_y;
slope_x = tmp;
/* 4. unstretch */
@ -230,6 +141,43 @@ ccl_device_forceinline float3 microfacet_sample_stretched(KernelGlobals kg,
return normalize(make_float3(-slope_x, -slope_y, 1.0f));
}
/* GGX VNDF importance sampling algorithm from:
* Sampling the GGX Distribution of Visible Normals.
* Eric Heitz, JCGT Vol. 7, No. 4, 2018.
* https://jcgt.org/published/0007/04/01/ */
ccl_device_forceinline float3 microfacet_ggx_sample_vndf(const float3 wi,
const float alpha_x,
const float alpha_y,
const float randu,
const float randv)
{
/* Section 3.2: Transforming the view direction to the hemisphere configuration. */
float3 wi_ = normalize(make_float3(alpha_x * wi.x, alpha_y * wi.y, wi.z));
/* Section 4.1: Orthonormal basis. */
float lensq = sqr(wi_.x) + sqr(wi_.y);
float3 T1, T2;
if (lensq > 1e-7f) {
T1 = make_float3(-wi_.y, wi_.x, 0.0f) * inversesqrtf(lensq);
T2 = cross(wi_, T1);
}
else {
/* Normal incidence, any basis is fine. */
T1 = make_float3(1.0f, 0.0f, 0.0f);
T2 = make_float3(0.0f, 1.0f, 0.0f);
}
/* Section 4.2: Parameterization of the projected area. */
float2 t = concentric_sample_disk(randu, randv);
t.y = mix(safe_sqrtf(1.0f - sqr(t.x)), t.y, 0.5f * (1.0f + wi_.z));
/* Section 4.3: Reprojection onto hemisphere. */
float3 H_ = t.x * T1 + t.y * T2 + safe_sqrtf(1.0f - len_squared(t)) * wi_;
/* Section 3.4: Transforming the normal back to the ellipsoid configuration. */
return normalize(make_float3(alpha_x * H_.x, alpha_y * H_.y, max(0.0f, H_.z)));
}
/* Calculate the reflection color
*
* If fresnel is used, the color is an interpolation of the F0 color and white
@ -476,10 +424,15 @@ ccl_device int bsdf_microfacet_sample(KernelGlobals kg,
/* Importance sampling with distribution of visible normals. Vectors are transformed to local
* space before and after sampling. */
float lambdaI;
const float3 local_I = make_float3(dot(X, wi), dot(Y, wi), cos_NI);
const float3 local_H = microfacet_sample_stretched<m_type>(
kg, local_I, alpha_x, alpha_y, randu, randv, &lambdaI);
float3 local_H;
if (m_type == MicrofacetType::GGX) {
local_H = microfacet_ggx_sample_vndf(local_I, alpha_x, alpha_y, randu, randv);
}
else {
/* m_type == MicrofacetType::BECKMANN */
local_H = microfacet_beckmann_sample_vndf(kg, local_I, alpha_x, alpha_y, randu, randv);
}
const float3 H = X * local_H.x + Y * local_H.y + N * local_H.z;
const float cos_NH = local_H.z;
@ -524,7 +477,7 @@ ccl_device int bsdf_microfacet_sample(KernelGlobals kg,
else {
label |= LABEL_GLOSSY;
float cos_NO = dot(N, *wo);
float D, lambdaO;
float D, lambdaI, lambdaO;
/* TODO: add support for anisotropic transmission. */
if (alpha_x == alpha_y || m_refractive) { /* Isotropic. */
@ -536,15 +489,13 @@ ccl_device int bsdf_microfacet_sample(KernelGlobals kg,
/* The masking-shadowing term for clearcoat has a fixed alpha of 0.25
* => alpha2 = 0.25 * 0.25 */
alpha2 = 0.0625f;
/* Recalculate lambdaI. */
lambdaI = bsdf_lambda<m_type>(alpha2, cos_NI);
}
else {
D = bsdf_D<m_type>(alpha2, cos_NH);
}
lambdaO = bsdf_lambda<m_type>(alpha2, cos_NO);
lambdaI = bsdf_lambda<m_type>(alpha2, cos_NI);
}
else { /* Anisotropic. */
const float3 local_O = make_float3(dot(X, *wo), dot(Y, *wo), cos_NO);
@ -552,6 +503,7 @@ ccl_device int bsdf_microfacet_sample(KernelGlobals kg,
D = bsdf_aniso_D<m_type>(alpha_x, alpha_y, local_H);
lambdaO = bsdf_aniso_lambda<m_type>(alpha_x, alpha_y, local_O);
lambdaI = bsdf_aniso_lambda<m_type>(alpha_x, alpha_y, local_I);
}
const float cos_HO = dot(H, *wo);

5
intern/cycles/util/math_float2.h
View File

@ -134,6 +134,11 @@ ccl_device_inline float len(const float2 a)
return sqrtf(dot(a, a));
}
ccl_device_inline float len_squared(const float2 a)
{
return dot(a, a);
}
#if !defined(__KERNEL_METAL__)
ccl_device_inline float distance(const float2 a, const float2 b)
{