ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/m_wgmdf.c
Revision: 2.6
Committed: Wed Dec 18 17:57:06 2024 UTC (4 months, 2 weeks ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.5: +14 -15 lines
Log Message:
perf: Minor tweaks

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.6 static const char RCSid[] = "$Id: m_wgmdf.c,v 2.5 2024/12/17 20:03:13 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * Shading function for programmable Ward-Geisler-Moroder-Duer material.
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     /* estimate of Fresnel function */
23     #define FRESNE(ci) (exp(-5.85*(ci)) - 0.00202943064)
24     #define FRESTHRESH 0.017999 /* minimum specularity for approx. */
25    
26     /*
27     * This routine implements the anisotropic Gaussian
28     * model described by Ward in a 1992 Siggraph article and updated by
29     * Geisler-Moroder and Duer in a 2010 article in High Performance Graphics.
30     * We do not reorient incoming ray, using side in part to determine
31     * reflectance values. Most parameters are programmable with their own
32     * modifiers and/or value expressions.
33     *
34     * Arguments for MAT_WGMDF are:
35     * 13+ rs_mod rs rs_urough rs_vrough
36     * ts_mod ts ts_urough ts_vrough
37     * td_mod
38     * ux uy uz funcfile transform
39     * 0
40     * 9+ rfdif gfdif bfdif
41     * rbdif gbdif bbdif
42     * rtdif gtdif btdif
43     * A10 ..
44     *
45     * Where the rs_urough or rs_vrough expression yields zero, mirror-Fresnel
46     * effects are computed, similar to MAT_PLASTIC and MAT_METAL. The
47     * rs* expressions should not vary with incident angle, or the material
48     * will not be physically valid. Similarly, the ts* expressions should
49     * give the same value for coincident direction vectors from either side.
50     * There are independent modifiers for specular reflection,
51     * transmission, and diffuse transmission. Diffuse reflection
52     * applies the material's main modifier, which doesn't apply to
53     * anything else by default. However, any of the modifiers may be
54     * ALIASMOD, which will use the main material modifier, or VOIDID,
55     * which will just be white.
56     * Diffuse reflection and transmission colors and patterns add to
57     * the specular components, and are only adjusted with mirror-Fresnel
58     * reflection if specular reflection is greater than FRESHTHRESH. The
59     * specular transmission is likewise adjusted in such cases. Specified
60     * values for all components should sum to less than 1, but like other
61     * Radiance materials, this is not enforced, nor is a warning issued.
62     */
63     /* specularity flags */
64     #define SP_REFL 01 /* has reflected specular component */
65     #define SP_TRAN 02 /* has transmitted specular */
66     #define SP_RPURE 04 /* mirror reflection */
67     #define SP_TPURE 010 /* has view component */
68     #define SP_FLAT 020 /* flat reflecting surface */
69     #define SP_RBLT 040 /* reflection below sample threshold */
70     #define SP_TBLT 0100 /* transmission below threshold */
71    
72     typedef struct {
73     char *nam; /* modifier name */
74     int hastexture; /* has a texture? */
75     FVECT pnorm; /* perturbed normal direction */
76     double pdot; /* perturbed dot product */
77     SCOLOR pcol; /* pattern color */
78     } MODVAL; /* modifier-derived values */
79    
80     typedef struct {
81     MODVAL mo; /* modifier parameters */
82     SCOLOR scol; /* modified diffuse color */
83     } DCOMP; /* diffuse component parameters */
84    
85     typedef struct {
86     MODVAL mo; /* modifier parameters */
87     SCOLOR scol; /* modified specular color */
88     FVECT u, v; /* u and v in-plane vectors */
89     double u_alpha; /* u roughness */
90     double v_alpha; /* v roughness */
91     } SCOMP; /* specular component parameters */
92    
93     typedef struct {
94     RAY *rp; /* ray pointer */
95     OBJREC *mtp; /* material pointer */
96     MFUNC *mf; /* pointer to expression list */
97     int specfl; /* specularity flags, defined above */
98     FVECT ulocal; /* u-vector in local coordinates */
99     DCOMP rd, td; /* diffuse component params */
100     SCOMP rs, ts; /* specular component params */
101     FVECT prdir; /* vector in transmitted direction */
102     } WGMDDAT; /* WGMD material data */
103    
104     #define clr_comps(wp) ((wp)->specfl = 0, \
105     (wp)->rd.mo.nam = (wp)->td.mo.nam = \
106     (wp)->rs.mo.nam = (wp)->ts.mo.nam = "")
107    
108     /* assign modifier values */
109     static int
110     set_modval(MODVAL *mp, OBJECT omod, const RAY *r)
111     {
112     RAY tr;
113    
114     if (!mp->nam[0])
115     mp->nam = (omod == OVOID) ? VOIDID : objptr(omod)->oname;
116     else if (!strcmp(mp->nam, VOIDID))
117     omod = OVOID;
118     else if (omod == OVOID)
119     return(0);
120     tr = *r; /* independent modifier */
121     raytexture(&tr, omod);
122     if (DOT(tr.pert,tr.pert) > FTINY*FTINY) {
123     mp->pdot = raynormal(mp->pnorm, &tr);
124     mp->hastexture = 1;
125     } else {
126     VCOPY(mp->pnorm, tr.ron);
127     mp->pdot = tr.rod;
128     mp->hastexture = 0;
129     }
130     copyscolor(mp->pcol, tr.pcol);
131     return(1);
132     }
133    
134     /* fill modifier values, using previous setting if found */
135     static int
136     fill_modval(MODVAL *mp, const WGMDDAT *wp)
137     {
138     if (mp == &wp->rd.mo) { /* special case (should be first) */
139     set_modval(mp, wp->mtp->omod, wp->rp);
140     return(1);
141     } /* use main modifier? */
142     if (!strcmp(mp->nam, ALIASMOD) || !strcmp(mp->nam, wp->rd.mo.nam)) {
143     *mp = wp->rd.mo;
144     return(1);
145     } /* check others */
146     if (mp != &wp->td.mo && !strcmp(mp->nam, wp->td.mo.nam)) {
147     *mp = wp->td.mo;
148     return(1);
149     }
150     if (mp != &wp->rs.mo && !strcmp(mp->nam, wp->rs.mo.nam)) {
151     *mp = wp->rs.mo;
152     return(1);
153     }
154     if (mp != &wp->ts.mo && !strcmp(mp->nam, wp->ts.mo.nam)) {
155     *mp = wp->ts.mo;
156     return(1);
157     } /* new modifier */
158     return(set_modval(mp, lastmod(objndx(wp->mtp), mp->nam), wp->rp));
159     }
160    
161 greg 2.6 /* set calculation context for given component of MAT_WGMDF */
162 greg 2.5 static int
163     setWGMDfunc(MODVAL *mp, const WGMDDAT *wp)
164     {
165 greg 2.6 static char lastMod[MAXSTR];
166 greg 2.5 double sf;
167     FVECT vec;
168    
169     if (setfunc(wp->mtp, wp->rp) == 0 &&
170     !strcmp(mp->nam, lastMod))
171     return(0); /* already set */
172 greg 2.6 strcpy(lastMod, mp->nam);
173 greg 2.5 /* else (re)assign special variables */
174 greg 2.6 sf = 1 - 2*(wp->rp->rod < 0);
175     varset("RdotP`", '=', mp->pdot*sf);
176     multv3(vec, mp->pnorm, funcxf.xfm);
177 greg 2.5 sf /= funcxf.sca;
178     varset("NxP`", '=', vec[0]*sf);
179     varset("NyP`", '=', vec[1]*sf);
180     varset("NzP`", '=', vec[2]*sf);
181     return(1);
182     }
183    
184 greg 2.1 /* assign indicated diffuse component (do !trans first) */
185     static void
186     set_dcomp(WGMDDAT *wp, int trans)
187     {
188     DCOMP *dp = trans ? &wp->td : &wp->rd;
189     const int offs = trans ? 6 : 3*(wp->rp->rod < 0);
190    
191     if (trans) { /* transmitted diffuse? */
192     if (intens(wp->mtp->oargs.farg+offs) <= FTINY) {
193     scolorblack(dp->scol);
194     return;
195     }
196     dp->mo.nam = wp->mtp->oargs.sarg[8];
197     if (!fill_modval(&dp->mo, wp)) {
198     sprintf(errmsg,
199     "unknown diffuse transmission modifier '%s'",
200     dp->mo.nam);
201     objerror(wp->mtp, USER, errmsg);
202     }
203     } else /* no priors for main mod */
204     fill_modval(&dp->mo, wp);
205    
206     setscolor(dp->scol, wp->mtp->oargs.farg[offs],
207     wp->mtp->oargs.farg[offs+1],
208     wp->mtp->oargs.farg[offs+2]);
209     smultscolor(dp->scol, dp->mo.pcol);
210     }
211    
212     /* assign indicated specular component */
213     static void
214     set_scomp(WGMDDAT *wp, int trans)
215     {
216 greg 2.6 SCOMP *sp = trans ? &wp->ts : &wp->rs;
217     EPNODE **exa = wp->mf->ep + 3*(trans != 0);
218     double coef;
219 greg 2.5 /* constant zero check */
220 greg 2.6 if (exa[0]->type == NUM && exa[0]->v.num <= FTINY) {
221 greg 2.5 scolorblack(sp->scol);
222     return;
223     } /* need modifier */
224     sp->mo.nam = wp->mtp->oargs.sarg[4*(trans != 0)];
225     if (!fill_modval(&sp->mo, wp)) {
226     sprintf(errmsg, "unknown specular %s modifier '%s'",
227     trans ? "transmission" : "reflection", sp->mo.nam);
228     objerror(wp->mtp, USER, errmsg);
229     }
230     setWGMDfunc(&sp->mo, wp);
231 greg 2.1 errno = 0;
232 greg 2.6 coef = evalue(exa[0]);
233 greg 2.1 if ((errno == EDOM) | (errno == ERANGE)) {
234     objerror(wp->mtp, WARNING, "specular compute error");
235     scolorblack(sp->scol);
236     return;
237     }
238     if (coef <= FTINY) { /* negligible value? */
239     scolorblack(sp->scol);
240     return;
241     }
242     copyscolor(sp->scol, sp->mo.pcol);
243     scalescolor(sp->scol, coef);
244     if (sintens(sp->scol) <= FTINY) {
245     scolorblack(sp->scol);
246     return; /* got black pattern */
247     }
248 greg 2.5 errno = 0; /* else get roughness */
249 greg 2.6 sp->u_alpha = evalue(exa[1]);
250     sp->v_alpha = (sp->u_alpha > FTINY) ? evalue(exa[2]) : 0.0;
251 greg 2.1 if ((errno == EDOM) | (errno == ERANGE)) {
252     objerror(wp->mtp, WARNING, "roughness compute error");
253     scolorblack(sp->scol);
254     return;
255     } /* we have something... */
256     wp->specfl |= trans ? SP_TRAN : SP_REFL;
257     if (sp->v_alpha <= FTINY) { /* is it pure specular? */
258     wp->specfl |= trans ? SP_TPURE : SP_RPURE;
259     sp->u_alpha = sp->v_alpha = 0.0;
260     return;
261     }
262     /* get anisotropic coordinates */
263     fcross(sp->v, sp->mo.pnorm, wp->ulocal);
264     if (normalize(sp->v) == 0.0) { /* orientation vector==normal? */
265     if (fabs(sp->u_alpha - sp->v_alpha) > 0.001)
266     objerror(wp->mtp, WARNING, "bad orientation vector");
267     getperpendicular(sp->u, sp->mo.pnorm, 1); /* punting */
268     fcross(sp->v, sp->mo.pnorm, sp->u);
269     sp->u_alpha = sp->v_alpha = sqrt( 0.5 *
270     (sp->u_alpha*sp->u_alpha + sp->v_alpha*sp->v_alpha) );
271     } else
272     fcross(sp->u, sp->v, sp->mo.pnorm);
273     }
274    
275     /* sample anisotropic Gaussian specular */
276     static void
277     agaussamp(WGMDDAT *wp)
278     {
279     RAY sr;
280     FVECT h;
281     double rv[2];
282     double d, sinp, cosp;
283     int maxiter, ntrials, nstarget, nstaken;
284     int i;
285     /* compute reflection */
286     if ((wp->specfl & (SP_REFL|SP_RPURE|SP_RBLT)) == SP_REFL &&
287     rayorigin(&sr, RSPECULAR, wp->rp, wp->rs.scol) == 0) {
288     SCOLOR scol;
289     nstarget = 1;
290     if (specjitter > 1.5) { /* multiple samples? */
291     nstarget = specjitter*wp->rp->rweight + .5;
292     if (sr.rweight <= minweight*nstarget)
293     nstarget = sr.rweight/minweight;
294     if (nstarget > 1) {
295     d = 1./nstarget;
296     scalescolor(sr.rcoef, d);
297     sr.rweight *= d;
298     } else
299     nstarget = 1;
300     }
301     scolorblack(scol);
302     dimlist[ndims++] = (int)(size_t)wp->mtp;
303     maxiter = MAXITER*nstarget;
304     for (nstaken = ntrials = 0; (nstaken < nstarget) &
305     (ntrials < maxiter); ntrials++) {
306     if (ntrials)
307     d = frandom();
308     else
309     d = urand(ilhash(dimlist,ndims)+samplendx);
310     multisamp(rv, 2, d);
311     d = 2.0*PI * rv[0];
312     cosp = tcos(d) * wp->rs.u_alpha;
313     sinp = tsin(d) * wp->rs.v_alpha;
314     d = 1./sqrt(cosp*cosp + sinp*sinp);
315     cosp *= d;
316     sinp *= d;
317     if ((0. <= specjitter) & (specjitter < 1.))
318     rv[1] = 1.0 - specjitter*rv[1];
319     d = (rv[1] <= FTINY) ? 1.0 : sqrt( -log(rv[1]) /
320     (cosp*cosp/(wp->rs.u_alpha*wp->rs.u_alpha) +
321     sinp*sinp/(wp->rs.v_alpha*wp->rs.v_alpha)) );
322     for (i = 0; i < 3; i++)
323     h[i] = wp->rs.mo.pnorm[i] +
324     d*(cosp*wp->rs.u[i] + sinp*wp->rs.v[i]);
325     d = -2.0 * DOT(h, wp->rp->rdir) / (1.0 + d*d);
326     VSUM(sr.rdir, wp->rp->rdir, h, d);
327     /* sample rejection test */
328     d = DOT(sr.rdir, wp->rp->ron);
329     if ((d > 0) ^ (wp->rp->rod > 0))
330     continue;
331     checknorm(sr.rdir);
332     if (nstarget > 1) { /* W-G-M-D adjustment */
333     if (nstaken) rayclear(&sr);
334     rayvalue(&sr);
335     d = 2./(1. + wp->rp->rod/d);
336     scalescolor(sr.rcol, d);
337     saddscolor(scol, sr.rcol);
338     } else {
339     rayvalue(&sr);
340     smultscolor(sr.rcol, sr.rcoef);
341     saddscolor(wp->rp->rcol, sr.rcol);
342     }
343     ++nstaken;
344     }
345     if (nstarget > 1) { /* final W-G-M-D weighting */
346     smultscolor(scol, sr.rcoef);
347     d = (double)nstarget/ntrials;
348     scalescolor(scol, d);
349     saddscolor(wp->rp->rcol, scol);
350     }
351     ndims--;
352     }
353     /* compute transmission */
354     if ((wp->specfl & (SP_TRAN|SP_TPURE|SP_TBLT)) == SP_TRAN &&
355     rayorigin(&sr, TSPECULAR, wp->rp, wp->ts.scol) == 0) {
356     nstarget = 1;
357     if (specjitter > 1.5) { /* multiple samples? */
358     nstarget = specjitter*wp->rp->rweight + .5;
359     if (sr.rweight <= minweight*nstarget)
360     nstarget = sr.rweight/minweight;
361     if (nstarget > 1) {
362     d = 1./nstarget;
363     scalescolor(sr.rcoef, d);
364     sr.rweight *= d;
365     } else
366     nstarget = 1;
367     }
368     dimlist[ndims++] = (int)(size_t)wp->mtp;
369     maxiter = MAXITER*nstarget;
370     for (nstaken = ntrials = 0; (nstaken < nstarget) &
371     (ntrials < maxiter); ntrials++) {
372     if (ntrials)
373     d = frandom();
374     else
375     d = urand(ilhash(dimlist,ndims)+1823+samplendx);
376     multisamp(rv, 2, d);
377     d = 2.0*PI * rv[0];
378     cosp = tcos(d) * wp->ts.u_alpha;
379     sinp = tsin(d) * wp->ts.v_alpha;
380     d = 1./sqrt(cosp*cosp + sinp*sinp);
381     cosp *= d;
382     sinp *= d;
383     if ((0. <= specjitter) & (specjitter < 1.))
384     rv[1] = 1.0 - specjitter*rv[1];
385     if (rv[1] <= FTINY)
386     d = 1.0;
387     else
388     d = sqrt(-log(rv[1]) /
389     (cosp*cosp/(wp->ts.u_alpha*wp->ts.u_alpha) +
390     sinp*sinp/(wp->ts.v_alpha*wp->ts.v_alpha)));
391     for (i = 0; i < 3; i++)
392     sr.rdir[i] = wp->prdir[i] +
393     d*(cosp*wp->ts.u[i] + sinp*wp->ts.v[i]);
394     /* rejection test */
395     if ((DOT(sr.rdir,wp->rp->ron) > 0) == (wp->rp->rod > 0))
396     continue;
397     normalize(sr.rdir); /* OK, normalize */
398     if (nstaken) /* multi-sampling? */
399     rayclear(&sr);
400     rayvalue(&sr);
401     smultscolor(sr.rcol, sr.rcoef);
402     saddscolor(wp->rp->rcol, sr.rcol);
403     ++nstaken;
404     }
405     ndims--;
406     }
407     }
408    
409     /* compute source contribution for MAT_WGMDF */
410     static void
411     dirwgmdf(SCOLOR scval, void *uwp, FVECT ldir, double omega)
412     {
413     WGMDDAT *wp = (WGMDDAT *)uwp;
414     const int hitfront = (wp->rp->rod > 0);
415     double fresadj = 1.;
416     double ldot;
417     double dtmp, dtmp1, dtmp2;
418     FVECT h;
419     double au2, av2;
420     SCOLOR sctmp;
421    
422     scolorblack(scval); /* will add component coefficients */
423    
424     /* XXX ignores which side is lit */
425     if (wp->specfl & SP_RPURE && pbright(wp->rs.scol) >= FRESTHRESH)
426     fresadj = 1. - FRESNE(fabs(DOT(wp->rs.mo.pnorm,ldir)));
427    
428     if (sintens(wp->rd.scol) > FTINY &&
429     ((ldot = DOT(wp->rd.mo.pnorm,ldir)) > 0) == hitfront) {
430     /*
431     * Compute diffuse reflection coefficient for source.
432     */
433     copyscolor(sctmp, wp->rd.scol);
434     dtmp = fabs(ldot) * omega * (1.0/PI) * fresadj;
435     scalescolor(sctmp, dtmp);
436     saddscolor(scval, sctmp);
437     }
438     if (sintens(wp->td.scol) > FTINY &&
439     ((ldot = DOT(wp->td.mo.pnorm,ldir)) > 0) ^ hitfront) {
440     /*
441     * Compute diffuse transmission coefficient for source.
442     */
443     copyscolor(sctmp, wp->td.scol);
444     dtmp = fabs(ldot) * omega * (1.0/PI) * fresadj;
445     scalescolor(sctmp, dtmp);
446     saddscolor(scval, sctmp);
447     }
448     #if 0 /* XXX not yet implemented */
449     if (ambRayInPmap(wp->rp))
450     return; /* specular accounted for in photon map */
451     #endif
452     if ((wp->specfl & (SP_REFL|SP_RPURE)) == SP_REFL &&
453     ((ldot = DOT(wp->rs.mo.pnorm,ldir)) > 0) == hitfront) {
454     /*
455     * Compute specular reflection coefficient for source using
456     * anisotropic Gaussian distribution model.
457     */
458     /* add source width if flat */
459     if (wp->specfl & SP_FLAT)
460     au2 = av2 = omega * (0.25/PI);
461     else
462     au2 = av2 = 0.0;
463     au2 += wp->rs.u_alpha*wp->rs.u_alpha;
464     av2 += wp->rs.v_alpha*wp->rs.v_alpha;
465     /* half vector */
466     VSUB(h, ldir, wp->rp->rdir);
467     /* ellipse */
468     dtmp1 = DOT(wp->rs.u, h);
469     dtmp1 *= dtmp1 / au2;
470     dtmp2 = DOT(wp->rs.v, h);
471     dtmp2 *= dtmp2 / av2;
472     /* W-G-M-D model */
473     dtmp = DOT(wp->rs.mo.pnorm, h);
474     dtmp *= dtmp;
475     dtmp1 = (dtmp1 + dtmp2) / dtmp;
476     dtmp = exp(-dtmp1) * DOT(h,h) /
477     (PI * dtmp*dtmp * sqrt(au2*av2));
478    
479     if (dtmp > FTINY) { /* worth using? */
480     copyscolor(sctmp, wp->rs.scol);
481     dtmp *= fabs(ldot) * omega;
482     scalescolor(sctmp, dtmp);
483     saddscolor(scval, sctmp);
484     }
485     }
486     if ((wp->specfl & (SP_TRAN|SP_TPURE)) == SP_TRAN &&
487     ((ldot = DOT(wp->ts.mo.pnorm,ldir)) > 0) ^ hitfront) {
488     /*
489     * Compute specular transmission coefficient for source.
490     */
491     /* roughness + source */
492     au2 = av2 = omega * (1.0/PI);
493     au2 += wp->ts.u_alpha*wp->ts.u_alpha;
494     av2 += wp->ts.v_alpha*wp->ts.v_alpha;
495     /* "half vector" */
496     VSUB(h, ldir, wp->prdir);
497     dtmp = DOT(h,h);
498     if (dtmp > FTINY*FTINY) {
499     dtmp1 = DOT(h,wp->ts.mo.pnorm);
500     dtmp = 1.0 - dtmp1*dtmp1/dtmp;
501     }
502     if (dtmp > FTINY*FTINY) {
503     dtmp1 = DOT(h,wp->ts.u);
504     dtmp1 *= dtmp1 / au2;
505     dtmp2 = DOT(h,wp->ts.v);
506     dtmp2 *= dtmp2 / av2;
507     dtmp = (dtmp1 + dtmp2) / dtmp;
508     dtmp = exp(-dtmp);
509     } else
510     dtmp = 1.0;
511     /* Gaussian */
512     dtmp *= (1.0/PI) * sqrt(-ldot/(wp->ts.mo.pdot*au2*av2));
513    
514     if (dtmp > FTINY) { /* worth using? */
515     copyscolor(sctmp, wp->ts.scol);
516     dtmp *= omega;
517     scalescolor(sctmp, dtmp);
518     saddscolor(scval, sctmp);
519     }
520     }
521     }
522    
523     /* color a ray that hit a programmable WGMD material */
524     int
525     m_wgmdf(OBJREC *m, RAY *r)
526     {
527     RAY lr;
528     WGMDDAT wd;
529     SCOLOR sctmp;
530     FVECT anorm;
531     int i;
532    
533     if (!backvis & (r->rod < 0.0)) {
534     raytrans(r);
535     return(1); /* backside invisible */
536     }
537     if ((m->oargs.nsargs < 13) | (m->oargs.nfargs < 9))
538     objerror(m, USER, "bad number of arguments");
539 greg 2.2
540     if (r->crtype & SHADOW && !strcmp(m->oargs.sarg[5], "0"))
541     return(1); /* first shadow test */
542 greg 2.1 clr_comps(&wd);
543     wd.rp = r;
544     wd.mtp = m;
545     wd.mf = getfunc(m, 12, 0xEEE, 1);
546 greg 2.5 set_dcomp(&wd, 0); /* gets main modifier */
547     setWGMDfunc(&wd.rd.mo, &wd); /* get local u vector */
548 greg 2.1 errno = 0;
549     for (i = 0; i < 3; i++)
550     wd.ulocal[i] = evalue(wd.mf->ep[6+i]);
551     if ((errno == EDOM) | (errno == ERANGE))
552     wd.ulocal[0] = wd.ulocal[1] = wd.ulocal[2] = 0.0;
553     else if (wd.mf->fxp != &unitxf)
554     multv3(wd.ulocal, wd.ulocal, wd.mf->fxp->xfm);
555    
556     set_scomp(&wd, 1); /* sets SP_TPURE */
557     if (r->crtype & SHADOW && !(wd.specfl & SP_TPURE))
558 greg 2.2 return(1); /* second shadow test */
559     set_dcomp(&wd, 1);
560 greg 2.1 set_scomp(&wd, 0);
561     wd.specfl |= SP_FLAT*(r->ro != NULL && isflat(r->ro->otype));
562     /* apply Fresnel adjustments? */
563     if (wd.specfl & SP_RPURE && pbright(wd.rs.scol) >= FRESTHRESH) {
564     const double fest = FRESNE(fabs(wd.rs.mo.pdot));
565     for (i = NCSAMP; i--; )
566     wd.rs.scol[i] += fest*(1. - wd.rs.scol[i]);
567 greg 2.4 scalescolor(wd.rd.scol, 1.-fest);
568 greg 2.1 scalescolor(wd.ts.scol, 1.-fest);
569     scalescolor(wd.td.scol, 1.-fest);
570     }
571     /* check specular thresholds */
572     wd.specfl |= SP_RBLT*((wd.specfl & (SP_REFL|SP_RPURE)) == SP_REFL &&
573     specthresh >= pbright(wd.rs.scol)-FTINY);
574     wd.specfl |= SP_TBLT*((wd.specfl & (SP_TRAN|SP_TPURE)) == SP_TRAN &&
575     specthresh >= pbright(wd.ts.scol)-FTINY);
576     /* get through direction */
577     if (wd.specfl & SP_TRAN && wd.ts.mo.hastexture &&
578     !(r->crtype & (SHADOW|AMBIENT))) {
579     for (i = 0; i < 3; i++) /* perturb */
580     wd.prdir[i] = r->rdir[i] - wd.ts.mo.pnorm[i] + r->ron[i];
581     if ((DOT(wd.prdir,r->ron) > 0) ^ (r->rod > 0))
582     normalize(wd.prdir); /* OK */
583     else /* too much */
584     VCOPY(wd.prdir, r->rdir);
585     } else
586     VCOPY(wd.prdir, r->rdir);
587     /* transmitted view ray? */
588     if ((wd.specfl & (SP_TRAN|SP_TPURE|SP_TBLT)) == (SP_TRAN|SP_TPURE) &&
589     rayorigin(&lr, TRANS, r, wd.ts.scol) == 0) {
590     VCOPY(lr.rdir, wd.prdir);
591     rayvalue(&lr);
592     smultscolor(lr.rcol, lr.rcoef);
593     saddscolor(r->rcol, lr.rcol);
594     if (scolor_mean(wd.ts.scol) >= 0.999) {
595     /* completely transparent */
596     smultscolor(lr.mcol, lr.rcoef);
597     copyscolor(r->mcol, lr.mcol);
598     r->rmt = r->rot + lr.rmt;
599     r->rxt = r->rot + lr.rxt;
600     } else if (pbright(wd.ts.scol) >
601     pbright(wd.td.scol) + pbright(wd.rd.scol))
602     r->rxt = r->rot + raydistance(&lr);
603     }
604     if (r->crtype & SHADOW)
605     return(1); /* the rest is shadow */
606     /* mirror ray? */
607     if ((wd.specfl & (SP_REFL|SP_RPURE|SP_RBLT)) == (SP_REFL|SP_RPURE) &&
608     rayorigin(&lr, REFLECTED, r, wd.rs.scol) == 0) {
609     VSUM(lr.rdir, r->rdir, wd.rs.mo.pnorm, 2.*wd.rs.mo.pdot);
610     /* fall back if would penetrate */
611     if (wd.rs.mo.hastexture &&
612     (DOT(lr.rdir,r->ron) > 0) ^ (r->rod > 0))
613     VSUM(lr.rdir, r->rdir, r->ron, 2.*r->rod);
614     checknorm(lr.rdir);
615     rayvalue(&lr);
616     smultscolor(lr.rcol, lr.rcoef);
617     copyscolor(r->mcol, lr.rcol);
618     saddscolor(r->rcol, lr.rcol);
619     r->rmt = r->rot;
620     if (wd.specfl & SP_FLAT &&
621     !wd.rs.mo.hastexture | (r->crtype & AMBIENT))
622     r->rmt += raydistance(&lr);
623     }
624     if (wd.specfl & (SP_REFL|SP_TRAN)) /* specularly scattered rays */
625     agaussamp(&wd); /* checks *BLT flags */
626    
627     if (sintens(wd.rd.scol) > FTINY) { /* ambient from this side */
628     if (r->rod > 0) {
629     VCOPY(anorm, wd.rd.mo.pnorm);
630     } else {
631     anorm[0] = -wd.rd.mo.pnorm[0];
632     anorm[1] = -wd.rd.mo.pnorm[1];
633     anorm[2] = -wd.rd.mo.pnorm[2];
634     }
635     copyscolor(sctmp, wd.rd.scol);
636     if (wd.specfl & SP_RBLT) /* add in specular as well? */
637     saddscolor(sctmp, wd.rs.scol);
638     multambient(sctmp, r, anorm);
639     saddscolor(r->rcol, sctmp); /* add to returned color */
640     }
641     if (sintens(wd.td.scol) > FTINY) { /* ambient from other side */
642     if (r->rod > 0) {
643     anorm[0] = -wd.td.mo.pnorm[0];
644     anorm[1] = -wd.td.mo.pnorm[1];
645     anorm[2] = -wd.td.mo.pnorm[2];
646     } else {
647     VCOPY(anorm, wd.td.mo.pnorm);
648     }
649     copyscolor(sctmp, wd.td.scol);
650     if (wd.specfl & SP_TBLT) /* add in specular as well? */
651     saddscolor(sctmp, wd.ts.scol)
652     multambient(sctmp, r, anorm);
653     saddscolor(r->rcol, sctmp);
654     }
655     direct(r, dirwgmdf, &wd); /* add direct component last */
656     return(1);
657     }