Blender  V3.3
subsurface_disk.h
Go to the documentation of this file.
1 /* SPDX-License-Identifier: Apache-2.0
2  * Copyright 2011-2022 Blender Foundation */
3 
5 
6 /* BSSRDF using disk based importance sampling.
7  *
8  * BSSRDF Importance Sampling, SIGGRAPH 2013
9  * http://library.imageworks.com/pdfs/imageworks-library-BSSRDF-sampling.pdf
10  */
11 
12 ccl_device_inline float3 subsurface_disk_eval(const float3 radius, float disk_r, float r)
13 {
14  const float3 eval = bssrdf_eval(radius, r);
15  const float pdf = bssrdf_pdf(radius, disk_r);
16  return (pdf > 0.0f) ? eval / pdf : zero_float3();
17 }
18 
19 /* Subsurface scattering step, from a point on the surface to other
20  * nearby points on the same object. */
23  RNGState rng_state,
24  ccl_private Ray &ray,
26 
27 {
28  float disk_u, disk_v;
29  path_state_rng_2D(kg, &rng_state, PRNG_BSDF_U, &disk_u, &disk_v);
30 
31  /* Read shading point info from integrator state. */
32  const float3 P = INTEGRATOR_STATE(state, ray, P);
33  const float ray_dP = INTEGRATOR_STATE(state, ray, dP);
34  const float time = INTEGRATOR_STATE(state, ray, time);
35  const float3 Ng = INTEGRATOR_STATE(state, subsurface, Ng);
36  const int object = INTEGRATOR_STATE(state, isect, object);
37  const uint32_t path_flag = INTEGRATOR_STATE(state, path, flag);
38 
39  /* Read subsurface scattering parameters. */
40  const float3 radius = INTEGRATOR_STATE(state, subsurface, radius);
41 
42  /* Pick random axis in local frame and point on disk. */
43  float3 disk_N, disk_T, disk_B;
44  float pick_pdf_N, pick_pdf_T, pick_pdf_B;
45 
46  disk_N = Ng;
47  make_orthonormals(disk_N, &disk_T, &disk_B);
48 
49  if (disk_v < 0.5f) {
50  pick_pdf_N = 0.5f;
51  pick_pdf_T = 0.25f;
52  pick_pdf_B = 0.25f;
53  disk_v *= 2.0f;
54  }
55  else if (disk_v < 0.75f) {
56  float3 tmp = disk_N;
57  disk_N = disk_T;
58  disk_T = tmp;
59  pick_pdf_N = 0.25f;
60  pick_pdf_T = 0.5f;
61  pick_pdf_B = 0.25f;
62  disk_v = (disk_v - 0.5f) * 4.0f;
63  }
64  else {
65  float3 tmp = disk_N;
66  disk_N = disk_B;
67  disk_B = tmp;
68  pick_pdf_N = 0.25f;
69  pick_pdf_T = 0.25f;
70  pick_pdf_B = 0.5f;
71  disk_v = (disk_v - 0.75f) * 4.0f;
72  }
73 
74  /* Sample point on disk. */
75  float phi = M_2PI_F * disk_v;
76  float disk_height, disk_r;
77 
78  bssrdf_sample(radius, disk_u, &disk_r, &disk_height);
79 
80  float3 disk_P = (disk_r * cosf(phi)) * disk_T + (disk_r * sinf(phi)) * disk_B;
81 
82  /* Create ray. */
83  ray.P = P + disk_N * disk_height + disk_P;
84  ray.D = -disk_N;
85  ray.tmin = 0.0f;
86  ray.tmax = 2.0f * disk_height;
87  ray.dP = ray_dP;
88  ray.dD = differential_zero_compact();
89  ray.time = time;
90  ray.self.object = OBJECT_NONE;
91  ray.self.prim = PRIM_NONE;
92  ray.self.light_object = OBJECT_NONE;
93  ray.self.light_prim = OBJECT_NONE;
94 
95  /* Intersect with the same object. if multiple intersections are found it
96  * will use at most BSSRDF_MAX_HITS hits, a random subset of all hits. */
97  uint lcg_state = lcg_state_init(
98  rng_state.rng_hash, rng_state.rng_offset, rng_state.sample, 0x68bc21eb);
99  const int max_hits = BSSRDF_MAX_HITS;
100 
101  scene_intersect_local(kg, &ray, &ss_isect, object, &lcg_state, max_hits);
102  const int num_eval_hits = min(ss_isect.num_hits, max_hits);
103  if (num_eval_hits == 0) {
104  return false;
105  }
106 
107  /* Sort for consistent renders between CPU and GPU, independent of the BVH
108  * traversal algorithm. */
109  sort_intersections_and_normals(ss_isect.hits, ss_isect.Ng, num_eval_hits);
110 
111  float3 weights[BSSRDF_MAX_HITS]; /* TODO: zero? */
112  float sum_weights = 0.0f;
113 
114  for (int hit = 0; hit < num_eval_hits; hit++) {
115  /* Get geometric normal. */
116  const int object = ss_isect.hits[hit].object;
117  const int object_flag = kernel_data_fetch(object_flag, object);
118  float3 hit_Ng = ss_isect.Ng[hit];
119  if (path_flag & PATH_RAY_SUBSURFACE_BACKFACING) {
120  hit_Ng = -hit_Ng;
121  }
122  if (object_flag & SD_OBJECT_NEGATIVE_SCALE_APPLIED) {
123  hit_Ng = -hit_Ng;
124  }
125 
126  if (!(object_flag & SD_OBJECT_TRANSFORM_APPLIED)) {
127  /* Transform normal to world space. */
128  Transform itfm;
129  object_fetch_transform_motion_test(kg, object, time, &itfm);
130  hit_Ng = normalize(transform_direction_transposed(&itfm, hit_Ng));
131  }
132 
133  /* Quickly retrieve P and Ng without setting up ShaderData. */
134  const float3 hit_P = ray.P + ray.D * ss_isect.hits[hit].t;
135 
136  /* Probability densities for local frame axes. */
137  const float pdf_N = pick_pdf_N * fabsf(dot(disk_N, hit_Ng));
138  const float pdf_T = pick_pdf_T * fabsf(dot(disk_T, hit_Ng));
139  const float pdf_B = pick_pdf_B * fabsf(dot(disk_B, hit_Ng));
140 
141  /* Multiple importance sample between 3 axes, power heuristic
142  * found to be slightly better than balance heuristic. pdf_N
143  * in the MIS weight and denominator cancelled out. */
144  float w = pdf_N / (sqr(pdf_N) + sqr(pdf_T) + sqr(pdf_B));
145  if (ss_isect.num_hits > max_hits) {
146  w *= ss_isect.num_hits / (float)max_hits;
147  }
148 
149  /* Real distance to sampled point. */
150  const float r = len(hit_P - P);
151 
152  /* Evaluate profiles. */
153  const float3 weight = subsurface_disk_eval(radius, disk_r, r) * w;
154 
155  /* Store result. */
156  ss_isect.Ng[hit] = hit_Ng;
157  weights[hit] = weight;
158  sum_weights += average(fabs(weight));
159  }
160 
161  if (sum_weights == 0.0f) {
162  return false;
163  }
164 
165  /* Use importance resampling, sampling one of the hits proportional to weight. */
166  const float r = lcg_step_float(&lcg_state) * sum_weights;
167  float partial_sum = 0.0f;
168 
169  for (int hit = 0; hit < num_eval_hits; hit++) {
170  const float3 weight = weights[hit];
171  const float sample_weight = average(fabs(weight));
172  float next_sum = partial_sum + sample_weight;
173 
174  if (r < next_sum) {
175  /* Return exit point. */
176  INTEGRATOR_STATE_WRITE(state, path, throughput) *= weight * sum_weights / sample_weight;
177 
178  ss_isect.hits[0] = ss_isect.hits[hit];
179  ss_isect.Ng[0] = ss_isect.Ng[hit];
180 
181  ray.P = ray.P + ray.D * ss_isect.hits[hit].t;
182  ray.D = ss_isect.Ng[hit];
183  ray.tmin = 0.0f;
184  ray.tmax = 1.0f;
185  return true;
186  }
187 
188  partial_sum = next_sum;
189  }
190 
191  return false;
192 }
193 
typedef float(TangentPoint)[2]
unsigned int uint
Definition: BLI_sys_types.h:67
_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
ccl_device_forceinline float3 bssrdf_eval(const float3 radius, float r)
Definition: bssrdf.h:244
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
SIMD_FORCE_INLINE const btScalar & w() const
Return the w value.
Definition: btQuadWord.h:119
#define sinf(x)
Definition: cuda/compat.h:102
#define cosf(x)
Definition: cuda/compat.h:101
#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
ccl_device_inline void sort_intersections_and_normals(ccl_private Intersection *hits, ccl_private float3 *Ng, uint num_hits)
double time
const KernelGlobalsCPU *ccl_restrict KernelGlobals
#define kernel_data_fetch(name, index)
ccl_device_forceinline float differential_zero_compact()
Definition: differential.h:112
int len
Definition: draw_manager.c:108
ccl_device_inline float3 transform_direction_transposed(ccl_private const Transform *t, const float3 a)
const int state
@ PRNG_BSDF_U
Definition: kernel/types.h:164
#define PRIM_NONE
Definition: kernel/types.h:41
@ PATH_RAY_SUBSURFACE_BACKFACING
Definition: kernel/types.h:261
#define OBJECT_NONE
Definition: kernel/types.h:40
@ SD_OBJECT_TRANSFORM_APPLIED
Definition: kernel/types.h:808
@ SD_OBJECT_NEGATIVE_SCALE_APPLIED
Definition: kernel/types.h:810
#define BSSRDF_MAX_HITS
Definition: kernel/types.h:31
ccl_device_inline uint lcg_state_init(const uint rng_hash, const uint rng_offset, const uint sample, const uint scramble)
Definition: lcg.h:33
ccl_device float lcg_step_float(T rng)
Definition: lcg.h:19
ccl_device_inline float average(const float2 &a)
Definition: math_float2.h:170
ccl_device_inline float2 fabs(const float2 &a)
Definition: math_float2.h:222
ccl_device_inline float3 zero_float3()
Definition: math_float3.h:80
static float P(float k)
Definition: math_interp.c:25
#define fabsf(x)
Definition: metal/compat.h:219
T dot(const vec_base< T, Size > &a, const vec_base< T, Size > &b)
vec_base< T, Size > normalize(const vec_base< T, Size > &v)
ccl_device_inline void path_state_rng_2D(KernelGlobals kg, ccl_private const RNGState *rng_state, int dimension, ccl_private float *fx, ccl_private float *fy)
Definition: path_state.h:307
#define min(a, b)
Definition: sort.c:35
IntegratorStateCPU *ccl_restrict IntegratorState
Definition: state.h:147
#define INTEGRATOR_STATE_WRITE(state, nested_struct, member)
Definition: state.h:155
#define INTEGRATOR_STATE(state, nested_struct, member)
Definition: state.h:154
unsigned int uint32_t
Definition: stdint.h:80
uint rng_hash
Definition: path_state.h:278
int sample
Definition: path_state.h:280
uint rng_offset
Definition: path_state.h:279
CCL_NAMESPACE_BEGIN ccl_device_inline float3 subsurface_disk_eval(const float3 radius, float disk_r, float r)
ccl_device_inline bool subsurface_disk(KernelGlobals kg, IntegratorState state, RNGState rng_state, ccl_private Ray &ray, ccl_private LocalIntersection &ss_isect)
ccl_device_inline float sqr(float a)
Definition: util/math.h:746
#define M_2PI_F
Definition: util/math.h:60
ccl_device_inline void make_orthonormals(const float3 N, ccl_private float3 *a, ccl_private float3 *b)
Definition: util/math.h:566