src/rt/normal.c

#ifndef lint
static const char RCSid[] = "$Id: normal.c,v 2.70 2015/05/21 05:54:54 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"
#include  "pmapmat.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 && 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_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 && (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;

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