Cycles: improved importance sampling for Beckmann and GGX glossy
[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 ccl_device int bsdf_microfacet_ggx_setup(ShaderClosure *sc)
346 {
347         sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ag */
348         
349         sc->type = CLOSURE_BSDF_MICROFACET_GGX_ID;
350
351         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
352 }
353
354 ccl_device int bsdf_microfacet_ggx_refraction_setup(ShaderClosure *sc)
355 {
356         sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ag */
357
358         sc->type = CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
359
360         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
361 }
362
363 ccl_device void bsdf_microfacet_ggx_blur(ShaderClosure *sc, float roughness)
364 {
365         sc->data0 = fmaxf(roughness, sc->data0); /* m_ag */
366 }
367
368 ccl_device float3 bsdf_microfacet_ggx_eval_reflect(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
369 {
370         float m_ag = max(sc->data0, 1e-4f);
371         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
372         float3 N = sc->N;
373
374         if(m_refractive || m_ag <= 1e-4f)
375                 return make_float3(0, 0, 0);
376
377         float cosNO = dot(N, I);
378         float cosNI = dot(N, omega_in);
379
380         if(cosNI > 0 && cosNO > 0) {
381                 /* get half vector */
382                 float3 m = normalize(omega_in + I);
383
384                 /* eq. 20: (F*G*D)/(4*in*on)
385                  * eq. 33: first we calculate D(m) */
386                 float alpha2 = m_ag * m_ag;
387                 float cosThetaM = dot(N, m);
388                 float cosThetaM2 = cosThetaM * cosThetaM;
389                 float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
390                 float cosThetaM4 = cosThetaM2 * cosThetaM2;
391                 float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
392                 /* eq. 34: now calculate G1(i,m) and G1(o,m) */
393                 float G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
394                 float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
395                 float G = G1o * G1i;
396
397                 /* eq. 20 */
398                 float common = D * 0.25f / cosNO;
399                 float out = G * common;
400
401                 /* eq. 2 in distribution of visible normals sampling
402                  * pm = Dw = G1o * dot(m, I) * D / dot(N, I); */
403
404                 /* eq. 38 - but see also:
405                  * eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
406                  * pdf = pm * 0.25 / dot(m, I); */
407                 *pdf = G1o * common;
408
409                 return make_float3(out, out, out);
410         }
411
412         return make_float3(0, 0, 0);
413 }
414
415 ccl_device float3 bsdf_microfacet_ggx_eval_transmit(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
416 {
417         float m_ag = max(sc->data0, 1e-4f);
418         float m_eta = sc->data1;
419         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
420         float3 N = sc->N;
421
422         if(!m_refractive || m_ag <= 1e-4f)
423                 return make_float3(0, 0, 0);
424
425         float cosNO = dot(N, I);
426         float cosNI = dot(N, omega_in);
427
428         if(cosNO <= 0 || cosNI >= 0)
429                 return make_float3(0, 0, 0); /* vectors on same side -- not possible */
430
431         /* compute half-vector of the refraction (eq. 16) */
432         float3 ht = -(m_eta * omega_in + I);
433         float3 Ht = normalize(ht);
434         float cosHO = dot(Ht, I);
435         float cosHI = dot(Ht, omega_in);
436
437         /* eq. 33: first we calculate D(m) with m=Ht: */
438         float alpha2 = m_ag * m_ag;
439         float cosThetaM = dot(N, Ht);
440         float cosThetaM2 = cosThetaM * cosThetaM;
441         float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
442         float cosThetaM4 = cosThetaM2 * cosThetaM2;
443         float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
444
445         /* eq. 34: now calculate G1(i,m) and G1(o,m) */
446         float G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
447         float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
448         float G = G1o * G1i;
449
450         /* probability */
451         float Ht2 = dot(ht, ht);
452
453         /* eq. 2 in distribution of visible normals sampling
454          * pm = Dw = G1o * dot(m, I) * D / dot(N, I); */
455
456         /* out = fabsf(cosHI * cosHO) * (m_eta * m_eta) * G * D / (cosNO * Ht2)
457          * pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2 */
458         float common = D * (m_eta * m_eta) / (cosNO * Ht2);
459         float out = G * fabsf(cosHI * cosHO) * common;
460         *pdf = G1o * cosHO * fabsf(cosHI) * common;
461
462         return make_float3(out, out, out);
463 }
464
465 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)
466 {
467         float m_ag = sc->data0;
468         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
469         float3 N = sc->N;
470
471         float cosNO = dot(N, I);
472         if(cosNO > 0) {
473                 float3 X, Y, Z = N;
474                 make_orthonormals(Z, &X, &Y);
475
476                 /* importance sampling with distribution of visible normals. vectors are
477                  * transformed to local space before and after */
478                 float3 local_I;
479
480                 local_I.x = dot(X, I);
481                 local_I.y = dot(Y, I);
482                 local_I.z = cosNO;
483
484                 float3 local_m;
485                 float G1o;
486
487                 local_m = microfacet_sample_stretched(local_I, m_ag, m_ag,
488                         randu, randv, false, &G1o);
489
490                 float3 m = X*local_m.x + Y*local_m.y + Z*local_m.z;
491                 float cosThetaM = local_m.z;
492
493                 /* reflection or refraction? */
494                 if(!m_refractive) {
495                         float cosMO = dot(m, I);
496
497                         if(cosMO > 0) {
498                                 /* eq. 39 - compute actual reflected direction */
499                                 *omega_in = 2 * cosMO * m - I;
500
501                                 if(dot(Ng, *omega_in) > 0) {
502                                         if(m_ag <= 1e-4f) {
503                                                 /* some high number for MIS */
504                                                 *pdf = 1e6f;
505                                                 *eval = make_float3(1e6f, 1e6f, 1e6f);
506                                         }
507                                         else {
508                                                 /* microfacet normal is visible to this ray */
509                                                 /* eq. 33 */
510                                                 float alpha2 = m_ag * m_ag;
511                                                 float cosThetaM2 = cosThetaM * cosThetaM;
512                                                 float cosThetaM4 = cosThetaM2 * cosThetaM2;
513                                                 float tanThetaM2 = 1/(cosThetaM2) - 1;
514                                                 float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
515
516                                                 /* eval BRDF*cosNI */
517                                                 float cosNI = dot(N, *omega_in);
518                                                 /* eq. 34: now calculate G1(i,m) */
519                                                 float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
520
521                                                 /* see eval function for derivation */
522                                                 float common = (G1o * D) * 0.25f / cosNO;
523                                                 float out = G1i * common;
524                                                 *pdf = common;
525
526                                                 *eval = make_float3(out, out, out);
527                                         }
528
529 #ifdef __RAY_DIFFERENTIALS__
530                                         *domega_in_dx = (2 * dot(m, dIdx)) * m - dIdx;
531                                         *domega_in_dy = (2 * dot(m, dIdy)) * m - dIdy;
532 #endif
533                                 }
534                         }
535                 }
536                 else {
537                         /* CAUTION: the i and o variables are inverted relative to the paper
538                          * eq. 39 - compute actual refractive direction */
539                         float3 R, T;
540 #ifdef __RAY_DIFFERENTIALS__
541                         float3 dRdx, dRdy, dTdx, dTdy;
542 #endif
543                         float m_eta = sc->data1;
544                         bool inside;
545
546                         fresnel_dielectric(m_eta, m, I, &R, &T,
547 #ifdef __RAY_DIFFERENTIALS__
548                                 dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,
549 #endif
550                                 &inside);
551                         
552                         if(!inside) {
553
554                                 *omega_in = T;
555 #ifdef __RAY_DIFFERENTIALS__
556                                 *domega_in_dx = dTdx;
557                                 *domega_in_dy = dTdy;
558 #endif
559
560                                 if(m_ag <= 1e-4f || fabsf(m_eta - 1.0f) < 1e-4f) {
561                                         /* some high number for MIS */
562                                         *pdf = 1e6f;
563                                         *eval = make_float3(1e6f, 1e6f, 1e6f);
564                                 }
565                                 else {
566                                         /* eq. 33 */
567                                         float alpha2 = m_ag * m_ag;
568                                         float cosThetaM2 = cosThetaM * cosThetaM;
569                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
570                                         float tanThetaM2 = 1/(cosThetaM2) - 1;
571                                         float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
572
573                                         /* eval BRDF*cosNI */
574                                         float cosNI = dot(N, *omega_in);
575
576                                         /* eq. 34: now calculate G1(i,m) */
577                                         float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
578
579                                         /* eq. 21 */
580                                         float cosHI = dot(m, *omega_in);
581                                         float cosHO = dot(m, I);
582                                         float Ht2 = m_eta * cosHI + cosHO;
583                                         Ht2 *= Ht2;
584
585                                         /* see eval function for derivation */
586                                         float common = (G1o * D) * (m_eta * m_eta) / (cosNO * Ht2);
587                                         float out = G1i * fabsf(cosHI * cosHO) * common;
588                                         *pdf = cosHO * fabsf(cosHI) * common;
589
590                                         *eval = make_float3(out, out, out);
591                                 }
592                         }
593                 }
594         }
595         return (m_refractive) ? LABEL_TRANSMIT|LABEL_GLOSSY : LABEL_REFLECT|LABEL_GLOSSY;
596 }
597
598 /* Beckmann microfacet with Smith shadow-masking from:
599  *
600  * Microfacet Models for Refraction through Rough Surfaces
601  * B. Walter, S. R. Marschner, H. Li, K. E. Torrance, EGSR 2007 */
602
603 ccl_device int bsdf_microfacet_beckmann_setup(ShaderClosure *sc)
604 {
605         sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ab */
606
607         sc->type = CLOSURE_BSDF_MICROFACET_BECKMANN_ID;
608         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
609 }
610
611 ccl_device int bsdf_microfacet_beckmann_refraction_setup(ShaderClosure *sc)
612 {
613         sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ab */
614
615         sc->type = CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
616         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
617 }
618
619 ccl_device void bsdf_microfacet_beckmann_blur(ShaderClosure *sc, float roughness)
620 {
621         sc->data0 = fmaxf(roughness, sc->data0); /* m_ab */
622 }
623
624 ccl_device float3 bsdf_microfacet_beckmann_eval_reflect(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
625 {
626         float m_ab = max(sc->data0, 1e-4f);
627         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
628         float3 N = sc->N;
629
630         if(m_refractive || m_ab <= 1e-4f)
631                 return make_float3(0, 0, 0);
632
633         float cosNO = dot(N, I);
634         float cosNI = dot(N, omega_in);
635
636         if(cosNO > 0 && cosNI > 0) {
637            /* get half vector */
638            float3 m = normalize(omega_in + I);
639
640            /* eq. 20: (F*G*D)/(4*in*on)
641             * eq. 25: first we calculate D(m) */
642            float alpha2 = m_ab * m_ab;
643            float cosThetaM = dot(N, m);
644            float cosThetaM2 = cosThetaM * cosThetaM;
645            float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
646            float cosThetaM4 = cosThetaM2 * cosThetaM2;
647            float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
648
649            /* eq. 26, 27: now calculate G1(i,m) and G1(o,m) */
650            float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
651            float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
652            float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
653            float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
654            float G = G1o * G1i;
655
656                 /* eq. 20 */
657                 float common = D * 0.25f / cosNO;
658                 float out = G * common;
659
660                 /* eq. 2 in distribution of visible normals sampling
661                  * pm = Dw = G1o * dot(m, I) * D / dot(N, I); */
662
663                 /* eq. 38 - but see also:
664                  * eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
665                  * pdf = pm * 0.25 / dot(m, I); */
666                 *pdf = G1o * common;
667
668                 return make_float3(out, out, out);
669         }
670
671         return make_float3(0, 0, 0);
672 }
673
674 ccl_device float3 bsdf_microfacet_beckmann_eval_transmit(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
675 {
676         float m_ab = max(sc->data0, 1e-4f);
677         float m_eta = sc->data1;
678         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
679         float3 N = sc->N;
680
681         if(!m_refractive || m_ab <= 1e-4f)
682                 return make_float3(0, 0, 0);
683
684         float cosNO = dot(N, I);
685         float cosNI = dot(N, omega_in);
686
687         if(cosNO <= 0 || cosNI >= 0)
688                 return make_float3(0, 0, 0);
689
690         /* compute half-vector of the refraction (eq. 16) */
691         float3 ht = -(m_eta * omega_in + I);
692         float3 Ht = normalize(ht);
693         float cosHO = dot(Ht, I);
694         float cosHI = dot(Ht, omega_in);
695
696         /* eq. 33: first we calculate D(m) with m=Ht: */
697         float alpha2 = m_ab * m_ab;
698         float cosThetaM = min(dot(N, Ht), 1.0f);
699         float cosThetaM2 = cosThetaM * cosThetaM;
700         float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
701         float cosThetaM4 = cosThetaM2 * cosThetaM2;
702         float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
703
704         /* eq. 26, 27: now calculate G1(i,m) and G1(o,m) */
705         float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
706         float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
707         float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
708         float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
709         float G = G1o * G1i;
710
711         /* probability */
712         float Ht2 = dot(ht, ht);
713
714         /* eq. 2 in distribution of visible normals sampling
715          * pm = Dw = G1o * dot(m, I) * D / dot(N, I); */
716
717         /* out = fabsf(cosHI * cosHO) * (m_eta * m_eta) * G * D / (cosNO * Ht2)
718          * pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2 */
719         float common = D * (m_eta * m_eta) / (cosNO * Ht2);
720         float out = G * fabsf(cosHI * cosHO) * common;
721         *pdf = G1o * cosHO * fabsf(cosHI) * common;
722
723         return make_float3(out, out, out);
724 }
725
726 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)
727 {
728         float m_ab = sc->data0;
729         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
730         float3 N = sc->N;
731
732         float cosNO = dot(N, I);
733         if(cosNO > 0) {
734                 float3 X, Y, Z = N;
735                 make_orthonormals(Z, &X, &Y);
736
737                 /* importance sampling with distribution of visible normals. vectors are
738                  * transformed to local space before and after */
739                 float3 local_I;
740
741                 local_I.x = dot(X, I);
742                 local_I.y = dot(Y, I);
743                 local_I.z = cosNO;
744
745                 float3 local_m;
746                 float G1o;
747
748                 local_m = microfacet_sample_stretched(local_I, m_ab, m_ab,
749                         randu, randv, true, &G1o);
750
751                 float3 m = X*local_m.x + Y*local_m.y + Z*local_m.z;
752                 float cosThetaM = local_m.z;
753
754                 /* reflection or refraction? */
755                 if(!m_refractive) {
756                         float cosMO = dot(m, I);
757
758                         if(cosMO > 0) {
759                                 /* eq. 39 - compute actual reflected direction */
760                                 *omega_in = 2 * cosMO * m - I;
761
762                                 if(dot(Ng, *omega_in) > 0) {
763                                         if(m_ab <= 1e-4f) {
764                                                 /* some high number for MIS */
765                                                 *pdf = 1e6f;
766                                                 *eval = make_float3(1e6f, 1e6f, 1e6f);
767                                         }
768                                         else {
769                                                 /* microfacet normal is visible to this ray
770                                                  * eq. 25 */
771                                                 float alpha2 = m_ab * m_ab;
772                                                 float cosThetaM2 = cosThetaM * cosThetaM;
773                                                 float cosThetaM4 = cosThetaM2 * cosThetaM2;
774                                                 float tanThetaM2 = 1/(cosThetaM2) - 1;
775                                                 float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
776
777                                                 /* eval BRDF*cosNI */
778                                                 float cosNI = dot(N, *omega_in);
779
780                                                 /* eq. 26, 27: now calculate G1(i,m) */
781                                                 float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
782                                                 float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
783                                                 float G = G1o * G1i;
784
785                                                 /* see eval function for derivation */
786                                                 float common = D * 0.25f / cosNO;
787                                                 float out = G * common;
788                                                 *pdf = G1o * common;
789
790                                                 *eval = make_float3(out, out, out);
791                                         }
792
793 #ifdef __RAY_DIFFERENTIALS__
794                                         *domega_in_dx = (2 * dot(m, dIdx)) * m - dIdx;
795                                         *domega_in_dy = (2 * dot(m, dIdy)) * m - dIdy;
796 #endif
797                                 }
798                         }
799                 }
800                 else {
801                         /* CAUTION: the i and o variables are inverted relative to the paper
802                          * eq. 39 - compute actual refractive direction */
803                         float3 R, T;
804 #ifdef __RAY_DIFFERENTIALS__
805                         float3 dRdx, dRdy, dTdx, dTdy;
806 #endif
807                         float m_eta = sc->data1;
808                         bool inside;
809
810                         fresnel_dielectric(m_eta, m, I, &R, &T,
811 #ifdef __RAY_DIFFERENTIALS__
812                                 dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,
813 #endif
814                                 &inside);
815
816                         if(!inside) {
817                                 *omega_in = T;
818
819 #ifdef __RAY_DIFFERENTIALS__
820                                 *domega_in_dx = dTdx;
821                                 *domega_in_dy = dTdy;
822 #endif
823
824                                 if(m_ab <= 1e-4f || fabsf(m_eta - 1.0f) < 1e-4f) {
825                                         /* some high number for MIS */
826                                         *pdf = 1e6f;
827                                         *eval = make_float3(1e6f, 1e6f, 1e6f);
828                                 }
829                                 else {
830                                         /* eq. 33 */
831                                         float alpha2 = m_ab * m_ab;
832                                         float cosThetaM2 = cosThetaM * cosThetaM;
833                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
834                                         float tanThetaM2 = 1/(cosThetaM2) - 1;
835                                         float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
836
837                                         /* eval BRDF*cosNI */
838                                         float cosNI = dot(N, *omega_in);
839
840                                         /* eq. 26, 27: now calculate G1(i,m) */
841                                         float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
842                                         float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
843                                         float G = G1o * G1i;
844
845                                         /* eq. 21 */
846                                         float cosHI = dot(m, *omega_in);
847                                         float cosHO = dot(m, I);
848                                         float Ht2 = m_eta * cosHI + cosHO;
849                                         Ht2 *= Ht2;
850
851                                         /* see eval function for derivation */
852                                         float common = D * (m_eta * m_eta) / (cosNO * Ht2);
853                                         float out = G * fabsf(cosHI * cosHO) * common;
854                                         *pdf = G1o * cosHO * fabsf(cosHI) * common;
855
856                                         *eval = make_float3(out, out, out);
857                                 }
858                         }
859                 }
860         }
861         return (m_refractive) ? LABEL_TRANSMIT|LABEL_GLOSSY : LABEL_REFLECT|LABEL_GLOSSY;
862 }
863
864 CCL_NAMESPACE_END
865
866 #endif /* __BSDF_MICROFACET_H__ */
867