ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/aniso.c
Revision: 2.59
Committed: Thu May 21 05:54:54 2015 UTC (9 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.58: +2 -2 lines
Log Message:
Made axis randomization optional in getperpendicular()

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: aniso.c,v 2.58 2015/02/24 19:39:26 greg Exp $";
3 #endif
4 /*
5 * Shading functions for anisotropic materials.
6 */
7
8 #include "copyright.h"
9
10 #include "ray.h"
11 #include "ambient.h"
12 #include "otypes.h"
13 #include "rtotypes.h"
14 #include "source.h"
15 #include "func.h"
16 #include "random.h"
17 #include "pmapmat.h"
18
19 #ifndef MAXITER
20 #define MAXITER 10 /* maximum # specular ray attempts */
21 #endif
22
23 /*
24 * This routine implements the anisotropic Gaussian
25 * model described by Ward in Siggraph `92 article, updated with
26 * normalization and sampling adjustments due to Geisler-Moroder and Duer.
27 * We orient the surface towards the incoming ray, so a single
28 * surface can be used to represent an infinitely thin object.
29 *
30 * Arguments for MAT_PLASTIC2 and MAT_METAL2 are:
31 * 4+ ux uy uz funcfile [transform...]
32 * 0
33 * 6 red grn blu specular-frac. u-rough v-rough
34 *
35 * Real arguments for MAT_TRANS2 are:
36 * 8 red grn blu rspec u-rough v-rough trans tspec
37 */
38
39 /* specularity flags */
40 #define SP_REFL 01 /* has reflected specular component */
41 #define SP_TRAN 02 /* has transmitted specular */
42 #define SP_FLAT 04 /* reflecting surface is flat */
43 #define SP_RBLT 010 /* reflection below sample threshold */
44 #define SP_TBLT 020 /* transmission below threshold */
45
46 typedef struct {
47 OBJREC *mp; /* material pointer */
48 RAY *rp; /* ray pointer */
49 short specfl; /* specularity flags, defined above */
50 COLOR mcolor; /* color of this material */
51 COLOR scolor; /* color of specular component */
52 FVECT vrefl; /* vector in reflected direction */
53 FVECT prdir; /* vector in transmitted direction */
54 FVECT u, v; /* u and v vectors orienting anisotropy */
55 double u_alpha; /* u roughness */
56 double v_alpha; /* v roughness */
57 double rdiff, rspec; /* reflected specular, diffuse */
58 double trans; /* transmissivity */
59 double tdiff, tspec; /* transmitted specular, diffuse */
60 FVECT pnorm; /* perturbed surface normal */
61 double pdot; /* perturbed dot product */
62 } ANISODAT; /* anisotropic material data */
63
64 static void getacoords(ANISODAT *np);
65 static void agaussamp(ANISODAT *np);
66
67
68 static void
69 diraniso( /* compute source contribution */
70 COLOR cval, /* returned coefficient */
71 void *nnp, /* material data */
72 FVECT ldir, /* light source direction */
73 double omega /* light source size */
74 )
75 {
76 ANISODAT *np = nnp;
77 double ldot;
78 double dtmp, dtmp1, dtmp2;
79 FVECT h;
80 double au2, av2;
81 COLOR ctmp;
82
83 setcolor(cval, 0.0, 0.0, 0.0);
84
85 ldot = DOT(np->pnorm, ldir);
86
87 if (ldot < 0.0 ? np->trans <= FTINY : np->trans >= 1.0-FTINY)
88 return; /* wrong side */
89
90 if ((ldot > FTINY) & (np->rdiff > FTINY)) {
91 /*
92 * Compute and add diffuse reflected component to returned
93 * color. The diffuse reflected component will always be
94 * modified by the color of the material.
95 */
96 copycolor(ctmp, np->mcolor);
97 dtmp = ldot * omega * np->rdiff * (1.0/PI);
98 scalecolor(ctmp, dtmp);
99 addcolor(cval, ctmp);
100 }
101
102 if ((ldot < -FTINY) & (np->tdiff > FTINY)) {
103 /*
104 * Compute diffuse transmission.
105 */
106 copycolor(ctmp, np->mcolor);
107 dtmp = -ldot * omega * np->tdiff * (1.0/PI);
108 scalecolor(ctmp, dtmp);
109 addcolor(cval, ctmp);
110 }
111
112 /* PMAP: skip direct specular refl via ambient bounce if already
113 * accounted for in photon map */
114 if (ambRayInPmap(np->rp))
115 return;
116
117 if (ldot > FTINY && np->specfl&SP_REFL) {
118 /*
119 * Compute specular reflection coefficient using
120 * anisotropic Gaussian distribution model.
121 */
122 /* add source width if flat */
123 if (np->specfl & SP_FLAT)
124 au2 = av2 = omega * (0.25/PI);
125 else
126 au2 = av2 = 0.0;
127 au2 += np->u_alpha*np->u_alpha;
128 av2 += np->v_alpha*np->v_alpha;
129 /* half vector */
130 VSUB(h, ldir, np->rp->rdir);
131 /* ellipse */
132 dtmp1 = DOT(np->u, h);
133 dtmp1 *= dtmp1 / au2;
134 dtmp2 = DOT(np->v, h);
135 dtmp2 *= dtmp2 / av2;
136 /* new W-G-M-D model */
137 dtmp = DOT(np->pnorm, h);
138 dtmp *= dtmp;
139 dtmp1 = (dtmp1 + dtmp2) / dtmp;
140 dtmp = exp(-dtmp1) * DOT(h,h) /
141 (PI * dtmp*dtmp * sqrt(au2*av2));
142 /* worth using? */
143 if (dtmp > FTINY) {
144 copycolor(ctmp, np->scolor);
145 dtmp *= ldot * omega;
146 scalecolor(ctmp, dtmp);
147 addcolor(cval, ctmp);
148 }
149 }
150
151 if (ldot < -FTINY && np->specfl&SP_TRAN) {
152 /*
153 * Compute specular transmission. Specular transmission
154 * is always modified by material color.
155 */
156 /* roughness + source */
157 au2 = av2 = omega * (1.0/PI);
158 au2 += np->u_alpha*np->u_alpha;
159 av2 += np->v_alpha*np->v_alpha;
160 /* "half vector" */
161 VSUB(h, ldir, np->prdir);
162 dtmp = DOT(h,h);
163 if (dtmp > FTINY*FTINY) {
164 dtmp1 = DOT(h,np->pnorm);
165 dtmp = 1.0 - dtmp1*dtmp1/dtmp;
166 if (dtmp > FTINY*FTINY) {
167 dtmp1 = DOT(h,np->u);
168 dtmp1 *= dtmp1 / au2;
169 dtmp2 = DOT(h,np->v);
170 dtmp2 *= dtmp2 / av2;
171 dtmp = (dtmp1 + dtmp2) / dtmp;
172 }
173 } else
174 dtmp = 0.0;
175 /* Gaussian */
176 dtmp = exp(-dtmp) * (1.0/PI) * sqrt(-ldot/(np->pdot*au2*av2));
177 /* worth using? */
178 if (dtmp > FTINY) {
179 copycolor(ctmp, np->mcolor);
180 dtmp *= np->tspec * omega;
181 scalecolor(ctmp, dtmp);
182 addcolor(cval, ctmp);
183 }
184 }
185 }
186
187
188 int
189 m_aniso( /* shade ray that hit something anisotropic */
190 OBJREC *m,
191 RAY *r
192 )
193 {
194 ANISODAT nd;
195 COLOR ctmp;
196 int i;
197 /* easy shadow test */
198 if (r->crtype & SHADOW)
199 return(1);
200
201 if (m->oargs.nfargs != (m->otype == MAT_TRANS2 ? 8 : 6))
202 objerror(m, USER, "bad number of real arguments");
203 /* check for back side */
204 if (r->rod < 0.0) {
205 if (!backvis) {
206 raytrans(r);
207 return(1);
208 }
209 raytexture(r, m->omod);
210 flipsurface(r); /* reorient if backvis */
211 } else
212 raytexture(r, m->omod);
213 /* get material color */
214 nd.mp = m;
215 nd.rp = r;
216 setcolor(nd.mcolor, m->oargs.farg[0],
217 m->oargs.farg[1],
218 m->oargs.farg[2]);
219 /* get roughness */
220 nd.specfl = 0;
221 nd.u_alpha = m->oargs.farg[4];
222 nd.v_alpha = m->oargs.farg[5];
223 if ((nd.u_alpha <= FTINY) | (nd.v_alpha <= FTINY))
224 objerror(m, USER, "roughness too small");
225
226 nd.pdot = raynormal(nd.pnorm, r); /* perturb normal */
227 if (nd.pdot < .001)
228 nd.pdot = .001; /* non-zero for diraniso() */
229 multcolor(nd.mcolor, r->pcol); /* modify material color */
230 /* get specular component */
231 if ((nd.rspec = m->oargs.farg[3]) > FTINY) {
232 nd.specfl |= SP_REFL;
233 /* compute specular color */
234 if (m->otype == MAT_METAL2)
235 copycolor(nd.scolor, nd.mcolor);
236 else
237 setcolor(nd.scolor, 1.0, 1.0, 1.0);
238 scalecolor(nd.scolor, nd.rspec);
239 /* check threshold */
240 if (specthresh >= nd.rspec-FTINY)
241 nd.specfl |= SP_RBLT;
242 /* compute refl. direction */
243 VSUM(nd.vrefl, r->rdir, nd.pnorm, 2.0*nd.pdot);
244 if (DOT(nd.vrefl, r->ron) <= FTINY) /* penetration? */
245 VSUM(nd.vrefl, r->rdir, r->ron, 2.0*r->rod);
246 }
247 /* compute transmission */
248 if (m->otype == MAT_TRANS2) {
249 nd.trans = m->oargs.farg[6]*(1.0 - nd.rspec);
250 nd.tspec = nd.trans * m->oargs.farg[7];
251 nd.tdiff = nd.trans - nd.tspec;
252 if (nd.tspec > FTINY) {
253 nd.specfl |= SP_TRAN;
254 /* check threshold */
255 if (specthresh >= nd.tspec-FTINY)
256 nd.specfl |= SP_TBLT;
257 if (DOT(r->pert,r->pert) <= FTINY*FTINY) {
258 VCOPY(nd.prdir, r->rdir);
259 } else {
260 for (i = 0; i < 3; i++) /* perturb */
261 nd.prdir[i] = r->rdir[i] - r->pert[i];
262 if (DOT(nd.prdir, r->ron) < -FTINY)
263 normalize(nd.prdir); /* OK */
264 else
265 VCOPY(nd.prdir, r->rdir);
266 }
267 }
268 } else
269 nd.tdiff = nd.tspec = nd.trans = 0.0;
270
271 /* diffuse reflection */
272 nd.rdiff = 1.0 - nd.trans - nd.rspec;
273
274 if (r->ro != NULL && isflat(r->ro->otype))
275 nd.specfl |= SP_FLAT;
276
277 getacoords(&nd); /* set up coordinates */
278
279 /* PMAP: skip indirect specular via ambient bounce if already accounted
280 * for in photon map */
281 if (!ambRayInPmap(r) && nd.specfl & (SP_REFL|SP_TRAN))
282 agaussamp(&nd);
283
284 if (nd.rdiff > FTINY) { /* ambient from this side */
285 copycolor(ctmp, nd.mcolor); /* modified by material color */
286 scalecolor(ctmp, nd.rdiff);
287 if (nd.specfl & SP_RBLT) /* add in specular as well? */
288 addcolor(ctmp, nd.scolor);
289 multambient(ctmp, r, nd.pnorm);
290 addcolor(r->rcol, ctmp); /* add to returned color */
291 }
292
293 if (nd.tdiff > FTINY) { /* ambient from other side */
294 FVECT bnorm;
295
296 flipsurface(r);
297 bnorm[0] = -nd.pnorm[0];
298 bnorm[1] = -nd.pnorm[1];
299 bnorm[2] = -nd.pnorm[2];
300 copycolor(ctmp, nd.mcolor); /* modified by color */
301 if (nd.specfl & SP_TBLT)
302 scalecolor(ctmp, nd.trans);
303 else
304 scalecolor(ctmp, nd.tdiff);
305 multambient(ctmp, r, bnorm);
306 addcolor(r->rcol, ctmp);
307 flipsurface(r);
308 }
309 /* add direct component */
310 direct(r, diraniso, &nd);
311
312 return(1);
313 }
314
315 static void
316 getacoords( /* set up coordinate system */
317 ANISODAT *np
318 )
319 {
320 MFUNC *mf;
321 int i;
322
323 mf = getfunc(np->mp, 3, 0x7, 1);
324 setfunc(np->mp, np->rp);
325 errno = 0;
326 for (i = 0; i < 3; i++)
327 np->u[i] = evalue(mf->ep[i]);
328 if ((errno == EDOM) | (errno == ERANGE))
329 np->u[0] = np->u[1] = np->u[2] = 0.0;
330 if (mf->fxp != &unitxf)
331 multv3(np->u, np->u, mf->fxp->xfm);
332 fcross(np->v, np->pnorm, np->u);
333 if (normalize(np->v) == 0.0) {
334 if (fabs(np->u_alpha - np->v_alpha) > 0.001)
335 objerror(np->mp, WARNING, "illegal orientation vector");
336 getperpendicular(np->u, np->pnorm, 1); /* punting */
337 fcross(np->v, np->pnorm, np->u);
338 np->u_alpha = np->v_alpha = sqrt( 0.5 *
339 (np->u_alpha*np->u_alpha + np->v_alpha*np->v_alpha) );
340 } else
341 fcross(np->u, np->v, np->pnorm);
342 }
343
344
345 static void
346 agaussamp( /* sample anisotropic Gaussian specular */
347 ANISODAT *np
348 )
349 {
350 RAY sr;
351 FVECT h;
352 double rv[2];
353 double d, sinp, cosp;
354 COLOR scol;
355 int maxiter, ntrials, nstarget, nstaken;
356 int i;
357 /* compute reflection */
358 if ((np->specfl & (SP_REFL|SP_RBLT)) == SP_REFL &&
359 rayorigin(&sr, SPECULAR, np->rp, np->scolor) == 0) {
360 nstarget = 1;
361 if (specjitter > 1.5) { /* multiple samples? */
362 nstarget = specjitter*np->rp->rweight + .5;
363 if (sr.rweight <= minweight*nstarget)
364 nstarget = sr.rweight/minweight;
365 if (nstarget > 1) {
366 d = 1./nstarget;
367 scalecolor(sr.rcoef, d);
368 sr.rweight *= d;
369 } else
370 nstarget = 1;
371 }
372 setcolor(scol, 0., 0., 0.);
373 dimlist[ndims++] = (int)(size_t)np->mp;
374 maxiter = MAXITER*nstarget;
375 for (nstaken = ntrials = 0; nstaken < nstarget &&
376 ntrials < maxiter; ntrials++) {
377 if (ntrials)
378 d = frandom();
379 else
380 d = urand(ilhash(dimlist,ndims)+samplendx);
381 multisamp(rv, 2, d);
382 d = 2.0*PI * rv[0];
383 cosp = tcos(d) * np->u_alpha;
384 sinp = tsin(d) * np->v_alpha;
385 d = 1./sqrt(cosp*cosp + sinp*sinp);
386 cosp *= d;
387 sinp *= d;
388 if ((0. <= specjitter) & (specjitter < 1.))
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/(np->u_alpha*np->u_alpha) +
395 sinp*sinp/(np->v_alpha*np->v_alpha)));
396 for (i = 0; i < 3; i++)
397 h[i] = np->pnorm[i] +
398 d*(cosp*np->u[i] + sinp*np->v[i]);
399 d = -2.0 * DOT(h, np->rp->rdir) / (1.0 + d*d);
400 VSUM(sr.rdir, np->rp->rdir, h, d);
401 /* sample rejection test */
402 if ((d = DOT(sr.rdir, np->rp->ron)) <= FTINY)
403 continue;
404 checknorm(sr.rdir);
405 if (nstarget > 1) { /* W-G-M-D adjustment */
406 if (nstaken) rayclear(&sr);
407 rayvalue(&sr);
408 d = 2./(1. + np->rp->rod/d);
409 scalecolor(sr.rcol, d);
410 addcolor(scol, sr.rcol);
411 } else {
412 rayvalue(&sr);
413 multcolor(sr.rcol, sr.rcoef);
414 addcolor(np->rp->rcol, sr.rcol);
415 }
416 ++nstaken;
417 }
418 if (nstarget > 1) { /* final W-G-M-D weighting */
419 multcolor(scol, sr.rcoef);
420 d = (double)nstarget/ntrials;
421 scalecolor(scol, d);
422 addcolor(np->rp->rcol, scol);
423 }
424 ndims--;
425 }
426 /* compute transmission */
427 copycolor(sr.rcoef, np->mcolor); /* modify by material color */
428 scalecolor(sr.rcoef, np->tspec);
429 if ((np->specfl & (SP_TRAN|SP_TBLT)) == SP_TRAN &&
430 rayorigin(&sr, SPECULAR, np->rp, sr.rcoef) == 0) {
431 nstarget = 1;
432 if (specjitter > 1.5) { /* multiple samples? */
433 nstarget = specjitter*np->rp->rweight + .5;
434 if (sr.rweight <= minweight*nstarget)
435 nstarget = sr.rweight/minweight;
436 if (nstarget > 1) {
437 d = 1./nstarget;
438 scalecolor(sr.rcoef, d);
439 sr.rweight *= d;
440 } else
441 nstarget = 1;
442 }
443 dimlist[ndims++] = (int)(size_t)np->mp;
444 maxiter = MAXITER*nstarget;
445 for (nstaken = ntrials = 0; nstaken < nstarget &&
446 ntrials < maxiter; ntrials++) {
447 if (ntrials)
448 d = frandom();
449 else
450 d = urand(ilhash(dimlist,ndims)+1823+samplendx);
451 multisamp(rv, 2, d);
452 d = 2.0*PI * rv[0];
453 cosp = tcos(d) * np->u_alpha;
454 sinp = tsin(d) * np->v_alpha;
455 d = 1./sqrt(cosp*cosp + sinp*sinp);
456 cosp *= d;
457 sinp *= d;
458 if ((0. <= specjitter) & (specjitter < 1.))
459 rv[1] = 1.0 - specjitter*rv[1];
460 if (rv[1] <= FTINY)
461 d = 1.0;
462 else
463 d = sqrt(-log(rv[1]) /
464 (cosp*cosp/(np->u_alpha*np->u_alpha) +
465 sinp*sinp/(np->v_alpha*np->v_alpha)));
466 for (i = 0; i < 3; i++)
467 sr.rdir[i] = np->prdir[i] +
468 d*(cosp*np->u[i] + sinp*np->v[i]);
469 if (DOT(sr.rdir, np->rp->ron) >= -FTINY)
470 continue;
471 normalize(sr.rdir); /* OK, normalize */
472 if (nstaken) /* multi-sampling */
473 rayclear(&sr);
474 rayvalue(&sr);
475 multcolor(sr.rcol, sr.rcoef);
476 addcolor(np->rp->rcol, sr.rcol);
477 ++nstaken;
478 }
479 ndims--;
480 }
481 }