[Bf-blender-cvs] [f5cb0cf] master: Cycles: improved importance sampling for Beckmann and GGX glossy

Brecht Van Lommel noreply at git.blender.org
Sat Jun 14 13:55:34 CEST 2014


Commit: f5cb0cf1a50350e32b6fec5056f23a20606c7ea0
Author: Brecht Van Lommel
Date:   Thu May 29 13:32:16 2014 +0200
https://developer.blender.org/rBf5cb0cf1a50350e32b6fec5056f23a20606c7ea0

Cycles: improved importance sampling for Beckmann and GGX glossy

Samples render slower than before, but hopefully this is made up for with
reduced noise in most cases. The main slowdown comes from samples that would
previously be wasted and turn out black, which are now continued.

GGX sampling is about the same speed as before, while for Beckmann it is slower
still. Perhaps optimizations are still possible there, but didn't find anything
easy.

Code from this paper, which comes with sample code:

Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals.
E. Heitz and E. d'Eon, EGSR 2014

Differential Revision: https://developer.blender.org/D572

===================================================================

M	intern/cycles/kernel/closure/bsdf_microfacet.h

===================================================================

diff --git a/intern/cycles/kernel/closure/bsdf_microfacet.h b/intern/cycles/kernel/closure/bsdf_microfacet.h
index 1ec35e4..ea0894e 100644
--- a/intern/cycles/kernel/closure/bsdf_microfacet.h
+++ b/intern/cycles/kernel/closure/bsdf_microfacet.h
@@ -35,7 +35,312 @@
 
 CCL_NAMESPACE_BEGIN
 
