src/rt/normal.c

#ifndef lint
static const char RCSid[] = "$Id: normal.c,v 2.75 2017/08/03 16:41:52 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 (ambRayInPmap(np->rp))
                return;         /* specular already in photon map */

        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 */
        /* No longer needed?
        if (shadowRayInPmap(r) || ambRayInPmap(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 {
                                                        /* perturb */
                                VSUB(nd.prdir, r->rdir, r->pert);
                                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) {
                if (mirtest > transtest+FTINY)
                        r->rt = mirdist;
                else if (transtest > FTINY)
                        r->rt = transdist;
                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.76
Committed:	Wed Jan 10 04:08:50 2018 UTC (6 years, 4 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R2
Changes since 2.75:	+3 -3 lines
Log Message:	Made behavior of BRTDfunc more consistent with other materials
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.76	static const char RCSid[] = "$Id: normal.c,v 2.75 2017/08/03 16:41:52 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.72
116	greg	2.69	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	greg	2.72
126			if (ambRayInPmap(np->rp))
127			return; /* specular already in photon map */
128
129	greg	2.2	if (ldot > FTINY && (np->specfl&(SP_REFL\|SP_PURE)) == SP_REFL) {
130	greg	1.3	/*
131			* Compute specular reflection coefficient using
132	greg	2.54	* Gaussian distribution model.
133	greg	1.3	*/
134	greg	2.3	/* roughness */
135	greg	2.16	dtmp = np->alpha2;
136	greg	2.3	/* + source if flat */
137			if (np->specfl & SP_FLAT)
138	greg	2.48	dtmp += omega * (0.25/PI);
139	greg	2.23	/* half vector */
140	greg	2.62	VSUB(vtmp, ldir, np->rp->rdir);
141	greg	2.16	d2 = DOT(vtmp, np->pnorm);
142	greg	2.23	d2 *= d2;
143	greg	2.54	d3 = DOT(vtmp,vtmp);
144			d4 = (d3 - d2) / d2;
145			/* new W-G-M-D model */
146			dtmp = exp(-d4/dtmp) * d3 / (PI * d2d2 dtmp);
147	greg	1.3	/* worth using? */
148			if (dtmp > FTINY) {
149			copycolor(ctmp, np->scolor);
150	greg	2.54	dtmp = ldot omega;
151	greg	1.3	scalecolor(ctmp, dtmp);
152			addcolor(cval, ctmp);
153			}
154			}
155	greg	2.69
156
157	greg	2.2	if (ldot < -FTINY && (np->specfl&(SP_TRAN\|SP_PURE)) == SP_TRAN) {
158	greg	1.3	/*
159	greg	1.4	* Compute specular transmission. Specular transmission
160	greg	1.13	* is always modified by material color.
161	greg	1.3	*/
162			/* roughness + source */
163	greg	2.48	dtmp = np->alpha2 + omega*(1.0/PI);
164	greg	2.54	/* Gaussian */
165	greg	2.53	dtmp = exp((2.DOT(np->prdir,ldir)-2.)/dtmp)/(PIdtmp);
166	greg	1.3	/* worth using? */
167			if (dtmp > FTINY) {
168	greg	1.13	copycolor(ctmp, np->mcolor);
169	greg	2.52	dtmp = np->tspec omega * sqrt(-ldot/np->pdot);
170	greg	1.13	scalecolor(ctmp, dtmp);
171	greg	1.3	addcolor(cval, ctmp);
172			}
173			}
174			}
175
176
177	greg	2.62	int
178	schorsch	2.47	m_normal( /* color a ray that hit something normal */
179	greg	2.62	OBJREC *m,
180			RAY *r
181	schorsch	2.47	)
182	greg	1.3	{
183			NORMDAT nd;
184	greg	2.38	double fest;
185	greg	1.9	double transtest, transdist;
186	greg	2.29	double mirtest, mirdist;
187			int hastexture;
188			double d;
189	greg	1.1	COLOR ctmp;
190	greg	2.62	int i;
191	greg	2.69
192			/* PMAP: skip transmitted shadow ray if accounted for in photon map */
193	rschregle	2.74	/* No longer needed?
194	greg	2.73	if (shadowRayInPmap(r) \|\| ambRayInPmap(r))
195	rschregle	2.74	return(1); */
196
197	greg	1.1	/* easy shadow test */
198			if (r->crtype & SHADOW && m->otype != MAT_TRANS)
199	greg	2.27	return(1);
200	greg	2.2
201			if (m->oargs.nfargs != (m->otype == MAT_TRANS ? 7 : 5))
202			objerror(m, USER, "bad number of arguments");
203	greg	2.29	/* check for back side */
204			if (r->rod < 0.0) {
205	greg	2.66	if (!backvis) {
206	greg	2.29	raytrans(r);
207			return(1);
208			}
209	greg	2.40	raytexture(r, m->omod);
210	greg	2.29	flipsurface(r); /* reorient if backvis */
211	greg	2.40	} else
212			raytexture(r, m->omod);
213	greg	1.3	nd.mp = m;
214	greg	2.16	nd.rp = r;
215	greg	1.1	/* get material color */
216	greg	1.3	setcolor(nd.mcolor, m->oargs.farg[0],
217	greg	1.1	m->oargs.farg[1],
218			m->oargs.farg[2]);
219			/* get roughness */
220	greg	2.2	nd.specfl = 0;
221	greg	1.3	nd.alpha2 = m->oargs.farg[4];
222	greg	2.2	if ((nd.alpha2 *= nd.alpha2) <= FTINY)
223			nd.specfl \|= SP_PURE;
224	greg	2.40
225	schorsch	2.45	if ( (hastexture = (DOT(r->pert,r->pert) > FTINY*FTINY)) ) {
226	greg	2.29	nd.pdot = raynormal(nd.pnorm, r); /* perturb normal */
227	greg	2.41	} else {
228	greg	2.29	VCOPY(nd.pnorm, r->ron);
229			nd.pdot = r->rod;
230			}
231	greg	2.42	if (r->ro != NULL && isflat(r->ro->otype))
232			nd.specfl \|= SP_FLAT;
233	greg	1.13	if (nd.pdot < .001)
234			nd.pdot = .001; /* non-zero for dirnorm() */
235	greg	1.3	multcolor(nd.mcolor, r->pcol); /* modify material color */
236	greg	2.29	mirtest = transtest = 0;
237			mirdist = transdist = r->rot;
238	greg	2.30	nd.rspec = m->oargs.farg[3];
239	greg	2.38	/* compute Fresnel approx. */
240	greg	2.51	if (nd.specfl & SP_PURE && nd.rspec >= FRESTHRESH) {
241	greg	2.62	fest = FRESNE(nd.pdot);
242	greg	2.38	nd.rspec += fest*(1. - nd.rspec);
243			} else
244			fest = 0.;
245	greg	1.3	/* compute transmission */
246	greg	1.1	if (m->otype == MAT_TRANS) {
247	greg	1.3	nd.trans = m->oargs.farg[5]*(1.0 - nd.rspec);
248			nd.tspec = nd.trans * m->oargs.farg[6];
249			nd.tdiff = nd.trans - nd.tspec;
250	greg	2.2	if (nd.tspec > FTINY) {
251			nd.specfl \|= SP_TRAN;
252	greg	2.5	/* check threshold */
253	greg	2.25	if (!(nd.specfl & SP_PURE) &&
254			specthresh >= nd.tspec-FTINY)
255	greg	2.5	nd.specfl \|= SP_TBLT;
256	greg	2.67	if (!hastexture \|\| r->crtype & (SHADOW\|AMBIENT)) {
257	greg	2.2	VCOPY(nd.prdir, r->rdir);
258			transtest = 2;
259			} else {
260	greg	2.76	/* perturb */
261			VSUB(nd.prdir, r->rdir, r->pert);
262	greg	2.7	if (DOT(nd.prdir, r->ron) < -FTINY)
263			normalize(nd.prdir); /* OK */
264			else
265			VCOPY(nd.prdir, r->rdir);
266	greg	2.2	}
267	greg	1.14	}
268	greg	1.1	} else
269	greg	1.3	nd.tdiff = nd.tspec = nd.trans = 0.0;
270	greg	1.1	/* transmitted ray */
271	greg	2.69
272	greg	2.71	if ((nd.specfl&(SP_TRAN\|SP_PURE\|SP_TBLT)) == (SP_TRAN\|SP_PURE)) {
273	greg	1.3	RAY lr;
274	greg	2.50	copycolor(lr.rcoef, nd.mcolor); /* modified by color */
275			scalecolor(lr.rcoef, nd.tspec);
276			if (rayorigin(&lr, TRANS, r, lr.rcoef) == 0) {
277	greg	1.14	VCOPY(lr.rdir, nd.prdir);
278	greg	1.1	rayvalue(&lr);
279	greg	2.50	multcolor(lr.rcol, lr.rcoef);
280	greg	1.1	addcolor(r->rcol, lr.rcol);
281	greg	1.9	transtest *= bright(lr.rcol);
282			transdist = r->rot + lr.rt;
283	greg	1.1	}
284	greg	2.11	} else
285			transtest = 0;
286	greg	2.2
287	greg	2.29	if (r->crtype & SHADOW) { /* the rest is shadow */
288			r->rt = transdist;
289	greg	2.27	return(1);
290	greg	2.30	}
291			/* get specular reflection */
292			if (nd.rspec > FTINY) {
293			nd.specfl \|= SP_REFL;
294			/* compute specular color */
295	greg	2.38	if (m->otype != MAT_METAL) {
296			setcolor(nd.scolor, nd.rspec, nd.rspec, nd.rspec);
297			} else if (fest > FTINY) {
298	greg	2.64	d = m->oargs.farg[3]*(1. - fest);
299	greg	2.38	for (i = 0; i < 3; i++)
300	greg	2.65	colval(nd.scolor,i) = fest +
301			colval(nd.mcolor,i)*d;
302	greg	2.38	} else {
303	greg	2.30	copycolor(nd.scolor, nd.mcolor);
304	greg	2.38	scalecolor(nd.scolor, nd.rspec);
305			}
306	greg	2.30	/* check threshold */
307			if (!(nd.specfl & SP_PURE) && specthresh >= nd.rspec-FTINY)
308			nd.specfl \|= SP_RBLT;
309			/* compute reflected ray */
310	greg	2.55	VSUM(nd.vrefl, r->rdir, nd.pnorm, 2.*nd.pdot);
311	greg	2.30	/* penetration? */
312			if (hastexture && DOT(nd.vrefl, r->ron) <= FTINY)
313	greg	2.55	VSUM(nd.vrefl, r->rdir, r->ron, 2.*r->rod);
314	greg	2.53	checknorm(nd.vrefl);
315	gregl	2.36	}
316			/* reflected ray */
317	greg	2.71	if ((nd.specfl&(SP_REFL\|SP_PURE\|SP_RBLT)) == (SP_REFL\|SP_PURE)) {
318	gregl	2.36	RAY lr;
319	greg	2.50	if (rayorigin(&lr, REFLECTED, r, nd.scolor) == 0) {
320	gregl	2.36	VCOPY(lr.rdir, nd.vrefl);
321			rayvalue(&lr);
322	greg	2.50	multcolor(lr.rcol, lr.rcoef);
323	gregl	2.36	addcolor(r->rcol, lr.rcol);
324	greg	2.67	if (nd.specfl & SP_FLAT &&
325			!hastexture \| (r->crtype & AMBIENT)) {
326	gregl	2.36	mirtest = 2.*bright(lr.rcol);
327			mirdist = r->rot + lr.rt;
328	greg	2.30	}
329			}
330	greg	2.29	}
331	greg	1.1	/* diffuse reflection */
332	greg	1.3	nd.rdiff = 1.0 - nd.trans - nd.rspec;
333	greg	1.1
334	greg	2.75	if (nd.specfl & SP_PURE && nd.rdiff <= FTINY && nd.tdiff <= FTINY) {
335			if (mirtest > transtest+FTINY)
336			r->rt = mirdist;
337			else if (transtest > FTINY)
338			r->rt = transdist;
339	greg	2.27	return(1); /* 100% pure specular */
340	greg	2.75	}
341	greg	2.71	if (!(nd.specfl & SP_PURE))
342			gaussamp(&nd); /* checks BLT flags /
343	greg	2.2
344	greg	1.3	if (nd.rdiff > FTINY) { /* ambient from this side */
345	greg	2.50	copycolor(ctmp, nd.mcolor); /* modified by material color */
346	greg	2.61	scalecolor(ctmp, nd.rdiff);
347			if (nd.specfl & SP_RBLT) /* add in specular as well? */
348			addcolor(ctmp, nd.scolor);
349	greg	2.50	multambient(ctmp, r, hastexture ? nd.pnorm : r->ron);
350	greg	1.2	addcolor(r->rcol, ctmp); /* add to returned color */
351			}
352	greg	1.3	if (nd.tdiff > FTINY) { /* ambient from other side */
353	greg	2.50	copycolor(ctmp, nd.mcolor); /* modified by color */
354			if (nd.specfl & SP_TBLT)
355			scalecolor(ctmp, nd.trans);
356			else
357			scalecolor(ctmp, nd.tdiff);
358	greg	1.1	flipsurface(r);
359	greg	2.32	if (hastexture) {
360			FVECT bnorm;
361			bnorm[0] = -nd.pnorm[0];
362			bnorm[1] = -nd.pnorm[1];
363			bnorm[2] = -nd.pnorm[2];
364	greg	2.50	multambient(ctmp, r, bnorm);
365	greg	2.32	} else
366	greg	2.50	multambient(ctmp, r, r->ron);
367	greg	1.1	addcolor(r->rcol, ctmp);
368			flipsurface(r);
369			}
370	greg	1.3	/* add direct component */
371			direct(r, dirnorm, &nd);
372	greg	1.9	/* check distance */
373	greg	2.29	d = bright(r->rcol);
374			if (transtest > d)
375	greg	1.9	r->rt = transdist;
376	greg	2.29	else if (mirtest > d)
377			r->rt = mirdist;
378	greg	2.27
379			return(1);
380	greg	2.2	}
381
382
383	greg	2.38	static void
384	greg	2.54	gaussamp( /* sample Gaussian specular */
385	greg	2.62	NORMDAT *np
386	schorsch	2.47	)
387	greg	2.2	{
388			RAY sr;
389			FVECT u, v, h;
390			double rv[2];
391			double d, sinp, cosp;
392	greg	2.62	COLOR scol;
393	greg	2.58	int maxiter, ntrials, nstarget, nstaken;
394	greg	2.62	int i;
395	greg	2.13	/* quick test */
396			if ((np->specfl & (SP_REFL\|SP_RBLT)) != SP_REFL &&
397			(np->specfl & (SP_TRAN\|SP_TBLT)) != SP_TRAN)
398			return;
399	greg	2.2	/* set up sample coordinates */
400	greg	2.70	getperpendicular(u, np->pnorm, rand_samp);
401	greg	2.2	fcross(v, np->pnorm, u);
402			/* compute reflection */
403	greg	2.5	if ((np->specfl & (SP_REFL\|SP_RBLT)) == SP_REFL &&
404	greg	2.63	rayorigin(&sr, SPECULAR, np->rp, np->scolor) == 0) {
405	greg	2.58	nstarget = 1;
406	greg	2.55	if (specjitter > 1.5) { /* multiple samples? */
407	greg	2.63	nstarget = specjitter*np->rp->rweight + .5;
408	greg	2.58	if (sr.rweight <= minweight*nstarget)
409			nstarget = sr.rweight/minweight;
410			if (nstarget > 1) {
411			d = 1./nstarget;
412			scalecolor(sr.rcoef, d);
413	greg	2.56	sr.rweight *= d;
414	greg	2.55	} else
415	greg	2.58	nstarget = 1;
416	greg	2.55	}
417	greg	2.58	setcolor(scol, 0., 0., 0.);
418	greg	2.60	dimlist[ndims++] = (int)(size_t)np->mp;
419	greg	2.58	maxiter = MAXITER*nstarget;
420			for (nstaken = ntrials = 0; nstaken < nstarget &&
421			ntrials < maxiter; ntrials++) {
422			if (ntrials)
423	greg	2.34	d = frandom();
424			else
425			d = urand(ilhash(dimlist,ndims)+samplendx);
426			multisamp(rv, 2, d);
427			d = 2.0PI rv[0];
428	gwlarson	2.37	cosp = tcos(d);
429			sinp = tsin(d);
430	greg	2.55	if ((0. <= specjitter) & (specjitter < 1.))
431			rv[1] = 1.0 - specjitter*rv[1];
432	greg	2.34	if (rv[1] <= FTINY)
433			d = 1.0;
434			else
435			d = sqrt( np->alpha2 * -log(rv[1]) );
436			for (i = 0; i < 3; i++)
437			h[i] = np->pnorm[i] + d(cospu[i] + sinp*v[i]);
438	greg	2.63	d = -2.0 * DOT(h, np->rp->rdir) / (1.0 + d*d);
439			VSUM(sr.rdir, np->rp->rdir, h, d);
440	greg	2.58	/* sample rejection test */
441	greg	2.63	if ((d = DOT(sr.rdir, np->rp->ron)) <= FTINY)
442	greg	2.55	continue;
443			checknorm(sr.rdir);
444	greg	2.58	if (nstarget > 1) { /* W-G-M-D adjustment */
445	greg	2.59	if (nstaken) rayclear(&sr);
446	greg	2.58	rayvalue(&sr);
447	greg	2.63	d = 2./(1. + np->rp->rod/d);
448	greg	2.58	scalecolor(sr.rcol, d);
449			addcolor(scol, sr.rcol);
450			} else {
451			rayvalue(&sr);
452			multcolor(sr.rcol, sr.rcoef);
453	greg	2.63	addcolor(np->rp->rcol, sr.rcol);
454	greg	2.34	}
455	greg	2.58	++nstaken;
456			}
457			if (nstarget > 1) { /* final W-G-M-D weighting */
458			multcolor(scol, sr.rcoef);
459			d = (double)nstarget/ntrials;
460			scalecolor(scol, d);
461	greg	2.63	addcolor(np->rp->rcol, scol);
462	greg	2.34	}
463	greg	2.2	ndims--;
464			}
465			/* compute transmission */
466	greg	2.50	copycolor(sr.rcoef, np->mcolor); /* modified by color */
467			scalecolor(sr.rcoef, np->tspec);
468	greg	2.8	if ((np->specfl & (SP_TRAN\|SP_TBLT)) == SP_TRAN &&
469	greg	2.63	rayorigin(&sr, SPECULAR, np->rp, sr.rcoef) == 0) {
470	greg	2.58	nstarget = 1;
471	greg	2.55	if (specjitter > 1.5) { /* multiple samples? */
472	greg	2.63	nstarget = specjitter*np->rp->rweight + .5;
473	greg	2.58	if (sr.rweight <= minweight*nstarget)
474			nstarget = sr.rweight/minweight;
475			if (nstarget > 1) {
476			d = 1./nstarget;
477	greg	2.56	scalecolor(sr.rcoef, d);
478			sr.rweight *= d;
479	greg	2.55	} else
480	greg	2.58	nstarget = 1;
481	greg	2.55	}
482	greg	2.60	dimlist[ndims++] = (int)(size_t)np->mp;
483	greg	2.58	maxiter = MAXITER*nstarget;
484			for (nstaken = ntrials = 0; nstaken < nstarget &&
485			ntrials < maxiter; ntrials++) {
486			if (ntrials)
487	greg	2.34	d = frandom();
488			else
489	greg	2.58	d = urand(ilhash(dimlist,ndims)+samplendx);
490	greg	2.34	multisamp(rv, 2, d);
491			d = 2.0PI rv[0];
492	gwlarson	2.37	cosp = tcos(d);
493			sinp = tsin(d);
494	greg	2.55	if ((0. <= specjitter) & (specjitter < 1.))
495			rv[1] = 1.0 - specjitter*rv[1];
496	greg	2.34	if (rv[1] <= FTINY)
497			d = 1.0;
498			else
499	gwlarson	2.37	d = sqrt( np->alpha2 * -log(rv[1]) );
500	greg	2.34	for (i = 0; i < 3; i++)
501			sr.rdir[i] = np->prdir[i] + d(cospu[i] + sinp*v[i]);
502	greg	2.58	/* sample rejection test */
503	greg	2.63	if (DOT(sr.rdir, np->rp->ron) >= -FTINY)
504	greg	2.55	continue;
505			normalize(sr.rdir); /* OK, normalize */
506	greg	2.59	if (nstaken) /* multi-sampling */
507	greg	2.55	rayclear(&sr);
508			rayvalue(&sr);
509			multcolor(sr.rcol, sr.rcoef);
510	greg	2.63	addcolor(np->rp->rcol, sr.rcol);
511	greg	2.58	++nstaken;
512	greg	2.34	}
513	greg	2.8	ndims--;
514			}
515	greg	1.1	}