src/rt/normal.c

#ifndef lint
static const char RCSid[] = "$Id: normal.c,v 2.86 2024/12/19 23:25:28 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.00202943064)
#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 */
        SCOLOR  mcolor;         /* color of this material */
        SCOLOR  scolor;         /* color of specular component */
        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 */
        SCOLOR  scval,                  /* 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;
        SCOLOR  sctmp;

        scolorblack(scval);

        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.
                 */
                copyscolor(sctmp, np->mcolor);
                dtmp = ldot * omega * lrdiff * (1.0/PI);
                scalescolor(sctmp, dtmp);
                saddscolor(scval, sctmp);
        }

        if ((ldot < -FTINY) & (ltdiff > FTINY)) {
                /*
                 *  Compute diffuse transmission.
                 */
                copyscolor(sctmp, np->mcolor);
                dtmp = -ldot * omega * ltdiff * (1.0/PI);
                scalescolor(sctmp, dtmp);
                saddscolor(scval, sctmp);
        }

        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 += (1. - dstrsrc) * 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) {
                        copyscolor(sctmp, np->scolor);
                        dtmp *= ldot * omega;
                        scalescolor(sctmp, dtmp);
                        saddscolor(scval, sctmp);
                }
        }
        

        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) {
                        copyscolor(sctmp, np->mcolor);
                        dtmp *= np->tspec * omega * sqrt(-ldot/np->pdot);
                        scalescolor(sctmp, dtmp);
                        saddscolor(scval, sctmp);
                }
        }
}


int
m_normal(                       /* color a ray that hit something normal */
        OBJREC  *m,
        RAY  *r
)
{
        NORMDAT  nd;
        double  fest;
        int     hastexture;
        double  d;
        SCOLOR  sctmp;
        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 */
        setscolor(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 (!hastexture && r->ro != NULL && isflat(r->ro->otype))
                nd.specfl |= SP_FLAT;
        if (nd.pdot < .001)
                nd.pdot = .001;                 /* non-zero for dirnorm() */
        smultscolor(nd.mcolor, r->pcol);        /* modify material color */
        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);
                        } 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;
                                                /* diffuse reflection */
        nd.rdiff = 1.0 - nd.trans - nd.rspec;
                                                /* transmitted ray */
        if ((nd.specfl&(SP_TRAN|SP_PURE|SP_TBLT)) == (SP_TRAN|SP_PURE)) {
                RAY  lr;
                copyscolor(lr.rcoef, nd.mcolor);        /* modified by color */
                scalescolor(lr.rcoef, nd.tspec);
                if (rayorigin(&lr, TRANS, r, lr.rcoef) == 0) {
                        VCOPY(lr.rdir, nd.prdir);
                        rayvalue(&lr);
                        smultscolor(lr.rcol, lr.rcoef);
                        saddscolor(r->rcol, lr.rcol);
                        if (nd.tspec >= 1.0-FTINY) {
                                                /* completely transparent */
                                smultscolor(lr.mcol, lr.rcoef);
                                copyscolor(r->mcol, lr.mcol);
                                r->rmt = r->rot + lr.rmt;
                                r->rxt = r->rot + lr.rxt;
                        } else if (nd.tspec > nd.tdiff + nd.rdiff)
                                r->rxt = r->rot + raydistance(&lr);
                }
        }

        if (r->crtype & SHADOW)                 /* the rest is shadow */
                return(1);
                                                /* get specular reflection */
        if (nd.rspec > FTINY) {
                nd.specfl |= SP_REFL;
                                                /* compute specular color */
                if (m->otype != MAT_METAL) {
                        setscolor(nd.scolor, nd.rspec, nd.rspec, nd.rspec);
                } else if (fest > FTINY) {
                        d = m->oargs.farg[3]*(1. - fest);
                        for (i = NCSAMP; i--; )
                                nd.scolor[i] = fest + nd.mcolor[i]*d;
                } else {
                        copyscolor(nd.scolor, nd.mcolor);
                        scalescolor(nd.scolor, nd.rspec);
                }
                                                /* check threshold */
                if (!(nd.specfl & SP_PURE) && specthresh >= nd.rspec-FTINY)
                        nd.specfl |= SP_RBLT;
        }
                                                /* 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) {
                                                /* compute reflected ray */
                        VSUM(lr.rdir, r->rdir, nd.pnorm, 2.*nd.pdot);
                                                /* penetration? */
                        if (hastexture && DOT(lr.rdir, r->ron) <= FTINY)
                                VSUM(lr.rdir, r->rdir, r->ron, 2.*r->rod);
                        checknorm(lr.rdir);
                        rayvalue(&lr);
                        smultscolor(lr.rcol, lr.rcoef);
                        copyscolor(r->mcol, lr.rcol);
                        saddscolor(r->rcol, lr.rcol);
                        r->rmt = r->rot;
                        if (nd.specfl & SP_FLAT && r->crtype & AMBIENT)
                                r->rmt += raydistance(&lr);
                }
        }

        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 */
                copyscolor(sctmp, nd.mcolor);   /* modified by material color */
                scalescolor(sctmp, nd.rdiff);
                if (nd.specfl & SP_RBLT)        /* add in specular as well? */
                        saddscolor(sctmp, nd.scolor);
                multambient(sctmp, r, nd.pnorm);
                saddscolor(r->rcol, sctmp);     /* add to returned color */
        }
        if (nd.tdiff > FTINY) {         /* ambient from other side */
                FVECT  bnorm;
                copyscolor(sctmp, nd.mcolor);   /* modified by color */
                if (nd.specfl & SP_TBLT) {
                        scalescolor(sctmp, nd.trans);
                } else {
                        scalescolor(sctmp, nd.tdiff);
                }
                bnorm[0] = -nd.pnorm[0];
                bnorm[1] = -nd.pnorm[1];
                bnorm[2] = -nd.pnorm[2];
                multambient(sctmp, r, bnorm);
                saddscolor(r->rcol, sctmp);
        }
                                        /* add direct component */
        direct(r, dirnorm, &nd);

        return(1);
}


