ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/m_wgmdf.c
Revision: 2.3
Committed: Wed Dec 11 17:34:03 2024 UTC (4 months, 2 weeks ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +2 -1 lines
Log Message:
fix: Was over-estimating diffuse reflection with Fresnel effect

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.3 static const char RCSid[] = "$Id: m_wgmdf.c,v 2.2 2024/12/10 03:16: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     /* assign indicated diffuse component (do !trans first) */
162     static void
163     set_dcomp(WGMDDAT *wp, int trans)
164     {
165     DCOMP *dp = trans ? &wp->td : &wp->rd;
166     const int offs = trans ? 6 : 3*(wp->rp->rod < 0);
167    
168     if (trans) { /* transmitted diffuse? */
169     if (intens(wp->mtp->oargs.farg+offs) <= FTINY) {
170     scolorblack(dp->scol);
171     return;
172     }
173     dp->mo.nam = wp->mtp->oargs.sarg[8];
174     if (!fill_modval(&dp->mo, wp)) {
175     sprintf(errmsg,
176     "unknown diffuse transmission modifier '%s'",
177     dp->mo.nam);
178     objerror(wp->mtp, USER, errmsg);
179     }
180     } else /* no priors for main mod */
181     fill_modval(&dp->mo, wp);
182    
183     setscolor(dp->scol, wp->mtp->oargs.farg[offs],
184     wp->mtp->oargs.farg[offs+1],
185     wp->mtp->oargs.farg[offs+2]);
186     smultscolor(dp->scol, dp->mo.pcol);
187     }
188    
189     /* assign indicated specular component */
190     static void
191     set_scomp(WGMDDAT *wp, int trans)
192     {
193     SCOMP *sp = trans ? &wp->ts : &wp->rs;
194     const int eoff = 3*(trans != 0);
195     double coef;
196    
197     setfunc(wp->mtp, wp->rp); /* get coefficient, first */
198     errno = 0;
199     coef = evalue(wp->mf->ep[eoff]);
200     if ((errno == EDOM) | (errno == ERANGE)) {
201     objerror(wp->mtp, WARNING, "specular compute error");
202     scolorblack(sp->scol);
203     return;
204     }
205     if (coef <= FTINY) { /* negligible value? */
206     scolorblack(sp->scol);
207     return;
208     } /* else get modifier */
209     sp->mo.nam = wp->mtp->oargs.sarg[4*(trans != 0)];
210     if (!fill_modval(&sp->mo, wp)) {
211     sprintf(errmsg, "unknown specular %s modifier '%s'",
212     trans ? "transmission" : "reflection", sp->mo.nam);
213     objerror(wp->mtp, USER, errmsg);
214     }
215     copyscolor(sp->scol, sp->mo.pcol);
216     scalescolor(sp->scol, coef);
217     if (sintens(sp->scol) <= FTINY) {
218     scolorblack(sp->scol);
219     return; /* got black pattern */
220     }
221     setfunc(wp->mtp, wp->rp); /* else get roughness */
222     errno = 0;
223     sp->u_alpha = evalue(wp->mf->ep[eoff+1]);
224     sp->v_alpha = (sp->u_alpha > FTINY) ? evalue(wp->mf->ep[eoff+2]) : 0.0;
225     if ((errno == EDOM) | (errno == ERANGE)) {
226     objerror(wp->mtp, WARNING, "roughness compute error");
227     scolorblack(sp->scol);
228     return;
229     } /* we have something... */
230     wp->specfl |= trans ? SP_TRAN : SP_REFL;
231     if (sp->v_alpha <= FTINY) { /* is it pure specular? */
232     wp->specfl |= trans ? SP_TPURE : SP_RPURE;
233     sp->u_alpha = sp->v_alpha = 0.0;
234     return;
235     }
236     /* get anisotropic coordinates */
237     fcross(sp->v, sp->mo.pnorm, wp->ulocal);
238     if (normalize(sp->v) == 0.0) { /* orientation vector==normal? */
239     if (fabs(sp->u_alpha - sp->v_alpha) > 0.001)
240     objerror(wp->mtp, WARNING, "bad orientation vector");
241     getperpendicular(sp->u, sp->mo.pnorm, 1); /* punting */
242     fcross(sp->v, sp->mo.pnorm, sp->u);
243     sp->u_alpha = sp->v_alpha = sqrt( 0.5 *
244     (sp->u_alpha*sp->u_alpha + sp->v_alpha*sp->v_alpha) );
245     } else
246     fcross(sp->u, sp->v, sp->mo.pnorm);
247     }
248    
249     /* sample anisotropic Gaussian specular */
250     static void
251     agaussamp(WGMDDAT *wp)
252     {
253     RAY sr;
254     FVECT h;
255     double rv[2];
256     double d, sinp, cosp;
257     int maxiter, ntrials, nstarget, nstaken;
258     int i;
259     /* compute reflection */
260     if ((wp->specfl & (SP_REFL|SP_RPURE|SP_RBLT)) == SP_REFL &&
261     rayorigin(&sr, RSPECULAR, wp->rp, wp->rs.scol) == 0) {
262     SCOLOR scol;
263     nstarget = 1;
264     if (specjitter > 1.5) { /* multiple samples? */
265     nstarget = specjitter*wp->rp->rweight + .5;
266     if (sr.rweight <= minweight*nstarget)
267     nstarget = sr.rweight/minweight;
268     if (nstarget > 1) {
269     d = 1./nstarget;
270     scalescolor(sr.rcoef, d);
271     sr.rweight *= d;
272     } else
273     nstarget = 1;
274     }
275     scolorblack(scol);
276     dimlist[ndims++] = (int)(size_t)wp->mtp;
277     maxiter = MAXITER*nstarget;
278     for (nstaken = ntrials = 0; (nstaken < nstarget) &
279     (ntrials < maxiter); ntrials++) {
280     if (ntrials)
281     d = frandom();
282     else
283     d = urand(ilhash(dimlist,ndims)+samplendx);
284     multisamp(rv, 2, d);
285     d = 2.0*PI * rv[0];
286     cosp = tcos(d) * wp->rs.u_alpha;
287     sinp = tsin(d) * wp->rs.v_alpha;
288     d = 1./sqrt(cosp*cosp + sinp*sinp);
289     cosp *= d;
290     sinp *= d;
291     if ((0. <= specjitter) & (specjitter < 1.))
292     rv[1] = 1.0 - specjitter*rv[1];
293     d = (rv[1] <= FTINY) ? 1.0 : sqrt( -log(rv[1]) /
294     (cosp*cosp/(wp->rs.u_alpha*wp->rs.u_alpha) +
295     sinp*sinp/(wp->rs.v_alpha*wp->rs.v_alpha)) );
296     for (i = 0; i < 3; i++)
297     h[i] = wp->rs.mo.pnorm[i] +
298     d*(cosp*wp->rs.u[i] + sinp*wp->rs.v[i]);
299     d = -2.0 * DOT(h, wp->rp->rdir) / (1.0 + d*d);
300     VSUM(sr.rdir, wp->rp->rdir, h, d);
301     /* sample rejection test */
302     d = DOT(sr.rdir, wp->rp->ron);
303     if ((d > 0) ^ (wp->rp->rod > 0))
304     continue;
305     checknorm(sr.rdir);
306     if (nstarget > 1) { /* W-G-M-D adjustment */
307     if (nstaken) rayclear(&sr);
308     rayvalue(&sr);
309     d = 2./(1. + wp->rp->rod/d);
310     scalescolor(sr.rcol, d);
311     saddscolor(scol, sr.rcol);
312     } else {
313     rayvalue(&sr);
314     smultscolor(sr.rcol, sr.rcoef);
315     saddscolor(wp->rp->rcol, sr.rcol);
316     }
317     ++nstaken;
318     }
319     if (nstarget > 1) { /* final W-G-M-D weighting */
320     smultscolor(scol, sr.rcoef);
321     d = (double)nstarget/ntrials;
322     scalescolor(scol, d);
323     saddscolor(wp->rp->rcol, scol);
324     }
325     ndims--;
326     }
327     /* compute transmission */
328     if ((wp->specfl & (SP_TRAN|SP_TPURE|SP_TBLT)) == SP_TRAN &&
329     rayorigin(&sr, TSPECULAR, wp->rp, wp->ts.scol) == 0) {
330     nstarget = 1;
331     if (specjitter > 1.5) { /* multiple samples? */
332     nstarget = specjitter*wp->rp->rweight + .5;
333     if (sr.rweight <= minweight*nstarget)
334     nstarget = sr.rweight/minweight;
335     if (nstarget > 1) {
336     d = 1./nstarget;
337     scalescolor(sr.rcoef, d);
338     sr.rweight *= d;
339     } else
340     nstarget = 1;
341     }
342     dimlist[ndims++] = (int)(size_t)wp->mtp;
343     maxiter = MAXITER*nstarget;
344     for (nstaken = ntrials = 0; (nstaken < nstarget) &
345     (ntrials < maxiter); ntrials++) {
346     if (ntrials)
347     d = frandom();
348     else
349     d = urand(ilhash(dimlist,ndims)+1823+samplendx);
350     multisamp(rv, 2, d);
351     d = 2.0*PI * rv[0];
352     cosp = tcos(d) * wp->ts.u_alpha;
353     sinp = tsin(d) * wp->ts.v_alpha;
354     d = 1./sqrt(cosp*cosp + sinp*sinp);
355     cosp *= d;
356     sinp *= d;
357     if ((0. <= specjitter) & (specjitter < 1.))
358     rv[1] = 1.0 - specjitter*rv[1];
359     if (rv[1] <= FTINY)
360     d = 1.0;
361     else
362     d = sqrt(-log(rv[1]) /
363     (cosp*cosp/(wp->ts.u_alpha*wp->ts.u_alpha) +
364     sinp*sinp/(wp->ts.v_alpha*wp->ts.v_alpha)));
365     for (i = 0; i < 3; i++)
366     sr.rdir[i] = wp->prdir[i] +
367     d*(cosp*wp->ts.u[i] + sinp*wp->ts.v[i]);
368     /* rejection test */
369     if ((DOT(sr.rdir,wp->rp->ron) > 0) == (wp->rp->rod > 0))
370     continue;
371     normalize(sr.rdir); /* OK, normalize */
372     if (nstaken) /* multi-sampling? */
373     rayclear(&sr);
374     rayvalue(&sr);
375     smultscolor(sr.rcol, sr.rcoef);
376     saddscolor(wp->rp->rcol, sr.rcol);
377     ++nstaken;
378     }
379     ndims--;
380     }
381     }
382    
383     /* compute source contribution for MAT_WGMDF */
384     static void
385     dirwgmdf(SCOLOR scval, void *uwp, FVECT ldir, double omega)
386     {
387     WGMDDAT *wp = (WGMDDAT *)uwp;
388     const int hitfront = (wp->rp->rod > 0);
389     double fresadj = 1.;
390     double ldot;
391     double dtmp, dtmp1, dtmp2;
392     FVECT h;
393     double au2, av2;
394     SCOLOR sctmp;
395    
396     scolorblack(scval); /* will add component coefficients */
397    
398     /* XXX ignores which side is lit */
399     if (wp->specfl & SP_RPURE && pbright(wp->rs.scol) >= FRESTHRESH)
400     fresadj = 1. - FRESNE(fabs(DOT(wp->rs.mo.pnorm,ldir)));
401    
402     if (sintens(wp->rd.scol) > FTINY &&
403     ((ldot = DOT(wp->rd.mo.pnorm,ldir)) > 0) == hitfront) {
404     /*
405     * Compute diffuse reflection coefficient for source.
406     */
407     copyscolor(sctmp, wp->rd.scol);
408     dtmp = fabs(ldot) * omega * (1.0/PI) * fresadj;
409     scalescolor(sctmp, dtmp);
410     saddscolor(scval, sctmp);
411     }
412     if (sintens(wp->td.scol) > FTINY &&
413     ((ldot = DOT(wp->td.mo.pnorm,ldir)) > 0) ^ hitfront) {
414     /*
415     * Compute diffuse transmission coefficient for source.
416     */
417     copyscolor(sctmp, wp->td.scol);
418     dtmp = fabs(ldot) * omega * (1.0/PI) * fresadj;
419     scalescolor(sctmp, dtmp);
420     saddscolor(scval, sctmp);
421     }
422     #if 0 /* XXX not yet implemented */
423     if (ambRayInPmap(wp->rp))
424     return; /* specular accounted for in photon map */
425     #endif
426     if ((wp->specfl & (SP_REFL|SP_RPURE)) == SP_REFL &&
427     ((ldot = DOT(wp->rs.mo.pnorm,ldir)) > 0) == hitfront) {
428     /*
429     * Compute specular reflection coefficient for source using
430     * anisotropic Gaussian distribution model.
431     */
432     /* add source width if flat */
433     if (wp->specfl & SP_FLAT)
434     au2 = av2 = omega * (0.25/PI);
435     else
436     au2 = av2 = 0.0;
437     au2 += wp->rs.u_alpha*wp->rs.u_alpha;
438     av2 += wp->rs.v_alpha*wp->rs.v_alpha;
439     /* half vector */
440     VSUB(h, ldir, wp->rp->rdir);
441     /* ellipse */
442     dtmp1 = DOT(wp->rs.u, h);
443     dtmp1 *= dtmp1 / au2;
444     dtmp2 = DOT(wp->rs.v, h);
445     dtmp2 *= dtmp2 / av2;
446     /* W-G-M-D model */
447     dtmp = DOT(wp->rs.mo.pnorm, h);
448     dtmp *= dtmp;
449     dtmp1 = (dtmp1 + dtmp2) / dtmp;
450     dtmp = exp(-dtmp1) * DOT(h,h) /
451     (PI * dtmp*dtmp * sqrt(au2*av2));
452    
453     if (dtmp > FTINY) { /* worth using? */
454     copyscolor(sctmp, wp->rs.scol);
455     dtmp *= fabs(ldot) * omega;
456     scalescolor(sctmp, dtmp);
457     saddscolor(scval, sctmp);
458     }
459     }
460     if ((wp->specfl & (SP_TRAN|SP_TPURE)) == SP_TRAN &&
461     ((ldot = DOT(wp->ts.mo.pnorm,ldir)) > 0) ^ hitfront) {
462     /*
463     * Compute specular transmission coefficient for source.
464     */
465     /* roughness + source */
466     au2 = av2 = omega * (1.0/PI);
467     au2 += wp->ts.u_alpha*wp->ts.u_alpha;
468     av2 += wp->ts.v_alpha*wp->ts.v_alpha;
469     /* "half vector" */
470     VSUB(h, ldir, wp->prdir);
471     dtmp = DOT(h,h);
472     if (dtmp > FTINY*FTINY) {
473     dtmp1 = DOT(h,wp->ts.mo.pnorm);
474     dtmp = 1.0 - dtmp1*dtmp1/dtmp;
475     }
476     if (dtmp > FTINY*FTINY) {
477     dtmp1 = DOT(h,wp->ts.u);
478     dtmp1 *= dtmp1 / au2;
479     dtmp2 = DOT(h,wp->ts.v);
480     dtmp2 *= dtmp2 / av2;
481     dtmp = (dtmp1 + dtmp2) / dtmp;
482     dtmp = exp(-dtmp);
483     } else
484     dtmp = 1.0;
485     /* Gaussian */
486     dtmp *= (1.0/PI) * sqrt(-ldot/(wp->ts.mo.pdot*au2*av2));
487    
488     if (dtmp > FTINY) { /* worth using? */
489     copyscolor(sctmp, wp->ts.scol);
490     dtmp *= omega;
491     scalescolor(sctmp, dtmp);
492     saddscolor(scval, sctmp);
493     }
494     }
495     }
496    
497     /* color a ray that hit a programmable WGMD material */
498     int
499     m_wgmdf(OBJREC *m, RAY *r)
500     {
501     RAY lr;
502     WGMDDAT wd;
503     SCOLOR sctmp;
504     FVECT anorm;
505     int i;
506    
507     if (!backvis & (r->rod < 0.0)) {
508     raytrans(r);
509     return(1); /* backside invisible */
510     }
511     if ((m->oargs.nsargs < 13) | (m->oargs.nfargs < 9))
512     objerror(m, USER, "bad number of arguments");
513 greg 2.2
514     if (r->crtype & SHADOW && !strcmp(m->oargs.sarg[5], "0"))
515     return(1); /* first shadow test */
516 greg 2.1 clr_comps(&wd);
517     wd.rp = r;
518     wd.mtp = m;
519     wd.mf = getfunc(m, 12, 0xEEE, 1);
520 greg 2.2 set_dcomp(&wd, 0); /* calls main modifier */
521 greg 2.1 setfunc(m, r); /* get local u vector */
522     errno = 0;
523     for (i = 0; i < 3; i++)
524     wd.ulocal[i] = evalue(wd.mf->ep[6+i]);
525     if ((errno == EDOM) | (errno == ERANGE))
526     wd.ulocal[0] = wd.ulocal[1] = wd.ulocal[2] = 0.0;
527     else if (wd.mf->fxp != &unitxf)
528     multv3(wd.ulocal, wd.ulocal, wd.mf->fxp->xfm);
529    
530     set_scomp(&wd, 1); /* sets SP_TPURE */
531     if (r->crtype & SHADOW && !(wd.specfl & SP_TPURE))
532 greg 2.2 return(1); /* second shadow test */
533     set_dcomp(&wd, 1);
534 greg 2.1 set_scomp(&wd, 0);
535     wd.specfl |= SP_FLAT*(r->ro != NULL && isflat(r->ro->otype));
536     /* apply Fresnel adjustments? */
537     if (wd.specfl & SP_RPURE && pbright(wd.rs.scol) >= FRESTHRESH) {
538     const double fest = FRESNE(fabs(wd.rs.mo.pdot));
539     for (i = NCSAMP; i--; )
540     wd.rs.scol[i] += fest*(1. - wd.rs.scol[i]);
541 greg 2.3 scalecolor(wd.rd.scol, 1.-fest);
542 greg 2.1 scalescolor(wd.ts.scol, 1.-fest);
543     scalescolor(wd.td.scol, 1.-fest);
544     }
545     /* check specular thresholds */
546     wd.specfl |= SP_RBLT*((wd.specfl & (SP_REFL|SP_RPURE)) == SP_REFL &&
547     specthresh >= pbright(wd.rs.scol)-FTINY);
548     wd.specfl |= SP_TBLT*((wd.specfl & (SP_TRAN|SP_TPURE)) == SP_TRAN &&
549     specthresh >= pbright(wd.ts.scol)-FTINY);
550     /* get through direction */
551     if (wd.specfl & SP_TRAN && wd.ts.mo.hastexture &&
552     !(r->crtype & (SHADOW|AMBIENT))) {
553     for (i = 0; i < 3; i++) /* perturb */
554     wd.prdir[i] = r->rdir[i] - wd.ts.mo.pnorm[i] + r->ron[i];
555     if ((DOT(wd.prdir,r->ron) > 0) ^ (r->rod > 0))
556     normalize(wd.prdir); /* OK */
557     else /* too much */
558     VCOPY(wd.prdir, r->rdir);
559     } else
560     VCOPY(wd.prdir, r->rdir);
561     /* transmitted view ray? */
562     if ((wd.specfl & (SP_TRAN|SP_TPURE|SP_TBLT)) == (SP_TRAN|SP_TPURE) &&
563     rayorigin(&lr, TRANS, r, wd.ts.scol) == 0) {
564     VCOPY(lr.rdir, wd.prdir);
565     rayvalue(&lr);
566     smultscolor(lr.rcol, lr.rcoef);
567     saddscolor(r->rcol, lr.rcol);
568     if (scolor_mean(wd.ts.scol) >= 0.999) {
569     /* completely transparent */
570     smultscolor(lr.mcol, lr.rcoef);
571     copyscolor(r->mcol, lr.mcol);
572     r->rmt = r->rot + lr.rmt;
573     r->rxt = r->rot + lr.rxt;
574     } else if (pbright(wd.ts.scol) >
575     pbright(wd.td.scol) + pbright(wd.rd.scol))
576     r->rxt = r->rot + raydistance(&lr);
577     }
578     if (r->crtype & SHADOW)
579     return(1); /* the rest is shadow */
580     /* mirror ray? */
581     if ((wd.specfl & (SP_REFL|SP_RPURE|SP_RBLT)) == (SP_REFL|SP_RPURE) &&
582     rayorigin(&lr, REFLECTED, r, wd.rs.scol) == 0) {
583     VSUM(lr.rdir, r->rdir, wd.rs.mo.pnorm, 2.*wd.rs.mo.pdot);
584     /* fall back if would penetrate */
585     if (wd.rs.mo.hastexture &&
586     (DOT(lr.rdir,r->ron) > 0) ^ (r->rod > 0))
587     VSUM(lr.rdir, r->rdir, r->ron, 2.*r->rod);
588     checknorm(lr.rdir);
589     rayvalue(&lr);
590     smultscolor(lr.rcol, lr.rcoef);
591     copyscolor(r->mcol, lr.rcol);
592     saddscolor(r->rcol, lr.rcol);
593     r->rmt = r->rot;
594     if (wd.specfl & SP_FLAT &&
595     !wd.rs.mo.hastexture | (r->crtype & AMBIENT))
596     r->rmt += raydistance(&lr);
597     }
598     if (wd.specfl & (SP_REFL|SP_TRAN)) /* specularly scattered rays */
599     agaussamp(&wd); /* checks *BLT flags */
600    
601     if (sintens(wd.rd.scol) > FTINY) { /* ambient from this side */
602     if (r->rod > 0) {
603     VCOPY(anorm, wd.rd.mo.pnorm);
604     } else {
605     anorm[0] = -wd.rd.mo.pnorm[0];
606     anorm[1] = -wd.rd.mo.pnorm[1];
607     anorm[2] = -wd.rd.mo.pnorm[2];
608     }
609     copyscolor(sctmp, wd.rd.scol);
610     if (wd.specfl & SP_RBLT) /* add in specular as well? */
611     saddscolor(sctmp, wd.rs.scol);
612     multambient(sctmp, r, anorm);
613     saddscolor(r->rcol, sctmp); /* add to returned color */
614     }
615     if (sintens(wd.td.scol) > FTINY) { /* ambient from other side */
616     if (r->rod > 0) {
617     anorm[0] = -wd.td.mo.pnorm[0];
618     anorm[1] = -wd.td.mo.pnorm[1];
619     anorm[2] = -wd.td.mo.pnorm[2];
620     } else {
621     VCOPY(anorm, wd.td.mo.pnorm);
622     }
623     copyscolor(sctmp, wd.td.scol);
624     if (wd.specfl & SP_TBLT) /* add in specular as well? */
625     saddscolor(sctmp, wd.ts.scol)
626     multambient(sctmp, r, anorm);
627     saddscolor(r->rcol, sctmp);
628     }
629     direct(r, dirwgmdf, &wd); /* add direct component last */
630     return(1);
631     }