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