Differential D572 Diff 1812 intern/cycles/kernel/closure/bsdf_microfacet.h

Changeset View

Standalone View

intern/cycles/kernel/closure/bsdf_microfacet.h

Show All 29 Lines
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.		* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/		*/

#ifndef __BSDF_MICROFACET_H__		#ifndef __BSDF_MICROFACET_H__
#define __BSDF_MICROFACET_H__		#define __BSDF_MICROFACET_H__

CCL_NAMESPACE_BEGIN		CCL_NAMESPACE_BEGIN

		// erf_inv (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 double evaluate_polynomial(const double poly[], size_t count, double z)
		{
		double sum = poly[count - 1];

		for(int i = (int)count - 2; i >= 0; --i)
		{
		sum *= z;
		sum += poly[i];
		}

		return sum;
		}

		ccl_device double erf_inv_impl(const double& p, const double& q)
		{
		double result = 0;

		if(p <= 0.5)
		{
		//
		// Evaluate inverse erf using the rational approximation:
		//
		// x = p(p+10)(Y+R(p))
		//
		// Where Y is a constant, and R(p) is optimised for a low
		// absolute error compared to \|Y\|.
		//
		// double: Max error found: 2.001849e-18
		// long double: Max error found: 1.017064e-20
		// Maximum Deviation Found (actual error term at infinite precision) 8.030e-21
		//
		static const float Y = 0.0891314744949340820313f;
		static const double P[] = {
		-0.000508781949658280665617,
		-0.00836874819741736770379,
		0.0334806625409744615033,
		-0.0126926147662974029034,
		-0.0365637971411762664006,
		0.0219878681111168899165,
		0.00822687874676915743155,
		-0.00538772965071242932965
		};
		static const double Q[] = {
		1,
		-0.970005043303290640362,
		-1.56574558234175846809,
		1.56221558398423026363,
		0.662328840472002992063,
		-0.71228902341542847553,
		-0.0527396382340099713954,
		0.0795283687341571680018,
		-0.00233393759374190016776,
		0.000886216390456424707504
		};
		double g = p * (p + 10);
		double r = evaluate_polynomial(P, sizeof(P)/sizeof(double), p) / evaluate_polynomial(Q, sizeof(Q)/sizeof(double), p);
		result = g * Y + g * r;
		}
		else if(q >= 0.25)
		{
		//
		// Rational approximation for 0.5 > q >= 0.25
		//
		// x = sqrt(-2*log(q)) / (Y + R(q))
		//
		// Where Y is a constant, and R(q) is optimised for a low
		// absolute error compared to Y.
		//
		// double : Max error found: 7.403372e-17
		// long double : Max error found: 6.084616e-20
		// Maximum Deviation Found (error term) 4.811e-20
		//
		static const float Y = 2.249481201171875f;
		static const double P[] = {
		-0.202433508355938759655,
		0.105264680699391713268,
		8.37050328343119927838,
		17.6447298408374015486,
		-18.8510648058714251895,
		-44.6382324441786960818,
		17.445385985570866523,
		21.1294655448340526258,
		-3.67192254707729348546
		};
		static const double Q[] = {
		1,
		6.24264124854247537712,
		3.9713437953343869095,
		-28.6608180499800029974,
		-20.1432634680485188801,
		48.5609213108739935468,
		10.8268667355460159008,
		-22.6436933413139721736,
		1.72114765761200282724
		};
		double g = sqrt(-2 * log(q));
		double xs = q - 0.25;
		double r = evaluate_polynomial(P, sizeof(P)/sizeof(double), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(double), xs);
		result = g / (Y + r);
		}
		else
		{
		//
		// For q < 0.25 we have a series of rational approximations all
		// of the general form:
		//
		// let: x = sqrt(-log(q))
		//
		// Then the result is given by:
		//
		// x(Y+R(x-B))
		//
		// where Y is a constant, B is the lowest value of x for which
		// the approximation is valid, and R(x-B) is optimised for a low
		// absolute error compared to Y.
		//
		// Note that almost all code will really go through the first
		// or maybe second approximation. After than we're dealing with very
		// small input values indeed: 80 and 128 bit long double's go all the
		// way down to ~ 1e-5000 so the "tail" is rather long...
		//
		double x = sqrt(-log(q));
		if(x < 3)
		{
		// Max error found: 1.089051e-20
		static const float Y = 0.807220458984375f;
		static const double P[] = {
		-0.131102781679951906451,
		-0.163794047193317060787,
		0.117030156341995252019,
		0.387079738972604337464,
		0.337785538912035898924,
		0.142869534408157156766,
		0.0290157910005329060432,
		0.00214558995388805277169,
		-0.679465575181126350155e-6,
		0.285225331782217055858e-7,
		-0.681149956853776992068e-9
		};
		static const double Q[] = {
		1,
		3.46625407242567245975,
		5.38168345707006855425,
		4.77846592945843778382,
		2.59301921623620271374,
		0.848854343457902036425,
		0.152264338295331783612,
		0.01105924229346489121
		};
		double xs = x - 1.125;
		double R = evaluate_polynomial(P, sizeof(P)/sizeof(double), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(double), xs);
		result = Y * x + R * x;
		}
		else if(x < 6)
		{
		// Max error found: 8.389174e-21
		static const float Y = 0.93995571136474609375f;
		static const double P[] = {
		-0.0350353787183177984712,
		-0.00222426529213447927281,
		0.0185573306514231072324,
		0.00950804701325919603619,
		0.00187123492819559223345,
		0.000157544617424960554631,
		0.460469890584317994083e-5,
		-0.230404776911882601748e-9,
		0.266339227425782031962e-11
		};
		static const double Q[] = {
		1,
		1.3653349817554063097,
		0.762059164553623404043,
		0.220091105764131249824,
		0.0341589143670947727934,
		0.00263861676657015992959,
		0.764675292302794483503e-4
		};
		double xs = x - 3;
		double R = evaluate_polynomial(P, sizeof(P)/sizeof(double), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(double), xs);
		result = Y * x + R * x;
		}
		else if(x < 18)
		{
		// Max error found: 1.481312e-19
		static const float Y = 0.98362827301025390625f;
		static const double P[] = {
		-0.0167431005076633737133,
		-0.00112951438745580278863,
		0.00105628862152492910091,
		0.000209386317487588078668,
		0.149624783758342370182e-4,
		0.449696789927706453732e-6,
		0.462596163522878599135e-8,
		-0.281128735628831791805e-13,
		0.99055709973310326855e-16
		};
		static const double Q[] = {
		1,
		0.591429344886417493481,
		0.138151865749083321638,
		0.0160746087093676504695,
		0.000964011807005165528527,
		0.275335474764726041141e-4,
		0.282243172016108031869e-6
		};
		double xs = x - 6;
		double R = evaluate_polynomial(P, sizeof(P)/sizeof(double), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(double), xs);
		result = Y * x + R * x;
		}
		else if(x < 44)
		{
		// Max error found: 5.697761e-20
		static const float Y = 0.99714565277099609375f;
		static const double P[] = {
		-0.0024978212791898131227,
		-0.779190719229053954292e-5,
		0.254723037413027451751e-4,
		0.162397777342510920873e-5,
		0.396341011304801168516e-7,
		0.411632831190944208473e-9,
		0.145596286718675035587e-11,
		-0.116765012397184275695e-17
		};
		static const double Q[] = {
		1,
		0.207123112214422517181,
		0.0169410838120975906478,
		0.000690538265622684595676,
		0.145007359818232637924e-4,
		0.144437756628144157666e-6,
		0.509761276599778486139e-9
		};
		double xs = x - 18;
		double R = evaluate_polynomial(P, sizeof(P)/sizeof(double), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(double), xs);
		result = Y * x + R * x;
		}
		else
		{
		// Max error found: 1.279746e-20
		static const float Y = 0.99941349029541015625f;
		static const double P[] = {
		-0.000539042911019078575891,
		-0.28398759004727721098e-6,
		0.899465114892291446442e-6,
		0.229345859265920864296e-7,
		0.225561444863500149219e-9,
		0.947846627503022684216e-12,
		0.135880130108924861008e-14,
		-0.348890393399948882918e-21
		};
		static const double Q[] = {
		1,
		0.0845746234001899436914,
		0.00282092984726264681981,
		0.468292921940894236786e-4,
		0.399968812193862100054e-6,
		0.161809290887904476097e-8,
		0.231558608310259605225e-11
		};
		double xs = x - 44;
		double R = evaluate_polynomial(P, sizeof(P)/sizeof(double), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(double), xs);
		result = Y * x + R * x;
		}
		}
		return result;
		}

		ccl_device double erf_inv(double z)
		{
		double p, q, s;

		if(z < 0)
		{
		p = -z;
		q = 1 - p;
		s = -1;
		}
		else
		{
		p = z;
		q = 1 - z;
		s = 1;
		}

		return s * erf_inv_impl(p, q);
		}

		ccl_device void beckmann_sample_slopes(
		// input
		const double theta_i, // incident direction
		double U1, double U2, // random numbers
		// output
		double& slope_x, double& slope_y // slopes
		)
		{
		// special case (normal incidence)
		if(theta_i < 0.0001)
		{
		const double r = sqrt(-log(U1));
		const double phi = 6.28318530718 * U2;
		slope_x = r * cos(phi);
		slope_y = r * sin(phi);
		return;
		}

		// precomputations
		const double sin_theta_i = sin(theta_i);
		const double cos_theta_i = cos(theta_i);
		const double tan_theta_i = sin_theta_i/cos_theta_i;
		const double a = 1.0 / tan_theta_i;
		const double erf_a = erf(a);
		const double exp_a2 = exp(-a*a);
		const double SQRT_PI_INV = 0.56418958354;
		const double Lambda = 0.5(erf_a-1) + 0.5SQRT_PI_INV*exp_a2/a;
		const double G1 = 1.0 / (1.0 + Lambda); // masking
		const double C= 1.0 - G1 * erf_a;

		// sample slope X
		if(U1 < C) {
		// rescale U1
		U1 = U1 / C;
		const double w_1 = 0.5 * SQRT_PI_INV * sin_theta_i * exp_a2;
		const double w_2 = cos_theta_i * (0.5 - 0.5*erf_a);
		const double p = w_1 / (w_1+w_2);

		if(U1 < p) {
		U1 = U1 / p;
		slope_x = -sqrt(-log(U1*exp_a2));
		}
		else {
		U1 = (U1-p) / (1.0-p);
		slope_x = erf_inv(U1-1.0-U1*erf_a);
		}
		}
		else {
		// rescale U1
		U1 = (U1-C) / (1.0-C);
		slope_x = erf_inv((-1.0+2.0U1)erf_a);

		const double p = (-slope_xsin_theta_i + cos_theta_i) / (2.0cos_theta_i);

		if (U2 > p) {
		slope_x = -slope_x;
		U2 = (U2-p) / (1.0-p);
		}
		else
		U2 = U2 / p;
		}

		// sample slope Y
		slope_y = erf_inv(2.0*U2-1.0);
		}

		ccl_device void ggx_sample_slopes(
		// input
		const double theta_i, // incident direction
		double U1, double U2, // random numbers
		// output
		double& slope_x, double& slope_y // slopes
		)
		{
		// special case (normal incidence)
		if(theta_i < 0.0001) {
		const double r = sqrt(U1/(1-U1));
		const double phi = 6.28318530718 * U2;
		slope_x = r * cos(phi);
		slope_y = r * sin(phi);
		return;
		}

		// precomputations
		const double tan_theta_i = tan(theta_i);
		const double a = 1 / (tan_theta_i);
		const double G1 = 2 / (1 + sqrt(1.0+1.0/(a*a)));
		// sample slope_x
		const double A = 2.0*U1/G1 - 1.0;
		const double tmp = 1.0 / (A*A-1.0);
		const double B = tan_theta_i;
		const double D = sqrt(BBtmptmp - (AA-BB)tmp);
		const double slope_x_1 = B*tmp - D;
		const double slope_x_2 = B*tmp + D;
		slope_x = (A < 0 \|\| slope_x_2 > 1.0/tan_theta_i) ? slope_x_1 : slope_x_2;

		// sample slope_y
		double S;

		if(U2 > 0.5) {
		S = 1.0;
		U2 = 2.0*(U2-0.5);
		}
		else {
		S = -1.0;
		U2 = 2.0*(0.5-U2);
		}

		const double z = (U2(U2(U20.27385-0.73369)+0.46341)) / (U2(U2(U20.093073+0.309420)-1.000000)+0.597999);
		slope_y = S * z * sqrt(1.0+slope_x*slope_x);
		}

		ccl_device void microfacet_sample_stretched(
		// input
		const float3 omega_i, // incident direction
		const double alpha_x, const double alpha_y, // anisotropic roughness
		const double U1, const double U2, // random numbers
		// output
		float3& omega_m, bool beckmann) // normal
		{
		// 1. stretch omega_i
		double omega_i_[3];
		omega_i_[0] = alpha_x * omega_i.x;
		omega_i_[1] = alpha_y * omega_i.y;
		omega_i_[2] = omega_i.z;

		// normalize
		double inv_omega_i = 1.0 / sqrt(omega_i_[0]omega_i_[0] + omega_i_[1]omega_i_[1] + omega_i_[2]*omega_i_[2]);
		omega_i_[0] *= inv_omega_i;
		omega_i_[1] *= inv_omega_i;
		omega_i_[2] *= inv_omega_i;

		// get polar coordinates of omega_i_
		double theta_ = 0.0;
		double phi_ = 0.0;
		if (omega_i_[2] < 0.99999)
		{
		theta_ = acos(omega_i_[2]);
		phi_ = atan2(omega_i_[1], omega_i_[0]);
		}

		// 2. sample P22_{omega_i}(x_slope, y_slope, 1, 1)
		double slope_x, slope_y;
		if(beckmann)
		beckmann_sample_slopes(theta_, U1, U2, slope_x, slope_y);
		else
		ggx_sample_slopes(theta_, U1, U2, slope_x, slope_y);

		// 3. rotate
		double tmp = cos(phi_)slope_x - sin(phi_)slope_y;
		slope_y = sin(phi_)slope_x + cos(phi_)slope_y;
		slope_x = tmp;

		// 4. unstretch
		slope_x = alpha_x * slope_x;
		slope_y = alpha_y * slope_y;

		// 5. compute normal
		double inv_omega_m = sqrt(slope_xslope_x + slope_yslope_y + 1.0);
		omega_m.x = (float)(-slope_x/inv_omega_m);
		omega_m.y = (float)(-slope_y/inv_omega_m);
		omega_m.z = (float)(1.0/inv_omega_m);
		}



