Page MenuHome
Differential D572 Diff 1840 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 */ /* Approximate erf and erfinv implementations
*
* Adapted from code (C) Copyright John Maddock 2006.
* Use, modification and distribution are subject to the
* Boost Software License, Version 1.0. (See accompanying file
* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) */
ccl_device float approx_erff_impl(float z)
{
float result;
if(z < 0.5f) {
if(z < 1e-10f) {
if(z == 0) {
result = 0;
}
else {
float c = 0.0033791670f;
result = z * 1.125f + z * c;
}
}
else {
float Y = 1.044948577f;
float zz = z * z;
float num = (((-0.007727583f * zz) + -0.050999073f)*zz + -0.338165134f)*zz + 0.083430589f;
float denom = (((0.000370900f * zz) + 0.008585719f)*zz + 0.087522260f)*zz + 0.455004033f;
result = z * (Y + num / denom);
}
}
else if(z < 2.5f) {
if(z < 1.5f) {
float Y = 0.4059357643f;
float fz = z - 0.5f;
float num = (((0.088890036f * fz) + 0.191003695f)*fz + 0.178114665f)*fz + -0.098090592f;
float denom = (((0.123850974f * fz) + 0.578052804f)*fz + 1.426280048f)*fz + 1.847590709f;
result = Y + num / denom;
result *= expf(-z * z) / z;
}
else {
float Y = 0.506728172f;
float fz = z - 1.5f;
float num = (((0.017567943f * fz) + 0.043948189f)*fz + 0.038654037f)*fz + -0.024350047f;
float denom = (((0.325732924f * fz) + 0.982403709f)*fz + 1.539914949f)*fz + 1;
result = Y + num / denom;
result *= expf(-z * z) / z;
}
result = 1 - result;
}
else {
result = 1;
}
return result;
}
ccl_device float approx_erff(float z)
{
float s = 1.0f;
if(z < 0.0f) {
s = -1.0f;
z = -z;
}
return s * approx_erff_impl(z);
}
ccl_device float approx_erfinvf_impl(float p, float q)
{
float result = 0;
if(p <= 0.5f) {
float Y = 0.089131474f;
float g = p * (p + 10);
float num = (((-0.012692614f * p) + 0.033480662f)*p + -0.008368748f)*p + -0.000508781f;
float denom = (((1.562215583f * p) + -1.565745582f)*p + -0.970005043f)*p + 1.0f;
float r = num / denom;
result = g * Y + g * r;
}
else if(q >= 0.25f) {
float Y = 2.249481201f;
float g = sqrtf(-2 * logf(q));
float xs = q - 0.25f;
float num = (((17.644729840f * xs) + 8.370503283f)*xs + 0.105264680f)*xs + -0.202433508f;
float denom = (((-28.660818049f * xs) + 3.971343795f)*xs + 6.242641248f)*xs + 1.0f;
float r = num / denom;
result = g / (Y + r);
}
else {
float x = sqrtf(-logf(q));
if(x < 3) {
float Y = 0.807220458f;
float xs = x - 1.125f;
float num = (((0.387079738f * xs) + 0.117030156f)*xs + -0.163794047f)*xs + -0.131102781f;
float denom = (((4.778465929f * xs) + 5.381683457f)*xs + 3.466254072f)*xs + 1.0f;
float R = num / denom;
result = Y * x + R * x;
}
else {
float Y = 0.939955711f;
float xs = x - 3;
float num = (((0.009508047f * xs) + 0.018557330f)*xs + -0.002224265f)*xs + -0.035035378f;
float denom = (((0.220091105f * xs) + 0.762059164f)*xs + 1.365334981f)*xs + 1.0f;
float R = num / denom;
result = Y * x + R * x;
}
}
return result;
}
ccl_device float approx_erfinvf(float z)
{
float p, q, s;
if(z < 0) {
p = -z;
q = 1 - p;
s = -1;
}
else {
p = z;
q = 1 - z;
s = 1;
}
return s * approx_erfinvf_impl(p, q);
}
/* Beckmann and GGX microfacet importance sampling from:
*
* Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals.
* E. Heitz and E. d'Eon, EGSR 2014 */
ccl_device_inline void microfacet_beckmann_sample_slopes(
const float cos_theta_i, const float sin_theta_i,
const float alpha_x, const float alpha_y,
float randu, float randv, float *slope_x, float *slope_y,
float *G1i)
{
/* special case (normal incidence) */
if(cos_theta_i >= 0.99999f) {
const float r = sqrtf(-logf(randu));
const float phi = M_2PI_F * randv;
*slope_x = r * cosf(phi);
*slope_y = r * sinf(phi);
*G1i = 1.0f;
return;
}
/* precomputations */
const float tan_theta_i = sin_theta_i/cos_theta_i;
const float inv_a = tan_theta_i;
const float a = 1.0f/inv_a;
const float erf_a = approx_erff(a);
const float exp_a2 = expf(-a*a);
const float SQRT_PI_INV = 0.56418958354f;
const float Lambda = 0.5f*(erf_a - 1.0f) + (0.5f*SQRT_PI_INV)*(exp_a2*inv_a);
const float G1 = 1.0f/(1.0f + Lambda); /* masking */
const float C = 1.0f - G1 * erf_a;
*G1i = G1;
/* sample slope X */
if(randu < C) {
/* rescale randu */
randu = randu / C;
const float w_1 = 0.5f * SQRT_PI_INV * sin_theta_i * exp_a2;
const float w_2 = cos_theta_i * (0.5f - 0.5f*erf_a);
const float p = w_1 / (w_1 + w_2);
if(randu < p) {
randu = randu / p;
*slope_x = -sqrtf(-logf(randu*exp_a2));
}
else {
randu = (randu - p) / (1.0f - p);
*slope_x = approx_erfinvf(randu - 1.0f - randu*erf_a);
}
}
else {
/* rescale randu */
randu = (randu - C) / (1.0f - C);
*slope_x = approx_erfinvf((-1.0f + 2.0f*randu)*erf_a);
const float p = (-(*slope_x)*sin_theta_i + cos_theta_i) / (2.0f*cos_theta_i);
if(randv > p) {
*slope_x = -(*slope_x);
randv = (randv - p) / (1.0f - p);
}
else
randv = randv / p;
}
/* sample slope Y */
*slope_y = approx_erfinvf(2.0f*randv - 1.0f);
}
ccl_device_inline void microfacet_ggx_sample_slopes(
const float cos_theta_i, const float sin_theta_i,
const float alpha_x, const float alpha_y,
float randu, float randv, float *slope_x, float *slope_y,
float *G1i)
{
/* special case (normal incidence) */
if(cos_theta_i >= 0.99999f) {
const float r = sqrtf(randu/(1.0f - randu));
const float phi = M_2PI_F * randv;
*slope_x = r * cosf(phi);
*slope_y = r * sinf(phi);
*G1i = 1.0f;
return;
}
/* precomputations */
const float tan_theta_i = sin_theta_i/cos_theta_i;
const float G1_inv = 0.5f * (1.0f + safe_sqrtf(1.0f + tan_theta_i*tan_theta_i));
*G1i = 1.0f/G1_inv;
/* sample slope_x */
const float A = 2.0f*randu*G1_inv - 1.0f;
const float AA = A*A;
const float tmp = 1.0f/(AA - 1.0f);
const float B = tan_theta_i;
const float BB = B*B;
const float D = safe_sqrtf(BB*(tmp*tmp) - (AA - BB)*tmp);
const float slope_x_1 = B*tmp - D;
const float slope_x_2 = B*tmp + D;
*slope_x = (A < 0.0f || slope_x_2*tan_theta_i > 1.0f)? slope_x_1: slope_x_2;
/* sample slope_y */
float S;
if(randv > 0.5f) {
S = 1.0f;
randv = 2.0f*(randv - 0.5f);
}
else {
S = -1.0f;
randv = 2.0f*(0.5f - randv);
}
const float z = (randv*(randv*(randv*0.27385f - 0.73369f) + 0.46341f)) / (randv*(randv*(randv*0.093073f + 0.309420f) - 1.000000f) + 0.597999f);
*slope_y = S * z * safe_sqrtf(1.0f + (*slope_x)*(*slope_x));
}
ccl_device_inline float3 microfacet_sample_stretched(const float3 omega_i,
const float alpha_x, const float alpha_y,
const float randu, const float randv,
bool beckmann, float *G1i)
{
/* 1. stretch omega_i */
float3 omega_i_ = make_float3(alpha_x * omega_i.x, alpha_y * omega_i.y, omega_i.z);
omega_i_ = normalize(omega_i_);
/* get polar coordinates of omega_i_ */
float costheta_ = 1.0f;
float sintheta_ = 0.0f;
float cosphi_ = 1.0f;
float sinphi_ = 0.0f;
if(omega_i_.z < 0.99999f) {
costheta_ = omega_i_.z;
sintheta_ = safe_sqrtf(1.0f - costheta_*costheta_);
float invlen = 1.0f/sintheta_;
cosphi_ = omega_i_.x * invlen;
sinphi_ = omega_i_.y * invlen;
}
/* 2. sample P22_{omega_i}(x_slope, y_slope, 1, 1) */
float slope_x, slope_y;
if(beckmann)
microfacet_beckmann_sample_slopes(costheta_, sintheta_,
alpha_x, alpha_y, randu, randv, &slope_x, &slope_y, G1i);
else
microfacet_ggx_sample_slopes(costheta_, sintheta_,
alpha_x, alpha_y, randu, randv, &slope_x, &slope_y, G1i);
/* 3. rotate */
float tmp = cosphi_*slope_x - sinphi_*slope_y;
slope_y = sinphi_*slope_x + cosphi_*slope_y;
slope_x = tmp;
/* 4. unstretch */
slope_x = alpha_x * slope_x;
slope_y = alpha_y * slope_y;
/* 5. compute normal */
return normalize(make_float3(-slope_x, -slope_y, 1.0f));
}
/* GGX microfacet with Smith shadow-masking from:
*
* Microfacet Models for Refraction through Rough Surfaces
* B. Walter, S. R. Marschner, H. Li, K. E. Torrance, EGSR 2007
*
* Anisotropic from:
*
* Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs
* E. Heitz, Research Report 2014
*
* Anisotropy is only supported for reflection currently, but adding it for
* tranmission is just a matter of copying code from reflection if needed. */
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); /* alpha_x */
sc->data1 = sc->data0; /* alpha_y */
sc->type = CLOSURE_BSDF_MICROFACET_GGX_ID;
return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
}
ccl_device int bsdf_microfacet_ggx_aniso_setup(ShaderClosure *sc)
{
sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* alpha_x */
sc->data1 = clamp(sc->data1, 0.0f, 1.0f); /* alpha_y */
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;
} }
ccl_device int bsdf_microfacet_ggx_refraction_setup(ShaderClosure *sc) ccl_device int bsdf_microfacet_ggx_refraction_setup(ShaderClosure *sc)
{ {
sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ag */ sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* alpha_x */
sc->data1 = sc->data1; /* alpha_y */
sc->type = CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID; sc->type = CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY; return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
} }
ccl_device void bsdf_microfacet_ggx_blur(ShaderClosure *sc, float roughness) ccl_device void bsdf_microfacet_ggx_blur(ShaderClosure *sc, float roughness)
{ {
sc->data0 = fmaxf(roughness, sc->data0); /* m_ag */ sc->data0 = fmaxf(roughness, sc->data0); /* alpha_x */
sc->data1 = fmaxf(roughness, sc->data1); /* alpha_y */
} }
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 alpha_x = sc->data0;
float alpha_y = 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 || fmaxf(alpha_x, alpha_y) <= 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) float alpha2 = alpha_x * alpha_y;
// eq. 33: first we calculate D(m) with m=Hr: float D, G1o, G1i;
float alpha2 = m_ag * m_ag;
float cosThetaM = dot(N, Hr); if(alpha_x == alpha_y) {
/* isotropic
* eq. 20: (F*G*D)/(4*in*on)
* eq. 33: first we calculate D(m) */
float cosThetaM = dot(N, m);
float cosThetaM2 = cosThetaM * cosThetaM; float cosThetaM2 = cosThetaM * cosThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2;
float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2; float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
/* eq. 34: now calculate G1(i,m) and G1(o,m) */
G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI)));
}
else {
/* anisotropic */
float3 X, Y, Z = N;
make_orthonormals_tangent(Z, sc->T, &X, &Y);
/* distribution */
float3 local_m = make_float3(dot(X, m), dot(Y, m), dot(Z, m));
float slope_x = -local_m.x/(local_m.z*alpha_x);
float slope_y = -local_m.y/(local_m.z*alpha_y);
float slope_len = 1 + slope_x*slope_x + slope_y*slope_y;
float cosThetaM = local_m.z;
float cosThetaM2 = cosThetaM * cosThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2; float cosThetaM4 = cosThetaM2 * cosThetaM2;
float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
// eq. 34: now calculate G1(i,m) and G1(o,m) D = 1 / ((slope_len * slope_len) * M_PI_F * alpha2 * cosThetaM4);
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))); /* G1(i,m) and G1(o,m) */
float tanThetaO2 = (1 - cosNO * cosNO) / (cosNO * cosNO);
float cosPhiO = dot(I, X);
float sinPhiO = dot(I, Y);
float alphaO2 = (cosPhiO*cosPhiO)*(alpha_x*alpha_x) + (sinPhiO*sinPhiO)*(alpha_y*alpha_y);
alphaO2 /= cosPhiO*cosPhiO + sinPhiO*sinPhiO;
G1o = 2 / (1 + safe_sqrtf(1 + alphaO2 * tanThetaO2));
float tanThetaI2 = (1 - cosNI * cosNI) / (cosNI * cosNI);
float cosPhiI = dot(omega_in, X);
float sinPhiI = dot(omega_in, Y);
float alphaI2 = (cosPhiI*cosPhiI)*(alpha_x*alpha_x) + (sinPhiI*sinPhiI)*(alpha_y*alpha_y);
alphaI2 /= cosPhiI*cosPhiI + sinPhiI*sinPhiI;
G1i = 2 / (1 + safe_sqrtf(1 + alphaI2 * tanThetaI2));
}
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 alpha_x = sc->data0;
float m_eta = sc->data1; float alpha_y = sc->data1;
float m_eta = sc->data2;
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 || fmaxf(alpha_x, alpha_y) <= 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:
float alpha2 = m_ag * m_ag; float D, G1o, G1i;
/* eq. 33: first we calculate D(m) with m=Ht: */
float alpha2 = alpha_x * alpha_y;
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)); D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
// eq. 34: now calculate G1(i,m) and G1(o,m)
float G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO))); /* eq. 34: now calculate G1(i,m) and G1(o,m) */
float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
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 alpha_x = sc->data0;
float alpha_y = 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;
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;
if(alpha_x == alpha_y)
make_orthonormals(Z, &X, &Y); make_orthonormals(Z, &X, &Y);
// generate a random microfacet normal m else
// eq. 35,36: make_orthonormals_tangent(Z, sc->T, &X, &Y);
// we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)
//tttt and sin(atan(x)) == x/sqrt(1+x^2) /* importance sampling with distribution of visible normals. vectors are
float alpha2 = m_ag * m_ag; * transformed to local space before and after */
float tanThetaM2 = alpha2 * randu / (1 - randu); float3 local_I = make_float3(dot(X, I), dot(Y, I), cosNO);
float cosThetaM = 1 / safe_sqrtf(1 + tanThetaM2); float3 local_m;
float sinThetaM = cosThetaM * safe_sqrtf(tanThetaM2); float G1o;
float phiM = M_2PI_F * randv;
float3 m = (cosf(phiM) * sinThetaM) * X + local_m = microfacet_sample_stretched(local_I, alpha_x, alpha_y,
(sinf(phiM) * sinThetaM) * Y + randu, randv, false, &G1o);
( cosThetaM) * Z;
float3 m = X*local_m.x + Y*local_m.y + Z*local_m.z;
float cosThetaM = local_m.z;
/* 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(fmaxf(alpha_x, alpha_y) <= 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 = alpha_x * alpha_y;
float D, G1i;
if(alpha_x == alpha_y) {
/* isotropic */
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 tanThetaM2 = 1/(cosThetaM2) - 1;
// eq. 24 D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
float pm = D * cosThetaM;
// convert into pdf of the sampled direction /* eval BRDF*cosNI */
// 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)
float G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO))); /* eq. 34: now calculate G1(i,m) */
float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI)));
float G = G1o * G1i; }
// eq. 20: (F*G*D)/(4*in*on) else {
float out = (G * D) * 0.25f / cosNO; /* anisotropic distribution */
float3 local_m = make_float3(dot(X, m), dot(Y, m), dot(Z, m));
float slope_x = -local_m.x/(local_m.z*alpha_x);
float slope_y = -local_m.y/(local_m.z*alpha_y);
float slope_len = 1 + slope_x*slope_x + slope_y*slope_y;
float cosThetaM = local_m.z;
float cosThetaM2 = cosThetaM * cosThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2;
D = 1 / ((slope_len * slope_len) * M_PI_F * alpha2 * cosThetaM4);
/* calculate G1(i,m) */
float cosNI = dot(N, *omega_in);
float tanThetaI2 = (1 - cosNI * cosNI) / (cosNI * cosNI);
float cosPhiI = dot(*omega_in, X);
float sinPhiI = dot(*omega_in, Y);
float alphaI2 = (cosPhiI*cosPhiI)*(alpha_x*alpha_x) + (sinPhiI*sinPhiI)*(alpha_y*alpha_y);
alphaI2 /= cosPhiI*cosPhiI + sinPhiI*sinPhiI;
G1i = 2 / (1 + safe_sqrtf(1 + alphaI2 * tanThetaI2));
}
/* see eval function for derivation */
float common = (G1o * D) * 0.25f / cosNO;
float out = G1i * common;
*pdf = 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->data2;
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(fmaxf(alpha_x, alpha_y) <= 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 = alpha_x * alpha_y;
float cosThetaM2 = cosThetaM * cosThetaM; float cosThetaM2 = cosThetaM * cosThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2; float cosThetaM4 = cosThetaM2 * cosThetaM2;
float tanThetaM2 = 1/(cosThetaM2) - 1;
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)
float G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO))); /* eq. 34: now calculate G1(i,m) */
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;
// 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 = (G1o * D) * (m_eta * m_eta) / (cosNO * Ht2);
float out = G1i * fabsf(cosHI * cosHO) * common;
*pdf = 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); /* alpha_x */
sc->data1 = sc->data0; /* alpha_y */
sc->type = CLOSURE_BSDF_MICROFACET_BECKMANN_ID;
return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
}
ccl_device int bsdf_microfacet_beckmann_aniso_setup(ShaderClosure *sc)
{
sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* alpha_x */
sc->data1 = clamp(sc->data1, 0.0f, 1.0f); /* alpha_y */
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;
} }
ccl_device int bsdf_microfacet_beckmann_refraction_setup(ShaderClosure *sc) ccl_device int bsdf_microfacet_beckmann_refraction_setup(ShaderClosure *sc)
{ {
sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ab */ sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* alpha_x */
sc->data1 = sc->data1; /* alpha_y */
sc->type = CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID; sc->type = CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY; return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
} }
ccl_device void bsdf_microfacet_beckmann_blur(ShaderClosure *sc, float roughness) ccl_device void bsdf_microfacet_beckmann_blur(ShaderClosure *sc, float roughness)
{ {
sc->data0 = fmaxf(roughness, sc->data0); /* m_ab */ sc->data0 = fmaxf(roughness, sc->data0); /* alpha_x */
sc->data1 = fmaxf(roughness, sc->data1); /* alpha_y */
} }
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 alpha_x = sc->data0;
float alpha_y = 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 || fmaxf(alpha_x, alpha_y) <= 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: float alpha2 = alpha_x * alpha_y;
float alpha2 = m_ab * m_ab; float D, G1o, G1i;
float cosThetaM = dot(N, Hr);
if(alpha_x == alpha_y) {
/* isotropic
* eq. 20: (F*G*D)/(4*in*on)
* eq. 25: first we calculate D(m) */
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); D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 * cosThetaM4);
// eq. 26, 27: now calculate G1(i,m) and G1(o,m)
float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO))); /* eq. 26, 27: now calculate G1(i,m) and G1(o,m) */
float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI))); float ao = 1 / (alpha_x * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f; float ai = 1 / (alpha_x * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f; G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
}
else {
/* anisotropic */
float3 X, Y, Z = N;
make_orthonormals_tangent(Z, sc->T, &X, &Y);
/* distribution */
float3 local_m = make_float3(dot(X, m), dot(Y, m), dot(Z, m));
float slope_x = -local_m.x/(local_m.z*alpha_x);
float slope_y = -local_m.y/(local_m.z*alpha_y);
float cosThetaM = local_m.z;
float cosThetaM2 = cosThetaM * cosThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2;
D = expf(-slope_x*slope_x - slope_y*slope_y) / (M_PI_F * alpha2 * cosThetaM4);
/* G1(i,m) and G1(o,m) */
float tanThetaO2 = (1 - cosNO * cosNO) / (cosNO * cosNO);
float cosPhiO = dot(I, X);
float sinPhiO = dot(I, Y);
float alphaO2 = (cosPhiO*cosPhiO)*(alpha_x*alpha_x) + (sinPhiO*sinPhiO)*(alpha_y*alpha_y);
alphaO2 /= cosPhiO*cosPhiO + sinPhiO*sinPhiO;
float tanThetaI2 = (1 - cosNI * cosNI) / (cosNI * cosNI);
float cosPhiI = dot(omega_in, X);
float sinPhiI = dot(omega_in, Y);
float alphaI2 = (cosPhiI*cosPhiI)*(alpha_x*alpha_x) + (sinPhiI*sinPhiI)*(alpha_y*alpha_y);
alphaI2 /= cosPhiI*cosPhiI + sinPhiI*sinPhiI;
float ao = 1 / (safe_sqrtf(alphaO2 * tanThetaO2));
float ai = 1 / (safe_sqrtf(alphaI2 * tanThetaI2));
G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
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 alpha_x = sc->data0;
float m_eta = sc->data1; float alpha_y = sc->data1;
float m_eta = sc->data2;
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 || fmaxf(alpha_x, alpha_y) <= 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:
float alpha2 = m_ab * m_ab; /* eq. 33: first we calculate D(m) with m=Ht: */
float alpha2 = alpha_x * alpha_y;
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)
float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO))); /* eq. 26, 27: now calculate G1(i,m) and G1(o,m) */
float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI))); float ao = 1 / (alpha_x * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
float ai = 1 / (alpha_x * 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 alpha_x = sc->data0;
float alpha_y = 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;
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;
if(alpha_x == alpha_y)
make_orthonormals(Z, &X, &Y); make_orthonormals(Z, &X, &Y);
// generate a random microfacet normal m else
// eq. 35,36: make_orthonormals_tangent(Z, sc->T, &X, &Y);
// we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)
//tttt and sin(atan(x)) == x/sqrt(1+x^2) /* importance sampling with distribution of visible normals. vectors are
float alpha2 = m_ab * m_ab; * transformed to local space before and after */
float tanThetaM, cosThetaM; float3 local_I = make_float3(dot(X, I), dot(Y, I), cosNO);
float3 local_m;
if(alpha2 == 0.0f) { float G1o;
tanThetaM = 0.0f;
cosThetaM = 1.0f; local_m = microfacet_sample_stretched(local_I, alpha_x, alpha_x,
} randu, randv, true, &G1o);
else {
tanThetaM = safe_sqrtf(-alpha2 * logf(1 - randu));
cosThetaM = 1 / safe_sqrtf(1 + tanThetaM * tanThetaM);
}
float sinThetaM = cosThetaM * tanThetaM;
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 = local_m.z;
/* 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(fmaxf(alpha_x, alpha_y) <= 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 = alpha_x * alpha_y;
float D, G1i;
if(alpha_x == alpha_y) {
/* istropic distribution */
float cosThetaM2 = cosThetaM * cosThetaM; float cosThetaM2 = cosThetaM * cosThetaM;
float tanThetaM2 = tanThetaM * tanThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2; float cosThetaM4 = cosThetaM2 * cosThetaM2;
float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 * cosThetaM4); float tanThetaM2 = 1/(cosThetaM2) - 1;
// eq. 24 D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 * cosThetaM4);
float pm = D * cosThetaM;
// convert into pdf of the sampled direction /* eval BRDF*cosNI */
// 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)
float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO))); /* eq. 26, 27: now calculate G1(i,m) */
float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI))); float ai = 1 / (alpha_x * 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; 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; }
else {
/* anisotropic distribution */
float3 local_m = make_float3(dot(X, m), dot(Y, m), dot(Z, m));
float slope_x = -local_m.x/(local_m.z*alpha_x);
float slope_y = -local_m.y/(local_m.z*alpha_y);
float cosThetaM = local_m.z;
float cosThetaM2 = cosThetaM * cosThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2;
D = expf(-slope_x*slope_x - slope_y*slope_y) / (M_PI_F * alpha2 * cosThetaM4);
/* G1(i,m) */
float cosNI = dot(N, *omega_in);
float tanThetaI2 = (1 - cosNI * cosNI) / (cosNI * cosNI);
float cosPhiI = dot(*omega_in, X);
float sinPhiI = dot(*omega_in, Y);
float alphaI2 = (cosPhiI*cosPhiI)*(alpha_x*alpha_x) + (sinPhiI*sinPhiI)*(alpha_y*alpha_y);
alphaI2 /= cosPhiI*cosPhiI + sinPhiI*sinPhiI;
float ai = 1 / (safe_sqrtf(alphaI2 * tanThetaI2));
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->data2;
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) {
// some high number for MIS if(fmaxf(alpha_x, alpha_y) <= 1e-4f || fabsf(m_eta - 1.0f) < 1e-4f) {
/* 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 = alpha_x * alpha_y;
float cosThetaM2 = cosThetaM * cosThetaM; float cosThetaM2 = cosThetaM * cosThetaM;
float tanThetaM2 = tanThetaM * tanThetaM;
float cosThetaM4 = cosThetaM2 * cosThetaM2; float cosThetaM4 = cosThetaM2 * cosThetaM2;
float tanThetaM2 = 1/(cosThetaM2) - 1;
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)
float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO))); /* eq. 26, 27: now calculate G1(i,m) */
float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI))); float ai = 1 / (alpha_x * 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 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__ */