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