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