--- ray/src/rt/srcsamp.c 1991/10/22 15:45:21 1.4 +++ ray/src/rt/srcsamp.c 1993/03/08 12:37:39 2.5 @@ -15,6 +15,9 @@ static char SCCSid[] = "$SunId$ LBL"; #include "random.h" +static int cyl_partit(), flt_partit(); + + double nextssamp(r, si) /* compute sample for source, rtn. distance */ register RAY *r; /* origin is read, direction is set */ @@ -24,11 +27,13 @@ register SRCINDEX *si; /* source index (modified to FVECT vpos; double d; register int i; -tryagain: +nextsample: while (++si->sp >= si->np) { /* get next sample */ if (++si->sn >= nsources) return(0.0); /* no more */ - if (srcsizerat <= FTINY) + if (source[si->sn].sflags & SSKIP) + si->np = 0; + else if (srcsizerat <= FTINY) nopart(si, r); else { for (i = si->sn; source[i].sflags & SVIRTUAL; @@ -74,26 +79,24 @@ tryagain: r->rdir[i] -= r->rorg[i]; /* compute distance */ if ((d = normalize(r->rdir)) == 0.0) - goto tryagain; /* 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); - if (si->dom <= FTINY) { /* behind source */ - si->np = 0; - goto tryagain; - } - si->dom *= (double)(size[SU]*size[SV])/(MAXSPART*MAXSPART); + si->dom *= size[SU]*size[SV]/(MAXSPART*(double)MAXSPART); } else if (source[si->sn].sflags & SCYL) { si->dom *= scylform(si->sn, r->rdir); - si->dom *= (double)size[SU]/MAXSPART; + si->dom *= size[SU]/(double)MAXSPART; } else { - si->dom *= (double)(size[SU]*size[SV]*size[SW]) / - (MAXSPART*MAXSPART*MAXSPART) ; + si->dom *= size[SU]*size[SV]*(double)size[SW] / + (MAXSPART*MAXSPART*(double)MAXSPART) ; } if (source[si->sn].sflags & SDISTANT) return(FHUGE); + if (si->dom <= 1e-4) + goto nextsample; /* behind source? */ si->dom /= d*d; return(d); /* sample OK, return distance */ } @@ -165,8 +168,8 @@ register RAY *r; return; } safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat); - if (dist2 <= 4.*rad2 || /* point too close to subdivide? */ - dist2cent >= safedist2) { /* too far? */ + if (dist2 <= 4.*rad2 || /* point too close to subdivide */ + dist2cent >= safedist2) { /* or too far */ setpart(si->spt, 0, S0); si->np = 1; return; @@ -218,19 +221,29 @@ double d2; flatpart(si, r) /* partition a flat source */ register SRCINDEX *si; -RAY *r; +register RAY *r; { - register double *vp; + register FLOAT *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) <= FTINY) { /* 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); - clrpart(si->spt); pi = 0; si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART, source[si->sn].sloc, @@ -294,7 +307,7 @@ scylform(sn, dir) /* compute cosine for cylinder's pr int sn; register FVECT dir; /* assume normalized */ { - register double *dv; + register FLOAT *dv; double d; dv = source[sn].ss[SU];