Fix Cycles CUDA performance on CUDA 8.0.
[blender-staging.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         float radius;
26         float sharpness;
27         float d;
28         float texture_blur;
29         float albedo;
30         float3 N;
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 ShaderClosure *sc, 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 Bssrdf *bssrdf = (const Bssrdf*)sc;
47         const float v = bssrdf->radius*bssrdf->radius*(0.25f*0.25f);
48         const float Rm = sqrtf(v*GAUSS_TRUNCATE);
49
50         if(r >= Rm)
51                 return 0.0f;
52
53         return expf(-r*r/(2.0f*v))/(2.0f*M_PI_F*v);
54 }
55
56 ccl_device float bssrdf_gaussian_pdf(const ShaderClosure *sc, float r)
57 {
58         /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
59         const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
60
61         return bssrdf_gaussian_eval(sc, r) * (1.0f/(area_truncated));
62 }
63
64 ccl_device void bssrdf_gaussian_sample(const ShaderClosure *sc, float xi, float *r, float *h)
65 {
66         /* xi = integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) = -exp(-r^2/(2*v))
67          * r = sqrt(-2*v*logf(xi)) */
68         const Bssrdf *bssrdf = (const Bssrdf*)sc;
69         const float v = bssrdf->radius*bssrdf->radius*(0.25f*0.25f);
70         const float Rm = sqrtf(v*GAUSS_TRUNCATE);
71
72         /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
73         const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
74
75         /* r(xi) */
76         const float r_squared = -2.0f*v*logf(1.0f - xi*area_truncated);
77         *r = sqrtf(r_squared);
78
79         /* h^2 + r^2 = Rm^2 */
80         *h = safe_sqrtf(Rm*Rm - r_squared);
81 }
82
83 /* Planar Cubic BSSRDF falloff
84  *
85  * This is basically (Rm - x)^3, with some factors to normalize it. For sampling
86  * we integrate 2*pi*x * (Rm - x)^3, which gives us a quintic equation that as
87  * far as I can tell has no closed form solution. So we get an iterative solution
88  * instead with newton-raphson. */
89
90 ccl_device float bssrdf_cubic_eval(const ShaderClosure *sc, float r)
91 {
92         const Bssrdf *bssrdf = (const Bssrdf*)sc;
93         const float sharpness = bssrdf->sharpness;
94
95         if(sharpness == 0.0f) {
96                 const float Rm = bssrdf->radius;
97
98                 if(r >= Rm)
99                         return 0.0f;
100
101                 /* integrate (2*pi*r * 10*(R - r)^3)/(pi * R^5) from 0 to R = 1 */
102                 const float Rm5 = (Rm*Rm) * (Rm*Rm) * Rm;
103                 const float f = Rm - r;
104                 const float num = f*f*f;
105
106                 return (10.0f * num) / (Rm5 * M_PI_F);
107
108         }
109         else {
110                 float Rm = bssrdf->radius*(1.0f + sharpness);
111
112                 if(r >= Rm)
113                         return 0.0f;
114
115                 /* custom variation with extra sharpness, to match the previous code */
116                 const float y = 1.0f/(1.0f + sharpness);
117                 float Rmy, ry, ryinv;
118
119                 if(sharpness == 1.0f) {
120                         Rmy = sqrtf(Rm);
121                         ry = sqrtf(r);
122                         ryinv = (ry > 0.0f)? 1.0f/ry: 0.0f;
123                 }
124                 else {
125                         Rmy = powf(Rm, y);
126                         ry = powf(r, y);
127                         ryinv = (r > 0.0f)? powf(r, 2.0f*y - 2.0f): 0.0f;
128                 }
129
130                 const float Rmy5 = (Rmy*Rmy) * (Rmy*Rmy) * Rmy;
131                 const float f = Rmy - ry;
132                 const float num = f*(f*f)*(y*ryinv);
133
134                 return (10.0f * num) / (Rmy5 * M_PI_F);
135         }
136 }
137
138 ccl_device float bssrdf_cubic_pdf(const ShaderClosure *sc, float r)
139 {
140         return bssrdf_cubic_eval(sc, r);
141 }
142
143 /* solve 10x^2 - 20x^3 + 15x^4 - 4x^5 - xi == 0 */
144 ccl_device_forceinline float bssrdf_cubic_quintic_root_find(float xi)
145 {
146         /* newton-raphson iteration, usually succeeds in 2-4 iterations, except
147          * outside 0.02 ... 0.98 where it can go up to 10, so overall performance
148          * should not be too bad */
149         const float tolerance = 1e-6f;
150         const int max_iteration_count = 10;
151         float x = 0.25f;
152         int i;
153
154         for(i = 0; i < max_iteration_count; i++) {
155                 float x2 = x*x;
156                 float x3 = x2*x;
157                 float nx = (1.0f - x);
158
159                 float f = 10.0f*x2 - 20.0f*x3 + 15.0f*x2*x2 - 4.0f*x2*x3 - xi;
160                 float f_ = 20.0f*(x*nx)*(nx*nx);
161
162                 if(fabsf(f) < tolerance || f_ == 0.0f)
163                         break;
164
165                 x = saturate(x - f/f_);
166         }
167
168         return x;
169 }
170
171 ccl_device void bssrdf_cubic_sample(const ShaderClosure *sc, float xi, float *r, float *h)
172 {
173         const Bssrdf *bssrdf = (const Bssrdf*)sc;
174         const float sharpness = bssrdf->sharpness;
175         float Rm = bssrdf->radius;
176         float r_ = bssrdf_cubic_quintic_root_find(xi);
177
178         if(sharpness != 0.0f) {
179                 r_ = powf(r_, 1.0f + sharpness);
180                 Rm *= (1.0f + sharpness);
181         }
182
183         r_ *= Rm;
184         *r = r_;
185
186         /* h^2 + r^2 = Rm^2 */
187         *h = safe_sqrtf(Rm*Rm - r_*r_);
188 }
189
190 /* Approximate Reflectance Profiles
191  * http://graphics.pixar.com/library/ApproxBSSRDF/paper.pdf
192  */
193
194 /* This is a bit arbitrary, just need big enough radius so it matches
195  * the mean free length, but still not too big so sampling is still
196  * effective. Might need some further tweaks.
197  */
198 #define BURLEY_TRUNCATE     16.0f
199 #define BURLEY_TRUNCATE_CDF 0.9963790093708328f // cdf(BURLEY_TRUNCATE)
200
201 ccl_device_inline float bssrdf_burley_fitting(float A)
202 {
203         /* Diffuse surface transmission, equation (6). */
204         return 1.9f - A + 3.5f * (A - 0.8f) * (A - 0.8f);
205 }
206
207 /* Scale mean free path length so it gives similar looking result
208  * to Cubic and Gaussian models.
209  */
210 ccl_device_inline float bssrdf_burley_compatible_mfp(float r)
211 {
212         return 0.25f * M_1_PI_F * r;
213 }
214
215 ccl_device void bssrdf_burley_setup(Bssrdf *bssrdf)
216 {
217         /* Mean free path length. */
218         const float l = bssrdf_burley_compatible_mfp(bssrdf->radius);
219         /* Surface albedo. */
220         const float A = bssrdf->albedo;
221         const float s = bssrdf_burley_fitting(A);
222         const float d = l / s;
223
224         bssrdf->d = d;
225 }
226
227 ccl_device float bssrdf_burley_eval(const ShaderClosure *sc, float r)
228 {
229         const Bssrdf *bssrdf = (const Bssrdf*)sc;
230         const float d = bssrdf->d;
231         const float Rm = BURLEY_TRUNCATE * d;
232
233         if(r >= Rm)
234                 return 0.0f;
235
236         /* Burley refletance profile, equation (3).
237          *
238          * NOTES:
239          * - Surface albedo is already included into sc->weight, no need to
240          *   multiply by this term here.
241          * - This is normalized diffuse model, so the equation is mutliplied
242          *   by 2*pi, which also matches cdf().
243          */
244         float exp_r_3_d = expf(-r / (3.0f * d));
245         float exp_r_d = exp_r_3_d * exp_r_3_d * exp_r_3_d;
246         return (exp_r_d + exp_r_3_d) / (4.0f*d);
247 }
248
249 ccl_device float bssrdf_burley_pdf(const ShaderClosure *sc, float r)
250 {
251         return bssrdf_burley_eval(sc, r) * (1.0f/BURLEY_TRUNCATE_CDF);
252 }
253
254 /* Find the radius for desired CDF value.
255  * Returns scaled radius, meaning the result is to be scaled up by d.
256  * Since there's no closed form solution we do Newton-Raphson method to find it.
257  */
258 ccl_device_forceinline float bssrdf_burley_root_find(float xi)
259 {
260         const float tolerance = 1e-6f;
261         const int max_iteration_count = 10;
262         /* Do initial guess based on manual curve fitting, this allows us to reduce
263          * number of iterations to maximum 4 across the [0..1] range. We keep maximum
264          * number of iteration higher just to be sure we didn't miss root in some
265          * corner case.
266          */
267         float r;
268         if(xi <= 0.9f) {
269                 r = expf(xi * xi * 2.4f) - 1.0f;
270         }
271         else {
272                 /* TODO(sergey): Some nicer curve fit is possible here. */
273                 r = 15.0f;
274         }
275         /* Solve against scaled radius. */
276         for(int i = 0; i < max_iteration_count; i++) {
277                 float exp_r_3 = expf(-r / 3.0f);
278                 float exp_r = exp_r_3 * exp_r_3 * exp_r_3;
279                 float f = 1.0f - 0.25f * exp_r - 0.75f * exp_r_3 - xi;
280                 float f_ = 0.25f * exp_r + 0.25f * exp_r_3;
281
282                 if(fabsf(f) < tolerance || f_ == 0.0f) {
283                         break;
284                 }
285
286                 r = r - f/f_;
287                 if(r < 0.0f) {
288                         r = 0.0f;
289                 }
290         }
291         return r;
292 }
293
294 ccl_device void bssrdf_burley_sample(const ShaderClosure *sc,
295                                      float xi,
296                                      float *r,
297                                      float *h)
298 {
299         const Bssrdf *bssrdf = (const Bssrdf*)sc;
300         const float d = bssrdf->d;
301         const float Rm = BURLEY_TRUNCATE * d;
302         const float r_ = bssrdf_burley_root_find(xi * BURLEY_TRUNCATE_CDF) * d;
303
304         *r = r_;
305
306         /* h^2 + r^2 = Rm^2 */
307         *h = safe_sqrtf(Rm*Rm - r_*r_);
308 }
309
310 /* None BSSRDF falloff
311  *
312  * Samples distributed over disk with no falloff, for reference. */
313
314 ccl_device float bssrdf_none_eval(const ShaderClosure *sc, float r)
315 {
316         const Bssrdf *bssrdf = (const Bssrdf*)sc;
317         const float Rm = bssrdf->radius;
318         return (r < Rm)? 1.0f: 0.0f;
319 }
320
321 ccl_device float bssrdf_none_pdf(const ShaderClosure *sc, float r)
322 {
323         /* integrate (2*pi*r)/(pi*Rm*Rm) from 0 to Rm = 1 */
324         const Bssrdf *bssrdf = (const Bssrdf*)sc;
325         const float Rm = bssrdf->radius;
326         const float area = (M_PI_F*Rm*Rm);
327
328         return bssrdf_none_eval(sc, r) / area;
329 }
330
331 ccl_device void bssrdf_none_sample(const ShaderClosure *sc, float xi, float *r, float *h)
332 {
333         /* xi = integrate (2*pi*r)/(pi*Rm*Rm) = r^2/Rm^2
334          * r = sqrt(xi)*Rm */
335         const Bssrdf *bssrdf = (const Bssrdf*)sc;
336         const float Rm = bssrdf->radius;
337         const float r_ = sqrtf(xi)*Rm;
338
339         *r = r_;
340
341         /* h^2 + r^2 = Rm^2 */
342         *h = safe_sqrtf(Rm*Rm - r_*r_);
343 }
344
345 /* Generic */
346
347 ccl_device_inline Bssrdf *bssrdf_alloc(ShaderData *sd, float3 weight)
348 {
349         Bssrdf *bssrdf = (Bssrdf*)closure_alloc(sd, sizeof(Bssrdf), CLOSURE_NONE_ID, weight);
350
351         if(!bssrdf)
352                 return NULL;
353
354         float sample_weight = fabsf(average(weight));
355         bssrdf->sample_weight = sample_weight;
356         return (sample_weight >= CLOSURE_WEIGHT_CUTOFF) ? bssrdf : NULL;
357 }
358
359 ccl_device int bssrdf_setup(Bssrdf *bssrdf, ClosureType type)
360 {
361         if(bssrdf->radius < BSSRDF_MIN_RADIUS) {
362                 /* revert to diffuse BSDF if radius too small */
363                 DiffuseBsdf *bsdf = (DiffuseBsdf*)bssrdf;
364                 bsdf->N = bssrdf->N;
365                 int flag = bsdf_diffuse_setup(bsdf);
366                 bsdf->type = CLOSURE_BSDF_BSSRDF_ID;
367                 return flag;
368         }
369         else {
370                 bssrdf->texture_blur = saturate(bssrdf->texture_blur);
371                 bssrdf->sharpness = saturate(bssrdf->sharpness);
372                 bssrdf->type = type;
373
374                 if(type == CLOSURE_BSSRDF_BURLEY_ID) {
375                         bssrdf_burley_setup(bssrdf);
376                 }
377
378                 return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSSRDF;
379         }
380 }
381
382 ccl_device void bssrdf_sample(const ShaderClosure *sc, float xi, float *r, float *h)
383 {
384         if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
385                 bssrdf_cubic_sample(sc, xi, r, h);
386         else if(sc->type == CLOSURE_BSSRDF_GAUSSIAN_ID)
387                 bssrdf_gaussian_sample(sc, xi, r, h);
388         else /*if(sc->type == CLOSURE_BSSRDF_BURLEY_ID)*/
389                 bssrdf_burley_sample(sc, xi, r, h);
390 }
391
392 ccl_device_forceinline float bssrdf_pdf(const ShaderClosure *sc, float r)
393 {
394         if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
395                 return bssrdf_cubic_pdf(sc, r);
396         else if(sc->type == CLOSURE_BSSRDF_GAUSSIAN_ID)
397                 return bssrdf_gaussian_pdf(sc, r);
398         else /*if(sc->type == CLOSURE_BSSRDF_BURLEY_ID)*/
399                 return bssrdf_burley_pdf(sc, r);
400 }
401
402 CCL_NAMESPACE_END
403
404 #endif /* __KERNEL_BSSRDF_H__ */
405