Page MenuHome
Differential D572 Diff 1819 intern/cycles/kernel/closure/bsdf_microfacet.h

Changeset View

Standalone View

View Options

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 void microfacet_beckmann_sample_slopes(const float theta_i,
float randu, float randv, float *slope_x, float *slope_y)
{
/* special case (normal incidence) */
if(theta_i < 0.0001f) {
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 sin_theta_i = sinf(theta_i);
const float cos_theta_i = cosf(theta_i);
const float tan_theta_i = sin_theta_i/cos_theta_i;
const float a = 1.0f / 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.5f*SQRT_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(randu*exp_a2));
}
else {
randu = (randu-p) / (1.0f - p);
*slope_x = erf_inv(randu - 1.0f - randu*erf_a);
}
}
else {
/* rescale randu */
randu = (randu - C) / (1.0f - C);
*slope_x = erf_inv((-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 = erf_inv(2.0f*randv - 1.0f);
}
ccl_device void microfacet_ggx_sample_slopes(const float theta_i,
float randu, float randv, float *slope_x, float *slope_y)
{
/* special case (normal incidence) */
if(theta_i < 0.0001f) {
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 = tanf(theta_i);
const float a = 1.0f / (tan_theta_i);
const float G1 = 2.0f / (1.0f + sqrtf(1.0f + 1.0f/(a*a)));
/* sample slope_x */
const float A = 2.0f*randu/G1 - 1.0f;
const float tmp = 1.0f / (A*A - 1.0f);
const float B = tan_theta_i;
const float D = sqrtf(B*B*tmp*tmp - (A*A-B*B)*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 > 1.0f/tan_theta_i) ? 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 * 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 theta_ = 0.0f;
float phi_ = 0.0f;
if(omega_i_.z < 0.99999f) {
theta_ = safe_acosf(omega_i_.z);
phi_ = atan2f(omega_i_.y, omega_i_.x);
}
/* 2. sample P22_{omega_i}(x_slope, y_slope, 1, 1) */
float slope_x, slope_y;
if(beckmann)
microfacet_beckmann_sample_slopes(theta_, randu, randv, &slope_x, &slope_y);
else
microfacet_ggx_sample_slopes(theta_, randu, randv, &slope_x, &slope_y);
/* 3. rotate */
float cosphi_ = cosf(phi_), sinphi_ = sinf(phi_);
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 = sqrtf(slope_x*slope_x + slope_y*slope_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: (F*G*D)/(4*in*on)
// eq. 33: first we calculate D(m) with m=Hr: /* eq. 20: (F*G*D)/(4*in*on)
* 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 = X*local_m.x + Y*local_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 BRDF*cosNI */
// 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: (F*G*D)/(4*in*on)
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 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;
// 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: (F*G*D)/(4*in*on)
// eq. 25: first we calculate D(m) with m=Hr: /* eq. 20: (F*G*D)/(4*in*on)
* 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 { local_I.x = dot(X, I);
tanThetaM = safe_sqrtf(-alpha2 * logf(1 - randu)); local_I.y = dot(Y, I);
cosThetaM = 1 / safe_sqrtf(1 + tanThetaM * tanThetaM); local_I.z = dot(Z, I);
}
float sinThetaM = cosThetaM * tanThetaM; microfacet_sample_stretched(local_I, m_ab, m_ab, randu, randv, true, &local_m);
float phiM = M_2PI_F * randv;
float3 m = (cosf(phiM) * sinThetaM) * X +
(sinf(phiM) * sinThetaM) * Y +
( cosThetaM) * Z;
float3 m = X*local_m.x + Y*local_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 BRDF*cosNI */
// 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: (F*G*D)/(4*in*on)
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 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;
// 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__ */