src/rt/normal.c

#ifndef lint
static const char RCSid[] = "$Id: normal.c,v 2.63 2012/07/29 21:56:16 greg Exp $";
#endif
/*
 *  normal.c - shading function for normal materials.
 *
 *     8/19/85
 *     12/19/85 - added stuff for metals.
 *     6/26/87 - improved specular model.
 *     9/28/87 - added model for translucent materials.
 *     Later changes described in delta comments.
 */

#include "copyright.h"

#include  "ray.h"
#include  "ambient.h"
#include  "source.h"
#include  "otypes.h"
#include  "rtotypes.h"
#include  "random.h"

#ifndef  MAXITER
#define  MAXITER        10              /* maximum # specular ray attempts */
#endif
                                        /* estimate of Fresnel function */
#define  FRESNE(ci)     (exp(-5.85*(ci)) - 0.00287989916)
#define  FRESTHRESH     0.017999        /* minimum specularity for approx. */


/*
 *      This routine implements the isotropic Gaussian
 *  model described by Ward in Siggraph `92 article.
 *      We orient the surface towards the incoming ray, so a single
 *  surface can be used to represent an infinitely thin object.
 *
 *  Arguments for MAT_PLASTIC and MAT_METAL are:
 *      red     grn     blu     specular-frac.  facet-slope
 *
 *  Arguments for MAT_TRANS are:
 *      red     grn     blu     rspec   rough   trans   tspec
 */

                                /* specularity flags */
#define  SP_REFL        01              /* has reflected specular component */
#define  SP_TRAN        02              /* has transmitted specular */
#define  SP_PURE        04              /* purely specular (zero roughness) */
#define  SP_FLAT        010             /* flat reflecting surface */
#define  SP_RBLT        020             /* reflection below sample threshold */
#define  SP_TBLT        040             /* transmission below threshold */

typedef struct {
        OBJREC  *mp;            /* material pointer */
        RAY  *rp;               /* ray pointer */
        short  specfl;          /* specularity flags, defined above */
        COLOR  mcolor;          /* color of this material */
        COLOR  scolor;          /* color of specular component */
        FVECT  vrefl;           /* vector in direction of reflected ray */
        FVECT  prdir;           /* vector in transmitted direction */
        double  alpha2;         /* roughness squared */
        double  rdiff, rspec;   /* reflected specular, diffuse */
        double  trans;          /* transmissivity */
        double  tdiff, tspec;   /* transmitted specular, diffuse */
        FVECT  pnorm;           /* perturbed surface normal */
        double  pdot;           /* perturbed dot product */
}  NORMDAT;             /* normal material data */

static void gaussamp(NORMDAT  *np);


static void
dirnorm(                /* compute source contribution */
        COLOR  cval,                    /* returned coefficient */
        void  *nnp,                     /* material data */
        FVECT  ldir,                    /* light source direction */
        double  omega                   /* light source size */
)
{
        NORMDAT *np = nnp;
        double  ldot;
        double  lrdiff, ltdiff;
        double  dtmp, d2, d3, d4;
        FVECT  vtmp;
        COLOR  ctmp;

        setcolor(cval, 0.0, 0.0, 0.0);

        ldot = DOT(np->pnorm, ldir);

        if (ldot < 0.0 ? np->trans <= FTINY : np->trans >= 1.0-FTINY)
                return;         /* wrong side */

                                /* Fresnel estimate */
        lrdiff = np->rdiff;
        ltdiff = np->tdiff;
        if (np->specfl & SP_PURE && np->rspec >= FRESTHRESH &&
                        (lrdiff > FTINY) | (ltdiff > FTINY)) {
                dtmp = 1. - FRESNE(fabs(ldot));
                lrdiff *= dtmp;
                ltdiff *= dtmp;
        }

        if (ldot > FTINY && lrdiff > FTINY) {
                /*
                 *  Compute and add diffuse reflected component to returned
                 *  color.  The diffuse reflected component will always be
                 *  modified by the color of the material.
                 */
                copycolor(ctmp, np->mcolor);
                dtmp = ldot * omega * lrdiff * (1.0/PI);
                scalecolor(ctmp, dtmp);
                addcolor(cval, ctmp);
        }
        if (ldot > FTINY && (np->specfl&(SP_REFL|SP_PURE)) == SP_REFL) {
                /*
                 *  Compute specular reflection coefficient using
                 *  Gaussian distribution model.
                 */
                                                /* roughness */
                dtmp = np->alpha2;
                                                /* + source if flat */
                if (np->specfl & SP_FLAT)
                        dtmp += omega * (0.25/PI);
                                                /* half vector */
                VSUB(vtmp, ldir, np->rp->rdir);
                d2 = DOT(vtmp, np->pnorm);
                d2 *= d2;
                d3 = DOT(vtmp,vtmp);
                d4 = (d3 - d2) / d2;
                                                /* new W-G-M-D model */
                dtmp = exp(-d4/dtmp) * d3 / (PI * d2*d2 * dtmp);
                                                /* worth using? */
                if (dtmp > FTINY) {
                        copycolor(ctmp, np->scolor);
                        dtmp *= ldot * omega;
                        scalecolor(ctmp, dtmp);
                        addcolor(cval, ctmp);
                }
        }
        if (ldot < -FTINY && ltdiff > FTINY) {
                /*
                 *  Compute diffuse transmission.
                 */
                copycolor(ctmp, np->mcolor);
                dtmp = -ldot * omega * ltdiff * (1.0/PI);
                scalecolor(ctmp, dtmp);
                addcolor(cval, ctmp);
        }
        if (ldot < -FTINY && (np->specfl&(SP_TRAN|SP_PURE)) == SP_TRAN) {
                /*
                 *  Compute specular transmission.  Specular transmission
                 *  is always modified by material color.
                 */
                                                /* roughness + source */
                dtmp = np->alpha2 + omega*(1.0/PI);
                                                /* Gaussian */
                dtmp = exp((2.*DOT(np->prdir,ldir)-2.)/dtmp)/(PI*dtmp);
                                                /* worth using? */
                if (dtmp > FTINY) {
                        copycolor(ctmp, np->mcolor);
                        dtmp *= np->tspec * omega * sqrt(-ldot/np->pdot);
                        scalecolor(ctmp, dtmp);
                        addcolor(cval, ctmp);
                }
        }
}


