intern/cycles/kernel/closure/bsdf_microfacet.h

   1 /*
   2  * Adapted from Open Shading Language with this license:
   3  *
   4  * Copyright (c) 2009-2010 Sony Pictures Imageworks Inc., et al.
   5  * All Rights Reserved.
   6  *
   7  * Modifications Copyright 2011, Blender Foundation.
   8  *
   9  * Redistribution and use in source and binary forms, with or without
  10  * modification, are permitted provided that the following conditions are
  11  * met:
  12  * * Redistributions of source code must retain the above copyright
  13  *   notice, this list of conditions and the following disclaimer.
  14  * * Redistributions in binary form must reproduce the above copyright
  15  *   notice, this list of conditions and the following disclaimer in the
  16  *   documentation and/or other materials provided with the distribution.
  17  * * Neither the name of Sony Pictures Imageworks nor the names of its
  18  *   contributors may be used to endorse or promote products derived from
  19  *   this software without specific prior written permission.
  20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  21  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  22  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  23  * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  24  * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  25  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  26  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  27  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  28  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  29  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  30  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  31  */
  32
  33 #ifndef __BSDF_MICROFACET_H__
  34 #define __BSDF_MICROFACET_H__
  35
  36 CCL_NAMESPACE_BEGIN
  37
  38 typedef ccl_addr_space struct MicrofacetExtra {
  39         float3 color, cspec0;
  40         float clearcoat;
  41 } MicrofacetExtra;
  42
  43 typedef ccl_addr_space struct MicrofacetBsdf {
  44         SHADER_CLOSURE_BASE;
  45
  46         float alpha_x, alpha_y, ior;
  47         MicrofacetExtra *extra;
  48         float3 T;
  49 } MicrofacetBsdf;
  50
  51 /* Beckmann and GGX microfacet importance sampling. */
  52
  53 ccl_device_inline void microfacet_beckmann_sample_slopes(
  54         KernelGlobals *kg,
  55         const float cos_theta_i, const float sin_theta_i,
  56         float randu, float randv, float *slope_x, float *slope_y,
  57         float *G1i)
  58 {
  59         /* special case (normal incidence) */
  60         if(cos_theta_i >= 0.99999f) {
  61                 const float r = sqrtf(-logf(randu));
  62                 const float phi = M_2PI_F * randv;
  63                 *slope_x = r * cosf(phi);
  64                 *slope_y = r * sinf(phi);
  65                 *G1i = 1.0f;
  66                 return;
  67         }
  68
  69         /* precomputations */
  70         const float tan_theta_i = sin_theta_i/cos_theta_i;
  71         const float inv_a = tan_theta_i;
  72         const float cot_theta_i = 1.0f/tan_theta_i;
  73         const float erf_a = fast_erff(cot_theta_i);
  74         const float exp_a2 = expf(-cot_theta_i*cot_theta_i);
  75         const float SQRT_PI_INV = 0.56418958354f;
  76         const float Lambda = 0.5f*(erf_a - 1.0f) + (0.5f*SQRT_PI_INV)*(exp_a2*inv_a);
  77         const float G1 = 1.0f/(1.0f + Lambda); /* masking */
  78
  79         *G1i = G1;
  80
  81 #if defined(__KERNEL_GPU__)
  82         /* Based on paper from Wenzel Jakob
  83          * An Improved Visible Normal Sampling Routine for the Beckmann Distribution
  84          *
  85          * http://www.mitsuba-renderer.org/~wenzel/files/visnormal.pdf
  86          *
  87          * Reformulation from OpenShadingLanguage which avoids using inverse
  88          * trigonometric functions.
  89          */
  90
  91         /* Sample slope X.
  92          *
  93          * Compute a coarse approximation using the approximation:
  94          *   exp(-ierf(x)^2) ~= 1 - x * x
  95          *   solve y = 1 + b + K * (1 - b * b)
  96          */
  97         float K = tan_theta_i * SQRT_PI_INV;
  98         float y_approx = randu * (1.0f + erf_a + K * (1 - erf_a * erf_a));
  99         float y_exact  = randu * (1.0f + erf_a + K * exp_a2);
 100         float b = K > 0 ? (0.5f - sqrtf(K * (K - y_approx + 1.0f) + 0.25f)) / K : y_approx - 1.0f;
 101
 102         /* Perform newton step to refine toward the true root. */
 103         float inv_erf = fast_ierff(b);
 104         float value  = 1.0f + b + K * expf(-inv_erf * inv_erf) - y_exact;
 105         /* Check if we are close enough already,
 106          * this also avoids NaNs as we get close to the root.
 107          */
 108         if(fabsf(value) > 1e-6f) {
 109                 b -= value / (1.0f - inv_erf * tan_theta_i); /* newton step 1. */
 110                 inv_erf = fast_ierff(b);
 111                 value  = 1.0f + b + K * expf(-inv_erf * inv_erf) - y_exact;
 112                 b -= value / (1.0f - inv_erf * tan_theta_i); /* newton step 2. */
 113                 /* Compute the slope from the refined value. */
 114                 *slope_x = fast_ierff(b);
 115         }
 116         else {
 117                 /* We are close enough already. */
 118                 *slope_x = inv_erf;
 119         }
 120         *slope_y = fast_ierff(2.0f*randv - 1.0f);
 121 #else
 122         /* Use precomputed table on CPU, it gives better perfomance. */
 123         int beckmann_table_offset = kernel_data.tables.beckmann_offset;
 124
 125         *slope_x = lookup_table_read_2D(kg, randu, cos_theta_i,
 126                 beckmann_table_offset, BECKMANN_TABLE_SIZE, BECKMANN_TABLE_SIZE);
 127         *slope_y = fast_ierff(2.0f*randv - 1.0f);
 128 #endif
 129 }
 130
 131 /* GGX microfacet importance sampling from:
 132  *
 133  * Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals.
 134  * E. Heitz and E. d'Eon, EGSR 2014
 135  */
 136
 137 ccl_device_inline void microfacet_ggx_sample_slopes(
 138         const float cos_theta_i, const float sin_theta_i,
 139         float randu, float randv, float *slope_x, float *slope_y,
 140         float *G1i)
 141 {
 142         /* special case (normal incidence) */
 143         if(cos_theta_i >= 0.99999f) {
 144                 const float r = sqrtf(randu/(1.0f - randu));
 145                 const float phi = M_2PI_F * randv;
 146                 *slope_x = r * cosf(phi);
 147                 *slope_y = r * sinf(phi);
 148                 *G1i = 1.0f;
 149
 150                 return;
 151         }
 152
 153         /* precomputations */
 154         const float tan_theta_i = sin_theta_i/cos_theta_i;
 155         const float G1_inv = 0.5f * (1.0f + safe_sqrtf(1.0f + tan_theta_i*tan_theta_i));
 156
 157         *G1i = 1.0f/G1_inv;
 158
 159         /* sample slope_x */
 160         const float A = 2.0f*randu*G1_inv - 1.0f;
 161         const float AA = A*A;
 162         const float tmp = 1.0f/(AA - 1.0f);
 163         const float B = tan_theta_i;
 164         const float BB = B*B;
 165         const float D = safe_sqrtf(BB*(tmp*tmp) - (AA - BB)*tmp);
 166         const float slope_x_1 = B*tmp - D;
 167         const float slope_x_2 = B*tmp + D;
 168         *slope_x = (A < 0.0f || slope_x_2*tan_theta_i > 1.0f)? slope_x_1: slope_x_2;
 169
 170         /* sample slope_y */
 171         float S;
 172
 173         if(randv > 0.5f) {
 174                 S = 1.0f;
 175                 randv = 2.0f*(randv - 0.5f);
 176         }
 177         else {
 178                 S = -1.0f;
 179                 randv = 2.0f*(0.5f - randv);
 180         }
 181
 182         const float z = (randv*(randv*(randv*0.27385f - 0.73369f) + 0.46341f)) / (randv*(randv*(randv*0.093073f + 0.309420f) - 1.000000f) + 0.597999f);
 183         *slope_y = S * z * safe_sqrtf(1.0f + (*slope_x)*(*slope_x));
 184 }
 185
 186 ccl_device_forceinline float3 microfacet_sample_stretched(
 187         KernelGlobals *kg, const float3 omega_i,
 188         const float alpha_x, const float alpha_y,
 189         const float randu, const float randv,
 190         bool beckmann, float *G1i)
 191 {
 192         /* 1. stretch omega_i */
 193         float3 omega_i_ = make_float3(alpha_x * omega_i.x, alpha_y * omega_i.y, omega_i.z);
 194         omega_i_ = normalize(omega_i_);
 195
 196         /* get polar coordinates of omega_i_ */
 197         float costheta_ = 1.0f;
 198         float sintheta_ = 0.0f;
 199         float cosphi_ = 1.0f;
 200         float sinphi_ = 0.0f;
 201
 202         if(omega_i_.z < 0.99999f) {
 203                 costheta_ = omega_i_.z;
 204                 sintheta_ = safe_sqrtf(1.0f - costheta_*costheta_);
 205
 206                 float invlen = 1.0f/sintheta_;
 207                 cosphi_ = omega_i_.x * invlen;
 208                 sinphi_ = omega_i_.y * invlen;
 209         }
 210
 211         /* 2. sample P22_{omega_i}(x_slope, y_slope, 1, 1) */
 212         float slope_x, slope_y;
 213
 214         if(beckmann) {
 215                 microfacet_beckmann_sample_slopes(kg, costheta_, sintheta_,
 216                         randu, randv, &slope_x, &slope_y, G1i);
 217         }
 218         else {
 219                 microfacet_ggx_sample_slopes(costheta_, sintheta_,
 220                         randu, randv, &slope_x, &slope_y, G1i);
 221         }
 222
 223         /* 3. rotate */
 224         float tmp = cosphi_*slope_x - sinphi_*slope_y;
 225         slope_y = sinphi_*slope_x + cosphi_*slope_y;
 226         slope_x = tmp;
 227
 228         /* 4. unstretch */
 229         slope_x = alpha_x * slope_x;
 230         slope_y = alpha_y * slope_y;
 231
 232         /* 5. compute normal */
 233         return normalize(make_float3(-slope_x, -slope_y, 1.0f));
 234 }
 235
 236 /* Calculate the reflection color
 237  *
 238  * If fresnel is used, the color is an interpolation of the F0 color and white
 239  * with respect to the fresnel
 240  *
 241  * Else it is simply white
 242  */
 243 ccl_device_forceinline float3 reflection_color(const MicrofacetBsdf *bsdf, float3 L, float3 H) {
 244         float3 F = make_float3(1.0f, 1.0f, 1.0f);
 245         bool use_fresnel = (bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_FRESNEL_ID
 246                            || bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_CLEARCOAT_ID
 247                            || bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_ANISO_FRESNEL_ID);
 248
 249         if(use_fresnel) {
 250                 float F0 = fresnel_dielectric_cos(1.0f, bsdf->ior);
 251
 252                 F = interpolate_fresnel_color(L, H, bsdf->ior, F0, bsdf->extra->cspec0);
 253         }
 254
 255         return F;
 256 }
 257
 258 ccl_device_forceinline float D_GTR1(float NdotH, float alpha)
 259 {
 260         if(alpha >= 1.0f) return M_1_PI_F;
 261         float alpha2 = alpha*alpha;
 262         float t = 1.0f + (alpha2 - 1.0f) * NdotH*NdotH;
 263         return (alpha2 - 1.0f) / (M_PI_F * logf(alpha2) * t);
 264 }
 265
 266 /* GGX microfacet with Smith shadow-masking from:
 267  *
 268  * Microfacet Models for Refraction through Rough Surfaces
 269  * B. Walter, S. R. Marschner, H. Li, K. E. Torrance, EGSR 2007
 270  *
 271  * Anisotropic from:
 272  *
 273  * Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs
 274  * E. Heitz, Research Report 2014
 275  *
 276  * Anisotropy is only supported for reflection currently, but adding it for
 277  * transmission is just a matter of copying code from reflection if needed. */
 278
 279 ccl_device int bsdf_microfacet_ggx_setup(MicrofacetBsdf *bsdf)
 280 {
 281         bsdf->extra = NULL;
 282
 283         bsdf->alpha_x = saturate(bsdf->alpha_x);
 284         bsdf->alpha_y = bsdf->alpha_x;
 285
 286         bsdf->type = CLOSURE_BSDF_MICROFACET_GGX_ID;
 287
 288         return SD_BSDF|SD_BSDF_HAS_EVAL;
 289 }
 290
 291 ccl_device int bsdf_microfacet_ggx_fresnel_setup(MicrofacetBsdf *bsdf, const ShaderData *sd)
 292 {
 293         bsdf->extra->cspec0.x = saturate(bsdf->extra->cspec0.x);
 294         bsdf->extra->cspec0.y = saturate(bsdf->extra->cspec0.y);
 295         bsdf->extra->cspec0.z = saturate(bsdf->extra->cspec0.z);
 296
 297         float F0 = fresnel_dielectric_cos(1.0f, bsdf->ior);
 298         float F = average(interpolate_fresnel_color(sd->I, bsdf->N, bsdf->ior, F0, bsdf->extra->cspec0));
 299         bsdf->sample_weight *= F;
 300
 301         bsdf->alpha_x = saturate(bsdf->alpha_x);
 302         bsdf->alpha_y = bsdf->alpha_x;
 303
 304         bsdf->type = CLOSURE_BSDF_MICROFACET_GGX_FRESNEL_ID;
 305
 306         return SD_BSDF|SD_BSDF_HAS_EVAL;
 307 }
 308
 309 ccl_device int bsdf_microfacet_ggx_clearcoat_setup(MicrofacetBsdf *bsdf, const ShaderData *sd)
 310 {
 311         bsdf->extra->cspec0.x = saturate(bsdf->extra->cspec0.x);
 312         bsdf->extra->cspec0.y = saturate(bsdf->extra->cspec0.y);
 313         bsdf->extra->cspec0.z = saturate(bsdf->extra->cspec0.z);
 314
 315         float F0 = fresnel_dielectric_cos(1.0f, bsdf->ior);
 316         float F = average(interpolate_fresnel_color(sd->I, bsdf->N, bsdf->ior, F0, bsdf->extra->cspec0));
 317         bsdf->sample_weight *= 0.25f * bsdf->extra->clearcoat * F;
 318
 319         bsdf->alpha_x = saturate(bsdf->alpha_x);
 320         bsdf->alpha_y = bsdf->alpha_x;
 321
 322         bsdf->type = CLOSURE_BSDF_MICROFACET_GGX_CLEARCOAT_ID;
 323
 324         return SD_BSDF|SD_BSDF_HAS_EVAL;
 325 }
 326
 327 ccl_device bool bsdf_microfacet_merge(const ShaderClosure *a, const ShaderClosure *b)
 328 {
 329         const MicrofacetBsdf *bsdf_a = (const MicrofacetBsdf*)a;
 330         const MicrofacetBsdf *bsdf_b = (const MicrofacetBsdf*)b;
 331
 332         return (isequal_float3(bsdf_a->N, bsdf_b->N)) &&
 333                (bsdf_a->alpha_x == bsdf_b->alpha_x) &&
 334                (bsdf_a->alpha_y == bsdf_b->alpha_y) &&
 335                (isequal_float3(bsdf_a->T, bsdf_b->T)) &&
 336                (bsdf_a->ior == bsdf_b->ior) &&
 337                ((bsdf_a->extra == NULL && bsdf_b->extra == NULL) ||
 338                 ((bsdf_a->extra && bsdf_b->extra) &&
 339                  (isequal_float3(bsdf_a->extra->color, bsdf_b->extra->color)) &&
 340                  (isequal_float3(bsdf_a->extra->cspec0, bsdf_b->extra->cspec0)) &&
 341                  (bsdf_a->extra->clearcoat == bsdf_b->extra->clearcoat)));
 342 }
 343
 344 ccl_device int bsdf_microfacet_ggx_aniso_setup(MicrofacetBsdf *bsdf)
 345 {
 346         bsdf->extra = NULL;
 347
 348         bsdf->alpha_x = saturate(bsdf->alpha_x);
 349         bsdf->alpha_y = saturate(bsdf->alpha_y);
 350
 351         bsdf->type = CLOSURE_BSDF_MICROFACET_GGX_ANISO_ID;
 352
 353         return SD_BSDF|SD_BSDF_HAS_EVAL;
 354 }
 355
 356 ccl_device int bsdf_microfacet_ggx_aniso_fresnel_setup(MicrofacetBsdf *bsdf, const ShaderData *sd)
 357 {
 358         bsdf->extra->cspec0.x = saturate(bsdf->extra->cspec0.x);
 359         bsdf->extra->cspec0.y = saturate(bsdf->extra->cspec0.y);
 360         bsdf->extra->cspec0.z = saturate(bsdf->extra->cspec0.z);
 361
 362         float F0 = fresnel_dielectric_cos(1.0f, bsdf->ior);
 363         float F = average(interpolate_fresnel_color(sd->I, bsdf->N, bsdf->ior, F0, bsdf->extra->cspec0));
 364         bsdf->sample_weight *= F;
 365
 366         bsdf->alpha_x = saturate(bsdf->alpha_x);
 367         bsdf->alpha_y = saturate(bsdf->alpha_y);
 368
 369         bsdf->type = CLOSURE_BSDF_MICROFACET_GGX_ANISO_FRESNEL_ID;
 370
 371         return SD_BSDF|SD_BSDF_HAS_EVAL;
 372 }
 373
 374 ccl_device int bsdf_microfacet_ggx_refraction_setup(MicrofacetBsdf *bsdf)
 375 {
 376         bsdf->extra = NULL;
 377
 378         bsdf->alpha_x = saturate(bsdf->alpha_x);
 379         bsdf->alpha_y = bsdf->alpha_x;
 380
 381         bsdf->type = CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
 382
 383         return SD_BSDF|SD_BSDF_HAS_EVAL;
 384 }
 385
 386 ccl_device void bsdf_microfacet_ggx_blur(ShaderClosure *sc, float roughness)
 387 {
 388         MicrofacetBsdf *bsdf = (MicrofacetBsdf*)sc;
 389
 390         bsdf->alpha_x = fmaxf(roughness, bsdf->alpha_x);
 391         bsdf->alpha_y = fmaxf(roughness, bsdf->alpha_y);
 392 }
 393
 394 ccl_device float3 bsdf_microfacet_ggx_eval_reflect(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
 395 {
 396         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
 397         float alpha_x = bsdf->alpha_x;
 398         float alpha_y = bsdf->alpha_y;
 399         bool m_refractive = bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
 400         float3 N = bsdf->N;
 401
 402         if(m_refractive || alpha_x*alpha_y <= 1e-7f)
 403                 return make_float3(0.0f, 0.0f, 0.0f);
 404
 405         float cosNO = dot(N, I);
 406         float cosNI = dot(N, omega_in);
 407
 408         if(cosNI > 0 && cosNO > 0) {
 409                 /* get half vector */
 410                 float3 m = normalize(omega_in + I);
 411                 float alpha2 = alpha_x * alpha_y;
 412                 float D, G1o, G1i;
 413
 414                 if(alpha_x == alpha_y) {
 415                         /* isotropic
 416                          * eq. 20: (F*G*D)/(4*in*on)
 417                          * eq. 33: first we calculate D(m) */
 418                         float cosThetaM = dot(N, m);
 419                         float cosThetaM2 = cosThetaM * cosThetaM;
 420                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
 421                         float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
 422
 423                         if(bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_CLEARCOAT_ID) {
 424                                 /* use GTR1 for clearcoat */
 425                                 D = D_GTR1(cosThetaM, bsdf->alpha_x);
 426
 427                                 /* the alpha value for clearcoat is a fixed 0.25 => alpha2 = 0.25 * 0.25 */
 428                                 alpha2 = 0.0625f;
 429                         }
 430                         else {
 431                                 /* use GTR2 otherwise */
 432                                 D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
 433                         }
 434
 435                         /* eq. 34: now calculate G1(i,m) and G1(o,m) */
 436                         G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
 437                         G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI)));
 438                 }
 439                 else {
 440                         /* anisotropic */
 441                         float3 X, Y, Z = N;
 442                         make_orthonormals_tangent(Z, bsdf->T, &X, &Y);
 443
 444                         /* distribution */
 445                         float3 local_m = make_float3(dot(X, m), dot(Y, m), dot(Z, m));
 446                         float slope_x = -local_m.x/(local_m.z*alpha_x);
 447                         float slope_y = -local_m.y/(local_m.z*alpha_y);
 448                         float slope_len = 1 + slope_x*slope_x + slope_y*slope_y;
 449
 450                         float cosThetaM = local_m.z;
 451                         float cosThetaM2 = cosThetaM * cosThetaM;
 452                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
 453
 454                         D = 1 / ((slope_len * slope_len) * M_PI_F * alpha2 * cosThetaM4);
 455
 456                         /* G1(i,m) and G1(o,m) */
 457                         float tanThetaO2 = (1 - cosNO * cosNO) / (cosNO * cosNO);
 458                         float cosPhiO = dot(I, X);
 459                         float sinPhiO = dot(I, Y);
 460
 461                         float alphaO2 = (cosPhiO*cosPhiO)*(alpha_x*alpha_x) + (sinPhiO*sinPhiO)*(alpha_y*alpha_y);
 462                         alphaO2 /= cosPhiO*cosPhiO + sinPhiO*sinPhiO;
 463
 464                         G1o = 2 / (1 + safe_sqrtf(1 + alphaO2 * tanThetaO2));
 465
 466                         float tanThetaI2 = (1 - cosNI * cosNI) / (cosNI * cosNI);
 467                         float cosPhiI = dot(omega_in, X);
 468                         float sinPhiI = dot(omega_in, Y);
 469
 470                         float alphaI2 = (cosPhiI*cosPhiI)*(alpha_x*alpha_x) + (sinPhiI*sinPhiI)*(alpha_y*alpha_y);
 471                         alphaI2 /= cosPhiI*cosPhiI + sinPhiI*sinPhiI;
 472
 473                         G1i = 2 / (1 + safe_sqrtf(1 + alphaI2 * tanThetaI2));
 474                 }
 475
 476                 float G = G1o * G1i;
 477
 478                 /* eq. 20 */
 479                 float common = D * 0.25f / cosNO;
 480
 481                 float3 F = reflection_color(bsdf, omega_in, m);
 482                 if(bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_CLEARCOAT_ID) {
 483                         F *= 0.25f * bsdf->extra->clearcoat;
 484                 }
 485
 486                 float3 out = F * G * common;
 487
 488                 /* eq. 2 in distribution of visible normals sampling
 489                  * pm = Dw = G1o * dot(m, I) * D / dot(N, I); */
 490
 491                 /* eq. 38 - but see also:
 492                  * eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
 493                  * pdf = pm * 0.25 / dot(m, I); */
 494                 *pdf = G1o * common;
 495
 496                 return out;
 497         }
 498
 499         return make_float3(0.0f, 0.0f, 0.0f);
 500 }
 501
 502 ccl_device float3 bsdf_microfacet_ggx_eval_transmit(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
 503 {
 504         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
 505         float alpha_x = bsdf->alpha_x;
 506         float alpha_y = bsdf->alpha_y;
 507         float m_eta = bsdf->ior;
 508         bool m_refractive = bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
 509         float3 N = bsdf->N;
 510
 511         if(!m_refractive || alpha_x*alpha_y <= 1e-7f)
 512                 return make_float3(0.0f, 0.0f, 0.0f);
 513
 514         float cosNO = dot(N, I);
 515         float cosNI = dot(N, omega_in);
 516
 517         if(cosNO <= 0 || cosNI >= 0)
 518                 return make_float3(0.0f, 0.0f, 0.0f); /* vectors on same side -- not possible */
 519
 520         /* compute half-vector of the refraction (eq. 16) */
 521         float3 ht = -(m_eta * omega_in + I);
 522         float3 Ht = normalize(ht);
 523         float cosHO = dot(Ht, I);
 524         float cosHI = dot(Ht, omega_in);
 525
 526         float D, G1o, G1i;
 527
 528         /* eq. 33: first we calculate D(m) with m=Ht: */
 529         float alpha2 = alpha_x * alpha_y;
 530         float cosThetaM = dot(N, Ht);
 531         float cosThetaM2 = cosThetaM * cosThetaM;
 532         float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
 533         float cosThetaM4 = cosThetaM2 * cosThetaM2;
 534         D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
 535
 536         /* eq. 34: now calculate G1(i,m) and G1(o,m) */
 537         G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
 538         G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI)));
 539
 540         float G = G1o * G1i;
 541
 542         /* probability */
 543         float Ht2 = dot(ht, ht);
 544
 545         /* eq. 2 in distribution of visible normals sampling
 546          * pm = Dw = G1o * dot(m, I) * D / dot(N, I); */
 547
 548         /* out = fabsf(cosHI * cosHO) * (m_eta * m_eta) * G * D / (cosNO * Ht2)
 549          * pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2 */
 550         float common = D * (m_eta * m_eta) / (cosNO * Ht2);
 551         float out = G * fabsf(cosHI * cosHO) * common;
 552         *pdf = G1o * fabsf(cosHO * cosHI) * common;
 553
 554         return make_float3(out, out, out);
 555 }
 556
 557 ccl_device int bsdf_microfacet_ggx_sample(KernelGlobals *kg, 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)
 558 {
 559         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
 560         float alpha_x = bsdf->alpha_x;
 561         float alpha_y = bsdf->alpha_y;
 562         bool m_refractive = bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
 563         float3 N = bsdf->N;
 564         int label;
 565
 566         float cosNO = dot(N, I);
 567         if(cosNO > 0) {
 568                 float3 X, Y, Z = N;
 569
 570                 if(alpha_x == alpha_y)
 571                         make_orthonormals(Z, &X, &Y);
 572                 else
 573                         make_orthonormals_tangent(Z, bsdf->T, &X, &Y);
 574
 575                 /* importance sampling with distribution of visible normals. vectors are
 576                  * transformed to local space before and after */
 577                 float3 local_I = make_float3(dot(X, I), dot(Y, I), cosNO);
 578                 float3 local_m;
 579                 float G1o;
 580
 581                 local_m = microfacet_sample_stretched(kg, local_I, alpha_x, alpha_y,
 582                         randu, randv, false, &G1o);
 583
 584                 float3 m = X*local_m.x + Y*local_m.y + Z*local_m.z;
 585                 float cosThetaM = local_m.z;
 586
 587                 /* reflection or refraction? */
 588                 if(!m_refractive) {
 589                         float cosMO = dot(m, I);
 590                         label = LABEL_REFLECT | LABEL_GLOSSY;
 591
 592                         if(cosMO > 0) {
 593                                 /* eq. 39 - compute actual reflected direction */
 594                                 *omega_in = 2 * cosMO * m - I;
 595
 596                                 if(dot(Ng, *omega_in) > 0) {
 597                                         if(alpha_x*alpha_y <= 1e-7f) {
 598                                                 /* some high number for MIS */
 599                                                 *pdf = 1e6f;
 600                                                 *eval = make_float3(1e6f, 1e6f, 1e6f);
 601
 602                                                 bool use_fresnel = (bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_FRESNEL_ID
 603                                                                    || bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_CLEARCOAT_ID
 604                                                                    || bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_ANISO_FRESNEL_ID);
 605
 606                                                 /* if fresnel is used, calculate the color with reflection_color(...) */
 607                                                 if(use_fresnel) {
 608                                                         *eval *= reflection_color(bsdf, *omega_in, m);
 609                                                 }
 610
 611                                                 label = LABEL_REFLECT | LABEL_SINGULAR;
 612                                         }
 613                                         else {
 614                                                 /* microfacet normal is visible to this ray */
 615                                                 /* eq. 33 */
 616                                                 float alpha2 = alpha_x * alpha_y;
 617                                                 float D, G1i;
 618
 619                                                 if(alpha_x == alpha_y) {
 620                                                         /* isotropic */
 621                                                         float cosThetaM2 = cosThetaM * cosThetaM;
 622                                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
 623                                                         float tanThetaM2 = 1/(cosThetaM2) - 1;
 624
 625                                                         /* eval BRDF*cosNI */
 626                                                         float cosNI = dot(N, *omega_in);
 627
 628                                                         if(bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_CLEARCOAT_ID) {
 629                                                                 /* use GTR1 for clearcoat */
 630                                                                 D = D_GTR1(cosThetaM, bsdf->alpha_x);
 631
 632                                                                 /* the alpha value for clearcoat is a fixed 0.25 => alpha2 = 0.25 * 0.25 */
 633                                                                 alpha2 = 0.0625f;
 634
 635                                                                 /* recalculate G1o */
 636                                                                 G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
 637                                                         }
 638                                                         else {
 639                                                                 /* use GTR2 otherwise */
 640                                                                 D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
 641                                                         }
 642
 643                                                         /* eq. 34: now calculate G1(i,m) */
 644                                                         G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI)));
 645                                                 }
 646                                                 else {
 647                                                         /* anisotropic distribution */
 648                                                         float3 local_m = make_float3(dot(X, m), dot(Y, m), dot(Z, m));
 649                                                         float slope_x = -local_m.x/(local_m.z*alpha_x);
 650                                                         float slope_y = -local_m.y/(local_m.z*alpha_y);
 651                                                         float slope_len = 1 + slope_x*slope_x + slope_y*slope_y;
 652
 653                                                         float cosThetaM = local_m.z;
 654                                                         float cosThetaM2 = cosThetaM * cosThetaM;
 655                                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
 656
 657                                                         D = 1 / ((slope_len * slope_len) * M_PI_F * alpha2 * cosThetaM4);
 658
 659                                                         /* calculate G1(i,m) */
 660                                                         float cosNI = dot(N, *omega_in);
 661
 662                                                         float tanThetaI2 = (1 - cosNI * cosNI) / (cosNI * cosNI);
 663                                                         float cosPhiI = dot(*omega_in, X);
 664                                                         float sinPhiI = dot(*omega_in, Y);
 665
 666                                                         float alphaI2 = (cosPhiI*cosPhiI)*(alpha_x*alpha_x) + (sinPhiI*sinPhiI)*(alpha_y*alpha_y);
 667                                                         alphaI2 /= cosPhiI*cosPhiI + sinPhiI*sinPhiI;
 668
 669                                                         G1i = 2 / (1 + safe_sqrtf(1 + alphaI2 * tanThetaI2));
 670                                                 }
 671
 672                                                 /* see eval function for derivation */
 673                                                 float common = (G1o * D) * 0.25f / cosNO;
 674                                                 *pdf = common;
 675
 676                                                 float3 F = reflection_color(bsdf, *omega_in, m);
 677
 678                                                 *eval = G1i * common * F;
 679                                         }
 680
 681                                         if(bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_CLEARCOAT_ID) {
 682                                                 *eval *= 0.25f * bsdf->extra->clearcoat;
 683                                         }
 684
 685 #ifdef __RAY_DIFFERENTIALS__
 686                                         *domega_in_dx = (2 * dot(m, dIdx)) * m - dIdx;
 687                                         *domega_in_dy = (2 * dot(m, dIdy)) * m - dIdy;
 688 #endif
 689                                 }
 690                         }
 691                 }
 692                 else {
 693                         label = LABEL_TRANSMIT | LABEL_GLOSSY;
 694
 695                         /* CAUTION: the i and o variables are inverted relative to the paper
 696                          * eq. 39 - compute actual refractive direction */
 697                         float3 R, T;
 698 #ifdef __RAY_DIFFERENTIALS__
 699                         float3 dRdx, dRdy, dTdx, dTdy;
 700 #endif
 701                         float m_eta = bsdf->ior, fresnel;
 702                         bool inside;
 703
 704                         fresnel = fresnel_dielectric(m_eta, m, I, &R, &T,
 705 #ifdef __RAY_DIFFERENTIALS__
 706                                 dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,
 707 #endif
 708                                 &inside);
 709
 710                         if(!inside && fresnel != 1.0f) {
 711
 712                                 *omega_in = T;
 713 #ifdef __RAY_DIFFERENTIALS__
 714                                 *domega_in_dx = dTdx;
 715                                 *domega_in_dy = dTdy;
 716 #endif
 717
 718                                 if(alpha_x*alpha_y <= 1e-7f || fabsf(m_eta - 1.0f) < 1e-4f) {
 719                                         /* some high number for MIS */
 720                                         *pdf = 1e6f;
 721                                         *eval = make_float3(1e6f, 1e6f, 1e6f);
 722                                         label = LABEL_TRANSMIT | LABEL_SINGULAR;
 723                                 }
 724                                 else {
 725                                         /* eq. 33 */
 726                                         float alpha2 = alpha_x * alpha_y;
 727                                         float cosThetaM2 = cosThetaM * cosThetaM;
 728                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
 729                                         float tanThetaM2 = 1/(cosThetaM2) - 1;
 730                                         float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
 731
 732                                         /* eval BRDF*cosNI */
 733                                         float cosNI = dot(N, *omega_in);
 734
 735                                         /* eq. 34: now calculate G1(i,m) */
 736                                         float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI)));
 737
 738                                         /* eq. 21 */
 739                                         float cosHI = dot(m, *omega_in);
 740                                         float cosHO = dot(m, I);
 741                                         float Ht2 = m_eta * cosHI + cosHO;
 742                                         Ht2 *= Ht2;
 743
 744                                         /* see eval function for derivation */
 745                                         float common = (G1o * D) * (m_eta * m_eta) / (cosNO * Ht2);
 746                                         float out = G1i * fabsf(cosHI * cosHO) * common;
 747                                         *pdf = cosHO * fabsf(cosHI) * common;
 748
 749                                         *eval = make_float3(out, out, out);
 750                                 }
 751                         }
 752                 }
 753         }
 754         else {
 755                 label = (m_refractive) ? LABEL_TRANSMIT|LABEL_GLOSSY : LABEL_REFLECT|LABEL_GLOSSY;
 756         }
 757         return label;
 758 }
 759
 760 /* Beckmann microfacet with Smith shadow-masking from:
 761  *
 762  * Microfacet Models for Refraction through Rough Surfaces
 763  * B. Walter, S. R. Marschner, H. Li, K. E. Torrance, EGSR 2007 */
 764
 765 ccl_device int bsdf_microfacet_beckmann_setup(MicrofacetBsdf *bsdf)
 766 {
 767         bsdf->alpha_x = saturate(bsdf->alpha_x);
 768         bsdf->alpha_y = bsdf->alpha_x;
 769
 770         bsdf->type = CLOSURE_BSDF_MICROFACET_BECKMANN_ID;
 771         return SD_BSDF|SD_BSDF_HAS_EVAL;
 772 }
 773
 774 ccl_device int bsdf_microfacet_beckmann_aniso_setup(MicrofacetBsdf *bsdf)
 775 {
 776         bsdf->alpha_x = saturate(bsdf->alpha_x);
 777         bsdf->alpha_y = saturate(bsdf->alpha_y);
 778
 779         bsdf->type = CLOSURE_BSDF_MICROFACET_BECKMANN_ANISO_ID;
 780         return SD_BSDF|SD_BSDF_HAS_EVAL;
 781 }
 782
 783 ccl_device int bsdf_microfacet_beckmann_refraction_setup(MicrofacetBsdf *bsdf)
 784 {
 785         bsdf->alpha_x = saturate(bsdf->alpha_x);
 786         bsdf->alpha_y = bsdf->alpha_x;
 787
 788         bsdf->type = CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
 789         return SD_BSDF|SD_BSDF_HAS_EVAL;
 790 }
 791
 792 ccl_device void bsdf_microfacet_beckmann_blur(ShaderClosure *sc, float roughness)
 793 {
 794         MicrofacetBsdf *bsdf = (MicrofacetBsdf*)sc;
 795
 796         bsdf->alpha_x = fmaxf(roughness, bsdf->alpha_x);
 797         bsdf->alpha_y = fmaxf(roughness, bsdf->alpha_y);
 798 }
 799
 800 ccl_device_inline float bsdf_beckmann_G1(float alpha, float cos_n)
 801 {
 802         cos_n *= cos_n;
 803         float invA = alpha * safe_sqrtf((1.0f - cos_n) / cos_n);
 804         if(invA < 0.625f) {
 805                 return 1.0f;
 806         }
 807
 808         float a = 1.0f / invA;
 809         return ((2.181f*a + 3.535f)*a) / ((2.577f*a + 2.276f)*a + 1.0f);
 810 }
 811
 812 ccl_device_inline float bsdf_beckmann_aniso_G1(float alpha_x, float alpha_y, float cos_n, float cos_phi, float sin_phi)
 813 {
 814         cos_n *= cos_n;
 815         sin_phi *= sin_phi;
 816         cos_phi *= cos_phi;
 817         alpha_x *= alpha_x;
 818         alpha_y *= alpha_y;
 819
 820         float alphaO2 = (cos_phi*alpha_x + sin_phi*alpha_y) / (cos_phi + sin_phi);
 821         float invA = safe_sqrtf(alphaO2 * (1 - cos_n) / cos_n);
 822         if(invA < 0.625f) {
 823                 return 1.0f;
 824         }
 825
 826         float a = 1.0f / invA;
 827         return ((2.181f*a + 3.535f)*a) / ((2.577f*a + 2.276f)*a + 1.0f);
 828 }
 829
 830 ccl_device float3 bsdf_microfacet_beckmann_eval_reflect(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
 831 {
 832         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
 833         float alpha_x = bsdf->alpha_x;
 834         float alpha_y = bsdf->alpha_y;
 835         bool m_refractive = bsdf->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
 836         float3 N = bsdf->N;
 837
 838         if(m_refractive || alpha_x*alpha_y <= 1e-7f)
 839                 return make_float3(0.0f, 0.0f, 0.0f);
 840
 841         float cosNO = dot(N, I);
 842         float cosNI = dot(N, omega_in);
 843
 844         if(cosNO > 0 && cosNI > 0) {
 845                 /* get half vector */
 846                 float3 m = normalize(omega_in + I);
 847
 848                 float alpha2 = alpha_x * alpha_y;
 849                 float D, G1o, G1i;
 850
 851                 if(alpha_x == alpha_y) {
 852                         /* isotropic
 853                          * eq. 20: (F*G*D)/(4*in*on)
 854                          * eq. 25: first we calculate D(m) */
 855                         float cosThetaM = dot(N, m);
 856                         float cosThetaM2 = cosThetaM * cosThetaM;
 857                         float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
 858                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
 859                         D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 * cosThetaM4);
 860
 861                         /* eq. 26, 27: now calculate G1(i,m) and G1(o,m) */
 862                         G1o = bsdf_beckmann_G1(alpha_x, cosNO);
 863                         G1i = bsdf_beckmann_G1(alpha_x, cosNI);
 864                 }
 865                 else {
 866                         /* anisotropic */
 867                         float3 X, Y, Z = N;
 868                         make_orthonormals_tangent(Z, bsdf->T, &X, &Y);
 869
 870                         /* distribution */
 871                         float3 local_m = make_float3(dot(X, m), dot(Y, m), dot(Z, m));
 872                         float slope_x = -local_m.x/(local_m.z*alpha_x);
 873                         float slope_y = -local_m.y/(local_m.z*alpha_y);
 874
 875                         float cosThetaM = local_m.z;
 876                         float cosThetaM2 = cosThetaM * cosThetaM;
 877                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
 878
 879                         D = expf(-slope_x*slope_x - slope_y*slope_y) / (M_PI_F * alpha2 * cosThetaM4);
 880
 881                         /* G1(i,m) and G1(o,m) */
 882                         G1o = bsdf_beckmann_aniso_G1(alpha_x, alpha_y, cosNO, dot(I, X), dot(I, Y));
 883                         G1i = bsdf_beckmann_aniso_G1(alpha_x, alpha_y, cosNI, dot(omega_in, X), dot(omega_in, Y));
 884                 }
 885
 886                 float G = G1o * G1i;
 887
 888                 /* eq. 20 */
 889                 float common = D * 0.25f / cosNO;
 890                 float out = G * common;
 891
 892                 /* eq. 2 in distribution of visible normals sampling
 893                  * pm = Dw = G1o * dot(m, I) * D / dot(N, I); */
 894
 895                 /* eq. 38 - but see also:
 896                  * eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
 897                  * pdf = pm * 0.25 / dot(m, I); */
 898                 *pdf = G1o * common;
 899
 900                 return make_float3(out, out, out);
 901         }
 902
 903         return make_float3(0.0f, 0.0f, 0.0f);
 904 }
 905
 906 ccl_device float3 bsdf_microfacet_beckmann_eval_transmit(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
 907 {
 908         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
 909         float alpha_x = bsdf->alpha_x;
 910         float alpha_y = bsdf->alpha_y;
 911         float m_eta = bsdf->ior;
 912         bool m_refractive = bsdf->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
 913         float3 N = bsdf->N;
 914
 915         if(!m_refractive || alpha_x*alpha_y <= 1e-7f)
 916                 return make_float3(0.0f, 0.0f, 0.0f);
 917
 918         float cosNO = dot(N, I);
 919         float cosNI = dot(N, omega_in);
 920
 921         if(cosNO <= 0 || cosNI >= 0)
 922                 return make_float3(0.0f, 0.0f, 0.0f);
 923
 924         /* compute half-vector of the refraction (eq. 16) */
 925         float3 ht = -(m_eta * omega_in + I);
 926         float3 Ht = normalize(ht);
 927         float cosHO = dot(Ht, I);
 928         float cosHI = dot(Ht, omega_in);
 929
 930         /* eq. 25: first we calculate D(m) with m=Ht: */
 931         float alpha2 = alpha_x * alpha_y;
 932         float cosThetaM = min(dot(N, Ht), 1.0f);
 933         float cosThetaM2 = cosThetaM * cosThetaM;
 934         float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
 935         float cosThetaM4 = cosThetaM2 * cosThetaM2;
 936         float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
 937
 938         /* eq. 26, 27: now calculate G1(i,m) and G1(o,m) */
 939         float G1o = bsdf_beckmann_G1(alpha_x, cosNO);
 940         float G1i = bsdf_beckmann_G1(alpha_x, cosNI);
 941         float G = G1o * G1i;
 942
 943         /* probability */
 944         float Ht2 = dot(ht, ht);
 945
 946         /* eq. 2 in distribution of visible normals sampling
 947          * pm = Dw = G1o * dot(m, I) * D / dot(N, I); */
 948
 949         /* out = fabsf(cosHI * cosHO) * (m_eta * m_eta) * G * D / (cosNO * Ht2)
 950          * pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2 */
 951         float common = D * (m_eta * m_eta) / (cosNO * Ht2);
 952         float out = G * fabsf(cosHI * cosHO) * common;
 953         *pdf = G1o * fabsf(cosHO * cosHI) * common;
 954
 955         return make_float3(out, out, out);
 956 }
 957
 958 ccl_device int bsdf_microfacet_beckmann_sample(KernelGlobals *kg, 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)
 959 {
 960         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
 961         float alpha_x = bsdf->alpha_x;
 962         float alpha_y = bsdf->alpha_y;
 963         bool m_refractive = bsdf->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
 964         float3 N = bsdf->N;
 965         int label;
 966
 967         float cosNO = dot(N, I);
 968         if(cosNO > 0) {
 969                 float3 X, Y, Z = N;
 970
 971                 if(alpha_x == alpha_y)
 972                         make_orthonormals(Z, &X, &Y);
 973                 else
 974                         make_orthonormals_tangent(Z, bsdf->T, &X, &Y);
 975
 976                 /* importance sampling with distribution of visible normals. vectors are
 977                  * transformed to local space before and after */
 978                 float3 local_I = make_float3(dot(X, I), dot(Y, I), cosNO);
 979                 float3 local_m;
 980                 float G1o;
 981
 982                 local_m = microfacet_sample_stretched(kg, local_I, alpha_x, alpha_x,
 983                         randu, randv, true, &G1o);
 984
 985                 float3 m = X*local_m.x + Y*local_m.y + Z*local_m.z;
 986                 float cosThetaM = local_m.z;
 987
 988                 /* reflection or refraction? */
 989                 if(!m_refractive) {
 990                         label = LABEL_REFLECT | LABEL_GLOSSY;
 991                         float cosMO = dot(m, I);
 992
 993                         if(cosMO > 0) {
 994                                 /* eq. 39 - compute actual reflected direction */
 995                                 *omega_in = 2 * cosMO * m - I;
 996
 997                                 if(dot(Ng, *omega_in) > 0) {
 998                                         if(alpha_x*alpha_y <= 1e-7f) {
 999                                                 /* some high number for MIS */
1000                                                 *pdf = 1e6f;
1001                                                 *eval = make_float3(1e6f, 1e6f, 1e6f);
1002                                                 label = LABEL_REFLECT | LABEL_SINGULAR;
1003                                         }
1004                                         else {
1005                                                 /* microfacet normal is visible to this ray
1006                                                  * eq. 25 */
1007                                                 float alpha2 = alpha_x * alpha_y;
1008                                                 float D, G1i;
1009
1010                                                 if(alpha_x == alpha_y) {
1011                                                         /* istropic distribution */
1012                                                         float cosThetaM2 = cosThetaM * cosThetaM;
1013                                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
1014                                                         float tanThetaM2 = 1/(cosThetaM2) - 1;
1015                                                         D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
1016
1017                                                         /* eval BRDF*cosNI */
1018                                                         float cosNI = dot(N, *omega_in);
1019
1020                                                         /* eq. 26, 27: now calculate G1(i,m) */
1021                                                         G1i = bsdf_beckmann_G1(alpha_x, cosNI);
1022                                                 }
1023                                                 else {
1024                                                         /* anisotropic distribution */
1025                                                         float3 local_m = make_float3(dot(X, m), dot(Y, m), dot(Z, m));
1026                                                         float slope_x = -local_m.x/(local_m.z*alpha_x);
1027                                                         float slope_y = -local_m.y/(local_m.z*alpha_y);
1028
1029                                                         float cosThetaM = local_m.z;
1030                                                         float cosThetaM2 = cosThetaM * cosThetaM;
1031                                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
1032
1033                                                         D = expf(-slope_x*slope_x - slope_y*slope_y) / (M_PI_F * alpha2 * cosThetaM4);
1034
1035                                                         /* G1(i,m) */
1036                                                         G1i = bsdf_beckmann_aniso_G1(alpha_x, alpha_y, dot(*omega_in, N), dot(*omega_in, X), dot(*omega_in, Y));
1037                                                 }
1038
1039                                                 float G = G1o * G1i;
1040
1041                                                 /* see eval function for derivation */
1042                                                 float common = D * 0.25f / cosNO;
1043                                                 float out = G * common;
1044                                                 *pdf = G1o * common;
1045
1046                                                 *eval = make_float3(out, out, out);
1047                                         }
1048
1049 #ifdef __RAY_DIFFERENTIALS__
1050                                         *domega_in_dx = (2 * dot(m, dIdx)) * m - dIdx;
1051                                         *domega_in_dy = (2 * dot(m, dIdy)) * m - dIdy;
1052 #endif
1053                                 }
1054                         }
1055                 }
1056                 else {
1057                         label = LABEL_TRANSMIT | LABEL_GLOSSY;
1058
1059                         /* CAUTION: the i and o variables are inverted relative to the paper
1060                          * eq. 39 - compute actual refractive direction */
1061                         float3 R, T;
1062 #ifdef __RAY_DIFFERENTIALS__
1063                         float3 dRdx, dRdy, dTdx, dTdy;
1064 #endif
1065                         float m_eta = bsdf->ior, fresnel;
1066                         bool inside;
1067
1068                         fresnel = fresnel_dielectric(m_eta, m, I, &R, &T,
1069 #ifdef __RAY_DIFFERENTIALS__
1070                                 dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,
1071 #endif
1072                                 &inside);
1073
1074                         if(!inside && fresnel != 1.0f) {
1075                                 *omega_in = T;
1076
1077 #ifdef __RAY_DIFFERENTIALS__
1078                                 *domega_in_dx = dTdx;
1079                                 *domega_in_dy = dTdy;
1080 #endif
1081
1082                                 if(alpha_x*alpha_y <= 1e-7f || fabsf(m_eta - 1.0f) < 1e-4f) {
1083                                         /* some high number for MIS */
1084                                         *pdf = 1e6f;
1085                                         *eval = make_float3(1e6f, 1e6f, 1e6f);
1086                                         label = LABEL_TRANSMIT | LABEL_SINGULAR;
1087                                 }
1088                                 else {
1089                                         /* eq. 33 */
1090                                         float alpha2 = alpha_x * alpha_y;
1091                                         float cosThetaM2 = cosThetaM * cosThetaM;
1092                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
1093                                         float tanThetaM2 = 1/(cosThetaM2) - 1;
1094                                         float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
1095
1096                                         /* eval BRDF*cosNI */
1097                                         float cosNI = dot(N, *omega_in);
1098
1099                                         /* eq. 26, 27: now calculate G1(i,m) */
1100                                         float G1i = bsdf_beckmann_G1(alpha_x, cosNI);
1101                                         float G = G1o * G1i;
1102
1103                                         /* eq. 21 */
1104                                         float cosHI = dot(m, *omega_in);
1105                                         float cosHO = dot(m, I);
1106                                         float Ht2 = m_eta * cosHI + cosHO;
1107                                         Ht2 *= Ht2;
1108
1109                                         /* see eval function for derivation */
1110                                         float common = D * (m_eta * m_eta) / (cosNO * Ht2);
1111                                         float out = G * fabsf(cosHI * cosHO) * common;
1112                                         *pdf = G1o * cosHO * fabsf(cosHI) * common;
1113
1114                                         *eval = make_float3(out, out, out);
1115                                 }
1116                         }
1117                 }
1118         }
1119         else {
1120                 label = (m_refractive) ? LABEL_TRANSMIT|LABEL_GLOSSY : LABEL_REFLECT|LABEL_GLOSSY;
1121         }
1122         return label;
1123 }
1124
1125 CCL_NAMESPACE_END
1126
1127 #endif /* __BSDF_MICROFACET_H__ */
1128