Cycles: better path termination for transparency.
[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         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