1 |
{ |
2 |
He-Torrance Reflectance Model (Siggraph 1991) |
3 |
|
4 |
This is the simplified version that doesn't account for |
5 |
changes in reflection due to changes in wavelength. Also, |
6 |
specular and directional-diffuse hightlights are left uncolored |
7 |
because coloring them requires multiple evaluations of some |
8 |
very expensive functions. |
9 |
|
10 |
The primitive for this function should look something like: |
11 |
|
12 |
void BRTDfunc name |
13 |
10 |
14 |
s s s |
15 |
0 0 0 |
16 |
dd dd dd |
17 |
He.cal |
18 |
0 |
19 |
13 amb_r amb_g amb_b |
20 |
amb_r amb_g amb_b |
21 |
0 0 0 |
22 |
sigma0 tau |
23 |
n_real n_imag |
24 |
|
25 |
For metals, the specular color may be modified like so: |
26 |
|
27 |
void BRTDfunc name |
28 |
10 |
29 |
s_r s_g s_b |
30 |
0 0 0 |
31 |
dd dd dd |
32 |
He.cal |
33 |
0 |
34 |
13 amb_r amb_g amb_b |
35 |
amb_r amb_g amb_b |
36 |
0 0 0 |
37 |
sigma0 tau |
38 |
n_real n_imag |
39 |
|
40 |
This doesn't work for the directional diffuse component, unfortunately. |
41 |
A second set of functions dd_r, dd_g and dd_b may be used, but they |
42 |
cost three times as much to compute! |
43 |
} |
44 |
|
45 |
{ Constants } |
46 |
lambda : .5; { wavelength (microns) } |
47 |
z0err : .0001; { accepted error in value of z0 } |
48 |
Dsumlim : .000001; { last term of D summation } |
49 |
Dsummax : 200; { maximum terms in D summation } |
50 |
|
51 |
{ Parameters } |
52 |
sigma0 = arg(10); { surface height deviation (microns) } |
53 |
tau = arg(11); { correlation distance (microns) } |
54 |
n_real = arg(12); { real part of index of refraction } |
55 |
n_imag = arg(13); { imaginary part of index of refraction } |
56 |
{ Derived parameters } |
57 |
n_k = n_imag/n_real; |
58 |
{ Constant functions } |
59 |
Exp(x) : if(-x-400, 0, exp(x)); { rayinit.cal version too timid for D() } |
60 |
|
61 |
{ Repeated formulas } |
62 |
cotexp(t) = tau/sigma0/2/tan(t); |
63 |
shadowf2(et,erfcet) = (1-.5*erfcet) / |
64 |
((Exp(-sq(et))/sqrt(PI)/et - erfcet)/2 + 1); |
65 |
shadowf1(t) = or(FTINY-sigma0, .01-abs(t)); |
66 |
shadowf0(t) = abs(t) - (PI/2-.0001); |
67 |
shadowf(t) = if(shadowf0(t), 0, if(shadowf1(t), 1, |
68 |
shadowf2(cotexp(t), erfc(cotexp(t))))); |
69 |
K(t) = if(abs(t)-FTINY, tan(t) * erfc(cotexp(t)), 0); |
70 |
fuvA(ct) = sq(n_real)*(1-sq(n_k)) - (1-sq(ct)); |
71 |
fuvB(ct) = sqrt(sq(fuvA(ct)) + 4*sq(sq(n_real)*n_k)); |
72 |
fu2(ct) = (fuvA(ct) + fuvB(ct))/2; |
73 |
fv2(ct) = (-fuvA(ct) + fuvB(ct))/2; |
74 |
fperp2(ct) = (sq(ct-sqrt(fu2(ct))) + fv2(ct)) / |
75 |
(sq(ct+sqrt(fu2(ct))) + fv2(ct)); |
76 |
fpara2(ct) = (sq(sq(n_real)*(1-sq(n_k))*ct - sqrt(fu2(ct))) + |
77 |
sq(2*sq(n_real)*n_k*ct - Sqrt(fv2(ct)))) / |
78 |
(sq(sq(n_real)*(1-sq(n_k))*ct + sqrt(fu2(ct))) + |
79 |
sq(2*sq(n_real)*n_k*ct + Sqrt(fv2(ct)))); |
80 |
fresnel2(ct) = (fperp2(ct) + fpara2(ct))/2; |
81 |
|
82 |
{ Formulas dependent only on reflected direction } |
83 |
theta_r = acos(RdotP); |
84 |
shadowf_r = shadowf(theta_r); |
85 |
K_r = K(theta_r); |
86 |
srx = Dy*NzP - Dz*NyP; sry = Dz*NxP - Dx*NzP; srz = Dx*NyP - Dy*NxP; |
87 |
srn2 = sq(srx) + sq(sry) + sq(srz); |
88 |
prx = sry*Dz - srz*Dy; |
89 |
pry = srz*Dx - srx*Dz; |
90 |
prz = srx*Dy - sry*Dx; |
91 |
s = fresnel2(RdotP)*Exp(-g(RdotP))*sq(shadowf_r); |
92 |
s_r = s*arg(1)*CrP; |
93 |
s_g = s*arg(2)*CgP; |
94 |
s_b = s*arg(3)*CbP; |
95 |
|
96 |
{ Formulas dependent on incident direction } |
97 |
{ z0 } |
98 |
z0d(Ki,z) = -(Ki+K_r)/(4*sigma0)*z*Exp(-sq(z/sigma0)/2) - sqrt(PI/2); |
99 |
z0lim(x) = if(x, max(x,z0err), min(x,-z0err)); |
100 |
z0off(Ki,z) = (sigma0/4*(Ki+K_r)*Exp(-sq(z/sigma0)/2)-sqrt(PI/2)*z)/ |
101 |
z0lim(z0d(Ki,z)); |
102 |
z0root(Ki, x0, x1, i) = if(i, |
103 |
if(z0err-abs(x1-x0), |
104 |
x1, |
105 |
z0root(Ki,x1,x1-z0off(Ki,x1),i-1)), |
106 |
0); |
107 |
z0(ti) = z0root(K(ti), .1, -z0off(K(ti),.1), 100); |
108 |
{ sigma } |
109 |
sigma(ti) = if( FTINY-sigma0, sigma0, |
110 |
sigma0/sqrt(1+sq(z0(ti)/sigma0)) ); |
111 |
{ g } |
112 |
g(cti) = sq(2*PI/lambda*sigma(Acos(cti))*(cti+RdotP)); |
113 |
{ |F|^2 } |
114 |
fresnel2dd(kix,kiy,kiz) = fresnel2(sqrt(sq(kix-Dx) + sq(kiy-Dy) + |
115 |
sq(kiz-Dz))/2); |
116 |
{ G } |
117 |
{ The bulk of G was found by Andrew Willmott to cancel. This is the original: |
118 |
G2( kix,kiy,kiz, six,siy,siz ) = |
119 |
sq( (sq(kix-Dx)+sq(kiy-Dy)+sq(kiz-Dz)) / |
120 |
(NxP*(kix-Dx)+NyP*(kiy-Dy)+NzP*(kiz-Dz)) ) / |
121 |
sq(sq(Dy*kiz-Dz*kiy)+sq(Dz*kix-Dx*kiz)+sq(Dx*kiy-Dy*kix)) * |
122 |
(sq(srx*kix+sry*kiy+srz*kiz) + |
123 |
sq(prx*kix+pry*kiy+prz*kiz)) * |
124 |
(sq(six*Dx+siy*Dy+siz*Dz) + |
125 |
sq((siy*kiz-siz*kiy)*Dx+(siz*kix-six*kiz)*Dy+(six*kiy-siy*kix)*Dz)) / |
126 |
srn2 / (sq(six)+sq(siy)+sq(siz)); |
127 |
G(kix,kiy,kiz) = G2(kix,kiy,kiz, |
128 |
kiy*NzP-kiz*NyP, kiz*NxP-kix*NzP, kix*NyP-kiy*NxP); |
129 |
-- Newer version below is much simpler: } |
130 |
G(kix,kiy,kiz) = sq( (sq(kix-Dx)+sq(kiy-Dy)+sq(kiz-Dz)) / |
131 |
(NxP*(kix-Dx)+NyP*(kiy-Dy)+NzP*(kiz-Dz)) ); |
132 |
{ D } |
133 |
Dsum2(m,lt,c,t,e,g) = if(or(m-Dsummax,and(lt-t,Dsumlim-t)),t, |
134 |
t+Dsum2(m+1,t,c*g/(m+1),c*g/(m+1)*Exp(-g-e/(m+1))/(m+1),e,g)); |
135 |
Dsum(e,g) = Dsum2(1,0,g,g*Exp(-g-e),e,g); |
136 |
D(kix,kiy,kiz) = sq(PI)/4/sq(lambda)*sq(tau) * |
137 |
Dsum(sq(2*PI/lambda)/4*sq(tau)* |
138 |
(sq(kix-Dx)+sq(kiy-Dy)+sq(kiz-Dz) - |
139 |
sq(NxP*(kix-Dx)+NyP*(kiy-Dy)+NzP*(kiz-Dz))), |
140 |
g(kix*NxP+kiy*NyP+kiz*NzP)); |
141 |
{ rho_dd } |
142 |
dd2(cti) = shadowf_r*shadowf(Acos(cti))/cti/RdotP; |
143 |
dd(kix,kiy,kiz) = dd2(kix*NxP+kiy*NyP+kiz*NzP)*G(kix,kiy,kiz)* |
144 |
fresnel2dd(kix,kiy,kiz)/PI*D(kix,kiy,kiz); |
145 |
{ Color version 3x as slow! } |
146 |
dd_r(kix,kiy,kiz) = dd(kix,kiy,kiz)*arg(1)*CrP; |
147 |
dd_g(kix,kiy,kiz) = dd(kix,kiy,kiz)*arg(2)*CgP; |
148 |
dd_b(kix,kiy,kiz) = dd(kix,kiy,kiz)*arg(3)*CbP; |