/* GGX */		/* GGX */

ccl_device int bsdf_microfacet_ggx_setup(ShaderClosure *sc)		ccl_device int bsdf_microfacet_ggx_setup(ShaderClosure *sc)
{		{
sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ag */		sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ag */

sc->type = CLOSURE_BSDF_MICROFACET_GGX_ID;		sc->type = CLOSURE_BSDF_MICROFACET_GGX_ID;

Show All 35 Lines	if(cosNI > 0 && cosNO > 0) {
float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;		float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
float cosThetaM4 = cosThetaM2 * cosThetaM2;		float cosThetaM4 = cosThetaM2 * cosThetaM2;
float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));		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 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 G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI)));
float G = G1o * G1i;		float G = G1o * G1i;
float out = (G * D) * 0.25f / cosNO;		float out = (G * D) * 0.25f / cosNO;

		#if 0
// eq. 24		// eq. 24
float pm = D * cosThetaM;		float pm = D * cosThetaM;
		#else
		// eq. 2 in Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals
		float Dw = G1o * dot(I, Hr) * D / cosNO;
		float pm = Dw;
		#endif

// convert into pdf of the sampled direction		// convert into pdf of the sampled direction
// eq. 38 - but see also:		// eq. 38 - but see also:
// eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf		// eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
pdf = pm 0.25f / dot(Hr, I);		pdf = pm 0.25f / dot(Hr, I);
return make_float3 (out, out, out);		return make_float3 (out, out, out);
}		}
return make_float3 (0, 0, 0);		return make_float3 (0, 0, 0);
}		}
▲ Show 20 Lines • Show All 45 Lines • ▼ Show 20 Lines	ccl_device int bsdf_microfacet_ggx_sample(const ShaderClosure sc, float3 Ng, float3 I, float3 dIdx, float3 dIdy, float randu, float randv, float3 eval, float3 omega_in, float3 domega_in_dx, float3 domega_in_dy, float pdf)
if(cosNO > 0) {		if(cosNO > 0) {
float3 X, Y, Z = N;		float3 X, Y, Z = N;
make_orthonormals(Z, &X, &Y);		make_orthonormals(Z, &X, &Y);
// generate a random microfacet normal m		// generate a random microfacet normal m
// eq. 35,36:		// eq. 35,36:
// we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)		// we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)
//tttt and sin(atan(x)) == x/sqrt(1+x^2)		//tttt and sin(atan(x)) == x/sqrt(1+x^2)
float alpha2 = m_ag * m_ag;		float alpha2 = m_ag * m_ag;

		#if 0
float tanThetaM2 = alpha2 * randu / (1 - randu);		float tanThetaM2 = alpha2 * randu / (1 - randu);
float cosThetaM = 1 / safe_sqrtf(1 + tanThetaM2);		float cosThetaM = 1 / safe_sqrtf(1 + tanThetaM2);
float sinThetaM = cosThetaM * safe_sqrtf(tanThetaM2);		float sinThetaM = cosThetaM * safe_sqrtf(tanThetaM2);
float phiM = M_2PI_F * randv;		float phiM = M_2PI_F * randv;
float3 m = (cosf(phiM) * sinThetaM) * X +		float3 m = (cosf(phiM) * sinThetaM) * X +
(sinf(phiM) * sinThetaM) * Y +		(sinf(phiM) * sinThetaM) * Y +
( cosThetaM) * Z;		( cosThetaM) * Z;
		#else
		float3 local_I, local_m;

		local_I.x = dot(X, I);
		local_I.y = dot(Y, I);
		local_I.z = dot(Z, I);

		microfacet_sample_stretched(local_I, m_ag, m_ag, randu, randv, local_m, false);

		float3 m = Xlocal_m.x + Ylocal_m.y + Z*local_m.z;
		float cosThetaM = dot(m, N);
		float tanThetaM2 = 1/(cosThetaM*cosThetaM) - 1;
		#endif

if(!m_refractive) {		if(!m_refractive) {
float cosMO = dot(m, I);		float cosMO = dot(m, I);
if(cosMO > 0) {		if(cosMO > 0) {
// eq. 39 - compute actual reflected direction		// eq. 39 - compute actual reflected direction
omega_in = 2 cosMO * m - I;		omega_in = 2 cosMO * m - I;
if(dot(Ng, *omega_in) > 0) {		if(dot(Ng, *omega_in) > 0) {
if (m_ag <= 1e-4f) {		if (m_ag <= 1e-4f) {
// some high number for MIS		// some high number for MIS
*pdf = 1e6f;		*pdf = 1e6f;
*eval = make_float3(1e6f, 1e6f, 1e6f);		*eval = make_float3(1e6f, 1e6f, 1e6f);
}		}
else {		else {
// microfacet normal is visible to this ray		// microfacet normal is visible to this ray
// eq. 33		// eq. 33
float cosThetaM2 = cosThetaM * cosThetaM;		float cosThetaM2 = cosThetaM * cosThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2;		float cosThetaM4 = cosThetaM2 * cosThetaM2;
float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));		float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
// 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 / cosMO;
// eval BRDF*cosNI		// eval BRDF*cosNI
float cosNI = dot(N, *omega_in);		float cosNI = dot(N, *omega_in);
// 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 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 G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI)));
float G = G1o * G1i;		float G = G1o * G1i;

		#if 0
		// eq. 24
		float pm = D * cosThetaM;
		#else
		// eq. 2 in Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals
		float Dw = G1o * dot(I, m) * D / cosNO;
		float pm = Dw;
		#endif

		// 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 / cosMO;
// eq. 20: (FGD)/(4inon)		// eq. 20: (FGD)/(4inon)
float out = (G * D) * 0.25f / cosNO;		float out = (G * D) * 0.25f / cosNO;
*eval = make_float3(out, out, out);		*eval = make_float3(out, out, out);
}		}

#ifdef __RAY_DIFFERENTIALS__		#ifdef __RAY_DIFFERENTIALS__
domega_in_dx = (2 dot(m, dIdx)) * m - dIdx;		domega_in_dx = (2 dot(m, dIdx)) * m - dIdx;
domega_in_dy = (2 dot(m, dIdy)) * m - dIdy;		domega_in_dy = (2 dot(m, dIdy)) * m - dIdy;
Show All 28 Lines	#endif
*pdf = 1e6f;		*pdf = 1e6f;
*eval = make_float3(1e6f, 1e6f, 1e6f);		*eval = make_float3(1e6f, 1e6f, 1e6f);
}		}
else {		else {
// eq. 33		// eq. 33
float cosThetaM2 = cosThetaM * cosThetaM;		float cosThetaM2 = cosThetaM * cosThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2;		float cosThetaM4 = cosThetaM2 * cosThetaM2;
float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));		float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
// eq. 24
float pm = D * cosThetaM;
// eval BRDF*cosNI		// eval BRDF*cosNI
float cosNI = dot(N, *omega_in);		float cosNI = dot(N, *omega_in);
// 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 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 G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI)));
float G = G1o * G1i;		float G = G1o * G1i;

		#if 0
		// eq. 24
		float pm = D * cosThetaM;
		#else
		// eq. 2 in Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals
		float Dw = G1o * dot(I, m) * D / cosNO;
		float pm = Dw;
		#endif

// eq. 21		// eq. 21
float cosHI = dot(m, *omega_in);		float cosHI = dot(m, *omega_in);
float cosHO = dot(m, I);		float cosHO = dot(m, I);
float Ht2 = m_eta * cosHI + cosHO;		float Ht2 = m_eta * cosHI + cosHO;
Ht2 *= Ht2;		Ht2 *= Ht2;
float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D)) / (cosNO * Ht2);		float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D)) / (cosNO * Ht2);
// eq. 38 and eq. 17		// eq. 38 and eq. 17
pdf = pm (m_eta * m_eta) * fabsf(cosHI) / Ht2;		pdf = pm (m_eta * m_eta) * fabsf(cosHI) / Ht2;
▲ Show 20 Lines • Show All 51 Lines • ▼ Show 20 Lines	if(cosNO > 0 && cosNI > 0) {
float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 * cosThetaM4);		float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 * cosThetaM4);
// eq. 26, 27: now calculate G1(i,m) and G1(o,m)		// eq. 26, 27: now calculate G1(i,m) and G1(o,m)
float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));		float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));		float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;		float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;		float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
float G = G1o * G1i;		float G = G1o * G1i;
float out = (G * D) * 0.25f / cosNO;		float out = (G * D) * 0.25f / cosNO;

		#if 0
