Fixed remaining syntax errors in OSL files. node_sepcomb_rgb.osl is split into 2...
[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                         }
217                         else {
218                                 // CAUTION: the i and o variables are inverted relative to the paper
219                                 // eq. 39 - compute actual refractive direction
220                                 Vec3 R, dRdx, dRdy;
221                                 Vec3 T, dTdx, dTdy;
222                                 bool inside;
223                                 fresnel_dielectric(m_eta, m, omega_out, domega_out_dx, domega_out_dy,
224                                                    R, dRdx, dRdy,
225                                                    T, dTdx, dTdy,
226                                                    inside);
227
228                                 if (!inside) {
229                                         omega_in = T;
230                                         domega_in_dx = dTdx;
231                                         domega_in_dy = dTdy;
232                                         // eq. 33
233                                         float cosThetaM2 = cosThetaM * cosThetaM;
234                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
235                                         float D = alpha2 / (float(M_PI) * cosThetaM4 * (alpha2 + tanThetaM2) * (alpha2 + tanThetaM2));
236                                         // eq. 24
237                                         float pm = D * cosThetaM;
238                                         // eval BRDF*cosNI
239                                         float cosNI = m_N.dot(omega_in);
240                                         // eq. 34: now calculate G1(i,m) and G1(o,m)
241                                         float G1o = 2 / (1 + sqrtf(1 + alpha2 * (1 - cosNO * cosNO) / (cosNO * cosNO)));
242                                         float G1i = 2 / (1 + sqrtf(1 + alpha2 * (1 - cosNI * cosNI) / (cosNI * cosNI)));
243                                         float G = G1o * G1i;
244                                         // eq. 21
245                                         float cosHI = m.dot(omega_in);
246                                         float cosHO = m.dot(omega_out);
247                                         float Ht2 = m_eta * cosHI + cosHO;
248                                         Ht2 *= Ht2;
249                                         float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D)) / (cosNO * Ht2);
250                                         // eq. 38 and eq. 17
251                                         pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2;
252                                         eval.setValue(out, out, out);
253
254                                         /* disabled for now - gives texture filtering problems */
255 #if 0
256                                         // Since there is some blur to this refraction, make the
257                                         // derivatives a bit bigger. In theory this varies with the
258                                         // roughness but the exact relationship is complex and
259                                         // requires more ops than are practical.
260                                         domega_in_dx *= 10;
261                                         domega_in_dy *= 10;
262 #endif
263                                 }
264                         }
265                 }
266                 return Refractive ? Labels::TRANSMIT : Labels::REFLECT;
267         }
268 };
269
270 // microfacet model with Beckmann facet distribution
271 // see http://www.graphics.cornell.edu/~bjw/microfacetbsdf.pdf
272 template <int Refractive = 0>
273 class MicrofacetBeckmannClosure : public BSDFClosure {
274 public:
275         Vec3 m_N;
276         float m_ab;   // width parameter (roughness)
277         float m_eta;  // index of refraction (for fresnel term)
278         MicrofacetBeckmannClosure() : BSDFClosure(Labels::GLOSSY, Refractive ? Back : Front) {
279         }
280
281         void setup()
282         {
283                 m_ab = clamp(m_ab, 1e-5f, 1.0f);
284         }
285
286         bool mergeable(const ClosurePrimitive *other) const {
287                 const MicrofacetBeckmannClosure *comp = (const MicrofacetBeckmannClosure *)other;
288                 return m_N == comp->m_N && m_ab == comp->m_ab &&
289                        m_eta == comp->m_eta && BSDFClosure::mergeable(other);
290         }
291
292         size_t memsize() const {
293                 return sizeof(*this);
294         }
295
296         const char *name() const {
297                 return Refractive ? "microfacet_beckmann_refraction"
298                            : "microfacet_beckmann";
299         }
300
301         void print_on(std::ostream &out) const
302         {
303                 out << name() << " (";
304                 out << "(" << m_N[0] << ", " << m_N[1] << ", " << m_N[2] << "), ";
305                 out << m_ab << ", ";
306                 out << m_eta;
307                 out << ")";
308         }
309
310         float albedo(const Vec3 &omega_out) const
311         {
312                 return 1.0f;
313         }
314
315         Color3 eval_reflect(const Vec3 &omega_out, const Vec3 &omega_in, float& pdf) const
316         {
317                 if (Refractive == 1) return Color3(0, 0, 0);
318                 float cosNO = m_N.dot(omega_out);
319                 float cosNI = m_N.dot(omega_in);
320                 if (cosNO > 0 && cosNI > 0) {
321                         // get half vector
322                         Vec3 Hr = omega_in + omega_out;
323                         Hr.normalize();
324                         // eq. 20: (F*G*D)/(4*in*on)
325                         // eq. 25: first we calculate D(m) with m=Hr:
326                         float alpha2 = m_ab * m_ab;
327                         float cosThetaM = m_N.dot(Hr);
328                         float cosThetaM2 = cosThetaM * cosThetaM;
329                         float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
330                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
331                         float D = expf(-tanThetaM2 / alpha2) / (float(M_PI) * alpha2 *  cosThetaM4);
332                         // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
333                         float ao = 1 / (m_ab * sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
334                         float ai = 1 / (m_ab * sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
335                         float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
336                         float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
337                         float G = G1o * G1i;
338                         float out = (G * D) * 0.25f / cosNO;
339                         // eq. 24
340                         float pm = D * cosThetaM;
341                         // convert into pdf of the sampled direction
342                         // eq. 38 - but see also:
343                         // eq. 17 in http://www.graphics.cornell.edu/~bjw/wardnotes.pdf
344                         pdf = pm * 0.25f / Hr.dot(omega_out);
345                         return Color3(out, out, out);
346                 }
347                 return Color3(0, 0, 0);
348         }
349
350         Color3 eval_transmit(const Vec3 &omega_out, const Vec3 &omega_in, float& pdf) const
351         {
352                 if (Refractive == 0) return Color3(0, 0, 0);
353                 float cosNO = m_N.dot(omega_out);
354                 float cosNI = m_N.dot(omega_in);
355                 if (cosNO <= 0 || cosNI >= 0)
356                         return Color3(0, 0, 0);
357                 // compute half-vector of the refraction (eq. 16)
358                 Vec3 ht = -(m_eta * omega_in + omega_out);
359                 Vec3 Ht = ht; Ht.normalize();
360                 float cosHO = Ht.dot(omega_out);
361
362                 float cosHI = Ht.dot(omega_in);
363                 // eq. 33: first we calculate D(m) with m=Ht:
364                 float alpha2 = m_ab * m_ab;
365                 float cosThetaM = m_N.dot(Ht);
366                 float cosThetaM2 = cosThetaM * cosThetaM;
367                 float tanThetaM2 = (1 - cosThetaM2) / cosThetaM2;
368                 float cosThetaM4 = cosThetaM2 * cosThetaM2;
369                 float D = expf(-tanThetaM2 / alpha2) / (float(M_PI) * alpha2 *  cosThetaM4);
370                 // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
371                 float ao = 1 / (m_ab * sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
372                 float ai = 1 / (m_ab * sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
373                 float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
374                 float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
375                 float G = G1o * G1i;
376                 // probability
377                 float invHt2 = 1 / ht.dot(ht);
378                 pdf = D * fabsf(cosThetaM) * (fabsf(cosHI) * (m_eta * m_eta)) * invHt2;
379                 float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D) * invHt2) / cosNO;
380                 return Color3(out, out, out);
381         }
382
383         ustring sample(const Vec3 &Ng,
384                        const Vec3 &omega_out, const Vec3 &domega_out_dx, const Vec3 &domega_out_dy,
385                        float randu, float randv,
386                        Vec3 &omega_in, Vec3 &domega_in_dx, Vec3 &domega_in_dy,
387                        float &pdf, Color3 &eval) const
388         {
389                 float cosNO = m_N.dot(omega_out);
390                 if (cosNO > 0) {
391                         Vec3 X, Y, Z = m_N;
392                         make_orthonormals(Z, X, Y);
393                         // generate a random microfacet normal m
394                         // eq. 35,36:
395                         // we take advantage of cos(atan(x)) == 1/sqrt(1+x^2)
396                         //                  and sin(atan(x)) == x/sqrt(1+x^2)
397                         float alpha2 = m_ab * m_ab;
398                         float tanThetaM = sqrtf(-alpha2 * logf(1 - randu));
399                         float cosThetaM = 1 / sqrtf(1 + tanThetaM * tanThetaM);
400                         float sinThetaM = cosThetaM * tanThetaM;
401                         float phiM = 2 * float(M_PI) * randv;
402                         Vec3 m = (cosf(phiM) * sinThetaM) * X +
403                                  (sinf(phiM) * sinThetaM) * Y +
404                                  cosThetaM  * Z;
405                         if (Refractive == 0) {
406                                 float cosMO = m.dot(omega_out);
407                                 if (cosMO > 0) {
408                                         // eq. 39 - compute actual reflected direction
409                                         omega_in = 2 * cosMO * m - omega_out;
410                                         if (Ng.dot(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) / (float(M_PI) * 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 = m_N.dot(omega_in);
425                                                 // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
426                                                 float ao = 1 / (m_ab * sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
427                                                 float ai = 1 / (m_ab * 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.setValue(out, out, out);
434                                                 domega_in_dx = (2 * m.dot(domega_out_dx)) * m - domega_out_dx;
435                                                 domega_in_dy = (2 * m.dot(domega_out_dy)) * m - domega_out_dy;
436
437                                                 /* disabled for now - gives texture filtering problems */
438 #if 0
439                                                 // Since there is some blur to this reflection, make the
440                                                 // derivatives a bit bigger. In theory this varies with the
441                                                 // roughness but the exact relationship is complex and
442                                                 // requires more ops than are practical.
443                                                 domega_in_dx *= 10;
444                                                 domega_in_dy *= 10;
445 #endif
446                                         }
447                                 }
448                         }
449                         else {
450                                 // CAUTION: the i and o variables are inverted relative to the paper
451                                 // eq. 39 - compute actual refractive direction
452                                 Vec3 R, dRdx, dRdy;
453                                 Vec3 T, dTdx, dTdy;
454                                 bool inside;
455                                 fresnel_dielectric(m_eta, m, omega_out, domega_out_dx, domega_out_dy,
456                                                    R, dRdx, dRdy,
457                                                    T, dTdx, dTdy,
458                                                    inside);
459                                 if (!inside) {
460                                         omega_in = T;
461                                         domega_in_dx = dTdx;
462                                         domega_in_dy = dTdy;
463                                         // eq. 33
464                                         float cosThetaM2 = cosThetaM * cosThetaM;
465                                         float tanThetaM2 = tanThetaM * tanThetaM;
466                                         float cosThetaM4 = cosThetaM2 * cosThetaM2;
467                                         float D = expf(-tanThetaM2 / alpha2) / (float(M_PI) * alpha2 *  cosThetaM4);
468                                         // eq. 24
469                                         float pm = D * cosThetaM;
470                                         // eval BRDF*cosNI
471                                         float cosNI = m_N.dot(omega_in);
472                                         // eq. 26, 27: now calculate G1(i,m) and G1(o,m)
473                                         float ao = 1 / (m_ab * sqrtf((1 - cosNO * cosNO) / (cosNO * cosNO)));
474                                         float ai = 1 / (m_ab * sqrtf((1 - cosNI * cosNI) / (cosNI * cosNI)));
475                                         float G1o = ao < 1.6f ? (3.535f * ao + 2.181f * ao * ao) / (1 + 2.276f * ao + 2.577f * ao * ao) : 1.0f;
476                                         float G1i = ai < 1.6f ? (3.535f * ai + 2.181f * ai * ai) / (1 + 2.276f * ai + 2.577f * ai * ai) : 1.0f;
477                                         float G = G1o * G1i;
478                                         // eq. 21
479                                         float cosHI = m.dot(omega_in);
480                                         float cosHO = m.dot(omega_out);
481                                         float Ht2 = m_eta * cosHI + cosHO;
482                                         Ht2 *= Ht2;
483                                         float out = (fabsf(cosHI * cosHO) * (m_eta * m_eta) * (G * D)) / (cosNO * Ht2);
484                                         // eq. 38 and eq. 17
485                                         pdf = pm * (m_eta * m_eta) * fabsf(cosHI) / Ht2;
486                                         eval.setValue(out, out, out);
487
488                                         /* disabled for now - gives texture filtering problems */
489 #if 0
490                                         // Since there is some blur to this refraction, make the
491                                         // derivatives a bit bigger. In theory this varies with the
492                                         // roughness but the exact relationship is complex and
493                                         // requires more ops than are practical.
494                                         domega_in_dx *= 10;
495                                         domega_in_dy *= 10;
496 #endif
497                                 }
498                         }
499                 }
500                 return Refractive ? Labels::TRANSMIT : Labels::REFLECT;
501         }
502 };
503
504
505
506 ClosureParam bsdf_microfacet_ggx_params[] = {
507         CLOSURE_VECTOR_PARAM(MicrofacetGGXClosure<0>, m_N),
508         CLOSURE_FLOAT_PARAM(MicrofacetGGXClosure<0>, m_ag),
509         CLOSURE_STRING_KEYPARAM("label"),
510         CLOSURE_FINISH_PARAM(MicrofacetGGXClosure<0>)
511 };
512
513 ClosureParam bsdf_microfacet_ggx_refraction_params[] = {
514         CLOSURE_VECTOR_PARAM(MicrofacetGGXClosure<1>, m_N),
515         CLOSURE_FLOAT_PARAM(MicrofacetGGXClosure<1>, m_ag),
516         CLOSURE_FLOAT_PARAM(MicrofacetGGXClosure<1>, m_eta),
517         CLOSURE_STRING_KEYPARAM("label"),
518         CLOSURE_FINISH_PARAM(MicrofacetGGXClosure<1>)
519 };
520
521 ClosureParam bsdf_microfacet_beckmann_params[] = {
522         CLOSURE_VECTOR_PARAM(MicrofacetBeckmannClosure<0>, m_N),
523         CLOSURE_FLOAT_PARAM(MicrofacetBeckmannClosure<0>, m_ab),
524         CLOSURE_STRING_KEYPARAM("label"),
525         CLOSURE_FINISH_PARAM(MicrofacetBeckmannClosure<0>)
526 };
527
528 ClosureParam bsdf_microfacet_beckmann_refraction_params[] = {
529         CLOSURE_VECTOR_PARAM(MicrofacetBeckmannClosure<1>, m_N),
530         CLOSURE_FLOAT_PARAM(MicrofacetBeckmannClosure<1>, m_ab),
531         CLOSURE_FLOAT_PARAM(MicrofacetBeckmannClosure<1>, m_eta),
532         CLOSURE_STRING_KEYPARAM("label"),
533         CLOSURE_FINISH_PARAM(MicrofacetBeckmannClosure<1>)
534 };
535
536 CLOSURE_PREPARE(bsdf_microfacet_ggx_prepare,                 MicrofacetGGXClosure<0>)
537 CLOSURE_PREPARE(bsdf_microfacet_ggx_refraction_prepare,      MicrofacetGGXClosure<1>)
538 CLOSURE_PREPARE(bsdf_microfacet_beckmann_prepare,            MicrofacetBeckmannClosure<0>)
539 CLOSURE_PREPARE(bsdf_microfacet_beckmann_refraction_prepare, MicrofacetBeckmannClosure<1>)
540
541 CCL_NAMESPACE_END
542