Blender  V3.3
bsdf_util.h
Go to the documentation of this file.
1 /* SPDX-License-Identifier: BSD-3-Clause
2  *
3  * Adapted from Open Shading Language
4  * Copyright (c) 2009-2010 Sony Pictures Imageworks Inc., et al.
5  * All Rights Reserved.
6  *
7  * Modifications Copyright 2011-2022 Blender Foundation. */
8 
9 #pragma once
10 
12 
14  const float3 N,
15  const float3 I,
19  const float3 dIdx,
20  const float3 dIdy,
21  ccl_private float3 *dRdx,
22  ccl_private float3 *dRdy,
23  ccl_private float3 *dTdx,
24  ccl_private float3 *dTdy,
25 #endif
26  ccl_private bool *is_inside)
27 {
28  float cos = dot(N, I), neta;
29  float3 Nn;
30 
31  // check which side of the surface we are on
32  if (cos > 0) {
33  // we are on the outside of the surface, going in
34  neta = 1 / eta;
35  Nn = N;
36  *is_inside = false;
37  }
38  else {
39  // we are inside the surface
40  cos = -cos;
41  neta = eta;
42  Nn = -N;
43  *is_inside = true;
44  }
45 
46  // compute reflection
47  *R = (2 * cos) * Nn - I;
48 #ifdef __RAY_DIFFERENTIALS__
49  *dRdx = (2 * dot(Nn, dIdx)) * Nn - dIdx;
50  *dRdy = (2 * dot(Nn, dIdy)) * Nn - dIdy;
51 #endif
52 
53  float arg = 1 - (neta * neta * (1 - (cos * cos)));
54  if (arg < 0) {
55  *T = make_float3(0.0f, 0.0f, 0.0f);
56 #ifdef __RAY_DIFFERENTIALS__
57  *dTdx = make_float3(0.0f, 0.0f, 0.0f);
58  *dTdy = make_float3(0.0f, 0.0f, 0.0f);
59 #endif
60  return 1; // total internal reflection
61  }
62  else {
63  float dnp = max(sqrtf(arg), 1e-7f);
64  float nK = (neta * cos) - dnp;
65  *T = -(neta * I) + (nK * Nn);
66 #ifdef __RAY_DIFFERENTIALS__
67  *dTdx = -(neta * dIdx) + ((neta - neta * neta * cos / dnp) * dot(dIdx, Nn)) * Nn;
68  *dTdy = -(neta * dIdy) + ((neta - neta * neta * cos / dnp) * dot(dIdy, Nn)) * Nn;
69 #endif
70  // compute Fresnel terms
71  float cosTheta1 = cos; // N.R
72  float cosTheta2 = -dot(Nn, *T);
73  float pPara = (cosTheta1 - eta * cosTheta2) / (cosTheta1 + eta * cosTheta2);
74  float pPerp = (eta * cosTheta1 - cosTheta2) / (eta * cosTheta1 + cosTheta2);
75  return 0.5f * (pPara * pPara + pPerp * pPerp);
76  }
77 }
78 
79 ccl_device float fresnel_dielectric_cos(float cosi, float eta)
80 {
81  // compute fresnel reflectance without explicitly computing
82  // the refracted direction
83  float c = fabsf(cosi);
84  float g = eta * eta - 1 + c * c;
85  if (g > 0) {
86  g = sqrtf(g);
87  float A = (g - c) / (g + c);
88  float B = (c * (g + c) - 1) / (c * (g - c) + 1);
89  return 0.5f * A * A * (1 + B * B);
90  }
91  return 1.0f; // TIR(no refracted component)
92 }
93 
94 ccl_device float3 fresnel_conductor(float cosi, const float3 eta, const float3 k)
95 {
96  float3 cosi2 = make_float3(cosi * cosi, cosi * cosi, cosi * cosi);
97  float3 one = make_float3(1.0f, 1.0f, 1.0f);
98  float3 tmp_f = eta * eta + k * k;
99  float3 tmp = tmp_f * cosi2;
100  float3 Rparl2 = (tmp - (2.0f * eta * cosi) + one) / (tmp + (2.0f * eta * cosi) + one);
101  float3 Rperp2 = (tmp_f - (2.0f * eta * cosi) + cosi2) / (tmp_f + (2.0f * eta * cosi) + cosi2);
102  return (Rparl2 + Rperp2) * 0.5f;
103 }
104 
106 {
107  float m = clamp(1.0f - u, 0.0f, 1.0f);
108  float m2 = m * m;
109  return m2 * m2 * m; // pow(m, 5)
110 }
111 
112 /* Calculate the fresnel color which is a blend between white and the F0 color (cspec0) */
114 interpolate_fresnel_color(float3 L, float3 H, float ior, float F0, float3 cspec0)
115 {
116  /* Calculate the fresnel interpolation factor
117  * The value from fresnel_dielectric_cos(...) has to be normalized because
118  * the cspec0 keeps the F0 color
119  */
120  float F0_norm = 1.0f / (1.0f - F0);
121  float FH = (fresnel_dielectric_cos(dot(L, H), ior) - F0) * F0_norm;
122 
123  /* Blend between white and a specular color with respect to the fresnel */
124  return cspec0 * (1.0f - FH) + make_float3(1.0f, 1.0f, 1.0f) * FH;
125 }
126 
128 {
129  float3 R = 2 * dot(N, I) * N - I;
130 
131  /* Reflection rays may always be at least as shallow as the incoming ray. */
132  float threshold = min(0.9f * dot(Ng, I), 0.01f);
133  if (dot(Ng, R) >= threshold) {
134  return N;
135  }
136 
137  /* Form coordinate system with Ng as the Z axis and N inside the X-Z-plane.
138  * The X axis is found by normalizing the component of N that's orthogonal to Ng.
139  * The Y axis isn't actually needed.
140  */
141  float NdotNg = dot(N, Ng);
142  float3 X = normalize(N - NdotNg * Ng);
143 
144  /* Keep math expressions. */
145  /* clang-format off */
146  /* Calculate N.z and N.x in the local coordinate system.
147  *
148  * The goal of this computation is to find a N' that is rotated towards Ng just enough
149  * to lift R' above the threshold (here called t), therefore dot(R', Ng) = t.
150  *
151  * According to the standard reflection equation,
152  * this means that we want dot(2*dot(N', I)*N' - I, Ng) = t.
153  *
154  * Since the Z axis of our local coordinate system is Ng, dot(x, Ng) is just x.z, so we get
155  * 2*dot(N', I)*N'.z - I.z = t.
156  *
157  * The rotation is simple to express in the coordinate system we formed -
158  * since N lies in the X-Z-plane, we know that N' will also lie in the X-Z-plane,
159  * so N'.y = 0 and therefore dot(N', I) = N'.x*I.x + N'.z*I.z .
160  *
161  * Furthermore, we want N' to be normalized, so N'.x = sqrt(1 - N'.z^2).
162  *
163  * With these simplifications,
164  * we get the final equation 2*(sqrt(1 - N'.z^2)*I.x + N'.z*I.z)*N'.z - I.z = t.
165  *
166  * The only unknown here is N'.z, so we can solve for that.
167  *
168  * The equation has four solutions in general:
169  *
170  * N'.z = +-sqrt(0.5*(+-sqrt(I.x^2*(I.x^2 + I.z^2 - t^2)) + t*I.z + I.x^2 + I.z^2)/(I.x^2 + I.z^2))
171  * We can simplify this expression a bit by grouping terms:
172  *
173  * a = I.x^2 + I.z^2
174  * b = sqrt(I.x^2 * (a - t^2))
175  * c = I.z*t + a
176  * N'.z = +-sqrt(0.5*(+-b + c)/a)
177  *
178  * Two solutions can immediately be discarded because they're negative so N' would lie in the
179  * lower hemisphere.
180  */
181  /* clang-format on */
182 
183  float Ix = dot(I, X), Iz = dot(I, Ng);
184  float Ix2 = sqr(Ix), Iz2 = sqr(Iz);
185  float a = Ix2 + Iz2;
186 
187  float b = safe_sqrtf(Ix2 * (a - sqr(threshold)));
188  float c = Iz * threshold + a;
189 
190  /* Evaluate both solutions.
191  * In many cases one can be immediately discarded (if N'.z would be imaginary or larger than
192  * one), so check for that first. If no option is viable (might happen in extreme cases like N
193  * being in the wrong hemisphere), give up and return Ng. */
194  float fac = 0.5f / a;
195  float N1_z2 = fac * (b + c), N2_z2 = fac * (-b + c);
196  bool valid1 = (N1_z2 > 1e-5f) && (N1_z2 <= (1.0f + 1e-5f));
197  bool valid2 = (N2_z2 > 1e-5f) && (N2_z2 <= (1.0f + 1e-5f));
198 
199  float2 N_new;
200  if (valid1 && valid2) {
201  /* If both are possible, do the expensive reflection-based check. */
202  float2 N1 = make_float2(safe_sqrtf(1.0f - N1_z2), safe_sqrtf(N1_z2));
203  float2 N2 = make_float2(safe_sqrtf(1.0f - N2_z2), safe_sqrtf(N2_z2));
204 
205  float R1 = 2 * (N1.x * Ix + N1.y * Iz) * N1.y - Iz;
206  float R2 = 2 * (N2.x * Ix + N2.y * Iz) * N2.y - Iz;
207 
208  valid1 = (R1 >= 1e-5f);
209  valid2 = (R2 >= 1e-5f);
210  if (valid1 && valid2) {
211  /* If both solutions are valid, return the one with the shallower reflection since it will be
212  * closer to the input (if the original reflection wasn't shallow, we would not be in this
213  * part of the function). */
214  N_new = (R1 < R2) ? N1 : N2;
215  }
216  else {
217  /* If only one reflection is valid (= positive), pick that one. */
218  N_new = (R1 > R2) ? N1 : N2;
219  }
220  }
221  else if (valid1 || valid2) {
222  /* Only one solution passes the N'.z criterium, so pick that one. */
223  float Nz2 = valid1 ? N1_z2 : N2_z2;
224  N_new = make_float2(safe_sqrtf(1.0f - Nz2), safe_sqrtf(Nz2));
225  }
226  else {
227  return Ng;
228  }
229 
230  return N_new.x * X + N_new.y * Ng;
231 }
232 
MINLINE float safe_sqrtf(float a)
in reality light always falls off quadratically Particle Retrieve the data of the particle that spawned the object for example to give variation to multiple instances of an object Point Retrieve information about points in a point cloud Retrieve the edges of an object as it appears to Cycles topology will always appear triangulated Convert a blackbody temperature to an RGB value Normal Generate a perturbed normal from an RGB normal map image Typically used for faking highly detailed surfaces Generate an OSL shader from a file or text data block Image Sample an image file as a texture Sky Generate a procedural sky texture Noise Generate fractal Perlin noise Wave Generate procedural bands or rings with noise Voronoi Generate Worley noise based on the distance to random points Typically used to generate textures such as or biological cells Brick Generate a procedural texture producing bricks Texture Retrieve multiple types of texture coordinates nTypically used as inputs for texture nodes Vector Convert a or normal between and object coordinate space Combine Create a color from its and value channels Color Retrieve a color or the default fallback if none is specified Separate Split a vector into its X
#define A
ccl_device float3 ensure_valid_reflection(float3 Ng, float3 I, float3 N)
Definition: bsdf_util.h:127
ccl_device float fresnel_dielectric_cos(float cosi, float eta)
Definition: bsdf_util.h:79
CCL_NAMESPACE_BEGIN ccl_device float fresnel_dielectric(float eta, const float3 N, const float3 I, ccl_private float3 *R, ccl_private float3 *T, ccl_private bool *is_inside)
Definition: bsdf_util.h:13
ccl_device_forceinline float3 interpolate_fresnel_color(float3 L, float3 H, float ior, float F0, float3 cspec0)
Definition: bsdf_util.h:114
ccl_device float schlick_fresnel(float u)
Definition: bsdf_util.h:105
ccl_device float3 fresnel_conductor(float cosi, const float3 eta, const float3 k)
Definition: bsdf_util.h:94
#define ccl_device_forceinline
Definition: cuda/compat.h:35
#define ccl_device
Definition: cuda/compat.h:32
#define ccl_private
Definition: cuda/compat.h:48
#define CCL_NAMESPACE_END
Definition: cuda/compat.h:9
static bool is_inside(int x, int y, int cols, int rows)
Definition: filesel.c:706
#define FH(b, c, d)
ccl_gpu_kernel_postfix ccl_global float int int int int float threshold
#define __RAY_DIFFERENTIALS__
Definition: kernel/types.h:62
#define N
#define T
#define B
#define R
#define L
#define H(x, y, z)
#define make_float2(x, y)
Definition: metal/compat.h:203
#define fabsf(x)
Definition: metal/compat.h:219
#define sqrtf(x)
Definition: metal/compat.h:243
#define make_float3(x, y, z)
Definition: metal/compat.h:204
static unsigned c
Definition: RandGen.cpp:83
static unsigned a[3]
Definition: RandGen.cpp:78
INLINE Rall1d< T, V, S > cos(const Rall1d< T, V, S > &arg)
Definition: rall1d.h:319
T dot(const vec_base< T, Size > &a, const vec_base< T, Size > &b)
T clamp(const T &a, const T &min, const T &max)
vec_base< T, Size > normalize(const vec_base< T, Size > &v)
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
static const pxr::TfToken ior("ior", pxr::TfToken::Immortal)
static const pxr::TfToken g("g", pxr::TfToken::Immortal)
#define I
#define min(a, b)
Definition: sort.c:35
float x
Definition: types_float2.h:15
float y
Definition: types_float2.h:15
float max
ccl_device_inline float sqr(float a)
Definition: util/math.h:746