1 |
{ |
2 |
He-Torrance Reflectance Model (Siggraph 1991) |
3 |
|
4 |
Complete spectral version (8/17/95) |
5 |
|
6 |
The primitive for this function should look something like: |
7 |
|
8 |
void BRTDfunc name |
9 |
10 |
10 |
s`r s`g s`b |
11 |
0 0 0 |
12 |
dd`r dd`g dd`b |
13 |
He3.cal |
14 |
0 |
15 |
17 amb`r amb`g amb`b |
16 |
amb`r amb`g amb`b |
17 |
0 0 0 |
18 |
sigma0 tau |
19 |
n_real`r n_imag`r |
20 |
n_real`g n_imag`g |
21 |
n_real`b n_imag`b |
22 |
} |
23 |
|
24 |
{ Constants } |
25 |
lambda`r : .67; { red wavelength (microns) } |
26 |
lambda`g : .55; { green wavelength (microns) } |
27 |
lambda`b : .43; { blue wavelength (microns) } |
28 |
z0err : .0001; { accepted error in value of z0 } |
29 |
Dsumlim : .000001; { last term of D summation } |
30 |
Dsummax : 200; { maximum terms in D summation } |
31 |
|
32 |
{ Parameters } |
33 |
sigma0 = arg(10); { surface height deviation (microns) } |
34 |
tau = arg(11); { correlation distance (microns) } |
35 |
n_real`r = arg(12); { red real part of index of refraction } |
36 |
n_imag`r = arg(13); { red imaginary part of index of refraction } |
37 |
n_real`g = arg(14); { green real part of index of refraction } |
38 |
n_imag`g = arg(15); { green imaginary part of index of refraction } |
39 |
n_real`b = arg(16); { blue real part of index of refraction } |
40 |
n_imag`b = arg(17); { blue imaginary part of index of refraction } |
41 |
{ Derived parameters } |
42 |
n_k`r = n_imag`r/n_real`r; |
43 |
n_k`g = n_imag`g/n_real`g; |
44 |
n_k`b = n_imag`b/n_real`b; |
45 |
|
46 |
{ Repeated formulas } |
47 |
cotexp(t) = tau/sigma0/2/tan(t); |
48 |
shadowf2(et,erfcet) = (1-.5*erfcet) / |
49 |
((Exp(-sq(et))/sqrt(PI)/et - erfcet)/2 + 1); |
50 |
shadowf1(t) = or(FTINY-sigma0, .01-abs(t)); |
51 |
shadowf0(t) = abs(t) - (PI/2-.0001); |
52 |
shadowf(t) = if(shadowf0(t), 0, if(shadowf1(t), 1, |
53 |
shadowf2(cotexp(t), erfc(cotexp(t))))); |
54 |
K(t) = if(abs(t)-FTINY, tan(t) * erfc(cotexp(t)), 0); |
55 |
fuvA`r(ct) = sq(n_real`r)*(1-sq(n_k`r)) - (1-sq(ct)); |
56 |
fuvB`r(ct) = sqrt(sq(fuvA`r(ct)) + 4*sq(sq(n_real`r)*n_k`r)); |
57 |
fu2`r(ct) = (fuvA`r(ct) + fuvB`r(ct))/2; |
58 |
fv2`r(ct) = (-fuvA`r(ct) + fuvB`r(ct))/2; |
59 |
fperp2`r(ct) = (sq(ct-sqrt(fu2`r(ct))) + fv2`r(ct)) / |
60 |
(sq(ct+sqrt(fu2`r(ct))) + fv2`r(ct)); |
61 |
fpara2`r(ct) = (sq(sq(n_real`r)*(1-sq(n_k`r))*ct - sqrt(fu2`r(ct))) + |
62 |
sq(2*sq(n_real`r)*n_k`r*ct - sqrt(fv2`r(ct)))) / |
63 |
(sq(sq(n_real`r)*(1-sq(n_k`r))*ct + sqrt(fu2`r(ct))) + |
64 |
sq(2*sq(n_real`r)*n_k`r*ct + sqrt(fv2`r(ct)))); |
65 |
fresnel2`r(ct) = (fperp2`r(ct) + fpara2`r(ct))/2; |
66 |
fuvA`g(ct) = sq(n_real`g)*(1-sq(n_k`g)) - (1-sq(ct)); |
67 |
fuvB`g(ct) = sqrt(sq(fuvA`g(ct)) + 4*sq(sq(n_real`g)*n_k`g)); |
68 |
fu2`g(ct) = (fuvA`g(ct) + fuvB`g(ct))/2; |
69 |
fv2`g(ct) = (-fuvA`g(ct) + fuvB`g(ct))/2; |
70 |
fperp2`g(ct) = (sq(ct-sqrt(fu2`g(ct))) + fv2`g(ct)) / |
71 |
(sq(ct+sqrt(fu2`g(ct))) + fv2`g(ct)); |
72 |
fpara2`g(ct) = (sq(sq(n_real`g)*(1-sq(n_k`g))*ct - sqrt(fu2`g(ct))) + |
73 |
sq(2*sq(n_real`g)*n_k`g*ct - sqrt(fv2`g(ct)))) / |
74 |
(sq(sq(n_real`g)*(1-sq(n_k`g))*ct + sqrt(fu2`g(ct))) + |
75 |
sq(2*sq(n_real`g)*n_k`g*ct + sqrt(fv2`g(ct)))); |
76 |
fresnel2`g(ct) = (fperp2`g(ct) + fpara2`g(ct))/2; |
77 |
fuvA`b(ct) = sq(n_real`b)*(1-sq(n_k`b)) - (1-sq(ct)); |
78 |
fuvB`b(ct) = sqrt(sq(fuvA`b(ct)) + 4*sq(sq(n_real`b)*n_k`b)); |
79 |
fu2`b(ct) = (fuvA`b(ct) + fuvB`b(ct))/2; |
80 |
fv2`b(ct) = (-fuvA`b(ct) + fuvB`b(ct))/2; |
81 |
fperp2`b(ct) = (sq(ct-sqrt(fu2`b(ct))) + fv2`b(ct)) / |
82 |
(sq(ct+sqrt(fu2`b(ct))) + fv2`b(ct)); |
83 |
fpara2`b(ct) = (sq(sq(n_real`b)*(1-sq(n_k`b))*ct - sqrt(fu2`b(ct))) + |
84 |
sq(2*sq(n_real`b)*n_k`b*ct - sqrt(fv2`b(ct)))) / |
85 |
(sq(sq(n_real`b)*(1-sq(n_k`b))*ct + sqrt(fu2`b(ct))) + |
86 |
sq(2*sq(n_real`b)*n_k`b*ct + sqrt(fv2`b(ct)))); |
87 |
fresnel2`b(ct) = (fperp2`b(ct) + fpara2`b(ct))/2; |
88 |
|
89 |
{ Formulas dependent only on reflected direction } |
90 |
theta_r = acos(RdotP); |
91 |
shadowf_r = shadowf(theta_r); |
92 |
K_r = K(theta_r); |
93 |
srx = Dy*NzP - Dz*NyP; sry = Dz*NxP - Dx*NzP; srz = Dx*NyP - Dy*NxP; |
94 |
srn2 = sq(srx) + sq(sry) + sq(srz); |
95 |
prx = sry*Dz - srz*Dy; |
96 |
pry = srz*Dx - srx*Dz; |
97 |
prz = srx*Dy - sry*Dx; |
98 |
s`r = fresnel2`r(RdotP)*Exp(-g`r(RdotP))*sq(shadowf_r); |
99 |
s`g = fresnel2`g(RdotP)*Exp(-g`g(RdotP))*sq(shadowf_r); |
100 |
s`b = fresnel2`b(RdotP)*Exp(-g`b(RdotP))*sq(shadowf_r); |
101 |
|
102 |
{ Formulas dependent on incident direction } |
103 |
{ z0 } |
104 |
z0d(Ki,z) = -(Ki+K_r)/(4*sigma0)*z*Exp(-sq(z/sigma0)/2) - sqrt(PI/2); |
105 |
z0lim(x) = if(x, max(x,z0err), min(x,-z0err)); |
106 |
z0off(Ki,z) = (sigma0/4*(Ki+K_r)*Exp(-sq(z/sigma0)/2)-sqrt(PI/2)*z)/ |
107 |
z0lim(z0d(Ki,z)); |
108 |
z0root(Ki, x0, x1, i) = if(i, |
109 |
if(z0err-abs(x1-x0), |
110 |
x1, |
111 |
z0root(Ki,x1,x1-z0off(Ki,x1),i-1)), |
112 |
0); |
113 |
z0(ti) = z0root(K(ti), .1, -z0off(K(ti),.1), 100); |
114 |
{ sigma } |
115 |
sigma(ti) = if( FTINY-sigma0, sigma0, |
116 |
sigma0/sqrt(1+sq(z0(ti)/sigma0)) ); |
117 |
{ g } |
118 |
g`r(cti) = sq(2*PI/lambda`r*sigma(Acos(cti))*(cti+RdotP)); |
119 |
g`g(cti) = sq(2*PI/lambda`g*sigma(Acos(cti))*(cti+RdotP)); |
120 |
g`b(cti) = sq(2*PI/lambda`b*sigma(Acos(cti))*(cti+RdotP)); |
121 |
{ |F|^2 } |
122 |
fresnel2dd`r(kix,kiy,kiz) = fresnel2`r(sqrt(sq(kix-Dx) + sq(kiy-Dy) + |
123 |
sq(kiz-Dz))/2); |
124 |
fresnel2dd`g(kix,kiy,kiz) = fresnel2`g(sqrt(sq(kix-Dx) + sq(kiy-Dy) + |
125 |
sq(kiz-Dz))/2); |
126 |
fresnel2dd`b(kix,kiy,kiz) = fresnel2`b(sqrt(sq(kix-Dx) + sq(kiy-Dy) + |
127 |
sq(kiz-Dz))/2); |
128 |
{ G } |
129 |
G(kix,kiy,kiz) = sq( (sq(kix-Dx)+sq(kiy-Dy)+sq(kiz-Dz)) / |
130 |
(NxP*(kix-Dx)+NyP*(kiy-Dy)+NzP*(kiz-Dz)) ); |
131 |
{ D } |
132 |
Dsum2(m,lt,c,t,e,g) = if(or(m-Dsummax,and(lt-t,Dsumlim-t)),0, |
133 |
t+Dsum2(m+1,t,c*g/(m+1),c*g/(m+1)*Exp(-g-e/(m+1))/(m+1),e,g)); |
134 |
Dsum(e,g) = Dsum2(1,0,g,g*Exp(-g-e),e,g); |
135 |
D`r(kix,kiy,kiz) = sq(PI)/4/sq(lambda`r)*sq(tau) * |
136 |
Dsum(sq(2*PI/lambda`r)/4*sq(tau)* |
137 |
(sq(kix-Dx)+sq(kiy-Dy)+sq(kiz-Dz) - |
138 |
sq(NxP*(kix-Dx)+NyP*(kiy-Dy)+NzP*(kiz-Dz))), |
139 |
g`r(kix*NxP+kiy*NyP+kiz*NzP)); |
140 |
D`g(kix,kiy,kiz) = sq(PI)/4/sq(lambda`g)*sq(tau) * |
141 |
Dsum(sq(2*PI/lambda`g)/4*sq(tau)* |
142 |
(sq(kix-Dx)+sq(kiy-Dy)+sq(kiz-Dz) - |
143 |
sq(NxP*(kix-Dx)+NyP*(kiy-Dy)+NzP*(kiz-Dz))), |
144 |
g`g(kix*NxP+kiy*NyP+kiz*NzP)); |
145 |
D`b(kix,kiy,kiz) = sq(PI)/4/sq(lambda`b)*sq(tau) * |
146 |
Dsum(sq(2*PI/lambda`b)/4*sq(tau)* |
147 |
(sq(kix-Dx)+sq(kiy-Dy)+sq(kiz-Dz) - |
148 |
sq(NxP*(kix-Dx)+NyP*(kiy-Dy)+NzP*(kiz-Dz))), |
149 |
g`b(kix*NxP+kiy*NyP+kiz*NzP)); |
150 |
{ rho_dd } |
151 |
dd2(cti) = shadowf_r*shadowf(Acos(cti))/cti/RdotP; |
152 |
dd`r(kix,kiy,kiz) = dd2(kix*NxP+kiy*NyP+kiz*NzP)*G(kix,kiy,kiz)* |
153 |
fresnel2dd`r(kix,kiy,kiz)/PI*D`r(kix,kiy,kiz); |
154 |
dd`g(kix,kiy,kiz) = dd2(kix*NxP+kiy*NyP+kiz*NzP)*G(kix,kiy,kiz)* |
155 |
fresnel2dd`g(kix,kiy,kiz)/PI*D`g(kix,kiy,kiz); |
156 |
dd`b(kix,kiy,kiz) = dd2(kix*NxP+kiy*NyP+kiz*NzP)*G(kix,kiy,kiz)* |
157 |
fresnel2dd`b(kix,kiy,kiz)/PI*D`b(kix,kiy,kiz); |