Differential D572 Diff 1820 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

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

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

				return sum;
				}

				ccl_device float erf_inv_impl(const float& p, const float& q)
				{
				float 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\|.
				//
				// float: Max error found: 2.001849e-18
				// long float: 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 float P[] = {
				-0.000508781949658280665617,
				-0.00836874819741736770379,
				0.0334806625409744615033,
				-0.0126926147662974029034,
				-0.0365637971411762664006,
				0.0219878681111168899165,
				0.00822687874676915743155,
				-0.00538772965071242932965
				};
				static const float Q[] = {
				1,
				-0.970005043303290640362,
				-1.56574558234175846809,
				1.56221558398423026363,
				0.662328840472002992063,
				-0.71228902341542847553,
				-0.0527396382340099713954,
				0.0795283687341571680018,
				-0.00233393759374190016776,
				0.000886216390456424707504
				};
				float g = p * (p + 10);
				float r = evaluate_polynomial(P, sizeof(P)/sizeof(float), p) / evaluate_polynomial(Q, sizeof(Q)/sizeof(float), p);
				result = g * Y + g * r;
				}
				else if(q >= 0.25)
				{
				//
				// Rational approximation for 0.5 > q >= 0.25
				//
				// x = sqrtf(-2*log(q)) / (Y + R(q))
				//
				// Where Y is a constant, and R(q) is optimised for a low
				// absolute error compared to Y.
				//
				// float : Max error found: 7.403372e-17
				// long float : Max error found: 6.084616e-20
				// Maximum Deviation Found (error term) 4.811e-20
				//
				static const float Y = 2.249481201171875f;
				static const float P[] = {
				-0.202433508355938759655,
				0.105264680699391713268,
				8.37050328343119927838,
				17.6447298408374015486,
				-18.8510648058714251895,
				-44.6382324441786960818,
				17.445385985570866523,
				21.1294655448340526258,
				-3.67192254707729348546
				};
				static const float Q[] = {
				1,
				6.24264124854247537712,
				3.9713437953343869095,
				-28.6608180499800029974,
				-20.1432634680485188801,
				48.5609213108739935468,
				10.8268667355460159008,
				-22.6436933413139721736,
				1.72114765761200282724
				};
				float g = sqrtf(-2 * log(q));
				float xs = q - 0.25;
				float r = evaluate_polynomial(P, sizeof(P)/sizeof(float), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(float), xs);
				result = g / (Y + r);
				}
				else
				{
				//
				// For q < 0.25 we have a series of rational approximations all
				// of the general form:
				//
				// let: x = sqrtf(-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 float's go all the
				// way down to ~ 1e-5000 so the "tail" is rather long...
				//
				float x = sqrtf(-log(q));
				if(x < 3)
				{
				// Max error found: 1.089051e-20
				static const float Y = 0.807220458984375f;
				static const float 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 float Q[] = {
				1,
				3.46625407242567245975,
				5.38168345707006855425,
				4.77846592945843778382,
				2.59301921623620271374,
				0.848854343457902036425,
				0.152264338295331783612,
				0.01105924229346489121
				};
				float xs = x - 1.125;
				float R = evaluate_polynomial(P, sizeof(P)/sizeof(float), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(float), xs);
				result = Y * x + R * x;
				}
				else if(x < 6)
				{
				// Max error found: 8.389174e-21
				static const float Y = 0.93995571136474609375f;
				static const float P[] = {
				-0.0350353787183177984712,
				-0.00222426529213447927281,
				0.0185573306514231072324,
				0.00950804701325919603619,
				0.00187123492819559223345,
				0.000157544617424960554631,
				0.460469890584317994083e-5,
				-0.230404776911882601748e-9,
				0.266339227425782031962e-11
				};
				static const float Q[] = {
				1,
				1.3653349817554063097,
				0.762059164553623404043,
				0.220091105764131249824,
				0.0341589143670947727934,
				0.00263861676657015992959,
				0.764675292302794483503e-4
				};
				float xs = x - 3;
				float R = evaluate_polynomial(P, sizeof(P)/sizeof(float), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(float), xs);
				result = Y * x + R * x;
				}
				else if(x < 18)
				{
				// Max error found: 1.481312e-19
				static const float Y = 0.98362827301025390625f;
				static const float 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 float Q[] = {
				1,
				0.591429344886417493481,
				0.138151865749083321638,
				0.0160746087093676504695,
				0.000964011807005165528527,
				0.275335474764726041141e-4,
				0.282243172016108031869e-6
				};
				float xs = x - 6;
				float R = evaluate_polynomial(P, sizeof(P)/sizeof(float), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(float), xs);
				result = Y * x + R * x;
				}
				else if(x < 44)
				{
				// Max error found: 5.697761e-20
				static const float Y = 0.99714565277099609375f;
				static const float 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 float Q[] = {
				1,
				0.207123112214422517181,
				0.0169410838120975906478,
				0.000690538265622684595676,
				0.145007359818232637924e-4,
				0.144437756628144157666e-6,
				0.509761276599778486139e-9
				};
				float xs = x - 18;
				float R = evaluate_polynomial(P, sizeof(P)/sizeof(float), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(float), xs);
				result = Y * x + R * x;
				}
				else
				{
				// Max error found: 1.279746e-20
				static const float Y = 0.99941349029541015625f;
				static const float 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 float Q[] = {
				1,
				0.0845746234001899436914,
				0.00282092984726264681981,
				0.468292921940894236786e-4,
				0.399968812193862100054e-6,
				0.161809290887904476097e-8,
				0.231558608310259605225e-11
				};
				float xs = x - 44;
				float R = evaluate_polynomial(P, sizeof(P)/sizeof(float), xs) / evaluate_polynomial(Q, sizeof(Q)/sizeof(float), xs);
				result = Y * x + R * x;
				}
				}
				return result;
				}

				ccl_device float erf_inv(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 * erf_inv_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,
				float randu, float randv, float slope_x, float slope_y)
				{
				/* 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);
				return;
				}

				/* precomputations */
				const float inv_tan_theta_i = cos_theta_i/sin_theta_i;
				const float a = inv_tan_theta_i;
				const float erf_a = erf(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.5fSQRT_PI_INV*exp_a2/a;
				const float G1 = 1.0f / (1.0f + Lambda); /* masking */
				const float C= 1.0f - G1 * erf_a;

				/* 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(randuexp_a2));
				}
				else {
				randu = (randu-p) / (1.0f - p);
				slope_x = erf_inv(randu - 1.0f - randuerf_a);
				}
				}
				else {
				/* rescale randu */
				randu = (randu - C) / (1.0f - C);
				slope_x = erf_inv((-1.0f + 2.0frandu)*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 = erf_inv(2.0frandv - 1.0f);
				}

				ccl_device_inline void microfacet_ggx_sample_slopes(
				const float cos_theta_i, const float sin_theta_i,
				float randu, float randv, float slope_x, float slope_y)
				{
				/* 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);
				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));

				/* sample slope_x */
				const float A = 2.0franduG1_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(tmptmp) - (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_2tan_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(randv0.27385f - 0.73369f) + 0.46341f)) / (randv(randv(randv0.093073f + 0.309420f) - 1.000000f) + 0.597999f);
				slope_y = S z * safe_sqrtf(1.0f + (slope_x)(*slope_x));
				}

				ccl_device void microfacet_sample_stretched(const float3 omega_i,
				const float alpha_x, const float alpha_y,
				const float randu, const float randv,
				bool beckmann, float3 *omega_m)
				{
				/* 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/hypotf(omega_i_.x, omega_i_.y);
				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_,
				randu, randv, &slope_x, &slope_y);
				else
				microfacet_ggx_sample_slopes(costheta_, sintheta_,
				randu, randv, &slope_x, &slope_y);

				/* 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 */
				float inv_omega_m = safe_sqrtf(slope_xslope_x + slope_yslope_y + 1.0f);
				*omega_m = make_float3(-slope_x, -slope_y, 1.0f) / inv_omega_m;
				}

				/* 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)			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;

	return SD_BSDF\|SD_BSDF_HAS_EVAL\|SD_BSDF_GLOSSY;			return SD_BSDF\|SD_BSDF_HAS_EVAL\|SD_BSDF_GLOSSY;
	Show All 15 Lines

	ccl_device float3 bsdf_microfacet_ggx_eval_reflect(const ShaderClosure sc, const float3 I, const float3 omega_in, float pdf)			ccl_device float3 bsdf_microfacet_ggx_eval_reflect(const ShaderClosure sc, const float3 I, const float3 omega_in, float pdf)
	{			{
	float m_ag = max(sc->data0, 1e-4f);			float m_ag = max(sc->data0, 1e-4f);
	int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;			int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
	float3 N = sc->N;			float3 N = sc->N;

	if(m_refractive \|\| m_ag <= 1e-4f)			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 cosNO = dot(N, I);
	float cosNI = dot(N, omega_in);			float cosNI = dot(N, omega_in);

	if(cosNI > 0 && cosNO > 0) {			if(cosNI > 0 && cosNO > 0) {
	// get half vector			/* get half vector */
	float3 Hr = normalize(omega_in + I);			float3 m = normalize(omega_in + I);
	// eq. 20: (FGD)/(4inon)
	// eq. 33: first we calculate D(m) with m=Hr:			/* eq. 20: (FGD)/(4inon)
				* eq. 33: first we calculate D(m) */
	float alpha2 = m_ag * m_ag;			float alpha2 = m_ag * m_ag;
	float cosThetaM = dot(N, Hr);			float cosThetaM = dot(N, m);
	float cosThetaM2 = cosThetaM * cosThetaM;			float cosThetaM2 = cosThetaM * cosThetaM;
	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;
	// eq. 24			/* eq. 20 */
	float pm = D * cosThetaM;			float common = D * 0.25f / cosNO;
	// convert into pdf of the sampled direction			float out = G * common;
	// eq. 38 - but see also:
	// eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf			/* eq. 2 in distribution of visible normals sampling
	pdf = pm 0.25f / dot(Hr, I);			* pm = Dw = G1o * dot(m, I) * D / dot(N, I); */

				/* eq. 38 - but see also:
				* eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
				* pdf = pm * 0.25 / dot(m, I); */
				pdf = G1o common;

	return make_float3 (out, out, out);			return make_float3(out, out, out);
	}			}

	return make_float3 (0, 0, 0);			return make_float3(0, 0, 0);
	}			}

	ccl_device float3 bsdf_microfacet_ggx_eval_transmit(const ShaderClosure sc, const float3 I, const float3 omega_in, float pdf)			ccl_device float3 bsdf_microfacet_ggx_eval_transmit(const ShaderClosure sc, const float3 I, const float3 omega_in, float pdf)
	{			{
	float m_ag = max(sc->data0, 1e-4f);			float m_ag = max(sc->data0, 1e-4f);
	float m_eta = sc->data1;			float m_eta = sc->data1;
	int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;			int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
	float3 N = sc->N;			float3 N = sc->N;

	if(!m_refractive \|\| m_ag <= 1e-4f)			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 cosNO = dot(N, I);
	float cosNI = dot(N, omega_in);			float cosNI = dot(N, omega_in);

	if(cosNO <= 0 \|\| cosNI >= 0)			if(cosNO <= 0 \|\| cosNI >= 0)
	return make_float3 (0, 0, 0); // vectors on same side -- not possible			return make_float3(0, 0, 0); /* vectors on same side -- not possible */
	// compute half-vector of the refraction (eq. 16)
				/* compute half-vector of the refraction (eq. 16) */
	float3 ht = -(m_eta * omega_in + I);			float3 ht = -(m_eta * omega_in + I);
	float3 Ht = normalize(ht);			float3 Ht = normalize(ht);
	float cosHO = dot(Ht, I);			float cosHO = dot(Ht, I);

	float cosHI = dot(Ht, omega_in);			float cosHI = dot(Ht, omega_in);
	// eq. 33: first we calculate D(m) with m=Ht:
				/* eq. 33: first we calculate D(m) with m=Ht: */
	float alpha2 = m_ag * m_ag;			float alpha2 = m_ag * m_ag;
	float cosThetaM = dot(N, Ht);			float cosThetaM = dot(N, Ht);
	float cosThetaM2 = cosThetaM * cosThetaM;			float cosThetaM2 = cosThetaM * cosThetaM;
	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;
	// probability
	float invHt2 = 1 / dot(ht, ht);			/* probability */
	pdf = D fabsf(cosThetaM) * (fabsf(cosHI) * (m_eta * m_eta)) * invHt2;			float Ht2 = dot(ht, ht);
	float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D) * invHt2) / cosNO;
				/* eq. 2 in distribution of visible normals sampling
				* pm = Dw = G1o * dot(m, I) * D / dot(N, I); */

				/* out = fabsf(cosHI * cosHO) * (m_eta * m_eta) * G * D / (cosNO * Ht2)
				* pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2 */
				float common = D * (m_eta * m_eta) / (cosNO * Ht2);
				float out = G * fabsf(cosHI * cosHO) * common;
				pdf = G1o cosHO * fabsf(cosHI) * common;

	return make_float3 (out, out, out);			return make_float3(out, out, out);
	}			}

	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)			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)
	{			{
	float m_ag = sc->data0;			float m_ag = sc->data0;
	int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;			int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_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
	// eq. 35,36:			/* importance sampling with distribution of visible normals. vectors are
	// we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)			* transformed to local space before and after */
	//tttt and sin(atan(x)) == x/sqrt(1+x^2)			float3 local_I, local_m;
	float alpha2 = m_ag * m_ag;
	float tanThetaM2 = alpha2 * randu / (1 - randu);			local_I.x = dot(X, I);
	float cosThetaM = 1 / safe_sqrtf(1 + tanThetaM2);			local_I.y = dot(Y, I);
	float sinThetaM = cosThetaM * safe_sqrtf(tanThetaM2);			local_I.z = dot(Z, I);
	float phiM = M_2PI_F * randv;
	float3 m = (cosf(phiM) * sinThetaM) * X +			microfacet_sample_stretched(local_I, m_ag, m_ag, randu, randv, false, &local_m);
	(sinf(phiM) * sinThetaM) * Y +
	( cosThetaM) * Z;			float3 m = Xlocal_m.x + Ylocal_m.y + Z*local_m.z;
				float cosThetaM = dot(m, N);
				float tanThetaM2 = 1/(cosThetaM*cosThetaM) - 1;

				/* reflection or refraction? */
	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 alpha2 = m_ag * m_ag;
	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 BRDFcosNI /
	// 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
	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;
	// eq. 20: (FGD)/(4inon)
	float out = (G * D) * 0.25f / cosNO;			/* see eval function for derivation */
				float common = D * 0.25f / cosNO;
				float out = G * common;
				pdf = G1o common;

	*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
	}			}
	}			}
	}			}
	else {			else {
	// CAUTION: the i and o variables are inverted relative to the paper			/* CAUTION: the i and o variables are inverted relative to the paper
	// eq. 39 - compute actual refractive direction			* eq. 39 - compute actual refractive direction */
	float3 R, T;			float3 R, T;
	#ifdef __RAY_DIFFERENTIALS__			#ifdef __RAY_DIFFERENTIALS__
	float3 dRdx, dRdy, dTdx, dTdy;			float3 dRdx, dRdy, dTdx, dTdy;
	#endif			#endif
	float m_eta = sc->data1;			float m_eta = sc->data1;
	bool inside;			bool inside;

	fresnel_dielectric(m_eta, m, I, &R, &T,			fresnel_dielectric(m_eta, m, I, &R, &T,
	#ifdef __RAY_DIFFERENTIALS__			#ifdef __RAY_DIFFERENTIALS__
	dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,			dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,
	#endif			#endif
	&inside);			&inside);

	if(!inside) {			if(!inside) {

	*omega_in = T;			*omega_in = T;
	#ifdef __RAY_DIFFERENTIALS__			#ifdef __RAY_DIFFERENTIALS__
	*domega_in_dx = dTdx;			*domega_in_dx = dTdx;
	*domega_in_dy = dTdy;			*domega_in_dy = dTdy;
	#endif			#endif

	if (m_ag <= 1e-4f \|\| fabsf(m_eta - 1.0f) < 1e-4f) {			if(m_ag <= 1e-4f \|\| fabsf(m_eta - 1.0f) < 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 {
	// eq. 33			/* eq. 33 */
				float alpha2 = m_ag * m_ag;
	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 BRDFcosNI /
	// 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;
	// 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);
	// eq. 38 and eq. 17			/* see eval function for derivation */
	pdf = pm (m_eta * m_eta) * fabsf(cosHI) / Ht2;			float common = D * (m_eta * m_eta) / (cosNO * Ht2);
				float out = G * fabsf(cosHI * cosHO) * common;
				pdf = G1o cosHO * fabsf(cosHI) * common;

	*eval = make_float3(out, out, out);			*eval = make_float3(out, out, out);
	}			}
	}			}
	}			}
	}			}
	return (m_refractive) ? LABEL_TRANSMIT\|LABEL_GLOSSY : LABEL_REFLECT\|LABEL_GLOSSY;			return (m_refractive) ? LABEL_TRANSMIT\|LABEL_GLOSSY : LABEL_REFLECT\|LABEL_GLOSSY;
	}			}

	/* BECKMANN */			/* Beckmann 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_beckmann_setup(ShaderClosure *sc)			ccl_device int bsdf_microfacet_beckmann_setup(ShaderClosure *sc)
	{			{
	sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ab */			sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ab */

	sc->type = CLOSURE_BSDF_MICROFACET_BECKMANN_ID;			sc->type = CLOSURE_BSDF_MICROFACET_BECKMANN_ID;
	return SD_BSDF\|SD_BSDF_HAS_EVAL\|SD_BSDF_GLOSSY;			return SD_BSDF\|SD_BSDF_HAS_EVAL\|SD_BSDF_GLOSSY;
	}			}
	Show All 13 Lines

	ccl_device float3 bsdf_microfacet_beckmann_eval_reflect(const ShaderClosure sc, const float3 I, const float3 omega_in, float pdf)			ccl_device float3 bsdf_microfacet_beckmann_eval_reflect(const ShaderClosure sc, const float3 I, const float3 omega_in, float pdf)
	{			{
	float m_ab = max(sc->data0, 1e-4f);			float m_ab = max(sc->data0, 1e-4f);
	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;

	if(m_refractive \|\| m_ab <= 1e-4f)			if(m_refractive \|\| m_ab <= 1e-4f)
	return make_float3 (0, 0, 0);			return make_float3(0, 0, 0);

	float cosNO = dot(N, I);			float cosNO = dot(N, I);
	float cosNI = dot(N, omega_in);			float cosNI = dot(N, omega_in);

	if(cosNO > 0 && cosNI > 0) {			if(cosNO > 0 && cosNI > 0) {
	// get half vector			/* get half vector */
	float3 Hr = normalize(omega_in + I);			float3 m = normalize(omega_in + I);
	// eq. 20: (FGD)/(4inon)
	// eq. 25: first we calculate D(m) with m=Hr:			/* eq. 20: (FGD)/(4inon)
				* eq. 25: first we calculate D(m) */
	float alpha2 = m_ab * m_ab;			float alpha2 = m_ab * m_ab;
	float cosThetaM = dot(N, Hr);			float cosThetaM = dot(N, m);
	float cosThetaM2 = cosThetaM * cosThetaM;			float cosThetaM2 = cosThetaM * cosThetaM;
	float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;			float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
	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. 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;
	// eq. 24			/* eq. 20 */
	float pm = D * cosThetaM;			float common = D * 0.25f / cosNO;
	// convert into pdf of the sampled direction			float out = G * common;
	// eq. 38 - but see also:
	// eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf			/* eq. 2 in distribution of visible normals sampling
	pdf = pm 0.25f / dot(Hr, I);			* pm = Dw = G1o * dot(m, I) * D / dot(N, I); */

				/* eq. 38 - but see also:
				* eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
				* pdf = pm * 0.25 / dot(m, I); */
				pdf = G1o common;

	return make_float3 (out, out, out);			return make_float3(out, out, out);
	}			}

	return make_float3 (0, 0, 0);			return make_float3(0, 0, 0);
	}			}

	ccl_device float3 bsdf_microfacet_beckmann_eval_transmit(const ShaderClosure sc, const float3 I, const float3 omega_in, float pdf)			ccl_device float3 bsdf_microfacet_beckmann_eval_transmit(const ShaderClosure sc, const float3 I, const float3 omega_in, float pdf)
	{			{
	float m_ab = max(sc->data0, 1e-4f);			float m_ab = max(sc->data0, 1e-4f);
	float m_eta = sc->data1;			float m_eta = sc->data1;
	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;

	if(!m_refractive \|\| m_ab <= 1e-4f)			if(!m_refractive \|\| m_ab <= 1e-4f)
	return make_float3 (0, 0, 0);			return make_float3(0, 0, 0);

	float cosNO = dot(N, I);			float cosNO = dot(N, I);
	float cosNI = dot(N, omega_in);			float cosNI = dot(N, omega_in);

	if(cosNO <= 0 \|\| cosNI >= 0)			if(cosNO <= 0 \|\| cosNI >= 0)
	return make_float3 (0, 0, 0);			return make_float3(0, 0, 0);
	// compute half-vector of the refraction (eq. 16)
				/* compute half-vector of the refraction (eq. 16) */
	float3 ht = -(m_eta * omega_in + I);			float3 ht = -(m_eta * omega_in + I);
	float3 Ht = normalize(ht);			float3 Ht = normalize(ht);
	float cosHO = dot(Ht, I);			float cosHO = dot(Ht, I);

	float cosHI = dot(Ht, omega_in);			float cosHI = dot(Ht, omega_in);
	// eq. 33: first we calculate D(m) with m=Ht:
				/* eq. 33: first we calculate D(m) with m=Ht: */
	float alpha2 = m_ab * m_ab;			float alpha2 = m_ab * m_ab;
	float cosThetaM = min(dot(N, Ht), 1.0f);			float cosThetaM = min(dot(N, Ht), 1.0f);
	float cosThetaM2 = cosThetaM * cosThetaM;			float cosThetaM2 = cosThetaM * cosThetaM;
	float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;			float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
	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. 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;
	// probability
	float invHt2 = 1 / dot(ht, ht);			/* probability */
	pdf = D fabsf(cosThetaM) * (fabsf(cosHI) * (m_eta * m_eta)) * invHt2;			float Ht2 = dot(ht, ht);
	float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D) * invHt2) / cosNO;
				/* eq. 2 in distribution of visible normals sampling
				* pm = Dw = G1o * dot(m, I) * D / dot(N, I); */

				/* out = fabsf(cosHI * cosHO) * (m_eta * m_eta) * G * D / (cosNO * Ht2)
				* pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2 */
				float common = D * (m_eta * m_eta) / (cosNO * Ht2);
				float out = G * fabsf(cosHI * cosHO) * common;
				pdf = G1o cosHO * fabsf(cosHI) * common;

	return make_float3 (out, out, out);			return make_float3(out, out, out);
	}			}

	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)			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
	// eq. 35,36:
	// we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)
	//tttt and sin(atan(x)) == x/sqrt(1+x^2)
	float alpha2 = m_ab * m_ab;
	float tanThetaM, cosThetaM;

	if(alpha2 == 0.0f) {			/* importance sampling with distribution of visible normals. vectors are
	tanThetaM = 0.0f;			* transformed to local space before and after */
	cosThetaM = 1.0f;			float3 local_I, local_m;
	}
	else {
	tanThetaM = safe_sqrtf(-alpha2 * logf(1 - randu));
	cosThetaM = 1 / safe_sqrtf(1 + tanThetaM * tanThetaM);
	}

	float sinThetaM = cosThetaM * tanThetaM;			local_I.x = dot(X, I);
	float phiM = M_2PI_F * randv;			local_I.y = dot(Y, I);
	float3 m = (cosf(phiM) * sinThetaM) * X +			local_I.z = dot(Z, I);
	(sinf(phiM) * sinThetaM) * Y +
	( cosThetaM) * Z;

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

				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);

				/* reflection or refraction? */
	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 alpha2 = m_ab * m_ab;
	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 BRDFcosNI /
	// 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
	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;
	// eq. 20: (FGD)/(4inon)
	float out = (G * D) * 0.25f / cosNO;			/* see eval function for derivation */
				float common = D * 0.25f / cosNO;
				float out = G * common;
				pdf = G1o common;

	*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
	}			}
	}			}
	}			}
	else {			else {
	// CAUTION: the i and o variables are inverted relative to the paper			/* CAUTION: the i and o variables are inverted relative to the paper
	// eq. 39 - compute actual refractive direction			* eq. 39 - compute actual refractive direction */
	float3 R, T;			float3 R, T;
	#ifdef __RAY_DIFFERENTIALS__			#ifdef __RAY_DIFFERENTIALS__
	float3 dRdx, dRdy, dTdx, dTdy;			float3 dRdx, dRdy, dTdx, dTdy;
	#endif			#endif
	float m_eta = sc->data1;			float m_eta = sc->data1;
	bool inside;			bool inside;

	fresnel_dielectric(m_eta, m, I, &R, &T,			fresnel_dielectric(m_eta, m, I, &R, &T,
	#ifdef __RAY_DIFFERENTIALS__			#ifdef __RAY_DIFFERENTIALS__
	dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,			dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,
	#endif			#endif
	&inside);			&inside);

	if(!inside) {			if(!inside) {
	*omega_in = T;			*omega_in = T;

	#ifdef __RAY_DIFFERENTIALS__			#ifdef __RAY_DIFFERENTIALS__
	*domega_in_dx = dTdx;			*domega_in_dx = dTdx;
	*domega_in_dy = dTdy;			*domega_in_dy = dTdy;
	#endif			#endif

	if (m_ab <= 1e-4f \|\| fabsf(m_eta - 1.0f) < 1e-4f) {			if(m_ab <= 1e-4f \|\| fabsf(m_eta - 1.0f) < 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 {
	// eq. 33			/* eq. 33 */
				float alpha2 = m_ab * m_ab;
	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 BRDFcosNI /
	// 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;
	// 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);
	// eq. 38 and eq. 17			/* see eval function for derivation */
	pdf = pm (m_eta * m_eta) * fabsf(cosHI) / Ht2;			float common = D * (m_eta * m_eta) / (cosNO * Ht2);
				float out = G * fabsf(cosHI * cosHO) * common;
				pdf = G1o cosHO * fabsf(cosHI) * common;

	*eval = make_float3(out, out, out);			*eval = make_float3(out, out, out);
	}			}
	}			}
	}			}
	}			}
	return (m_refractive) ? LABEL_TRANSMIT\|LABEL_GLOSSY : LABEL_REFLECT\|LABEL_GLOSSY;			return (m_refractive) ? LABEL_TRANSMIT\|LABEL_GLOSSY : LABEL_REFLECT\|LABEL_GLOSSY;
	}			}

	CCL_NAMESPACE_END			CCL_NAMESPACE_END

	#endif /* __BSDF_MICROFACET_H__ */			#endif /* __BSDF_MICROFACET_H__ */