19 |
|
extern double specthresh; /* specular sampling threshold */ |
20 |
|
extern double specjitter; /* specular sampling jitter */ |
21 |
|
|
22 |
+ |
static agaussamp(); |
23 |
+ |
|
24 |
|
/* |
25 |
< |
* This anisotropic reflection model uses a variant on the |
26 |
< |
* exponential Gaussian used in normal.c. |
25 |
> |
* This routine implements the anisotropic Gaussian |
26 |
> |
* model described by Ward in Siggraph `92 article. |
27 |
|
* We orient the surface towards the incoming ray, so a single |
28 |
|
* surface can be used to represent an infinitely thin object. |
29 |
|
* |
55 |
|
FVECT vrefl; /* vector in reflected direction */ |
56 |
|
FVECT prdir; /* vector in transmitted direction */ |
57 |
|
FVECT u, v; /* u and v vectors orienting anisotropy */ |
58 |
< |
double u_alpha2; /* u roughness squared */ |
59 |
< |
double v_alpha2; /* v roughness squared */ |
58 |
> |
double u_alpha; /* u roughness */ |
59 |
> |
double v_alpha; /* v roughness */ |
60 |
|
double rdiff, rspec; /* reflected specular, diffuse */ |
61 |
|
double trans; /* transmissivity */ |
62 |
|
double tdiff, tspec; /* transmitted specular, diffuse */ |
105 |
|
au2 = av2 = omega/(4.0*PI); |
106 |
|
else |
107 |
|
au2 = av2 = 0.0; |
108 |
< |
au2 += np->u_alpha2; |
109 |
< |
av2 += np->v_alpha2; |
108 |
> |
au2 += np->u_alpha*np->u_alpha; |
109 |
> |
av2 += np->v_alpha*np->v_alpha; |
110 |
|
/* half vector */ |
111 |
|
h[0] = ldir[0] - np->rp->rdir[0]; |
112 |
|
h[1] = ldir[1] - np->rp->rdir[1]; |
113 |
|
h[2] = ldir[2] - np->rp->rdir[2]; |
112 |
– |
normalize(h); |
114 |
|
/* ellipse */ |
115 |
|
dtmp1 = DOT(np->u, h); |
116 |
|
dtmp1 *= dtmp1 / au2; |
117 |
|
dtmp2 = DOT(np->v, h); |
118 |
|
dtmp2 *= dtmp2 / av2; |
119 |
|
/* gaussian */ |
120 |
< |
dtmp = (dtmp1 + dtmp2) / (1.0 + DOT(np->pnorm, h)); |
121 |
< |
dtmp = exp(-2.0*dtmp) * 1.0/(4.0*PI) |
120 |
> |
dtmp = DOT(np->pnorm, h); |
121 |
> |
dtmp = (dtmp1 + dtmp2) / (dtmp*dtmp); |
122 |
> |
dtmp = exp(-dtmp) * (0.25/PI) |
123 |
|
* sqrt(ldot/(np->pdot*au2*av2)); |
124 |
|
/* worth using? */ |
125 |
|
if (dtmp > FTINY) { |
145 |
|
*/ |
146 |
|
/* roughness + source */ |
147 |
|
au2 = av2 = omega / PI; |
148 |
< |
au2 += np->u_alpha2; |
149 |
< |
av2 += np->v_alpha2; |
148 |
> |
au2 += np->u_alpha*np->u_alpha; |
149 |
> |
av2 += np->v_alpha*np->v_alpha; |
150 |
|
/* "half vector" */ |
151 |
|
h[0] = ldir[0] - np->prdir[0]; |
152 |
|
h[1] = ldir[1] - np->prdir[1]; |
153 |
|
h[2] = ldir[2] - np->prdir[2]; |
154 |
< |
dtmp = DOT(h,np->pnorm); |
153 |
< |
dtmp = DOT(h,h) - dtmp*dtmp; |
154 |
> |
dtmp = DOT(h,h); |
155 |
|
if (dtmp > FTINY*FTINY) { |
156 |
< |
dtmp1 = DOT(h,np->u); |
157 |
< |
dtmp1 = dtmp1*dtmp1 / (au2*dtmp); |
158 |
< |
dtmp2 = DOT(h,np->v); |
159 |
< |
dtmp2 = dtmp2*dtmp2 / (av2*dtmp); |
160 |
< |
dtmp = 2. - 2.*DOT(ldir,np->prdir); |
161 |
< |
dtmp *= dtmp1 + dtmp2; |
156 |
> |
dtmp1 = DOT(h,np->pnorm); |
157 |
> |
dtmp = 1.0 - dtmp1*dtmp1/dtmp; |
158 |
> |
if (dtmp > FTINY*FTINY) { |
159 |
> |
dtmp1 = DOT(h,np->u); |
160 |
> |
dtmp1 *= dtmp1 / au2; |
161 |
> |
dtmp2 = DOT(h,np->v); |
162 |
> |
dtmp2 *= dtmp2 / av2; |
163 |
> |
dtmp = (dtmp1 + dtmp2) / dtmp; |
164 |
> |
} |
165 |
|
} else |
166 |
|
dtmp = 0.0; |
167 |
|
/* gaussian */ |
168 |
< |
dtmp = exp(-dtmp) * 1.0/(4.0*PI) |
168 |
> |
dtmp = exp(-dtmp) * (1.0/PI) |
169 |
|
* sqrt(-ldot/(np->pdot*au2*av2)); |
170 |
|
/* worth using? */ |
171 |
|
if (dtmp > FTINY) { |
200 |
|
m->oargs.farg[2]); |
201 |
|
/* get roughness */ |
202 |
|
nd.specfl = 0; |
203 |
< |
nd.u_alpha2 = m->oargs.farg[4]; |
204 |
< |
nd.u_alpha2 *= nd.u_alpha2; |
205 |
< |
nd.v_alpha2 = m->oargs.farg[5]; |
202 |
< |
nd.v_alpha2 *= nd.v_alpha2; |
203 |
< |
if (nd.u_alpha2 < FTINY*FTINY || nd.v_alpha2 <= FTINY*FTINY) |
203 |
> |
nd.u_alpha = m->oargs.farg[4]; |
204 |
> |
nd.v_alpha = m->oargs.farg[5]; |
205 |
> |
if (nd.u_alpha < FTINY || nd.v_alpha <= FTINY) |
206 |
|
objerror(m, USER, "roughness too small"); |
207 |
|
/* reorient if necessary */ |
208 |
|
if (r->rod < 0.0) |
349 |
|
d = urand(ilhash(dimlist,ndims)+samplendx); |
350 |
|
multisamp(rv, 2, d); |
351 |
|
d = 2.0*PI * rv[0]; |
352 |
< |
cosp = cos(d); |
353 |
< |
sinp = sin(d); |
354 |
< |
d = sqrt(np->u_alpha2*cosp*cosp + np->v_alpha2*sinp*sinp); |
352 |
> |
cosp = cos(d) * np->u_alpha; |
353 |
> |
sinp = sin(d) * np->v_alpha; |
354 |
> |
d = sqrt(cosp*cosp + sinp*sinp); |
355 |
|
cosp /= d; |
356 |
|
sinp /= d; |
357 |
|
rv[1] = 1.0 - specjitter*rv[1]; |
359 |
|
d = 1.0; |
360 |
|
else |
361 |
|
d = sqrt(-log(rv[1]) / |
362 |
< |
(cosp*cosp/np->u_alpha2 + |
363 |
< |
sinp*sinp/np->v_alpha2)); |
362 |
> |
(cosp*cosp/(np->u_alpha*np->u_alpha) + |
363 |
> |
sinp*sinp/(np->v_alpha*np->v_alpha))); |
364 |
|
for (i = 0; i < 3; i++) |
365 |
|
h[i] = np->pnorm[i] + |
366 |
|
d*(cosp*np->u[i] + sinp*np->v[i]); |
381 |
|
d = urand(ilhash(dimlist,ndims)+1823+samplendx); |
382 |
|
multisamp(rv, 2, d); |
383 |
|
d = 2.0*PI * rv[0]; |
384 |
< |
cosp = cos(d); |
385 |
< |
sinp = sin(d); |
384 |
> |
cosp = cos(d) * np->u_alpha; |
385 |
> |
sinp = sin(d) * np->v_alpha; |
386 |
> |
d = sqrt(cosp*cosp + sinp*sinp); |
387 |
> |
cosp /= d; |
388 |
> |
sinp /= d; |
389 |
|
rv[1] = 1.0 - specjitter*rv[1]; |
390 |
|
if (rv[1] <= FTINY) |
391 |
|
d = 1.0; |
392 |
|
else |
393 |
|
d = sqrt(-log(rv[1]) / |
394 |
< |
(cosp*cosp*4./np->u_alpha2 + |
395 |
< |
sinp*sinp*4./np->v_alpha2)); |
394 |
> |
(cosp*cosp/(np->u_alpha*np->u_alpha) + |
395 |
> |
sinp*sinp/(np->v_alpha*np->u_alpha))); |
396 |
|
for (i = 0; i < 3; i++) |
397 |
|
sr.rdir[i] = np->prdir[i] + |
398 |
|
d*(cosp*np->u[i] + sinp*np->v[i]); |