Blender  V3.3
bssrdf.h
Go to the documentation of this file.
1 /* SPDX-License-Identifier: Apache-2.0
2  * Copyright 2011-2022 Blender Foundation */
3 
4 #pragma once
5 
7 
8 typedef struct Bssrdf {
10 
13  float roughness;
14  float anisotropy;
16 
17 static_assert(sizeof(ShaderClosure) >= sizeof(Bssrdf), "Bssrdf is too large!");
18 
19 /* Random Walk BSSRDF */
20 
21 ccl_device float bssrdf_dipole_compute_Rd(float alpha_prime, float fourthirdA)
22 {
23  float s = sqrtf(3.0f * (1.0f - alpha_prime));
24  return 0.5f * alpha_prime * (1.0f + expf(-fourthirdA * s)) * expf(-s);
25 }
26 
27 ccl_device float bssrdf_dipole_compute_alpha_prime(float rd, float fourthirdA)
28 {
29  /* Little Newton solver. */
30  if (rd < 1e-4f) {
31  return 0.0f;
32  }
33  if (rd >= 0.995f) {
34  return 0.999999f;
35  }
36 
37  float x0 = 0.0f;
38  float x1 = 1.0f;
39  float xmid, fmid;
40 
41  constexpr const int max_num_iterations = 12;
42  for (int i = 0; i < max_num_iterations; ++i) {
43  xmid = 0.5f * (x0 + x1);
44  fmid = bssrdf_dipole_compute_Rd(xmid, fourthirdA);
45  if (fmid < rd) {
46  x0 = xmid;
47  }
48  else {
49  x1 = xmid;
50  }
51  }
52 
53  return xmid;
54 }
55 
57  const ClosureType type,
58  const float eta)
59 {
61  /* Scale mean free path length so it gives similar looking result to older
62  * Cubic, Gaussian and Burley models. */
63  bssrdf->radius *= 0.25f * M_1_PI_F;
64  }
65  else {
66  /* Adjust radius based on IOR and albedo. */
67  const float inv_eta = 1.0f / eta;
68  const float F_dr = inv_eta * (-1.440f * inv_eta + 0.710f) + 0.668f + 0.0636f * eta;
69  const float fourthirdA = (4.0f / 3.0f) * (1.0f + F_dr) /
70  (1.0f - F_dr); /* From Jensen's `Fdr` ratio formula. */
71 
72  const float3 alpha_prime = make_float3(
73  bssrdf_dipole_compute_alpha_prime(bssrdf->albedo.x, fourthirdA),
74  bssrdf_dipole_compute_alpha_prime(bssrdf->albedo.y, fourthirdA),
75  bssrdf_dipole_compute_alpha_prime(bssrdf->albedo.z, fourthirdA));
76 
77  bssrdf->radius *= sqrt(3.0f * (one_float3() - alpha_prime));
78  }
79 }
80 
81 /* Christensen-Burley BSSRDF.
82  *
83  * Approximate Reflectance Profiles from
84  * http://graphics.pixar.com/library/ApproxBSSRDF/paper.pdf
85  */
86 
87 /* This is a bit arbitrary, just need big enough radius so it matches
88  * the mean free length, but still not too big so sampling is still
89  * effective. */
90 #define BURLEY_TRUNCATE 16.0f
91 #define BURLEY_TRUNCATE_CDF 0.9963790093708328f // cdf(BURLEY_TRUNCATE)
92 
94 {
95  /* Diffuse surface transmission, equation (6). */
96  return 1.9f - A + 3.5f * (A - 0.8f) * (A - 0.8f);
97 }
98 
99 /* Scale mean free path length so it gives similar looking result
100  * to Cubic and Gaussian models. */
102 {
103  return 0.25f * M_1_PI_F * r;
104 }
105 
107 {
108  /* Mean free path length. */
109  const float3 l = bssrdf_burley_compatible_mfp(bssrdf->radius);
110  /* Surface albedo. */
111  const float3 A = bssrdf->albedo;
112  const float3 s = make_float3(
114 
115  bssrdf->radius = l / s;
116 }
117 
118 ccl_device float bssrdf_burley_eval(const float d, float r)
119 {
120  const float Rm = BURLEY_TRUNCATE * d;
121 
122  if (r >= Rm)
123  return 0.0f;
124 
125  /* Burley reflectance profile, equation (3).
126  *
127  * NOTES:
128  * - Surface albedo is already included into `sc->weight`, no need to
129  * multiply by this term here.
130  * - This is normalized diffuse model, so the equation is multiplied
131  * by `2*pi`, which also matches `cdf()`.
132  */
133  float exp_r_3_d = expf(-r / (3.0f * d));
134  float exp_r_d = exp_r_3_d * exp_r_3_d * exp_r_3_d;
135  return (exp_r_d + exp_r_3_d) / (4.0f * d);
136 }
137 
138 ccl_device float bssrdf_burley_pdf(const float d, float r)
139 {
140  if (r == 0.0f) {
141  return 0.0f;
142  }
143 
144  return bssrdf_burley_eval(d, r) * (1.0f / BURLEY_TRUNCATE_CDF);
145 }
146 
147 /* Find the radius for desired CDF value.
148  * Returns scaled radius, meaning the result is to be scaled up by d.
149  * Since there's no closed form solution we do Newton-Raphson method to find it.
150  */
152 {
153  const float tolerance = 1e-6f;
154  const int max_iteration_count = 10;
155  /* Do initial guess based on manual curve fitting, this allows us to reduce
156  * number of iterations to maximum 4 across the [0..1] range. We keep maximum
157  * number of iteration higher just to be sure we didn't miss root in some
158  * corner case.
159  */
160  float r;
161  if (xi <= 0.9f) {
162  r = expf(xi * xi * 2.4f) - 1.0f;
163  }
164  else {
165  /* TODO(sergey): Some nicer curve fit is possible here. */
166  r = 15.0f;
167  }
168  /* Solve against scaled radius. */
169  for (int i = 0; i < max_iteration_count; i++) {
170  float exp_r_3 = expf(-r / 3.0f);
171  float exp_r = exp_r_3 * exp_r_3 * exp_r_3;
172  float f = 1.0f - 0.25f * exp_r - 0.75f * exp_r_3 - xi;
173  float f_ = 0.25f * exp_r + 0.25f * exp_r_3;
174 
175  if (fabsf(f) < tolerance || f_ == 0.0f) {
176  break;
177  }
178 
179  r = r - f / f_;
180  if (r < 0.0f) {
181  r = 0.0f;
182  }
183  }
184  return r;
185 }
186 
187 ccl_device void bssrdf_burley_sample(const float d,
188  float xi,
189  ccl_private float *r,
190  ccl_private float *h)
191 {
192  const float Rm = BURLEY_TRUNCATE * d;
193  const float r_ = bssrdf_burley_root_find(xi * BURLEY_TRUNCATE_CDF) * d;
194 
195  *r = r_;
196 
197  /* h^2 + r^2 = Rm^2 */
198  *h = safe_sqrtf(Rm * Rm - r_ * r_);
199 }
200 
202 {
203  float channels = 0;
204  if (radius.x > 0.0f) {
205  channels += 1.0f;
206  }
207  if (radius.y > 0.0f) {
208  channels += 1.0f;
209  }
210  if (radius.z > 0.0f) {
211  channels += 1.0f;
212  }
213  return channels;
214 }
215 
216 ccl_device void bssrdf_sample(const float3 radius,
217  float xi,
218  ccl_private float *r,
219  ccl_private float *h)
220 {
221  const float num_channels = bssrdf_num_channels(radius);
222  float sampled_radius;
223 
224  /* Sample color channel and reuse random number. Only a subset of channels
225  * may be used if their radius was too small to handle as BSSRDF. */
226  xi *= num_channels;
227 
228  if (xi < 1.0f) {
229  sampled_radius = (radius.x > 0.0f) ? radius.x : (radius.y > 0.0f) ? radius.y : radius.z;
230  }
231  else if (xi < 2.0f) {
232  xi -= 1.0f;
233  sampled_radius = (radius.x > 0.0f && radius.y > 0.0f) ? radius.y : radius.z;
234  }
235  else {
236  xi -= 2.0f;
237  sampled_radius = radius.z;
238  }
239 
240  /* Sample BSSRDF. */
241  bssrdf_burley_sample(sampled_radius, xi, r, h);
242 }
243 
245 {
246  return make_float3(bssrdf_burley_pdf(radius.x, r),
247  bssrdf_burley_pdf(radius.y, r),
248  bssrdf_burley_pdf(radius.z, r));
249 }
250 
251 ccl_device_forceinline float bssrdf_pdf(const float3 radius, float r)
252 {
253  float3 pdf = bssrdf_eval(radius, r);
254  return (pdf.x + pdf.y + pdf.z) / bssrdf_num_channels(radius);
255 }
256 
257 /* Setup */
258 
260 {
262  sd, sizeof(Bssrdf), CLOSURE_NONE_ID, weight);
263 
264  if (bssrdf == NULL) {
265  return NULL;
266  }
267 
268  float sample_weight = fabsf(average(weight));
269  bssrdf->sample_weight = sample_weight;
270  return (sample_weight >= CLOSURE_WEIGHT_CUTOFF) ? bssrdf : NULL;
271 }
272 
276  const float ior)
277 {
278  int flag = 0;
279 
280  /* Add retro-reflection component as separate diffuse BSDF. */
281  if (bssrdf->roughness != FLT_MAX) {
283  sd, sizeof(PrincipledDiffuseBsdf), bssrdf->weight);
284 
285  if (bsdf) {
286  bsdf->N = bssrdf->N;
287  bsdf->roughness = bssrdf->roughness;
289 
290  /* Ad-hoc weight adjustment to avoid retro-reflection taking away half the
291  * samples from BSSRDF. */
292  bsdf->sample_weight *= bsdf_principled_diffuse_retro_reflection_sample_weight(bsdf, sd->I);
293  }
294  }
295 
296  /* Verify if the radii are large enough to sample without precision issues. */
297  int bssrdf_channels = 3;
298  float3 diffuse_weight = make_float3(0.0f, 0.0f, 0.0f);
299 
300  if (bssrdf->radius.x < BSSRDF_MIN_RADIUS) {
301  diffuse_weight.x = bssrdf->weight.x;
302  bssrdf->weight.x = 0.0f;
303  bssrdf->radius.x = 0.0f;
304  bssrdf_channels--;
305  }
306  if (bssrdf->radius.y < BSSRDF_MIN_RADIUS) {
307  diffuse_weight.y = bssrdf->weight.y;
308  bssrdf->weight.y = 0.0f;
309  bssrdf->radius.y = 0.0f;
310  bssrdf_channels--;
311  }
312  if (bssrdf->radius.z < BSSRDF_MIN_RADIUS) {
313  diffuse_weight.z = bssrdf->weight.z;
314  bssrdf->weight.z = 0.0f;
315  bssrdf->radius.z = 0.0f;
316  bssrdf_channels--;
317  }
318 
319  if (bssrdf_channels < 3) {
320  /* Add diffuse BSDF if any radius too small. */
321 #ifdef __PRINCIPLED__
322  if (bssrdf->roughness != FLT_MAX) {
324  sd, sizeof(PrincipledDiffuseBsdf), diffuse_weight);
325 
326  if (bsdf) {
327  bsdf->N = bssrdf->N;
328  bsdf->roughness = bssrdf->roughness;
330  }
331  }
332  else
333 #endif /* __PRINCIPLED__ */
334  {
336  sd, sizeof(DiffuseBsdf), diffuse_weight);
337 
338  if (bsdf) {
339  bsdf->N = bssrdf->N;
340  flag |= bsdf_diffuse_setup(bsdf);
341  }
342  }
343  }
344 
345  /* Setup BSSRDF if radius is large enough. */
346  if (bssrdf_channels > 0) {
347  bssrdf->type = type;
348  bssrdf->sample_weight = fabsf(average(bssrdf->weight)) * bssrdf_channels;
349 
351 
352  flag |= SD_BSSRDF;
353  }
354  else {
355  bssrdf->type = type;
356  bssrdf->sample_weight = 0.0f;
357  }
358 
359  return flag;
360 }
361 
sqrt(x)+1/max(0
MINLINE float safe_sqrtf(float a)
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble w _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat w _GL_VOID_RET _GL_VOID GLint GLint GLint w _GL_VOID_RET _GL_VOID GLshort GLshort GLshort w _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble y2 _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat y2 _GL_VOID_RET _GL_VOID GLint GLint GLint y2 _GL_VOID_RET _GL_VOID GLshort GLshort GLshort y2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLuint *buffer _GL_VOID_RET _GL_VOID GLdouble t _GL_VOID_RET _GL_VOID GLfloat t _GL_VOID_RET _GL_VOID GLint t _GL_VOID_RET _GL_VOID GLshort t _GL_VOID_RET _GL_VOID GLdouble GLdouble r _GL_VOID_RET _GL_VOID GLfloat GLfloat r _GL_VOID_RET _GL_VOID GLint GLint r _GL_VOID_RET _GL_VOID GLshort GLshort r _GL_VOID_RET _GL_VOID GLdouble GLdouble r
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum type
Group Output data from inside of a node group A color picker Mix two input colors RGB to Convert a color s luminance to a grayscale value Generate a normal vector and a dot product Bright Control the brightness and contrast of the input color Vector Map an input vectors to used to fine tune the interpolation of the input Camera Retrieve information about the camera and how it relates to the current shading point s position Clamp a value between a minimum and a maximum Vector Perform vector math operation Invert a producing a negative Combine Generate a color from its and blue channels(Deprecated)") DefNode(ShaderNode
CCL_NAMESPACE_BEGIN ccl_device ccl_private ShaderClosure * closure_alloc(ccl_private ShaderData *sd, int size, ClosureType type, float3 weight)
Definition: alloc.h:8
ccl_device_inline ccl_private ShaderClosure * bsdf_alloc(ccl_private ShaderData *sd, int size, float3 weight)
Definition: alloc.h:50
ATTR_WARN_UNUSED_RESULT const BMLoop * l
#define A
ccl_device int bsdf_diffuse_setup(ccl_private DiffuseBsdf *bsdf)
Definition: bsdf_diffuse.h:23
ccl_device_inline float bsdf_principled_diffuse_retro_reflection_sample_weight(ccl_private PrincipledDiffuseBsdf *bsdf, const float3 I)
ccl_device int bsdf_principled_diffuse_setup(ccl_private PrincipledDiffuseBsdf *bsdf)
@ PRINCIPLED_DIFFUSE_RETRO_REFLECTION
@ PRINCIPLED_DIFFUSE_LAMBERT
CCL_NAMESPACE_BEGIN struct Bssrdf Bssrdf
ccl_device_forceinline float bssrdf_burley_root_find(float xi)
Definition: bssrdf.h:151
ccl_device float bssrdf_num_channels(const float3 radius)
Definition: bssrdf.h:201
ccl_device_inline ccl_private Bssrdf * bssrdf_alloc(ccl_private ShaderData *sd, float3 weight)
Definition: bssrdf.h:259
ccl_device float bssrdf_burley_pdf(const float d, float r)
Definition: bssrdf.h:138
#define BURLEY_TRUNCATE
Definition: bssrdf.h:90
ccl_device_inline float bssrdf_burley_fitting(float A)
Definition: bssrdf.h:93
ccl_device float bssrdf_dipole_compute_alpha_prime(float rd, float fourthirdA)
Definition: bssrdf.h:27
#define BURLEY_TRUNCATE_CDF
Definition: bssrdf.h:91
ccl_device_inline float3 bssrdf_burley_compatible_mfp(float3 r)
Definition: bssrdf.h:101
ccl_device void bssrdf_setup_radius(ccl_private Bssrdf *bssrdf, const ClosureType type, const float eta)
Definition: bssrdf.h:56
ccl_device_forceinline float3 bssrdf_eval(const float3 radius, float r)
Definition: bssrdf.h:244
ccl_device void bssrdf_burley_setup(ccl_private Bssrdf *bssrdf)
Definition: bssrdf.h:106
ccl_device float bssrdf_dipole_compute_Rd(float alpha_prime, float fourthirdA)
Definition: bssrdf.h:21
ccl_device void bssrdf_burley_sample(const float d, float xi, ccl_private float *r, ccl_private float *h)
Definition: bssrdf.h:187
ccl_device int bssrdf_setup(ccl_private ShaderData *sd, ccl_private Bssrdf *bssrdf, ClosureType type, const float ior)
Definition: bssrdf.h:273
ccl_device float bssrdf_burley_eval(const float d, float r)
Definition: bssrdf.h:118
ccl_device void bssrdf_sample(const float3 radius, float xi, ccl_private float *r, ccl_private float *h)
Definition: bssrdf.h:216
ccl_device_forceinline float bssrdf_pdf(const float3 radius, float r)
Definition: bssrdf.h:251
#define ccl_device_forceinline
Definition: cuda/compat.h:35
#define ccl_device
Definition: cuda/compat.h:32
#define expf(x)
Definition: cuda/compat.h:106
#define ccl_private
Definition: cuda/compat.h:48
#define ccl_device_inline
Definition: cuda/compat.h:34
#define CCL_NAMESPACE_END
Definition: cuda/compat.h:9
ClosureType
@ CLOSURE_BSSRDF_RANDOM_WALK_FIXED_RADIUS_ID
@ CLOSURE_BSSRDF_BURLEY_ID
@ CLOSURE_NONE_ID
#define CLOSURE_WEIGHT_CUTOFF
@ SD_BSSRDF
Definition: kernel/types.h:746
ShaderData
Definition: kernel/types.h:925
ShaderClosure
Definition: kernel/types.h:726
#define BSSRDF_MIN_RADIUS
Definition: kernel/types.h:30
ccl_device_inline float average(const float2 &a)
Definition: math_float2.h:170
ccl_device_inline float3 one_float3()
Definition: math_float3.h:89
#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 const pxr::TfToken ior("ior", pxr::TfToken::Immortal)
closure color bssrdf(string method, normal N, vector radius, color albedo) BUILTIN
Definition: bssrdf.h:8
float anisotropy
Definition: bssrdf.h:14
float roughness
Definition: bssrdf.h:13
float3 albedo
Definition: bssrdf.h:12
SHADER_CLOSURE_BASE
Definition: bssrdf.h:9
float3 radius
Definition: bssrdf.h:11
float z
float y
float x
#define M_1_PI_F
Definition: util/math.h:43