// eq. 24		// eq. 24
float pm = D * cosThetaM;		float pm = D * cosThetaM;
		#else
		// eq. 2 in Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals
		float Dw = G1o * dot(I, Hr) * D / cosNO;
		float pm = Dw;
		#endif

// convert into pdf of the sampled direction		// convert into pdf of the sampled direction
// eq. 38 - but see also:		// eq. 38 - but see also:
// eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf		// eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
pdf = pm 0.25f / dot(Hr, I);		pdf = pm 0.25f / dot(Hr, I);
return make_float3 (out, out, out);		return make_float3 (out, out, out);
}		}
return make_float3 (0, 0, 0);		return make_float3 (0, 0, 0);
}		}
▲ Show 20 Lines • Show All 42 Lines • ▼ Show 20 Lines	ccl_device int bsdf_microfacet_beckmann_sample(const ShaderClosure sc, float3 Ng, float3 I, float3 dIdx, float3 dIdy, float randu, float randv, float3 eval, float3 omega_in, float3 domega_in_dx, float3 domega_in_dy, float pdf)
float m_ab = sc->data0;		float m_ab = sc->data0;
int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;		int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
float3 N = sc->N;		float3 N = sc->N;

float cosNO = dot(N, I);		float cosNO = dot(N, I);
if(cosNO > 0) {		if(cosNO > 0) {
float3 X, Y, Z = N;		float3 X, Y, Z = N;
make_orthonormals(Z, &X, &Y);		make_orthonormals(Z, &X, &Y);

// generate a random microfacet normal m		// generate a random microfacet normal m
// eq. 35,36:		// eq. 35,36:
// we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)		// we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)
//tttt and sin(atan(x)) == x/sqrt(1+x^2)		//tttt and sin(atan(x)) == x/sqrt(1+x^2)
float alpha2 = m_ab * m_ab;		float alpha2 = m_ab * m_ab;

		#if 0
float tanThetaM, cosThetaM;		float tanThetaM, cosThetaM;

if(alpha2 == 0.0f) {		if(alpha2 == 0.0f) {
tanThetaM = 0.0f;		tanThetaM = 0.0f;
cosThetaM = 1.0f;		cosThetaM = 1.0f;
}		}
else {		else {
tanThetaM = safe_sqrtf(-alpha2 * logf(1 - randu));		tanThetaM = safe_sqrtf(-alpha2 * logf(1 - randu));
cosThetaM = 1 / safe_sqrtf(1 + tanThetaM * tanThetaM);		cosThetaM = 1 / safe_sqrtf(1 + tanThetaM * tanThetaM);
}		}

float sinThetaM = cosThetaM * tanThetaM;		float sinThetaM = cosThetaM * tanThetaM;
float phiM = M_2PI_F * randv;		float phiM = M_2PI_F * randv;
float3 m = (cosf(phiM) * sinThetaM) * X +		float3 m = (cosf(phiM) * sinThetaM) * X +
(sinf(phiM) * sinThetaM) * Y +		(sinf(phiM) * sinThetaM) * Y +
( cosThetaM) * Z;		( cosThetaM) * Z;
		#else
		float3 local_I, local_m;

		local_I.x = dot(X, I);
		local_I.y = dot(Y, I);
		local_I.z = dot(Z, I);

		microfacet_sample_stretched(local_I, m_ab, m_ab, randu, randv, local_m, true);

		float3 m = Xlocal_m.x + Ylocal_m.y + Z*local_m.z;
		float cosThetaM = dot(m, N);
		float tanThetaM = safe_sqrtf(1/(cosThetaM*cosThetaM) - 1);
		#endif

if(!m_refractive) {		if(!m_refractive) {
float cosMO = dot(m, I);		float cosMO = dot(m, I);
if(cosMO > 0) {		if(cosMO > 0) {
// eq. 39 - compute actual reflected direction		// eq. 39 - compute actual reflected direction
omega_in = 2 cosMO * m - I;		omega_in = 2 cosMO * m - I;
if(dot(Ng, *omega_in) > 0) {		if(dot(Ng, *omega_in) > 0) {
if (m_ab <= 1e-4f) {		if (m_ab <= 1e-4f) {
// some high number for MIS		// some high number for MIS
*pdf = 1e6f;		*pdf = 1e6f;
*eval = make_float3(1e6f, 1e6f, 1e6f);		*eval = make_float3(1e6f, 1e6f, 1e6f);
}		}
else {		else {



// microfacet normal is visible to this ray		// microfacet normal is visible to this ray
// eq. 25		// eq. 25
float cosThetaM2 = cosThetaM * cosThetaM;		float cosThetaM2 = cosThetaM * cosThetaM;
float tanThetaM2 = tanThetaM * tanThetaM;		float tanThetaM2 = tanThetaM * tanThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2;		float cosThetaM4 = cosThetaM2 * cosThetaM2;
float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 * cosThetaM4);		float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 * cosThetaM4);
// 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 / cosMO;
// Eval BRDF*cosNI		// Eval BRDF*cosNI
float cosNI = dot(N, *omega_in);		float cosNI = dot(N, *omega_in);
// eq. 26, 27: now calculate G1(i,m) and G1(o,m)		// eq. 26, 27: now calculate G1(i,m) and G1(o,m)
float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));		float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));		float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;		float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;		float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
float G = G1o * G1i;		float G = G1o * G1i;

		#if 0
		// eq. 24
		float pm = D * cosThetaM;
		#else
		// eq. 2 in Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals
		float Dw = G1o * dot(I, m) * D / cosNO;
		float pm = Dw;
		#endif

		// 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 / cosMO;
// eq. 20: (FGD)/(4inon)		// eq. 20: (FGD)/(4inon)
float out = (G * D) * 0.25f / cosNO;		float out = (G * D) * 0.25f / cosNO;
*eval = make_float3(out, out, out);		*eval = make_float3(out, out, out);
}		}
#ifdef __RAY_DIFFERENTIALS__		#ifdef __RAY_DIFFERENTIALS__
domega_in_dx = (2 dot(m, dIdx)) * m - dIdx;		domega_in_dx = (2 dot(m, dIdx)) * m - dIdx;
domega_in_dy = (2 dot(m, dIdy)) * m - dIdy;		domega_in_dy = (2 dot(m, dIdy)) * m - dIdy;
#endif		#endif
Show All 27 Lines	#endif
*eval = make_float3(1e6f, 1e6f, 1e6f);		*eval = make_float3(1e6f, 1e6f, 1e6f);
}		}
else {		else {
// eq. 33		// eq. 33
float cosThetaM2 = cosThetaM * cosThetaM;		float cosThetaM2 = cosThetaM * cosThetaM;
float tanThetaM2 = tanThetaM * tanThetaM;		float tanThetaM2 = tanThetaM * tanThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2;		float cosThetaM4 = cosThetaM2 * cosThetaM2;
float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 * cosThetaM4);		float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 * cosThetaM4);
// eq. 24
float pm = D * cosThetaM;
// eval BRDF*cosNI		// eval BRDF*cosNI
float cosNI = dot(N, *omega_in);		float cosNI = dot(N, *omega_in);
// eq. 26, 27: now calculate G1(i,m) and G1(o,m)		// eq. 26, 27: now calculate G1(i,m) and G1(o,m)
float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));		float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));		float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;		float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;		float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
float G = G1o * G1i;		float G = G1o * G1i;

		#if 0
		// eq. 24
		float pm = D * cosThetaM;
		#else
		// eq. 2 in Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals
		float Dw = G1o * dot(I, m) * D / cosNO;
		float pm = Dw;
		#endif

// eq. 21		// eq. 21
float cosHI = dot(m, *omega_in);		float cosHI = dot(m, *omega_in);
float cosHO = dot(m, I);		float cosHO = dot(m, I);
float Ht2 = m_eta * cosHI + cosHO;		float Ht2 = m_eta * cosHI + cosHO;
Ht2 *= Ht2;		Ht2 *= Ht2;
float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D)) / (cosNO * Ht2);		float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D)) / (cosNO * Ht2);
// eq. 38 and eq. 17		// eq. 38 and eq. 17
pdf = pm (m_eta * m_eta) * fabsf(cosHI) / Ht2;		pdf = pm (m_eta * m_eta) * fabsf(cosHI) / Ht2;
Show All 12 Lines