int
m_normal(                       /* color a ray that hit something normal */
        OBJREC  *m,
        RAY  *r
)
{
        NORMDAT  nd;
        double  fest;
        double  transtest, transdist;
        double  mirtest, mirdist;
        int     hastexture;
        double  d;
        COLOR  ctmp;
        int  i;
                                                /* easy shadow test */
        if (r->crtype & SHADOW && m->otype != MAT_TRANS)
                return(1);

        if (m->oargs.nfargs != (m->otype == MAT_TRANS ? 7 : 5))
                objerror(m, USER, "bad number of arguments");
                                                /* check for back side */
        if (r->rod < 0.0) {
                if (!backvis && m->otype != MAT_TRANS) {
                        raytrans(r);
                        return(1);
                }
                raytexture(r, m->omod);
                flipsurface(r);                 /* reorient if backvis */
        } else
                raytexture(r, m->omod);
        nd.mp = m;
        nd.rp = r;
                                                /* get material color */
        setcolor(nd.mcolor, m->oargs.farg[0],
                           m->oargs.farg[1],
                           m->oargs.farg[2]);
                                                /* get roughness */
        nd.specfl = 0;
        nd.alpha2 = m->oargs.farg[4];
        if ((nd.alpha2 *= nd.alpha2) <= FTINY)
                nd.specfl |= SP_PURE;

        if ( (hastexture = (DOT(r->pert,r->pert) > FTINY*FTINY)) ) {
                nd.pdot = raynormal(nd.pnorm, r);       /* perturb normal */
        } else {
                VCOPY(nd.pnorm, r->ron);
                nd.pdot = r->rod;
        }
        if (r->ro != NULL && isflat(r->ro->otype))
                nd.specfl |= SP_FLAT;
        if (nd.pdot < .001)
                nd.pdot = .001;                 /* non-zero for dirnorm() */
        multcolor(nd.mcolor, r->pcol);          /* modify material color */
        mirtest = transtest = 0;
        mirdist = transdist = r->rot;
        nd.rspec = m->oargs.farg[3];
                                                /* compute Fresnel approx. */
        if (nd.specfl & SP_PURE && nd.rspec >= FRESTHRESH) {
                fest = FRESNE(nd.pdot);
                nd.rspec += fest*(1. - nd.rspec);
        } else
                fest = 0.;
                                                /* compute transmission */
        if (m->otype == MAT_TRANS) {
                nd.trans = m->oargs.farg[5]*(1.0 - nd.rspec);
                nd.tspec = nd.trans * m->oargs.farg[6];
                nd.tdiff = nd.trans - nd.tspec;
                if (nd.tspec > FTINY) {
                        nd.specfl |= SP_TRAN;
                                                        /* check threshold */
                        if (!(nd.specfl & SP_PURE) &&
                                        specthresh >= nd.tspec-FTINY)
                                nd.specfl |= SP_TBLT;
                        if (!hastexture || r->crtype & SHADOW) {
                                VCOPY(nd.prdir, r->rdir);
                                transtest = 2;
                        } else {
                                for (i = 0; i < 3; i++)         /* perturb */
                                        nd.prdir[i] = r->rdir[i] - r->pert[i];
                                if (DOT(nd.prdir, r->ron) < -FTINY)
                                        normalize(nd.prdir);    /* OK */
                                else
                                        VCOPY(nd.prdir, r->rdir);
                        }
                }
        } else
                nd.tdiff = nd.tspec = nd.trans = 0.0;
                                                /* transmitted ray */
        if ((nd.specfl&(SP_TRAN|SP_PURE|SP_TBLT)) == (SP_TRAN|SP_PURE)) {
                RAY  lr;
                copycolor(lr.rcoef, nd.mcolor); /* modified by color */
                scalecolor(lr.rcoef, nd.tspec);
                if (rayorigin(&lr, TRANS, r, lr.rcoef) == 0) {
                        VCOPY(lr.rdir, nd.prdir);
                        rayvalue(&lr);
                        multcolor(lr.rcol, lr.rcoef);
                        addcolor(r->rcol, lr.rcol);
                        transtest *= bright(lr.rcol);
                        transdist = r->rot + lr.rt;
                }
        } else
                transtest = 0;

        if (r->crtype & SHADOW) {               /* the rest is shadow */
                r->rt = transdist;
                return(1);
        }
                                                /* get specular reflection */
        if (nd.rspec > FTINY) {
                nd.specfl |= SP_REFL;
                                                /* compute specular color */
                if (m->otype != MAT_METAL) {
                        setcolor(nd.scolor, nd.rspec, nd.rspec, nd.rspec);
                } else if (fest > FTINY) {
                        d = m->oargs.farg[3]*(1. - fest);
                        for (i = 0; i < 3; i++)
                                nd.scolor[i] = fest + nd.mcolor[i]*d;
                } else {
                        copycolor(nd.scolor, nd.mcolor);
                        scalecolor(nd.scolor, nd.rspec);
                }
                                                /* check threshold */
                if (!(nd.specfl & SP_PURE) && specthresh >= nd.rspec-FTINY)
                        nd.specfl |= SP_RBLT;
                                                /* compute reflected ray */
                VSUM(nd.vrefl, r->rdir, nd.pnorm, 2.*nd.pdot);
                                                /* penetration? */
                if (hastexture && DOT(nd.vrefl, r->ron) <= FTINY)
                        VSUM(nd.vrefl, r->rdir, r->ron, 2.*r->rod);
                checknorm(nd.vrefl);
        }
                                                /* reflected ray */
        if ((nd.specfl&(SP_REFL|SP_PURE|SP_RBLT)) == (SP_REFL|SP_PURE)) {
                RAY  lr;
                if (rayorigin(&lr, REFLECTED, r, nd.scolor) == 0) {
                        VCOPY(lr.rdir, nd.vrefl);
                        rayvalue(&lr);
                        multcolor(lr.rcol, lr.rcoef);
                        addcolor(r->rcol, lr.rcol);
                        if (!hastexture && nd.specfl & SP_FLAT) {
                                mirtest = 2.*bright(lr.rcol);
                                mirdist = r->rot + lr.rt;
                        }
                }
        }
                                                /* diffuse reflection */
        nd.rdiff = 1.0 - nd.trans - nd.rspec;

        if (nd.specfl & SP_PURE && nd.rdiff <= FTINY && nd.tdiff <= FTINY)
                return(1);                      /* 100% pure specular */

        if (!(nd.specfl & SP_PURE))
                gaussamp(&nd);          /* checks *BLT flags */

        if (nd.rdiff > FTINY) {         /* ambient from this side */
                copycolor(ctmp, nd.mcolor);     /* modified by material color */
                scalecolor(ctmp, nd.rdiff);
                if (nd.specfl & SP_RBLT)        /* add in specular as well? */
                        addcolor(ctmp, nd.scolor);
                multambient(ctmp, r, hastexture ? nd.pnorm : r->ron);
                addcolor(r->rcol, ctmp);        /* add to returned color */
        }
        if (nd.tdiff > FTINY) {         /* ambient from other side */
                copycolor(ctmp, nd.mcolor);     /* modified by color */
                if (nd.specfl & SP_TBLT)
                        scalecolor(ctmp, nd.trans);
                else
                        scalecolor(ctmp, nd.tdiff);
                flipsurface(r);
                if (hastexture) {
                        FVECT  bnorm;
                        bnorm[0] = -nd.pnorm[0];
                        bnorm[1] = -nd.pnorm[1];
                        bnorm[2] = -nd.pnorm[2];
                        multambient(ctmp, r, bnorm);
                } else
                        multambient(ctmp, r, r->ron);
                addcolor(r->rcol, ctmp);
                flipsurface(r);
        }
                                        /* add direct component */
        direct(r, dirnorm, &nd);
                                        /* check distance */
        d = bright(r->rcol);
        if (transtest > d)
                r->rt = transdist;
        else if (mirtest > d)
                r->rt = mirdist;

        return(1);
}