static void
gaussamp(                       /* sample Gaussian specular */
        NORMDAT  *np
)
{
        RAY  sr;
        FVECT  u, v, h;
        double  rv[2];
        double  d, sinp, cosp;
        SCOLOR  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, RSPECULAR, 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;
                                scalescolor(sr.rcoef, d);
                                sr.rweight *= d;
                        } else
                                nstarget = 1;
                }
                scolorblack(scol);
                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);
                                scalescolor(sr.rcol, d);
                                saddscolor(scol, sr.rcol);
                        } else {
                                rayvalue(&sr);
                                smultscolor(sr.rcol, sr.rcoef);
                                saddscolor(np->rp->rcol, sr.rcol);
                        }
                        ++nstaken;
                }
                if (nstarget > 1) {             /* final W-G-M-D weighting */
                        smultscolor(scol, sr.rcoef);
                        d = (double)nstarget/ntrials;
                        scalescolor(scol, d);
                        saddscolor(np->rp->rcol, scol);
                }
                ndims--;
        }
                                        /* compute transmission */
        copyscolor(sr.rcoef, np->mcolor);       /* modified by color */
        scalescolor(sr.rcoef, np->tspec);
        if ((np->specfl & (SP_TRAN|SP_TBLT)) == SP_TRAN &&
                        rayorigin(&sr, TSPECULAR, 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;
                                scalescolor(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);
                        smultscolor(sr.rcol, sr.rcoef);
                        saddscolor(np->rp->rcol, sr.rcol);
                        ++nstaken;
                }
                ndims--;
        }
}
Revision:	2.87
Committed:	Fri Dec 20 16:29:50 2024 UTC (4 months, 1 week ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 2.86:	+2 -2 lines
Log Message:	perf: Adjustment to source spread in lobe speculars that accounts for -dj
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.87	static const char RCSid[] = "$Id: normal.c,v 2.86 2024/12/19 23:25:28 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.82	#define FRESNE(ci) (exp(-5.85*(ci)) - 0.00202943064)
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	2.83	SCOLOR mcolor; /* color of this material */
58			SCOLOR scolor; /* color of specular component */
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	greg	2.83	SCOLOR scval, /* 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	2.83	SCOLOR sctmp;
85	greg	1.3
86	greg	2.83	scolorblack(scval);
87	greg	1.3
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.83	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	greg	2.83	copyscolor(sctmp, np->mcolor);
110	greg	2.49	dtmp = ldot * omega * lrdiff * (1.0/PI);
111	greg	2.83	scalescolor(sctmp, dtmp);
112			saddscolor(scval, sctmp);
113	greg	1.3	}
114	greg	2.72
115	greg	2.83	if ((ldot < -FTINY) & (ltdiff > FTINY)) {
116	greg	2.69	/*
117			* Compute diffuse transmission.
118			*/
119	greg	2.83	copyscolor(sctmp, np->mcolor);
120	greg	2.69	dtmp = -ldot * omega * ltdiff * (1.0/PI);
121	greg	2.83	scalescolor(sctmp, dtmp);
122			saddscolor(scval, sctmp);
123	greg	2.69	}
124	greg	2.72
125			if (ambRayInPmap(np->rp))
126			return; /* specular already in photon map */
127
128	greg	2.83	if ((ldot > FTINY) & ((np->specfl&(SP_REFL\|SP_PURE)) == SP_REFL)) {
129	greg	1.3	/*
130			* Compute specular reflection coefficient using
131	greg	2.54	* Gaussian distribution model.
132	greg	1.3	*/
133	greg	2.3	/* roughness */
134	greg	2.16	dtmp = np->alpha2;
135	greg	2.3	/* + source if flat */
136			if (np->specfl & SP_FLAT)
137	greg	2.87	dtmp += (1. - dstrsrc) * omega * (0.25/PI);
138	greg	2.23	/* half vector */
139	greg	2.62	VSUB(vtmp, ldir, np->rp->rdir);
140	greg	2.16	d2 = DOT(vtmp, np->pnorm);
141	greg	2.23	d2 *= d2;
142	greg	2.54	d3 = DOT(vtmp,vtmp);
143			d4 = (d3 - d2) / d2;
144			/* new W-G-M-D model */
145			dtmp = exp(-d4/dtmp) * d3 / (PI * d2d2 dtmp);
146	greg	1.3	/* worth using? */
147			if (dtmp > FTINY) {
148	greg	2.83	copyscolor(sctmp, np->scolor);
149	greg	2.54	dtmp = ldot omega;
150	greg	2.83	scalescolor(sctmp, dtmp);
151			saddscolor(scval, sctmp);
152	greg	1.3	}
153			}
154	greg	2.69
155
156	greg	2.83	if ((ldot < -FTINY) & ((np->specfl&(SP_TRAN\|SP_PURE)) == SP_TRAN)) {
157	greg	1.3	/*
158	greg	1.4	* Compute specular transmission. Specular transmission
159	greg	1.13	* is always modified by material color.
160	greg	1.3	*/
161			/* roughness + source */
162	greg	2.48	dtmp = np->alpha2 + omega*(1.0/PI);
163	greg	2.54	/* Gaussian */
164	greg	2.53	dtmp = exp((2.DOT(np->prdir,ldir)-2.)/dtmp)/(PIdtmp);
165	greg	1.3	/* worth using? */
166			if (dtmp > FTINY) {
167	greg	2.83	copyscolor(sctmp, np->mcolor);
168	greg	2.52	dtmp = np->tspec omega * sqrt(-ldot/np->pdot);
169	greg	2.83	scalescolor(sctmp, dtmp);
170			saddscolor(scval, sctmp);
171	greg	1.3	}
172			}
173			}
174
175
176	greg	2.62	int
177	schorsch	2.47	m_normal( /* color a ray that hit something normal */
178	greg	2.62	OBJREC *m,
179			RAY *r
180	schorsch	2.47	)
181	greg	1.3	{
182			NORMDAT nd;
183	greg	2.38	double fest;
184	greg	2.29	int hastexture;
185			double d;
186	greg	2.83	SCOLOR sctmp;
187	greg	2.62	int i;
188	greg	2.69
189			/* PMAP: skip transmitted shadow ray if accounted for in photon map */
190	rschregle	2.74	/* No longer needed?
191	greg	2.73	if (shadowRayInPmap(r) \|\| ambRayInPmap(r))
192	rschregle	2.74	return(1); */
193
194	greg	1.1	/* easy shadow test */
195			if (r->crtype & SHADOW && m->otype != MAT_TRANS)
196	greg	2.27	return(1);
197	greg	2.2
198			if (m->oargs.nfargs != (m->otype == MAT_TRANS ? 7 : 5))
199			objerror(m, USER, "bad number of arguments");
200	greg	2.29	/* check for back side */
201			if (r->rod < 0.0) {
202	greg	2.66	if (!backvis) {
203	greg	2.29	raytrans(r);
204			return(1);
205			}
206	greg	2.40	raytexture(r, m->omod);
207	greg	2.29	flipsurface(r); /* reorient if backvis */
208	greg	2.40	} else
209			raytexture(r, m->omod);
210	greg	1.3	nd.mp = m;
211	greg	2.16	nd.rp = r;
212	greg	1.1	/* get material color */
213	greg	2.83	setscolor(nd.mcolor, m->oargs.farg[0],
214	greg	1.1	m->oargs.farg[1],
215			m->oargs.farg[2]);
216			/* get roughness */
217	greg	2.2	nd.specfl = 0;
218	greg	1.3	nd.alpha2 = m->oargs.farg[4];
219	greg	2.2	if ((nd.alpha2 *= nd.alpha2) <= FTINY)
220			nd.specfl \|= SP_PURE;
221	greg	2.40
222	schorsch	2.45	if ( (hastexture = (DOT(r->pert,r->pert) > FTINY*FTINY)) ) {
223	greg	2.29	nd.pdot = raynormal(nd.pnorm, r); /* perturb normal */
224	greg	2.41	} else {
225	greg	2.29	VCOPY(nd.pnorm, r->ron);
226			nd.pdot = r->rod;
227			}
228	greg	2.86	if (!hastexture && r->ro != NULL && isflat(r->ro->otype))
229	greg	2.42	nd.specfl \|= SP_FLAT;
230	greg	1.13	if (nd.pdot < .001)
231			nd.pdot = .001; /* non-zero for dirnorm() */
232	greg	2.83	smultscolor(nd.mcolor, r->pcol); /* modify material color */
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			} else {
254	greg	2.76	/* perturb */
255			VSUB(nd.prdir, r->rdir, r->pert);
256	greg	2.7	if (DOT(nd.prdir, r->ron) < -FTINY)
257			normalize(nd.prdir); /* OK */
258			else
259			VCOPY(nd.prdir, r->rdir);
260	greg	2.2	}
261	greg	1.14	}
262	greg	1.1	} else
263	greg	1.3	nd.tdiff = nd.tspec = nd.trans = 0.0;
264	greg	2.79	/* diffuse reflection */
265			nd.rdiff = 1.0 - nd.trans - nd.rspec;
266	greg	1.1	/* transmitted ray */
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.83	copyscolor(lr.rcoef, nd.mcolor); /* modified by color */
270			scalescolor(lr.rcoef, nd.tspec);
271	greg	2.50	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.83	smultscolor(lr.rcol, lr.rcoef);
275			saddscolor(r->rcol, lr.rcol);
276	greg	2.78	if (nd.tspec >= 1.0-FTINY) {
277			/* completely transparent */
278	greg	2.83	smultscolor(lr.mcol, lr.rcoef);
279			copyscolor(r->mcol, lr.mcol);
280	greg	2.78	r->rmt = r->rot + lr.rmt;
281			r->rxt = r->rot + lr.rxt;
282	greg	2.79	} else if (nd.tspec > nd.tdiff + nd.rdiff)
283	greg	2.78	r->rxt = r->rot + raydistance(&lr);
284	greg	1.1	}
285	greg	2.77	}
286	greg	2.2
287	greg	2.77	if (r->crtype & SHADOW) /* the rest is shadow */
288	greg	2.27	return(1);
289	greg	2.30	/* get specular reflection */
290			if (nd.rspec > FTINY) {
291			nd.specfl \|= SP_REFL;
292			/* compute specular color */
293	greg	2.38	if (m->otype != MAT_METAL) {
294	greg	2.83	setscolor(nd.scolor, nd.rspec, nd.rspec, nd.rspec);
295	greg	2.38	} else if (fest > FTINY) {
296	greg	2.64	d = m->oargs.farg[3]*(1. - fest);
297	greg	2.83	for (i = NCSAMP; i--; )
298			nd.scolor[i] = fest + nd.mcolor[i]*d;
299	greg	2.38	} else {
300	greg	2.83	copyscolor(nd.scolor, nd.mcolor);
301			scalescolor(nd.scolor, nd.rspec);
302	greg	2.38	}
303	greg	2.30	/* check threshold */
304			if (!(nd.specfl & SP_PURE) && specthresh >= nd.rspec-FTINY)
305			nd.specfl \|= SP_RBLT;
306	gregl	2.36	}
307			/* reflected ray */
308	greg	2.71	if ((nd.specfl&(SP_REFL\|SP_PURE\|SP_RBLT)) == (SP_REFL\|SP_PURE)) {
309	gregl	2.36	RAY lr;
310	greg	2.50	if (rayorigin(&lr, REFLECTED, r, nd.scolor) == 0) {
311	greg	2.85	/* compute reflected ray */
312			VSUM(lr.rdir, r->rdir, nd.pnorm, 2.*nd.pdot);
313			/* penetration? */
314			if (hastexture && DOT(lr.rdir, r->ron) <= FTINY)
315			VSUM(lr.rdir, r->rdir, r->ron, 2.*r->rod);
316			checknorm(lr.rdir);
317	gregl	2.36	rayvalue(&lr);
318	greg	2.83	smultscolor(lr.rcol, lr.rcoef);
319			copyscolor(r->mcol, lr.rcol);
320			saddscolor(r->rcol, lr.rcol);
321	greg	2.80	r->rmt = r->rot;
322	greg	2.86	if (nd.specfl & SP_FLAT && r->crtype & AMBIENT)
323	greg	2.80	r->rmt += raydistance(&lr);
324	greg	2.30	}
325	greg	2.29	}
326	greg	1.1
327	greg	2.77	if (nd.specfl & SP_PURE && nd.rdiff <= FTINY && nd.tdiff <= FTINY)
328	greg	2.27	return(1); /* 100% pure specular */
329	greg	2.77
330	greg	2.71	if (!(nd.specfl & SP_PURE))
331			gaussamp(&nd); /* checks BLT flags /
332	greg	2.2
333	greg	1.3	if (nd.rdiff > FTINY) { /* ambient from this side */
334	greg	2.83	copyscolor(sctmp, nd.mcolor); /* modified by material color */
335			scalescolor(sctmp, nd.rdiff);
336	greg	2.61	if (nd.specfl & SP_RBLT) /* add in specular as well? */
337	greg	2.83	saddscolor(sctmp, nd.scolor);
338	greg	2.84	multambient(sctmp, r, nd.pnorm);
339	greg	2.83	saddscolor(r->rcol, sctmp); /* add to returned color */
340	greg	1.2	}
341	greg	1.3	if (nd.tdiff > FTINY) { /* ambient from other side */
342	greg	2.84	FVECT bnorm;
343	greg	2.83	copyscolor(sctmp, nd.mcolor); /* modified by color */
344			if (nd.specfl & SP_TBLT) {
345			scalescolor(sctmp, nd.trans);
346			} else {
347			scalescolor(sctmp, nd.tdiff);
348			}
349	greg	2.84	bnorm[0] = -nd.pnorm[0];
350			bnorm[1] = -nd.pnorm[1];
351			bnorm[2] = -nd.pnorm[2];
352			multambient(sctmp, r, bnorm);
353	greg	2.83	saddscolor(r->rcol, sctmp);
354	greg	1.1	}
355	greg	1.3	/* add direct component */
356			direct(r, dirnorm, &nd);
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.83	SCOLOR 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	greg	2.70	getperpendicular(u, np->pnorm, rand_samp);
380	greg	2.2	fcross(v, np->pnorm, u);
381			/* compute reflection */
382	greg	2.5	if ((np->specfl & (SP_REFL\|SP_RBLT)) == SP_REFL &&
383	greg	2.84	rayorigin(&sr, RSPECULAR, np->rp, np->scolor) == 0) {
384	greg	2.58	nstarget = 1;
385	greg	2.55	if (specjitter > 1.5) { /* multiple samples? */
386	greg	2.63	nstarget = specjitter*np->rp->rweight + .5;
387	greg	2.58	if (sr.rweight <= minweight*nstarget)
388			nstarget = sr.rweight/minweight;
389			if (nstarget > 1) {
390			d = 1./nstarget;
391	greg	2.83	scalescolor(sr.rcoef, d);
392	greg	2.56	sr.rweight *= d;
393	greg	2.55	} else
394	greg	2.58	nstarget = 1;
395	greg	2.55	}
396	greg	2.83	scolorblack(scol);
397	greg	2.60	dimlist[ndims++] = (int)(size_t)np->mp;
398	greg	2.58	maxiter = MAXITER*nstarget;
399			for (nstaken = ntrials = 0; nstaken < nstarget &&
400			ntrials < maxiter; ntrials++) {
401			if (ntrials)
402	greg	2.34	d = frandom();
403			else
404			d = urand(ilhash(dimlist,ndims)+samplendx);
405			multisamp(rv, 2, d);
406			d = 2.0PI rv[0];
407	gwlarson	2.37	cosp = tcos(d);
408			sinp = tsin(d);
409	greg	2.55	if ((0. <= specjitter) & (specjitter < 1.))
410			rv[1] = 1.0 - specjitter*rv[1];
411	greg	2.34	if (rv[1] <= FTINY)
412			d = 1.0;
413			else
414			d = sqrt( np->alpha2 * -log(rv[1]) );
415			for (i = 0; i < 3; i++)
416			h[i] = np->pnorm[i] + d(cospu[i] + sinp*v[i]);
417	greg	2.63	d = -2.0 * DOT(h, np->rp->rdir) / (1.0 + d*d);
418			VSUM(sr.rdir, np->rp->rdir, h, d);
419	greg	2.58	/* sample rejection test */
420	greg	2.63	if ((d = DOT(sr.rdir, np->rp->ron)) <= FTINY)
421	greg	2.55	continue;
422			checknorm(sr.rdir);
423	greg	2.58	if (nstarget > 1) { /* W-G-M-D adjustment */
424	greg	2.59	if (nstaken) rayclear(&sr);
425	greg	2.58	rayvalue(&sr);
426	greg	2.63	d = 2./(1. + np->rp->rod/d);
427	greg	2.83	scalescolor(sr.rcol, d);
428			saddscolor(scol, sr.rcol);
429	greg	2.58	} else {
430			rayvalue(&sr);
431	greg	2.83	smultscolor(sr.rcol, sr.rcoef);
432			saddscolor(np->rp->rcol, sr.rcol);
433	greg	2.34	}
434	greg	2.58	++nstaken;
435			}
436			if (nstarget > 1) { /* final W-G-M-D weighting */
437	greg	2.83	smultscolor(scol, sr.rcoef);
438	greg	2.58	d = (double)nstarget/ntrials;
439	greg	2.83	scalescolor(scol, d);
440			saddscolor(np->rp->rcol, scol);
441	greg	2.34	}
442	greg	2.2	ndims--;
443			}
444			/* compute transmission */
445	greg	2.83	copyscolor(sr.rcoef, np->mcolor); /* modified by color */
446			scalescolor(sr.rcoef, np->tspec);
447	greg	2.8	if ((np->specfl & (SP_TRAN\|SP_TBLT)) == SP_TRAN &&
448	greg	2.84	rayorigin(&sr, TSPECULAR, np->rp, sr.rcoef) == 0) {
449	greg	2.58	nstarget = 1;
450	greg	2.55	if (specjitter > 1.5) { /* multiple samples? */
451	greg	2.63	nstarget = specjitter*np->rp->rweight + .5;
452	greg	2.58	if (sr.rweight <= minweight*nstarget)
453			nstarget = sr.rweight/minweight;
454			if (nstarget > 1) {
455			d = 1./nstarget;
456	greg	2.83	scalescolor(sr.rcoef, d);
457	greg	2.56	sr.rweight *= d;
458	greg	2.55	} else
459	greg	2.58	nstarget = 1;
460	greg	2.55	}
461	greg	2.60	dimlist[ndims++] = (int)(size_t)np->mp;
462	greg	2.58	maxiter = MAXITER*nstarget;
463			for (nstaken = ntrials = 0; nstaken < nstarget &&
464			ntrials < maxiter; ntrials++) {
465			if (ntrials)
466	greg	2.34	d = frandom();
467			else
468	greg	2.58	d = urand(ilhash(dimlist,ndims)+samplendx);
469	greg	2.34	multisamp(rv, 2, d);
470			d = 2.0PI rv[0];
471	gwlarson	2.37	cosp = tcos(d);
472			sinp = tsin(d);
473	greg	2.55	if ((0. <= specjitter) & (specjitter < 1.))
474			rv[1] = 1.0 - specjitter*rv[1];
475	greg	2.34	if (rv[1] <= FTINY)
476			d = 1.0;
477			else
478	gwlarson	2.37	d = sqrt( np->alpha2 * -log(rv[1]) );
479	greg	2.34	for (i = 0; i < 3; i++)
480			sr.rdir[i] = np->prdir[i] + d(cospu[i] + sinp*v[i]);
481	greg	2.58	/* sample rejection test */
482	greg	2.63	if (DOT(sr.rdir, np->rp->ron) >= -FTINY)
483	greg	2.55	continue;
484			normalize(sr.rdir); /* OK, normalize */
485	greg	2.59	if (nstaken) /* multi-sampling */
486	greg	2.55	rayclear(&sr);
487			rayvalue(&sr);
488	greg	2.83	smultscolor(sr.rcol, sr.rcoef);
489			saddscolor(np->rp->rcol, sr.rcol);
490	greg	2.58	++nstaken;
491	greg	2.34	}
492	greg	2.8	ndims--;
493			}
494	greg	1.1	}