--- ray/src/rt/srcsamp.c 1991/10/21 12:57:13 1.1 +++ ray/src/rt/srcsamp.c 2003/02/22 02:07:29 2.7 @@ -1,47 +1,101 @@ -/* 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.7 2003/02/22 02:07:29 greg Exp $"; #endif - /* * Source sampling routines + * + * External symbols declared in source.h */ -#include "standard.h" +/* ==================================================================== + * The Radiance Software License, Version 1.0 + * + * Copyright (c) 1990 - 2002 The Regents of the University of California, + * through Lawrence Berkeley National Laboratory. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * + * 3. The end-user documentation included with the redistribution, + * if any, must include the following acknowledgment: + * "This product includes Radiance software + * (http://radsite.lbl.gov/) + * developed by the Lawrence Berkeley National Laboratory + * (http://www.lbl.gov/)." + * Alternately, this acknowledgment may appear in the software itself, + * if and wherever such third-party acknowledgments normally appear. + * + * 4. The names "Radiance," "Lawrence Berkeley National Laboratory" + * and "The Regents of the University of California" must + * not be used to endorse or promote products derived from this + * software without prior written permission. For written + * permission, please contact radiance@radsite.lbl.gov. + * + * 5. Products derived from this software may not be called "Radiance", + * nor may "Radiance" appear in their name, without prior written + * permission of Lawrence Berkeley National Laboratory. + * + * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED + * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR + * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF + * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT + * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * ==================================================================== + * + * This software consists of voluntary contributions made by many + * individuals on behalf of Lawrence Berkeley National Laboratory. For more + * information on Lawrence Berkeley National Laboratory, please see + * . + */ -#include "object.h" +#include "ray.h" #include "source.h" #include "random.h" -extern int dimlist[]; /* dimension list for distribution */ -extern int ndims; /* number of dimensions so far */ -extern int samplendx; /* index for this sample */ +static int cyl_partit(), flt_partit(); double -nextssamp(org, dir, si) /* compute sample for source, rtn. distance */ -FVECT org, dir; /* origin is read only, direction is set */ +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) */ { int cent[3], size[3], parr[2]; 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) - nopart(si, org); + 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; i = source[i].sa.sv.sn) ; /* partition source */ - (*sfun[source[i].so->otype].of->partit)(si, org); + (*sfun[source[i].so->otype].of->partit)(si, r); } si->sp = -1; } @@ -71,41 +125,41 @@ tryagain: vpos[i] += (double)cent[i]/MAXSPART; /* compute direction */ for (i = 0; i < 3; i++) - dir[i] = source[si->sn].sloc[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]; if (!(source[si->sn].sflags & SDISTANT)) for (i = 0; i < 3; i++) - dir[i] -= org[i]; + r->rdir[i] -= r->rorg[i]; /* compute distance */ - if ((d = normalize(dir)) == 0.0) - goto tryagain; /* at source! */ + if ((d = normalize(r->rdir)) == 0.0) + goto nextsample; /* at source! */ /* compute sample size */ - si->dom = source[si->sn].ss2; if (source[si->sn].sflags & SFLAT) { - si->dom *= sflatform(si->sn, dir); - if (si->dom <= FTINY) { /* behind source */ - si->sp = si->np; - goto tryagain; - } - si->dom *= (double)(size[SU]*size[SV])/(MAXSPART*MAXSPART); + si->dom = sflatform(si->sn, r->rdir); + si->dom *= size[SU]*size[SV]/(MAXSPART*(double)MAXSPART); } else if (source[si->sn].sflags & SCYL) { - si->dom *= scylform(si->sn, dir); - si->dom *= (double)size[SU]/MAXSPART; + si->dom = scylform(si->sn, r->rdir); + 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) + if (source[si->sn].sflags & SDISTANT) { + si->dom *= source[si->sn].ss2; return(FHUGE); - si->dom /= d*d; + } + if (si->dom <= 1e-4) + goto nextsample; /* behind source? */ + si->dom *= source[si->sn].ss2/(d*d); return(d); /* sample OK, return distance */ } +int 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) */ @@ -137,9 +191,10 @@ unsigned char *pt; /* partition array */ } -nopart(si, ro) /* single source partition */ +void +nopart(si, r) /* single source partition */ register SRCINDEX *si; -FVECT ro; +RAY *r; { clrpart(si->spt); setpart(si->spt, 0, S0); @@ -147,9 +202,10 @@ FVECT ro; } -cylpart(si, ro) /* partition a cylinder */ +void +cylpart(si, r) /* partition a cylinder */ SRCINDEX *si; -FVECT ro; +register RAY *r; { double dist2, safedist2, dist2cent, rad2; FVECT v; @@ -157,11 +213,11 @@ FVECT ro; int pi; /* first check point location */ clrpart(si->spt); - sp = &source[si->sn]; - rad2 = 1.273 * DOT(sp->ss[SV],sp->ss[SV]); - v[0] = ro[0] - sp->sloc[0]; - v[1] = ro[1] - sp->sloc[1]; - v[2] = ro[2] - sp->sloc[2]; + 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; @@ -171,15 +227,15 @@ FVECT ro; si->np = 0; return; } - safedist2 *= 4./(srcsizerat*srcsizerat); - if (dist2 <= 4.*rad2 || /* point too close to subdivide? */ - dist2cent >= safedist2) { + 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(ro, si->spt, &pi, MAXSPART, + si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART, sp->sloc, sp->ss[SU], safedist2); } @@ -212,32 +268,46 @@ double d2; newct[0] = cent[0] - newax[0]; newct[1] = cent[1] - newax[1]; newct[2] = cent[2] - newax[2]; - npl = cyl_partit(ro, pt, pi, mp*3/4, newct, newax, d2); + npl = cyl_partit(ro, pt, pi, mp/2, newct, newax, d2); /* upper half */ newct[0] = cent[0] + newax[0]; newct[1] = cent[1] + newax[1]; newct[2] = cent[2] + newax[2]; - npu = cyl_partit(ro, pt, pi, mp-npl, newct, newax, d2); + npu = cyl_partit(ro, pt, pi, mp/2, newct, newax, d2); /* return total */ return(npl + npu); } -flatpart(si, ro) /* partition a flat source */ +void +flatpart(si, r) /* partition a flat source */ register SRCINDEX *si; -FVECT ro; +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 = 4./(srcsizerat*srcsizerat) * DOT(vp,vp); + du2 = dv2 * DOT(vp,vp); vp = source[si->sn].ss[SV]; - dv2 = 4./(srcsizerat*srcsizerat) * DOT(vp,vp); - clrpart(si->spt); + dv2 *= DOT(vp,vp); pi = 0; - si->np = flt_partit(ro, si->spt, &pi, MAXSPART, source[si->sn].sloc, + 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); } @@ -282,12 +352,12 @@ double du2, dv2; newct[0] = cent[0] - newax[0]; newct[1] = cent[1] - newax[1]; newct[2] = cent[2] - newax[2]; - npl = flt_partit(ro, pt, pi, mp*3/4, newct, u, v, du2, dv2); + npl = flt_partit(ro, pt, pi, mp/2, newct, u, v, du2, dv2); /* upper half */ newct[0] = cent[0] + newax[0]; newct[1] = cent[1] + newax[1]; newct[2] = cent[2] + newax[2]; - npu = flt_partit(ro, pt, pi, mp-npl, newct, u, v, du2, dv2); + npu = flt_partit(ro, pt, pi, mp/2, newct, u, v, du2, dv2); /* return total */ return(npl + npu); } @@ -298,7 +368,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];