790368ee888b0cb3e330c07575b944691054881e
[blender.git] / intern / cycles / kernel / closure / bssrdf.h
1 /*
2  * Copyright 2011-2013 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 #ifndef __KERNEL_BSSRDF_H__
18 #define __KERNEL_BSSRDF_H__
19
20 CCL_NAMESPACE_BEGIN
21
22 typedef ccl_addr_space struct Bssrdf {
23         SHADER_CLOSURE_BASE;
24
25         float3 radius;
26         float3 albedo;
27         float sharpness;
28         float texture_blur;
29         float roughness;
30         float channels;
31 } Bssrdf;
32
33 /* Planar Truncated Gaussian
34  *
35  * Note how this is different from the typical gaussian, this one integrates
36  * to 1 over the plane (where you get an extra 2*pi*x factor). We are lucky
37  * that integrating x*exp(-x) gives a nice closed form solution. */
38
39 /* paper suggests 1/12.46 which is much too small, suspect it's *12.46 */
40 #define GAUSS_TRUNCATE 12.46f
41
42 ccl_device float bssrdf_gaussian_eval(const float radius, float r)
43 {
44         /* integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) from 0 to Rm
45          * = 1 - exp(-Rm*Rm/(2*v)) */
46         const float v = radius*radius*(0.25f*0.25f);
47         const float Rm = sqrtf(v*GAUSS_TRUNCATE);
48
49         if(r >= Rm)
50                 return 0.0f;
51
52         return expf(-r*r/(2.0f*v))/(2.0f*M_PI_F*v);
53 }
54
55 ccl_device float bssrdf_gaussian_pdf(const float radius, float r)
56 {
57         /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
58         const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
59
60         return bssrdf_gaussian_eval(radius, r) * (1.0f/(area_truncated));
61 }
62
63 ccl_device void bssrdf_gaussian_sample(const float radius, float xi, float *r, float *h)
64 {
65         /* xi = integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) = -exp(-r^2/(2*v))
66          * r = sqrt(-2*v*logf(xi)) */
67         const float v = radius*radius*(0.25f*0.25f);
68         const float Rm = sqrtf(v*GAUSS_TRUNCATE);
69
70         /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
71         const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
72
73         /* r(xi) */
74         const float r_squared = -2.0f*v*logf(1.0f - xi*area_truncated);
75         *r = sqrtf(r_squared);
76
77         /* h^2 + r^2 = Rm^2 */
78         *h = safe_sqrtf(Rm*Rm - r_squared);
79 }
80
81 /* Planar Cubic BSSRDF falloff
82  *
83  * This is basically (Rm - x)^3, with some factors to normalize it. For sampling
84  * we integrate 2*pi*x * (Rm - x)^3, which gives us a quintic equation that as
85  * far as I can tell has no closed form solution. So we get an iterative solution
86  * instead with newton-raphson. */
87
88 ccl_device float bssrdf_cubic_eval(const float radius, const float sharpness, float r)
89 {
90         if(sharpness == 0.0f) {
91                 const float Rm = radius;
92
93                 if(r >= Rm)
94                         return 0.0f;
95
96                 /* integrate (2*pi*r * 10*(R - r)^3)/(pi * R^5) from 0 to R = 1 */
97                 const float Rm5 = (Rm*Rm) * (Rm*Rm) * Rm;
98                 const float f = Rm - r;
99                 const float num = f*f*f;
100
101                 return (10.0f * num) / (Rm5 * M_PI_F);
102
103         }
104         else {
105                 float Rm = radius*(1.0f + sharpness);
106
107                 if(r >= Rm)
108                         return 0.0f;
109
110                 /* custom variation with extra sharpness, to match the previous code */
111                 const float y = 1.0f/(1.0f + sharpness);
112                 float Rmy, ry, ryinv;
113
114                 if(sharpness == 1.0f) {
115                         Rmy = sqrtf(Rm);
116                         ry = sqrtf(r);
117                         ryinv = (ry > 0.0f)? 1.0f/ry: 0.0f;
118                 }
119                 else {
120                         Rmy = powf(Rm, y);
121                         ry = powf(r, y);
122                         ryinv = (r > 0.0f)? powf(r, y - 1.0f): 0.0f;
123                 }
124
125                 const float Rmy5 = (Rmy*Rmy) * (Rmy*Rmy) * Rmy;
126                 const float f = Rmy - ry;
127                 const float num = f*(f*f)*(y*ryinv);
128
129                 return (10.0f * num) / (Rmy5 * M_PI_F);
130         }
131 }
132
133 ccl_device float bssrdf_cubic_pdf(const float radius, const float sharpness, float r)
134 {
135         return bssrdf_cubic_eval(radius, sharpness, r);
136 }
137
138 /* solve 10x^2 - 20x^3 + 15x^4 - 4x^5 - xi == 0 */
139 ccl_device_forceinline float bssrdf_cubic_quintic_root_find(float xi)
140 {
141         /* newton-raphson iteration, usually succeeds in 2-4 iterations, except
142          * outside 0.02 ... 0.98 where it can go up to 10, so overall performance
143          * should not be too bad */
144         const float tolerance = 1e-6f;
145         const int max_iteration_count = 10;
146         float x = 0.25f;
147         int i;
148
149         for(i = 0; i < max_iteration_count; i++) {
150                 float x2 = x*x;
151                 float x3 = x2*x;
152                 float nx = (1.0f - x);
153
154                 float f = 10.0f*x2 - 20.0f*x3 + 15.0f*x2*x2 - 4.0f*x2*x3 - xi;
155                 float f_ = 20.0f*(x*nx)*(nx*nx);
156
157                 if(fabsf(f) < tolerance || f_ == 0.0f)
158                         break;
159
160                 x = saturate(x - f/f_);
161         }
162
163         return x;
164 }
165
166 ccl_device void bssrdf_cubic_sample(const float radius, const float sharpness, float xi, float *r, float *h)
167 {
168         float Rm = radius;
169         float r_ = bssrdf_cubic_quintic_root_find(xi);
170
171         if(sharpness != 0.0f) {
172                 r_ = powf(r_, 1.0f + sharpness);
173                 Rm *= (1.0f + sharpness);
174         }
175
176         r_ *= Rm;
177         *r = r_;
178
179         /* h^2 + r^2 = Rm^2 */
180         *h = safe_sqrtf(Rm*Rm - r_*r_);
181 }
182
183 /* Approximate Reflectance Profiles
184  * http://graphics.pixar.com/library/ApproxBSSRDF/paper.pdf
185  */
186
187 /* This is a bit arbitrary, just need big enough radius so it matches
188  * the mean free length, but still not too big so sampling is still
189  * effective. Might need some further tweaks.
190  */
191 #define BURLEY_TRUNCATE     16.0f
192 #define BURLEY_TRUNCATE_CDF 0.9963790093708328f // cdf(BURLEY_TRUNCATE)
193
194 ccl_device_inline float bssrdf_burley_fitting(float A)
195 {
196         /* Diffuse surface transmission, equation (6). */
197         return 1.9f - A + 3.5f * (A - 0.8f) * (A - 0.8f);
198 }
199
200 /* Scale mean free path length so it gives similar looking result
201  * to Cubic and Gaussian models.
202  */
203 ccl_device_inline float3 bssrdf_burley_compatible_mfp(float3 r)
204 {
205         return 0.25f * M_1_PI_F * r;
206 }
207
208 ccl_device void bssrdf_burley_setup(Bssrdf *bssrdf)
209 {
210         /* Mean free path length. */
211         const float3 l = bssrdf_burley_compatible_mfp(bssrdf->radius);
212         /* Surface albedo. */
213         const float3 A = bssrdf->albedo;
214         const float3 s = make_float3(bssrdf_burley_fitting(A.x),
215                                  bssrdf_burley_fitting(A.y),
216                                  bssrdf_burley_fitting(A.z));
217
218         bssrdf->radius = l / s;
219 }
220
221 ccl_device float bssrdf_burley_eval(const float d, float r)
222 {
223         const float Rm = BURLEY_TRUNCATE * d;
224
225         if(r >= Rm)
226                 return 0.0f;
227
228         /* Burley refletance profile, equation (3).
229          *
230          * NOTES:
231          * - Surface albedo is already included into sc->weight, no need to
232          *   multiply by this term here.
233          * - This is normalized diffuse model, so the equation is mutliplied
234          *   by 2*pi, which also matches cdf().
235          */
236         float exp_r_3_d = expf(-r / (3.0f * d));
237         float exp_r_d = exp_r_3_d * exp_r_3_d * exp_r_3_d;
238         return (exp_r_d + exp_r_3_d) / (4.0f*d);
239 }
240
241 ccl_device float bssrdf_burley_pdf(const float d, float r)
242 {
243         return bssrdf_burley_eval(d, r) * (1.0f/BURLEY_TRUNCATE_CDF);
244 }
245
246 /* Find the radius for desired CDF value.
247  * Returns scaled radius, meaning the result is to be scaled up by d.
248  * Since there's no closed form solution we do Newton-Raphson method to find it.
249  */
250 ccl_device_forceinline float bssrdf_burley_root_find(float xi)
251 {
252         const float tolerance = 1e-6f;
253         const int max_iteration_count = 10;
254         /* Do initial guess based on manual curve fitting, this allows us to reduce
255          * number of iterations to maximum 4 across the [0..1] range. We keep maximum
256          * number of iteration higher just to be sure we didn't miss root in some
257          * corner case.
258          */
259         float r;
260         if(xi <= 0.9f) {
261                 r = expf(xi * xi * 2.4f) - 1.0f;
262         }
263         else {
264                 /* TODO(sergey): Some nicer curve fit is possible here. */
265                 r = 15.0f;
266         }
267         /* Solve against scaled radius. */
268         for(int i = 0; i < max_iteration_count; i++) {
269                 float exp_r_3 = expf(-r / 3.0f);
270                 float exp_r = exp_r_3 * exp_r_3 * exp_r_3;
271                 float f = 1.0f - 0.25f * exp_r - 0.75f * exp_r_3 - xi;
272                 float f_ = 0.25f * exp_r + 0.25f * exp_r_3;
273
274                 if(fabsf(f) < tolerance || f_ == 0.0f) {
275                         break;
276                 }
277
278                 r = r - f/f_;
279                 if(r < 0.0f) {
280                         r = 0.0f;
281                 }
282         }
283         return r;
284 }
285
286 ccl_device void bssrdf_burley_sample(const float d,
287                                      float xi,
288                                      float *r,
289                                      float *h)
290 {
291         const float Rm = BURLEY_TRUNCATE * d;
292         const float r_ = bssrdf_burley_root_find(xi * BURLEY_TRUNCATE_CDF) * d;
293
294         *r = r_;
295
296         /* h^2 + r^2 = Rm^2 */
297         *h = safe_sqrtf(Rm*Rm - r_*r_);
298 }
299
300 /* None BSSRDF falloff
301  *
302  * Samples distributed over disk with no falloff, for reference. */
303
304 ccl_device float bssrdf_none_eval(const float radius, float r)
305 {
306         const float Rm = radius;
307         return (r < Rm)? 1.0f: 0.0f;
308 }
309
310 ccl_device float bssrdf_none_pdf(const float radius, float r)
311 {
312         /* integrate (2*pi*r)/(pi*Rm*Rm) from 0 to Rm = 1 */
313         const float Rm = radius;
314         const float area = (M_PI_F*Rm*Rm);
315
316         return bssrdf_none_eval(radius, r) / area;
317 }
318
319 ccl_device void bssrdf_none_sample(const float radius, float xi, float *r, float *h)
320 {
321         /* xi = integrate (2*pi*r)/(pi*Rm*Rm) = r^2/Rm^2
322          * r = sqrt(xi)*Rm */
323         const float Rm = radius;
324         const float r_ = sqrtf(xi)*Rm;
325
326         *r = r_;
327
328         /* h^2 + r^2 = Rm^2 */
329         *h = safe_sqrtf(Rm*Rm - r_*r_);
330 }
331
332 /* Generic */
333
334 ccl_device_inline Bssrdf *bssrdf_alloc(ShaderData *sd, float3 weight)
335 {
336         Bssrdf *bssrdf = (Bssrdf*)closure_alloc(sd, sizeof(Bssrdf), CLOSURE_NONE_ID, weight);
337
338         if(bssrdf == NULL) {
339                 return NULL;
340         }
341
342         float sample_weight = fabsf(average(weight));
343         bssrdf->sample_weight = sample_weight;
344         return (sample_weight >= CLOSURE_WEIGHT_CUTOFF) ? bssrdf : NULL;
345 }
346
347 ccl_device int bssrdf_setup(ShaderData *sd, Bssrdf *bssrdf, ClosureType type)
348 {
349         int flag = 0;
350         int bssrdf_channels = 3;
351         float3 diffuse_weight = make_float3(0.0f, 0.0f, 0.0f);
352
353         /* Verify if the radii are large enough to sample without precision issues. */
354         if(bssrdf->radius.x < BSSRDF_MIN_RADIUS) {
355                 diffuse_weight.x = bssrdf->weight.x;
356                 bssrdf->weight.x = 0.0f;
357                 bssrdf->radius.x = 0.0f;
358                 bssrdf_channels--;
359         }
360         if(bssrdf->radius.y < BSSRDF_MIN_RADIUS) {
361                 diffuse_weight.y = bssrdf->weight.y;
362                 bssrdf->weight.y = 0.0f;
363                 bssrdf->radius.y = 0.0f;
364                 bssrdf_channels--;
365         }
366         if(bssrdf->radius.z < BSSRDF_MIN_RADIUS) {
367                 diffuse_weight.z = bssrdf->weight.z;
368                 bssrdf->weight.z = 0.0f;
369                 bssrdf->radius.z = 0.0f;
370                 bssrdf_channels--;
371         }
372
373         if(bssrdf_channels < 3) {
374                 /* Add diffuse BSDF if any radius too small. */
375 #ifdef __PRINCIPLED__
376                 if(type == CLOSURE_BSSRDF_PRINCIPLED_ID) {
377                         float roughness = bssrdf->roughness;
378                         float3 N = bssrdf->N;
379
380                         PrincipledDiffuseBsdf *bsdf = (PrincipledDiffuseBsdf*)bsdf_alloc(sd, sizeof(PrincipledDiffuseBsdf), diffuse_weight);
381
382                         if(bsdf) {
383                                 bsdf->type = CLOSURE_BSDF_BSSRDF_PRINCIPLED_ID;
384                                 bsdf->N = N;
385                                 bsdf->roughness = roughness;
386                                 flag |= bsdf_principled_diffuse_setup(bsdf);
387                         }
388                 }
389                 else
390 #endif  /* __PRINCIPLED__ */
391                 {
392                         DiffuseBsdf *bsdf = (DiffuseBsdf*)bsdf_alloc(sd, sizeof(DiffuseBsdf), diffuse_weight);
393
394                         if(bsdf) {
395                                 bsdf->type = CLOSURE_BSDF_BSSRDF_ID;
396                                 bsdf->N = bssrdf->N;
397                                 flag |= bsdf_diffuse_setup(bsdf);
398                         }
399                 }
400         }
401
402         /* Setup BSSRDF if radius is large enough. */
403         if(bssrdf_channels > 0) {
404                 bssrdf->type = type;
405                 bssrdf->channels = bssrdf_channels;
406                 bssrdf->sample_weight = fabsf(average(bssrdf->weight)) * bssrdf->channels;
407                 bssrdf->texture_blur = saturate(bssrdf->texture_blur);
408                 bssrdf->sharpness = saturate(bssrdf->sharpness);
409
410                 if(type == CLOSURE_BSSRDF_BURLEY_ID ||
411                    type == CLOSURE_BSSRDF_PRINCIPLED_ID ||
412                    type == CLOSURE_BSSRDF_RANDOM_WALK_ID)
413                 {
414                         bssrdf_burley_setup(bssrdf);
415                 }
416
417                 flag |= SD_BSSRDF;
418         }
419         else {
420                 bssrdf->type = type;
421                 bssrdf->sample_weight = 0.0f;
422         }
423
424         return flag;
425 }
426
427 ccl_device void bssrdf_sample(const ShaderClosure *sc, float xi, float *r, float *h)
428 {
429         const Bssrdf *bssrdf = (const Bssrdf*)sc;
430         float radius;
431
432         /* Sample color channel and reuse random number. Only a subset of channels
433          * may be used if their radius was too small to handle as BSSRDF. */
434         xi *= bssrdf->channels;
435
436         if(xi < 1.0f) {
437                 radius = (bssrdf->radius.x > 0.0f)? bssrdf->radius.x:
438                          (bssrdf->radius.y > 0.0f)? bssrdf->radius.y:
439                                                     bssrdf->radius.z;
440         }
441         else if(xi < 2.0f) {
442                 xi -= 1.0f;
443                 radius = (bssrdf->radius.x > 0.0f)? bssrdf->radius.y:
444                                                     bssrdf->radius.z;
445         }
446         else {
447                 xi -= 2.0f;
448                 radius = bssrdf->radius.z;
449         }
450
451         /* Sample BSSRDF. */
452         if(bssrdf->type == CLOSURE_BSSRDF_CUBIC_ID) {
453                 bssrdf_cubic_sample(radius, bssrdf->sharpness, xi, r, h);
454         }
455         else if(bssrdf->type == CLOSURE_BSSRDF_GAUSSIAN_ID){
456                 bssrdf_gaussian_sample(radius, xi, r, h);
457         }
458         else { /*if(bssrdf->type == CLOSURE_BSSRDF_BURLEY_ID || bssrdf->type == CLOSURE_BSSRDF_PRINCIPLED_ID)*/
459                 bssrdf_burley_sample(radius, xi, r, h);
460         }
461 }
462
463 ccl_device float bssrdf_channel_pdf(const Bssrdf *bssrdf, float radius, float r)
464 {
465         if(radius == 0.0f) {
466                 return 0.0f;
467         }
468         else if(bssrdf->type == CLOSURE_BSSRDF_CUBIC_ID) {
469                 return bssrdf_cubic_pdf(radius, bssrdf->sharpness, r);
470         }
471         else if(bssrdf->type == CLOSURE_BSSRDF_GAUSSIAN_ID) {
472                 return bssrdf_gaussian_pdf(radius, r);
473         }
474         else { /*if(bssrdf->type == CLOSURE_BSSRDF_BURLEY_ID || bssrdf->type == CLOSURE_BSSRDF_PRINCIPLED_ID)*/
475                 return bssrdf_burley_pdf(radius, r);
476         }
477 }
478
479 ccl_device_forceinline float3 bssrdf_eval(const ShaderClosure *sc, float r)
480 {
481         const Bssrdf *bssrdf = (const Bssrdf*)sc;
482
483         return make_float3(
484                 bssrdf_channel_pdf(bssrdf, bssrdf->radius.x, r),
485                 bssrdf_channel_pdf(bssrdf, bssrdf->radius.y, r),
486                 bssrdf_channel_pdf(bssrdf, bssrdf->radius.z, r));
487 }
488
489 ccl_device_forceinline float bssrdf_pdf(const ShaderClosure *sc, float r)
490 {
491         const Bssrdf *bssrdf = (const Bssrdf*)sc;
492         float3 pdf = bssrdf_eval(sc, r);
493
494         return (pdf.x + pdf.y + pdf.z) / bssrdf->channels;
495 }
496
497 CCL_NAMESPACE_END
498
499 #endif /* __KERNEL_BSSRDF_H__ */
500