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