ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/aniso.c
Revision: 2.64
Committed: Fri Apr 5 01:10:26 2024 UTC (5 weeks, 5 days ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.63: +3 -6 lines
Log Message:
fix: Improved tracking of reflected vs. transmitted rays for antimatter

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: aniso.c,v 2.63 2023/11/17 20:02: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 SCOLOR mcolor; /* color of this material */
51 SCOLOR 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 SCOLOR scval, /* 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 SCOLOR sctmp;
82
83 scolorblack(scval);
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 copyscolor(sctmp, np->mcolor);
97 dtmp = ldot * omega * np->rdiff * (1.0/PI);
98 scalescolor(sctmp, dtmp);
99 saddscolor(scval, sctmp);
100 }
101
102 if ((ldot < -FTINY) & (np->tdiff > FTINY)) {
103 /*
104 * Compute diffuse transmission.
105 */
106 copyscolor(sctmp, np->mcolor);
107 dtmp = -ldot * omega * np->tdiff * (1.0/PI);
108 scalescolor(sctmp, dtmp);
109 saddscolor(scval, sctmp);
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 copyscolor(sctmp, np->scolor);
143 dtmp *= ldot * omega;
144 scalescolor(sctmp, dtmp);
145 saddscolor(scval, sctmp);
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 copyscolor(sctmp, np->mcolor);
178 dtmp *= np->tspec * omega;
179 scalescolor(sctmp, dtmp);
180 saddscolor(scval, sctmp);
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 SCOLOR sctmp;
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 setscolor(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 smultscolor(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 copyscolor(nd.scolor, nd.mcolor);
234 else
235 setscolor(nd.scolor, 1.0, 1.0, 1.0);
236 scalescolor(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 copyscolor(sctmp, nd.mcolor); /* modified by material color */
282 scalescolor(sctmp, nd.rdiff);
283 if (nd.specfl & SP_RBLT) /* add in specular as well? */
284 saddscolor(sctmp, nd.scolor);
285 multambient(sctmp, r, nd.pnorm);
286 saddscolor(r->rcol, sctmp); /* add to returned color */
287 }
288
289 if (nd.tdiff > FTINY) { /* ambient from other side */
290 FVECT bnorm;
291 bnorm[0] = -nd.pnorm[0];
292 bnorm[1] = -nd.pnorm[1];
293 bnorm[2] = -nd.pnorm[2];
294 copyscolor(sctmp, nd.mcolor); /* modified by color */
295 if (nd.specfl & SP_TBLT) {
296 scalescolor(sctmp, nd.trans);
297 } else {
298 scalescolor(sctmp, nd.tdiff);
299 }
300 multambient(sctmp, r, bnorm);
301 saddscolor(r->rcol, sctmp);
302 }
303 /* add direct component */
304 direct(r, diraniso, &nd);
305
306 return(1);
307 }
308
309 static void
310 getacoords( /* set up coordinate system */
311 ANISODAT *np
312 )
313 {
314 MFUNC *mf;
315 int i;
316
317 mf = getfunc(np->mp, 3, 0x7, 1);
318 setfunc(np->mp, np->rp);
319 errno = 0;
320 for (i = 0; i < 3; i++)
321 np->u[i] = evalue(mf->ep[i]);
322 if ((errno == EDOM) | (errno == ERANGE))
323 np->u[0] = np->u[1] = np->u[2] = 0.0;
324 if (mf->fxp != &unitxf)
325 multv3(np->u, np->u, mf->fxp->xfm);
326 fcross(np->v, np->pnorm, np->u);
327 if (normalize(np->v) == 0.0) {
328 if (fabs(np->u_alpha - np->v_alpha) > 0.001)
329 objerror(np->mp, WARNING, "illegal orientation vector");
330 getperpendicular(np->u, np->pnorm, 1); /* punting */
331 fcross(np->v, np->pnorm, np->u);
332 np->u_alpha = np->v_alpha = sqrt( 0.5 *
333 (np->u_alpha*np->u_alpha + np->v_alpha*np->v_alpha) );
334 } else
335 fcross(np->u, np->v, np->pnorm);
336 }
337
338
339 static void
340 agaussamp( /* sample anisotropic Gaussian specular */
341 ANISODAT *np
342 )
343 {
344 RAY sr;
345 FVECT h;
346 double rv[2];
347 double d, sinp, cosp;
348 SCOLOR scol;
349 int maxiter, ntrials, nstarget, nstaken;
350 int i;
351 /* compute reflection */
352 if ((np->specfl & (SP_REFL|SP_RBLT)) == SP_REFL &&
353 rayorigin(&sr, RSPECULAR, np->rp, np->scolor) == 0) {
354 nstarget = 1;
355 if (specjitter > 1.5) { /* multiple samples? */
356 nstarget = specjitter*np->rp->rweight + .5;
357 if (sr.rweight <= minweight*nstarget)
358 nstarget = sr.rweight/minweight;
359 if (nstarget > 1) {
360 d = 1./nstarget;
361 scalescolor(sr.rcoef, d);
362 sr.rweight *= d;
363 } else
364 nstarget = 1;
365 }
366 scolorblack(scol);
367 dimlist[ndims++] = (int)(size_t)np->mp;
368 maxiter = MAXITER*nstarget;
369 for (nstaken = ntrials = 0; nstaken < nstarget &&
370 ntrials < maxiter; ntrials++) {
371 if (ntrials)
372 d = frandom();
373 else
374 d = urand(ilhash(dimlist,ndims)+samplendx);
375 multisamp(rv, 2, d);
376 d = 2.0*PI * rv[0];
377 cosp = tcos(d) * np->u_alpha;
378 sinp = tsin(d) * np->v_alpha;
379 d = 1./sqrt(cosp*cosp + sinp*sinp);
380 cosp *= d;
381 sinp *= d;
382 if ((0. <= specjitter) & (specjitter < 1.))
383 rv[1] = 1.0 - specjitter*rv[1];
384 if (rv[1] <= FTINY)
385 d = 1.0;
386 else
387 d = sqrt(-log(rv[1]) /
388 (cosp*cosp/(np->u_alpha*np->u_alpha) +
389 sinp*sinp/(np->v_alpha*np->v_alpha)));
390 for (i = 0; i < 3; i++)
391 h[i] = np->pnorm[i] +
392 d*(cosp*np->u[i] + sinp*np->v[i]);
393 d = -2.0 * DOT(h, np->rp->rdir) / (1.0 + d*d);
394 VSUM(sr.rdir, np->rp->rdir, h, d);
395 /* sample rejection test */
396 if ((d = DOT(sr.rdir, np->rp->ron)) <= FTINY)
397 continue;
398 checknorm(sr.rdir);
399 if (nstarget > 1) { /* W-G-M-D adjustment */
400 if (nstaken) rayclear(&sr);
401 rayvalue(&sr);
402 d = 2./(1. + np->rp->rod/d);
403 scalescolor(sr.rcol, d);
404 saddscolor(scol, sr.rcol);
405 } else {
406 rayvalue(&sr);
407 smultscolor(sr.rcol, sr.rcoef);
408 saddscolor(np->rp->rcol, sr.rcol);
409 }
410 ++nstaken;
411 }
412 if (nstarget > 1) { /* final W-G-M-D weighting */
413 smultscolor(scol, sr.rcoef);
414 d = (double)nstarget/ntrials;
415 scalescolor(scol, d);
416 saddscolor(np->rp->rcol, scol);
417 }
418 ndims--;
419 }
420 /* compute transmission */
421 copyscolor(sr.rcoef, np->mcolor); /* modify by material color */
422 scalescolor(sr.rcoef, np->tspec);
423 if ((np->specfl & (SP_TRAN|SP_TBLT)) == SP_TRAN &&
424 rayorigin(&sr, TSPECULAR, np->rp, sr.rcoef) == 0) {
425 nstarget = 1;
426 if (specjitter > 1.5) { /* multiple samples? */
427 nstarget = specjitter*np->rp->rweight + .5;
428 if (sr.rweight <= minweight*nstarget)
429 nstarget = sr.rweight/minweight;
430 if (nstarget > 1) {
431 d = 1./nstarget;
432 scalescolor(sr.rcoef, d);
433 sr.rweight *= d;
434 } else
435 nstarget = 1;
436 }
437 dimlist[ndims++] = (int)(size_t)np->mp;
438 maxiter = MAXITER*nstarget;
439 for (nstaken = ntrials = 0; nstaken < nstarget &&
440 ntrials < maxiter; ntrials++) {
441 if (ntrials)
442 d = frandom();
443 else
444 d = urand(ilhash(dimlist,ndims)+1823+samplendx);
445 multisamp(rv, 2, d);
446 d = 2.0*PI * rv[0];
447 cosp = tcos(d) * np->u_alpha;
448 sinp = tsin(d) * np->v_alpha;
449 d = 1./sqrt(cosp*cosp + sinp*sinp);
450 cosp *= d;
451 sinp *= d;
452 if ((0. <= specjitter) & (specjitter < 1.))
453 rv[1] = 1.0 - specjitter*rv[1];
454 if (rv[1] <= FTINY)
455 d = 1.0;
456 else
457 d = sqrt(-log(rv[1]) /
458 (cosp*cosp/(np->u_alpha*np->u_alpha) +
459 sinp*sinp/(np->v_alpha*np->v_alpha)));
460 for (i = 0; i < 3; i++)
461 sr.rdir[i] = np->prdir[i] +
462 d*(cosp*np->u[i] + sinp*np->v[i]);
463 if (DOT(sr.rdir, np->rp->ron) >= -FTINY)
464 continue;
465 normalize(sr.rdir); /* OK, normalize */
466 if (nstaken) /* multi-sampling */
467 rayclear(&sr);
468 rayvalue(&sr);
469 smultscolor(sr.rcol, sr.rcoef);
470 saddscolor(np->rp->rcol, sr.rcol);
471 ++nstaken;
472 }
473 ndims--;
474 }
475 }