| 1 |
greg |
1.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;
|
