Cycles: change __device and similar qualifiers to ccl_device in kernel code.
[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 ccl_device int bssrdf_setup(ShaderClosure *sc, ClosureType type)
23 {
24         if(sc->data0 < BSSRDF_MIN_RADIUS) {
25                 /* revert to diffuse BSDF if radius too small */
26                 sc->data0 = 0.0f;
27                 sc->data1 = 0.0f;
28                 int flag = bsdf_diffuse_setup(sc);
29                 sc->type = CLOSURE_BSDF_BSSRDF_ID;
30                 return flag;
31         }
32         else {
33                 sc->data1 = clamp(sc->data1, 0.0f, 1.0f); /* texture blur */
34                 sc->T.x = clamp(sc->T.x, 0.0f, 1.0f); /* sharpness */
35                 sc->type = type;
36
37                 return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSSRDF;
38         }
39 }
40
41 /* Planar Truncated Gaussian
42  *
43  * Note how this is different from the typical gaussian, this one integrates
44  * to 1 over the plane (where you get an extra 2*pi*x factor). We are lucky
45  * that integrating x*exp(-x) gives a nice closed form solution. */
46
47 /* paper suggests 1/12.46 which is much too small, suspect it's *12.46 */
48 #define GAUSS_TRUNCATE 12.46f
49
50 ccl_device float bssrdf_gaussian_eval(ShaderClosure *sc, float r)
51 {
52         /* integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) from 0 to Rm
53          * = 1 - exp(-Rm*Rm/(2*v)) */
54         const float v = sc->data0*sc->data0*(0.25f*0.25f);
55         const float Rm = sqrtf(v*GAUSS_TRUNCATE);
56
57         if(r >= Rm)
58                 return 0.0f;
59
60         return expf(-r*r/(2.0f*v))/(2.0f*M_PI_F*v);
61 }
62
63 ccl_device float bssrdf_gaussian_pdf(ShaderClosure *sc, float r)
64 {
65         /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
66         const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
67
68         return bssrdf_gaussian_eval(sc, r) * (1.0f/(area_truncated));
69 }
70
71 ccl_device void bssrdf_gaussian_sample(ShaderClosure *sc, float xi, float *r, float *h)
72 {
73         /* xi = integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) = -exp(-r^2/(2*v))
74          * r = sqrt(-2*v*logf(xi)) */
75
76         const float v = sc->data0*sc->data0*(0.25f*0.25f);
77         const float Rm = sqrtf(v*GAUSS_TRUNCATE);
78
79         /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
80         const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
81
82         /* r(xi) */
83         const float r_squared = -2.0f*v*logf(1.0f - xi*area_truncated);
84         *r = sqrtf(r_squared);
85
86          /* h^2 + r^2 = Rm^2 */
87          *h = sqrtf(Rm*Rm - r_squared);
88 }
89
90 /* Planar Cubic BSSRDF falloff
91  *
92  * This is basically (Rm - x)^3, with some factors to normalize it. For sampling
93  * we integrate 2*pi*x * (Rm - x)^3, which gives us a quintic equation that as
94  * far as I can tell has no closed form solution. So we get an iterative solution
95  * instead with newton-raphson. */
96
97 ccl_device float bssrdf_cubic_eval(ShaderClosure *sc, float r)
98 {
99         const float sharpness = sc->T.x;
100
101         if(sharpness == 0.0f) {
102                 const float Rm = sc->data0;
103
104                 if(r >= Rm)
105                         return 0.0f;
106
107                 /* integrate (2*pi*r * 10*(R - r)^3)/(pi * R^5) from 0 to R = 1 */
108                 const float Rm5 = (Rm*Rm) * (Rm*Rm) * Rm;
109                 const float f = Rm - r;
110                 const float num = f*f*f;
111
112                 return (10.0f * num) / (Rm5 * M_PI_F);
113
114         }
115         else {
116                 float Rm = sc->data0*(1.0f + sharpness);
117
118                 if(r >= Rm)
119                         return 0.0f;
120
121                 /* custom variation with extra sharpness, to match the previous code */
122                 const float y = 1.0f/(1.0f + sharpness);
123                 float Rmy, ry, ryinv;
124
125                 if(sharpness == 1.0f) {
126                         Rmy = sqrtf(Rm);
127                         ry = sqrtf(r);
128                         ryinv = (ry > 0.0f)? 1.0f/ry: 0.0f;
129                 }
130                 else {
131                         Rmy = powf(Rm, y);
132                         ry = powf(r, y);
133                         ryinv = (r > 0.0f)? powf(r, 2.0f*y - 2.0f): 0.0f;
134                 }
135
136                 const float Rmy5 = (Rmy*Rmy) * (Rmy*Rmy) * Rmy;
137                 const float f = Rmy - ry;
138                 const float num = f*(f*f)*(y*ryinv);
139
140                 return (10.0f * num) / (Rmy5 * M_PI_F);
141         }
142 }
143
144 ccl_device float bssrdf_cubic_pdf(ShaderClosure *sc, float r)
145 {
146         return bssrdf_cubic_eval(sc, r);
147 }
148
149 /* solve 10x^2 - 20x^3 + 15x^4 - 4x^5 - xi == 0 */
150 ccl_device float bssrdf_cubic_quintic_root_find(float xi)
151 {
152         /* newton-raphson iteration, usually succeeds in 2-4 iterations, except
153          * outside 0.02 ... 0.98 where it can go up to 10, so overall performance
154          * should not be too bad */
155         const float tolerance = 1e-6f;
156         const int max_iteration_count = 10;
157         float x = 0.25f;
158         int i;
159
160         for (i = 0; i < max_iteration_count; i++) {
161                 float x2 = x*x;
162                 float x3 = x2*x;
163                 float nx = (1.0f - x);
164
165                 float f = 10.0f*x2 - 20.0f*x3 + 15.0f*x2*x2 - 4.0f*x2*x3 - xi;
166                 float f_ = 20.0f*(x*nx)*(nx*nx);
167
168                 if(fabsf(f) < tolerance || f_ == 0.0f)
169                         break;
170
171                 x = clamp(x - f/f_, 0.0f, 1.0f);
172         }
173
174         return x;
175 }
176
177 ccl_device void bssrdf_cubic_sample(ShaderClosure *sc, float xi, float *r, float *h)
178 {
179         float Rm = sc->data0;
180         float r_ = bssrdf_cubic_quintic_root_find(xi);
181
182         const float sharpness = sc->T.x;
183         if(sharpness != 0.0f) {
184                 r_ = powf(r_, 1.0f + sharpness);
185                 Rm *= (1.0f + sharpness);
186         }
187         
188         r_ *= Rm;
189         *r = r_;
190
191         /* h^2 + r^2 = Rm^2 */
192         *h = sqrtf(Rm*Rm - r_*r_);
193 }
194
195 /* None BSSRDF falloff
196  * 
197  * Samples distributed over disk with no falloff, for reference. */
198
199 ccl_device float bssrdf_none_eval(ShaderClosure *sc, float r)
200 {
201         const float Rm = sc->data0;
202         return (r < Rm)? 1.0f: 0.0f;
203 }
204
205 ccl_device float bssrdf_none_pdf(ShaderClosure *sc, float r)
206 {
207         /* integrate (2*pi*r)/(pi*Rm*Rm) from 0 to Rm = 1 */
208         const float Rm = sc->data0;
209         const float area = (M_PI_F*Rm*Rm);
210
211         return bssrdf_none_eval(sc, r) / area;
212 }
213
214 ccl_device void bssrdf_none_sample(ShaderClosure *sc, float xi, float *r, float *h)
215 {
216         /* xi = integrate (2*pi*r)/(pi*Rm*Rm) = r^2/Rm^2
217          * r = sqrt(xi)*Rm */
218         const float Rm = sc->data0;
219         const float r_ = sqrtf(xi)*Rm;
220
221         *r = r_;
222
223         /* h^2 + r^2 = Rm^2 */
224         *h = sqrtf(Rm*Rm - r_*r_);
225 }
226
227 /* Generic */
228
229 ccl_device void bssrdf_sample(ShaderClosure *sc, float xi, float *r, float *h)
230 {
231         if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
232                 bssrdf_cubic_sample(sc, xi, r, h);
233         else
234                 bssrdf_gaussian_sample(sc, xi, r, h);
235 }
236
237 ccl_device float bssrdf_pdf(ShaderClosure *sc, float r)
238 {
239         if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
240                 return bssrdf_cubic_pdf(sc, r);
241         else
242                 return bssrdf_gaussian_pdf(sc, r);
243 }
244
245 CCL_NAMESPACE_END
246
247 #endif /* __KERNEL_BSSRDF_H__ */
248