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