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 /* Planar Truncated Gaussian
23  *
24  * Note how this is different from the typical gaussian, this one integrates
25  * to 1 over the plane (where you get an extra 2*pi*x factor). We are lucky
26  * that integrating x*exp(-x) gives a nice closed form solution. */
28 /* paper suggests 1/12.46 which is much too small, suspect it's *12.46 */
29 #define GAUSS_TRUNCATE 12.46f
31 ccl_device float bssrdf_gaussian_eval(ShaderClosure *sc, float r)
32 {
33         /* integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) from 0 to Rm
34          * = 1 - exp(-Rm*Rm/(2*v)) */
35         const float v = sc->data0*sc->data0*(0.25f*0.25f);
36         const float Rm = sqrtf(v*GAUSS_TRUNCATE);
38         if(r >= Rm)
39                 return 0.0f;
41         return expf(-r*r/(2.0f*v))/(2.0f*M_PI_F*v);
42 }
44 ccl_device float bssrdf_gaussian_pdf(ShaderClosure *sc, float r)
45 {
46         /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
47         const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
49         return bssrdf_gaussian_eval(sc, r) * (1.0f/(area_truncated));
50 }
52 ccl_device void bssrdf_gaussian_sample(ShaderClosure *sc, float xi, float *r, float *h)
53 {
54         /* xi = integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) = -exp(-r^2/(2*v))
55          * r = sqrt(-2*v*logf(xi)) */
57         const float v = sc->data0*sc->data0*(0.25f*0.25f);
58         const float Rm = sqrtf(v*GAUSS_TRUNCATE);
60         /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
61         const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
63         /* r(xi) */
64         const float r_squared = -2.0f*v*logf(1.0f - xi*area_truncated);
65         *r = sqrtf(r_squared);
67         /* h^2 + r^2 = Rm^2 */
68         *h = safe_sqrtf(Rm*Rm - r_squared);
69 }
71 /* Planar Cubic BSSRDF falloff
72  *
73  * This is basically (Rm - x)^3, with some factors to normalize it. For sampling
74  * we integrate 2*pi*x * (Rm - x)^3, which gives us a quintic equation that as
75  * far as I can tell has no closed form solution. So we get an iterative solution
76  * instead with newton-raphson. */
78 ccl_device float bssrdf_cubic_eval(ShaderClosure *sc, float r)
79 {
80         const float sharpness = sc->T.x;
82         if(sharpness == 0.0f) {
83                 const float Rm = sc->data0;
85                 if(r >= Rm)
86                         return 0.0f;
88                 /* integrate (2*pi*r * 10*(R - r)^3)/(pi * R^5) from 0 to R = 1 */
89                 const float Rm5 = (Rm*Rm) * (Rm*Rm) * Rm;
90                 const float f = Rm - r;
91                 const float num = f*f*f;
93                 return (10.0f * num) / (Rm5 * M_PI_F);
95         }
96         else {
97                 float Rm = sc->data0*(1.0f + sharpness);
99                 if(r >= Rm)
100                         return 0.0f;
102                 /* custom variation with extra sharpness, to match the previous code */
103                 const float y = 1.0f/(1.0f + sharpness);
104                 float Rmy, ry, ryinv;
106                 if(sharpness == 1.0f) {
107                         Rmy = sqrtf(Rm);
108                         ry = sqrtf(r);
109                         ryinv = (ry > 0.0f)? 1.0f/ry: 0.0f;
110                 }
111                 else {
112                         Rmy = powf(Rm, y);
113                         ry = powf(r, y);
114                         ryinv = (r > 0.0f)? powf(r, 2.0f*y - 2.0f): 0.0f;
115                 }
117                 const float Rmy5 = (Rmy*Rmy) * (Rmy*Rmy) * Rmy;
118                 const float f = Rmy - ry;
119                 const float num = f*(f*f)*(y*ryinv);
121                 return (10.0f * num) / (Rmy5 * M_PI_F);
122         }
123 }
125 ccl_device float bssrdf_cubic_pdf(ShaderClosure *sc, float r)
126 {
127         return bssrdf_cubic_eval(sc, r);
128 }
130 /* solve 10x^2 - 20x^3 + 15x^4 - 4x^5 - xi == 0 */
131 ccl_device float bssrdf_cubic_quintic_root_find(float xi)
132 {
133         /* newton-raphson iteration, usually succeeds in 2-4 iterations, except
134          * outside 0.02 ... 0.98 where it can go up to 10, so overall performance
135          * should not be too bad */
136         const float tolerance = 1e-6f;
137         const int max_iteration_count = 10;
138         float x = 0.25f;
139         int i;
141         for(i = 0; i < max_iteration_count; i++) {
142                 float x2 = x*x;
143                 float x3 = x2*x;
144                 float nx = (1.0f - x);
146                 float f = 10.0f*x2 - 20.0f*x3 + 15.0f*x2*x2 - 4.0f*x2*x3 - xi;
147                 float f_ = 20.0f*(x*nx)*(nx*nx);
149                 if(fabsf(f) < tolerance || f_ == 0.0f)
150                         break;
152                 x = saturate(x - f/f_);
153         }
155         return x;
156 }
158 ccl_device void bssrdf_cubic_sample(ShaderClosure *sc, float xi, float *r, float *h)
159 {
160         float Rm = sc->data0;
161         float r_ = bssrdf_cubic_quintic_root_find(xi);
163         const float sharpness = sc->T.x;
164         if(sharpness != 0.0f) {
165                 r_ = powf(r_, 1.0f + sharpness);
166                 Rm *= (1.0f + sharpness);
167         }
169         r_ *= Rm;
170         *r = r_;
172         /* h^2 + r^2 = Rm^2 */
173         *h = safe_sqrtf(Rm*Rm - r_*r_);
174 }
176 /* Approximate Reflectance Profiles
177  * http://graphics.pixar.com/library/ApproxBSSRDF/paper.pdf
178  */
180 /* This is a bit arbitrary, just need big enough radius so it matches
181  * the mean free length, but still not too big so sampling is still
182  * effective. Might need some further tweaks.
183  */
184 #define BURLEY_TRUNCATE     16.0f
185 #define BURLEY_TRUNCATE_CDF 0.9963790093708328f // cdf(BURLEY_TRUNCATE)
187 ccl_device_inline float bssrdf_burley_fitting(float A)
188 {
189         /* Diffuse surface transmission, equation (6). */
190         return 1.9f - A + 3.5f * (A - 0.8f) * (A - 0.8f);
191 }
193 /* Scale mean free path length so it gives similar looking result
194  * to Cubic and Gaussian models.
195  */
196 ccl_device_inline float bssrdf_burley_compatible_mfp(float r)
197 {
198         return 0.25f * M_1_PI_F * r;
199 }
202 {
203         /* Mean free path length. */
204         const float l = bssrdf_burley_compatible_mfp(sc->data0);
205         /* Surface albedo. */
206         const float A = sc->data2;
207         const float s = bssrdf_burley_fitting(A);
208         const float d = l / s;
210         sc->custom1 = d;
211 }
213 ccl_device float bssrdf_burley_eval(ShaderClosure *sc, float r)
214 {
215         const float d = sc->custom1;
216         const float Rm = BURLEY_TRUNCATE * d;
218         if(r >= Rm)
219                 return 0.0f;
221         /* Burley refletance profile, equation (3).
222          *
223          * NOTES:
224          * - Surface albedo is already included into sc->weight, no need to
225          *   multiply by this term here.
226          * - This is normalized diffuse model, so the equation is mutliplied
227          *   by 2*pi, which also matches cdf().
228          */
229         float exp_r_3_d = expf(-r / (3.0f * d));
230         float exp_r_d = exp_r_3_d * exp_r_3_d * exp_r_3_d;
231         return (exp_r_d + exp_r_3_d) / (4.0f*d);
232 }
234 ccl_device float bssrdf_burley_pdf(ShaderClosure *sc, float r)
235 {
236         return bssrdf_burley_eval(sc, r) * (1.0f/BURLEY_TRUNCATE_CDF);
237 }
239 /* Find the radius for desired CDF value.
240  * Returns scaled radius, meaning the result is to be scaled up by d.
241  * Since there's no closed form solution we do Newton-Raphson method to find it.
242  */
243 ccl_device float bssrdf_burley_root_find(float xi)
244 {
245         const float tolerance = 1e-6f;
246         const int max_iteration_count = 10;
247         /* Do initial guess based on manual curve fitting, this allows us to reduce
248          * number of iterations to maximum 4 across the [0..1] range. We keep maximum
249          * number of iteration higher just to be sure we didn't miss root in some
250          * corner case.
251          */
252         float r;
253         if (xi <= 0.9f) {
254                 r = expf(xi * xi * 2.4f) - 1.0f;
255         }
256         else {
257                 /* TODO(sergey): Some nicer curve fit is possible here. */
258                 r = 15.0f;
259         }
260         /* Solve against scaled radius. */
261         for(int i = 0; i < max_iteration_count; i++) {
262                 float exp_r_3 = expf(-r / 3.0f);
263                 float exp_r = exp_r_3 * exp_r_3 * exp_r_3;
264                 float f = 1.0f - 0.25f * exp_r - 0.75f * exp_r_3 - xi;
265                 float f_ = 0.25f * exp_r + 0.25f * exp_r_3;
267                 if(fabsf(f) < tolerance || f_ == 0.0f) {
268                         break;
269                 }
271                 r = r - f/f_;
272                 if(r < 0.0f) {
273                         r = 0.0f;
274                 }
275         }
276         return r;
277 }
280                                      float xi,
281                                      float *r,
282                                      float *h)
283 {
284         const float d = sc->custom1;
285         const float Rm = BURLEY_TRUNCATE * d;
286         const float r_ = bssrdf_burley_root_find(xi * BURLEY_TRUNCATE_CDF) * d;
288         *r = r_;
290         /* h^2 + r^2 = Rm^2 */
291         *h = safe_sqrtf(Rm*Rm - r_*r_);
292 }
294 /* None BSSRDF falloff
295  *
296  * Samples distributed over disk with no falloff, for reference. */
298 ccl_device float bssrdf_none_eval(ShaderClosure *sc, float r)
299 {
300         const float Rm = sc->data0;
301         return (r < Rm)? 1.0f: 0.0f;
302 }
304 ccl_device float bssrdf_none_pdf(ShaderClosure *sc, float r)
305 {
306         /* integrate (2*pi*r)/(pi*Rm*Rm) from 0 to Rm = 1 */
307         const float Rm = sc->data0;
308         const float area = (M_PI_F*Rm*Rm);
310         return bssrdf_none_eval(sc, r) / area;
311 }
313 ccl_device void bssrdf_none_sample(ShaderClosure *sc, float xi, float *r, float *h)
314 {
315         /* xi = integrate (2*pi*r)/(pi*Rm*Rm) = r^2/Rm^2
316          * r = sqrt(xi)*Rm */
317         const float Rm = sc->data0;
318         const float r_ = sqrtf(xi)*Rm;
320         *r = r_;
322         /* h^2 + r^2 = Rm^2 */
323         *h = safe_sqrtf(Rm*Rm - r_*r_);
324 }
326 /* Generic */
328 ccl_device int bssrdf_setup(ShaderClosure *sc, ClosureType type)
329 {
331                 /* revert to diffuse BSDF if radius too small */
332                 sc->data0 = 0.0f;
333                 sc->data1 = 0.0f;
334                 int flag = bsdf_diffuse_setup(sc);
335                 sc->type = CLOSURE_BSDF_BSSRDF_ID;
336                 return flag;
337         }
338         else {
339                 sc->data1 = saturate(sc->data1); /* texture blur */
340                 sc->T.x = saturate(sc->T.x); /* sharpness */
341                 sc->type = type;
343                 if(type == CLOSURE_BSSRDF_BURLEY_ID) {
344                         bssrdf_burley_setup(sc);
345                 }
347                 return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSSRDF;
348         }
349 }
351 ccl_device void bssrdf_sample(ShaderClosure *sc, float xi, float *r, float *h)
352 {
353         if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
354                 bssrdf_cubic_sample(sc, xi, r, h);
355         else if(sc->type == CLOSURE_BSSRDF_GAUSSIAN_ID)
356                 bssrdf_gaussian_sample(sc, xi, r, h);
357         else /*if(sc->type == CLOSURE_BSSRDF_BURLEY_ID)*/
358                 bssrdf_burley_sample(sc, xi, r, h);
359 }
361 ccl_device float bssrdf_pdf(ShaderClosure *sc, float r)
362 {
363         if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
364                 return bssrdf_cubic_pdf(sc, r);
365         else if(sc->type == CLOSURE_BSSRDF_GAUSSIAN_ID)
366                 return bssrdf_gaussian_pdf(sc, r);
367         else /*if(sc->type == CLOSURE_BSSRDF_BURLEY_ID)*/
368                 return bssrdf_burley_pdf(sc, r);
369 }
371 CCL_NAMESPACE_END
373 #endif /* __KERNEL_BSSRDF_H__ */