--- ray/src/rt/raytrace.c 2003/03/11 17:08:55 2.37 +++ ray/src/rt/raytrace.c 2023/11/15 18:02:53 2.88 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: raytrace.c,v 2.37 2003/03/11 17:08:55 greg Exp $"; +static const char RCSid[] = "$Id: raytrace.c,v 2.88 2023/11/15 18:02:53 greg Exp $"; #endif /* * raytrace.c - routines for tracing and shading rays. @@ -10,42 +10,52 @@ static const char RCSid[] = "$Id: raytrace.c,v 2.37 20 #include "copyright.h" #include "ray.h" - +#include "source.h" #include "otypes.h" - #include "otspecial.h" +#include "random.h" +#include "pmap.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 FLOAT Lambfa[5] = {PI, PI, PI, 0.0, 0.0}; +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 */ -static int raymove(), checkhit(); -static void checkset(); - -#ifndef MAXLOOP -#define MAXLOOP 0 /* modifier loop detection */ -#endif - #define RAYHIT (-1) /* return value for intercepted ray */ +static int raymove(FVECT pos, OBJECT *cxs, int dirf, RAY *r, CUBE *cu); +static int checkhit(RAY *r, CUBE *cu, OBJECT *cxs); +static void checkset(OBJECT *os, OBJECT *cs); + int -rayorigin(r, ro, rt, rw) /* start new ray from old one */ -register RAY *r, *ro; -int rt; -double rw; +rayorigin( /* start new ray from old one */ + RAY *r, + int rt, + const RAY *ro, + const SCOLOR rc +) { - double re; - + double rw, re; + /* assign coefficient/weight */ + if (rc == NULL) { + rw = 1.0; + setscolor(r->rcoef, 1., 1., 1.); + } else { + rw = sintens(rc); + if (rw > 1.0) + rw = 1.0; /* avoid calculation growth */ + if (rc != r->rcoef) + copyscolor(r->rcoef, rc); + } if ((r->parent = ro) == NULL) { /* primary ray */ r->rlvl = 0; r->rweight = rw; @@ -58,6 +68,10 @@ double rw; r->gecc = seccg; r->slights = NULL; } else { /* spawned ray */ + if (ro->rot >= FHUGE*.99) { + memset(r, 0, sizeof(RAY)); + return(-1); /* illegal continuation */ + } r->rlvl = ro->rlvl; if (rt & RAYREFL) { r->rlvl++; @@ -77,21 +91,53 @@ double rw; 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); + /* ambient in photon map? */ + if (ro != NULL && ro->crtype & AMBIENT) { + if (causticPhotonMapping) + return(-1); + if (photonMapping && rt != TRANS) + return(-1); + } + 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 */ + scalescolor(r->rcoef, rw); + r->rweight = minweight; + return(0); + } + return((r->rweight >= minweight) & (r->rlvl <= abs(maxdepth)) ? 0 : -1); } void -rayclear(r) /* clear a ray for (re)evaluation */ -register RAY *r; +rayclear( /* clear a ray for (re)evaluation */ + RAY *r +) { r->rno = raynum++; r->newcset = r->clipset; @@ -99,17 +145,23 @@ register RAY *r; r->robj = OVOID; r->ro = NULL; r->rox = NULL; - r->rt = r->rot = FHUGE; + r->rxt = r->rmt = r->rot = FHUGE; + VCOPY(r->rop, r->rorg); + r->ron[0] = -r->rdir[0]; r->ron[1] = -r->rdir[1]; r->ron[2] = -r->rdir[2]; + r->rod = 1.0; r->pert[0] = r->pert[1] = r->pert[2] = 0.0; + r->rflips = 0; r->uv[0] = r->uv[1] = 0.0; - setcolor(r->pcol, 1.0, 1.0, 1.0); - setcolor(r->rcol, 0.0, 0.0, 0.0); + setscolor(r->pcol, 1.0, 1.0, 1.0); + scolorblack(r->mcol); + scolorblack(r->rcol); } void -raytrace(r) /* trace a ray and compute its value */ -RAY *r; +raytrace( /* trace a ray and compute its value */ + RAY *r +) { if (localhit(r, &thescene)) raycont(r); /* hit local surface, evaluate */ @@ -119,16 +171,17 @@ RAY *r; } 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 */ } void -raycont(r) /* check for clipped object and continue */ -register RAY *r; +raycont( /* check for clipped object and continue */ + RAY *r +) { if ((r->clipset != NULL && inset(r->clipset, r->ro->omod)) || !rayshade(r, r->ro->omod)) @@ -137,35 +190,53 @@ register RAY *r; void -raytrans(r) /* transmit ray as is */ -register RAY *r; +raytrans( /* transmit ray as is */ + RAY *r +) { RAY tr; - if (rayorigin(&tr, r, TRANS, 1.0) == 0) { - VCOPY(tr.rdir, r->rdir); - rayvalue(&tr); - copycolor(r->rcol, tr.rcol); - r->rt = r->rot + tr.rt; + rayorigin(&tr, TRANS, r, NULL); /* always continue */ + VCOPY(tr.rdir, r->rdir); + rayvalue(&tr); + copyscolor(r->mcol, tr.mcol); + copyscolor(r->rcol, tr.rcol); + r->rmt = r->rot + tr.rmt; + r->rxt = r->rot + tr.rxt; +} + + +int +raytirrad( /* irradiance hack */ + OBJREC *m, + RAY *r +) +{ + if (ofun[m->otype].flags & (T_M|T_X) && m->otype != MAT_CLIP) { + if (istransp(m->otype) || isBSDFproxy(m)) { + raytrans(r); + return(1); + } + if (!islight(m->otype)) { + setscolor(r->pcol, 1.0, 1.0, 1.0); + return((*ofun[Lamb.otype].funp)(&Lamb, r)); + } } + return(0); /* not a qualifying surface */ } int -rayshade(r, mod) /* shade ray r with material mod */ -register RAY *r; -int mod; +rayshade( /* shade ray r with material mod */ + RAY *r, + int mod +) { - int gotmat; - register OBJREC *m; -#if MAXLOOP - static int depth = 0; - /* check for infinite loop */ - if (depth++ >= MAXLOOP) - objerror(r->ro, USER, "possible modifier loop"); -#endif - r->rt = r->rot; /* set effective ray length */ - for (gotmat = 0; !gotmat && mod != OVOID; mod = m->omod) { + int tst_irrad = do_irrad && !(r->crtype & ~(PRIMARY|TRANS)); + OBJREC *m; + + r->rxt = r->rot; /* preset effective ray length */ + for ( ; mod != OVOID; mod = m->omod) { m = objptr(mod); /****** unnecessary test since modifier() is always called if (!ismodifier(m->otype)) { @@ -174,32 +245,22 @@ int mod; } ******/ /* hack for irradiance calculation */ - if (do_irrad && !(r->crtype & ~(PRIMARY|TRANS))) { - if (irr_ignore(m->otype)) { -#if MAXLOOP - depth--; -#endif - raytrans(r); - return(1); - } - if (!islight(m->otype)) - m = &Lamb; - } - /* materials call raytexture */ - gotmat = (*ofun[m->otype].funp)(m, r); + if (tst_irrad && raytirrad(m, r)) + return(1); + + if ((*ofun[m->otype].funp)(m, r)) + return(1); /* materials call raytexture() */ } -#if MAXLOOP - depth--; -#endif - return(gotmat); + return(0); /* no material! */ } void -rayparticipate(r) /* compute ray medium participation */ -register RAY *r; +rayparticipate( /* compute ray medium participation */ + RAY *r +) { - COLOR ce, ca; + SCOLOR ce, ca; double re, ge, be; if (intens(r->cext) <= 1./FHUGE) @@ -212,32 +273,33 @@ register RAY *r; 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 */ + setscolor(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)); + smultscolor(r->rcol, ce); /* path extinction */ if (r->crtype & SHADOW || intens(r->albedo) <= FTINY) return; /* no scattering */ - setcolor(ca, - colval(r->albedo,RED)*colval(ambval,RED)*(1.-colval(ce,RED)), - colval(r->albedo,GRN)*colval(ambval,GRN)*(1.-colval(ce,GRN)), - colval(r->albedo,BLU)*colval(ambval,BLU)*(1.-colval(ce,BLU))); - addcolor(r->rcol, ca); /* ambient in scattering */ + + /* PMAP: indirect inscattering accounted for by volume photons? */ + if (!volumePhotonMapping) { + setscolor(ca, + colval(r->albedo,RED)*colval(ambval,RED)*(1.-colval(ce,RED)), + colval(r->albedo,GRN)*colval(ambval,GRN)*(1.-colval(ce,GRN)), + colval(r->albedo,BLU)*colval(ambval,BLU)*(1.-colval(ce,BLU))); + saddscolor(r->rcol, ca); /* ambient in scattering */ + } + srcscatter(r); /* source in scattering */ } -raytexture(r, mod) /* get material modifiers */ -RAY *r; -int mod; +void +raytexture( /* get material modifiers */ + RAY *r, + OBJECT mod +) { - register OBJREC *m; -#if MAXLOOP - static int depth = 0; - /* check for infinite loop */ - if (depth++ >= MAXLOOP) - objerror(r->ro, USER, "modifier loop"); -#endif + OBJREC *m; /* execute textures and patterns */ for ( ; mod != OVOID; mod = m->omod) { m = objptr(mod); @@ -253,21 +315,21 @@ int mod; objerror(r->ro, USER, errmsg); } } -#if MAXLOOP - depth--; /* end here */ -#endif } int -raymixture(r, fore, back, coef) /* mix modifiers */ -register RAY *r; -OBJECT fore, back; -double coef; +raymixture( /* mix modifiers */ + RAY *r, + OBJECT fore, + OBJECT back, + double coef +) { RAY fr, br; + double mfore, mback; int foremat, backmat; - register int i; + int i; /* bound coefficient */ if (coef > 1.0) coef = 1.0; @@ -276,44 +338,59 @@ double coef; /* compute foreground and background */ foremat = backmat = 0; /* foreground */ - copystruct(&fr, r); - if (coef > FTINY) + fr = *r; + if (coef > FTINY) { + fr.rweight *= coef; + scalescolor(fr.rcoef, coef); foremat = rayshade(&fr, fore); + } /* background */ - copystruct(&br, r); - if (coef < 1.0-FTINY) + br = *r; + if (coef < 1.0-FTINY) { + br.rweight *= 1.0-coef; + scalescolor(br.rcoef, 1.0-coef); backmat = rayshade(&br, back); + } /* check for transparency */ - if (backmat ^ foremat) + if (backmat ^ foremat) { if (backmat && coef > FTINY) raytrans(&fr); else if (foremat && coef < 1.0-FTINY) raytrans(&br); + } /* mix perturbations */ for (i = 0; i < 3; i++) r->pert[i] = coef*fr.pert[i] + (1.0-coef)*br.pert[i]; /* mix pattern colors */ - scalecolor(fr.pcol, coef); - scalecolor(br.pcol, 1.0-coef); - copycolor(r->pcol, fr.pcol); - addcolor(r->pcol, br.pcol); + scalescolor(fr.pcol, coef); + scalescolor(br.pcol, 1.0-coef); + copyscolor(r->pcol, fr.pcol); + saddscolor(r->pcol, br.pcol); /* return value tells if material */ if (!foremat & !backmat) return(0); /* mix returned ray values */ - scalecolor(fr.rcol, coef); - scalecolor(br.rcol, 1.0-coef); - copycolor(r->rcol, fr.rcol); - addcolor(r->rcol, br.rcol); - r->rt = bright(fr.rcol) > bright(br.rcol) ? fr.rt : br.rt; + scalescolor(fr.rcol, coef); + scalescolor(br.rcol, 1.0-coef); + copyscolor(r->rcol, fr.rcol); + saddscolor(r->rcol, br.rcol); + scalescolor(fr.mcol, coef); + scalescolor(br.mcol, 1.0-coef); + copyscolor(r->mcol, fr.mcol); + saddscolor(r->mcol, br.mcol); + mfore = pbright(fr.mcol); mback = pbright(br.mcol); + r->rmt = mfore > mback ? fr.rmt : br.rmt; + r->rxt = pbright(fr.rcol)-mfore > pbright(br.rcol)-mback ? + fr.rxt : br.rxt; return(1); } double -raydist(r, flags) /* compute (cumulative) ray distance */ -register RAY *r; -register int flags; +raydist( /* compute (cumulative) ray distance */ + const RAY *r, + int flags +) { double sum = 0.0; @@ -325,13 +402,39 @@ register int flags; } +void +raycontrib( /* compute (cumulative) ray contribution */ + SCOLOR rc, + const RAY *r, + int flags +) +{ + static int warnedPM = 0; + + setscolor(rc, 1., 1., 1.); + + while (r != NULL && r->crtype&flags) { + smultscolor(rc, r->rcoef); + /* check for participating medium */ + if (!warnedPM && (bright(r->cext) > FTINY) | + (bright(r->albedo) > FTINY)) { + error(WARNING, + "ray contribution calculation does not support participating media"); + warnedPM++; + } + r = r->parent; + } +} + + double -raynormal(norm, r) /* compute perturbed normal for ray */ -FVECT norm; -register RAY *r; +raynormal( /* compute perturbed normal for ray */ + FVECT norm, + RAY *r +) { double newdot; - register int i; + int i; /* The perturbation is added to the surface normal to obtain * the new normal. If the new normal would affect the surface @@ -361,15 +464,16 @@ register RAY *r; void -newrayxf(r) /* get new tranformation matrix for ray */ -RAY *r; +newrayxf( /* get new tranformation matrix for ray */ + RAY *r +) { static struct xfn { struct xfn *next; FULLXF xf; } xfseed = { &xfseed }, *xflast = &xfseed; - register struct xfn *xp; - register RAY *rp; + struct xfn *xp; + const RAY *rp; /* * Search for transform in circular list that @@ -380,7 +484,7 @@ RAY *r; if (rp->rox == &xp->xf) { /* xp in use */ xp = xp->next; /* move to next */ if (xp == xflast) { /* need new one */ - xp = (struct xfn *)malloc(sizeof(struct xfn)); + xp = (struct xfn *)bmalloc(sizeof(struct xfn)); if (xp == NULL) error(SYSTEM, "out of memory in newrayxf"); @@ -398,8 +502,9 @@ RAY *r; void -flipsurface(r) /* reverse surface orientation */ -register RAY *r; +flipsurface( /* reverse surface orientation */ + RAY *r +) { r->rod = -r->rod; r->ron[0] = -r->ron[0]; @@ -408,13 +513,57 @@ register RAY *r; r->pert[0] = -r->pert[0]; r->pert[1] = -r->pert[1]; r->pert[2] = -r->pert[2]; + r->rflips++; } +int +rayreject( /* check if candidate hit is worse than current */ + OBJREC *o, + RAY *r, + double t, + double rod +) +{ + OBJREC *mnew, *mray; + + if ((t <= FTINY) | (t > r->rot + FTINY)) + return(1); + if (t < r->rot - FTINY) /* is new hit significantly closer? */ + return(0); + /* coincident point, so decide... */ + if (o == r->ro) + return(1); /* shouldn't happen */ + if (r->ro == NULL) + return(0); /* ditto */ + mnew = findmaterial(o); + mray = findmaterial(r->ro); /* check material transparencies */ + if (mnew == NULL) { + if (mray != NULL) + return(1); /* old has material, new does not */ + } else if (mray == NULL) { + return(0); /* new has material, old does not */ + } else if (istransp(mnew->otype)) { + if (!istransp(mray->otype)) + return(1); /* new is transparent, old is not */ + } else if (istransp(mray->otype)) { + return(0); /* old is transparent, new is not */ + } + if (rod <= 0) { /* check which side we hit */ + if (r->rod > 0) + return(1); /* old hit front, new did not */ + } else if (r->rod <= 0) { + return(0); /* new hit front, old did not */ + } + /* earlier modifier definition wins tie */ + return (r->ro->omod >= o->omod); +} + void -rayhit(oset, r) /* standard ray hit test */ -OBJECT *oset; -RAY *r; +rayhit( /* standard ray hit test */ + OBJECT *oset, + RAY *r +) { OBJREC *o; int i; @@ -428,15 +577,16 @@ RAY *r; int -localhit(r, scene) /* check for hit in the octree */ -register RAY *r; -register CUBE *scene; +localhit( /* check for hit in the octree */ + RAY *r, + CUBE *scene +) { OBJECT cxset[MAXCSET+1]; /* set of checked objects */ FVECT curpos; /* current cube position */ int sflags; /* sign flags */ double t, dt; - register int i; + int i; nrays++; /* increment trace counter */ sflags = 0; @@ -447,14 +597,15 @@ register CUBE *scene; 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; @@ -477,32 +628,32 @@ register CUBE *scene; 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); } cxset[0] = 0; raymove(curpos, cxset, sflags, r, scene); - return(r->ro != NULL & r->ro != &Aftplane); + return((r->ro != NULL) & (r->ro != &Aftplane)); } static int -raymove(pos, cxs, dirf, r, cu) /* check for hit as we move */ -FVECT pos; /* current position, modified herein */ -OBJECT *cxs; /* checked objects, modified by checkhit */ -int dirf; /* direction indicators to speed tests */ -register RAY *r; -register CUBE *cu; +raymove( /* check for hit as we move */ + FVECT pos, /* current position, modified herein */ + OBJECT *cxs, /* checked objects, modified by checkhit */ + int dirf, /* direction indicators to speed tests */ + RAY *r, + CUBE *cu +) { int ax; double dt, t; if (istree(cu->cutree)) { /* recurse on subcubes */ CUBE cukid; - register int br, sgn; + int br, sgn; cukid.cusize = cu->cusize * 0.5; /* find subcube */ VCOPY(cukid.cuorg, cu->cuorg); @@ -568,18 +719,17 @@ register CUBE *cu; 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); } static int -checkhit(r, cu, cxs) /* check for hit in full cube */ -register RAY *r; -CUBE *cu; -OBJECT *cxs; +checkhit( /* check for hit in full cube */ + RAY *r, + CUBE *cu, + OBJECT *cxs +) { OBJECT oset[MAXSET+1]; @@ -596,12 +746,13 @@ OBJECT *cxs; static void -checkset(os, cs) /* modify checked set and set to check */ -register OBJECT *os; /* os' = os - cs */ -register OBJECT *cs; /* cs' = cs + os */ +checkset( /* modify checked set and set to check */ + OBJECT *os, /* os' = os - cs */ + OBJECT *cs /* cs' = cs + os */ +) { OBJECT cset[MAXCSET+MAXSET+1]; - register int i, j; + int i, j; int k; /* copy os in place, cset <- cs */ cset[0] = 0;