ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/aniso.c
Revision: 2.61
Committed: Wed Sep 2 18:59:01 2015 UTC (8 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, rad5R2, rad5R0, rad5R1, rad5R3
Changes since 2.60: +5 -1 lines
Log Message:
Had to reinstate ambRayInPmap() macro to avoid over-counting bug

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: aniso.c,v 2.60 2015/05/26 13:21:07 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 if (ambRayInPmap(np->rp))
113 return; /* specular accounted for in photon map */
114
115 if (ldot > FTINY && np->specfl&SP_REFL) {
116 /*
117 * Compute specular reflection coefficient using
118 * anisotropic Gaussian distribution model.
119 */
120 /* add source width if flat */
121 if (np->specfl & SP_FLAT)
122 au2 = av2 = omega * (0.25/PI);
123 else
124 au2 = av2 = 0.0;
125 au2 += np->u_alpha*np->u_alpha;
126 av2 += np->v_alpha*np->v_alpha;
127 /* half vector */
128 VSUB(h, ldir, np->rp->rdir);
129 /* ellipse */
130 dtmp1 = DOT(np->u, h);
131 dtmp1 *= dtmp1 / au2;
132 dtmp2 = DOT(np->v, h);
133 dtmp2 *= dtmp2 / av2;
134 /* new W-G-M-D model */
135 dtmp = DOT(np->pnorm, h);
136 dtmp *= dtmp;
137 dtmp1 = (dtmp1 + dtmp2) / dtmp;
138 dtmp = exp(-dtmp1) * DOT(h,h) /
139 (PI * dtmp*dtmp * sqrt(au2*av2));
140 /* worth using? */
141 if (dtmp > FTINY) {
142 copycolor(ctmp, np->scolor);
143 dtmp *= ldot * omega;
144 scalecolor(ctmp, dtmp);
145 addcolor(cval, ctmp);
146 }
147 }
148
149 if (ldot < -FTINY && np->specfl&SP_TRAN) {
150 /*
151 * Compute specular transmission. Specular transmission
152 * is always modified by material color.
153 */
154 /* roughness + source */
155 au2 = av2 = omega * (1.0/PI);
156 au2 += np->u_alpha*np->u_alpha;
157 av2 += np->v_alpha*np->v_alpha;
158 /* "half vector" */
159 VSUB(h, ldir, np->prdir);
160 dtmp = DOT(h,h);
161 if (dtmp > FTINY*FTINY) {
162 dtmp1 = DOT(h,np->pnorm);
163 dtmp = 1.0 - dtmp1*dtmp1/dtmp;
164 if (dtmp > FTINY*FTINY) {
165 dtmp1 = DOT(h,np->u);
166 dtmp1 *= dtmp1 / au2;
167 dtmp2 = DOT(h,np->v);
168 dtmp2 *= dtmp2 / av2;
169 dtmp = (dtmp1 + dtmp2) / dtmp;
170 }
171 } else
172 dtmp = 0.0;
173 /* Gaussian */
174 dtmp = exp(-dtmp) * (1.0/PI) * sqrt(-ldot/(np->pdot*au2*av2));
175 /* worth using? */
176 if (dtmp > FTINY) {
177 copycolor(ctmp, np->mcolor);
178 dtmp *= np->tspec * omega;
179 scalecolor(ctmp, dtmp);
180 addcolor(cval, ctmp);
181 }
182 }
183 }
184
185
186 int
187 m_aniso( /* shade ray that hit something anisotropic */
188 OBJREC *m,
189 RAY *r
190 )
191 {
192 ANISODAT nd;
193 COLOR ctmp;
194 int i;
195 /* easy shadow test */
196 if (r->crtype & SHADOW)
197 return(1);
198
199 if (m->oargs.nfargs != (m->otype == MAT_TRANS2 ? 8 : 6))
200 objerror(m, USER, "bad number of real arguments");
201 /* check for back side */
202 if (r->rod < 0.0) {
203 if (!backvis) {
204 raytrans(r);
205 return(1);
206 }
207 raytexture(r, m->omod);
208 flipsurface(r); /* reorient if backvis */
209 } else
210 raytexture(r, m->omod);
211 /* get material color */
212 nd.mp = m;
213 nd.rp = r;
214 setcolor(nd.mcolor, m->oargs.farg[0],
215 m->oargs.farg[1],
216 m->oargs.farg[2]);
217 /* get roughness */
218 nd.specfl = 0;
219 nd.u_alpha = m->oargs.farg[4];
220 nd.v_alpha = m->oargs.farg[5];
221 if ((nd.u_alpha <= FTINY) | (nd.v_alpha <= FTINY))
222 objerror(m, USER, "roughness too small");
223
224 nd.pdot = raynormal(nd.pnorm, r); /* perturb normal */
225 if (nd.pdot < .001)
226 nd.pdot = .001; /* non-zero for diraniso() */
227 multcolor(nd.mcolor, r->pcol); /* modify material color */
228 /* get specular component */
229 if ((nd.rspec = m->oargs.farg[3]) > FTINY) {
230 nd.specfl |= SP_REFL;
231 /* compute specular color */
232 if (m->otype == MAT_METAL2)
233 copycolor(nd.scolor, nd.mcolor);
234 else
235 setcolor(nd.scolor, 1.0, 1.0, 1.0);
236 scalecolor(nd.scolor, nd.rspec);
237 /* check threshold */
238 if (specthresh >= nd.rspec-FTINY)
239 nd.specfl |= SP_RBLT;
240 /* compute refl. direction */
241 VSUM(nd.vrefl, r->rdir, nd.pnorm, 2.0*nd.pdot);
242 if (DOT(nd.vrefl, r->ron) <= FTINY) /* penetration? */
243 VSUM(nd.vrefl, r->rdir, r->ron, 2.0*r->rod);
244 }
245 /* compute transmission */
246 if (m->otype == MAT_TRANS2) {
247 nd.trans = m->oargs.farg[6]*(1.0 - nd.rspec);
248 nd.tspec = nd.trans * m->oargs.farg[7];
249 nd.tdiff = nd.trans - nd.tspec;
250 if (nd.tspec > FTINY) {
251 nd.specfl |= SP_TRAN;
252 /* check threshold */
253 if (specthresh >= nd.tspec-FTINY)
254 nd.specfl |= SP_TBLT;
255 if (DOT(r->pert,r->pert) <= FTINY*FTINY) {
256 VCOPY(nd.prdir, r->rdir);
257 } else {
258 for (i = 0; i < 3; i++) /* perturb */
259 nd.prdir[i] = r->rdir[i] - r->pert[i];
260 if (DOT(nd.prdir, r->ron) < -FTINY)
261 normalize(nd.prdir); /* OK */
262 else
263 VCOPY(nd.prdir, r->rdir);
264 }
265 }
266 } else
267 nd.tdiff = nd.tspec = nd.trans = 0.0;
268
269 /* diffuse reflection */
270 nd.rdiff = 1.0 - nd.trans - nd.rspec;
271
272 if (r->ro != NULL && isflat(r->ro->otype))
273 nd.specfl |= SP_FLAT;
274
275 getacoords(&nd); /* set up coordinates */
276
277 if (nd.specfl & (SP_REFL|SP_TRAN))
278 agaussamp(&nd);
279
280 if (nd.rdiff > FTINY) { /* ambient from this side */
281 copycolor(ctmp, nd.mcolor); /* modified by material color */
282 scalecolor(ctmp, nd.rdiff);
283 if (nd.specfl & SP_RBLT) /* add in specular as well? */
284 addcolor(ctmp, nd.scolor);
285 multambient(ctmp, r, nd.pnorm);
286 addcolor(r->rcol, ctmp); /* add to returned color */
287 }
288
289 if (nd.tdiff > FTINY) { /* ambient from other side */
290 FVECT bnorm;
291
292 flipsurface(r);
293 bnorm[0] = -nd.pnorm[0];
294 bnorm[1] = -nd.pnorm[1];
295 bnorm[2] = -nd.pnorm[2];
296 copycolor(ctmp, nd.mcolor); /* modified by color */
297 if (nd.specfl & SP_TBLT)
298 scalecolor(ctmp, nd.trans);
299 else
300 scalecolor(ctmp, nd.tdiff);
301 multambient(ctmp, r, bnorm);
302 addcolor(r->rcol, ctmp);
303 flipsurface(r);
304 }
305 /* add direct component */
306 direct(r, diraniso, &nd);
307
308 return(1);
309 }
310
311 static void
312 getacoords( /* set up coordinate system */
313 ANISODAT *np
314 )
315 {
316 MFUNC *mf;
317 int i;
318
319 mf = getfunc(np->mp, 3, 0x7, 1);
320 setfunc(np->mp, np->rp);
321 errno = 0;
322 for (i = 0; i < 3; i++)
323 np->u[i] = evalue(mf->ep[i]);
324 if ((errno == EDOM) | (errno == ERANGE))
325 np->u[0] = np->u[1] = np->u[2] = 0.0;
326 if (mf->fxp != &unitxf)
327 multv3(np->u, np->u, mf->fxp->xfm);
328 fcross(np->v, np->pnorm, np->u);
329 if (normalize(np->v) == 0.0) {
330 if (fabs(np->u_alpha - np->v_alpha) > 0.001)
331 objerror(np->mp, WARNING, "illegal orientation vector");
332 getperpendicular(np->u, np->pnorm, 1); /* punting */
333 fcross(np->v, np->pnorm, np->u);
334 np->u_alpha = np->v_alpha = sqrt( 0.5 *
335 (np->u_alpha*np->u_alpha + np->v_alpha*np->v_alpha) );
336 } else
337 fcross(np->u, np->v, np->pnorm);
338 }
339
340
341 static void
342 agaussamp( /* sample anisotropic Gaussian specular */
343 ANISODAT *np
344 )
345 {
346 RAY sr;
347 FVECT h;
348 double rv[2];
349 double d, sinp, cosp;
350 COLOR scol;
351 int maxiter, ntrials, nstarget, nstaken;
352 int i;
353 /* compute reflection */
354 if ((np->specfl & (SP_REFL|SP_RBLT)) == SP_REFL &&
355 rayorigin(&sr, SPECULAR, np->rp, np->scolor) == 0) {
356 nstarget = 1;
357 if (specjitter > 1.5) { /* multiple samples? */
358 nstarget = specjitter*np->rp->rweight + .5;
359 if (sr.rweight <= minweight*nstarget)
360 nstarget = sr.rweight/minweight;
361 if (nstarget > 1) {
362 d = 1./nstarget;
363 scalecolor(sr.rcoef, d);
364 sr.rweight *= d;
365 } else
366 nstarget = 1;
367 }
368 setcolor(scol, 0., 0., 0.);
369 dimlist[ndims++] = (int)(size_t)np->mp;
370 maxiter = MAXITER*nstarget;
371 for (nstaken = ntrials = 0; nstaken < nstarget &&
372 ntrials < maxiter; ntrials++) {
373 if (ntrials)
374 d = frandom();
375 else
376 d = urand(ilhash(dimlist,ndims)+samplendx);
377 multisamp(rv, 2, d);
378 d = 2.0*PI * rv[0];
379 cosp = tcos(d) * np->u_alpha;
380 sinp = tsin(d) * np->v_alpha;
381 d = 1./sqrt(cosp*cosp + sinp*sinp);
382 cosp *= d;
383 sinp *= d;
384 if ((0. <= specjitter) & (specjitter < 1.))
385 rv[1] = 1.0 - specjitter*rv[1];
386 if (rv[1] <= FTINY)
387 d = 1.0;
388 else
389 d = sqrt(-log(rv[1]) /
390 (cosp*cosp/(np->u_alpha*np->u_alpha) +
391 sinp*sinp/(np->v_alpha*np->v_alpha)));
392 for (i = 0; i < 3; i++)
393 h[i] = np->pnorm[i] +
394 d*(cosp*np->u[i] + sinp*np->v[i]);
395 d = -2.0 * DOT(h, np->rp->rdir) / (1.0 + d*d);
396 VSUM(sr.rdir, np->rp->rdir, h, d);
397 /* sample rejection test */
398 if ((d = DOT(sr.rdir, np->rp->ron)) <= FTINY)
399 continue;
400 checknorm(sr.rdir);
401 if (nstarget > 1) { /* W-G-M-D adjustment */
402 if (nstaken) rayclear(&sr);
403 rayvalue(&sr);
404 d = 2./(1. + np->rp->rod/d);
405 scalecolor(sr.rcol, d);
406 addcolor(scol, sr.rcol);
407 } else {
408 rayvalue(&sr);
409 multcolor(sr.rcol, sr.rcoef);
410 addcolor(np->rp->rcol, sr.rcol);
411 }
412 ++nstaken;
413 }
414 if (nstarget > 1) { /* final W-G-M-D weighting */
415 multcolor(scol, sr.rcoef);
416 d = (double)nstarget/ntrials;
417 scalecolor(scol, d);
418 addcolor(np->rp->rcol, scol);
419 }
420 ndims--;
421 }
422 /* compute transmission */
423 copycolor(sr.rcoef, np->mcolor); /* modify by material color */
424 scalecolor(sr.rcoef, np->tspec);
425 if ((np->specfl & (SP_TRAN|SP_TBLT)) == SP_TRAN &&
426 rayorigin(&sr, SPECULAR, np->rp, sr.rcoef) == 0) {
427 nstarget = 1;
428 if (specjitter > 1.5) { /* multiple samples? */
429 nstarget = specjitter*np->rp->rweight + .5;
430 if (sr.rweight <= minweight*nstarget)
431 nstarget = sr.rweight/minweight;
432 if (nstarget > 1) {
433 d = 1./nstarget;
434 scalecolor(sr.rcoef, d);
435 sr.rweight *= d;
436 } else
437 nstarget = 1;
438 }
439 dimlist[ndims++] = (int)(size_t)np->mp;
440 maxiter = MAXITER*nstarget;
441 for (nstaken = ntrials = 0; nstaken < nstarget &&
442 ntrials < maxiter; ntrials++) {
443 if (ntrials)
444 d = frandom();
445 else
446 d = urand(ilhash(dimlist,ndims)+1823+samplendx);
447 multisamp(rv, 2, d);
448 d = 2.0*PI * rv[0];
449 cosp = tcos(d) * np->u_alpha;
450 sinp = tsin(d) * np->v_alpha;
451 d = 1./sqrt(cosp*cosp + sinp*sinp);
452 cosp *= d;
453 sinp *= d;
454 if ((0. <= specjitter) & (specjitter < 1.))
455 rv[1] = 1.0 - specjitter*rv[1];
456 if (rv[1] <= FTINY)
457 d = 1.0;
458 else
459 d = sqrt(-log(rv[1]) /
460 (cosp*cosp/(np->u_alpha*np->u_alpha) +
461 sinp*sinp/(np->v_alpha*np->v_alpha)));
462 for (i = 0; i < 3; i++)
463 sr.rdir[i] = np->prdir[i] +
464 d*(cosp*np->u[i] + sinp*np->v[i]);
465 if (DOT(sr.rdir, np->rp->ron) >= -FTINY)
466 continue;
467 normalize(sr.rdir); /* OK, normalize */
468 if (nstaken) /* multi-sampling */
469 rayclear(&sr);
470 rayvalue(&sr);
471 multcolor(sr.rcol, sr.rcoef);
472 addcolor(np->rp->rcol, sr.rcol);
473 ++nstaken;
474 }
475 ndims--;
476 }
477 }