home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 8 / CDASC08.ISO / NEWS / RADIANCE / SRC / RT / SRCSAMP.C < prev    next >
C/C++ Source or Header  |  1993-10-07  |  8KB  |  315 lines
  1. /* Copyright (c) 1991 Regents of the University of California */
  2.  
  3. #ifndef lint
  4. static char SCCSid[] = "@(#)srcsamp.c 2.4 9/22/92 LBL";
  5. #endif
  6.  
  7. /*
  8.  * Source sampling routines
  9.  */
  10.  
  11. #include  "ray.h"
  12.  
  13. #include  "source.h"
  14.  
  15. #include  "random.h"
  16.  
  17.  
  18. double
  19. nextssamp(r, si)        /* compute sample for source, rtn. distance */
  20. register RAY  *r;        /* origin is read, direction is set */
  21. register SRCINDEX  *si;        /* source index (modified to current) */
  22. {
  23.     int  cent[3], size[3], parr[2];
  24.     FVECT  vpos;
  25.     double  d;
  26.     register int  i;
  27. nextsample:
  28.     while (++si->sp >= si->np) {    /* get next sample */
  29.         if (++si->sn >= nsources)
  30.             return(0.0);    /* no more */
  31.         if (source[si->sn].sflags & SSKIP)
  32.             si->np = 0;
  33.         else if (srcsizerat <= FTINY)
  34.             nopart(si, r);
  35.         else {
  36.             for (i = si->sn; source[i].sflags & SVIRTUAL;
  37.                     i = source[i].sa.sv.sn)
  38.                 ;        /* partition source */
  39.             (*sfun[source[i].so->otype].of->partit)(si, r);
  40.         }
  41.         si->sp = -1;
  42.     }
  43.                     /* get partition */
  44.     cent[0] = cent[1] = cent[2] = 0;
  45.     size[0] = size[1] = size[2] = MAXSPART;
  46.     parr[0] = 0; parr[1] = si->sp;
  47.     if (!skipparts(cent, size, parr, si->spt))
  48.         error(CONSISTENCY, "bad source partition in nextssamp");
  49.                     /* compute sample */
  50.     if (dstrsrc > FTINY) {            /* jitter sample */
  51.         dimlist[ndims] = si->sn + 8831;
  52.         dimlist[ndims+1] = si->sp + 3109;
  53.         d = urand(ilhash(dimlist,ndims+2)+samplendx);
  54.         if (source[si->sn].sflags & SFLAT) {
  55.             multisamp(vpos, 2, d);
  56.             vpos[2] = 0.5;
  57.         } else
  58.             multisamp(vpos, 3, d);
  59.         for (i = 0; i < 3; i++)
  60.             vpos[i] = dstrsrc * (1. - 2.*vpos[i]) *
  61.                     (double)size[i]/MAXSPART;
  62.     } else
  63.         vpos[0] = vpos[1] = vpos[2] = 0.0;
  64.  
  65.     for (i = 0; i < 3; i++)
  66.         vpos[i] += (double)cent[i]/MAXSPART;
  67.                     /* compute direction */
  68.     for (i = 0; i < 3; i++)
  69.         r->rdir[i] = source[si->sn].sloc[i] +
  70.                 vpos[SU]*source[si->sn].ss[SU][i] +
  71.                 vpos[SV]*source[si->sn].ss[SV][i] +
  72.                 vpos[SW]*source[si->sn].ss[SW][i];
  73.  
  74.     if (!(source[si->sn].sflags & SDISTANT))
  75.         for (i = 0; i < 3; i++)
  76.             r->rdir[i] -= r->rorg[i];
  77.                     /* compute distance */
  78.     if ((d = normalize(r->rdir)) == 0.0)
  79.         goto nextsample;        /* at source! */
  80.  
  81.                     /* compute sample size */
  82.     si->dom  = source[si->sn].ss2;
  83.     if (source[si->sn].sflags & SFLAT) {
  84.         si->dom *= sflatform(si->sn, r->rdir);
  85.         si->dom *= size[SU]*size[SV]/(MAXSPART*(double)MAXSPART);
  86.     } else if (source[si->sn].sflags & SCYL) {
  87.         si->dom *= scylform(si->sn, r->rdir);
  88.         si->dom *= size[SU]/(double)MAXSPART;
  89.     } else {
  90.         si->dom *= size[SU]*size[SV]*(double)size[SW] /
  91.                 (MAXSPART*MAXSPART*(double)MAXSPART) ;
  92.     }
  93.     if (source[si->sn].sflags & SDISTANT)
  94.         return(FHUGE);
  95.     if (si->dom <= 1e-4)
  96.         goto nextsample;        /* behind source? */
  97.     si->dom /= d*d;
  98.     return(d);        /* sample OK, return distance */
  99. }
  100.  
  101.  
  102. skipparts(ct, sz, pp, pt)        /* skip to requested partition */
  103. int  ct[3], sz[3];        /* center and size of partition (returned) */
  104. register int  pp[2];        /* current index, number to skip (modified) */
  105. unsigned char  *pt;        /* partition array */
  106. {
  107.     register int  p;
  108.                     /* check this partition */
  109.     p = spart(pt, pp[0]);
  110.     pp[0]++;
  111.     if (p == S0)            /* leaf partition */
  112.         if (pp[1]) {
  113.             pp[1]--;
  114.             return(0);    /* not there yet */
  115.         } else
  116.             return(1);    /* we've arrived */
  117.                 /* else check lower */
  118.     sz[p] >>= 1;
  119.     ct[p] -= sz[p];
  120.     if (skipparts(ct, sz, pp, pt))
  121.         return(1);    /* return hit */
  122.                 /* else check upper */
  123.     ct[p] += sz[p] << 1;
  124.     if (skipparts(ct, sz, pp, pt))
  125.         return(1);    /* return hit */
  126.                 /* else return to starting position */
  127.     ct[p] -= sz[p];
  128.     sz[p] <<= 1;
  129.     return(0);        /* return miss */
  130. }
  131.  
  132.  
  133. nopart(si, r)            /* single source partition */
  134. register SRCINDEX  *si;
  135. RAY  *r;
  136. {
  137.     clrpart(si->spt);
  138.     setpart(si->spt, 0, S0);
  139.     si->np = 1;
  140. }
  141.  
  142.  
  143. cylpart(si, r)            /* partition a cylinder */
  144. SRCINDEX  *si;
  145. register RAY  *r;
  146. {
  147.     double  dist2, safedist2, dist2cent, rad2;
  148.     FVECT  v;
  149.     register SRCREC  *sp;
  150.     int  pi;
  151.                     /* first check point location */
  152.     clrpart(si->spt);
  153.     sp = source + si->sn;
  154.     rad2 = 1.365 * DOT(sp->ss[SV],sp->ss[SV]);
  155.     v[0] = r->rorg[0] - sp->sloc[0];
  156.     v[1] = r->rorg[1] - sp->sloc[1];
  157.     v[2] = r->rorg[2] - sp->sloc[2];
  158.     dist2 = DOT(v,sp->ss[SU]);
  159.     safedist2 = DOT(sp->ss[SU],sp->ss[SU]);
  160.     dist2 *= dist2 / safedist2;
  161.     dist2cent = DOT(v,v);
  162.     dist2 = dist2cent - dist2;
  163.     if (dist2 <= rad2) {        /* point inside extended cylinder */
  164.         si->np = 0;
  165.         return;
  166.     }
  167.     safedist2 *= 4.*r->rweight*r->rweight/(srcsizerat*srcsizerat);
  168.     if (dist2 <= 4.*rad2 ||        /* point too close to subdivide */
  169.             dist2cent >= safedist2) {    /* or too far */
  170.         setpart(si->spt, 0, S0);
  171.         si->np = 1;
  172.         return;
  173.     }
  174.     pi = 0;
  175.     si->np = cyl_partit(r->rorg, si->spt, &pi, MAXSPART,
  176.             sp->sloc, sp->ss[SU], safedist2);
  177. }
  178.  
  179.  
  180. static int
  181. cyl_partit(ro, pt, pi, mp, cent, axis, d2)    /* slice a cylinder */
  182. FVECT  ro;
  183. unsigned char  *pt;
  184. register int  *pi;
  185. int  mp;
  186. FVECT  cent, axis;
  187. double  d2;
  188. {
  189.     FVECT  newct, newax;
  190.     int  npl, npu;
  191.  
  192.     if (mp < 2 || dist2(ro, cent) >= d2) {    /* hit limit? */
  193.         setpart(pt, *pi, S0);
  194.         (*pi)++;
  195.         return(1);
  196.     }
  197.                     /* subdivide */
  198.     setpart(pt, *pi, SU);
  199.     (*pi)++;
  200.     newax[0] = .5*axis[0];
  201.     newax[1] = .5*axis[1];
  202.     newax[2] = .5*axis[2];
  203.     d2 *= 0.25;
  204.                     /* lower half */
  205.     newct[0] = cent[0] - newax[0];
  206.     newct[1] = cent[1] - newax[1];
  207.     newct[2] = cent[2] - newax[2];
  208.     npl = cyl_partit(ro, pt, pi, mp/2, newct, newax, d2);
  209.                     /* upper half */
  210.     newct[0] = cent[0] + newax[0];
  211.     newct[1] = cent[1] + newax[1];
  212.     newct[2] = cent[2] + newax[2];
  213.     npu = cyl_partit(ro, pt, pi, mp/2, newct, newax, d2);
  214.                     /* return total */
  215.     return(npl + npu);
  216. }
  217.  
  218.  
  219. flatpart(si, r)                /* partition a flat source */
  220. register SRCINDEX  *si;
  221. register RAY  *r;
  222. {
  223.     register FLOAT  *vp;
  224.     FVECT  v;
  225.     double  du2, dv2;
  226.     int  pi;
  227.  
  228.     clrpart(si->spt);
  229.     vp = source[si->sn].sloc;
  230.     v[0] = r->rorg[0] - vp[0];
  231.     v[1] = r->rorg[1] - vp[1];
  232.     v[2] = r->rorg[2] - vp[2];
  233.     vp = source[si->sn].snorm;
  234.     if (DOT(v,vp) <= FTINY) {    /* behind source */
  235.         si->np = 0;
  236.         return;
  237.     }
  238.     dv2 = 2.*r->rweight/srcsizerat;
  239.     dv2 *= dv2;
  240.     vp = source[si->sn].ss[SU];
  241.     du2 = dv2 * DOT(vp,vp);
  242.     vp = source[si->sn].ss[SV];
  243.     dv2 *= DOT(vp,vp);
  244.     pi = 0;
  245.     si->np = flt_partit(r->rorg, si->spt, &pi, MAXSPART,
  246.         source[si->sn].sloc,
  247.         source[si->sn].ss[SU], source[si->sn].ss[SV], du2, dv2);
  248. }
  249.  
  250.  
  251. static int
  252. flt_partit(ro, pt, pi, mp, cent, u, v, du2, dv2)    /* partition flatty */
  253. FVECT  ro;
  254. unsigned char  *pt;
  255. register int  *pi;
  256. int  mp;
  257. FVECT  cent, u, v;
  258. double  du2, dv2;
  259. {
  260.     double  d2;
  261.     FVECT  newct, newax;
  262.     int  npl, npu;
  263.  
  264.     if (mp < 2 || ((d2 = dist2(ro, cent)) >= du2
  265.             && d2 >= dv2)) {    /* hit limit? */
  266.         setpart(pt, *pi, S0);
  267.         (*pi)++;
  268.         return(1);
  269.     }
  270.     if (du2 > dv2) {            /* subdivide in U */
  271.         setpart(pt, *pi, SU);
  272.         (*pi)++;
  273.         newax[0] = .5*u[0];
  274.         newax[1] = .5*u[1];
  275.         newax[2] = .5*u[2];
  276.         u = newax;
  277.         du2 *= 0.25;
  278.     } else {                /* subdivide in V */
  279.         setpart(pt, *pi, SV);
  280.         (*pi)++;
  281.         newax[0] = .5*v[0];
  282.         newax[1] = .5*v[1];
  283.         newax[2] = .5*v[2];
  284.         v = newax;
  285.         dv2 *= 0.25;
  286.     }
  287.                     /* lower half */
  288.     newct[0] = cent[0] - newax[0];
  289.     newct[1] = cent[1] - newax[1];
  290.     newct[2] = cent[2] - newax[2];
  291.     npl = flt_partit(ro, pt, pi, mp/2, newct, u, v, du2, dv2);
  292.                     /* upper half */
  293.     newct[0] = cent[0] + newax[0];
  294.     newct[1] = cent[1] + newax[1];
  295.     newct[2] = cent[2] + newax[2];
  296.     npu = flt_partit(ro, pt, pi, mp/2, newct, u, v, du2, dv2);
  297.                 /* return total */
  298.     return(npl + npu);
  299. }
  300.  
  301.  
  302. double
  303. scylform(sn, dir)        /* compute cosine for cylinder's projection */
  304. int  sn;
  305. register FVECT  dir;        /* assume normalized */
  306. {
  307.     register FLOAT  *dv;
  308.     double  d;
  309.  
  310.     dv = source[sn].ss[SU];
  311.     d = DOT(dir, dv);
  312.     d *= d / DOT(dv,dv);
  313.     return(sqrt(1. - d));
  314. }
  315.