ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/aniso.c
Revision: 2.60
Committed: Tue May 26 13:21:07 2015 UTC (8 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.59: +2 -10 lines
Log Message:
Removed deprecated ambRayInPmap() macro from code

File Contents

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