-/* GGX */
+/* Approximate erf and erfinv implementations
+ *
+ * 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;
+}
+
+ccl_device float approx_erff(float z)
+{
+	float s = 1.0f;
+
+	if(z < 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;
+	}
+	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);
+	}
+	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;
+}
+
+ccl_device float approx_erfinvf(float z)
+{
+   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);
+}
+
+/* Beckmann and GGX microfacet importance sampling from:
+ * 
+ * Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals.
+ * E. Heitz and E. d'Eon, EGSR 2014 */
+
+ccl_device_inline void microfacet_beckmann_sample_slopes(
+	const float cos_theta_i, const float sin_theta_i,
+	const float alpha_x, const float alpha_y,
+	float randu, float randv, float *slope_x, float *slope_y,
+	float *G1i)
+{
+	/* 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);
+		*G1i = 1.0f;
+		return;
+	}
+
+	/* precomputations */
+	const float tan_theta_i = sin_theta_i/cos_theta_i;
+	const float inv_a = tan_theta_i;
+	const float a = 1.0f/inv_a;
+	const float erf_a = approx_erff(a);
+	const float exp_a2 = expf(-a*a);
+	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);
+	const float G1 = 1.0f/(1.0f + Lambda); /* masking */
+	const float C = 1.0f - G1 * erf_a;
+
+	*G1i = G1;
+
+	/* sample slope X */
+	if(randu < C) {
+		/* rescale randu */
+		randu = randu / C;
+		const float w_1 = 0.5f * SQRT_PI_INV * sin_theta_i * exp_a2;
+		const float w_2 = cos_theta_i * (0.5f - 0.5f*erf_a);
+		const float p = w_1 / (w_1 + w_2);
+
+		if(randu < p) {
+			randu = randu / p;
+			*slope_x = -sqrtf(-logf(randu*exp_a2));
+		}
+		else {
+			randu = (randu - p) / (1.0f - p);
+			*slope_x = approx_erfinvf(randu - 1.0f - randu*erf_a);
+		}
+	}
+	else {
+		/* rescale randu */
+		randu = (randu - C) / (1.0f - C);
+		*slope_x = approx_erfinvf((-1.0f + 2.0f*randu)*erf_a);
+
+		const float p = (-(*slope_x)*sin_theta_i + cos_theta_i) / (2.0f*cos_theta_i);
+
+		if(randv > p) {
+			*slope_x = -(*slope_x);
+			randv = (randv - p) / (1.0f - p);
+		}
+		else
+			randv = randv / p;
+	}
+
+	/* sample slope Y */
+	*slope_y = approx_erfinvf(2.0f*randv - 1.0f);
+}
+
+ccl_device_inline void microfacet_ggx_sample_slopes(
+	const float cos_theta_i, const float sin_theta_i,
+	const float alpha_x, const float alpha_y,
+	float randu, float randv, float *slope_x, float *slope_y,
+	float *G1i)
+{
+	/* 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);
+		*G1i = 1.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));
+
+	*G1i = 1.0f/G1_inv;
+
+	/* 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));
+}
+
+ccl_device_inline float3 microfacet_sample_stretched(const float3 omega_i,
+	const float alpha_x, const float alpha_y,
+	const float randu, const float randv,
+	bool beckmann, float *G1i)
+{
+	/* 1. stretch omega_i */
+	float3 omega_i_ = make_float3(alpha_x * omega_i.x, alpha_y * omega_i.y, omega_i.z);
+	omega_i_ = normalize(omega_i_);
+
+	/* get polar coordinates of omega_i_ */
+	float costheta_ = 1.0f;
+	float sintheta_ = 0.0f;
+	float cosphi_ = 1.0f;
+	float sinphi_ = 0.0f;
+
+	if(omega_i_.z < 0.99999f) {
+		costheta_ = omega_i_.z;
+		sintheta_ = safe_sqrtf(1.0f - costheta_*costheta_);
+
+		float invlen = 1.0f/sintheta_;
+		cosphi_ = omega_i_.x * invlen;
+		sinphi_ = omega_i_.y * invlen;
+	}
+
+	/* 2. sample P22_{omega_i}(x_slope, y_slope, 1, 1) */
+	float slope_x, slope_y;
+
+	if(beckmann)
+		microfacet_beckmann_sample_slopes(costheta_, sintheta_,
+			alpha_x, alpha_y, randu, randv, &slope_x, &slope_y, G1i);
+	else
+		microfacet_ggx_sample_slopes(costheta_, sintheta_,
+			alpha_x, alpha_y, randu, randv, &slope_x, &slope_y, G1i);
+
+	/* 3. rotate */
+	float tmp = cosphi_*slope_x - sinphi_*slope_y;
+	slope_y = sinphi_*slope_x + cosphi_*slope_y;
+	slope_x = tmp;
+
+	/* 4. unstretch */
+	slope_x = alpha_x * slope_x;
+	slope_y = alpha_y * slope_y;
+
+	/* 5. compute normal */
+	return normalize(make_float3(-slope_x, -slope_y, 1.0f));
+} 
+
+/* GGX microfacet with Smith shadow-masking from:
+ *
+ * Microfacet Models for Refraction through Rough Surfaces
+ * B. Walter, S. R. Marschner, H. Li, K. E. Torrance, EGSR 2007 */
 
 ccl_device int bsdf_microfacet_ggx_setup(ShaderClosure *sc)
 {
@@ -67,34 +372,44 @@ ccl_device float3 bsdf_microfacet_ggx_eval_reflect(const ShaderClosure *sc, cons
 	float3 N = sc->N;
 
 	if(m_refractive || m_ag <= 1e-4f)
-		return make_float3 (0, 0, 0);
+		return make_float3(0, 0, 0);
+
 	float cosNO = dot(N, I);
 	float cosNI = dot(N, omega_in);
+
 	if(cosNI > 0 && cosNO > 0) {
-		// get half vector
-		float3 Hr = normalize(omega_in + I);
-		// eq. 20: (F*G*D)/(4*in*on)
-		// eq. 33: first we calculate D(m) with m=Hr:
+		/* get half vector */
+		float3 m = normalize(omega_in + I);
+
+		/* eq. 20: (F*G*D)/(4*in*on)
+		 * eq. 33: first we calculate D(m) */
 		float alpha2 = m_ag * m_ag;
-		float cosThetaM = dot(N, Hr);
+		float cosThetaM = dot(N, m);
 		float cosThetaM2 = cosThetaM * cosThetaM;
 		float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
 		float cosThetaM4 = cosThetaM2 * cosThetaM2;
 		float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
-		// eq. 34: now calculate G1(i,m) and G1(o,m)
+		/* eq. 34: now calculate G1(i,m) and G1(o,m) */
 		float G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
 		float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
 		float G = G1o * G1i;
-		float out = (G * D) * 0.25f / cosNO;
-		// eq. 24
-		float pm = D * cosThetaM;
-		// convert into pdf of the sampled direction
-		// eq. 38 - but see also:
-		// eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
-		*pdf = pm * 0.25f / dot(Hr, I);
-		return make_float3 (out, out, out);
+
+		/*

@@ Diff output truncated at 10240 characters. @@




More information about the Bf-blender-cvs mailing list