home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 8 / CDASC08.ISO / NEWS / RADIANCE / SRC / RT / SRCSUPP.C < prev    next >
C/C++ Source or Header  |  1993-10-07  |  12KB  |  499 lines
  1. /* Copyright (c) 1992 Regents of the University of California */
  2.  
  3. #ifndef lint
  4. static char SCCSid[] = "@(#)srcsupp.c 2.6 8/28/92 LBL";
  5. #endif
  6.  
  7. /*
  8.  *  Support routines for source objects and materials
  9.  */
  10.  
  11. #include  "ray.h"
  12.  
  13. #include  "otypes.h"
  14.  
  15. #include  "source.h"
  16.  
  17. #include  "cone.h"
  18.  
  19. #include  "face.h"
  20.  
  21. #define SRCINC        4        /* realloc increment for array */
  22.  
  23. SRCREC  *source = NULL;            /* our list of sources */
  24. int  nsources = 0;            /* the number of sources */
  25.  
  26. SRCFUNC  sfun[NUMOTYPE];        /* source dispatch table */
  27.  
  28.  
  29. initstypes()            /* initialize source dispatch table */
  30. {
  31.     extern VSMATERIAL  mirror_vs, direct1_vs, direct2_vs;
  32.     extern int  fsetsrc(), ssetsrc(), sphsetsrc(), cylsetsrc(), rsetsrc();
  33.     extern int  nopart(), flatpart(), cylpart();
  34.     extern double  fgetplaneq(), rgetplaneq();
  35.     extern double  fgetmaxdisk(), rgetmaxdisk();
  36.     static SOBJECT  fsobj = {fsetsrc, flatpart, fgetplaneq, fgetmaxdisk};
  37.     static SOBJECT  ssobj = {ssetsrc, nopart};
  38.     static SOBJECT  sphsobj = {sphsetsrc, nopart};
  39.     static SOBJECT  cylsobj = {cylsetsrc, cylpart};
  40.     static SOBJECT  rsobj = {rsetsrc, flatpart, rgetplaneq, rgetmaxdisk};
  41.  
  42.     sfun[MAT_MIRROR].mf = &mirror_vs;
  43.     sfun[MAT_DIRECT1].mf = &direct1_vs;
  44.     sfun[MAT_DIRECT2].mf = &direct2_vs;
  45.     sfun[OBJ_FACE].of = &fsobj;
  46.     sfun[OBJ_SOURCE].of = &ssobj;
  47.     sfun[OBJ_SPHERE].of = &sphsobj;
  48.     sfun[OBJ_CYLINDER].of = &cylsobj;
  49.     sfun[OBJ_RING].of = &rsobj;
  50. }
  51.  
  52.  
  53. int
  54. newsource()            /* allocate new source in our array */
  55. {
  56.     if (nsources == 0)
  57.         source = (SRCREC *)malloc(SRCINC*sizeof(SRCREC));
  58.     else if (nsources%SRCINC == 0)
  59.         source = (SRCREC *)realloc((char *)source,
  60.                 (unsigned)(nsources+SRCINC)*sizeof(SRCREC));
  61.     if (source == NULL)
  62.         return(-1);
  63.     source[nsources].sflags = 0;
  64.     source[nsources].nhits = 1;
  65.     source[nsources].ntests = 2;    /* initial hit probability = 1/2 */
  66.     return(nsources++);
  67. }
  68.  
  69.  
  70. setflatss(src)                /* set sampling for a flat source */
  71. register SRCREC  *src;
  72. {
  73.     double  mult;
  74.     register int  i;
  75.  
  76.     src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
  77.     for (i = 0; i < 3; i++)
  78.         if (src->snorm[i] < 0.6 && src->snorm[i] > -0.6)
  79.             break;
  80.     src->ss[SV][i] = 1.0;
  81.     fcross(src->ss[SU], src->ss[SV], src->snorm);
  82.     mult = .5 * sqrt( src->ss2 / DOT(src->ss[SU],src->ss[SU]) );
  83.     for (i = 0; i < 3; i++)
  84.         src->ss[SU][i] *= mult;
  85.     fcross(src->ss[SV], src->snorm, src->ss[SU]);
  86. }
  87.  
  88.  
  89. fsetsrc(src, so)            /* set a face as a source */
  90. register SRCREC  *src;
  91. OBJREC  *so;
  92. {
  93.     register FACE  *f;
  94.     register int  i, j;
  95.     double  d;
  96.     
  97.     src->sa.success = 2*AIMREQT-1;        /* bitch on second failure */
  98.     src->so = so;
  99.                         /* get the face */
  100.     f = getface(so);
  101.                         /* find the center */
  102.     for (j = 0; j < 3; j++) {
  103.         src->sloc[j] = 0.0;
  104.         for (i = 0; i < f->nv; i++)
  105.             src->sloc[j] += VERTEX(f,i)[j];
  106.         src->sloc[j] /= (double)f->nv;
  107.     }
  108.     if (!inface(src->sloc, f))
  109.         objerror(so, USER, "cannot hit center");
  110.     src->sflags |= SFLAT;
  111.     VCOPY(src->snorm, f->norm);
  112.     src->ss2 = f->area;
  113.                         /* find maximum radius */
  114.     src->srad = 0.;
  115.     for (i = 0; i < f->nv; i++) {
  116.         d = dist2(VERTEX(f,i), src->sloc);
  117.         if (d > src->srad)
  118.             src->srad = d;
  119.     }
  120.     src->srad = sqrt(src->srad);
  121.                         /* compute size vectors */
  122.     if (f->nv == 4 || (f->nv == 5 &&    /* parallelogram case */
  123.             dist2(VERTEX(f,0),VERTEX(f,4)) <= FTINY*FTINY))
  124.         for (j = 0; j < 3; j++) {
  125.             src->ss[SU][j] = .5*(VERTEX(f,1)[j]-VERTEX(f,0)[j]);
  126.             src->ss[SV][j] = .5*(VERTEX(f,3)[j]-VERTEX(f,0)[j]);
  127.         }
  128.     else
  129.         setflatss(src);
  130. }
  131.  
  132.  
  133. ssetsrc(src, so)            /* set a source as a source */
  134. register SRCREC  *src;
  135. register OBJREC  *so;
  136. {
  137.     double  theta;
  138.     
  139.     src->sa.success = 2*AIMREQT-1;        /* bitch on second failure */
  140.     src->so = so;
  141.     if (so->oargs.nfargs != 4)
  142.         objerror(so, USER, "bad arguments");
  143.     src->sflags |= SDISTANT;
  144.     VCOPY(src->sloc, so->oargs.farg);
  145.     if (normalize(src->sloc) == 0.0)
  146.         objerror(so, USER, "zero direction");
  147.     theta = PI/180.0/2.0 * so->oargs.farg[3];
  148.     if (theta <= FTINY)
  149.         objerror(so, USER, "zero size");
  150.     src->ss2 = 2.0*PI * (1.0 - cos(theta));
  151.                     /* the following is approximate */
  152.     src->srad = sqrt(src->ss2/PI);
  153.     VCOPY(src->snorm, src->sloc);
  154.     setflatss(src);            /* hey, whatever works */
  155.     src->ss[SW][0] = src->ss[SW][1] = src->ss[SW][2] = 0.0;
  156. }
  157.  
  158.  
  159. sphsetsrc(src, so)            /* set a sphere as a source */
  160. register SRCREC  *src;
  161. register OBJREC  *so;
  162. {
  163.     register int  i;
  164.  
  165.     src->sa.success = 2*AIMREQT-1;        /* bitch on second failure */
  166.     src->so = so;
  167.     if (so->oargs.nfargs != 4)
  168.         objerror(so, USER, "bad # arguments");
  169.     if (so->oargs.farg[3] <= FTINY)
  170.         objerror(so, USER, "illegal radius");
  171.     VCOPY(src->sloc, so->oargs.farg);
  172.     src->srad = so->oargs.farg[3];
  173.     src->ss2 = PI * src->srad * src->srad;
  174.     for (i = 0; i < 3; i++)
  175.         src->ss[SU][i] = src->ss[SV][i] = src->ss[SW][i] = 0.0;
  176.     for (i = 0; i < 3; i++)
  177.         src->ss[i][i] = .7236 * so->oargs.farg[3];
  178. }
  179.  
  180.  
  181. rsetsrc(src, so)            /* set a ring (disk) as a source */
  182. register SRCREC  *src;
  183. OBJREC  *so;
  184. {
  185.     register CONE  *co;
  186.     
  187.     src->sa.success = 2*AIMREQT-1;        /* bitch on second failure */
  188.     src->so = so;
  189.                         /* get the ring */
  190.     co = getcone(so, 0);
  191.     VCOPY(src->sloc, CO_P0(co));
  192.     if (CO_R0(co) > 0.0)
  193.         objerror(so, USER, "cannot hit center");
  194.     src->sflags |= SFLAT;
  195.     VCOPY(src->snorm, co->ad);
  196.     src->srad = CO_R1(co);
  197.     src->ss2 = PI * src->srad * src->srad;
  198.     setflatss(src);
  199. }
  200.  
  201.  
  202. cylsetsrc(src, so)            /* set a cylinder as a source */
  203. register SRCREC  *src;
  204. OBJREC  *so;
  205. {
  206.     register CONE  *co;
  207.     register int  i;
  208.     
  209.     src->sa.success = 4*AIMREQT-1;        /* bitch on fourth failure */
  210.     src->so = so;
  211.                         /* get the cylinder */
  212.     co = getcone(so, 0);
  213.     if (CO_R0(co) > .2*co->al)        /* heuristic constraint */
  214.         objerror(so, WARNING, "source aspect too small");
  215.     src->sflags |= SCYL;
  216.     for (i = 0; i < 3; i++)
  217.         src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
  218.     src->srad = .5*co->al;
  219.     src->ss2 = 2.*CO_R0(co)*co->al;
  220.                         /* set sampling vectors */
  221.     for (i = 0; i < 3; i++)
  222.         src->ss[SU][i] = .5 * co->al * co->ad[i];
  223.     src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
  224.     for (i = 0; i < 3; i++)
  225.         if (co->ad[i] < 0.6 && co->ad[i] > -0.6)
  226.             break;
  227.     src->ss[SV][i] = 1.0;
  228.     fcross(src->ss[SW], src->ss[SV], co->ad);
  229.     normalize(src->ss[SW]);
  230.     for (i = 0; i < 3; i++)
  231.         src->ss[SW][i] *= .8559 * CO_R0(co);
  232.     fcross(src->ss[SV], src->ss[SW], co->ad);
  233. }
  234.  
  235.  
  236. SPOT *
  237. makespot(m)            /* make a spotlight */
  238. register OBJREC  *m;
  239. {
  240.     register SPOT  *ns;
  241.  
  242.     if ((ns = (SPOT *)m->os) != NULL)
  243.         return(ns);
  244.     if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
  245.         return(NULL);
  246.     ns->siz = 2.0*PI * (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
  247.     VCOPY(ns->aim, m->oargs.farg+4);
  248.     if ((ns->flen = normalize(ns->aim)) == 0.0)
  249.         objerror(m, USER, "zero focus vector");
  250.     m->os = (char *)ns;
  251.     return(ns);
  252. }
  253.  
  254.  
  255. spotout(r, s, dist)        /* check if we're outside spot region */
  256. register RAY  *r;
  257. register SPOT  *s;
  258. int  dist;
  259. {
  260.     double  d;
  261.     FVECT  vd;
  262.     
  263.     if (s == NULL)
  264.         return(0);
  265.     if (dist) {            /* distant source */
  266.         vd[0] = s->aim[0] - r->rorg[0];
  267.         vd[1] = s->aim[1] - r->rorg[1];
  268.         vd[2] = s->aim[2] - r->rorg[2];
  269.         d = DOT(r->rdir,vd);
  270.         /*            wrong side?
  271.         if (d <= FTINY)
  272.             return(1);    */
  273.         d = DOT(vd,vd) - d*d;
  274.         if (PI*d > s->siz)
  275.             return(1);    /* out */
  276.         return(0);    /* OK */
  277.     }
  278.                     /* local source */
  279.     if (s->siz < 2.0*PI * (1.0 + DOT(s->aim,r->rdir)))
  280.         return(1);    /* out */
  281.     return(0);    /* OK */
  282. }
  283.  
  284.  
  285. double
  286. fgetmaxdisk(ocent, op)        /* get center and squared radius of face */
  287. FVECT  ocent;
  288. OBJREC  *op;
  289. {
  290.     double  maxrad2;
  291.     double  d;
  292.     register int  i, j;
  293.     register FACE  *f;
  294.     
  295.     f = getface(op);
  296.     if (f->area == 0.)
  297.         return(0.);
  298.     for (i = 0; i < 3; i++) {
  299.         ocent[i] = 0.;
  300.         for (j = 0; j < f->nv; j++)
  301.             ocent[i] += VERTEX(f,j)[i];
  302.         ocent[i] /= (double)f->nv;
  303.     }
  304.     d = DOT(ocent,f->norm);
  305.     for (i = 0; i < 3; i++)
  306.         ocent[i] += (f->offset - d)*f->norm[i];
  307.     maxrad2 = 0.;
  308.     for (j = 0; j < f->nv; j++) {
  309.         d = dist2(VERTEX(f,j), ocent);
  310.         if (d > maxrad2)
  311.             maxrad2 = d;
  312.     }
  313.     return(maxrad2);
  314. }
  315.  
  316.  
  317. double
  318. rgetmaxdisk(ocent, op)        /* get center and squared radius of ring */
  319. FVECT  ocent;
  320. OBJREC  *op;
  321. {
  322.     register CONE  *co;
  323.     
  324.     co = getcone(op, 0);
  325.     VCOPY(ocent, CO_P0(co));
  326.     return(CO_R1(co)*CO_R1(co));
  327. }
  328.  
  329.  
  330. double
  331. fgetplaneq(nvec, op)            /* get plane equation for face */
  332. FVECT  nvec;
  333. OBJREC  *op;
  334. {
  335.     register FACE  *fo;
  336.  
  337.     fo = getface(op);
  338.     VCOPY(nvec, fo->norm);
  339.     return(fo->offset);
  340. }
  341.  
  342.  
  343. double
  344. rgetplaneq(nvec, op)            /* get plane equation for ring */
  345. FVECT  nvec;
  346. OBJREC  *op;
  347. {
  348.     register CONE  *co;
  349.  
  350.     co = getcone(op, 0);
  351.     VCOPY(nvec, co->ad);
  352.     return(DOT(nvec, CO_P0(co)));
  353. }
  354.  
  355.  
  356. commonspot(sp1, sp2, org)    /* set sp1 to intersection of sp1 and sp2 */
  357. register SPOT  *sp1, *sp2;
  358. FVECT  org;
  359. {
  360.     FVECT  cent;
  361.     double  rad2, cos1, cos2;
  362.  
  363.     cos1 = 1. - sp1->siz/(2.*PI);
  364.     cos2 = 1. - sp2->siz/(2.*PI);
  365.     if (sp2->siz >= 2.*PI-FTINY)        /* BIG, just check overlap */
  366.         return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
  367.                     sqrt((1.-cos1*cos1)*(1.-cos2*cos2)));
  368.                 /* compute and check disks */
  369.     rad2 = intercircle(cent, sp1->aim, sp2->aim,
  370.             1./(cos1*cos1) - 1.,  1./(cos2*cos2) - 1.);
  371.     if (rad2 <= FTINY || normalize(cent) == 0.)
  372.         return(0);
  373.     VCOPY(sp1->aim, cent);
  374.     sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2));
  375.     return(1);
  376. }
  377.  
  378.  
  379. commonbeam(sp1, sp2, dir)    /* set sp1 to intersection of sp1 and sp2 */
  380. register SPOT  *sp1, *sp2;
  381. FVECT  dir;
  382. {
  383.     FVECT  cent, c1, c2;
  384.     double  rad2, d;
  385.     register int  i;
  386.                     /* move centers to common plane */
  387.     d = DOT(sp1->aim, dir);
  388.     for (i = 0; i < 3; i++)
  389.         c1[i] = sp1->aim[i] - d*dir[i];
  390.     d = DOT(sp2->aim, dir);
  391.     for (i = 0; i < 3; i++)
  392.         c2[i] = sp2->aim[i] - d*dir[i];
  393.                     /* compute overlap */
  394.     rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
  395.     if (rad2 <= FTINY)
  396.         return(0);
  397.     VCOPY(sp1->aim, cent);
  398.     sp1->siz = PI*rad2;
  399.     return(1);
  400. }
  401.  
  402.  
  403. checkspot(sp, nrm)        /* check spotlight for behind source */
  404. register SPOT  *sp;    /* spotlight */
  405. FVECT  nrm;        /* source surface normal */
  406. {
  407.     double  d, d1;
  408.  
  409.     d = DOT(sp->aim, nrm);
  410.     if (d > FTINY)            /* center in front? */
  411.         return(1);
  412.                     /* else check horizon */
  413.     d1 = 1. - sp->siz/(2.*PI);
  414.     return(1.-FTINY-d*d < d1*d1);
  415. }
  416.  
  417.  
  418. double
  419. spotdisk(oc, op, sp, pos)    /* intersect spot with object op */
  420. FVECT  oc;
  421. OBJREC  *op;
  422. register SPOT  *sp;
  423. FVECT  pos;
  424. {
  425.     FVECT  onorm;
  426.     double  offs, d, dist;
  427.     register int  i;
  428.  
  429.     offs = getplaneq(onorm, op);
  430.     d = -DOT(onorm, sp->aim);
  431.     if (d >= -FTINY && d <= FTINY)
  432.         return(0.);
  433.     dist = (DOT(pos, onorm) - offs)/d;
  434.     if (dist < 0.)
  435.         return(0.);
  436.     for (i = 0; i < 3; i++)
  437.         oc[i] = pos[i] + dist*sp->aim[i];
  438.     return(sp->siz*dist*dist/PI/(d*d));
  439. }
  440.  
  441.  
  442. double
  443. beamdisk(oc, op, sp, dir)    /* intersect beam with object op */
  444. FVECT  oc;
  445. OBJREC  *op;
  446. register SPOT  *sp;
  447. FVECT  dir;
  448. {
  449.     FVECT  onorm;
  450.     double  offs, d, dist;
  451.     register int  i;
  452.  
  453.     offs = getplaneq(onorm, op);
  454.     d = -DOT(onorm, dir);
  455.     if (d >= -FTINY && d <= FTINY)
  456.         return(0.);
  457.     dist = (DOT(sp->aim, onorm) - offs)/d;
  458.     for (i = 0; i < 3; i++)
  459.         oc[i] = sp->aim[i] + dist*dir[i];
  460.     return(sp->siz/PI/(d*d));
  461. }
  462.  
  463.  
  464. double
  465. intercircle(cc, c1, c2, r1s, r2s)    /* intersect two circles */
  466. FVECT  cc;            /* midpoint (return value) */
  467. FVECT  c1, c2;            /* circle centers */
  468. double  r1s, r2s;        /* radii squared */
  469. {
  470.     double  a2, d2, l;
  471.     FVECT  disp;
  472.     register int  i;
  473.  
  474.     for (i = 0; i < 3; i++)
  475.         disp[i] = c2[i] - c1[i];
  476.     d2 = DOT(disp,disp);
  477.                     /* circle within overlap? */
  478.     if (r1s < r2s) {
  479.         if (r2s >= r1s + d2) {
  480.             VCOPY(cc, c1);
  481.             return(r1s);
  482.         }
  483.     } else {
  484.         if (r1s >= r2s + d2) {
  485.             VCOPY(cc, c2);
  486.             return(r2s);
  487.         }
  488.     }
  489.     a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
  490.                     /* no overlap? */
  491.     if (a2 <= 0.)
  492.         return(0.);
  493.                     /* overlap, compute center */
  494.     l = sqrt((r1s - a2)/d2);
  495.     for (i = 0; i < 3; i++)
  496.         cc[i] = c1[i] + l*disp[i];
  497.     return(a2);
  498. }
  499.