1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
|
/* SPDX-License-Identifier: Apache-2.0
* Copyright 2011-2022 Blender Foundation */
#include "kernel/camera/projection.h"
#include "kernel/bvh/bvh.h"
CCL_NAMESPACE_BEGIN
/* Random walk subsurface scattering.
*
* "Practical and Controllable Subsurface Scattering for Production Path
* Tracing". Matt Jen-Yuan Chiang, Peter Kutz, Brent Burley. SIGGRAPH 2016. */
/* Support for anisotropy from:
* "Path Traced Subsurface Scattering using Anisotropic Phase Functions
* and Non-Exponential Free Flights".
* Magnus Wrenninge, Ryusuke Villemin, Christophe Hery.
* https://graphics.pixar.com/library/PathTracedSubsurface/ */
ccl_device void subsurface_random_walk_remap(const float albedo,
const float d,
float g,
ccl_private float *sigma_t,
ccl_private float *alpha)
{
/* Compute attenuation and scattering coefficients from albedo. */
const float g2 = g * g;
const float g3 = g2 * g;
const float g4 = g3 * g;
const float g5 = g4 * g;
const float g6 = g5 * g;
const float g7 = g6 * g;
const float A = 1.8260523782f + -1.28451056436f * g + -1.79904629312f * g2 +
9.19393289202f * g3 + -22.8215585862f * g4 + 32.0234874259f * g5 +
-23.6264803333f * g6 + 7.21067002658f * g7;
const float B = 4.98511194385f +
0.127355959438f *
expf(31.1491581433f * g + -201.847017512f * g2 + 841.576016723f * g3 +
-2018.09288505f * g4 + 2731.71560286f * g5 + -1935.41424244f * g6 +
559.009054474f * g7);
const float C = 1.09686102424f + -0.394704063468f * g + 1.05258115941f * g2 +
-8.83963712726f * g3 + 28.8643230661f * g4 + -46.8802913581f * g5 +
38.5402837518f * g6 + -12.7181042538f * g7;
const float D = 0.496310210422f + 0.360146581622f * g + -2.15139309747f * g2 +
17.8896899217f * g3 + -55.2984010333f * g4 + 82.065982243f * g5 +
-58.5106008578f * g6 + 15.8478295021f * g7;
const float E = 4.23190299701f +
0.00310603949088f *
expf(76.7316253952f * g + -594.356773233f * g2 + 2448.8834203f * g3 +
-5576.68528998f * g4 + 7116.60171912f * g5 + -4763.54467887f * g6 +
1303.5318055f * g7);
const float F = 2.40602999408f + -2.51814844609f * g + 9.18494908356f * g2 +
-79.2191708682f * g3 + 259.082868209f * g4 + -403.613804597f * g5 +
302.85712436f * g6 + -87.4370473567f * g7;
const float blend = powf(albedo, 0.25f);
*alpha = (1.0f - blend) * A * powf(atanf(B * albedo), C) +
blend * D * powf(atanf(E * albedo), F);
*alpha = clamp(*alpha, 0.0f, 0.999999f); // because of numerical precision
float sigma_t_prime = 1.0f / fmaxf(d, 1e-16f);
*sigma_t = sigma_t_prime / (1.0f - g);
}
ccl_device void subsurface_random_walk_coefficients(const float3 albedo,
const float3 radius,
const float anisotropy,
ccl_private float3 *sigma_t,
ccl_private float3 *alpha,
ccl_private float3 *throughput)
{
float sigma_t_x, sigma_t_y, sigma_t_z;
float alpha_x, alpha_y, alpha_z;
subsurface_random_walk_remap(albedo.x, radius.x, anisotropy, &sigma_t_x, &alpha_x);
subsurface_random_walk_remap(albedo.y, radius.y, anisotropy, &sigma_t_y, &alpha_y);
subsurface_random_walk_remap(albedo.z, radius.z, anisotropy, &sigma_t_z, &alpha_z);
/* Throughput already contains closure weight at this point, which includes the
* albedo, as well as closure mixing and Fresnel weights. Divide out the albedo
* which will be added through scattering. */
*throughput = safe_divide_color(*throughput, albedo);
/* With low albedo values (like 0.025) we get diffusion_length 1.0 and
* infinite phase functions. To avoid a sharp discontinuity as we go from
* such values to 0.0, increase alpha and reduce the throughput to compensate. */
const float min_alpha = 0.2f;
if (alpha_x < min_alpha) {
(*throughput).x *= alpha_x / min_alpha;
alpha_x = min_alpha;
}
if (alpha_y < min_alpha) {
(*throughput).y *= alpha_y / min_alpha;
alpha_y = min_alpha;
}
if (alpha_z < min_alpha) {
(*throughput).z *= alpha_z / min_alpha;
alpha_z = min_alpha;
}
*sigma_t = make_float3(sigma_t_x, sigma_t_y, sigma_t_z);
*alpha = make_float3(alpha_x, alpha_y, alpha_z);
}
/* References for Dwivedi sampling:
*
* [1] "A Zero-variance-based Sampling Scheme for Monte Carlo Subsurface Scattering"
* by Jaroslav Křivánek and Eugene d'Eon (SIGGRAPH 2014)
* https://cgg.mff.cuni.cz/~jaroslav/papers/2014-zerovar/
*
* [2] "Improving the Dwivedi Sampling Scheme"
* by Johannes Meng, Johannes Hanika, and Carsten Dachsbacher (EGSR 2016)
* https://cg.ivd.kit.edu/1951.php
*
* [3] "Zero-Variance Theory for Efficient Subsurface Scattering"
* by Eugene d'Eon and Jaroslav Křivánek (SIGGRAPH 2020)
* https://iliyan.com/publications/RenderingCourse2020
*/
ccl_device_forceinline float eval_phase_dwivedi(float v, float phase_log, float cos_theta)
{
/* Eq. 9 from [2] using precomputed log((v + 1) / (v - 1)) */
return 1.0f / ((v - cos_theta) * phase_log);
}
ccl_device_forceinline float sample_phase_dwivedi(float v, float phase_log, float rand)
{
/* Based on Eq. 10 from [2]: `v - (v + 1) * pow((v - 1) / (v + 1), rand)`
* Since we're already pre-computing `phase_log = log((v + 1) / (v - 1))` for the evaluation,
* we can implement the power function like this. */
return v - (v + 1.0f) * expf(-rand * phase_log);
}
ccl_device_forceinline float diffusion_length_dwivedi(float alpha)
{
/* Eq. 67 from [3] */
return 1.0f / sqrtf(1.0f - powf(alpha, 2.44294f - 0.0215813f * alpha + 0.578637f / alpha));
}
ccl_device_forceinline float3 direction_from_cosine(float3 D, float cos_theta, float randv)
{
float sin_theta = safe_sqrtf(1.0f - cos_theta * cos_theta);
float phi = M_2PI_F * randv;
float3 dir = make_float3(sin_theta * cosf(phi), sin_theta * sinf(phi), cos_theta);
float3 T, B;
make_orthonormals(D, &T, &B);
return dir.x * T + dir.y * B + dir.z * D;
}
ccl_device_forceinline float3 subsurface_random_walk_pdf(float3 sigma_t,
float t,
bool hit,
ccl_private float3 *transmittance)
{
float3 T = volume_color_transmittance(sigma_t, t);
if (transmittance) {
*transmittance = T;
}
return hit ? T : sigma_t * T;
}
/* Define the below variable to get the similarity code active,
* and the value represents the cutoff level */
#define SUBSURFACE_RANDOM_WALK_SIMILARITY_LEVEL 9
ccl_device_inline bool subsurface_random_walk(KernelGlobals kg,
IntegratorState state,
RNGState rng_state,
ccl_private Ray &ray,
ccl_private LocalIntersection &ss_isect)
{
float bssrdf_u, bssrdf_v;
path_state_rng_2D(kg, &rng_state, PRNG_BSDF_U, &bssrdf_u, &bssrdf_v);
const float3 P = INTEGRATOR_STATE(state, ray, P);
const float3 N = INTEGRATOR_STATE(state, ray, D);
const float ray_dP = INTEGRATOR_STATE(state, ray, dP);
const float time = INTEGRATOR_STATE(state, ray, time);
const float3 Ng = INTEGRATOR_STATE(state, subsurface, Ng);
const int object = INTEGRATOR_STATE(state, isect, object);
const int prim = INTEGRATOR_STATE(state, isect, prim);
/* Sample diffuse surface scatter into the object. */
float3 D;
float pdf;
sample_cos_hemisphere(-N, bssrdf_u, bssrdf_v, &D, &pdf);
if (dot(-Ng, D) <= 0.0f) {
return false;
}
/* Setup ray. */
ray.P = P;
ray.D = D;
ray.tmin = 0.0f;
ray.tmax = FLT_MAX;
ray.time = time;
ray.dP = ray_dP;
ray.dD = differential_zero_compact();
ray.self.object = object;
ray.self.prim = prim;
ray.self.light_object = OBJECT_NONE;
ray.self.light_prim = PRIM_NONE;
/* Convert subsurface to volume coefficients.
* The single-scattering albedo is named alpha to avoid confusion with the surface albedo. */
const float3 albedo = INTEGRATOR_STATE(state, subsurface, albedo);
const float3 radius = INTEGRATOR_STATE(state, subsurface, radius);
const float anisotropy = INTEGRATOR_STATE(state, subsurface, anisotropy);
float3 sigma_t, alpha;
float3 throughput = INTEGRATOR_STATE_WRITE(state, path, throughput);
subsurface_random_walk_coefficients(albedo, radius, anisotropy, &sigma_t, &alpha, &throughput);
float3 sigma_s = sigma_t * alpha;
/* Theoretically it should be better to use the exact alpha for the channel we're sampling at
* each bounce, but in practice there doesn't seem to be a noticeable difference in exchange
* for making the code significantly more complex and slower (if direction sampling depends on
* the sampled channel, we need to compute its PDF per-channel and consider it for MIS later on).
*
* Since the strength of the guided sampling increases as alpha gets lower, using a value that
* is too low results in fireflies while one that's too high just gives a bit more noise.
* Therefore, the code here uses the highest of the three albedos to be safe. */
const float diffusion_length = diffusion_length_dwivedi(reduce_max(alpha));
if (diffusion_length == 1.0f) {
/* With specific values of alpha the length might become 1, which in asymptotic makes phase to
* be infinite. After first bounce it will cause throughput to be 0. Do early output, avoiding
* numerical issues and extra unneeded work. */
return false;
}
/* Precompute term for phase sampling. */
const float phase_log = logf((diffusion_length + 1.0f) / (diffusion_length - 1.0f));
/* Modify state for RNGs, decorrelated from other paths. */
rng_state.rng_hash = cmj_hash(rng_state.rng_hash + rng_state.rng_offset, 0xdeadbeef);
/* Random walk until we hit the surface again. */
bool hit = false;
bool have_opposite_interface = false;
float opposite_distance = 0.0f;
/* TODO: Disable for `alpha > 0.999` or so? */
/* Our heuristic, a compromise between guiding and classic. */
const float guided_fraction = 1.0f - fmaxf(0.5f, powf(fabsf(anisotropy), 0.125f));
#ifdef SUBSURFACE_RANDOM_WALK_SIMILARITY_LEVEL
float3 sigma_s_star = sigma_s * (1.0f - anisotropy);
float3 sigma_t_star = sigma_t - sigma_s + sigma_s_star;
float3 sigma_t_org = sigma_t;
float3 sigma_s_org = sigma_s;
const float anisotropy_org = anisotropy;
const float guided_fraction_org = guided_fraction;
#endif
for (int bounce = 0; bounce < BSSRDF_MAX_BOUNCES; bounce++) {
/* Advance random number offset. */
rng_state.rng_offset += PRNG_BOUNCE_NUM;
#ifdef SUBSURFACE_RANDOM_WALK_SIMILARITY_LEVEL
// shadow with local variables according to depth
float anisotropy, guided_fraction;
float3 sigma_s, sigma_t;
if (bounce <= SUBSURFACE_RANDOM_WALK_SIMILARITY_LEVEL) {
anisotropy = anisotropy_org;
guided_fraction = guided_fraction_org;
sigma_t = sigma_t_org;
sigma_s = sigma_s_org;
}
else {
anisotropy = 0.0f;
guided_fraction = 0.75f; // back to isotropic heuristic from Blender
sigma_t = sigma_t_star;
sigma_s = sigma_s_star;
}
#endif
/* Sample color channel, use MIS with balance heuristic. */
float rphase = path_state_rng_1D(kg, &rng_state, PRNG_PHASE_CHANNEL);
float3 channel_pdf;
int channel = volume_sample_channel(alpha, throughput, rphase, &channel_pdf);
float sample_sigma_t = volume_channel_get(sigma_t, channel);
float randt = path_state_rng_1D(kg, &rng_state, PRNG_SCATTER_DISTANCE);
/* We need the result of the ray-cast to compute the full guided PDF, so just remember the
* relevant terms to avoid recomputing them later. */
float backward_fraction = 0.0f;
float forward_pdf_factor = 0.0f;
float forward_stretching = 1.0f;
float backward_pdf_factor = 0.0f;
float backward_stretching = 1.0f;
/* For the initial ray, we already know the direction, so just do classic distance sampling. */
if (bounce > 0) {
/* Decide whether we should use guided or classic sampling. */
bool guided = (path_state_rng_1D(kg, &rng_state, PRNG_LIGHT_TERMINATE) < guided_fraction);
/* Determine if we want to sample away from the incoming interface.
* This only happens if we found a nearby opposite interface, and the probability for it
* depends on how close we are to it already.
* This probability term comes from the recorded presentation of [3]. */
bool guide_backward = false;
if (have_opposite_interface) {
/* Compute distance of the random walk between the tangent plane at the starting point
* and the assumed opposite interface (the parallel plane that contains the point we
* found in our ray query for the opposite side). */
float x = clamp(dot(ray.P - P, -N), 0.0f, opposite_distance);
backward_fraction = 1.0f /
(1.0f + expf((opposite_distance - 2.0f * x) / diffusion_length));
guide_backward = path_state_rng_1D(kg, &rng_state, PRNG_TERMINATE) < backward_fraction;
}
/* Sample scattering direction. */
float scatter_u, scatter_v;
path_state_rng_2D(kg, &rng_state, PRNG_BSDF_U, &scatter_u, &scatter_v);
float cos_theta;
float hg_pdf;
if (guided) {
cos_theta = sample_phase_dwivedi(diffusion_length, phase_log, scatter_u);
/* The backwards guiding distribution is just mirrored along `sd->N`, so swapping the
* sign here is enough to sample from that instead. */
if (guide_backward) {
cos_theta = -cos_theta;
}
float3 newD = direction_from_cosine(N, cos_theta, scatter_v);
hg_pdf = single_peaked_henyey_greenstein(dot(ray.D, newD), anisotropy);
ray.D = newD;
}
else {
float3 newD = henyey_greenstrein_sample(ray.D, anisotropy, scatter_u, scatter_v, &hg_pdf);
cos_theta = dot(newD, N);
ray.D = newD;
}
/* Compute PDF factor caused by phase sampling (as the ratio of guided / classic).
* Since phase sampling is channel-independent, we can get away with applying a factor
* to the guided PDF, which implicitly means pulling out the classic PDF term and letting
* it cancel with an equivalent term in the numerator of the full estimator.
* For the backward PDF, we again reuse the same probability distribution with a sign swap.
*/
forward_pdf_factor = M_1_2PI_F * eval_phase_dwivedi(diffusion_length, phase_log, cos_theta) /
hg_pdf;
backward_pdf_factor = M_1_2PI_F *
eval_phase_dwivedi(diffusion_length, phase_log, -cos_theta) / hg_pdf;
/* Prepare distance sampling.
* For the backwards case, this also needs the sign swapped since now directions against
* `sd->N` (and therefore with negative cos_theta) are preferred. */
forward_stretching = (1.0f - cos_theta / diffusion_length);
backward_stretching = (1.0f + cos_theta / diffusion_length);
if (guided) {
sample_sigma_t *= guide_backward ? backward_stretching : forward_stretching;
}
}
/* Sample direction along ray. */
float t = -logf(1.0f - randt) / sample_sigma_t;
/* On the first bounce, we use the ray-cast to check if the opposite side is nearby.
* If yes, we will later use backwards guided sampling in order to have a decent
* chance of connecting to it.
* TODO: Maybe use less than 10 times the mean free path? */
if (bounce == 0) {
ray.tmax = max(t, 10.0f / (reduce_min(sigma_t)));
}
else {
ray.tmax = t;
/* After the first bounce the object can intersect the same surface again */
ray.self.object = OBJECT_NONE;
ray.self.prim = PRIM_NONE;
}
scene_intersect_local(kg, &ray, &ss_isect, object, NULL, 1);
hit = (ss_isect.num_hits > 0);
if (hit) {
ray.tmax = ss_isect.hits[0].t;
}
if (bounce == 0) {
/* Check if we hit the opposite side. */
if (hit) {
have_opposite_interface = true;
opposite_distance = dot(ray.P + ray.tmax * ray.D - P, -N);
}
/* Apart from the opposite side check, we were supposed to only trace up to distance t,
* so check if there would have been a hit in that case. */
hit = ray.tmax < t;
}
/* Use the distance to the exit point for the throughput update if we found one. */
if (hit) {
t = ray.tmax;
}
/* Advance to new scatter location. */
ray.P += t * ray.D;
float3 transmittance;
float3 pdf = subsurface_random_walk_pdf(sigma_t, t, hit, &transmittance);
if (bounce > 0) {
/* Compute PDF just like we do for classic sampling, but with the stretched sigma_t. */
float3 guided_pdf = subsurface_random_walk_pdf(forward_stretching * sigma_t, t, hit, NULL);
if (have_opposite_interface) {
/* First step of MIS: Depending on geometry we might have two methods for guided
* sampling, so perform MIS between them. */
float3 back_pdf = subsurface_random_walk_pdf(backward_stretching * sigma_t, t, hit, NULL);
guided_pdf = mix(
guided_pdf * forward_pdf_factor, back_pdf * backward_pdf_factor, backward_fraction);
}
else {
/* Just include phase sampling factor otherwise. */
guided_pdf *= forward_pdf_factor;
}
/* Now we apply the MIS balance heuristic between the classic and guided sampling. */
pdf = mix(pdf, guided_pdf, guided_fraction);
}
/* Finally, we're applying MIS again to combine the three color channels.
* Altogether, the MIS computation combines up to nine different estimators:
* {classic, guided, backward_guided} x {r, g, b} */
throughput *= (hit ? transmittance : sigma_s * transmittance) / dot(channel_pdf, pdf);
if (hit) {
/* If we hit the surface, we are done. */
break;
}
else if (throughput.x < VOLUME_THROUGHPUT_EPSILON &&
throughput.y < VOLUME_THROUGHPUT_EPSILON &&
throughput.z < VOLUME_THROUGHPUT_EPSILON) {
/* Avoid unnecessary work and precision issue when throughput gets really small. */
break;
}
}
if (hit) {
kernel_assert(isfinite_safe(throughput));
INTEGRATOR_STATE_WRITE(state, path, throughput) = throughput;
}
return hit;
}
CCL_NAMESPACE_END
|