src/rt/normal.c

#ifndef lint
static const char RCSid[] = "$Id: normal.c,v 2.82 2023/03/07 02:24:44 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  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 */
        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 += 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 (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;
                                                /* 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);
                        smultscolor(lr.rcol, lr.rcoef);
                        copyscolor(r->mcol, lr.rcol);
                        saddscolor(r->rcol, lr.rcol);
                        r->rmt = r->rot;
                        if (nd.specfl & SP_FLAT &&
                                        !hastexture | (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, hastexture ? nd.pnorm : r->ron);
                saddscolor(r->rcol, sctmp);     /* add to returned color */
        }
        if (nd.tdiff > FTINY) {         /* ambient from other side */
                copyscolor(sctmp, nd.mcolor);   /* modified by color */
                if (nd.specfl & SP_TBLT) {
                        scalescolor(sctmp, nd.trans);
                } else {
                        scalescolor(sctmp, nd.tdiff);
                }
                flipsurface(r);
                if (hastexture) {
                        FVECT  bnorm;
                        bnorm[0] = -nd.pnorm[0];
                        bnorm[1] = -nd.pnorm[1];
                        bnorm[2] = -nd.pnorm[2];
                        multambient(sctmp, r, bnorm);
                } else
                        multambient(sctmp, r, r->ron);
                saddscolor(r->rcol, sctmp);
                flipsurface(r);
        }
                                        /* 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, 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;
                                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, 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;
                                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.83
Committed:	Wed Nov 15 18:02:53 2023 UTC (5 months, 3 weeks ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.82:	+68 -68 lines
Log Message:	feat(rpict,rtrace,rcontrib,rtpict): Hyperspectral rendering (except photon map)
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.83	static const char RCSid[] = "$Id: normal.c,v 2.82 2023/03/07 02:24:44 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.1	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	greg	2.83	SCOLOR scval, /* 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	2.83	SCOLOR sctmp;
86	greg	1.3
87	greg	2.83	scolorblack(scval);
88	greg	1.3
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.83	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	greg	2.83	copyscolor(sctmp, np->mcolor);
111	greg	2.49	dtmp = ldot * omega * lrdiff * (1.0/PI);
112	greg	2.83	scalescolor(sctmp, dtmp);
113			saddscolor(scval, sctmp);
114	greg	1.3	}
115	greg	2.72
116	greg	2.83	if ((ldot < -FTINY) & (ltdiff > FTINY)) {
117	greg	2.69	/*
118			* Compute diffuse transmission.
119			*/
120	greg	2.83	copyscolor(sctmp, np->mcolor);
121	greg	2.69	dtmp = -ldot * omega * ltdiff * (1.0/PI);
122	greg	2.83	scalescolor(sctmp, dtmp);
123			saddscolor(scval, sctmp);
124	greg	2.69	}
125	greg	2.72
126			if (ambRayInPmap(np->rp))
127			return; /* specular already in photon map */
128
129	greg	2.83	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	greg	2.83	copyscolor(sctmp, np->scolor);
150	greg	2.54	dtmp = ldot omega;
151	greg	2.83	scalescolor(sctmp, dtmp);
152			saddscolor(scval, sctmp);
153	greg	1.3	}
154			}
155	greg	2.69
156
157	greg	2.83	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	2.83	copyscolor(sctmp, np->mcolor);
169	greg	2.52	dtmp = np->tspec omega * sqrt(-ldot/np->pdot);
170	greg	2.83	scalescolor(sctmp, dtmp);
171			saddscolor(scval, sctmp);
172	greg	1.3	}
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	2.29	int hastexture;
186			double d;
187	greg	2.83	SCOLOR sctmp;
188	greg	2.62	int i;
189	greg	2.69
190			/* PMAP: skip transmitted shadow ray if accounted for in photon map */
191	rschregle	2.74	/* No longer needed?
192	greg	2.73	if (shadowRayInPmap(r) \|\| ambRayInPmap(r))
193	rschregle	2.74	return(1); */
194
195	greg	1.1	/* easy shadow test */
196			if (r->crtype & SHADOW && m->otype != MAT_TRANS)
197	greg	2.27	return(1);
198	greg	2.2
199			if (m->oargs.nfargs != (m->otype == MAT_TRANS ? 7 : 5))
200			objerror(m, USER, "bad number of arguments");
201	greg	2.29	/* check for back side */
202			if (r->rod < 0.0) {
203	greg	2.66	if (!backvis) {
204	greg	2.29	raytrans(r);
205			return(1);
206			}
207	greg	2.40	raytexture(r, m->omod);
208	greg	2.29	flipsurface(r); /* reorient if backvis */
209	greg	2.40	} else
210			raytexture(r, m->omod);
211	greg	1.3	nd.mp = m;
212	greg	2.16	nd.rp = r;
213	greg	1.1	/* get material color */
214	greg	2.83	setscolor(nd.mcolor, m->oargs.farg[0],
215	greg	1.1	m->oargs.farg[1],
216			m->oargs.farg[2]);
217			/* get roughness */
218	greg	2.2	nd.specfl = 0;
219	greg	1.3	nd.alpha2 = m->oargs.farg[4];
220	greg	2.2	if ((nd.alpha2 *= nd.alpha2) <= FTINY)
221			nd.specfl \|= SP_PURE;
222	greg	2.40
223	schorsch	2.45	if ( (hastexture = (DOT(r->pert,r->pert) > FTINY*FTINY)) ) {
224	greg	2.29	nd.pdot = raynormal(nd.pnorm, r); /* perturb normal */
225	greg	2.41	} else {
226	greg	2.29	VCOPY(nd.pnorm, r->ron);
227			nd.pdot = r->rod;
228			}
229	greg	2.42	if (r->ro != NULL && isflat(r->ro->otype))
230			nd.specfl \|= SP_FLAT;
231	greg	1.13	if (nd.pdot < .001)
232			nd.pdot = .001; /* non-zero for dirnorm() */
233	greg	2.83	smultscolor(nd.mcolor, r->pcol); /* modify material color */
234	greg	2.30	nd.rspec = m->oargs.farg[3];
235	greg	2.38	/* compute Fresnel approx. */
236	greg	2.51	if (nd.specfl & SP_PURE && nd.rspec >= FRESTHRESH) {
237	greg	2.62	fest = FRESNE(nd.pdot);
238	greg	2.38	nd.rspec += fest*(1. - nd.rspec);
239			} else
240			fest = 0.;
241	greg	1.3	/* compute transmission */
242	greg	1.1	if (m->otype == MAT_TRANS) {
243	greg	1.3	nd.trans = m->oargs.farg[5]*(1.0 - nd.rspec);
244			nd.tspec = nd.trans * m->oargs.farg[6];
245			nd.tdiff = nd.trans - nd.tspec;
246	greg	2.2	if (nd.tspec > FTINY) {
247			nd.specfl \|= SP_TRAN;
248	greg	2.5	/* check threshold */
249	greg	2.25	if (!(nd.specfl & SP_PURE) &&
250			specthresh >= nd.tspec-FTINY)
251	greg	2.5	nd.specfl \|= SP_TBLT;
252	greg	2.67	if (!hastexture \|\| r->crtype & (SHADOW\|AMBIENT)) {
253	greg	2.2	VCOPY(nd.prdir, r->rdir);
254			} else {
255	greg	2.76	/* perturb */
256			VSUB(nd.prdir, r->rdir, r->pert);
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	2.79	/* diffuse reflection */
266			nd.rdiff = 1.0 - nd.trans - nd.rspec;
267	greg	1.1	/* transmitted ray */
268	greg	2.71	if ((nd.specfl&(SP_TRAN\|SP_PURE\|SP_TBLT)) == (SP_TRAN\|SP_PURE)) {
269	greg	1.3	RAY lr;
270	greg	2.83	copyscolor(lr.rcoef, nd.mcolor); /* modified by color */
271			scalescolor(lr.rcoef, nd.tspec);
272	greg	2.50	if (rayorigin(&lr, TRANS, r, lr.rcoef) == 0) {
273	greg	1.14	VCOPY(lr.rdir, nd.prdir);
274	greg	1.1	rayvalue(&lr);
275	greg	2.83	smultscolor(lr.rcol, lr.rcoef);
276			saddscolor(r->rcol, lr.rcol);
277	greg	2.78	if (nd.tspec >= 1.0-FTINY) {
278			/* completely transparent */
279	greg	2.83	smultscolor(lr.mcol, lr.rcoef);
280			copyscolor(r->mcol, lr.mcol);
281	greg	2.78	r->rmt = r->rot + lr.rmt;
282			r->rxt = r->rot + lr.rxt;
283	greg	2.79	} else if (nd.tspec > nd.tdiff + nd.rdiff)
284	greg	2.78	r->rxt = r->rot + raydistance(&lr);
285	greg	1.1	}
286	greg	2.77	}
287	greg	2.2
288	greg	2.77	if (r->crtype & SHADOW) /* the rest is shadow */
289	greg	2.27	return(1);
290	greg	2.30	/* get specular reflection */
291			if (nd.rspec > FTINY) {
292			nd.specfl \|= SP_REFL;
293			/* compute specular color */
294	greg	2.38	if (m->otype != MAT_METAL) {
295	greg	2.83	setscolor(nd.scolor, nd.rspec, nd.rspec, nd.rspec);
296	greg	2.38	} else if (fest > FTINY) {
297	greg	2.64	d = m->oargs.farg[3]*(1. - fest);
298	greg	2.83	for (i = NCSAMP; i--; )
299			nd.scolor[i] = fest + nd.mcolor[i]*d;
300	greg	2.38	} else {
301	greg	2.83	copyscolor(nd.scolor, nd.mcolor);
302			scalescolor(nd.scolor, nd.rspec);
303	greg	2.38	}
304	greg	2.30	/* check threshold */
305			if (!(nd.specfl & SP_PURE) && specthresh >= nd.rspec-FTINY)
306			nd.specfl \|= SP_RBLT;
307			/* compute reflected ray */
308	greg	2.55	VSUM(nd.vrefl, r->rdir, nd.pnorm, 2.*nd.pdot);
309	greg	2.30	/* penetration? */
310			if (hastexture && DOT(nd.vrefl, r->ron) <= FTINY)
311	greg	2.55	VSUM(nd.vrefl, r->rdir, r->ron, 2.*r->rod);
312	greg	2.53	checknorm(nd.vrefl);
313	gregl	2.36	}
314			/* reflected ray */
315	greg	2.71	if ((nd.specfl&(SP_REFL\|SP_PURE\|SP_RBLT)) == (SP_REFL\|SP_PURE)) {
316	gregl	2.36	RAY lr;
317	greg	2.50	if (rayorigin(&lr, REFLECTED, r, nd.scolor) == 0) {
318	gregl	2.36	VCOPY(lr.rdir, nd.vrefl);
319			rayvalue(&lr);
320	greg	2.83	smultscolor(lr.rcol, lr.rcoef);
321			copyscolor(r->mcol, lr.rcol);
322			saddscolor(r->rcol, lr.rcol);
323	greg	2.80	r->rmt = r->rot;
324	greg	2.67	if (nd.specfl & SP_FLAT &&
325	greg	2.77	!hastexture \| (r->crtype & AMBIENT))
326	greg	2.80	r->rmt += raydistance(&lr);
327	greg	2.30	}
328	greg	2.29	}
329	greg	1.1
330	greg	2.77	if (nd.specfl & SP_PURE && nd.rdiff <= FTINY && nd.tdiff <= FTINY)
331	greg	2.27	return(1); /* 100% pure specular */
332	greg	2.77
333	greg	2.71	if (!(nd.specfl & SP_PURE))
334			gaussamp(&nd); /* checks BLT flags /
335	greg	2.2
336	greg	1.3	if (nd.rdiff > FTINY) { /* ambient from this side */
337	greg	2.83	copyscolor(sctmp, nd.mcolor); /* modified by material color */
338			scalescolor(sctmp, nd.rdiff);
339	greg	2.61	if (nd.specfl & SP_RBLT) /* add in specular as well? */
340	greg	2.83	saddscolor(sctmp, nd.scolor);
341			multambient(sctmp, r, hastexture ? nd.pnorm : r->ron);
342			saddscolor(r->rcol, sctmp); /* add to returned color */
343	greg	1.2	}
344	greg	1.3	if (nd.tdiff > FTINY) { /* ambient from other side */
345	greg	2.83	copyscolor(sctmp, nd.mcolor); /* modified by color */
346			if (nd.specfl & SP_TBLT) {
347			scalescolor(sctmp, nd.trans);
348			} else {
349			scalescolor(sctmp, nd.tdiff);
350			}
351	greg	1.1	flipsurface(r);
352	greg	2.32	if (hastexture) {
353			FVECT bnorm;
354			bnorm[0] = -nd.pnorm[0];
355			bnorm[1] = -nd.pnorm[1];
356			bnorm[2] = -nd.pnorm[2];
357	greg	2.83	multambient(sctmp, r, bnorm);
358	greg	2.32	} else
359	greg	2.83	multambient(sctmp, r, r->ron);
360			saddscolor(r->rcol, sctmp);
361	greg	1.1	flipsurface(r);
362			}
363	greg	1.3	/* add direct component */
364			direct(r, dirnorm, &nd);
365	greg	2.27
366			return(1);
367	greg	2.2	}
368
369
370	greg	2.38	static void
371	greg	2.54	gaussamp( /* sample Gaussian specular */
372	greg	2.62	NORMDAT *np
373	schorsch	2.47	)
374	greg	2.2	{
375			RAY sr;
376			FVECT u, v, h;
377			double rv[2];
378			double d, sinp, cosp;
379	greg	2.83	SCOLOR scol;
380	greg	2.58	int maxiter, ntrials, nstarget, nstaken;
381	greg	2.62	int i;
382	greg	2.13	/* quick test */
383			if ((np->specfl & (SP_REFL\|SP_RBLT)) != SP_REFL &&
384			(np->specfl & (SP_TRAN\|SP_TBLT)) != SP_TRAN)
385			return;
386	greg	2.2	/* set up sample coordinates */
387	greg	2.70	getperpendicular(u, np->pnorm, rand_samp);
388	greg	2.2	fcross(v, np->pnorm, u);
389			/* compute reflection */
390	greg	2.5	if ((np->specfl & (SP_REFL\|SP_RBLT)) == SP_REFL &&
391	greg	2.63	rayorigin(&sr, SPECULAR, np->rp, np->scolor) == 0) {
392	greg	2.58	nstarget = 1;
393	greg	2.55	if (specjitter > 1.5) { /* multiple samples? */
394	greg	2.63	nstarget = specjitter*np->rp->rweight + .5;
395	greg	2.58	if (sr.rweight <= minweight*nstarget)
396			nstarget = sr.rweight/minweight;
397			if (nstarget > 1) {
398			d = 1./nstarget;
399	greg	2.83	scalescolor(sr.rcoef, d);
400	greg	2.56	sr.rweight *= d;
401	greg	2.55	} else
402	greg	2.58	nstarget = 1;
403	greg	2.55	}
404	greg	2.83	scolorblack(scol);
405	greg	2.60	dimlist[ndims++] = (int)(size_t)np->mp;
406	greg	2.58	maxiter = MAXITER*nstarget;
407			for (nstaken = ntrials = 0; nstaken < nstarget &&
408			ntrials < maxiter; ntrials++) {
409			if (ntrials)
410	greg	2.34	d = frandom();
411			else
412			d = urand(ilhash(dimlist,ndims)+samplendx);
413			multisamp(rv, 2, d);
414			d = 2.0PI rv[0];
415	gwlarson	2.37	cosp = tcos(d);
416			sinp = tsin(d);
417	greg	2.55	if ((0. <= specjitter) & (specjitter < 1.))
418			rv[1] = 1.0 - specjitter*rv[1];
419	greg	2.34	if (rv[1] <= FTINY)
420			d = 1.0;
421			else
422			d = sqrt( np->alpha2 * -log(rv[1]) );
423			for (i = 0; i < 3; i++)
424			h[i] = np->pnorm[i] + d(cospu[i] + sinp*v[i]);
425	greg	2.63	d = -2.0 * DOT(h, np->rp->rdir) / (1.0 + d*d);
426			VSUM(sr.rdir, np->rp->rdir, h, d);
427	greg	2.58	/* sample rejection test */
428	greg	2.63	if ((d = DOT(sr.rdir, np->rp->ron)) <= FTINY)
429	greg	2.55	continue;
430			checknorm(sr.rdir);
431	greg	2.58	if (nstarget > 1) { /* W-G-M-D adjustment */
432	greg	2.59	if (nstaken) rayclear(&sr);
433	greg	2.58	rayvalue(&sr);
434	greg	2.63	d = 2./(1. + np->rp->rod/d);
435	greg	2.83	scalescolor(sr.rcol, d);
436			saddscolor(scol, sr.rcol);
437	greg	2.58	} else {
438			rayvalue(&sr);
439	greg	2.83	smultscolor(sr.rcol, sr.rcoef);
440			saddscolor(np->rp->rcol, sr.rcol);
441	greg	2.34	}
442	greg	2.58	++nstaken;
443			}
444			if (nstarget > 1) { /* final W-G-M-D weighting */
445	greg	2.83	smultscolor(scol, sr.rcoef);
446	greg	2.58	d = (double)nstarget/ntrials;
447	greg	2.83	scalescolor(scol, d);
448			saddscolor(np->rp->rcol, scol);
449	greg	2.34	}
450	greg	2.2	ndims--;
451			}
452			/* compute transmission */
453	greg	2.83	copyscolor(sr.rcoef, np->mcolor); /* modified by color */
454			scalescolor(sr.rcoef, np->tspec);
455	greg	2.8	if ((np->specfl & (SP_TRAN\|SP_TBLT)) == SP_TRAN &&
456	greg	2.63	rayorigin(&sr, SPECULAR, np->rp, sr.rcoef) == 0) {
457	greg	2.58	nstarget = 1;
458	greg	2.55	if (specjitter > 1.5) { /* multiple samples? */
459	greg	2.63	nstarget = specjitter*np->rp->rweight + .5;
460	greg	2.58	if (sr.rweight <= minweight*nstarget)
461			nstarget = sr.rweight/minweight;
462			if (nstarget > 1) {
463			d = 1./nstarget;
464	greg	2.83	scalescolor(sr.rcoef, d);
465	greg	2.56	sr.rweight *= d;
466	greg	2.55	} else
467	greg	2.58	nstarget = 1;
468	greg	2.55	}
469	greg	2.60	dimlist[ndims++] = (int)(size_t)np->mp;
470	greg	2.58	maxiter = MAXITER*nstarget;
471			for (nstaken = ntrials = 0; nstaken < nstarget &&
472			ntrials < maxiter; ntrials++) {
473			if (ntrials)
474	greg	2.34	d = frandom();
475			else
476	greg	2.58	d = urand(ilhash(dimlist,ndims)+samplendx);
477	greg	2.34	multisamp(rv, 2, d);
478			d = 2.0PI rv[0];
479	gwlarson	2.37	cosp = tcos(d);
480			sinp = tsin(d);
481	greg	2.55	if ((0. <= specjitter) & (specjitter < 1.))
482			rv[1] = 1.0 - specjitter*rv[1];
483	greg	2.34	if (rv[1] <= FTINY)
484			d = 1.0;
485			else
486	gwlarson	2.37	d = sqrt( np->alpha2 * -log(rv[1]) );
487	greg	2.34	for (i = 0; i < 3; i++)
488			sr.rdir[i] = np->prdir[i] + d(cospu[i] + sinp*v[i]);
489	greg	2.58	/* sample rejection test */
490	greg	2.63	if (DOT(sr.rdir, np->rp->ron) >= -FTINY)
491	greg	2.55	continue;
492			normalize(sr.rdir); /* OK, normalize */
493	greg	2.59	if (nstaken) /* multi-sampling */
494	greg	2.55	rayclear(&sr);
495			rayvalue(&sr);
496	greg	2.83	smultscolor(sr.rcol, sr.rcoef);
497			saddscolor(np->rp->rcol, sr.rcol);
498	greg	2.58	++nstaken;
499	greg	2.34	}
500	greg	2.8	ndims--;
501			}
502	greg	1.1	}