1 /*
2  * Copyright 2011-2013 Blender Foundation
3  *
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
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  */
17 #ifndef __KERNEL_BSSRDF_H__
18 #define __KERNEL_BSSRDF_H__
20 CCL_NAMESPACE_BEGIN
22 typedef ccl_addr_space struct Bssrdf {
26   float3 albedo;
27   float sharpness;
28   float texture_blur;
29   float roughness;
30   float channels;
31 } Bssrdf;
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. */
39 /* paper suggests 1/12.46 which is much too small, suspect it's *12.46 */
40 #define GAUSS_TRUNCATE 12.46f
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);
49   if (r >= Rm)
50     return 0.0f;
52   return expf(-r * r / (2.0f * v)) / (2.0f * M_PI_F * v);
53 }
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);
60   return bssrdf_gaussian_eval(radius, r) * (1.0f / (area_truncated));
61 }
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);
70   /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
71   const float area_truncated = 1.0f - expf(-0.5f * GAUSS_TRUNCATE);
73   /* r(xi) */
74   const float r_squared = -2.0f * v * logf(1.0f - xi * area_truncated);
75   *r = sqrtf(r_squared);
77   /* h^2 + r^2 = Rm^2 */
78   *h = safe_sqrtf(Rm * Rm - r_squared);
79 }
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. */
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;
93     if (r >= Rm)
94       return 0.0f;
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;
101     return (10.0f * num) / (Rm5 * M_PI_F);
102   }
103   else {
104     float Rm = radius * (1.0f + sharpness);
106     if (r >= Rm)
107       return 0.0f;
109     /* custom variation with extra sharpness, to match the previous code */
110     const float y = 1.0f / (1.0f + sharpness);
111     float Rmy, ry, ryinv;
113     if (sharpness == 1.0f) {
114       Rmy = sqrtf(Rm);
115       ry = sqrtf(r);
116       ryinv = (ry > 0.0f) ? 1.0f / ry : 0.0f;
117     }
118     else {
119       Rmy = powf(Rm, y);
120       ry = powf(r, y);
121       ryinv = (r > 0.0f) ? powf(r, y - 1.0f) : 0.0f;
122     }
124     const float Rmy5 = (Rmy * Rmy) * (Rmy * Rmy) * Rmy;
125     const float f = Rmy - ry;
126     const float num = f * (f * f) * (y * ryinv);
128     return (10.0f * num) / (Rmy5 * M_PI_F);
129   }
130 }
132 ccl_device float bssrdf_cubic_pdf(const float radius, const float sharpness, float r)
133 {
135 }
137 /* solve 10x^2 - 20x^3 + 15x^4 - 4x^5 - xi == 0 */
138 ccl_device_forceinline float bssrdf_cubic_quintic_root_find(float xi)
139 {
140   /* newton-raphson iteration, usually succeeds in 2-4 iterations, except
141    * outside 0.02 ... 0.98 where it can go up to 10, so overall performance
142    * should not be too bad */
143   const float tolerance = 1e-6f;
144   const int max_iteration_count = 10;
145   float x = 0.25f;
146   int i;
148   for (i = 0; i < max_iteration_count; i++) {
149     float x2 = x * x;
150     float x3 = x2 * x;
151     float nx = (1.0f - x);
153     float f = 10.0f * x2 - 20.0f * x3 + 15.0f * x2 * x2 - 4.0f * x2 * x3 - xi;
154     float f_ = 20.0f * (x * nx) * (nx * nx);
156     if (fabsf(f) < tolerance || f_ == 0.0f)
157       break;
159     x = saturate(x - f / f_);
160   }
162   return x;
163 }
165 ccl_device void bssrdf_cubic_sample(
166     const float radius, const float sharpness, float xi, float *r, float *h)
167 {
169   float r_ = bssrdf_cubic_quintic_root_find(xi);
171   if (sharpness != 0.0f) {
172     r_ = powf(r_, 1.0f + sharpness);
173     Rm *= (1.0f + sharpness);
174   }
176   r_ *= Rm;
177   *r = r_;
179   /* h^2 + r^2 = Rm^2 */
180   *h = safe_sqrtf(Rm * Rm - r_ * r_);
181 }
183 /* Approximate Reflectance Profiles
184  * http://graphics.pixar.com/library/ApproxBSSRDF/paper.pdf
185  */
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)
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 }
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 }
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(
215       bssrdf_burley_fitting(A.x), bssrdf_burley_fitting(A.y), bssrdf_burley_fitting(A.z));
217   bssrdf->radius = l / s;
218 }
220 ccl_device float bssrdf_burley_eval(const float d, float r)
221 {
222   const float Rm = BURLEY_TRUNCATE * d;
224   if (r >= Rm)
225     return 0.0f;
227   /* Burley refletance profile, equation (3).
228    *
229    * NOTES:
230    * - Surface albedo is already included into sc->weight, no need to
231    *   multiply by this term here.
232    * - This is normalized diffuse model, so the equation is mutliplied
233    *   by 2*pi, which also matches cdf().
234    */
235   float exp_r_3_d = expf(-r / (3.0f * d));
236   float exp_r_d = exp_r_3_d * exp_r_3_d * exp_r_3_d;
237   return (exp_r_d + exp_r_3_d) / (4.0f * d);
238 }
240 ccl_device float bssrdf_burley_pdf(const float d, float r)
241 {
242   return bssrdf_burley_eval(d, r) * (1.0f / BURLEY_TRUNCATE_CDF);
243 }
245 /* Find the radius for desired CDF value.
246  * Returns scaled radius, meaning the result is to be scaled up by d.
247  * Since there's no closed form solution we do Newton-Raphson method to find it.
248  */
249 ccl_device_forceinline float bssrdf_burley_root_find(float xi)
250 {
251   const float tolerance = 1e-6f;
252   const int max_iteration_count = 10;
253   /* Do initial guess based on manual curve fitting, this allows us to reduce
254    * number of iterations to maximum 4 across the [0..1] range. We keep maximum
255    * number of iteration higher just to be sure we didn't miss root in some
256    * corner case.
257    */
258   float r;
259   if (xi <= 0.9f) {
260     r = expf(xi * xi * 2.4f) - 1.0f;
261   }
262   else {
263     /* TODO(sergey): Some nicer curve fit is possible here. */
264     r = 15.0f;
265   }
266   /* Solve against scaled radius. */
267   for (int i = 0; i < max_iteration_count; i++) {
268     float exp_r_3 = expf(-r / 3.0f);
269     float exp_r = exp_r_3 * exp_r_3 * exp_r_3;
270     float f = 1.0f - 0.25f * exp_r - 0.75f * exp_r_3 - xi;
271     float f_ = 0.25f * exp_r + 0.25f * exp_r_3;
273     if (fabsf(f) < tolerance || f_ == 0.0f) {
274       break;
275     }
277     r = r - f / f_;
278     if (r < 0.0f) {
279       r = 0.0f;
280     }
281   }
282   return r;
283 }
285 ccl_device void bssrdf_burley_sample(const float d, float xi, float *r, float *h)
286 {
287   const float Rm = BURLEY_TRUNCATE * d;
288   const float r_ = bssrdf_burley_root_find(xi * BURLEY_TRUNCATE_CDF) * d;
290   *r = r_;
292   /* h^2 + r^2 = Rm^2 */
293   *h = safe_sqrtf(Rm * Rm - r_ * r_);
294 }
296 /* None BSSRDF falloff
297  *
298  * Samples distributed over disk with no falloff, for reference. */
300 ccl_device float bssrdf_none_eval(const float radius, float r)
301 {
302   const float Rm = radius;
303   return (r < Rm) ? 1.0f : 0.0f;
304 }
306 ccl_device float bssrdf_none_pdf(const float radius, float r)
307 {
308   /* integrate (2*pi*r)/(pi*Rm*Rm) from 0 to Rm = 1 */
309   const float Rm = radius;
310   const float area = (M_PI_F * Rm * Rm);
312   return bssrdf_none_eval(radius, r) / area;
313 }
315 ccl_device void bssrdf_none_sample(const float radius, float xi, float *r, float *h)
316 {
317   /* xi = integrate (2*pi*r)/(pi*Rm*Rm) = r^2/Rm^2
318    * r = sqrt(xi)*Rm */
319   const float Rm = radius;
320   const float r_ = sqrtf(xi) * Rm;
322   *r = r_;
324   /* h^2 + r^2 = Rm^2 */
325   *h = safe_sqrtf(Rm * Rm - r_ * r_);
326 }
328 /* Generic */
330 ccl_device_inline Bssrdf *bssrdf_alloc(ShaderData *sd, float3 weight)
331 {
332   Bssrdf *bssrdf = (Bssrdf *)closure_alloc(sd, sizeof(Bssrdf), CLOSURE_NONE_ID, weight);
334   if (bssrdf == NULL) {
335     return NULL;
336   }
338   float sample_weight = fabsf(average(weight));
339   bssrdf->sample_weight = sample_weight;
340   return (sample_weight >= CLOSURE_WEIGHT_CUTOFF) ? bssrdf : NULL;
341 }
343 ccl_device int bssrdf_setup(ShaderData *sd, Bssrdf *bssrdf, ClosureType type)
344 {
345   int flag = 0;
346   int bssrdf_channels = 3;
347   float3 diffuse_weight = make_float3(0.0f, 0.0f, 0.0f);
349   /* Verify if the radii are large enough to sample without precision issues. */
351     diffuse_weight.x = bssrdf->weight.x;
352     bssrdf->weight.x = 0.0f;
354     bssrdf_channels--;
355   }
357     diffuse_weight.y = bssrdf->weight.y;
358     bssrdf->weight.y = 0.0f;
360     bssrdf_channels--;
361   }
363     diffuse_weight.z = bssrdf->weight.z;
364     bssrdf->weight.z = 0.0f;
366     bssrdf_channels--;
367   }
369   if (bssrdf_channels < 3) {
371 #ifdef __PRINCIPLED__
372     if (type == CLOSURE_BSSRDF_PRINCIPLED_ID || type == CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID) {
373       float roughness = bssrdf->roughness;
374       float3 N = bssrdf->N;
376       PrincipledDiffuseBsdf *bsdf = (PrincipledDiffuseBsdf *)bsdf_alloc(
377           sd, sizeof(PrincipledDiffuseBsdf), diffuse_weight);
379       if (bsdf) {
380         bsdf->type = CLOSURE_BSDF_BSSRDF_PRINCIPLED_ID;
381         bsdf->N = N;
382         bsdf->roughness = roughness;
383         flag |= bsdf_principled_diffuse_setup(bsdf);
384       }
385     }
386     else
387 #endif /* __PRINCIPLED__ */
388     {
389       DiffuseBsdf *bsdf = (DiffuseBsdf *)bsdf_alloc(sd, sizeof(DiffuseBsdf), diffuse_weight);
391       if (bsdf) {
392         bsdf->type = CLOSURE_BSDF_BSSRDF_ID;
393         bsdf->N = bssrdf->N;
394         flag |= bsdf_diffuse_setup(bsdf);
395       }
396     }
397   }
399   /* Setup BSSRDF if radius is large enough. */
400   if (bssrdf_channels > 0) {
401     bssrdf->type = type;
402     bssrdf->channels = bssrdf_channels;
403     bssrdf->sample_weight = fabsf(average(bssrdf->weight)) * bssrdf->channels;
404     bssrdf->texture_blur = saturate(bssrdf->texture_blur);
405     bssrdf->sharpness = saturate(bssrdf->sharpness);
407     if (type == CLOSURE_BSSRDF_BURLEY_ID || type == CLOSURE_BSSRDF_PRINCIPLED_ID ||
408         type == CLOSURE_BSSRDF_RANDOM_WALK_ID ||
409         type == CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID) {
410       bssrdf_burley_setup(bssrdf);
411     }
413     flag |= SD_BSSRDF;
414   }
415   else {
416     bssrdf->type = type;
417     bssrdf->sample_weight = 0.0f;
418   }
420   return flag;
421 }
423 ccl_device void bssrdf_sample(const ShaderClosure *sc, float xi, float *r, float *h)
424 {
425   const Bssrdf *bssrdf = (const Bssrdf *)sc;
428   /* Sample color channel and reuse random number. Only a subset of channels
429    * may be used if their radius was too small to handle as BSSRDF. */
430   xi *= bssrdf->channels;
432   if (xi < 1.0f) {
436   }
437   else if (xi < 2.0f) {
438     xi -= 1.0f;
440   }
441   else {
442     xi -= 2.0f;
444   }
446   /* Sample BSSRDF. */
447   if (bssrdf->type == CLOSURE_BSSRDF_CUBIC_ID) {
448     bssrdf_cubic_sample(radius, bssrdf->sharpness, xi, r, h);
449   }
450   else if (bssrdf->type == CLOSURE_BSSRDF_GAUSSIAN_ID) {
452   }
453   else { /*if(bssrdf->type == CLOSURE_BSSRDF_BURLEY_ID || bssrdf->type == CLOSURE_BSSRDF_PRINCIPLED_ID)*/
455   }
456 }
458 ccl_device float bssrdf_channel_pdf(const Bssrdf *bssrdf, float radius, float r)
459 {
460   if (radius == 0.0f) {
461     return 0.0f;
462   }
463   else if (bssrdf->type == CLOSURE_BSSRDF_CUBIC_ID) {
465   }
466   else if (bssrdf->type == CLOSURE_BSSRDF_GAUSSIAN_ID) {
468   }
469   else { /*if(bssrdf->type == CLOSURE_BSSRDF_BURLEY_ID || bssrdf->type == CLOSURE_BSSRDF_PRINCIPLED_ID)*/
471   }
472 }
474 ccl_device_forceinline float3 bssrdf_eval(const ShaderClosure *sc, float r)
475 {
476   const Bssrdf *bssrdf = (const Bssrdf *)sc;