intern/cycles/kernel/closure/bsdf_hair_principled.h

   1 /*
   2  * Copyright 2018 Blender Foundation
   3  *
   4  * Licensed under the Apache License, Version 2.0 (the "License");
   5  * you may not use this file except in compliance with the License.
   6  * You may obtain a copy of the License at
   7  *
   8  * http://www.apache.org/licenses/LICENSE-2.0
   9  *
  10  * Unless required by applicable law or agreed to in writing, software
  11  * distributed under the License is distributed on an "AS IS" BASIS,
  12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  13  * See the License for the specific language governing permissions and
  14  * limitations under the License.
  15  */
  16
  17 #ifdef __KERNEL_CPU__
  18 #include <fenv.h>
  19 #endif
  20
  21 #include "kernel/kernel_color.h"
  22
  23 #ifndef __BSDF_HAIR_PRINCIPLED_H__
  24 #define __BSDF_HAIR_PRINCIPLED_H__
  25
  26 CCL_NAMESPACE_BEGIN
  27
  28 typedef ccl_addr_space struct PrincipledHairExtra {
  29         /* Geometry data. */
  30         float4 geom;
  31 } PrincipledHairExtra;
  32
  33 typedef ccl_addr_space struct PrincipledHairBSDF {
  34         SHADER_CLOSURE_BASE;
  35
  36         /* Absorption coefficient. */
  37         float3 sigma;
  38         /* Variance of the underlying logistic distribution. */
  39         float v;
  40         /* Scale factor of the underlying logistic distribution. */
  41         float s;
  42         /* Cuticle tilt angle. */
  43         float alpha;
  44         /* IOR. */
  45         float eta;
  46         /* Effective variance for the diffuse bounce only. */
  47         float m0_roughness;
  48
  49         /* Extra closure. */
  50         PrincipledHairExtra *extra;
  51 } PrincipledHairBSDF;
  52
  53 static_assert(sizeof(ShaderClosure) >= sizeof(PrincipledHairBSDF), "PrincipledHairBSDF is too large!");
  54 static_assert(sizeof(ShaderClosure) >= sizeof(PrincipledHairExtra), "PrincipledHairExtra is too large!");
  55
  56 ccl_device_inline float cos_from_sin(const float s)
  57 {
  58         return safe_sqrtf(1.0f - s*s);
  59 }
  60
  61 /* Gives the change in direction in the normal plane for the given angles and p-th-order scattering. */
  62 ccl_device_inline float delta_phi(int p, float gamma_o, float gamma_t)
  63 {
  64         return 2.0f * p * gamma_t - 2.0f * gamma_o + p * M_PI_F;
  65 }
  66
  67 /* Remaps the given angle to [-pi, pi]. */
  68 ccl_device_inline float wrap_angle(float a)
  69 {
  70         while(a > M_PI_F) {
  71                 a -= M_2PI_F;
  72         }
  73         while(a < -M_PI_F) {
  74                 a += M_2PI_F;
  75         }
  76         return a;
  77 }
  78
  79 /* Logistic distribution function. */
  80 ccl_device_inline float logistic(float x, float s)
  81 {
  82         float v = expf(-fabsf(x)/s);
  83         return v / (s * sqr(1.0f + v));
  84 }
  85
  86 /* Logistic cumulative density function. */
  87 ccl_device_inline float logistic_cdf(float x, float s)
  88 {
  89         float arg = -x/s;
  90         /* expf() overflows if arg >= 89.0. */
  91         if(arg > 88.0f) {
  92                 return 0.0f;
  93         }
  94         else {
  95                 return 1.0f / (1.0f + expf(arg));
  96         }
  97 }
  98
  99 /* Numerical approximation to the Bessel function of the first kind. */
 100 ccl_device_inline float bessel_I0(float x)
 101 {
 102         x = sqr(x);
 103         float val = 1.0f + 0.25f*x;
 104         float pow_x_2i = sqr(x);
 105         uint64_t i_fac_2 = 1;
 106         int pow_4_i = 16;
 107         for(int i = 2; i < 10; i++) {
 108                 i_fac_2 *= i*i;
 109                 float newval = val + pow_x_2i / (pow_4_i * i_fac_2);
 110                 if(val == newval) {
 111                         return val;
 112                 }
 113                 val = newval;
 114                 pow_x_2i *= x;
 115                 pow_4_i *= 4;
 116         }
 117         return val;
 118 }
 119
 120 /* Logarithm of the Bessel function of the first kind. */
 121 ccl_device_inline float log_bessel_I0(float x)
 122 {
 123         if(x > 12.0f) {
 124                 /* log(1/x) == -log(x) iff x > 0.
 125                  * This is only used with positive cosines */
 126                 return x + 0.5f * (1.f / (8.0f * x) - M_LN_2PI_F - logf(x));
 127         }
 128         else {
 129                 return logf(bessel_I0(x));
 130         }
 131 }
 132
 133 /* Logistic distribution limited to the interval [-pi, pi]. */
 134 ccl_device_inline float trimmed_logistic(float x, float s)
 135 {
 136         /* The logistic distribution is symmetric and centered around zero,
 137          * so logistic_cdf(x, s) = 1 - logistic_cdf(-x, s).
 138          * Therefore, logistic_cdf(x, s)-logistic_cdf(-x, s) = 1 - 2*logistic_cdf(-x, s) */
 139         float scaling_fac = 1.0f - 2.0f*logistic_cdf(-M_PI_F, s);
 140         float val = logistic(x, s);
 141         return safe_divide(val, scaling_fac);
 142 }
 143
 144 /* Sampling function for the trimmed logistic function. */
 145 ccl_device_inline float sample_trimmed_logistic(float u, float s)
 146 {
 147         float cdf_minuspi = logistic_cdf(-M_PI_F, s);
 148         float x = -s*logf(1.0f / (u*(1.0f - 2.0f*cdf_minuspi) + cdf_minuspi) - 1.0f);
 149         return clamp(x, -M_PI_F, M_PI_F);
 150 }
 151
 152 /* Azimuthal scattering function Np. */
 153 ccl_device_inline float azimuthal_scattering(float phi,
 154                                              int p,
 155                                              float s,
 156                                              float gamma_o,
 157                                              float gamma_t)
 158 {
 159         float phi_o = wrap_angle(phi - delta_phi(p, gamma_o, gamma_t));
 160         float val = trimmed_logistic(phi_o, s);
 161         return val;
 162 }
 163
 164 /* Longitudinal scattering function Mp. */
 165 ccl_device_inline float longitudinal_scattering(float sin_theta_i,
 166                                                 float cos_theta_i,
 167                                                 float sin_theta_o,
 168                                                 float cos_theta_o,
 169                                                 float v)
 170 {
 171         float inv_v = 1.0f/v;
 172         float cos_arg = cos_theta_i * cos_theta_o * inv_v;
 173         float sin_arg = sin_theta_i * sin_theta_o * inv_v;
 174         if(v <= 0.1f) {
 175                 float i0 = log_bessel_I0(cos_arg);
 176                 float val = expf(i0 - sin_arg - inv_v + 0.6931f + logf(0.5f*inv_v));
 177                 return val;
 178         }
 179         else {
 180                 float i0 = bessel_I0(cos_arg);
 181                 float val = (expf(-sin_arg) * i0) / (sinhf(inv_v) * 2.0f * v);
 182                 return val;
 183         }
 184 }
 185
 186 /* Combine the three values using their luminances. */
 187 ccl_device_inline float4 combine_with_energy(KernelGlobals *kg, float3 c)
 188 {
 189         return make_float4(c.x, c.y, c.z, linear_rgb_to_gray(kg, c));
 190 }
 191
 192 #ifdef __HAIR__
 193 /* Set up the hair closure. */
 194 ccl_device int bsdf_principled_hair_setup(ShaderData *sd, PrincipledHairBSDF *bsdf)
 195 {
 196         bsdf->type = CLOSURE_BSDF_HAIR_PRINCIPLED_ID;
 197         bsdf->v = clamp(bsdf->v, 0.001f, 1.0f);
 198         bsdf->s = clamp(bsdf->s, 0.001f, 1.0f);
 199         /* Apply Primary Reflection Roughness modifier. */
 200         bsdf->m0_roughness = clamp(bsdf->m0_roughness*bsdf->v, 0.001f, 1.0f);
 201
 202         /* Map from roughness_u and roughness_v to variance and scale factor. */
 203         bsdf->v = sqr(0.726f*bsdf->v + 0.812f*sqr(bsdf->v) + 3.700f*pow20(bsdf->v));
 204         bsdf->s =    (0.265f*bsdf->s + 1.194f*sqr(bsdf->s) + 5.372f*pow22(bsdf->s))*M_SQRT_PI_8_F;
 205         bsdf->m0_roughness = sqr(0.726f*bsdf->m0_roughness + 0.812f*sqr(bsdf->m0_roughness) + 3.700f*pow20(bsdf->m0_roughness));
 206
 207         /* Compute local frame, aligned to curve tangent and ray direction. */
 208         float3 X = safe_normalize(sd->dPdu);
 209         float3 Y = safe_normalize(cross(X, sd->I));
 210         float3 Z = safe_normalize(cross(X, Y));
 211         /* TODO: the solution below works where sd->Ng is the normal
 212          * pointing from the center of the curve to the shading point.
 213          * It doesn't work for triangles, see https://developer.blender.org/T43625 */
 214
 215         /* h -1..0..1 means the rays goes from grazing the hair, to hitting it at
 216          * the center, to grazing the other edge. This is the sine of the angle
 217          * between sd->Ng and Z, as seen from the tangent X. */
 218
 219         /* TODO: we convert this value to a cosine later and discard the sign, so
 220          * we could probably save some operations. */
 221         float h = dot(cross(sd->Ng, X), Z);
 222
 223         kernel_assert(fabsf(h) < 1.0f + 1e-4f);
 224         kernel_assert(isfinite3_safe(Y));
 225         kernel_assert(isfinite_safe(h));
 226
 227         bsdf->extra->geom = make_float4(Y.x, Y.y, Y.z, h);
 228
 229         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_NEEDS_LCG;
 230 }
 231
 232 #endif /* __HAIR__ */
 233
 234 /* Given the Fresnel term and transmittance, generate the attenuation terms for each bounce. */
 235 ccl_device_inline void hair_attenuation(KernelGlobals *kg,
 236                                         float f,
 237                                         float3 T,
 238                                         float4 *Ap)
 239 {
 240         /* Primary specular (R). */
 241         Ap[0] = make_float4(f, f, f, f);
 242
 243         /* Transmission (TT). */
 244         float3 col = sqr(1.0f - f) * T;
 245         Ap[1] = combine_with_energy(kg, col);
 246
 247         /* Secondary specular (TRT). */
 248         col *= T*f;
 249         Ap[2] = combine_with_energy(kg, col);
 250
 251         /* Residual component (TRRT+). */
 252         col *= safe_divide_color(T*f, make_float3(1.0f, 1.0f, 1.0f) - T*f);
 253         Ap[3] = combine_with_energy(kg, col);
 254
 255         /* Normalize sampling weights. */
 256         float totweight = Ap[0].w + Ap[1].w + Ap[2].w + Ap[3].w;
 257         float fac = safe_divide(1.0f, totweight);
 258
 259         Ap[0].w *= fac;
 260         Ap[1].w *= fac;
 261         Ap[2].w *= fac;
 262         Ap[3].w *= fac;
 263 }
 264
 265 /* Given the tilt angle, generate the rotated theta_i for the different bounces. */
 266 ccl_device_inline void hair_alpha_angles(float sin_theta_i,
 267                                          float cos_theta_i,
 268                                          float alpha,
 269                                          float *angles)
 270 {
 271         float sin_1alpha = sinf(alpha);
 272         float cos_1alpha = cos_from_sin(sin_1alpha);
 273         float sin_2alpha = 2.0f*sin_1alpha*cos_1alpha;
 274         float cos_2alpha = sqr(cos_1alpha) - sqr(sin_1alpha);
 275         float sin_4alpha = 2.0f*sin_2alpha*cos_2alpha;
 276         float cos_4alpha = sqr(cos_2alpha) - sqr(sin_2alpha);
 277
 278         angles[0] = sin_theta_i*cos_2alpha + cos_theta_i*sin_2alpha;
 279         angles[1] = fabsf(cos_theta_i*cos_2alpha - sin_theta_i*sin_2alpha);
 280         angles[2] = sin_theta_i*cos_1alpha - cos_theta_i*sin_1alpha;
 281         angles[3] = fabsf(cos_theta_i*cos_1alpha + sin_theta_i*sin_1alpha);
 282         angles[4] = sin_theta_i*cos_4alpha - cos_theta_i*sin_4alpha;
 283         angles[5] = fabsf(cos_theta_i*cos_4alpha + sin_theta_i*sin_4alpha);
 284 }
 285
 286 /* Evaluation function for our shader. */
 287 ccl_device float3 bsdf_principled_hair_eval(KernelGlobals *kg,
 288                                             const ShaderData *sd,
 289                                             const ShaderClosure *sc,
 290                                             const float3 omega_in,
 291                                             float *pdf)
 292 {
 293         kernel_assert(isfinite3_safe(sd->P) && isfinite_safe(sd->ray_length));
 294
 295         const PrincipledHairBSDF *bsdf = (const PrincipledHairBSDF*) sc;
 296         float3 Y = float4_to_float3(bsdf->extra->geom);
 297
 298         float3 X = safe_normalize(sd->dPdu);
 299         kernel_assert(fabsf(dot(X, Y)) < 1e-3f);
 300         float3 Z = safe_normalize(cross(X, Y));
 301
 302         float3 wo = make_float3(dot(sd->I, X), dot(sd->I, Y), dot(sd->I, Z));
 303         float3 wi = make_float3(dot(omega_in, X), dot(omega_in, Y), dot(omega_in, Z));
 304
 305         float sin_theta_o = wo.x;
 306         float cos_theta_o = cos_from_sin(sin_theta_o);
 307         float phi_o = atan2f(wo.z, wo.y);
 308
 309         float sin_theta_t = sin_theta_o / bsdf->eta;
 310         float cos_theta_t = cos_from_sin(sin_theta_t);
 311
 312         float sin_gamma_o = bsdf->extra->geom.w;
 313         float cos_gamma_o = cos_from_sin(sin_gamma_o);
 314         float gamma_o = safe_asinf(sin_gamma_o);
 315
 316         float sin_gamma_t = sin_gamma_o * cos_theta_o / sqrtf(sqr(bsdf->eta) - sqr(sin_theta_o));
 317         float cos_gamma_t = cos_from_sin(sin_gamma_t);
 318         float gamma_t = safe_asinf(sin_gamma_t);
 319
 320         float3 T = exp3(-bsdf->sigma * (2.0f * cos_gamma_t / cos_theta_t));
 321         float4 Ap[4];
 322         hair_attenuation(kg, fresnel_dielectric_cos(cos_theta_o * cos_gamma_o, bsdf->eta), T, Ap);
 323
 324         float sin_theta_i = wi.x;
 325         float cos_theta_i = cos_from_sin(sin_theta_i);
 326         float phi_i = atan2f(wi.z, wi.y);
 327
 328         float phi = phi_i - phi_o;
 329
 330         float angles[6];
 331         hair_alpha_angles(sin_theta_i, cos_theta_i, bsdf->alpha, angles);
 332
 333         float4 F;
 334         float Mp, Np;
 335
 336         /* Primary specular (R). */
 337         Mp = longitudinal_scattering(angles[0], angles[1], sin_theta_o, cos_theta_o, bsdf->m0_roughness);
 338         Np = azimuthal_scattering(phi, 0, bsdf->s, gamma_o, gamma_t);
 339         F  = Ap[0] * Mp * Np;
 340         kernel_assert(isfinite3_safe(float4_to_float3(F)));
 341
 342         /* Transmission (TT). */
 343         Mp = longitudinal_scattering(angles[2], angles[3], sin_theta_o, cos_theta_o, 0.25f*bsdf->v);
 344         Np = azimuthal_scattering(phi, 1, bsdf->s, gamma_o, gamma_t);
 345         F += Ap[1] * Mp * Np;
 346         kernel_assert(isfinite3_safe(float4_to_float3(F)));
 347
 348         /* Secondary specular (TRT). */
 349         Mp = longitudinal_scattering(angles[4], angles[5], sin_theta_o, cos_theta_o, 4.0f*bsdf->v);
 350         Np = azimuthal_scattering(phi, 2, bsdf->s, gamma_o, gamma_t);
 351         F += Ap[2] * Mp * Np;
 352         kernel_assert(isfinite3_safe(float4_to_float3(F)));
 353
 354         /* Residual component (TRRT+). */
 355         Mp = longitudinal_scattering(sin_theta_i, cos_theta_i, sin_theta_o, cos_theta_o, 4.0f*bsdf->v);
 356         Np = M_1_2PI_F;
 357         F += Ap[3] * Mp * Np;
 358         kernel_assert(isfinite3_safe(float4_to_float3(F)));
 359
 360         *pdf = F.w;
 361         return float4_to_float3(F);
 362 }
 363
 364 /* Sampling function for the hair shader. */
 365 ccl_device int bsdf_principled_hair_sample(KernelGlobals *kg,
 366                                            const ShaderClosure *sc,
 367                                            ShaderData *sd,
 368                                            float randu,
 369                                            float randv,
 370                                            float3 *eval,
 371                                            float3 *omega_in,
 372                                            float3 *domega_in_dx,
 373                                            float3 *domega_in_dy,
 374                                            float *pdf)
 375 {
 376         PrincipledHairBSDF *bsdf = (PrincipledHairBSDF*) sc;
 377
 378         float3 Y = float4_to_float3(bsdf->extra->geom);
 379
 380         float3 X = safe_normalize(sd->dPdu);
 381         kernel_assert(fabsf(dot(X, Y)) < 1e-3f);
 382         float3 Z = safe_normalize(cross(X, Y));
 383
 384         float3 wo = make_float3(dot(sd->I, X), dot(sd->I, Y), dot(sd->I, Z));
 385
 386         float2 u[2];
 387         u[0] = make_float2(randu, randv);
 388         u[1].x = lcg_step_float_addrspace(&sd->lcg_state);
 389         u[1].y = lcg_step_float_addrspace(&sd->lcg_state);
 390
 391         float sin_theta_o = wo.x;
 392         float cos_theta_o = cos_from_sin(sin_theta_o);
 393         float phi_o = atan2f(wo.z, wo.y);
 394
 395         float sin_theta_t = sin_theta_o / bsdf->eta;
 396         float cos_theta_t = cos_from_sin(sin_theta_t);
 397
 398         float sin_gamma_o = bsdf->extra->geom.w;
 399         float cos_gamma_o = cos_from_sin(sin_gamma_o);
 400         float gamma_o = safe_asinf(sin_gamma_o);
 401
 402         float sin_gamma_t = sin_gamma_o * cos_theta_o / sqrtf(sqr(bsdf->eta) - sqr(sin_theta_o));
 403         float cos_gamma_t = cos_from_sin(sin_gamma_t);
 404         float gamma_t = safe_asinf(sin_gamma_t);
 405
 406         float3 T = exp3(-bsdf->sigma * (2.0f * cos_gamma_t / cos_theta_t));
 407         float4 Ap[4];
 408         hair_attenuation(kg, fresnel_dielectric_cos(cos_theta_o * cos_gamma_o, bsdf->eta), T, Ap);
 409
 410         int p = 0;
 411         for(; p < 3; p++) {
 412                 if(u[0].x < Ap[p].w) {
 413                         break;
 414                 }
 415                 u[0].x -= Ap[p].w;
 416         }
 417
 418         float v = bsdf->v;
 419         if(p == 1) {
 420                 v *= 0.25f;
 421         }
 422         if(p >= 2) {
 423                 v *= 4.0f;
 424         }
 425
 426         u[1].x = max(u[1].x, 1e-5f);
 427         float fac = 1.0f + v*logf(u[1].x + (1.0f - u[1].x)*expf(-2.0f/v));
 428         float sin_theta_i = -fac * sin_theta_o + cos_from_sin(fac) * cosf(M_2PI_F * u[1].y) * cos_theta_o;
 429         float cos_theta_i = cos_from_sin(sin_theta_i);
 430
 431         float angles[6];
 432         if(p < 3) {
 433                 hair_alpha_angles(sin_theta_i, cos_theta_i, -bsdf->alpha, angles);
 434                 sin_theta_i = angles[2*p];
 435                 cos_theta_i = angles[2*p+1];
 436         }
 437
 438         float phi;
 439         if(p < 3) {
 440                 phi = delta_phi(p, gamma_o, gamma_t) + sample_trimmed_logistic(u[0].y, bsdf->s);
 441         }
 442         else {
 443                 phi = M_2PI_F*u[0].y;
 444         }
 445         float phi_i = phi_o + phi;
 446
 447         hair_alpha_angles(sin_theta_i, cos_theta_i, bsdf->alpha, angles);
 448
 449         float4 F;
 450         float Mp, Np;
 451
 452         /* Primary specular (R). */
 453         Mp = longitudinal_scattering(angles[0], angles[1], sin_theta_o, cos_theta_o, bsdf->m0_roughness);
 454         Np = azimuthal_scattering(phi, 0, bsdf->s, gamma_o, gamma_t);
 455         F  = Ap[0] * Mp * Np;
 456         kernel_assert(isfinite3_safe(float4_to_float3(F)));
 457
 458         /* Transmission (TT). */
 459         Mp = longitudinal_scattering(angles[2], angles[3], sin_theta_o, cos_theta_o, 0.25f*bsdf->v);
 460         Np = azimuthal_scattering(phi, 1, bsdf->s, gamma_o, gamma_t);
 461         F += Ap[1] * Mp * Np;
 462         kernel_assert(isfinite3_safe(float4_to_float3(F)));
 463
 464         /* Secondary specular (TRT). */
 465         Mp = longitudinal_scattering(angles[4], angles[5], sin_theta_o, cos_theta_o, 4.0f*bsdf->v);
 466         Np = azimuthal_scattering(phi, 2, bsdf->s, gamma_o, gamma_t);
 467         F += Ap[2] * Mp * Np;
 468         kernel_assert(isfinite3_safe(float4_to_float3(F)));
 469
 470         /* Residual component (TRRT+). */
 471         Mp = longitudinal_scattering(sin_theta_i, cos_theta_i, sin_theta_o, cos_theta_o, 4.0f*bsdf->v);
 472         Np = M_1_2PI_F;
 473         F += Ap[3] * Mp * Np;
 474         kernel_assert(isfinite3_safe(float4_to_float3(F)));
 475
 476         *eval = float4_to_float3(F);
 477         *pdf = F.w;
 478
 479         *omega_in = X*sin_theta_i + Y*cos_theta_i*cosf(phi_i) + Z*cos_theta_i*sinf(phi_i);
 480
 481 #ifdef __RAY_DIFFERENTIALS__
 482         float3 N = safe_normalize(sd->I + *omega_in);
 483         *domega_in_dx = (2 * dot(N, sd->dI.dx)) * N - sd->dI.dx;
 484         *domega_in_dy = (2 * dot(N, sd->dI.dy)) * N - sd->dI.dy;
 485 #endif
 486
 487         return LABEL_GLOSSY|((p == 0)? LABEL_REFLECT : LABEL_TRANSMIT);
 488 }
 489
 490 /* Implements Filter Glossy by capping the effective roughness. */
 491 ccl_device void bsdf_principled_hair_blur(ShaderClosure *sc, float roughness)
 492 {
 493         PrincipledHairBSDF *bsdf = (PrincipledHairBSDF*)sc;
 494
 495         bsdf->v = fmaxf(roughness, bsdf->v);
 496         bsdf->s = fmaxf(roughness, bsdf->s);
 497         bsdf->m0_roughness = fmaxf(roughness, bsdf->m0_roughness);
 498 }
 499
 500 CCL_NAMESPACE_END
 501
 502 #endif /* __BSDF_HAIR_PRINCIPLED_H__ */