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