--- ray/src/rt/source.c 1991/06/26 14:10:29 1.38 +++ ray/src/rt/source.c 2003/07/27 22:12:03 2.35 @@ -1,53 +1,60 @@ -/* Copyright (c) 1990 Regents of the University of California */ - #ifndef lint -static char SCCSid[] = "$SunId$ LBL"; +static const char RCSid[] = "$Id: source.c,v 2.35 2003/07/27 22:12:03 schorsch Exp $"; #endif - /* * source.c - routines dealing with illumination sources. * - * 8/20/85 + * External symbols declared in source.h */ +#include "copyright.h" + #include "ray.h" -#include "octree.h" - #include "otypes.h" #include "source.h" #include "random.h" +extern double ssampdist; /* scatter sampling distance */ + +#ifndef MAXSSAMP +#define MAXSSAMP 16 /* maximum samples per ray */ +#endif + /* * Structures used by direct() */ typedef struct { + int sno; /* source number */ FVECT dir; /* source direction */ COLOR coef; /* material coefficient */ COLOR val; /* contribution */ } CONTRIB; /* direct contribution */ typedef struct { - int sno; /* source number */ + int sndx; /* source index (to CONTRIB array) */ float brt; /* brightness (for comparison) */ } CNTPTR; /* contribution pointer */ static CONTRIB *srccnt; /* source contributions in direct() */ static CNTPTR *cntord; /* source ordering in direct() */ +static int maxcntr = 0; /* size of contribution arrays */ +void marksources() /* find and mark source objects */ { + int foundsource = 0; int i; register OBJREC *o, *m; register int ns; /* initialize dispatch table */ initstypes(); /* find direct sources */ - for (i = 0; i < nobjects; i++) { + for (i = 0; i < nsceneobjs; i++) { o = objptr(i); @@ -67,6 +74,9 @@ marksources() /* find and mark source objects */ o->otype != OBJ_SOURCE && m->oargs.farg[3] <= FTINY) continue; /* don't bother */ + if (m->oargs.farg[0] <= FTINY && m->oargs.farg[1] <= FTINY && + m->oargs.farg[2] <= FTINY) + continue; /* don't bother */ if (sfun[o->otype].of == NULL || sfun[o->otype].of->setsrc == NULL) @@ -80,7 +90,7 @@ marksources() /* find and mark source objects */ if (m->otype == MAT_GLOW) { source[ns].sflags |= SPROX; source[ns].sl.prox = m->oargs.farg[3]; - if (o->otype == OBJ_SOURCE) + if (source[ns].sflags & SDISTANT) source[ns].sflags |= SSKIP; } else if (m->otype == MAT_SPOT) { source[ns].sflags |= SSPOT; @@ -93,15 +103,19 @@ marksources() /* find and mark source objects */ source[ns].sflags |= SSKIP; } } + if (!(source[ns].sflags & SSKIP)) + foundsource++; } - if (nsources <= 0) { + if (!foundsource) { error(WARNING, "no light sources found"); return; } markvirtuals(); /* find and add virtual sources */ - srccnt = (CONTRIB *)malloc(nsources*sizeof(CONTRIB)); - cntord = (CNTPTR *)malloc(nsources*sizeof(CNTPTR)); - if (srccnt == NULL || cntord == NULL) + /* allocate our contribution arrays */ + maxcntr = nsources + MAXSPART; /* start with this many */ + srccnt = (CONTRIB *)malloc(maxcntr*sizeof(CONTRIB)); + cntord = (CNTPTR *)malloc(maxcntr*sizeof(CNTPTR)); + if ((srccnt == NULL) | (cntord == NULL)) goto memerr; return; memerr: @@ -109,95 +123,65 @@ memerr: } -double -srcray(sr, r, sn) /* send a ray to a source, return domega */ +void +freesources() /* free all source structures */ +{ + if (nsources > 0) { + free((void *)source); + source = NULL; + nsources = 0; + } + if (maxcntr <= 0) + return; + free((void *)srccnt); + srccnt = NULL; + free((void *)cntord); + cntord = NULL; + maxcntr = 0; +} + + +int +srcray(sr, r, si) /* send a ray to a source, return domega */ register RAY *sr; /* returned source ray */ RAY *r; /* ray which hit object */ -register int sn; /* source number */ +SRCINDEX *si; /* source sample index */ { - double ddot; /* (distance times) cosine */ - FVECT vd; - double d; - register int i; + double d; /* distance to source */ + register SRCREC *srcp; - if (source[sn].sflags & SSKIP) - return(0.0); /* skip this source */ + rayorigin(sr, r, SHADOW, 1.0); /* ignore limits */ - rayorigin(sr, r, SHADOW, 1.0); /* ignore limits */ - - sr->rsrc = sn; /* remember source */ - /* get source direction */ - if (source[sn].sflags & SDISTANT) { - if (source[sn].sflags & SSPOT) { /* check location */ - for (i = 0; i < 3; i++) - vd[i] = sr->rorg[i] - source[sn].sl.s->aim[i]; - d = DOT(source[sn].sloc,vd); - d = DOT(vd,vd) - d*d; - if (PI*d > source[sn].sl.s->siz) - return(0.0); - } - /* constant direction */ - VCOPY(sr->rdir, source[sn].sloc); - } else { /* compute direction */ - for (i = 0; i < 3; i++) - sr->rdir[i] = source[sn].sloc[i] - sr->rorg[i]; - - if (source[sn].sflags & SFLAT && - (ddot = -DOT(sr->rdir, source[sn].snorm)) <= FTINY) - return(0.0); /* behind surface! */ + while ((d = nextssamp(sr, si)) != 0.0) { + sr->rsrc = si->sn; /* remember source */ + srcp = source + si->sn; + if (srcp->sflags & SDISTANT) { + if (srcp->sflags & SSPOT && spotout(sr, srcp->sl.s)) + continue; + return(1); /* sample OK */ } - if (dstrsrc > FTINY) { - /* distribute source direction */ - dimlist[ndims++] = sn; - for (i = 0; i < 3; i++) { - dimlist[ndims] = i + 8831; - vd[i] = dstrsrc * source[sn].ss * - (1.0 - 2.0*urand(ilhash(dimlist,ndims+1)+samplendx)); - } - ndims--; - if (source[sn].sflags & SFLAT) { /* project offset */ - d = DOT(vd, source[sn].snorm); - for (i = 0; i < 3; i++) - vd[i] -= d * source[sn].snorm[i]; - } - for (i = 0; i < 3; i++) /* offset source direction */ - sr->rdir[i] += vd[i]; - - } else if (source[sn].sflags & SDISTANT) - /* already normalized */ - return(source[sn].ss2); - - if ((d = normalize(sr->rdir)) == 0.0) - /* at source! */ - return(0.0); - - if (source[sn].sflags & SDISTANT) - /* domega constant */ - return(source[sn].ss2); - + /* local source */ /* check proximity */ - if (source[sn].sflags & SPROX && - d > source[sn].sl.prox) - return(0.0); - /* compute dot product */ - if (source[sn].sflags & SFLAT) - ddot /= d; - else - ddot = 1.0; + if (srcp->sflags & SPROX && d > srcp->sl.prox) + continue; /* check angle */ - if (source[sn].sflags & SSPOT) { - if (source[sn].sl.s->siz < 2.0*PI * - (1.0 + DOT(source[sn].sl.s->aim,sr->rdir))) - return(0.0); - d += source[sn].sl.s->flen; /* adjust length */ + if (srcp->sflags & SSPOT) { + if (spotout(sr, srcp->sl.s)) + continue; + /* adjust solid angle */ + si->dom *= d*d; + d += srcp->sl.s->flen; + si->dom /= d*d; } - /* compute domega */ - return(ddot*source[sn].ss2/(d*d)); + return(1); /* sample OK */ + } + return(0); /* no more samples */ } +void srcvalue(r) /* punch ray to source and compute value */ -RAY *r; +register RAY *r; { register SRCREC *sp; @@ -206,7 +190,9 @@ RAY *r; /* check intersection */ if (!(*ofun[sp->so->otype].funp)(sp->so, r)) return; - raycont(r); /* compute contribution */ + if (!rayshade(r, r->ro->omod)) /* compute contribution */ + goto nomat; + rayparticipate(r); return; } /* compute intersection */ @@ -214,9 +200,12 @@ RAY *r; (*ofun[sp->so->otype].funp)(sp->so, r)) { if (sp->sa.success >= 0) sp->sa.success++; - raycont(r); /* compute contribution */ + if (!rayshade(r, r->ro->omod)) /* compute contribution */ + goto nomat; + rayparticipate(r); return; } + /* we missed our mark! */ if (sp->sa.success < 0) return; /* bitched already */ sp->sa.success -= AIMREQT; @@ -225,9 +214,51 @@ RAY *r; sprintf(errmsg, "aiming failure for light source \"%s\"", sp->so->oname); error(WARNING, errmsg); /* issue warning */ + return; +nomat: + objerror(r->ro, USER, "material not found"); } +int +sourcehit(r) /* check to see if ray hit distant source */ +register RAY *r; +{ + int first, last; + register int i; + + if (r->rsrc >= 0) { /* check only one if aimed */ + first = last = r->rsrc; + } else { /* otherwise check all */ + first = 0; last = nsources-1; + } + for (i = first; i <= last; i++) + if ((source[i].sflags & (SDISTANT|SVIRTUAL)) == SDISTANT) + /* + * Check to see if ray is within + * solid angle of source. + */ + if (2.0*PI * (1.0 - DOT(source[i].sloc,r->rdir)) + <= source[i].ss2) { + r->ro = source[i].so; + if (!(source[i].sflags & SSKIP)) + break; + } + + if (r->ro != NULL) { + r->robj = objndx(r->ro); + for (i = 0; i < 3; i++) + r->ron[i] = -r->rdir[i]; + r->rod = 1.0; + r->pert[0] = r->pert[1] = r->pert[2] = 0.0; + r->uv[0] = r->uv[1] = 0.0; + r->rox = NULL; + return(1); + } + return(0); +} + + static int cntcmp(sc1, sc2) /* contribution compare (descending) */ register CNTPTR *sc1, *sc2; @@ -240,50 +271,58 @@ register CNTPTR *sc1, *sc2; } +void direct(r, f, p) /* add direct component */ RAY *r; /* ray that hit surface */ -int (*f)(); /* direct component coefficient function */ +void (*f)(); /* direct component coefficient function */ char *p; /* data for f */ { - extern double pow(); + extern void (*trace)(); register int sn; + register CONTRIB *scp; + SRCINDEX si; int nshadcheck, ncnts; int nhits; - double dom, prob, ourthresh, hwt; + double prob, ourthresh, hwt; RAY sr; /* NOTE: srccnt and cntord global so no recursion */ if (nsources <= 0) return; /* no sources?! */ - /* compute number to check */ - nshadcheck = pow((double)nsources, shadcert) + .5; - /* modify threshold */ - ourthresh = shadthresh / r->rweight; /* potential contributions */ - for (sn = 0; sn < nsources; sn++) { - cntord[sn].sno = sn; - cntord[sn].brt = 0.0; - /* get source ray */ - if ((dom = srcray(&sr, r, sn)) == 0.0) - continue; - VCOPY(srccnt[sn].dir, sr.rdir); + initsrcindex(&si); + for (sn = 0; srcray(&sr, r, &si); sn++) { + if (sn >= maxcntr) { + maxcntr = sn + MAXSPART; + srccnt = (CONTRIB *)realloc((void *)srccnt, + maxcntr*sizeof(CONTRIB)); + cntord = (CNTPTR *)realloc((void *)cntord, + maxcntr*sizeof(CNTPTR)); + if ((srccnt == NULL) | (cntord == NULL)) + error(SYSTEM, "out of memory in direct"); + } + cntord[sn].sndx = sn; + scp = srccnt + sn; + scp->sno = sr.rsrc; /* compute coefficient */ - (*f)(srccnt[sn].coef, p, srccnt[sn].dir, dom); - cntord[sn].brt = bright(srccnt[sn].coef); + (*f)(scp->coef, p, sr.rdir, si.dom); + cntord[sn].brt = bright(scp->coef); if (cntord[sn].brt <= 0.0) continue; + VCOPY(scp->dir, sr.rdir); /* compute potential */ sr.revf = srcvalue; rayvalue(&sr); - copycolor(srccnt[sn].val, sr.rcol); - multcolor(srccnt[sn].val, srccnt[sn].coef); - cntord[sn].brt = bright(srccnt[sn].val); + copycolor(scp->val, sr.rcol); + multcolor(scp->val, scp->coef); + cntord[sn].brt = bright(scp->val); } /* sort contributions */ - qsort(cntord, nsources, sizeof(CNTPTR), cntcmp); + qsort(cntord, sn, sizeof(CNTPTR), cntcmp); { /* find last */ register int l, m; - sn = 0; ncnts = l = nsources; + ncnts = l = sn; + sn = 0; while ((m = (sn + ncnts) >> 1) != l) { if (cntord[m].brt > 0.0) sn = m; @@ -292,54 +331,262 @@ char *p; /* data for f */ l = m; } } + if (ncnts == 0) + return; /* no contributions! */ /* accumulate tail */ for (sn = ncnts-1; sn > 0; sn--) cntord[sn-1].brt += cntord[sn].brt; + /* compute number to check */ + nshadcheck = pow((double)ncnts, shadcert) + .5; + /* modify threshold */ + ourthresh = shadthresh / r->rweight; /* test for shadows */ - nhits = 0; - for (sn = 0; sn < ncnts; sn++) { + for (nhits = 0, hwt = 0.0, sn = 0; sn < ncnts; + hwt += (double)source[scp->sno].nhits / + (double)source[scp->sno].ntests, + sn++) { /* check threshold */ if ((sn+nshadcheck>=ncnts ? cntord[sn].brt : cntord[sn].brt-cntord[sn+nshadcheck].brt) < ourthresh*bright(r->rcol)) break; - /* get statistics */ - source[cntord[sn].sno].ntests++; + scp = srccnt + cntord[sn].sndx; /* test for hit */ rayorigin(&sr, r, SHADOW, 1.0); - VCOPY(sr.rdir, srccnt[cntord[sn].sno].dir); - sr.rsrc = cntord[sn].sno; + VCOPY(sr.rdir, scp->dir); + sr.rsrc = scp->sno; + /* keep statistics */ + if (source[scp->sno].ntests++ > 0xfffffff0) { + source[scp->sno].ntests >>= 1; + source[scp->sno].nhits >>= 1; + } if (localhit(&sr, &thescene) && - ( sr.ro != source[cntord[sn].sno].so || - source[cntord[sn].sno].sflags & SFOLLOW )) { + ( sr.ro != source[scp->sno].so || + source[scp->sno].sflags & SFOLLOW )) { /* follow entire path */ raycont(&sr); + rayparticipate(&sr); + if (trace != NULL) + (*trace)(&sr); /* trace execution */ if (bright(sr.rcol) <= FTINY) continue; /* missed! */ - copycolor(srccnt[cntord[sn].sno].val, sr.rcol); - multcolor(srccnt[cntord[sn].sno].val, - srccnt[cntord[sn].sno].coef); + copycolor(scp->val, sr.rcol); + multcolor(scp->val, scp->coef); } /* add contribution if hit */ - addcolor(r->rcol, srccnt[cntord[sn].sno].val); + addcolor(r->rcol, scp->val); nhits++; - source[cntord[sn].sno].nhits++; + source[scp->sno].nhits++; } - /* surface hit rate */ - if (sn > 0) - hwt = (double)nhits / (double)sn; + /* source hit rate */ + if (hwt > FTINY) + hwt = (double)nhits / hwt; else hwt = 0.5; #ifdef DEBUG - sprintf(errmsg, "%d tested, %d untested, %f hit rate\n", + sprintf(errmsg, "%d tested, %d untested, %f conditional hit rate\n", sn, ncnts-sn, hwt); eputs(errmsg); #endif /* add in untested sources */ for ( ; sn < ncnts; sn++) { - prob = hwt * (double)source[cntord[sn].sno].nhits / - (double)source[cntord[sn].sno].ntests; - scalecolor(srccnt[cntord[sn].sno].val, prob); - addcolor(r->rcol, srccnt[cntord[sn].sno].val); + scp = srccnt + cntord[sn].sndx; + prob = hwt * (double)source[scp->sno].nhits / + (double)source[scp->sno].ntests; + if (prob > 1.0) + prob = 1.0; + scalecolor(scp->val, prob); + addcolor(r->rcol, scp->val); } +} + + +void +srcscatter(r) /* compute source scattering into ray */ +register RAY *r; +{ + int oldsampndx; + int nsamps; + RAY sr; + SRCINDEX si; + double t, d; + double re, ge, be; + COLOR cvext; + int i, j; + + if (r->slights == NULL || r->slights[0] == 0 + || r->gecc >= 1.-FTINY || r->rot >= FHUGE) + return; + if (ssampdist <= FTINY || (nsamps = r->rot/ssampdist + .5) < 1) + nsamps = 1; +#if MAXSSAMP + else if (nsamps > MAXSSAMP) + nsamps = MAXSSAMP; +#endif + oldsampndx = samplendx; + samplendx = random()&0x7fff; /* randomize */ + for (i = r->slights[0]; i > 0; i--) { /* for each source */ + for (j = 0; j < nsamps; j++) { /* for each sample position */ + samplendx++; + t = r->rot * (j+frandom())/nsamps; + /* extinction */ + re = t*colval(r->cext,RED); + ge = t*colval(r->cext,GRN); + be = t*colval(r->cext,BLU); + setcolor(cvext, re > 92. ? 0. : exp(-re), + ge > 92. ? 0. : exp(-ge), + be > 92. ? 0. : exp(-be)); + if (intens(cvext) <= FTINY) + break; /* too far away */ + sr.rorg[0] = r->rorg[0] + r->rdir[0]*t; + sr.rorg[1] = r->rorg[1] + r->rdir[1]*t; + sr.rorg[2] = r->rorg[2] + r->rdir[2]*t; + sr.rmax = 0.; + initsrcindex(&si); /* sample ray to this source */ + si.sn = r->slights[i]; + nopart(&si, &sr); + if (!srcray(&sr, NULL, &si) || + sr.rsrc != r->slights[i]) + continue; /* no path */ + copycolor(sr.cext, r->cext); + copycolor(sr.albedo, r->albedo); + sr.gecc = r->gecc; + sr.slights = r->slights; + rayvalue(&sr); /* eval. source ray */ + if (bright(sr.rcol) <= FTINY) + continue; + if (r->gecc <= FTINY) /* compute P(theta) */ + d = 1.; + else { + d = DOT(r->rdir, sr.rdir); + d = 1. + r->gecc*r->gecc - 2.*r->gecc*d; + d = (1. - r->gecc*r->gecc) / (d*sqrt(d)); + } + /* other factors */ + d *= si.dom * r->rot / (4.*PI*nsamps); + multcolor(sr.rcol, r->cext); + multcolor(sr.rcol, r->albedo); + scalecolor(sr.rcol, d); + multcolor(sr.rcol, cvext); + addcolor(r->rcol, sr.rcol); /* add it in */ + } + } + samplendx = oldsampndx; +} + + +/**************************************************************** + * The following macros were separated from the m_light() routine + * because they are very nasty and difficult to understand. + */ + +/* illumblock * + * + * We cannot allow an illum to pass to another illum, because that + * would almost certainly constitute overcounting. + * However, we do allow an illum to pass to another illum + * that is actually going to relay to a virtual light source. + * We also prevent an illum from passing to a glow; this provides a + * convenient mechanism for defining detailed light source + * geometry behind (or inside) an effective radiator. + */ + +static int weaksrcmod(obj) int obj; /* efficiency booster function */ +{register OBJREC *o = objptr(obj); +return((o->otype==MAT_ILLUM)|(o->otype==MAT_GLOW));} + +#define illumblock(m, r) (!(source[r->rsrc].sflags&SVIRTUAL) && \ + r->rod > 0.0 && \ + weaksrcmod(source[r->rsrc].so->omod)) + +/* wrongsource * + * + * This source is the wrong source (ie. overcounted) if we are + * aimed to a different source than the one we hit and the one + * we hit is not an illum that should be passed. + */ + +#define wrongsource(m, r) (r->rsrc>=0 && source[r->rsrc].so!=r->ro && \ + (m->otype!=MAT_ILLUM || illumblock(m,r))) + +/* distglow * + * + * A distant glow is an object that sometimes acts as a light source, + * but is too far away from the test point to be one in this case. + * (Glows with negative radii should NEVER participate in illumination.) + */ + +#define distglow(m, r, d) (m->otype==MAT_GLOW && \ + m->oargs.farg[3] >= -FTINY && \ + d > m->oargs.farg[3]) + +/* badcomponent * + * + * We must avoid counting light sources in the ambient calculation, + * since the direct component is handled separately. Therefore, any + * ambient ray which hits an active light source must be discarded. + * The same is true for stray specular samples, since the specular + * contribution from light sources is calculated separately. + */ + +#define badcomponent(m, r) (r->crtype&(AMBIENT|SPECULAR) && \ + !(r->crtype&SHADOW || r->rod < 0.0 || \ + /* not 100% correct */ distglow(m, r, r->rot))) + +/* passillum * + * + * An illum passes to another material type when we didn't hit it + * on purpose (as part of a direct calculation), or it is relaying + * a virtual light source. + */ + +#define passillum(m, r) (m->otype==MAT_ILLUM && \ + (r->rsrc<0 || source[r->rsrc].so!=r->ro || \ + source[r->rsrc].sflags&SVIRTUAL)) + +/* srcignore * + * + * The -dv flag is normally on for sources to be visible. + */ + +#define srcignore(m, r) !(directvis || r->crtype&SHADOW || \ + distglow(m, r, raydist(r,PRIMARY))) + + +int +m_light(m, r) /* ray hit a light source */ +register OBJREC *m; +register RAY *r; +{ + /* check for over-counting */ + if (badcomponent(m, r)) + return(1); + if (wrongsource(m, r)) + return(1); + /* check for passed illum */ + if (passillum(m, r)) { + if (m->oargs.nsargs && strcmp(m->oargs.sarg[0], VOIDID)) + return(rayshade(r,lastmod(objndx(m),m->oargs.sarg[0]))); + raytrans(r); + return(1); + } + /* otherwise treat as source */ + /* check for behind */ + if (r->rod < 0.0) + return(1); + /* check for invisibility */ + if (srcignore(m, r)) + return(1); + /* check for outside spot */ + if (m->otype==MAT_SPOT && spotout(r, makespot(m))) + return(1); + /* get distribution pattern */ + raytexture(r, m->omod); + /* get source color */ + setcolor(r->rcol, m->oargs.farg[0], + m->oargs.farg[1], + m->oargs.farg[2]); + /* modify value */ + multcolor(r->rcol, r->pcol); + return(1); }