static void
gaussamp(                       /* sample Gaussian specular */
        NORMDAT  *np
)
{
        RAY  sr;
        FVECT  u, v, h;
        double  rv[2];
        double  d, sinp, cosp;
        COLOR  scol;
        int  maxiter, ntrials, nstarget, nstaken;
        int  i;
                                        /* quick test */
        if ((np->specfl & (SP_REFL|SP_RBLT)) != SP_REFL &&
                        (np->specfl & (SP_TRAN|SP_TBLT)) != SP_TRAN)
                return;
                                        /* set up sample coordinates */
        v[0] = v[1] = v[2] = 0.0;
        for (i = 0; i < 3; i++)
                if (np->pnorm[i] < 0.6 && np->pnorm[i] > -0.6)
                        break;
        v[i] = 1.0;
        fcross(u, v, np->pnorm);
        normalize(u);
        fcross(v, np->pnorm, u);
                                        /* compute reflection */
        if ((np->specfl & (SP_REFL|SP_RBLT)) == SP_REFL &&
                        rayorigin(&sr, SPECULAR, np->rp, np->scolor) == 0) {
                nstarget = 1;
                if (specjitter > 1.5) { /* multiple samples? */
                        nstarget = specjitter*np->rp->rweight + .5;
                        if (sr.rweight <= minweight*nstarget)
                                nstarget = sr.rweight/minweight;
                        if (nstarget > 1) {
                                d = 1./nstarget;
                                scalecolor(sr.rcoef, d);
                                sr.rweight *= d;
                        } else
                                nstarget = 1;
                }
                setcolor(scol, 0., 0., 0.);
                dimlist[ndims++] = (int)(size_t)np->mp;
                maxiter = MAXITER*nstarget;
                for (nstaken = ntrials = 0; nstaken < nstarget &&
                                                ntrials < maxiter; ntrials++) {
                        if (ntrials)
                                d = frandom();
                        else
                                d = urand(ilhash(dimlist,ndims)+samplendx);
                        multisamp(rv, 2, d);
                        d = 2.0*PI * rv[0];
                        cosp = tcos(d);
                        sinp = tsin(d);
                        if ((0. <= specjitter) & (specjitter < 1.))
                                rv[1] = 1.0 - specjitter*rv[1];
                        if (rv[1] <= FTINY)
                                d = 1.0;
                        else
                                d = sqrt( np->alpha2 * -log(rv[1]) );
                        for (i = 0; i < 3; i++)
                                h[i] = np->pnorm[i] + d*(cosp*u[i] + sinp*v[i]);
                        d = -2.0 * DOT(h, np->rp->rdir) / (1.0 + d*d);
                        VSUM(sr.rdir, np->rp->rdir, h, d);
                                                /* sample rejection test */
                        if ((d = DOT(sr.rdir, np->rp->ron)) <= FTINY)
                                continue;
                        checknorm(sr.rdir);
                        if (nstarget > 1) {     /* W-G-M-D adjustment */
                                if (nstaken) rayclear(&sr);
                                rayvalue(&sr);
                                d = 2./(1. + np->rp->rod/d);
                                scalecolor(sr.rcol, d);
                                addcolor(scol, sr.rcol);
                        } else {
                                rayvalue(&sr);
                                multcolor(sr.rcol, sr.rcoef);
                                addcolor(np->rp->rcol, sr.rcol);
                        }
                        ++nstaken;
                }
                if (nstarget > 1) {             /* final W-G-M-D weighting */
                        multcolor(scol, sr.rcoef);
                        d = (double)nstarget/ntrials;
                        scalecolor(scol, d);
                        addcolor(np->rp->rcol, scol);
                }
                ndims--;
        }
                                        /* compute transmission */
        copycolor(sr.rcoef, np->mcolor);        /* modified by color */
        scalecolor(sr.rcoef, np->tspec);
        if ((np->specfl & (SP_TRAN|SP_TBLT)) == SP_TRAN &&
                        rayorigin(&sr, SPECULAR, np->rp, sr.rcoef) == 0) {
                nstarget = 1;
                if (specjitter > 1.5) { /* multiple samples? */
                        nstarget = specjitter*np->rp->rweight + .5;
                        if (sr.rweight <= minweight*nstarget)
                                nstarget = sr.rweight/minweight;
                        if (nstarget > 1) {
                                d = 1./nstarget;
                                scalecolor(sr.rcoef, d);
                                sr.rweight *= d;
                        } else
                                nstarget = 1;
                }
                dimlist[ndims++] = (int)(size_t)np->mp;
                maxiter = MAXITER*nstarget;
                for (nstaken = ntrials = 0; nstaken < nstarget &&
                                                ntrials < maxiter; ntrials++) {
                        if (ntrials)
                                d = frandom();
                        else
                                d = urand(ilhash(dimlist,ndims)+samplendx);
                        multisamp(rv, 2, d);
                        d = 2.0*PI * rv[0];
                        cosp = tcos(d);
                        sinp = tsin(d);
                        if ((0. <= specjitter) & (specjitter < 1.))
                                rv[1] = 1.0 - specjitter*rv[1];
                        if (rv[1] <= FTINY)
                                d = 1.0;
                        else
                                d = sqrt( np->alpha2 * -log(rv[1]) );
                        for (i = 0; i < 3; i++)
                                sr.rdir[i] = np->prdir[i] + d*(cosp*u[i] + sinp*v[i]);
                                                /* sample rejection test */
                        if (DOT(sr.rdir, np->rp->ron) >= -FTINY)
                                continue;
                        normalize(sr.rdir);     /* OK, normalize */
                        if (nstaken)            /* multi-sampling */
                                rayclear(&sr);
                        rayvalue(&sr);
                        multcolor(sr.rcol, sr.rcoef);
                        addcolor(np->rp->rcol, sr.rcol);
                        ++nstaken;
                }
                ndims--;
        }
}
Revision:	2.64
Committed:	Mon Jul 30 17:46:50 2012 UTC (11 years, 9 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.63:	+2 -2 lines
Log Message:	Fixed bug in Fresnel approximation for metals
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.64	static const char RCSid[] = "$Id: normal.c,v 2.63 2012/07/29 21:56:16 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* normal.c - shading function for normal materials.
6			*
7			* 8/19/85
8			* 12/19/85 - added stuff for metals.
9			* 6/26/87 - improved specular model.
10			* 9/28/87 - added model for translucent materials.
11	greg	2.2	* Later changes described in delta comments.
12	greg	1.1	*/
13
14	greg	2.39	#include "copyright.h"
15	greg	2.38
16	greg	1.1	#include "ray.h"
17	greg	2.46	#include "ambient.h"
18	schorsch	2.47	#include "source.h"
19	greg	1.1	#include "otypes.h"
20	schorsch	2.47	#include "rtotypes.h"
21	greg	2.2	#include "random.h"
22
23	greg	2.34	#ifndef MAXITER
24			#define MAXITER 10 /* maximum # specular ray attempts */
25			#endif
26	greg	2.38	/* estimate of Fresnel function */
27	greg	2.44	#define FRESNE(ci) (exp(-5.85*(ci)) - 0.00287989916)
28	greg	2.51	#define FRESTHRESH 0.017999 /* minimum specularity for approx. */
29	greg	2.34
30	greg	2.24
31	greg	1.1	/*
32	greg	2.22	* This routine implements the isotropic Gaussian
33			* model described by Ward in Siggraph `92 article.
34	greg	1.1	* We orient the surface towards the incoming ray, so a single
35			* surface can be used to represent an infinitely thin object.
36			*
37			* Arguments for MAT_PLASTIC and MAT_METAL are:
38			* red grn blu specular-frac. facet-slope
39			*
40			* Arguments for MAT_TRANS are:
41			* red grn blu rspec rough trans tspec
42			*/
43
44	greg	2.2	/* specularity flags */
45			#define SP_REFL 01 /* has reflected specular component */
46			#define SP_TRAN 02 /* has transmitted specular */
47	greg	2.11	#define SP_PURE 04 /* purely specular (zero roughness) */
48			#define SP_FLAT 010 /* flat reflecting surface */
49			#define SP_RBLT 020 /* reflection below sample threshold */
50			#define SP_TBLT 040 /* transmission below threshold */
51	greg	1.1
52	greg	1.3	typedef struct {
53			OBJREC mp; / material pointer */
54	greg	2.16	RAY rp; / ray pointer */
55	greg	2.2	short specfl; /* specularity flags, defined above */
56	greg	1.1	COLOR mcolor; /* color of this material */
57			COLOR scolor; /* color of specular component */
58			FVECT vrefl; /* vector in direction of reflected ray */
59	greg	1.14	FVECT prdir; /* vector in transmitted direction */
60	greg	2.2	double alpha2; /* roughness squared */
61	greg	1.1	double rdiff, rspec; /* reflected specular, diffuse */
62			double trans; /* transmissivity */
63			double tdiff, tspec; /* transmitted specular, diffuse */
64			FVECT pnorm; /* perturbed surface normal */
65			double pdot; /* perturbed dot product */
66	greg	1.3	} NORMDAT; /* normal material data */
67
68	greg	2.63	static void gaussamp(NORMDAT *np);
69	schorsch	2.47
70	greg	1.3
71	greg	2.38	static void
72	schorsch	2.47	dirnorm( /* compute source contribution */
73			COLOR cval, /* returned coefficient */
74	greg	2.62	void nnp, / material data */
75	schorsch	2.47	FVECT ldir, /* light source direction */
76			double omega /* light source size */
77			)
78	greg	1.3	{
79	greg	2.62	NORMDAT *np = nnp;
80	greg	1.1	double ldot;
81	greg	2.49	double lrdiff, ltdiff;
82	greg	2.54	double dtmp, d2, d3, d4;
83	greg	2.16	FVECT vtmp;
84	greg	1.3	COLOR ctmp;
85
86			setcolor(cval, 0.0, 0.0, 0.0);
87
88			ldot = DOT(np->pnorm, ldir);
89
90			if (ldot < 0.0 ? np->trans <= FTINY : np->trans >= 1.0-FTINY)
91			return; /* wrong side */
92
93	greg	2.38	/* Fresnel estimate */
94	greg	2.49	lrdiff = np->rdiff;
95			ltdiff = np->tdiff;
96	greg	2.51	if (np->specfl & SP_PURE && np->rspec >= FRESTHRESH &&
97	greg	2.49	(lrdiff > FTINY) \| (ltdiff > FTINY)) {
98			dtmp = 1. - FRESNE(fabs(ldot));
99			lrdiff *= dtmp;
100			ltdiff *= dtmp;
101			}
102	greg	2.38
103	greg	2.49	if (ldot > FTINY && lrdiff > FTINY) {
104	greg	1.3	/*
105	greg	1.4	* Compute and add diffuse reflected component to returned
106			* color. The diffuse reflected component will always be
107			* modified by the color of the material.
108	greg	1.3	*/
109			copycolor(ctmp, np->mcolor);
110	greg	2.49	dtmp = ldot * omega * lrdiff * (1.0/PI);
111	greg	1.3	scalecolor(ctmp, dtmp);
112			addcolor(cval, ctmp);
113			}
114	greg	2.2	if (ldot > FTINY && (np->specfl&(SP_REFL\|SP_PURE)) == SP_REFL) {
115	greg	1.3	/*
116			* Compute specular reflection coefficient using
117	greg	2.54	* Gaussian distribution model.
118	greg	1.3	*/
119	greg	2.3	/* roughness */
120	greg	2.16	dtmp = np->alpha2;
121	greg	2.3	/* + source if flat */
122			if (np->specfl & SP_FLAT)
123	greg	2.48	dtmp += omega * (0.25/PI);
124	greg	2.23	/* half vector */
125	greg	2.62	VSUB(vtmp, ldir, np->rp->rdir);
126	greg	2.16	d2 = DOT(vtmp, np->pnorm);
127	greg	2.23	d2 *= d2;
128	greg	2.54	d3 = DOT(vtmp,vtmp);
129			d4 = (d3 - d2) / d2;
130			/* new W-G-M-D model */
131			dtmp = exp(-d4/dtmp) * d3 / (PI * d2d2 dtmp);
132	greg	1.3	/* worth using? */
133			if (dtmp > FTINY) {
134			copycolor(ctmp, np->scolor);
135	greg	2.54	dtmp = ldot omega;
136	greg	1.3	scalecolor(ctmp, dtmp);
137			addcolor(cval, ctmp);
138			}
139			}
140	greg	2.49	if (ldot < -FTINY && ltdiff > FTINY) {
141	greg	1.3	/*
142			* Compute diffuse transmission.
143			*/
144			copycolor(ctmp, np->mcolor);
145	greg	2.49	dtmp = -ldot * omega * ltdiff * (1.0/PI);
146	greg	1.3	scalecolor(ctmp, dtmp);
147			addcolor(cval, ctmp);
148			}
149	greg	2.2	if (ldot < -FTINY && (np->specfl&(SP_TRAN\|SP_PURE)) == SP_TRAN) {
150	greg	1.3	/*
151	greg	1.4	* Compute specular transmission. Specular transmission
152	greg	1.13	* is always modified by material color.
153	greg	1.3	*/
154			/* roughness + source */
155	greg	2.48	dtmp = np->alpha2 + omega*(1.0/PI);
156	greg	2.54	/* Gaussian */
157	greg	2.53	dtmp = exp((2.DOT(np->prdir,ldir)-2.)/dtmp)/(PIdtmp);
158	greg	1.3	/* worth using? */
159			if (dtmp > FTINY) {
160	greg	1.13	copycolor(ctmp, np->mcolor);
161	greg	2.52	dtmp = np->tspec omega * sqrt(-ldot/np->pdot);
162	greg	1.13	scalecolor(ctmp, dtmp);
163	greg	1.3	addcolor(cval, ctmp);
164			}
165			}
166			}
167
168
169	greg	2.62	int
170	schorsch	2.47	m_normal( /* color a ray that hit something normal */
171	greg	2.62	OBJREC *m,
172			RAY *r
173	schorsch	2.47	)
174	greg	1.3	{
175			NORMDAT nd;
176	greg	2.38	double fest;
177	greg	1.9	double transtest, transdist;
178	greg	2.29	double mirtest, mirdist;
179			int hastexture;
180			double d;
181	greg	1.1	COLOR ctmp;
182	greg	2.62	int i;
183	greg	1.1	/* easy shadow test */
184			if (r->crtype & SHADOW && m->otype != MAT_TRANS)
185	greg	2.27	return(1);
186	greg	2.2
187			if (m->oargs.nfargs != (m->otype == MAT_TRANS ? 7 : 5))
188			objerror(m, USER, "bad number of arguments");
189	greg	2.29	/* check for back side */
190			if (r->rod < 0.0) {
191			if (!backvis && m->otype != MAT_TRANS) {
192			raytrans(r);
193			return(1);
194			}
195	greg	2.40	raytexture(r, m->omod);
196	greg	2.29	flipsurface(r); /* reorient if backvis */
197	greg	2.40	} else
198			raytexture(r, m->omod);
199	greg	1.3	nd.mp = m;
200	greg	2.16	nd.rp = r;
201	greg	1.1	/* get material color */
202	greg	1.3	setcolor(nd.mcolor, m->oargs.farg[0],
203	greg	1.1	m->oargs.farg[1],
204			m->oargs.farg[2]);
205			/* get roughness */
206	greg	2.2	nd.specfl = 0;
207	greg	1.3	nd.alpha2 = m->oargs.farg[4];
208	greg	2.2	if ((nd.alpha2 *= nd.alpha2) <= FTINY)
209			nd.specfl \|= SP_PURE;
210	greg	2.40
211	schorsch	2.45	if ( (hastexture = (DOT(r->pert,r->pert) > FTINY*FTINY)) ) {
212	greg	2.29	nd.pdot = raynormal(nd.pnorm, r); /* perturb normal */
213	greg	2.41	} else {
214	greg	2.29	VCOPY(nd.pnorm, r->ron);
215			nd.pdot = r->rod;
216			}
217	greg	2.42	if (r->ro != NULL && isflat(r->ro->otype))
218			nd.specfl \|= SP_FLAT;
219	greg	1.13	if (nd.pdot < .001)
220			nd.pdot = .001; /* non-zero for dirnorm() */
221	greg	1.3	multcolor(nd.mcolor, r->pcol); /* modify material color */
222	greg	2.29	mirtest = transtest = 0;
223			mirdist = transdist = r->rot;
224	greg	2.30	nd.rspec = m->oargs.farg[3];
225	greg	2.38	/* compute Fresnel approx. */
226	greg	2.51	if (nd.specfl & SP_PURE && nd.rspec >= FRESTHRESH) {
227	greg	2.62	fest = FRESNE(nd.pdot);
228	greg	2.38	nd.rspec += fest*(1. - nd.rspec);
229			} else
230			fest = 0.;
231	greg	1.3	/* compute transmission */
232	greg	1.1	if (m->otype == MAT_TRANS) {
233	greg	1.3	nd.trans = m->oargs.farg[5]*(1.0 - nd.rspec);
234			nd.tspec = nd.trans * m->oargs.farg[6];
235			nd.tdiff = nd.trans - nd.tspec;
236	greg	2.2	if (nd.tspec > FTINY) {
237			nd.specfl \|= SP_TRAN;
238	greg	2.5	/* check threshold */
239	greg	2.25	if (!(nd.specfl & SP_PURE) &&
240			specthresh >= nd.tspec-FTINY)
241	greg	2.5	nd.specfl \|= SP_TBLT;
242	greg	2.29	if (!hastexture \|\| r->crtype & SHADOW) {
243	greg	2.2	VCOPY(nd.prdir, r->rdir);
244			transtest = 2;
245			} else {
246			for (i = 0; i < 3; i++) /* perturb */
247	greg	2.19	nd.prdir[i] = r->rdir[i] - r->pert[i];
248	greg	2.7	if (DOT(nd.prdir, r->ron) < -FTINY)
249			normalize(nd.prdir); /* OK */
250			else
251			VCOPY(nd.prdir, r->rdir);
252	greg	2.2	}
253	greg	1.14	}
254	greg	1.1	} else
255	greg	1.3	nd.tdiff = nd.tspec = nd.trans = 0.0;
256	greg	1.1	/* transmitted ray */
257	gregl	2.36	if ((nd.specfl&(SP_TRAN\|SP_PURE\|SP_TBLT)) == (SP_TRAN\|SP_PURE)) {
258	greg	1.3	RAY lr;
259	greg	2.50	copycolor(lr.rcoef, nd.mcolor); /* modified by color */
260			scalecolor(lr.rcoef, nd.tspec);
261			if (rayorigin(&lr, TRANS, r, lr.rcoef) == 0) {
262	greg	1.14	VCOPY(lr.rdir, nd.prdir);
263	greg	1.1	rayvalue(&lr);
264	greg	2.50	multcolor(lr.rcol, lr.rcoef);
265	greg	1.1	addcolor(r->rcol, lr.rcol);
266	greg	1.9	transtest *= bright(lr.rcol);
267			transdist = r->rot + lr.rt;
268	greg	1.1	}
269	greg	2.11	} else
270			transtest = 0;
271	greg	2.2
272	greg	2.29	if (r->crtype & SHADOW) { /* the rest is shadow */
273			r->rt = transdist;
274	greg	2.27	return(1);
275	greg	2.30	}
276			/* get specular reflection */
277			if (nd.rspec > FTINY) {
278			nd.specfl \|= SP_REFL;
279			/* compute specular color */
280	greg	2.38	if (m->otype != MAT_METAL) {
281			setcolor(nd.scolor, nd.rspec, nd.rspec, nd.rspec);
282			} else if (fest > FTINY) {
283	greg	2.64	d = m->oargs.farg[3]*(1. - fest);
284	greg	2.38	for (i = 0; i < 3; i++)
285			nd.scolor[i] = fest + nd.mcolor[i]*d;
286			} else {
287	greg	2.30	copycolor(nd.scolor, nd.mcolor);
288	greg	2.38	scalecolor(nd.scolor, nd.rspec);
289			}
290	greg	2.30	/* check threshold */
291			if (!(nd.specfl & SP_PURE) && specthresh >= nd.rspec-FTINY)
292			nd.specfl \|= SP_RBLT;
293			/* compute reflected ray */
294	greg	2.55	VSUM(nd.vrefl, r->rdir, nd.pnorm, 2.*nd.pdot);
295	greg	2.30	/* penetration? */
296			if (hastexture && DOT(nd.vrefl, r->ron) <= FTINY)
297	greg	2.55	VSUM(nd.vrefl, r->rdir, r->ron, 2.*r->rod);
298	greg	2.53	checknorm(nd.vrefl);
299	gregl	2.36	}
300			/* reflected ray */
301			if ((nd.specfl&(SP_REFL\|SP_PURE\|SP_RBLT)) == (SP_REFL\|SP_PURE)) {
302			RAY lr;
303	greg	2.50	if (rayorigin(&lr, REFLECTED, r, nd.scolor) == 0) {
304	gregl	2.36	VCOPY(lr.rdir, nd.vrefl);
305			rayvalue(&lr);
306	greg	2.50	multcolor(lr.rcol, lr.rcoef);
307	gregl	2.36	addcolor(r->rcol, lr.rcol);
308			if (!hastexture && nd.specfl & SP_FLAT) {
309			mirtest = 2.*bright(lr.rcol);
310			mirdist = r->rot + lr.rt;
311	greg	2.30	}
312			}
313	greg	2.29	}
314	greg	1.1	/* diffuse reflection */
315	greg	1.3	nd.rdiff = 1.0 - nd.trans - nd.rspec;
316	greg	1.1
317	greg	2.2	if (nd.specfl & SP_PURE && nd.rdiff <= FTINY && nd.tdiff <= FTINY)
318	greg	2.27	return(1); /* 100% pure specular */
319	greg	2.3
320	gregl	2.36	if (!(nd.specfl & SP_PURE))
321	greg	2.63	gaussamp(&nd); /* checks BLT flags /
322	greg	2.2
323	greg	1.3	if (nd.rdiff > FTINY) { /* ambient from this side */
324	greg	2.50	copycolor(ctmp, nd.mcolor); /* modified by material color */
325	greg	2.61	scalecolor(ctmp, nd.rdiff);
326			if (nd.specfl & SP_RBLT) /* add in specular as well? */
327			addcolor(ctmp, nd.scolor);
328	greg	2.50	multambient(ctmp, r, hastexture ? nd.pnorm : r->ron);
329	greg	1.2	addcolor(r->rcol, ctmp); /* add to returned color */
330			}
331	greg	1.3	if (nd.tdiff > FTINY) { /* ambient from other side */
332	greg	2.50	copycolor(ctmp, nd.mcolor); /* modified by color */
333			if (nd.specfl & SP_TBLT)
334			scalecolor(ctmp, nd.trans);
335			else
336			scalecolor(ctmp, nd.tdiff);
337	greg	1.1	flipsurface(r);
338	greg	2.32	if (hastexture) {
339			FVECT bnorm;
340			bnorm[0] = -nd.pnorm[0];
341			bnorm[1] = -nd.pnorm[1];
342			bnorm[2] = -nd.pnorm[2];
343	greg	2.50	multambient(ctmp, r, bnorm);
344	greg	2.32	} else
345	greg	2.50	multambient(ctmp, r, r->ron);
346	greg	1.1	addcolor(r->rcol, ctmp);
347			flipsurface(r);
348			}
349	greg	1.3	/* add direct component */
350			direct(r, dirnorm, &nd);
351	greg	1.9	/* check distance */
352	greg	2.29	d = bright(r->rcol);
353			if (transtest > d)
354	greg	1.9	r->rt = transdist;
355	greg	2.29	else if (mirtest > d)
356			r->rt = mirdist;
357	greg	2.27
358			return(1);
359	greg	2.2	}
360
361
362	greg	2.38	static void
363	greg	2.54	gaussamp( /* sample Gaussian specular */
364	greg	2.62	NORMDAT *np
365	schorsch	2.47	)
366	greg	2.2	{
367			RAY sr;
368			FVECT u, v, h;
369			double rv[2];
370			double d, sinp, cosp;
371	greg	2.62	COLOR scol;
372	greg	2.58	int maxiter, ntrials, nstarget, nstaken;
373	greg	2.62	int i;
374	greg	2.13	/* quick test */
375			if ((np->specfl & (SP_REFL\|SP_RBLT)) != SP_REFL &&
376			(np->specfl & (SP_TRAN\|SP_TBLT)) != SP_TRAN)
377			return;
378	greg	2.2	/* set up sample coordinates */
379			v[0] = v[1] = v[2] = 0.0;
380			for (i = 0; i < 3; i++)
381			if (np->pnorm[i] < 0.6 && np->pnorm[i] > -0.6)
382			break;
383			v[i] = 1.0;
384			fcross(u, v, np->pnorm);
385			normalize(u);
386			fcross(v, np->pnorm, u);
387			/* compute reflection */
388	greg	2.5	if ((np->specfl & (SP_REFL\|SP_RBLT)) == SP_REFL &&
389	greg	2.63	rayorigin(&sr, SPECULAR, np->rp, np->scolor) == 0) {
390	greg	2.58	nstarget = 1;
391	greg	2.55	if (specjitter > 1.5) { /* multiple samples? */
392	greg	2.63	nstarget = specjitter*np->rp->rweight + .5;
393	greg	2.58	if (sr.rweight <= minweight*nstarget)
394			nstarget = sr.rweight/minweight;
395			if (nstarget > 1) {
396			d = 1./nstarget;
397			scalecolor(sr.rcoef, d);
398	greg	2.56	sr.rweight *= d;
399	greg	2.55	} else
400	greg	2.58	nstarget = 1;
401	greg	2.55	}
402	greg	2.58	setcolor(scol, 0., 0., 0.);
403	greg	2.60	dimlist[ndims++] = (int)(size_t)np->mp;
404	greg	2.58	maxiter = MAXITER*nstarget;
405			for (nstaken = ntrials = 0; nstaken < nstarget &&
406			ntrials < maxiter; ntrials++) {
407			if (ntrials)
408	greg	2.34	d = frandom();
409			else
410			d = urand(ilhash(dimlist,ndims)+samplendx);
411			multisamp(rv, 2, d);
412			d = 2.0PI rv[0];
413	gwlarson	2.37	cosp = tcos(d);
414			sinp = tsin(d);
415	greg	2.55	if ((0. <= specjitter) & (specjitter < 1.))
416			rv[1] = 1.0 - specjitter*rv[1];
417	greg	2.34	if (rv[1] <= FTINY)
418			d = 1.0;
419			else
420			d = sqrt( np->alpha2 * -log(rv[1]) );
421			for (i = 0; i < 3; i++)
422			h[i] = np->pnorm[i] + d(cospu[i] + sinp*v[i]);
423	greg	2.63	d = -2.0 * DOT(h, np->rp->rdir) / (1.0 + d*d);
424			VSUM(sr.rdir, np->rp->rdir, h, d);
425	greg	2.58	/* sample rejection test */
426	greg	2.63	if ((d = DOT(sr.rdir, np->rp->ron)) <= FTINY)
427	greg	2.55	continue;
428			checknorm(sr.rdir);
429	greg	2.58	if (nstarget > 1) { /* W-G-M-D adjustment */
430	greg	2.59	if (nstaken) rayclear(&sr);
431	greg	2.58	rayvalue(&sr);
432	greg	2.63	d = 2./(1. + np->rp->rod/d);
433	greg	2.58	scalecolor(sr.rcol, d);
434			addcolor(scol, sr.rcol);
435			} else {
436			rayvalue(&sr);
437			multcolor(sr.rcol, sr.rcoef);
438	greg	2.63	addcolor(np->rp->rcol, sr.rcol);
439	greg	2.34	}
440	greg	2.58	++nstaken;
441			}
442			if (nstarget > 1) { /* final W-G-M-D weighting */
443			multcolor(scol, sr.rcoef);
444			d = (double)nstarget/ntrials;
445			scalecolor(scol, d);
446	greg	2.63	addcolor(np->rp->rcol, scol);
447	greg	2.34	}
448	greg	2.2	ndims--;
449			}
450			/* compute transmission */
451	greg	2.50	copycolor(sr.rcoef, np->mcolor); /* modified by color */
452			scalecolor(sr.rcoef, np->tspec);
453	greg	2.8	if ((np->specfl & (SP_TRAN\|SP_TBLT)) == SP_TRAN &&
454	greg	2.63	rayorigin(&sr, SPECULAR, np->rp, sr.rcoef) == 0) {
455	greg	2.58	nstarget = 1;
456	greg	2.55	if (specjitter > 1.5) { /* multiple samples? */
457	greg	2.63	nstarget = specjitter*np->rp->rweight + .5;
458	greg	2.58	if (sr.rweight <= minweight*nstarget)
459			nstarget = sr.rweight/minweight;
460			if (nstarget > 1) {
461			d = 1./nstarget;
462	greg	2.56	scalecolor(sr.rcoef, d);
463			sr.rweight *= d;
464	greg	2.55	} else
465	greg	2.58	nstarget = 1;
466	greg	2.55	}
467	greg	2.60	dimlist[ndims++] = (int)(size_t)np->mp;
468	greg	2.58	maxiter = MAXITER*nstarget;
469			for (nstaken = ntrials = 0; nstaken < nstarget &&
470			ntrials < maxiter; ntrials++) {
471			if (ntrials)
472	greg	2.34	d = frandom();
473			else
474	greg	2.58	d = urand(ilhash(dimlist,ndims)+samplendx);
475	greg	2.34	multisamp(rv, 2, d);
476			d = 2.0PI rv[0];
477	gwlarson	2.37	cosp = tcos(d);
478			sinp = tsin(d);
479	greg	2.55	if ((0. <= specjitter) & (specjitter < 1.))
480			rv[1] = 1.0 - specjitter*rv[1];
481	greg	2.34	if (rv[1] <= FTINY)
482			d = 1.0;
483			else
484	gwlarson	2.37	d = sqrt( np->alpha2 * -log(rv[1]) );
485	greg	2.34	for (i = 0; i < 3; i++)
486			sr.rdir[i] = np->prdir[i] + d(cospu[i] + sinp*v[i]);
487	greg	2.58	/* sample rejection test */
488	greg	2.63	if (DOT(sr.rdir, np->rp->ron) >= -FTINY)
489	greg	2.55	continue;
490			normalize(sr.rdir); /* OK, normalize */
491	greg	2.59	if (nstaken) /* multi-sampling */
492	greg	2.55	rayclear(&sr);
493			rayvalue(&sr);
494			multcolor(sr.rcol, sr.rcoef);
495	greg	2.63	addcolor(np->rp->rcol, sr.rcol);
496	greg	2.58	++nstaken;
497	greg	2.34	}
498	greg	2.8	ndims--;
499			}
500	greg	1.1	}