Back to index

radiance  4R0+20100331
srcsamp.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: srcsamp.c,v 2.16 2009/06/06 02:11:44 greg Exp $";
00003 #endif
00004 /*
00005  * Source sampling routines
00006  *
00007  *  External symbols declared in source.h
00008  */
00009 
00010 #include "copyright.h"
00011 
00012 #include  "ray.h"
00013 
00014 #include  "source.h"
00015 
00016 #include  "random.h"
00017 
00018 
00019 static int  cyl_partit(), flt_partit();
00020 
00021 
00022 double
00023 nextssamp(r, si)            /* compute sample for source, rtn. distance */
00024 register RAY  *r;           /* origin is read, direction is set */
00025 register SRCINDEX  *si;            /* source index (modified to current) */
00026 {
00027        int  cent[3], size[3], parr[2];
00028        SRCREC  *srcp;
00029        FVECT  vpos;
00030        double  d;
00031        register int  i;
00032 nextsample:
00033        while (++si->sp >= si->np) {       /* get next sample */
00034               if (++si->sn >= nsources)
00035                      return(0.0);  /* no more */
00036               if (source[si->sn].sflags & SSKIP)
00037                      si->np = 0;
00038               else if (srcsizerat <= FTINY)
00039                      nopart(si, r);
00040               else {
00041                      for (i = si->sn; source[i].sflags & SVIRTUAL;
00042                                    i = source[i].sa.sv.sn)
00043                             ;             /* partition source */
00044                      (*sfun[source[i].so->otype].of->partit)(si, r);
00045               }
00046               si->sp = -1;
00047        }
00048                                    /* get partition */
00049        cent[0] = cent[1] = cent[2] = 0;
00050        size[0] = size[1] = size[2] = MAXSPART;
00051        parr[0] = 0; parr[1] = si->sp;
00052        if (!skipparts(cent, size, parr, si->spt))
00053               error(CONSISTENCY, "bad source partition in nextssamp");
00054                                    /* compute sample */
00055        srcp = source + si->sn;
00056        if (dstrsrc > FTINY) {                    /* jitter sample */
00057               dimlist[ndims] = si->sn + 8831;
00058               dimlist[ndims+1] = si->sp + 3109;
00059               d = urand(ilhash(dimlist,ndims+2)+samplendx);
00060               if (srcp->sflags & SFLAT) {
00061                      multisamp(vpos, 2, d);
00062                      vpos[SW] = 0.5;
00063               } else
00064                      multisamp(vpos, 3, d);
00065               for (i = 0; i < 3; i++)
00066                      vpos[i] = dstrsrc * (1. - 2.*vpos[i]) *
00067                                    (double)size[i]/MAXSPART;
00068        } else
00069               vpos[0] = vpos[1] = vpos[2] = 0.0;
00070 
00071        for (i = 0; i < 3; i++)
00072               vpos[i] += (double)cent[i]/MAXSPART;
00073                                    /* avoid circular aiming failures */
00074        if ((srcp->sflags & SCIR) && (si->np > 1 || dstrsrc > 0.7)) {
00075               FVECT  trim;
00076               if (srcp->sflags & (SFLAT|SDISTANT)) {
00077                      d = 1.12837917;             /* correct setflatss() */
00078                      trim[SU] = d*sqrt(1.0 - 0.5*vpos[SV]*vpos[SV]);
00079                      trim[SV] = d*sqrt(1.0 - 0.5*vpos[SU]*vpos[SU]);
00080                      trim[SW] = 0.0;
00081               } else {
00082                      trim[SW] = trim[SU] = vpos[SU]*vpos[SU];
00083                      d = vpos[SV]*vpos[SV];
00084                      if (d > trim[SW]) trim[SW] = d;
00085                      trim[SU] += d;
00086                      d = vpos[SW]*vpos[SW];
00087                      if (d > trim[SW]) trim[SW] = d;
00088                      trim[SU] += d;
00089                      if (trim[SU] > FTINY*FTINY) {
00090                             d = 1.0/0.7236;      /* correct sphsetsrc() */
00091                             trim[SW] = trim[SV] = trim[SU] =
00092                                           d*sqrt(trim[SW]/trim[SU]);
00093                      } else
00094                             trim[SW] = trim[SV] = trim[SU] = 0.0;
00095               }
00096               for (i = 0; i < 3; i++)
00097                      vpos[i] *= trim[i];
00098        }
00099                                    /* compute direction */
00100        for (i = 0; i < 3; i++)
00101               r->rdir[i] = srcp->sloc[i] +
00102                             vpos[SU]*srcp->ss[SU][i] +
00103                             vpos[SV]*srcp->ss[SV][i] +
00104                             vpos[SW]*srcp->ss[SW][i];
00105 
00106        if (!(srcp->sflags & SDISTANT))
00107               for (i = 0; i < 3; i++)
00108                      r->rdir[i] -= r->rorg[i];
00109                                    /* compute distance */
00110        if ((d = normalize(r->rdir)) == 0.0)
00111               goto nextsample;            /* at source! */
00112 
00113                                    /* compute sample size */
00114        if (srcp->sflags & SFLAT) {
00115               si->dom = sflatform(si->sn, r->rdir);
00116               si->dom *= size[SU]*size[SV]/(MAXSPART*(double)MAXSPART);
00117        } else if (srcp->sflags & SCYL) {
00118               si->dom = scylform(si->sn, r->rdir);
00119               si->dom *= size[SU]/(double)MAXSPART;
00120        } else {
00121               si->dom = size[SU]*size[SV]*(double)size[SW] /
00122                             (MAXSPART*MAXSPART*(double)MAXSPART) ;
00123        }
00124        if (srcp->sflags & SDISTANT) {
00125               si->dom *= srcp->ss2;
00126               return(FHUGE);
00127        }
00128        if (si->dom <= 1e-4)
00129               goto nextsample;            /* behind source? */
00130        si->dom *= srcp->ss2/(d*d);
00131        return(d);           /* sample OK, return distance */
00132 }
00133 
00134 
00135 int
00136 skipparts(ct, sz, pp, pt)          /* skip to requested partition */
00137 int  ct[3], sz[3];          /* center and size of partition (returned) */
00138 register int  pp[2];        /* current index, number to skip (modified) */
00139 unsigned char  *pt;         /* partition array */
00140 {
00141        register int  p;
00142                                    /* check this partition */
00143        p = spart(pt, pp[0]);
00144        pp[0]++;
00145        if (p == S0) {                     /* leaf partition */
00146               if (pp[1]) {
00147                      pp[1]--;
00148                      return(0);    /* not there yet */
00149               } else
00150                      return(1);    /* we've arrived */
00151        }
00152                             /* else check lower */
00153        sz[p] >>= 1;
00154        ct[p] -= sz[p];
00155        if (skipparts(ct, sz, pp, pt))
00156               return(1);    /* return hit */
00157                             /* else check upper */
00158        ct[p] += sz[p] << 1;
00159        if (skipparts(ct, sz, pp, pt))
00160               return(1);    /* return hit */
00161                             /* else return to starting position */
00162        ct[p] -= sz[p];
00163        sz[p] <<= 1;
00164        return(0);           /* return miss */
00165 }
00166 
00167 
00168 void
00169 nopart(si, r)               /* single source partition */
00170 register SRCINDEX  *si;
00171 RAY  *r;
00172 {
00173        clrpart(si->spt);
00174        setpart(si->spt, 0, S0);
00175        si->np = 1;
00176 }
00177 
00178 
00179 void
00180 cylpart(si, r)                     /* partition a cylinder */
00181 SRCINDEX  *si;
00182 register RAY  *r;
00183 {
00184        double  dist2, safedist2, dist2cent, rad2;
00185        FVECT  v;
00186        register SRCREC  *sp;
00187        int  pi;
00188                                    /* first check point location */
00189        clrpart(si->spt);
00190        sp = source + si->sn;
00191        rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
00192        v[0] = r->rorg[0] - sp->sloc[0];
00193        v[1] = r->rorg[1] - sp->sloc[1];
00194        v[2] = r->rorg[2] - sp->sloc[2];
00195        dist2 = DOT(v,sp->ss[SU]);
00196        safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
00197        dist2 *= dist2 / safedist2;
00198        dist2cent = DOT(v,v);
00199        dist2 = dist2cent - dist2;
00200        if (dist2 <= rad2) {        /* point inside extended cylinder */
00201               si->np = 0;
00202               return;
00203        }
00204        safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
00205        if (dist2 <= 4.*rad2 ||            /* point too close to subdivide */
00206                      dist2cent >= safedist2) {   /* or too far */
00207               setpart(si->spt, 0, S0);
00208               si->np = 1;
00209               return;
00210        }
00211        pi = 0;
00212        si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
00213                      sp->sloc, sp->ss[SU], safedist2);
00214 }
00215 
00216 
00217 static int
00218 cyl_partit(ro, pt, pi, mp, cent, axis, d2)       /* slice a cylinder */
00219 FVECT  ro;
00220 unsigned char  *pt;
00221 register int  *pi;
00222 int  mp;
00223 FVECT  cent, axis;
00224 double  d2;
00225 {
00226        FVECT  newct, newax;
00227        int  npl, npu;
00228 
00229        if (mp < 2 || dist2(ro, cent) >= d2) {    /* hit limit? */
00230               setpart(pt, *pi, S0);
00231               (*pi)++;
00232               return(1);
00233        }
00234                                    /* subdivide */
00235        setpart(pt, *pi, SU);
00236        (*pi)++;
00237        newax[0] = .5*axis[0];
00238        newax[1] = .5*axis[1];
00239        newax[2] = .5*axis[2];
00240        d2 *= 0.25;
00241                                    /* lower half */
00242        newct[0] = cent[0] - newax[0];
00243        newct[1] = cent[1] - newax[1];
00244        newct[2] = cent[2] - newax[2];
00245        npl = cyl_partit(ro, pt, pi, mp/2, newct, newax, d2);
00246                                    /* upper half */
00247        newct[0] = cent[0] + newax[0];
00248        newct[1] = cent[1] + newax[1];
00249        newct[2] = cent[2] + newax[2];
00250        npu = cyl_partit(ro, pt, pi, mp/2, newct, newax, d2);
00251                                    /* return total */
00252        return(npl + npu);
00253 }
00254 
00255 
00256 void
00257 flatpart(si, r)                           /* partition a flat source */
00258 register SRCINDEX  *si;
00259 register RAY  *r;
00260 {
00261        register RREAL  *vp;
00262        FVECT  v;
00263        double  du2, dv2;
00264        int  pi;
00265 
00266        clrpart(si->spt);
00267        vp = source[si->sn].sloc;
00268        v[0] = r->rorg[0] - vp[0];
00269        v[1] = r->rorg[1] - vp[1];
00270        v[2] = r->rorg[2] - vp[2];
00271        vp = source[si->sn].snorm;
00272        if (DOT(v,vp) <= 0.) {             /* behind source */
00273               si->np = 0;
00274               return;
00275        }
00276        dv2 = 2.*r->rweight/srcsizerat;
00277        dv2 *= dv2;
00278        vp = source[si->sn].ss[SU];
00279        du2 = dv2 * DOT(vp,vp);
00280        vp = source[si->sn].ss[SV];
00281        dv2 *= DOT(vp,vp);
00282        pi = 0;
00283        si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
00284               source[si->sn].sloc,
00285               source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
00286 }
00287 
00288 
00289 static int
00290 flt_partit(ro, pt, pi, mp, cent, u, v, du2, dv2) /* partition flatty */
00291 FVECT  ro;
00292 unsigned char  *pt;
00293 register int  *pi;
00294 int  mp;
00295 FVECT  cent, u, v;
00296 double  du2, dv2;
00297 {
00298        double  d2;
00299        FVECT  newct, newax;
00300        int  npl, npu;
00301 
00302        if (mp < 2 || ((d2 = dist2(ro, cent)) >= du2
00303                      && d2 >= dv2)) {     /* hit limit? */
00304               setpart(pt, *pi, S0);
00305               (*pi)++;
00306               return(1);
00307        }
00308        if (du2 > dv2) {                   /* subdivide in U */
00309               setpart(pt, *pi, SU);
00310               (*pi)++;
00311               newax[0] = .5*u[0];
00312               newax[1] = .5*u[1];
00313               newax[2] = .5*u[2];
00314               u = newax;
00315               du2 *= 0.25;
00316        } else {                           /* subdivide in V */
00317               setpart(pt, *pi, SV);
00318               (*pi)++;
00319               newax[0] = .5*v[0];
00320               newax[1] = .5*v[1];
00321               newax[2] = .5*v[2];
00322               v = newax;
00323               dv2 *= 0.25;
00324        }
00325                                    /* lower half */
00326        newct[0] = cent[0] - newax[0];
00327        newct[1] = cent[1] - newax[1];
00328        newct[2] = cent[2] - newax[2];
00329        npl = flt_partit(ro, pt, pi, mp/2, newct, u, v, du2, dv2);
00330                                    /* upper half */
00331        newct[0] = cent[0] + newax[0];
00332        newct[1] = cent[1] + newax[1];
00333        newct[2] = cent[2] + newax[2];
00334        npu = flt_partit(ro, pt, pi, mp/2, newct, u, v, du2, dv2);
00335                             /* return total */
00336        return(npl + npu);
00337 }
00338 
00339 
00340 double
00341 scylform(sn, dir)           /* compute cosine for cylinder's projection */
00342 int  sn;
00343 register FVECT  dir;        /* assume normalized */
00344 {
00345        register RREAL  *dv;
00346        double  d;
00347 
00348        dv = source[sn].ss[SU];
00349        d = DOT(dir, dv);
00350        d *= d / DOT(dv,dv);
00351        return(sqrt(1. - d));
00352 }