1ec35e444fe42bf88ca153b13158cb8e71ee8d44
[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 /* GGX */
39
40 ccl_device int bsdf_microfacet_ggx_setup(ShaderClosure *sc)
41 {
42         sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ag */
43         
44         sc->type = CLOSURE_BSDF_MICROFACET_GGX_ID;
45
46         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
47 }
48
49 ccl_device int bsdf_microfacet_ggx_refraction_setup(ShaderClosure *sc)
50 {
51         sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ag */
52
53         sc->type = CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
54
55         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
56 }
57
58 ccl_device void bsdf_microfacet_ggx_blur(ShaderClosure *sc, float roughness)
59 {
60         sc->data0 = fmaxf(roughness, sc->data0); /* m_ag */
61 }
62
63 ccl_device float3 bsdf_microfacet_ggx_eval_reflect(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
64 {
65         float m_ag = max(sc->data0, 1e-4f);
66         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
67         float3 N = sc->N;
68
69         if(m_refractive || m_ag <= 1e-4f)
70                 return make_float3 (0, 0, 0);
71         float cosNO = dot(N, I);
72         float cosNI = dot(N, omega_in);
73         if(cosNI > 0 && cosNO > 0) {
74                 // get half vector
75                 float3 Hr = normalize(omega_in + I);
76                 // eq. 20: (F*G*D)/(4*in*on)
77                 // eq. 33: first we calculate D(m) with m=Hr:
78                 float alpha2 = m_ag * m_ag;
79                 float cosThetaM = dot(N, Hr);
80                 float cosThetaM2 = cosThetaM * cosThetaM;
81                 float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
82                 float cosThetaM4 = cosThetaM2 * cosThetaM2;
83                 float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
84                 // eq. 34: now calculate G1(i,m) and G1(o,m)
85                 float G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
86                 float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
87                 float G = G1o * G1i;
88                 float out = (G * D) * 0.25f / cosNO;
89                 // eq. 24
90                 float pm = D * cosThetaM;
91                 // convert into pdf of the sampled direction
92                 // eq. 38 - but see also:
93                 // eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
94                 *pdf = pm * 0.25f / dot(Hr, I);
95                 return make_float3 (out, out, out);
96         }
97         return make_float3 (0, 0, 0);
98 }
99
100 ccl_device float3 bsdf_microfacet_ggx_eval_transmit(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
101 {
102         float m_ag = max(sc->data0, 1e-4f);
103         float m_eta = sc->data1;
104         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
105         float3 N = sc->N;
106
107         if(!m_refractive || m_ag <= 1e-4f)
108                 return make_float3 (0, 0, 0);
109         float cosNO = dot(N, I);
110         float cosNI = dot(N, omega_in);
111         if(cosNO <= 0 || cosNI >= 0)
112                 return make_float3 (0, 0, 0); // vectors on same side -- not possible
113         // compute half-vector of the refraction (eq. 16)
114         float3 ht = -(m_eta * omega_in + I);
115         float3 Ht = normalize(ht);
116         float cosHO = dot(Ht, I);
117
118         float cosHI = dot(Ht, omega_in);
119         // eq. 33: first we calculate D(m) with m=Ht:
120         float alpha2 = m_ag * m_ag;
121         float cosThetaM = dot(N, Ht);
122         float cosThetaM2 = cosThetaM * cosThetaM;
123         float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
124         float cosThetaM4 = cosThetaM2 * cosThetaM2;
125         float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
126         // eq. 34: now calculate G1(i,m) and G1(o,m)
127         float G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
128         float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
129         float G = G1o * G1i;
130         // probability
131         float invHt2 = 1 / dot(ht, ht);
132         *pdf = D * fabsf(cosThetaM) * (fabsf(cosHI) * (m_eta * m_eta)) * invHt2;
133         float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D) * invHt2) / cosNO;
134         return make_float3 (out, out, out);
135 }
136
137 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)
138 {
139         float m_ag = sc->data0;
140         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID;
141         float3 N = sc->N;
142
143         float cosNO = dot(N, I);
144         if(cosNO > 0) {
145                 float3 X, Y, Z = N;
146                 make_orthonormals(Z, &X, &Y);
147                 // generate a random microfacet normal m
148                 // eq. 35,36:
149                 // we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)
150                 //tttt  and sin(atan(x)) == x/sqrt(1+x^2)
151                 float alpha2 = m_ag * m_ag;
152                 float tanThetaM2 = alpha2 * randu / (1 - randu);
153                 float cosThetaM  = 1 / safe_sqrtf(1 + tanThetaM2);
154                 float sinThetaM  = cosThetaM * safe_sqrtf(tanThetaM2);
155                 float phiM = M_2PI_F * randv;
156                 float3 m = (cosf(phiM) * sinThetaM) * X +
157                            (sinf(phiM) * sinThetaM) * Y +
158                            (             cosThetaM) * Z;
159                 if(!m_refractive) {
160                         float cosMO = dot(m, I);
161                         if(cosMO > 0) {
162                                 // eq. 39 - compute actual reflected direction
163                                 *omega_in = 2 * cosMO * m - I;
164                                 if(dot(Ng, *omega_in) > 0) {
165                                         if (m_ag <= 1e-4f) {
166                                                 // some high number for MIS
167                                                 *pdf = 1e6f;
168                                                 *eval = make_float3(1e6f, 1e6f, 1e6f);
169                                         }
170                                         else {
171                                                 // microfacet normal is visible to this ray
172                                                 // eq. 33
173                                                 float cosThetaM2 = cosThetaM * cosThetaM;
174                                                 float cosThetaM4 = cosThetaM2 * cosThetaM2;
175                                                 float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
176                                                 // eq. 24
177                                                 float pm = D * cosThetaM;
178                                                 // convert into pdf of the sampled direction
179                                                 // eq. 38 - but see also:
180                                                 // eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
181                                                 *pdf = pm * 0.25f / cosMO;
182                                                 // eval BRDF*cosNI
183                                                 float cosNI = dot(N, *omega_in);
184                                                 // eq. 34: now calculate G1(i,m) and G1(o,m)
185                                                 float G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
186                                                 float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
187                                                 float G = G1o * G1i;
188                                                 // eq. 20: (F*G*D)/(4*in*on)
189                                                 float out = (G * D) * 0.25f / cosNO;
190                                                 *eval = make_float3(out, out, out);
191                                         }
192
193 #ifdef __RAY_DIFFERENTIALS__
194                                         *domega_in_dx = (2 * dot(m, dIdx)) * m - dIdx;
195                                         *domega_in_dy = (2 * dot(m, dIdy)) * m - dIdy;
196 #endif
197                                 }
198                         }
199                 }
200                 else {
201                         // CAUTION: the i and o variables are inverted relative to the paper
202                         // eq. 39 - compute actual refractive direction
203                         float3 R, T;
204 #ifdef __RAY_DIFFERENTIALS__
205                         float3 dRdx, dRdy, dTdx, dTdy;
206 #endif
207                         float m_eta = sc->data1;
208                         bool inside;
209                         fresnel_dielectric(m_eta, m, I, &R, &T,
210 #ifdef __RAY_DIFFERENTIALS__
211                                 dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,
212 #endif
213                                 &inside);
214                         
215                         if(!inside) {
216                                 *omega_in = T;
217 #ifdef __RAY_DIFFERENTIALS__
218                                 *domega_in_dx = dTdx;
219                                 *domega_in_dy = dTdy;
220 #endif
221
222                                 if (m_ag <= 1e-4f || fabsf(m_eta - 1.0f) < 1e-4f) {
223                                         // some high number for MIS
224                                         *pdf = 1e6f;
225                                         *eval = make_float3(1e6f, 1e6f, 1e6f);
226                                 }
227                                 else {
228                                         // eq. 33
229                                         float cosThetaM2 = cosThetaM * cosThetaM;
230                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
231                                         float D = alpha2 / (M_PI_F * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
232                                         // eq. 24
233                                         float pm = D * cosThetaM;
234                                         // eval BRDF*cosNI
235                                         float cosNI = dot(N, *omega_in);
236                                         // eq. 34: now calculate G1(i,m) and G1(o,m)
237                                         float G1o = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
238                                         float G1i = 2 / (1 + safe_sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
239                                         float G = G1o * G1i;
240                                         // eq. 21
241                                         float cosHI = dot(m, *omega_in);
242                                         float cosHO = dot(m, I);
243                                         float Ht2 = m_eta * cosHI + cosHO;
244                                         Ht2 *= Ht2;
245                                         float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D)) / (cosNO * Ht2);
246                                         // eq. 38 and eq. 17
247                                         *pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2;
248                                         *eval = make_float3(out, out, out);
249                                 }
250                         }
251                 }
252         }
253         return (m_refractive) ? LABEL_TRANSMIT|LABEL_GLOSSY : LABEL_REFLECT|LABEL_GLOSSY;
254 }
255
256 /* BECKMANN */
257
258 ccl_device int bsdf_microfacet_beckmann_setup(ShaderClosure *sc)
259 {
260         sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ab */
261
262         sc->type = CLOSURE_BSDF_MICROFACET_BECKMANN_ID;
263         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
264 }
265
266 ccl_device int bsdf_microfacet_beckmann_refraction_setup(ShaderClosure *sc)
267 {
268         sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ab */
269
270         sc->type = CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
271         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_GLOSSY;
272 }
273
274 ccl_device void bsdf_microfacet_beckmann_blur(ShaderClosure *sc, float roughness)
275 {
276         sc->data0 = fmaxf(roughness, sc->data0); /* m_ab */
277 }
278
279 ccl_device float3 bsdf_microfacet_beckmann_eval_reflect(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
280 {
281         float m_ab = max(sc->data0, 1e-4f);
282         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
283         float3 N = sc->N;
284
285         if(m_refractive || m_ab <= 1e-4f)
286                 return make_float3 (0, 0, 0);
287         float cosNO = dot(N, I);
288         float cosNI = dot(N, omega_in);
289         if(cosNO > 0 && cosNI > 0) {
290            // get half vector
291            float3 Hr = normalize(omega_in + I);
292            // eq. 20: (F*G*D)/(4*in*on)
293            // eq. 25: first we calculate D(m) with m=Hr:
294            float alpha2 = m_ab * m_ab;
295            float cosThetaM = dot(N, Hr);
296            float cosThetaM2 = cosThetaM * cosThetaM;
297            float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
298            float cosThetaM4 = cosThetaM2 * cosThetaM2;
299            float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
300            // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
301            float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
302            float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
303            float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
304            float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
305            float G = G1o * G1i;
306            float out = (G * D) * 0.25f / cosNO;
307            // eq. 24
308            float pm = D * cosThetaM;
309            // convert into pdf of the sampled direction
310            // eq. 38 - but see also:
311            // eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
312            *pdf = pm * 0.25f / dot(Hr, I);
313            return make_float3 (out, out, out);
314         }
315         return make_float3 (0, 0, 0);
316 }
317
318 ccl_device float3 bsdf_microfacet_beckmann_eval_transmit(const ShaderClosure *sc, const float3 I, const float3 omega_in, float *pdf)
319 {
320         float m_ab = max(sc->data0, 1e-4f);
321         float m_eta = sc->data1;
322         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
323         float3 N = sc->N;
324
325         if(!m_refractive || m_ab <= 1e-4f)
326                 return make_float3 (0, 0, 0);
327         float cosNO = dot(N, I);
328         float cosNI = dot(N, omega_in);
329         if(cosNO <= 0 || cosNI >= 0)
330                 return make_float3 (0, 0, 0);
331         // compute half-vector of the refraction (eq. 16)
332         float3 ht = -(m_eta * omega_in + I);
333         float3 Ht = normalize(ht);
334         float cosHO = dot(Ht, I);
335
336         float cosHI = dot(Ht, omega_in);
337         // eq. 33: first we calculate D(m) with m=Ht:
338         float alpha2 = m_ab * m_ab;
339         float cosThetaM = min(dot(N, Ht), 1.0f);
340         float cosThetaM2 = cosThetaM * cosThetaM;
341         float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
342         float cosThetaM4 = cosThetaM2 * cosThetaM2;
343         float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
344         // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
345         float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
346         float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
347         float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
348         float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
349         float G = G1o * G1i;
350         // probability
351         float invHt2 = 1 / dot(ht, ht);
352         *pdf = D * fabsf(cosThetaM) * (fabsf(cosHI) * (m_eta * m_eta)) * invHt2;
353         float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D) * invHt2) / cosNO;
354         return make_float3 (out, out, out);
355 }
356
357 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)
358 {
359         float m_ab = sc->data0;
360         int m_refractive = sc->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID;
361         float3 N = sc->N;
362
363         float cosNO = dot(N, I);
364         if(cosNO > 0) {
365                 float3 X, Y, Z = N;
366                 make_orthonormals(Z, &X, &Y);
367                 // generate a random microfacet normal m
368                 // eq. 35,36:
369                 // we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)
370                 //tttt  and sin(atan(x)) == x/sqrt(1+x^2)
371                 float alpha2 = m_ab * m_ab;
372                 float tanThetaM, cosThetaM;
373
374                 if(alpha2 == 0.0f) {
375                         tanThetaM = 0.0f;
376                         cosThetaM = 1.0f;
377                 }
378                 else {
379                         tanThetaM = safe_sqrtf(-alpha2 * logf(1 - randu));
380                         cosThetaM = 1 / safe_sqrtf(1 + tanThetaM * tanThetaM);
381                 }
382
383                 float sinThetaM = cosThetaM * tanThetaM;
384                 float phiM = M_2PI_F * randv;
385                 float3 m = (cosf(phiM) * sinThetaM) * X +
386                            (sinf(phiM) * sinThetaM) * Y +
387                            (             cosThetaM) * Z;
388
389                 if(!m_refractive) {
390                         float cosMO = dot(m, I);
391                         if(cosMO > 0) {
392                                 // eq. 39 - compute actual reflected direction
393                                 *omega_in = 2 * cosMO * m - I;
394                                 if(dot(Ng, *omega_in) > 0) {
395                                         if (m_ab <= 1e-4f) {
396                                                 // some high number for MIS
397                                                 *pdf = 1e6f;
398                                                 *eval = make_float3(1e6f, 1e6f, 1e6f);
399                                         }
400                                         else {
401                                                 // microfacet normal is visible to this ray
402                                                 // eq. 25
403                                                 float cosThetaM2 = cosThetaM * cosThetaM;
404                                                 float tanThetaM2 = tanThetaM * tanThetaM;
405                                                 float cosThetaM4 = cosThetaM2 * cosThetaM2;
406                                                 float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
407                                                 // eq. 24
408                                                 float pm = D * cosThetaM;
409                                                 // convert into pdf of the sampled direction
410                                                 // eq. 38 - but see also:
411                                                 // eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
412                                                 *pdf = pm * 0.25f / cosMO;
413                                                 // Eval BRDF*cosNI
414                                                 float cosNI = dot(N, *omega_in);
415                                                 // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
416                                                 float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
417                                                 float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
418                                                 float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
419                                                 float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
420                                                 float G = G1o * G1i;
421                                                 // eq. 20: (F*G*D)/(4*in*on)
422                                                 float out = (G * D) * 0.25f / cosNO;
423                                                 *eval = make_float3(out, out, out);
424                                         }
425 #ifdef __RAY_DIFFERENTIALS__
426                                         *domega_in_dx = (2 * dot(m, dIdx)) * m - dIdx;
427                                         *domega_in_dy = (2 * dot(m, dIdy)) * m - dIdy;
428 #endif
429                                 }
430                         }
431                 }
432                 else {
433                         // CAUTION: the i and o variables are inverted relative to the paper
434                         // eq. 39 - compute actual refractive direction
435                         float3 R, T;
436 #ifdef __RAY_DIFFERENTIALS__
437                         float3 dRdx, dRdy, dTdx, dTdy;
438 #endif
439                         float m_eta = sc->data1;
440                         bool inside;
441                         fresnel_dielectric(m_eta, m, I, &R, &T,
442 #ifdef __RAY_DIFFERENTIALS__
443                                 dIdx, dIdy, &dRdx, &dRdy, &dTdx, &dTdy,
444 #endif
445                                 &inside);
446
447                         if(!inside) {
448                                 *omega_in = T;
449 #ifdef __RAY_DIFFERENTIALS__
450                                 *domega_in_dx = dTdx;
451                                 *domega_in_dy = dTdy;
452 #endif
453                                 if (m_ab <= 1e-4f || fabsf(m_eta - 1.0f) < 1e-4f) {
454                                         // some high number for MIS
455                                         *pdf = 1e6f;
456                                         *eval = make_float3(1e6f, 1e6f, 1e6f);
457                                 }
458                                 else {
459                                         // eq. 33
460                                         float cosThetaM2 = cosThetaM * cosThetaM;
461                                         float tanThetaM2 = tanThetaM * tanThetaM;
462                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
463                                         float D = expf(-tanThetaM2 / alpha2) / (M_PI_F * alpha2 *  cosThetaM4);
464                                         // eq. 24
465                                         float pm = D * cosThetaM;
466                                         // eval BRDF*cosNI
467                                         float cosNI = dot(N, *omega_in);
468                                         // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
469                                         float ao = 1 / (m_ab * safe_sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
470                                         float ai = 1 / (m_ab * safe_sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
471                                         float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
472                                         float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
473                                         float G = G1o * G1i;
474                                         // eq. 21
475                                         float cosHI = dot(m, *omega_in);
476                                         float cosHO = dot(m, I);
477                                         float Ht2 = m_eta * cosHI + cosHO;
478                                         Ht2 *= Ht2;
479                                         float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D)) / (cosNO * Ht2);
480                                         // eq. 38 and eq. 17
481                                         *pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2;
482                                         *eval = make_float3(out, out, out);
483                                 }
484                         }
485                 }
486         }
487         return (m_refractive) ? LABEL_TRANSMIT|LABEL_GLOSSY : LABEL_REFLECT|LABEL_GLOSSY;
488 }
489
490 CCL_NAMESPACE_END
491
492 #endif /* __BSDF_MICROFACET_H__ */
493