--- ray/src/rt/raytrace.c 2004/03/30 16:13:01 2.45 +++ ray/src/rt/raytrace.c 2010/10/25 22:57:45 2.62 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: raytrace.c,v 2.45 2004/03/30 16:13:01 schorsch Exp $"; +static const char RCSid[] = "$Id: raytrace.c,v 2.62 2010/10/25 22:57:45 greg Exp $"; #endif /* * raytrace.c - routines for tracing and shading rays. @@ -13,16 +13,17 @@ static const char RCSid[] = "$Id: raytrace.c,v 2.45 20 #include "source.h" #include "otypes.h" #include "otspecial.h" +#include "random.h" #define MAXCSET ((MAXSET+1)*2-1) /* maximum check set size */ -unsigned long raynum = 0; /* next unique ray number */ -unsigned long nrays = 0; /* number of calls to localhit */ +RNUMBER raynum = 0; /* next unique ray number */ +RNUMBER nrays = 0; /* number of calls to localhit */ static RREAL Lambfa[5] = {PI, PI, PI, 0.0, 0.0}; OBJREC Lamb = { OVOID, MAT_PLASTIC, "Lambertian", - {0, 5, NULL, Lambfa}, NULL, + {NULL, Lambfa, 0, 5}, NULL }; /* a Lambertian surface */ OBJREC Aftplane; /* aft clipping plane object */ @@ -36,14 +37,22 @@ static void checkset(OBJECT *os, OBJECT *cs); extern int rayorigin( /* start new ray from old one */ - register RAY *r, - register RAY *ro, + RAY *r, int rt, - double rw + const RAY *ro, + const COLOR rc ) { - double re; - + double rw, re; + /* assign coefficient/weight */ + if (rc == NULL) { + rw = 1.0; + setcolor(r->rcoef, 1., 1., 1.); + } else { + rw = intens(rc); + if (rc != r->rcoef) + copycolor(r->rcoef, rc); + } if ((r->parent = ro) == NULL) { /* primary ray */ r->rlvl = 0; r->rweight = rw; @@ -56,6 +65,10 @@ rayorigin( /* start new ray from old one */ r->gecc = seccg; r->slights = NULL; } else { /* spawned ray */ + if (ro->rot >= FHUGE) { + memset(r, 0, sizeof(RAY)); + return(-1); /* illegal continuation */ + } r->rlvl = ro->rlvl; if (rt & RAYREFL) { r->rlvl++; @@ -75,15 +88,39 @@ rayorigin( /* start new ray from old one */ r->crtype = ro->crtype | (r->rtype = rt); VCOPY(r->rorg, ro->rop); r->rweight = ro->rweight * rw; - /* estimate absorption */ + /* estimate extinction */ re = colval(ro->cext,RED) < colval(ro->cext,GRN) ? colval(ro->cext,RED) : colval(ro->cext,GRN); if (colval(ro->cext,BLU) < re) re = colval(ro->cext,BLU); - if (re > 0.) - r->rweight *= exp(-re*ro->rot); + re *= ro->rot; + if (re > 0.1) { + if (re > 92.) { + r->rweight = 0.0; + } else { + r->rweight *= exp(-re); + } + } } rayclear(r); - return(r->rlvl <= maxdepth && r->rweight >= minweight ? 0 : -1); + if (r->rweight <= 0.0) /* check for expiration */ + return(-1); + if (r->crtype & SHADOW) /* shadow commitment */ + return(0); + if (maxdepth <= 0 && rc != NULL) { /* Russian roulette */ + if (minweight <= 0.0) + error(USER, "zero ray weight in Russian roulette"); + if (maxdepth < 0 && r->rlvl > -maxdepth) + return(-1); /* upper reflection limit */ + if (r->rweight >= minweight) + return(0); + if (frandom() > r->rweight/minweight) + return(-1); + rw = minweight/r->rweight; /* promote survivor */ + scalecolor(r->rcoef, rw); + r->rweight = minweight; + return(0); + } + return(r->rlvl <= abs(maxdepth) && r->rweight >= minweight ? 0 : -1); } @@ -119,10 +156,10 @@ raytrace( /* trace a ray and compute its value */ } else if (sourcehit(r)) rayshade(r, r->ro->omod); /* distant source */ - rayparticipate(r); /* for participating medium */ - if (trace != NULL) (*trace)(r); /* trace execution */ + + rayparticipate(r); /* for participating medium */ } @@ -144,7 +181,7 @@ raytrans( /* transmit ray as is */ { RAY tr; - if (rayorigin(&tr, r, TRANS, 1.0) == 0) { + if (rayorigin(&tr, TRANS, r, NULL) == 0) { VCOPY(tr.rdir, r->rdir); rayvalue(&tr); copycolor(r->rcol, tr.rcol); @@ -159,10 +196,10 @@ rayshade( /* shade ray r with material mod */ int mod ) { - int gotmat; register OBJREC *m; + r->rt = r->rot; /* set effective ray length */ - for (gotmat = 0; !gotmat && mod != OVOID; mod = m->omod) { + for ( ; mod != OVOID; mod = m->omod) { m = objptr(mod); /****** unnecessary test since modifier() is always called if (!ismodifier(m->otype)) { @@ -181,10 +218,10 @@ rayshade( /* shade ray r with material mod */ if (!islight(m->otype)) m = &Lamb; } - /* materials call raytexture */ - gotmat = (*ofun[m->otype].funp)(m, r); + if ((*ofun[m->otype].funp)(m, r)) + return(1); /* materials call raytexture() */ } - return(gotmat); + return(0); /* no material! */ } @@ -206,10 +243,10 @@ rayparticipate( /* compute ray medium participation ge *= 1. - colval(r->albedo,GRN); be *= 1. - colval(r->albedo,BLU); } - setcolor(ce, re<=0. ? 1. : re>92. ? 0. : exp(-re), - ge<=0. ? 1. : ge>92. ? 0. : exp(-ge), - be<=0. ? 1. : be>92. ? 0. : exp(-be)); - multcolor(r->rcol, ce); /* path absorption */ + setcolor(ce, re<=FTINY ? 1. : re>92. ? 0. : exp(-re), + ge<=FTINY ? 1. : ge>92. ? 0. : exp(-ge), + be<=FTINY ? 1. : be>92. ? 0. : exp(-be)); + multcolor(r->rcol, ce); /* path extinction */ if (r->crtype & SHADOW || intens(r->albedo) <= FTINY) return; /* no scattering */ setcolor(ca, @@ -266,12 +303,18 @@ raymixture( /* mix modifiers */ foremat = backmat = 0; /* foreground */ fr = *r; - if (coef > FTINY) + if (coef > FTINY) { + fr.rweight *= coef; + scalecolor(fr.rcoef, coef); foremat = rayshade(&fr, fore); + } /* background */ br = *r; - if (coef < 1.0-FTINY) + if (coef < 1.0-FTINY) { + br.rweight *= 1.0-coef; + scalecolor(br.rcoef, 1.0-coef); backmat = rayshade(&br, back); + } /* check for transparency */ if (backmat ^ foremat) { if (backmat && coef > FTINY) @@ -302,7 +345,7 @@ raymixture( /* mix modifiers */ extern double raydist( /* compute (cumulative) ray distance */ - register RAY *r, + register const RAY *r, register int flags ) { @@ -316,6 +359,32 @@ raydist( /* compute (cumulative) ray distance */ } +extern void +raycontrib( /* compute (cumulative) ray contribution */ + double rc[3], + const RAY *r, + int flags +) +{ + double eext[3]; + int i; + + eext[0] = eext[1] = eext[2] = 0.; + rc[0] = rc[1] = rc[2] = 1.; + + while (r != NULL && r->crtype&flags) { + for (i = 3; i--; ) { + rc[i] *= colval(r->rcoef,i); + eext[i] += r->rot * colval(r->cext,i); + } + r = r->parent; + } + for (i = 3; i--; ) + rc[i] *= (eext[i] <= FTINY) ? 1. : + (eext[i] > 92.) ? 0. : exp(-eext[i]); +} + + extern double raynormal( /* compute perturbed normal for ray */ FVECT norm, @@ -362,7 +431,7 @@ newrayxf( /* get new tranformation matrix for ray */ FULLXF xf; } xfseed = { &xfseed }, *xflast = &xfseed; register struct xfn *xp; - register RAY *rp; + register const RAY *rp; /* * Search for transform in circular list that @@ -443,14 +512,15 @@ localhit( /* check for hit in the octree */ else if (r->rdir[i] < -1e-7) sflags |= 0x10 << i; } - if (sflags == 0) - error(CONSISTENCY, "zero ray direction in localhit"); + if (!sflags) { + error(WARNING, "zero ray direction in localhit"); + return(0); + } /* start off assuming nothing hit */ if (r->rmax > FTINY) { /* except aft plane if one */ r->ro = &Aftplane; r->rot = r->rmax; - for (i = 0; i < 3; i++) - r->rop[i] = r->rorg[i] + r->rot*r->rdir[i]; + VSUM(r->rop, r->rorg, r->rdir, r->rot); } /* find global cube entrance point */ t = 0.0; @@ -473,8 +543,7 @@ localhit( /* check for hit in the octree */ if (t >= r->rot) /* clipped already */ return(0); /* advance position */ - for (i = 0; i < 3; i++) - curpos[i] += r->rdir[i]*t; + VSUM(curpos, curpos, r->rdir, t); if (!incube(scene, curpos)) /* non-intersecting ray */ return(0); @@ -565,9 +634,7 @@ raymove( /* check for hit as we move */ ax = 2; } } - pos[0] += r->rdir[0]*t; - pos[1] += r->rdir[1]*t; - pos[2] += r->rdir[2]*t; + VSUM(pos, pos, r->rdir, t); return(ax); }