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