1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
/*
* Copyright 2013, Blender Foundation.
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software Foundation,
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#ifndef __KERNEL_BSSRDF_H__
#define __KERNEL_BSSRDF_H__
CCL_NAMESPACE_BEGIN
__device int bssrdf_setup(ShaderClosure *sc)
{
if(sc->data0 < BSSRDF_MIN_RADIUS) {
/* revert to diffuse BSDF if radius too small */
sc->data0 = 0.0f;
sc->data1 = 0.0f;
return bsdf_diffuse_setup(sc);
}
else {
/* IOR param */
sc->data1 = max(sc->data1, 1.0f);
sc->type = CLOSURE_BSSRDF_ID;
return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSSRDF;
}
}
/* Simple Cubic BSSRDF falloff */
__device float bssrdf_cubic(float ld, float r)
{
if(ld == 0.0f)
return (r == 0.0f)? 1.0f: 0.0f;
return powf(ld - min(r, ld), 3.0f) * 4.0f/powf(ld, 4.0f);
}
/* Original BSSRDF fallof function */
typedef struct BSSRDFParams {
float eta; /* index of refraction */
float sigma_t_; /* reduced extinction coefficient */
float sigma_tr; /* effective extinction coefficient */
float Fdr; /* diffuse fresnel reflectance */
float D; /* diffusion constant */
float A;
float alpha_; /* reduced albedo */
float zr; /* distance of virtual lightsource above surface */
float zv; /* distance of virtual lightsource below surface */
float ld; /* mean free path */
float ro; /* diffuse reflectance */
} BSSRDFParams;
__device float bssrdf_reduced_albedo_Rd(float alpha_, float A, float ro)
{
float sq;
sq = sqrtf(3.0f*(1.0f - alpha_));
return (alpha_/2.0f)*(1.0f + expf((-4.0f/3.0f)*A*sq))*expf(-sq) - ro;
}
__device float bssrdf_compute_reduced_albedo(float A, float ro)
{
const float tolerance = 1e-8f;
const int max_iteration_count = 20;
float d, fsub, xn_1 = 0.0f, xn = 1.0f, fxn, fxn_1;
int i;
/* use secant method to compute reduced albedo using Rd function inverse
* with a given reflectance */
fxn = bssrdf_reduced_albedo_Rd(xn, A, ro);
fxn_1 = bssrdf_reduced_albedo_Rd(xn_1, A, ro);
for (i= 0; i < max_iteration_count; i++) {
fsub = (fxn - fxn_1);
if (fabsf(fsub) < tolerance)
break;
d = ((xn - xn_1)/fsub)*fxn;
if (fabsf(d) < tolerance)
break;
xn_1 = xn;
fxn_1 = fxn;
xn = xn - d;
if (xn > 1.0f) xn = 1.0f;
if (xn_1 > 1.0f) xn_1 = 1.0f;
fxn = bssrdf_reduced_albedo_Rd(xn, A, ro);
}
/* avoid division by zero later */
if (xn <= 0.0f)
xn = 0.00001f;
return xn;
}
__device void bssrdf_setup_params(BSSRDFParams *ss, float refl, float radius, float ior)
{
ss->eta = ior;
ss->Fdr = -1.440f/ior*ior + 0.710f/ior + 0.668f + 0.0636f*ior;
ss->A = (1.0f + ss->Fdr)/(1.0f - ss->Fdr);
ss->ld = radius;
ss->ro = min(refl, 0.999f);
ss->alpha_ = bssrdf_compute_reduced_albedo(ss->A, ss->ro);
ss->sigma_tr = 1.0f/ss->ld;
ss->sigma_t_ = ss->sigma_tr/sqrtf(3.0f*(1.0f - ss->alpha_));
ss->D = 1.0f/(3.0f*ss->sigma_t_);
ss->zr = 1.0f/ss->sigma_t_;
ss->zv = ss->zr + 4.0f*ss->A*ss->D;
}
/* exponential falloff function */
__device float bssrdf_original(const BSSRDFParams *ss, float r)
{
if(ss->ld == 0.0f)
return (r == 0.0f)? 1.0f: 0.0f;
float rr = r*r;
float sr, sv, Rdr, Rdv;
sr = sqrtf(rr + ss->zr*ss->zr);
sv = sqrtf(rr + ss->zv*ss->zv);
Rdr = ss->zr*(1.0f + ss->sigma_tr*sr)*expf(-ss->sigma_tr*sr)/(sr*sr*sr);
Rdv = ss->zv*(1.0f + ss->sigma_tr*sv)*expf(-ss->sigma_tr*sv)/(sv*sv*sv);
return ss->alpha_*(1.0f/M_4PI_F)*(Rdr + Rdv);
}
CCL_NAMESPACE_END
#endif /* __KERNEL_BSSRDF_H__ */
|
