Fix assert rendering hair tests on some systems.
[blender.git] / intern / cycles / kernel / closure / bsdf_hair_principled.h
1 /*
2  * Copyright 2018 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 #ifdef __KERNEL_CPU__
18 #include <fenv.h>
19 #endif
20
21 #include "kernel/kernel_color.h"
22
23 #ifndef __BSDF_HAIR_PRINCIPLED_H__
24 #define __BSDF_HAIR_PRINCIPLED_H__
25
26 CCL_NAMESPACE_BEGIN
27
28 typedef ccl_addr_space struct PrincipledHairExtra {
29         /* Geometry data. */
30         float4 geom;
31 } PrincipledHairExtra;
32
33 typedef ccl_addr_space struct PrincipledHairBSDF {
34         SHADER_CLOSURE_BASE;
35
36         /* Absorption coefficient. */
37         float3 sigma;
38         /* Variance of the underlying logistic distribution. */
39         float v;
40         /* Scale factor of the underlying logistic distribution. */
41         float s;
42         /* Cuticle tilt angle. */
43         float alpha;
44         /* IOR. */
45         float eta;
46         /* Effective variance for the diffuse bounce only. */
47         float m0_roughness;
48
49         /* Extra closure. */
50         PrincipledHairExtra *extra;
51 } PrincipledHairBSDF;
52
53 static_assert(sizeof(ShaderClosure) >= sizeof(PrincipledHairBSDF), "PrincipledHairBSDF is too large!");
54 static_assert(sizeof(ShaderClosure) >= sizeof(PrincipledHairExtra), "PrincipledHairExtra is too large!");
55
56 ccl_device_inline float cos_from_sin(const float s)
57 {
58         return safe_sqrtf(1.0f - s*s);
59 }
60
61 /* Gives the change in direction in the normal plane for the given angles and p-th-order scattering. */
62 ccl_device_inline float delta_phi(int p, float gamma_o, float gamma_t)
63 {
64         return 2.0f * p * gamma_t - 2.0f * gamma_o + p * M_PI_F;
65 }
66
67 /* Remaps the given angle to [-pi, pi]. */
68 ccl_device_inline float wrap_angle(float a)
69 {
70         while(a > M_PI_F) {
71                 a -= M_2PI_F;
72         }
73         while(a < -M_PI_F) {
74                 a += M_2PI_F;
75         }
76         return a;
77 }
78
79 /* Logistic distribution function. */
80 ccl_device_inline float logistic(float x, float s)
81 {
82         float v = expf(-fabsf(x)/s);
83         return v / (s * sqr(1.0f + v));
84 }
85
86 /* Logistic cumulative density function. */
87 ccl_device_inline float logistic_cdf(float x, float s)
88 {
89         float arg = -x/s;
90         /* expf() overflows if arg >= 89.0. */
91         if(arg > 88.0f) {
92                 return 0.0f;
93         }
94         else {
95                 return 1.0f / (1.0f + expf(arg));
96         }
97 }
98
99 /* Numerical approximation to the Bessel function of the first kind. */
100 ccl_device_inline float bessel_I0(float x)
101 {
102         x = sqr(x);
103         float val = 1.0f + 0.25f*x;
104         float pow_x_2i = sqr(x);
105         uint64_t i_fac_2 = 1;
106         int pow_4_i = 16;
107         for(int i = 2; i < 10; i++) {
108                 i_fac_2 *= i*i;
109                 float newval = val + pow_x_2i / (pow_4_i * i_fac_2);
110                 if(val == newval) {
111                         return val;
112                 }
113                 val = newval;
114                 pow_x_2i *= x;
115                 pow_4_i *= 4;
116         }
117         return val;
118 }
119
120 /* Logarithm of the Bessel function of the first kind. */
121 ccl_device_inline float log_bessel_I0(float x)
122 {
123         if(x > 12.0f) {
124                 /* log(1/x) == -log(x) iff x > 0.
125                  * This is only used with positive cosines */
126                 return x + 0.5f * (1.f / (8.0f * x) - M_LN_2PI_F - logf(x));
127         }
128         else {
129                 return logf(bessel_I0(x));
130         }
131 }
132
133 /* Logistic distribution limited to the interval [-pi, pi]. */
134 ccl_device_inline float trimmed_logistic(float x, float s)
135 {
136         /* The logistic distribution is symmetric and centered around zero,
137          * so logistic_cdf(x, s) = 1 - logistic_cdf(-x, s).
138          * Therefore, logistic_cdf(x, s)-logistic_cdf(-x, s) = 1 - 2*logistic_cdf(-x, s) */
139         float scaling_fac = 1.0f - 2.0f*logistic_cdf(-M_PI_F, s);
140         float val = logistic(x, s);
141         return safe_divide(val, scaling_fac);
142 }
143
144 /* Sampling function for the trimmed logistic function. */
145 ccl_device_inline float sample_trimmed_logistic(float u, float s)
146 {
147         float cdf_minuspi = logistic_cdf(-M_PI_F, s);
148         float x = -s*logf(1.0f / (u*(1.0f - 2.0f*cdf_minuspi) + cdf_minuspi) - 1.0f);
149         return clamp(x, -M_PI_F, M_PI_F);
150 }
151
152 /* Azimuthal scattering function Np. */
153 ccl_device_inline float azimuthal_scattering(float phi,
154                                              int p,
155                                              float s,
156                                              float gamma_o,
157                                              float gamma_t)
158 {
159         float phi_o = wrap_angle(phi - delta_phi(p, gamma_o, gamma_t));
160         float val = trimmed_logistic(phi_o, s);
161         return val;
162 }
163
164 /* Longitudinal scattering function Mp. */
165 ccl_device_inline float longitudinal_scattering(float sin_theta_i,
166                                                 float cos_theta_i,
167                                                 float sin_theta_o,
168                                                 float cos_theta_o,
169                                                 float v)
170 {
171         float inv_v = 1.0f/v;
172         float cos_arg = cos_theta_i * cos_theta_o * inv_v;
173         float sin_arg = sin_theta_i * sin_theta_o * inv_v;
174         if(v <= 0.1f) {
175                 float i0 = log_bessel_I0(cos_arg);
176                 float val = expf(i0 - sin_arg - inv_v + 0.6931f + logf(0.5f*inv_v));
177                 return val;
178         }
179         else {
180                 float i0 = bessel_I0(cos_arg);
181                 float val = (expf(-sin_arg) * i0) / (sinhf(inv_v) * 2.0f * v);
182                 return val;
183         }
184 }
185
186 /* Combine the three values using their luminances. */
187 ccl_device_inline float4 combine_with_energy(KernelGlobals *kg, float3 c)
188 {
189         return make_float4(c.x, c.y, c.z, linear_rgb_to_gray(kg, c));
190 }
191
192 #ifdef __HAIR__
193 /* Set up the hair closure. */
194 ccl_device int bsdf_principled_hair_setup(ShaderData *sd, PrincipledHairBSDF *bsdf)
195 {
196         bsdf->type = CLOSURE_BSDF_HAIR_PRINCIPLED_ID;
197         bsdf->v = clamp(bsdf->v, 0.001f, 1.0f);
198         bsdf->s = clamp(bsdf->s, 0.001f, 1.0f);
199         /* Apply Primary Reflection Roughness modifier. */
200         bsdf->m0_roughness = clamp(bsdf->m0_roughness*bsdf->v, 0.001f, 1.0f);
201
202         /* Map from roughness_u and roughness_v to variance and scale factor. */
203         bsdf->v = sqr(0.726f*bsdf->v + 0.812f*sqr(bsdf->v) + 3.700f*pow20(bsdf->v));
204         bsdf->s =    (0.265f*bsdf->s + 1.194f*sqr(bsdf->s) + 5.372f*pow22(bsdf->s))*M_SQRT_PI_8_F;
205         bsdf->m0_roughness = sqr(0.726f*bsdf->m0_roughness + 0.812f*sqr(bsdf->m0_roughness) + 3.700f*pow20(bsdf->m0_roughness));
206
207         /* Compute local frame, aligned to curve tangent and ray direction. */
208         float3 X = safe_normalize(sd->dPdu);
209         float3 Y = safe_normalize(cross(X, sd->I));
210         float3 Z = safe_normalize(cross(X, Y));
211         /* TODO: the solution below works where sd->Ng is the normal
212          * pointing from the center of the curve to the shading point.
213          * It doesn't work for triangles, see https://developer.blender.org/T43625 */
214
215         /* h -1..0..1 means the rays goes from grazing the hair, to hitting it at
216          * the center, to grazing the other edge. This is the sine of the angle
217          * between sd->Ng and Z, as seen from the tangent X. */
218
219         /* TODO: we convert this value to a cosine later and discard the sign, so
220          * we could probably save some operations. */
221         float h = dot(cross(sd->Ng, X), Z);
222
223         kernel_assert(fabsf(h) < 1.0f + 1e-4f);
224         kernel_assert(isfinite3_safe(Y));
225         kernel_assert(isfinite_safe(h));
226
227         bsdf->extra->geom = make_float4(Y.x, Y.y, Y.z, h);
228
229         return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSDF_NEEDS_LCG;
230 }
231
232 #endif /* __HAIR__ */
233
234 /* Given the Fresnel term and transmittance, generate the attenuation terms for each bounce. */
235 ccl_device_inline void hair_attenuation(KernelGlobals *kg,
236                                         float f,
237                                         float3 T,
238                                         float4 *Ap)
239 {
240         /* Primary specular (R). */
241         Ap[0] = make_float4(f, f, f, f);
242
243         /* Transmission (TT). */
244         float3 col = sqr(1.0f - f) * T;
245         Ap[1] = combine_with_energy(kg, col);
246
247         /* Secondary specular (TRT). */
248         col *= T*f;
249         Ap[2] = combine_with_energy(kg, col);
250
251         /* Residual component (TRRT+). */
252         col *= safe_divide_color(T*f, make_float3(1.0f, 1.0f, 1.0f) - T*f);
253         Ap[3] = combine_with_energy(kg, col);
254
255         /* Normalize sampling weights. */
256         float totweight = Ap[0].w + Ap[1].w + Ap[2].w + Ap[3].w;
257         float fac = safe_divide(1.0f, totweight);
258
259         Ap[0].w *= fac;
260         Ap[1].w *= fac;
261         Ap[2].w *= fac;
262         Ap[3].w *= fac;
263 }
264
265 /* Given the tilt angle, generate the rotated theta_i for the different bounces. */
266 ccl_device_inline void hair_alpha_angles(float sin_theta_i,
267                                          float cos_theta_i,
268                                          float alpha,
269                                          float *angles)
270 {
271         float sin_1alpha = sinf(alpha);
272         float cos_1alpha = cos_from_sin(sin_1alpha);
273         float sin_2alpha = 2.0f*sin_1alpha*cos_1alpha;
274         float cos_2alpha = sqr(cos_1alpha) - sqr(sin_1alpha);
275         float sin_4alpha = 2.0f*sin_2alpha*cos_2alpha;
276         float cos_4alpha = sqr(cos_2alpha) - sqr(sin_2alpha);
277
278         angles[0] = sin_theta_i*cos_2alpha + cos_theta_i*sin_2alpha;
279         angles[1] = fabsf(cos_theta_i*cos_2alpha - sin_theta_i*sin_2alpha);
280         angles[2] = sin_theta_i*cos_1alpha - cos_theta_i*sin_1alpha;
281         angles[3] = fabsf(cos_theta_i*cos_1alpha + sin_theta_i*sin_1alpha);
282         angles[4] = sin_theta_i*cos_4alpha - cos_theta_i*sin_4alpha;
283         angles[5] = fabsf(cos_theta_i*cos_4alpha + sin_theta_i*sin_4alpha);
284 }
285
286 /* Evaluation function for our shader. */
287 ccl_device float3 bsdf_principled_hair_eval(KernelGlobals *kg,
288                                             const ShaderData *sd,
289                                             const ShaderClosure *sc,
290                                             const float3 omega_in,
291                                             float *pdf)
292 {
293         kernel_assert(isfinite3_safe(sd->P) && isfinite_safe(sd->ray_length));
294
295         const PrincipledHairBSDF *bsdf = (const PrincipledHairBSDF*) sc;
296         float3 Y = float4_to_float3(bsdf->extra->geom);
297
298         float3 X = safe_normalize(sd->dPdu);
299         kernel_assert(fabsf(dot(X, Y)) < 1e-3f);
300         float3 Z = safe_normalize(cross(X, Y));
301
302         float3 wo = make_float3(dot(sd->I, X), dot(sd->I, Y), dot(sd->I, Z));
303         float3 wi = make_float3(dot(omega_in, X), dot(omega_in, Y), dot(omega_in, Z));
304
305         float sin_theta_o = wo.x;
306         float cos_theta_o = cos_from_sin(sin_theta_o);
307         float phi_o = atan2f(wo.z, wo.y);
308
309         float sin_theta_t = sin_theta_o / bsdf->eta;
310         float cos_theta_t = cos_from_sin(sin_theta_t);
311
312         float sin_gamma_o = bsdf->extra->geom.w;
313         float cos_gamma_o = cos_from_sin(sin_gamma_o);
314         float gamma_o = safe_asinf(sin_gamma_o);
315
316         float sin_gamma_t = sin_gamma_o * cos_theta_o / sqrtf(sqr(bsdf->eta) - sqr(sin_theta_o));
317         float cos_gamma_t = cos_from_sin(sin_gamma_t);
318         float gamma_t = safe_asinf(sin_gamma_t);
319
320         float3 T = exp3(-bsdf->sigma * (2.0f * cos_gamma_t / cos_theta_t));
321         float4 Ap[4];
322         hair_attenuation(kg, fresnel_dielectric_cos(cos_theta_o * cos_gamma_o, bsdf->eta), T, Ap);
323
324         float sin_theta_i = wi.x;
325         float cos_theta_i = cos_from_sin(sin_theta_i);
326         float phi_i = atan2f(wi.z, wi.y);
327
328         float phi = phi_i - phi_o;
329
330         float angles[6];
331         hair_alpha_angles(sin_theta_i, cos_theta_i, bsdf->alpha, angles);
332
333         float4 F;
334         float Mp, Np;
335
336         /* Primary specular (R). */
337         Mp = longitudinal_scattering(angles[0], angles[1], sin_theta_o, cos_theta_o, bsdf->m0_roughness);
338         Np = azimuthal_scattering(phi, 0, bsdf->s, gamma_o, gamma_t);
339         F  = Ap[0] * Mp * Np;
340         kernel_assert(isfinite3_safe(float4_to_float3(F)));
341
342         /* Transmission (TT). */
343         Mp = longitudinal_scattering(angles[2], angles[3], sin_theta_o, cos_theta_o, 0.25f*bsdf->v);
344         Np = azimuthal_scattering(phi, 1, bsdf->s, gamma_o, gamma_t);
345         F += Ap[1] * Mp * Np;
346         kernel_assert(isfinite3_safe(float4_to_float3(F)));
347
348         /* Secondary specular (TRT). */
349         Mp = longitudinal_scattering(angles[4], angles[5], sin_theta_o, cos_theta_o, 4.0f*bsdf->v);
350         Np = azimuthal_scattering(phi, 2, bsdf->s, gamma_o, gamma_t);
351         F += Ap[2] * Mp * Np;
352         kernel_assert(isfinite3_safe(float4_to_float3(F)));
353
354         /* Residual component (TRRT+). */
355         Mp = longitudinal_scattering(sin_theta_i, cos_theta_i, sin_theta_o, cos_theta_o, 4.0f*bsdf->v);
356         Np = M_1_2PI_F;
357         F += Ap[3] * Mp * Np;
358         kernel_assert(isfinite3_safe(float4_to_float3(F)));
359
360         *pdf = F.w;
361         return float4_to_float3(F);
362 }
363
364 /* Sampling function for the hair shader. */
365 ccl_device int bsdf_principled_hair_sample(KernelGlobals *kg,
366                                            const ShaderClosure *sc,
367                                            ShaderData *sd,
368                                            float randu,
369                                            float randv,
370                                            float3 *eval,
371                                            float3 *omega_in,
372                                            float3 *domega_in_dx,
373                                            float3 *domega_in_dy,
374                                            float *pdf)
375 {
376         PrincipledHairBSDF *bsdf = (PrincipledHairBSDF*) sc;
377
378         float3 Y = float4_to_float3(bsdf->extra->geom);
379
380         float3 X = safe_normalize(sd->dPdu);
381         kernel_assert(fabsf(dot(X, Y)) < 1e-3f);
382         float3 Z = safe_normalize(cross(X, Y));
383
384         float3 wo = make_float3(dot(sd->I, X), dot(sd->I, Y), dot(sd->I, Z));
385
386         float2 u[2];
387         u[0] = make_float2(randu, randv);
388         u[1].x = lcg_step_float_addrspace(&sd->lcg_state);
389         u[1].y = lcg_step_float_addrspace(&sd->lcg_state);
390
391         float sin_theta_o = wo.x;
392         float cos_theta_o = cos_from_sin(sin_theta_o);
393         float phi_o = atan2f(wo.z, wo.y);
394
395         float sin_theta_t = sin_theta_o / bsdf->eta;
396         float cos_theta_t = cos_from_sin(sin_theta_t);
397
398         float sin_gamma_o = bsdf->extra->geom.w;
399         float cos_gamma_o = cos_from_sin(sin_gamma_o);
400         float gamma_o = safe_asinf(sin_gamma_o);
401
402         float sin_gamma_t = sin_gamma_o * cos_theta_o / sqrtf(sqr(bsdf->eta) - sqr(sin_theta_o));
403         float cos_gamma_t = cos_from_sin(sin_gamma_t);
404         float gamma_t = safe_asinf(sin_gamma_t);
405
406         float3 T = exp3(-bsdf->sigma * (2.0f * cos_gamma_t / cos_theta_t));
407         float4 Ap[4];
408         hair_attenuation(kg, fresnel_dielectric_cos(cos_theta_o * cos_gamma_o, bsdf->eta), T, Ap);
409
410         int p = 0;
411         for(; p < 3; p++) {
412                 if(u[0].x < Ap[p].w) {
413                         break;
414                 }
415                 u[0].x -= Ap[p].w;
416         }
417
418         float v = bsdf->v;
419         if(p == 1) {
420                 v *= 0.25f;
421         }
422         if(p >= 2) {
423                 v *= 4.0f;
424         }
425
426         u[1].x = max(u[1].x, 1e-5f);
427         float fac = 1.0f + v*logf(u[1].x + (1.0f - u[1].x)*expf(-2.0f/v));
428         float sin_theta_i = -fac * sin_theta_o + cos_from_sin(fac) * cosf(M_2PI_F * u[1].y) * cos_theta_o;
429         float cos_theta_i = cos_from_sin(sin_theta_i);
430
431         float angles[6];
432         if(p < 3) {
433                 hair_alpha_angles(sin_theta_i, cos_theta_i, -bsdf->alpha, angles);
434                 sin_theta_i = angles[2*p];
435                 cos_theta_i = angles[2*p+1];
436         }
437
438         float phi;
439         if(p < 3) {
440                 phi = delta_phi(p, gamma_o, gamma_t) + sample_trimmed_logistic(u[0].y, bsdf->s);
441         }
442         else {
443                 phi = M_2PI_F*u[0].y;
444         }
445         float phi_i = phi_o + phi;
446
447         hair_alpha_angles(sin_theta_i, cos_theta_i, bsdf->alpha, angles);
448
449         float4 F;
450         float Mp, Np;
451
452         /* Primary specular (R). */
453         Mp = longitudinal_scattering(angles[0], angles[1], sin_theta_o, cos_theta_o, bsdf->m0_roughness);
454         Np = azimuthal_scattering(phi, 0, bsdf->s, gamma_o, gamma_t);
455         F  = Ap[0] * Mp * Np;
456         kernel_assert(isfinite3_safe(float4_to_float3(F)));
457
458         /* Transmission (TT). */
459         Mp = longitudinal_scattering(angles[2], angles[3], sin_theta_o, cos_theta_o, 0.25f*bsdf->v);
460         Np = azimuthal_scattering(phi, 1, bsdf->s, gamma_o, gamma_t);
461         F += Ap[1] * Mp * Np;
462         kernel_assert(isfinite3_safe(float4_to_float3(F)));
463
464         /* Secondary specular (TRT). */
465         Mp = longitudinal_scattering(angles[4], angles[5], sin_theta_o, cos_theta_o, 4.0f*bsdf->v);
466         Np = azimuthal_scattering(phi, 2, bsdf->s, gamma_o, gamma_t);
467         F += Ap[2] * Mp * Np;
468         kernel_assert(isfinite3_safe(float4_to_float3(F)));
469
470         /* Residual component (TRRT+). */
471         Mp = longitudinal_scattering(sin_theta_i, cos_theta_i, sin_theta_o, cos_theta_o, 4.0f*bsdf->v);
472         Np = M_1_2PI_F;
473         F += Ap[3] * Mp * Np;
474         kernel_assert(isfinite3_safe(float4_to_float3(F)));
475
476         *eval = float4_to_float3(F);
477         *pdf = F.w;
478
479         *omega_in = X*sin_theta_i + Y*cos_theta_i*cosf(phi_i) + Z*cos_theta_i*sinf(phi_i);
480
481 #ifdef __RAY_DIFFERENTIALS__
482         float3 N = safe_normalize(sd->I + *omega_in);
483         *domega_in_dx = (2 * dot(N, sd->dI.dx)) * N - sd->dI.dx;
484         *domega_in_dy = (2 * dot(N, sd->dI.dy)) * N - sd->dI.dy;
485 #endif
486
487         return LABEL_GLOSSY|((p == 0)? LABEL_REFLECT : LABEL_TRANSMIT);
488 }
489
490 /* Implements Filter Glossy by capping the effective roughness. */
491 ccl_device void bsdf_principled_hair_blur(ShaderClosure *sc, float roughness)
492 {
493         PrincipledHairBSDF *bsdf = (PrincipledHairBSDF*)sc;
494
495         bsdf->v = fmaxf(roughness, bsdf->v);
496         bsdf->s = fmaxf(roughness, bsdf->s);
497         bsdf->m0_roughness = fmaxf(roughness, bsdf->m0_roughness);
498 }
499
500 CCL_NAMESPACE_END
501
502 #endif /* __BSDF_HAIR_PRINCIPLED_H__ */