Cycles: svn merge -r41225:41232 ^/trunk/blender
[blender.git] / intern / cycles / kernel / osl / bsdf_microfacet.cpp
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 #include <OpenImageIO/fmath.h>
34
35 #include <OSL/genclosure.h>
36
37 #include "osl_closures.h"
38
39 #include "util_math.h"
40
41 using namespace OSL;
42
43 CCL_NAMESPACE_BEGIN
44
45 // TODO: fresnel_dielectric is only used for derivatives, could be optimized
46
47 // TODO: refactor these two classes so they share everything by the microfacet
48 //       distribution terms
49
50 // microfacet model with GGX facet distribution
51 // see http://www.graphics.cornell.edu/~bjw/microfacetbsdf.pdf
52 template <int Refractive = 0>
53 class MicrofacetGGXClosure : public BSDFClosure {
54 public:
55     Vec3 m_N;
56     float m_ag;   // width parameter (roughness)
57     float m_eta;  // index of refraction (for fresnel term)
58     MicrofacetGGXClosure() : BSDFClosure(Labels::GLOSSY, Refractive ? Back : Front) { m_eta = 1.0f; }
59
60     void setup()
61         {
62                 m_ag = clamp(m_ag, 1e-5f, 1.0f);
63         }
64
65     bool mergeable (const ClosurePrimitive *other) const {
66         const MicrofacetGGXClosure *comp = (const MicrofacetGGXClosure *)other;
67         return m_N == comp->m_N && m_ag == comp->m_ag &&
68             m_eta == comp->m_eta && BSDFClosure::mergeable(other);
69     }
70
71     size_t memsize () const { return sizeof(*this); }
72
73     const char *name () const {
74         return Refractive ? "microfacet_ggx_refraction" : "microfacet_ggx";
75     }
76
77     void print_on (std::ostream &out) const {
78         out << name() << " (";
79         out << "(" << m_N[0] << ", " << m_N[1] << ", " << m_N[2] << "), ";
80         out << m_ag << ", ";
81         out << m_eta;
82         out << ")";
83     }
84
85     float albedo (const Vec3 &omega_out) const
86     {
87                 return 1.0f;
88     }
89
90     Color3 eval_reflect (const Vec3 &omega_out, const Vec3 &omega_in, float& pdf) const
91     {
92         if (Refractive == 1) return Color3 (0, 0, 0);
93         float cosNO = m_N.dot(omega_out);
94         float cosNI = m_N.dot(omega_in);
95         if (cosNI > 0 && cosNO > 0) {
96             // get half vector
97             Vec3 Hr = omega_in + omega_out;
98             Hr.normalize();
99             // eq. 20: (F*G*D)/(4*in*on)
100             // eq. 33: first we calculate D(m) with m=Hr:
101             float alpha2 = m_ag * m_ag;
102             float cosThetaM = m_N.dot(Hr);
103             float cosThetaM2 = cosThetaM * cosThetaM;
104             float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
105             float cosThetaM4 = cosThetaM2 * cosThetaM2;
106             float D = alpha2 / ((float) M_PI * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
107             // eq. 34: now calculate G1(i,m) and G1(o,m)
108             float G1o = 2 / (1 + sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
109             float G1i = 2 / (1 + sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
110             float G = G1o * G1i;
111             float out = (G * D) * 0.25f / cosNO;
112             // eq. 24
113             float pm = D * cosThetaM;
114             // convert into pdf of the sampled direction
115             // eq. 38 - but see also:
116             // eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
117             pdf = pm * 0.25f / Hr.dot(omega_out);
118             return Color3 (out, out, out);
119         }
120         return Color3 (0, 0, 0);
121     }
122
123     Color3 eval_transmit (const Vec3 &omega_out, const Vec3 &omega_in, float& pdf) const
124     {
125         if (Refractive == 0) return Color3 (0, 0, 0);
126         float cosNO = m_N.dot(omega_out);
127         float cosNI = m_N.dot(omega_in);
128         if (cosNO <= 0 || cosNI >= 0)
129             return Color3 (0, 0, 0); // vectors on same side -- not possible
130         // compute half-vector of the refraction (eq. 16)
131         Vec3 ht = -(m_eta * omega_in + omega_out);
132         Vec3 Ht = ht; Ht.normalize();
133         float cosHO = Ht.dot(omega_out);
134
135                 float cosHI = Ht.dot(omega_in);
136                 // eq. 33: first we calculate D(m) with m=Ht:
137                 float alpha2 = m_ag * m_ag;
138                 float cosThetaM = m_N.dot(Ht);
139                 float cosThetaM2 = cosThetaM * cosThetaM;
140                 float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
141                 float cosThetaM4 = cosThetaM2 * cosThetaM2;
142                 float D = alpha2 / ((float) M_PI * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
143                 // eq. 34: now calculate G1(i,m) and G1(o,m)
144                 float G1o = 2 / (1 + sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
145                 float G1i = 2 / (1 + sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
146                 float G = G1o * G1i;
147                 // probability
148                 float invHt2 = 1 / ht.dot(ht);
149                 pdf = D * fabsf(cosThetaM) * (fabsf(cosHI) * (m_eta * m_eta)) * invHt2;
150                 float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D) * invHt2) / cosNO;
151                 return Color3 (out, out, out);
152     }
153
154     ustring sample (const Vec3 &Ng,
155                  const Vec3 &omega_out, const Vec3 &domega_out_dx, const Vec3 &domega_out_dy,
156                  float randu, float randv,
157                  Vec3 &omega_in, Vec3 &domega_in_dx, Vec3 &domega_in_dy,
158                  float &pdf, Color3 &eval) const
159     {
160         float cosNO = m_N.dot(omega_out);
161         if (cosNO > 0) {
162             Vec3 X, Y, Z = m_N;
163             make_orthonormals(Z, X, Y);
164             // generate a random microfacet normal m
165             // eq. 35,36:
166             // we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)
167             //                  and sin(atan(x)) == x/sqrt(1+x^2)
168             float alpha2 = m_ag * m_ag;
169             float tanThetaM2 = alpha2 * randu / (1 - randu);
170             float cosThetaM  = 1 / sqrtf(1 + tanThetaM2);
171             float sinThetaM  = cosThetaM * sqrtf(tanThetaM2);
172             float phiM = 2 * float(M_PI) * randv;
173             Vec3 m = (cosf(phiM) * sinThetaM) * X +
174                      (sinf(phiM) * sinThetaM) * Y +
175                                    cosThetaM  * Z;
176             if (Refractive == 0) {
177                 float cosMO = m.dot(omega_out);
178                 if (cosMO > 0) {
179                     // eq. 39 - compute actual reflected direction
180                     omega_in = 2 * cosMO * m - omega_out;
181                     if (Ng.dot(omega_in) > 0) {
182                         // microfacet normal is visible to this ray
183                         // eq. 33
184                         float cosThetaM2 = cosThetaM * cosThetaM;
185                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
186                         float D = alpha2 / (float(M_PI) * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
187                         // eq. 24
188                         float pm = D * cosThetaM;
189                         // convert into pdf of the sampled direction
190                         // eq. 38 - but see also:
191                         // eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
192                         pdf = pm * 0.25f / cosMO;
193                         // eval BRDF*cosNI
194                         float cosNI = m_N.dot(omega_in);
195                         // eq. 34: now calculate G1(i,m) and G1(o,m)
196                         float G1o = 2 / (1 + sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
197                         float G1i = 2 / (1 + sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
198                         float G = G1o * G1i;
199                         // eq. 20: (F*G*D)/(4*in*on)
200                         float out = (G * D) * 0.25f / cosNO;
201                         eval.setValue(out, out, out);
202                         domega_in_dx = (2 * m.dot(domega_out_dx)) * m - domega_out_dx;
203                         domega_in_dy = (2 * m.dot(domega_out_dy)) * m - domega_out_dy;
204
205                                                 /* disabled for now - gives texture filtering problems */
206 #if 0
207                         // Since there is some blur to this reflection, make the
208                         // derivatives a bit bigger. In theory this varies with the
209                         // roughness but the exact relationship is complex and
210                         // requires more ops than are practical.
211                         domega_in_dx *= 10;
212                         domega_in_dy *= 10;
213 #endif
214                     }
215                 }
216             } else {
217                 // CAUTION: the i and o variables are inverted relative to the paper
218                 // eq. 39 - compute actual refractive direction
219                 Vec3 R, dRdx, dRdy;
220                 Vec3 T, dTdx, dTdy;
221                 bool inside;
222                 fresnel_dielectric(m_eta, m, omega_out, domega_out_dx, domega_out_dy,
223                                    R, dRdx, dRdy,
224                                    T, dTdx, dTdy,
225                                    inside);
226
227                 if (!inside) {
228                     omega_in = T;
229                     domega_in_dx = dTdx;
230                     domega_in_dy = dTdy;
231                     // eq. 33
232                     float cosThetaM2 = cosThetaM * cosThetaM;
233                     float cosThetaM4 = cosThetaM2 * cosThetaM2;
234                     float D = alpha2 / (float(M_PI) * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
235                     // eq. 24
236                     float pm = D * cosThetaM;
237                     // eval BRDF*cosNI
238                     float cosNI = m_N.dot(omega_in);
239                     // eq. 34: now calculate G1(i,m) and G1(o,m)
240                     float G1o = 2 / (1 + sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
241                     float G1i = 2 / (1 + sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI))); 
242                     float G = G1o * G1i;
243                     // eq. 21
244                     float cosHI = m.dot(omega_in);
245                     float cosHO = m.dot(omega_out);
246                     float Ht2 = m_eta * cosHI + cosHO;
247                     Ht2 *= Ht2;
248                     float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D)) / (cosNO * Ht2);
249                     // eq. 38 and eq. 17
250                     pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2;
251                     eval.setValue(out, out, out);
252
253                                         /* disabled for now - gives texture filtering problems */
254 #if 0
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;
260                     domega_in_dy *= 10;
261 #endif
262                 }
263             }
264         }
265         return Refractive ? Labels::TRANSMIT : Labels::REFLECT;
266     }
267 };
268
269 // microfacet model with Beckmann facet distribution
270 // see http://www.graphics.cornell.edu/~bjw/microfacetbsdf.pdf
271 template <int Refractive = 0>
272 class MicrofacetBeckmannClosure : public BSDFClosure {
273 public:
274     Vec3 m_N;
275     float m_ab;   // width parameter (roughness)
276     float m_eta;  // index of refraction (for fresnel term)
277     MicrofacetBeckmannClosure() : BSDFClosure(Labels::GLOSSY, Refractive ? Back : Front) { }
278
279     void setup()
280         {
281                 m_ab = clamp(m_ab, 1e-5f, 1.0f);
282         }
283
284     bool mergeable (const ClosurePrimitive *other) const {
285         const MicrofacetBeckmannClosure *comp = (const MicrofacetBeckmannClosure *)other;
286         return m_N == comp->m_N && m_ab == comp->m_ab &&
287             m_eta == comp->m_eta && BSDFClosure::mergeable(other);
288     }
289
290     size_t memsize () const { return sizeof(*this); }
291
292     const char * name () const {
293         return Refractive ? "microfacet_beckmann_refraction" 
294                           : "microfacet_beckmann";
295     }
296
297     void print_on (std::ostream &out) const
298     {
299         out << name() << " (";
300         out << "(" << m_N[0] << ", " << m_N[1] << ", " << m_N[2] << "), ";
301         out << m_ab << ", ";
302         out << m_eta;
303         out << ")";
304     }
305
306     float albedo (const Vec3 &omega_out) const
307     {
308                 return 1.0f;
309     }
310
311     Color3 eval_reflect (const Vec3 &omega_out, const Vec3 &omega_in, float& pdf) const
312     {
313         if (Refractive == 1) return Color3 (0, 0, 0);
314         float cosNO = m_N.dot(omega_out);
315         float cosNI = m_N.dot(omega_in);
316         if (cosNO > 0 && cosNI > 0) {
317            // get half vector
318            Vec3 Hr = omega_in + omega_out;
319            Hr.normalize();
320            // eq. 20: (F*G*D)/(4*in*on)
321            // eq. 25: first we calculate D(m) with m=Hr:
322            float alpha2 = m_ab * m_ab;
323            float cosThetaM = m_N.dot(Hr);
324            float cosThetaM2 = cosThetaM * cosThetaM;
325            float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
326            float cosThetaM4 = cosThetaM2 * cosThetaM2;
327            float D = expf(-tanThetaM2 / alpha2) / (float(M_PI) * alpha2 *  cosThetaM4);
328            // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
329            float ao = 1 / (m_ab * sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
330            float ai = 1 / (m_ab * sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
331            float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
332            float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
333            float G = G1o * G1i;
334            float out = (G * D) * 0.25f / cosNO;
335            // eq. 24
336            float pm = D * cosThetaM;
337            // convert into pdf of the sampled direction
338            // eq. 38 - but see also:
339            // eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
340            pdf = pm * 0.25f / Hr.dot(omega_out);
341            return Color3 (out, out, out);
342         }
343         return Color3 (0, 0, 0);
344     }
345
346     Color3 eval_transmit (const Vec3 &omega_out, const Vec3 &omega_in, float& pdf) const
347     {
348         if (Refractive == 0) return Color3 (0, 0, 0);
349         float cosNO = m_N.dot(omega_out);
350         float cosNI = m_N.dot(omega_in);
351         if (cosNO <= 0 || cosNI >= 0)
352             return Color3 (0, 0, 0);
353         // compute half-vector of the refraction (eq. 16)
354         Vec3 ht = -(m_eta * omega_in + omega_out);
355         Vec3 Ht = ht; Ht.normalize();
356         float cosHO = Ht.dot(omega_out);
357
358                 float cosHI = Ht.dot(omega_in);
359                 // eq. 33: first we calculate D(m) with m=Ht:
360                 float alpha2 = m_ab * m_ab;
361                 float cosThetaM = m_N.dot(Ht);
362                 float cosThetaM2 = cosThetaM * cosThetaM;
363                 float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
364                 float cosThetaM4 = cosThetaM2 * cosThetaM2;
365                 float D = expf(-tanThetaM2 / alpha2) / (float(M_PI) * alpha2 *  cosThetaM4);
366                 // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
367                 float ao = 1 / (m_ab * sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
368                 float ai = 1 / (m_ab * sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
369                 float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
370                 float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
371                 float G = G1o * G1i;
372                 // probability
373                 float invHt2 = 1 / ht.dot(ht);
374                 pdf = D * fabsf(cosThetaM) * (fabsf(cosHI) * (m_eta * m_eta)) * invHt2;
375                 float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D) * invHt2) / cosNO;
376                 return Color3 (out, out, out);
377     }
378
379     ustring sample (const Vec3 &Ng,
380                  const Vec3 &omega_out, const Vec3 &domega_out_dx, const Vec3 &domega_out_dy,
381                  float randu, float randv,
382                  Vec3 &omega_in, Vec3 &domega_in_dx, Vec3 &domega_in_dy,
383                  float &pdf, Color3 &eval) const
384     {
385         float cosNO = m_N.dot(omega_out);
386         if (cosNO > 0) {
387             Vec3 X, Y, Z = m_N;
388             make_orthonormals(Z, X, Y);
389             // generate a random microfacet normal m
390             // eq. 35,36:
391             // we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)
392             //                  and sin(atan(x)) == x/sqrt(1+x^2)
393             float alpha2 = m_ab * m_ab;
394             float tanThetaM = sqrtf(-alpha2 * logf(1 - randu));
395             float cosThetaM = 1 / sqrtf(1 + tanThetaM * tanThetaM);
396             float sinThetaM = cosThetaM * tanThetaM;
397             float phiM = 2 * float(M_PI) * randv;
398             Vec3 m = (cosf(phiM) * sinThetaM) * X +
399                      (sinf(phiM) * sinThetaM) * Y +
400                                    cosThetaM  * Z;
401             if (Refractive == 0) {
402                 float cosMO = m.dot(omega_out);
403                 if (cosMO > 0) {
404                     // eq. 39 - compute actual reflected direction
405                     omega_in = 2 * cosMO * m - omega_out;
406                     if (Ng.dot(omega_in) > 0) {
407                         // microfacet normal is visible to this ray
408                         // eq. 25
409                         float cosThetaM2 = cosThetaM * cosThetaM;
410                         float tanThetaM2 = tanThetaM * tanThetaM;
411                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
412                         float D = expf(-tanThetaM2 / alpha2) / (float(M_PI) * alpha2 *  cosThetaM4);
413                         // eq. 24
414                         float pm = D * cosThetaM;
415                         // convert into pdf of the sampled direction
416                         // eq. 38 - but see also:
417                         // eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
418                         pdf = pm * 0.25f / cosMO;
419                         // Eval BRDF*cosNI
420                         float cosNI = m_N.dot(omega_in);
421                         // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
422                         float ao = 1 / (m_ab * sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
423                         float ai = 1 / (m_ab * sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
424                         float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
425                         float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
426                         float G = G1o * G1i;
427                         // eq. 20: (F*G*D)/(4*in*on)
428                         float out = (G * D) * 0.25f / cosNO;
429                         eval.setValue(out, out, out);
430                         domega_in_dx = (2 * m.dot(domega_out_dx)) * m - domega_out_dx;
431                         domega_in_dy = (2 * m.dot(domega_out_dy)) * m - domega_out_dy;
432
433                                                 /* disabled for now - gives texture filtering problems */
434 #if 0
435                         // Since there is some blur to this reflection, make the
436                         // derivatives a bit bigger. In theory this varies with the
437                         // roughness but the exact relationship is complex and
438                         // requires more ops than are practical.
439                         domega_in_dx *= 10;
440                         domega_in_dy *= 10;
441 #endif
442                     }
443                 }
444             } else {
445                 // CAUTION: the i and o variables are inverted relative to the paper
446                 // eq. 39 - compute actual refractive direction
447                 Vec3 R, dRdx, dRdy;
448                 Vec3 T, dTdx, dTdy;
449                 bool inside;
450                 fresnel_dielectric(m_eta, m, omega_out, domega_out_dx, domega_out_dy,
451                                    R, dRdx, dRdy,
452                                    T, dTdx, dTdy,
453                                    inside);
454                 if (!inside) {
455                     omega_in = T;
456                     domega_in_dx = dTdx;
457                     domega_in_dy = dTdy;
458                     // eq. 33
459                     float cosThetaM2 = cosThetaM * cosThetaM;
460                     float tanThetaM2 = tanThetaM * tanThetaM;
461                     float cosThetaM4 = cosThetaM2 * cosThetaM2;
462                     float D = expf(-tanThetaM2 / alpha2) / (float(M_PI) * alpha2 *  cosThetaM4);
463                     // eq. 24
464                     float pm = D * cosThetaM;
465                     // eval BRDF*cosNI
466                     float cosNI = m_N.dot(omega_in);
467                     // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
468                     float ao = 1 / (m_ab * sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
469                     float ai = 1 / (m_ab * sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
470                     float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
471                     float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
472                     float G = G1o * G1i;
473                     // eq. 21
474                     float cosHI = m.dot(omega_in);
475                     float cosHO = m.dot(omega_out);
476                     float Ht2 = m_eta * cosHI + cosHO;
477                     Ht2 *= Ht2;
478                     float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D)) / (cosNO * Ht2);
479                     // eq. 38 and eq. 17
480                     pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2;
481                     eval.setValue(out, out, out);
482
483                                         /* disabled for now - gives texture filtering problems */
484 #if 0
485                     // Since there is some blur to this refraction, make the
486                     // derivatives a bit bigger. In theory this varies with the
487                     // roughness but the exact relationship is complex and
488                     // requires more ops than are practical.
489                     domega_in_dx *= 10;
490                     domega_in_dy *= 10;
491 #endif
492                 }
493             }
494         }
495         return Refractive ? Labels::TRANSMIT : Labels::REFLECT;
496     }
497 };
498
499
500
501 ClosureParam bsdf_microfacet_ggx_params[] = {
502     CLOSURE_VECTOR_PARAM(MicrofacetGGXClosure<0>, m_N),
503     CLOSURE_FLOAT_PARAM (MicrofacetGGXClosure<0>, m_ag),
504     CLOSURE_STRING_KEYPARAM("label"),
505     CLOSURE_FINISH_PARAM(MicrofacetGGXClosure<0>) };
506
507 ClosureParam bsdf_microfacet_ggx_refraction_params[] = {
508     CLOSURE_VECTOR_PARAM(MicrofacetGGXClosure<1>, m_N),
509     CLOSURE_FLOAT_PARAM (MicrofacetGGXClosure<1>, m_ag),
510     CLOSURE_FLOAT_PARAM (MicrofacetGGXClosure<1>, m_eta),
511     CLOSURE_STRING_KEYPARAM("label"),
512     CLOSURE_FINISH_PARAM(MicrofacetGGXClosure<1>) };
513
514 ClosureParam bsdf_microfacet_beckmann_params[] = {
515     CLOSURE_VECTOR_PARAM(MicrofacetBeckmannClosure<0>, m_N),
516     CLOSURE_FLOAT_PARAM (MicrofacetBeckmannClosure<0>, m_ab),
517     CLOSURE_STRING_KEYPARAM("label"),
518     CLOSURE_FINISH_PARAM(MicrofacetBeckmannClosure<0>) };
519
520 ClosureParam bsdf_microfacet_beckmann_refraction_params[] = {
521     CLOSURE_VECTOR_PARAM(MicrofacetBeckmannClosure<1>, m_N),
522     CLOSURE_FLOAT_PARAM (MicrofacetBeckmannClosure<1>, m_ab),
523     CLOSURE_FLOAT_PARAM (MicrofacetBeckmannClosure<1>, m_eta),
524     CLOSURE_STRING_KEYPARAM("label"),
525     CLOSURE_FINISH_PARAM(MicrofacetBeckmannClosure<1>) };
526
527 CLOSURE_PREPARE(bsdf_microfacet_ggx_prepare,                 MicrofacetGGXClosure<0>)
528 CLOSURE_PREPARE(bsdf_microfacet_ggx_refraction_prepare,      MicrofacetGGXClosure<1>)
529 CLOSURE_PREPARE(bsdf_microfacet_beckmann_prepare,            MicrofacetBeckmannClosure<0>)
530 CLOSURE_PREPARE(bsdf_microfacet_beckmann_refraction_prepare, MicrofacetBeckmannClosure<1>)
531
532 CCL_NAMESPACE_END
533