0a6dd4dcbdff722ee92c45bc67cfc3e43779768c
[blender-staging.git] / intern / cycles / kernel / closure / bsdf_microfacet_multi.h
1 /*
2  * Copyright 2011-2016 Blender Foundation
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16
17 CCL_NAMESPACE_BEGIN
18
19 /* Most of the code is based on the supplemental implementations from https://eheitzresearch.wordpress.com/240-2/. */
20
21 /* === GGX Microfacet distribution functions === */
22
23 /* Isotropic GGX microfacet distribution */
24 ccl_device_inline float D_ggx(float3 wm, float alpha)
25 {
26         wm.z *= wm.z;
27         alpha *= alpha;
28         float tmp = (1.0f - wm.z) + alpha * wm.z;
29         return alpha / max(M_PI_F * tmp*tmp, 1e-7f);
30 }
31
32 /* Anisotropic GGX microfacet distribution */
33 ccl_device_inline float D_ggx_aniso(const float3 wm, const float2 alpha)
34 {
35         float slope_x = -wm.x/alpha.x;
36         float slope_y = -wm.y/alpha.y;
37         float tmp = wm.z*wm.z + slope_x*slope_x + slope_y*slope_y;
38
39         return 1.0f / max(M_PI_F * tmp*tmp * alpha.x*alpha.y, 1e-7f);
40 }
41
42 /* Sample slope distribution (based on page 14 of the supplemental implementation). */
43 ccl_device_inline float2 mf_sampleP22_11(const float cosI, const float2 randU)
44 {
45         if(cosI > 0.9999f || cosI < 1e-6f) {
46                 const float r = sqrtf(randU.x / (1.0f - randU.x));
47                 const float phi = M_2PI_F * randU.y;
48                 return make_float2(r*cosf(phi), r*sinf(phi));
49         }
50
51         const float sinI = sqrtf(1.0f - cosI*cosI);
52         const float tanI = sinI/cosI;
53         const float projA = 0.5f * (cosI + 1.0f);
54         if(projA < 0.0001f)
55                 return make_float2(0.0f, 0.0f);
56         const float A = 2.0f*randU.x*projA / cosI - 1.0f;
57         float tmp = A*A-1.0f;
58         if(fabsf(tmp) < 1e-7f)
59                 return make_float2(0.0f, 0.0f);
60         tmp = 1.0f / tmp;
61         const float D = safe_sqrtf(tanI*tanI*tmp*tmp - (A*A-tanI*tanI)*tmp);
62
63         const float slopeX2 = tanI*tmp + D;
64         const float slopeX = (A < 0.0f || slopeX2 > 1.0f/tanI)? (tanI*tmp - D) : slopeX2;
65
66         float U2;
67         if(randU.y >= 0.5f)
68                 U2 = 2.0f*(randU.y - 0.5f);
69         else
70                 U2 = 2.0f*(0.5f - randU.y);
71         const float z = (U2*(U2*(U2*0.27385f-0.73369f)+0.46341f)) / (U2*(U2*(U2*0.093073f+0.309420f)-1.0f)+0.597999f);
72         const float slopeY = z * sqrtf(1.0f + slopeX*slopeX);
73
74         if(randU.y >= 0.5f)
75                 return make_float2(slopeX, slopeY);
76         else
77                 return make_float2(slopeX, -slopeY);
78 }
79
80 /* Visible normal sampling for the GGX distribution (based on page 7 of the supplemental implementation). */
81 ccl_device_inline float3 mf_sample_vndf(const float3 wi, const float2 alpha, const float2 randU)
82 {
83         const float3 wi_11 = normalize(make_float3(alpha.x*wi.x, alpha.y*wi.y, wi.z));
84         const float2 slope_11 = mf_sampleP22_11(wi_11.z, randU);
85
86         const float2 cossin_phi = normalize(make_float2(wi_11.x, wi_11.y));
87         const float slope_x = alpha.x*(cossin_phi.x * slope_11.x - cossin_phi.y * slope_11.y);
88         const float slope_y = alpha.y*(cossin_phi.y * slope_11.x + cossin_phi.x * slope_11.y);
89
90         kernel_assert(isfinite(slope_x));
91         return normalize(make_float3(-slope_x, -slope_y, 1.0f));
92 }
93
94 /* === Phase functions: Glossy, Diffuse and Glass === */
95
96 /* Phase function for reflective materials, either without a fresnel term (for compatibility) or with the conductive fresnel term. */
97 ccl_device_inline float3 mf_sample_phase_glossy(const float3 wi, float3 *n, float3 *k, float3 *weight, const float3 wm)
98 {
99         if(n && k)
100                 *weight *= fresnel_conductor(dot(wi, wm), *n, *k);
101
102         return -wi + 2.0f * wm * dot(wi, wm);
103 }
104
105 ccl_device_inline float3 mf_eval_phase_glossy(const float3 w, const float lambda, const float3 wo, const float2 alpha, float3 *n, float3 *k)
106 {
107         if(w.z > 0.9999f)
108                 return make_float3(0.0f, 0.0f, 0.0f);
109
110         const float3 wh = normalize(wo - w);
111         if(wh.z < 0.0f)
112                 return make_float3(0.0f, 0.0f, 0.0f);
113
114         float pArea = (w.z < -0.9999f)? 1.0f: lambda*w.z;
115
116         const float dotW_WH = dot(-w, wh);
117         if(dotW_WH < 0.0f)
118                 return make_float3(0.0f, 0.0f, 0.0f);
119
120         float phase = max(0.0f, dotW_WH) * 0.25f / max(pArea * dotW_WH, 1e-7f);
121         if(alpha.x == alpha.y)
122                 phase *= D_ggx(wh, alpha.x);
123         else
124                 phase *= D_ggx_aniso(wh, alpha);
125
126         if(n && k) {
127                 /* Apply conductive fresnel term. */
128                 return phase * fresnel_conductor(dotW_WH, *n, *k);
129         }
130
131         return make_float3(phase, phase, phase);
132 }
133
134 /* Phase function for rough lambertian diffuse surfaces. */
135 ccl_device_inline float3 mf_sample_phase_diffuse(const float3 wm, const float randu, const float randv)
136 {
137         float3 tm, bm;
138         make_orthonormals(wm, &tm, &bm);
139
140         float2 disk = concentric_sample_disk(randu, randv);
141         return disk.x*tm + disk.y*bm + safe_sqrtf(1.0f - disk.x*disk.x - disk.y*disk.y)*wm;
142 }
143
144 ccl_device_inline float3 mf_eval_phase_diffuse(const float3 w, const float3 wm)
145 {
146         const float v = max(0.0f, dot(w, wm)) * M_1_PI_F;
147         return make_float3(v, v, v);
148 }
149
150 /* Phase function for dielectric transmissive materials, including both reflection and refraction according to the dielectric fresnel term. */
151 ccl_device_inline float3 mf_sample_phase_glass(const float3 wi, const float eta, const float3 wm, const float randV, bool *outside)
152 {
153         float cosI = dot(wi, wm);
154         float f = fresnel_dielectric_cos(cosI, eta);
155         if(randV < f) {
156                 *outside = true;
157                 return -wi + 2.0f * wm * cosI;
158         }
159         *outside = false;
160         float inv_eta = 1.0f/eta;
161         float cosT = -safe_sqrtf(1.0f - (1.0f - cosI*cosI) * inv_eta*inv_eta);
162         return normalize(wm*(cosI*inv_eta + cosT) - wi*inv_eta);
163 }
164
165 ccl_device_inline float3 mf_eval_phase_glass(const float3 w, const float lambda, const float3 wo, const bool wo_outside, const float2 alpha, const float eta)
166 {
167         if(w.z > 0.9999f)
168                 return make_float3(0.0f, 0.0f, 0.0f);
169
170         float pArea = (w.z < -0.9999f)? 1.0f: lambda*w.z;
171         float v;
172         if(wo_outside) {
173                 const float3 wh = normalize(wo - w);
174                 if(wh.z < 0.0f)
175                         return make_float3(0.0f, 0.0f, 0.0f);
176
177                 const float dotW_WH = dot(-w, wh);
178                 v = fresnel_dielectric_cos(dotW_WH, eta) * max(0.0f, dotW_WH) * D_ggx(wh, alpha.x) * 0.25f / (pArea * dotW_WH);
179         }
180         else {
181                 float3 wh = normalize(wo*eta - w);
182                 if(wh.z < 0.0f)
183                         wh = -wh;
184                 const float dotW_WH = dot(-w, wh), dotWO_WH = dot(wo, wh);
185                 if(dotW_WH < 0.0f)
186                         return make_float3(0.0f, 0.0f, 0.0f);
187
188                 float temp = dotW_WH + eta*dotWO_WH;
189                 v = (1.0f - fresnel_dielectric_cos(dotW_WH, eta)) * max(0.0f, dotW_WH) * max(0.0f, -dotWO_WH) * D_ggx(wh, alpha.x) / (pArea * temp * temp);
190         }
191
192         return make_float3(v, v, v);
193 }
194
195 /* === Utility functions for the random walks === */
196
197 /* Smith Lambda function for GGX (based on page 12 of the supplemental implementation). */
198 ccl_device_inline float mf_lambda(const float3 w, const float2 alpha)
199 {
200         if(w.z > 0.9999f)
201                 return 0.0f;
202         else if(w.z < -0.9999f)
203                 return -0.9999f;
204
205         const float inv_wz2 = 1.0f / max(w.z*w.z, 1e-7f);
206         const float2 wa = make_float2(w.x, w.y)*alpha;
207         float v = sqrtf(1.0f + dot(wa, wa) * inv_wz2);
208         if(w.z <= 0.0f)
209                 v = -v;
210
211         return 0.5f*(v - 1.0f);
212 }
213
214 /* Height distribution CDF (based on page 4 of the supplemental implementation). */
215 ccl_device_inline float mf_invC1(const float h)
216 {
217         return 2.0f * saturate(h) - 1.0f;
218 }
219
220 ccl_device_inline float mf_C1(const float h)
221 {
222         return saturate(0.5f * (h + 1.0f));
223 }
224
225 /* Masking function (based on page 16 of the supplemental implementation). */
226 ccl_device_inline float mf_G1(const float3 w, const float C1, const float lambda)
227 {
228         if(w.z > 0.9999f)
229                 return 1.0f;
230         if(w.z < 1e-5f)
231                 return 0.0f;
232         return powf(C1, lambda);
233 }
234
235 /* Sampling from the visible height distribution (based on page 17 of the supplemental implementation). */
236 ccl_device_inline bool mf_sample_height(const float3 w, float *h, float *C1, float *G1, float *lambda, const float U)
237 {
238         if(w.z > 0.9999f)
239                 return false;
240         if(w.z < -0.9999f) {
241                 *C1 *= U;
242                 *h = mf_invC1(*C1);
243                 *G1 = mf_G1(w, *C1, *lambda);
244         }
245         else if(fabsf(w.z) >= 0.0001f) {
246                 if(U > 1.0f - *G1)
247                         return false;
248                 if(*lambda >= 0.0f) {
249                         *C1 = 1.0f;
250                 }
251                 else {
252                         *C1 *= powf(1.0f-U, -1.0f / *lambda);
253                 }
254                 *h = mf_invC1(*C1);
255                 *G1 = mf_G1(w, *C1, *lambda);
256         }
257         return true;
258 }
259
260 /* === PDF approximations for the different phase functions. ===
261  * As explained in bsdf_microfacet_multi_impl.h, using approximations with MIS still produces an unbiased result. */
262
263 /* Approximation for the albedo of the single-scattering GGX distribution,
264  * the missing energy is then approximated as a diffuse reflection for the PDF. */
265 ccl_device_inline float mf_ggx_albedo(float r)
266 {
267         float albedo = 0.806495f*expf(-1.98712f*r*r) + 0.199531f;
268         albedo -= ((((((1.76741f*r - 8.43891f)*r + 15.784f)*r - 14.398f)*r + 6.45221f)*r - 1.19722f)*r + 0.027803f)*r + 0.00568739f;
269         return saturate(albedo);
270 }
271
272 ccl_device_inline float mf_ggx_pdf(const float3 wi, const float3 wo, const float alpha)
273 {
274         float D = D_ggx(normalize(wi+wo), alpha);
275         float lambda = mf_lambda(wi, make_float2(alpha, alpha));
276         float albedo = mf_ggx_albedo(alpha);
277         return 0.25f * D / max((1.0f + lambda) * wi.z, 1e-7f) + (1.0f - albedo) * wo.z;
278 }
279
280 ccl_device_inline float mf_ggx_aniso_pdf(const float3 wi, const float3 wo, const float2 alpha)
281 {
282         return 0.25f * D_ggx_aniso(normalize(wi+wo), alpha) / ((1.0f + mf_lambda(wi, alpha)) * wi.z) + (1.0f - mf_ggx_albedo(sqrtf(alpha.x*alpha.y))) * wo.z;
283 }
284
285 ccl_device_inline float mf_diffuse_pdf(const float3 wo)
286 {
287         return M_1_PI_F * wo.z;
288 }
289
290 ccl_device_inline float mf_glass_pdf(const float3 wi, const float3 wo, const float alpha, const float eta)
291 {
292         float3 wh;
293         float fresnel;
294         if(wi.z*wo.z > 0.0f) {
295                 wh = normalize(wi + wo);
296                 fresnel = fresnel_dielectric_cos(dot(wi, wh), eta);
297         }
298         else {
299                 wh = normalize(wi + wo*eta);
300                 fresnel = 1.0f - fresnel_dielectric_cos(dot(wi, wh), eta);
301         }
302         if(wh.z < 0.0f)
303                 wh = -wh;
304         float3 r_wi = (wi.z < 0.0f)? -wi: wi;
305         return fresnel * max(0.0f, dot(r_wi, wh)) * D_ggx(wh, alpha) / ((1.0f + mf_lambda(r_wi, make_float2(alpha, alpha))) * r_wi.z) + fabsf(wo.z);
306 }
307
308 /* === Actual random walk implementations, one version of mf_eval and mf_sample per phase function. === */
309
310 #define MF_NAME_JOIN(x,y) x ## _ ## y
311 #define MF_NAME_EVAL(x,y) MF_NAME_JOIN(x,y)
312 #define MF_FUNCTION_FULL_NAME(prefix) MF_NAME_EVAL(prefix, MF_PHASE_FUNCTION)
313
314 #define MF_PHASE_FUNCTION glass
315 #define MF_MULTI_GLASS
316 #include "bsdf_microfacet_multi_impl.h"
317
318 /* The diffuse phase function is not implemented as a node yet. */
319 #if 0
320 #define MF_PHASE_FUNCTION diffuse
321 #define MF_MULTI_DIFFUSE
322 #include "bsdf_microfacet_multi_impl.h"
323 #endif
324
325 #define MF_PHASE_FUNCTION glossy
326 #define MF_MULTI_GLOSSY
327 #include "bsdf_microfacet_multi_impl.h"
328
329 ccl_device void bsdf_microfacet_multi_ggx_blur(ShaderClosure *sc, float roughness)
330 {
331         MicrofacetBsdf *bsdf = (MicrofacetBsdf*)sc;
332
333         bsdf->alpha_x = fmaxf(roughness, bsdf->alpha_x);
334         bsdf->alpha_y = fmaxf(roughness, bsdf->alpha_y);
335 }
336
337 /* === Closure implementations === */
338
339 /* Multiscattering GGX Glossy closure */
340
341 ccl_device int bsdf_microfacet_multi_ggx_common_setup(MicrofacetBsdf *bsdf)
342 {
343         bsdf->alpha_x = clamp(bsdf->alpha_x, 1e-4f, 1.0f);
344         bsdf->alpha_y = clamp(bsdf->alpha_y, 1e-4f, 1.0f);
345         bsdf->extra->color.x = saturate(bsdf->extra->color.x);
346         bsdf->extra->color.y = saturate(bsdf->extra->color.y);
347         bsdf->extra->color.z = saturate(bsdf->extra->color.z);
348
349         bsdf->type = CLOSURE_BSDF_MICROFACET_MULTI_GGX_ID;
350
351         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_NEEDS_LCG;
352 }
353
354 ccl_device int bsdf_microfacet_multi_ggx_aniso_setup(MicrofacetBsdf *bsdf)
355 {
356         if(is_zero(bsdf->T))
357                 bsdf->T = make_float3(1.0f, 0.0f, 0.0f);
358
359         return bsdf_microfacet_multi_ggx_common_setup(bsdf);
360 }
361
362 ccl_device int bsdf_microfacet_multi_ggx_setup(MicrofacetBsdf *bsdf)
363 {
364         bsdf->alpha_y = bsdf->alpha_x;
365
366         return bsdf_microfacet_multi_ggx_common_setup(bsdf);
367 }
368
369 ccl_device float3 bsdf_microfacet_multi_ggx_eval_transmit(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf, ccl_addr_space uint *lcg_state) {
370         *pdf = 0.0f;
371         return make_float3(0.0f, 0.0f, 0.0f);
372 }
373
374 ccl_device float3 bsdf_microfacet_multi_ggx_eval_reflect(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf, ccl_addr_space uint *lcg_state) {
375         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
376
377         if(bsdf->alpha_x*bsdf->alpha_y < 1e-7f) {
378                 return make_float3(0.0f, 0.0f, 0.0f);
379         }
380
381         bool is_aniso = (bsdf->alpha_x != bsdf->alpha_y);
382         float3 X, Y, Z;
383         Z = bsdf->N;
384         if(is_aniso)
385                 make_orthonormals_tangent(Z, bsdf->T, &X, &Y);
386         else
387                 make_orthonormals(Z, &X, &Y);
388
389         float3 localI = make_float3(dot(I, X), dot(I, Y), dot(I, Z));
390         float3 localO = make_float3(dot(omega_in, X), dot(omega_in, Y), dot(omega_in, Z));
391
392         if(is_aniso)
393                 *pdf = mf_ggx_aniso_pdf(localI, localO, make_float2(bsdf->alpha_x, bsdf->alpha_y));
394         else
395                 *pdf = mf_ggx_pdf(localI, localO, bsdf->alpha_x);
396         return mf_eval_glossy(localI, localO, true, bsdf->extra->color, bsdf->alpha_x, bsdf->alpha_y, lcg_state, NULL, NULL);
397 }
398
399 ccl_device int bsdf_microfacet_multi_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, ccl_addr_space uint *lcg_state)
400 {
401         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
402
403         float3 X, Y, Z;
404         Z = bsdf->N;
405
406         if(bsdf->alpha_x*bsdf->alpha_y < 1e-7f) {
407                 *omega_in = 2*dot(Z, I)*Z - I;
408                 *pdf = 1e6f;
409                 *eval = make_float3(1e6f, 1e6f, 1e6f);
410                 return LABEL_REFLECT|LABEL_SINGULAR;
411         }
412
413         bool is_aniso = (bsdf->alpha_x != bsdf->alpha_y);
414         if(is_aniso)
415                 make_orthonormals_tangent(Z, bsdf->T, &X, &Y);
416         else
417                 make_orthonormals(Z, &X, &Y);
418
419         float3 localI = make_float3(dot(I, X), dot(I, Y), dot(I, Z));
420         float3 localO;
421
422         *eval = mf_sample_glossy(localI, &localO, bsdf->extra->color, bsdf->alpha_x, bsdf->alpha_y, lcg_state, NULL, NULL);
423         if(is_aniso)
424                 *pdf = mf_ggx_aniso_pdf(localI, localO, make_float2(bsdf->alpha_x, bsdf->alpha_y));
425         else
426                 *pdf = mf_ggx_pdf(localI, localO, bsdf->alpha_x);
427         *eval *= *pdf;
428
429         *omega_in = X*localO.x + Y*localO.y + Z*localO.z;
430 #ifdef __RAY_DIFFERENTIALS__
431         *domega_in_dx = (2 * dot(Z, dIdx)) * Z - dIdx;
432         *domega_in_dy = (2 * dot(Z, dIdy)) * Z - dIdy;
433 #endif
434         return LABEL_REFLECT|LABEL_GLOSSY;
435 }
436
437 /* Multiscattering GGX Glass closure */
438
439 ccl_device int bsdf_microfacet_multi_ggx_glass_setup(MicrofacetBsdf *bsdf)
440 {
441         bsdf->alpha_x = clamp(bsdf->alpha_x, 1e-4f, 1.0f);
442         bsdf->alpha_y = bsdf->alpha_x;
443         bsdf->ior = max(0.0f, bsdf->ior);
444         bsdf->extra->color.x = saturate(bsdf->extra->color.x);
445         bsdf->extra->color.y = saturate(bsdf->extra->color.y);
446         bsdf->extra->color.z = saturate(bsdf->extra->color.z);
447
448         bsdf->type = CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_ID;
449
450         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_NEEDS_LCG;
451 }
452
453 ccl_device float3 bsdf_microfacet_multi_ggx_glass_eval_transmit(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf, ccl_addr_space uint *lcg_state) {
454         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
455
456         if(bsdf->alpha_x*bsdf->alpha_y < 1e-7f) {
457                 return make_float3(0.0f, 0.0f, 0.0f);
458         }
459
460         float3 X, Y, Z;
461         Z = bsdf->N;
462         make_orthonormals(Z, &X, &Y);
463
464         float3 localI = make_float3(dot(I, X), dot(I, Y), dot(I, Z));
465         float3 localO = make_float3(dot(omega_in, X), dot(omega_in, Y), dot(omega_in, Z));
466
467         *pdf = mf_glass_pdf(localI, localO, bsdf->alpha_x, bsdf->ior);
468         return mf_eval_glass(localI, localO, false, bsdf->extra->color, bsdf->alpha_x, bsdf->alpha_y, lcg_state, bsdf->ior);
469 }
470
471 ccl_device float3 bsdf_microfacet_multi_ggx_glass_eval_reflect(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf, ccl_addr_space uint *lcg_state) {
472         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
473
474         if(bsdf->alpha_x*bsdf->alpha_y < 1e-7f) {
475                 return make_float3(0.0f, 0.0f, 0.0f);
476         }
477
478         float3 X, Y, Z;
479         Z = bsdf->N;
480         make_orthonormals(Z, &X, &Y);
481
482         float3 localI = make_float3(dot(I, X), dot(I, Y), dot(I, Z));
483         float3 localO = make_float3(dot(omega_in, X), dot(omega_in, Y), dot(omega_in, Z));
484
485         *pdf = mf_glass_pdf(localI, localO, bsdf->alpha_x, bsdf->ior);
486         return mf_eval_glass(localI, localO, true, bsdf->extra->color, bsdf->alpha_x, bsdf->alpha_y, lcg_state, bsdf->ior);
487 }
488
489 ccl_device int bsdf_microfacet_multi_ggx_glass_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, ccl_addr_space uint *lcg_state)
490 {
491         const MicrofacetBsdf *bsdf = (const MicrofacetBsdf*)sc;
492
493         float3 X, Y, Z;
494         Z = bsdf->N;
495
496         if(bsdf->alpha_x*bsdf->alpha_y < 1e-7f) {
497                 float3 R, T;
498 #ifdef __RAY_DIFFERENTIALS__
499                 float3 dRdx, dRdy, dTdx, dTdy;
500 #endif
501                 bool inside;
502                 float fresnel = fresnel_dielectric(bsdf->ior, Z, I, &R, &T,
503 #ifdef __RAY_DIFFERENTIALS__
504                                 dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,
505 #endif
506                                 &inside);
507
508                 *pdf = 1e6f;
509                 *eval = make_float3(1e6f, 1e6f, 1e6f);
510                 if(randu < fresnel) {
511                         *omega_in = R;
512 #ifdef __RAY_DIFFERENTIALS__
513                         *domega_in_dx = dRdx;
514                         *domega_in_dy = dRdy;
515 #endif
516                         return LABEL_REFLECT|LABEL_SINGULAR;
517                 }
518                 else {
519                         *omega_in = T;
520 #ifdef __RAY_DIFFERENTIALS__
521                         *domega_in_dx = dTdx;
522                         *domega_in_dy = dTdy;
523 #endif
524                         return LABEL_TRANSMIT|LABEL_SINGULAR;
525                 }
526         }
527
528         make_orthonormals(Z, &X, &Y);
529
530         float3 localI = make_float3(dot(I, X), dot(I, Y), dot(I, Z));
531         float3 localO;
532
533         *eval = mf_sample_glass(localI, &localO, bsdf->extra->color, bsdf->alpha_x, bsdf->alpha_y, lcg_state, bsdf->ior);
534         *pdf = mf_glass_pdf(localI, localO, bsdf->alpha_x, bsdf->ior);
535         *eval *= *pdf;
536
537         *omega_in = X*localO.x + Y*localO.y + Z*localO.z;
538         if(localO.z*localI.z > 0.0f) {
539 #ifdef __RAY_DIFFERENTIALS__
540                 *domega_in_dx = (2 * dot(Z, dIdx)) * Z - dIdx;
541                 *domega_in_dy = (2 * dot(Z, dIdy)) * Z - dIdy;
542 #endif
543                 return LABEL_REFLECT|LABEL_GLOSSY;
544         }
545         else {
546 #ifdef __RAY_DIFFERENTIALS__
547                 float cosI = dot(Z, I);
548                 float dnp = max(sqrtf(1.0f - (bsdf->ior * bsdf->ior * (1.0f - cosI*cosI))), 1e-7f);
549                 *domega_in_dx = -(bsdf->ior * dIdx) + ((bsdf->ior - bsdf->ior * bsdf->ior * cosI / dnp) * dot(dIdx, Z)) * Z;
550                 *domega_in_dy = -(bsdf->ior * dIdy) + ((bsdf->ior - bsdf->ior * bsdf->ior * cosI / dnp) * dot(dIdy, Z)) * Z;
551 #endif
552
553                 return LABEL_TRANSMIT|LABEL_GLOSSY;
554         }
555 }
556
557 CCL_NAMESPACE_END