--- ray/src/rt/srcsamp.c 1991/10/23 13:43:52 1.6 +++ ray/src/rt/srcsamp.c 2014/07/15 23:44:53 2.19 @@ -1,13 +1,14 @@ -/* Copyright (c) 1991 Regents of the University of California */ - #ifndef lint -static char SCCSid[] = "$SunId$ LBL"; +static const char RCSid[] = "$Id: srcsamp.c,v 2.19 2014/07/15 23:44:53 greg Exp $"; #endif - /* * Source sampling routines + * + * External symbols declared in source.h */ +#include "copyright.h" + #include "ray.h" #include "source.h" @@ -15,20 +16,40 @@ static char SCCSid[] = "$SunId$ LBL"; #include "random.h" +static int +srcskip( /* pre-emptive test for out-of-range glow */ + SRCREC *sp, + FVECT orig +) +{ + if (sp->sflags & SSKIP) + return(1); + + if ((sp->sflags & (SPROX|SDISTANT)) != SPROX) + return(0); + + return(dist2(orig, sp->sloc) > + (sp->sl.prox + sp->srad)*(sp->sl.prox + sp->srad)); +} + double -nextssamp(r, si) /* compute sample for source, rtn. distance */ -register RAY *r; /* origin is read, direction is set */ -register SRCINDEX *si; /* source index (modified to current) */ +nextssamp( /* compute sample for source, rtn. distance */ + RAY *r, /* origin is read, direction is set */ + SRCINDEX *si /* source index (modified to current) */\ +) { int cent[3], size[3], parr[2]; - FVECT vpos; + SRCREC *srcp; + double vpos[3]; double d; - register int i; - + int i; +nextsample: while (++si->sp >= si->np) { /* get next sample */ if (++si->sn >= nsources) return(0.0); /* no more */ - if (srcsizerat <= FTINY) + if (srcskip(source+si->sn, r->rorg)) + si->np = 0; + else if (srcsizerat <= FTINY) nopart(si, r); else { for (i = si->sn; source[i].sflags & SVIRTUAL; @@ -45,71 +66,103 @@ register SRCINDEX *si; /* source index (modified to if (!skipparts(cent, size, parr, si->spt)) error(CONSISTENCY, "bad source partition in nextssamp"); /* compute sample */ + srcp = source + si->sn; if (dstrsrc > FTINY) { /* jitter sample */ dimlist[ndims] = si->sn + 8831; dimlist[ndims+1] = si->sp + 3109; d = urand(ilhash(dimlist,ndims+2)+samplendx); - if (source[si->sn].sflags & SFLAT) { + if (srcp->sflags & SFLAT) { multisamp(vpos, 2, d); - vpos[2] = 0.5; + vpos[SW] = 0.5; } else multisamp(vpos, 3, d); for (i = 0; i < 3; i++) vpos[i] = dstrsrc * (1. - 2.*vpos[i]) * - (double)size[i]/MAXSPART; + (double)size[i]*(1.0/MAXSPART); } else vpos[0] = vpos[1] = vpos[2] = 0.0; - for (i = 0; i < 3; i++) - vpos[i] += (double)cent[i]/MAXSPART; + VSUM(vpos, vpos, cent, 1.0/MAXSPART); + /* avoid circular aiming failures */ + if ((srcp->sflags & SCIR) && (si->np > 1 || dstrsrc > 0.7)) { + FVECT trim; + if (srcp->sflags & (SFLAT|SDISTANT)) { + d = 1.12837917; /* correct setflatss() */ + trim[SU] = d*sqrt(1.0 - 0.5*vpos[SV]*vpos[SV]); + trim[SV] = d*sqrt(1.0 - 0.5*vpos[SU]*vpos[SU]); + trim[SW] = 0.0; + } else { + trim[SW] = trim[SU] = vpos[SU]*vpos[SU]; + d = vpos[SV]*vpos[SV]; + if (d > trim[SW]) trim[SW] = d; + trim[SU] += d; + d = vpos[SW]*vpos[SW]; + if (d > trim[SW]) trim[SW] = d; + trim[SU] += d; + if (trim[SU] > FTINY*FTINY) { + d = 1.0/0.7236; /* correct sphsetsrc() */ + trim[SW] = trim[SV] = trim[SU] = + d*sqrt(trim[SW]/trim[SU]); + } else + trim[SW] = trim[SV] = trim[SU] = 0.0; + } + for (i = 0; i < 3; i++) + vpos[i] *= trim[i]; + } /* compute direction */ for (i = 0; i < 3; i++) - r->rdir[i] = source[si->sn].sloc[i] + - vpos[SU]*source[si->sn].ss[SU][i] + - vpos[SV]*source[si->sn].ss[SV][i] + - vpos[SW]*source[si->sn].ss[SW][i]; + r->rdir[i] = srcp->sloc[i] + + vpos[SU]*srcp->ss[SU][i] + + vpos[SV]*srcp->ss[SV][i] + + vpos[SW]*srcp->ss[SW][i]; - if (!(source[si->sn].sflags & SDISTANT)) - for (i = 0; i < 3; i++) - r->rdir[i] -= r->rorg[i]; + if (!(srcp->sflags & SDISTANT)) + VSUB(r->rdir, r->rdir, r->rorg); /* compute distance */ if ((d = normalize(r->rdir)) == 0.0) - return(nextssamp(r, si)); /* at source! */ + goto nextsample; /* at source! */ /* compute sample size */ - si->dom = source[si->sn].ss2; - if (source[si->sn].sflags & SFLAT) { - si->dom *= sflatform(si->sn, r->rdir); - si->dom *= (double)(size[SU]*size[SV])/(MAXSPART*MAXSPART); - } else if (source[si->sn].sflags & SCYL) { - si->dom *= scylform(si->sn, r->rdir); - si->dom *= (double)size[SU]/MAXSPART; + if (srcp->sflags & SFLAT) { + si->dom = sflatform(si->sn, r->rdir); + si->dom *= size[SU]*size[SV]*(1.0/MAXSPART/MAXSPART); + } else if (srcp->sflags & SCYL) { + si->dom = scylform(si->sn, r->rdir); + si->dom *= size[SU]*(1.0/MAXSPART); } else { - si->dom *= (double)(size[SU]*size[SV]*size[SW]) / - (MAXSPART*MAXSPART*MAXSPART) ; + si->dom = size[SU]*size[SV]*(double)size[SW] * + (1.0/MAXSPART/MAXSPART/MAXSPART) ; } - if (source[si->sn].sflags & SDISTANT) + if (srcp->sflags & SDISTANT) { + si->dom *= srcp->ss2; return(FHUGE); - si->dom /= d*d; + } + if (si->dom <= 1e-4) + goto nextsample; /* behind source? */ + si->dom *= srcp->ss2/(d*d); return(d); /* sample OK, return distance */ } -skipparts(ct, sz, pp, pt) /* skip to requested partition */ -int ct[3], sz[3]; /* center and size of partition (returned) */ -register int pp[2]; /* current index, number to skip (modified) */ -unsigned char *pt; /* partition array */ +int +skipparts( /* skip to requested partition */ + int ct[3], + int sz[3], /* center and size of partition (returned) */ + int pp[2], /* current index, number to skip (modified) */ + unsigned char *pt /* partition array */ +) { - register int p; + int p; /* check this partition */ p = spart(pt, pp[0]); pp[0]++; - if (p == S0) /* leaf partition */ + if (p == S0) { /* leaf partition */ if (pp[1]) { pp[1]--; return(0); /* not there yet */ } else return(1); /* we've arrived */ + } /* else check lower */ sz[p] >>= 1; ct[p] -= sz[p]; @@ -126,9 +179,11 @@ unsigned char *pt; /* partition array */ } -nopart(si, r) /* single source partition */ -register SRCINDEX *si; -RAY *r; +void +nopart( /* single source partition */ + SRCINDEX *si, + RAY *r +) { clrpart(si->spt); setpart(si->spt, 0, S0); @@ -136,51 +191,16 @@ RAY *r; } -cylpart(si, r) /* partition a cylinder */ -SRCINDEX *si; -register RAY *r; -{ - double dist2, safedist2, dist2cent, rad2; - FVECT v; - register SRCREC *sp; - int pi; - /* first check point location */ - clrpart(si->spt); - sp = source + si->sn; - rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]); - v[0] = r->rorg[0] - sp->sloc[0]; - v[1] = r->rorg[1] - sp->sloc[1]; - v[2] = r->rorg[2] - sp->sloc[2]; - dist2 = DOT(v,sp->ss[SU]); - safedist2 = DOT(sp->ss[SU],sp->ss[SU]); - dist2 *= dist2 / safedist2; - dist2cent = DOT(v,v); - dist2 = dist2cent - dist2; - if (dist2 <= rad2) { /* point inside extended cylinder */ - si->np = 0; - return; - } - safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat); - if (dist2 <= 4.*rad2 || /* point too close to subdivide */ - dist2cent >= safedist2) { /* or too far */ - setpart(si->spt, 0, S0); - si->np = 1; - return; - } - pi = 0; - si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART, - sp->sloc, sp->ss[SU], safedist2); -} - - static int -cyl_partit(ro, pt, pi, mp, cent, axis, d2) /* slice a cylinder */ -FVECT ro; -unsigned char *pt; -register int *pi; -int mp; -FVECT cent, axis; -double d2; +cyl_partit( /* slice a cylinder */ + FVECT ro, + unsigned char *pt, + int *pi, + int mp, + FVECT cent, + FVECT axis, + double d2 +) { FVECT newct, newax; int npl, npu; @@ -212,46 +232,57 @@ double d2; } -flatpart(si, r) /* partition a flat source */ -register SRCINDEX *si; -register RAY *r; +void +cylpart( /* partition a cylinder */ + SRCINDEX *si, + RAY *r +) { - register FLOAT *vp; + double dist2, safedist2, dist2cent, rad2; FVECT v; - double du2, dv2; + SRCREC *sp; int pi; - + /* first check point location */ clrpart(si->spt); - vp = source[si->sn].sloc; - v[0] = r->rorg[0] - vp[0]; - v[1] = r->rorg[1] - vp[1]; - v[2] = r->rorg[2] - vp[2]; - vp = source[si->sn].snorm; - if (DOT(v,vp) <= FTINY) { /* behind source */ + sp = source + si->sn; + rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]); + v[0] = r->rorg[0] - sp->sloc[0]; + v[1] = r->rorg[1] - sp->sloc[1]; + v[2] = r->rorg[2] - sp->sloc[2]; + dist2 = DOT(v,sp->ss[SU]); + safedist2 = DOT(sp->ss[SU],sp->ss[SU]); + dist2 *= dist2 / safedist2; + dist2cent = DOT(v,v); + dist2 = dist2cent - dist2; + if (dist2 <= rad2) { /* point inside extended cylinder */ si->np = 0; return; } - dv2 = 2.*r->rweight/srcsizerat; - dv2 *= dv2; - vp = source[si->sn].ss[SU]; - du2 = dv2 * DOT(vp,vp); - vp = source[si->sn].ss[SV]; - dv2 *= DOT(vp,vp); + safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat); + if (dist2 <= 4.*rad2 || /* point too close to subdivide */ + dist2cent >= safedist2) { /* or too far */ + setpart(si->spt, 0, S0); + si->np = 1; + return; + } pi = 0; - si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART, - source[si->sn].sloc, - source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2); + si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART, + sp->sloc, sp->ss[SU], safedist2); } static int -flt_partit(ro, pt, pi, mp, cent, u, v, du2, dv2) /* partition flatty */ -FVECT ro; -unsigned char *pt; -register int *pi; -int mp; -FVECT cent, u, v; -double du2, dv2; +flt_partit( /* partition flatty */ + FVECT ro, + unsigned char *pt, + int *pi, + int mp, + FVECT cent, + FVECT u, + FVECT v, + double du2, + double dv2 +) { double d2; FVECT newct, newax; @@ -295,12 +326,47 @@ double du2, dv2; } +void +flatpart( /* partition a flat source */ + SRCINDEX *si, + RAY *r +) +{ + RREAL *vp; + FVECT v; + double du2, dv2; + int pi; + + clrpart(si->spt); + vp = source[si->sn].sloc; + v[0] = r->rorg[0] - vp[0]; + v[1] = r->rorg[1] - vp[1]; + v[2] = r->rorg[2] - vp[2]; + vp = source[si->sn].snorm; + if (DOT(v,vp) <= 0.) { /* behind source */ + si->np = 0; + return; + } + dv2 = 2.*r->rweight/srcsizerat; + dv2 *= dv2; + vp = source[si->sn].ss[SU]; + du2 = dv2 * DOT(vp,vp); + vp = source[si->sn].ss[SV]; + dv2 *= DOT(vp,vp); + pi = 0; + si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART, + source[si->sn].sloc, + source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2); +} + + double -scylform(sn, dir) /* compute cosine for cylinder's projection */ -int sn; -register FVECT dir; /* assume normalized */ +scylform( /* compute cosine for cylinder's projection */ + int sn, + FVECT dir /* assume normalized */ +) { - register FLOAT *dv; + RREAL *dv; double d; dv = source[sn].ss[SU];