home *** CD-ROM | disk | FTP | other *** search
/ Encyclopedia of Graphics File Formats Companion / GFF_CD.ISO / formats / uray / code / intersec.c < prev    next >
C/C++ Source or Header  |  1994-06-20  |  9KB  |  294 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *            Copyright (c) 1988, David B. Wecker            *
  4.  *                All Rights Reserved                *
  5.  *                                    *
  6.  * This file is part of DBW_uRAY                    *
  7.  *                                    *
  8.  * DBW_uRAY is distributed in the hope that it will be useful, but    *
  9.  * WITHOUT ANY WARRANTY. No author or distributor accepts        *
  10.  * responsibility to anyone for the consequences of using it or for    *
  11.  * whether it serves any particular purpose or works at all, unless    *
  12.  * he says so in writing. Refer to the DBW_uRAY General Public        *
  13.  * License for full details.                        *
  14.  *                                    *
  15.  * Everyone is granted permission to copy, modify and redistribute    *
  16.  * DBW_uRAY, but only under the conditions described in the        *
  17.  * DBW_uRAY General Public License. A copy of this license is        *
  18.  * supposed to have been given to you along with DBW_uRAY so you    *
  19.  * can know your rights and responsibilities. It should be in a file    *
  20.  * named LICENSE. Among other things, the copyright notice and this    *
  21.  * notice must be preserved on all copies.                *
  22.  ************************************************************************
  23.  *                                    *
  24.  * Authors:                                *
  25.  *    DBW - David B. Wecker                        *
  26.  *                                    *
  27.  * Versions:                                *
  28.  *    V1.0 881023 DBW    - First released version            *
  29.  *    V1.1 881110 DBW - Fixed scan coherence code            *
  30.  *    V1.2 881125 DBW - Removed ALL scan coherence code (useless)    *
  31.  *              added "fat" extent boxes            *
  32.  *    V1.3 881203 DBW - Fixed single precision TOLerances        *
  33.  *                                    *
  34.  ************************************************************************/
  35.  
  36. #include "uray.h"
  37.  
  38. /*************************************************************************/
  39. /******** This module computes all ray object intersections **************/
  40. /*************************************************************************/
  41.  
  42. /* see if we hit a bounding box (n1 = near side, f1 = far side if good) */
  43. /* near point = from + n1 * direction                    */
  44. /* far  point = from + f1 * direction                    */
  45.  
  46. NODE    *hitextent(s,n,f)
  47. EXTENT    *s;
  48. FLTDBL    *n,*f;
  49.     {
  50.  
  51.     /* Use registers to be as efficient as possible */
  52.     register FLTDBL  n1 = -MYHUGE,f1 = MYHUGE,n2,f2,tval,sval;
  53.     register int     i;
  54.  
  55.     /* check each dimension of the bounding box */
  56.     for (i=0; i<3; i++) {
  57.  
  58.     /* see which side of the box that we're on */
  59.     if ((tval=dirinv[i]) < 0.0) {
  60.         n2  = (s->max[i] - (sval=from[i]))    * tval;
  61.         f2  = (s->min[i] -  sval)        * tval;
  62.         }
  63.     else {
  64.         n2  = (s->min[i] - (sval=from[i]))    * tval;
  65.         f2  = (s->max[i] -  sval)        * tval;
  66.         }
  67.  
  68.     /* if new near is worse than the old near (and old far) then bad */
  69.     if (n2 > n1 && (n1=n2) > f1) return NULL;
  70.  
  71.  
  72.     /* if new far is worse than old far and (and old near) or if we're
  73.         too close to 0 then return bad */
  74.     if (f2 < f1 && ((f1=f2) < n1 || f1 < TOL)) return NULL;
  75.     }
  76.  
  77.     /* are we worse than the far we ENTERED the routine with? */
  78.     if (n1 > *f) return NULL;
  79.  
  80.     /* set the return values */
  81.     *n = n1;
  82.     *f = f1;
  83.  
  84.     /* return the object hit */
  85.     return (NODE *)s;
  86.     }
  87.  
  88. /* see if we hit the sphere in question */
  89. NODE    *hitsphere(s,n1)
  90. SPHERE    *s;
  91. FLTDBL    *n1;
  92.     {
  93.     register FLTDBL    r_r,d_r,radical,t1,t2;
  94.     VEC            r;
  95.  
  96.     vcomb(-1.0, from, s->cen, r);
  97.     r_r = vdot(r,r);
  98.     d_r = vdot(direction,r);
  99.  
  100.     if ((radical = (d_r*d_r) + s->rad - r_r) < 0.0) return NULL;
  101.     radical = sqrt(radical);
  102.  
  103.     if (d_r < radical)  { t1 = d_r + radical; t2 = d_r - radical; }
  104.     else        { t1 = d_r - radical; t2 = d_r + radical; }
  105.  
  106.     /* since there are 2 possible intersection points, pick the best */
  107.     if (t1 <= TOL)    t1 = t2;
  108.  
  109.     /* are we too close to the intersection point */
  110.     if (t1 <= TOL)    return NULL;
  111.  
  112.     /* are we further than the bounding box intersection */
  113.     if (t1  > *n1)    return NULL;
  114.  
  115.     /* return the intersection point */
  116.     *n1 = t1;
  117.  
  118.     /* return the object that we hit */
  119.     return (NODE *)s;
  120.     }
  121.  
  122. /* take care of all planar intersections (quad, triangle, ring) */
  123. NODE    *hitplane(q,n1)
  124. RING        *q;
  125. FLTDBL        *n1;
  126.     {
  127.     VEC        sfs,h;
  128.     FLTDBL  det,r;
  129.  
  130.     /* First get the determinant of the interseection plane */
  131.     det  = q->v1[0] *      ((q->v2[1] * direction[2])-(q->v2[2] * direction[1]));
  132.     det -= q->v2[0] *      ((q->v1[1] * direction[2])-(q->v1[2] * direction[1]));
  133.     det += direction[0] * ((q->v1[1] * q->v2[2])    -(q->v1[2] * q->v2[1]));
  134.  
  135.     /* 0 determinant says that we missed the plane of intersection */
  136.     if (det == 0.0) return NULL;
  137.  
  138.     vcomb(-1.0,q->p0,from,h);
  139.  
  140.     /* now build the exact intersection point */
  141.     sfs[2]  = h[0] * ((q->v1[1] * q->v2[2]) - (q->v1[2] * q->v2[1]));
  142.     sfs[2] -= h[1] * ((q->v1[0] * q->v2[2]) - (q->v1[2] * q->v2[0]));
  143.     sfs[2] += h[2] * ((q->v1[0] * q->v2[1]) - (q->v1[1] * q->v2[0]));
  144.     sfs[2] /= det;
  145.     sfs[2]  = -sfs[2];
  146.  
  147.     /* Too close to 0? */
  148.     if (sfs[2] <= TOL)    return NULL;
  149.  
  150.     /* Further than bounding box? */
  151.     if (sfs[2] >  *n1)    return NULL;
  152.  
  153.     /* Now build comparison points for each object type */
  154.     sfs[0]  = h[0] * ((q->v2[1] * direction[2]) - (q->v2[2] * direction[1]));
  155.     sfs[0] -= h[1] * ((q->v2[0] * direction[2]) - (q->v2[2] * direction[0]));
  156.     sfs[0] += h[2] * ((q->v2[0] * direction[1]) - (q->v2[1] * direction[0]));
  157.     sfs[0] /= det;
  158.  
  159.     sfs[1]  = h[0] * ((q->v1[2] * direction[1]) - (q->v1[1] * direction[2]));
  160.     sfs[1] += h[1] * ((q->v1[0] * direction[2]) - (q->v1[2] * direction[0]));
  161.     sfs[1] -= h[2] * ((q->v1[0] * direction[1]) - (q->v1[1] * direction[0]));
  162.     sfs[1] /= det;
  163.  
  164.     switch (q->typ) {
  165.  
  166.     /* Make sure that we aren't off the triangles surface */
  167.     case TYP_T:
  168.     if (sfs[0] < 0.0 || sfs[1] < 0.0 || sfs[0]+sfs[1] > 1.0) return NULL;
  169.     break;
  170.  
  171.     /* make sure that we're within the bounds of the quad */
  172.     case TYP_Q:
  173.     if (sfs[0] < 0.0 || sfs[1] < 0.0 || sfs[0] > 1.0 || sfs[1] > 1.0)
  174.         return NULL;
  175.     break;
  176.  
  177.     /* make sure that we're between the ring's radii */
  178.     case TYP_R:
  179.     r = sfs[0] * sfs[0] + sfs[1] * sfs[1];
  180.     if (r > q->rad2 || r < q->rad1) return NULL;
  181.     break;
  182.     }
  183.  
  184.     /* return the intersection point */
  185.     *n1 = sfs[2];
  186.     return (NODE *)q;
  187.     }
  188.  
  189. /* recursive part of the intersection() routine */
  190. NODE *intersection2(cur,n1,f1)
  191. NODE    *cur;
  192. FLTDBL    *n1,*f1;
  193.     {
  194.  
  195.     register NODE    *best = NULL,*hit;
  196.     register int    i;
  197.     FLTDBL        ln1,lf1,ln2,lf2;
  198.  
  199.     /* save copies of the near and far intersection points (for extents) */
  200.     ln1 = *n1;
  201.     lf1 = *f1;
  202.  
  203.     /* check each object in the hierarchy */
  204.     while (cur) {
  205.  
  206.     /* re-init the intersection points */
  207.     hit = NULL;
  208.     ln2 = ln1;
  209.     lf2 = lf1;
  210.  
  211.     /* do the specific object */
  212.     switch (cur->typ) {
  213.  
  214.         /* if we hit an extent, call ourselves recursively */
  215.         case TYP_E:    
  216.         if (hitextent(cur,&ln2,&lf2)) {
  217.         gehits++;
  218.         hit = intersection2(((EXTENT *)cur)->child,&ln2,&lf2);
  219.         }
  220.         else fehits++;
  221.         break;
  222.  
  223.         /* if we hit a sphere, save the intersection info */
  224.         case TYP_S:
  225.         hit = hitsphere(cur,&lf2);
  226.         ln2 = lf2;
  227.         break;
  228.  
  229.         /* if we hit a planar, save the intersection info */
  230.         case TYP_T:
  231.         case TYP_Q:
  232.         case TYP_R:
  233.         hit = hitplane(cur,&lf2);
  234.         ln2 = lf2;
  235.         break;
  236.         }
  237.  
  238.     /* if NOT an extent, update the node statistics */
  239.     if (cur->typ != TYP_E) {
  240.         if (hit)    gnhits++;
  241.         else    fnhits++;
  242.         }
  243.  
  244.     /* if hit a non extent, then save this is the best so far */
  245.     if (hit && hit->typ != TYP_E) {
  246.         best = hit;
  247.         ln1  = ln2;
  248.         lf1  = lf2;
  249.         }
  250.  
  251.     /* go get next object */
  252.     cur  = cur->next;
  253.     }
  254.  
  255.     /* if we got something... return in */
  256.     if (best) {
  257.     *n1 = ln1;
  258.     *f1 = lf1;
  259.     }
  260.     return best;
  261.     }
  262.  
  263. /* return closest node (frm + n1 * dir) or NULL */
  264. NODE *intersection(frm, dir, n1)
  265. VEC    frm, dir;
  266. FLTDBL    *n1;
  267.     {
  268.     register int    i;
  269.     NODE        *hit = NULL;
  270.     FLTDBL        f1;
  271.  
  272.     /* initialize near and far */
  273.     *n1    = -MYHUGE;
  274.     f1  =  MYHUGE;
  275.  
  276.     /* set up global vectors FROM, DIRECTION and DIRINV (inverse of dir) */
  277.     for (i=0; i<3; i++) {
  278.     from[i]        = frm[i];
  279.     direction[i]    = dir[i];
  280.     if (FABS(direction[i]) < TOL)  {
  281.         if (direction[i] < 0.0)   dirinv[i] = -MYHUGE;
  282.         else              dirinv[i] =  MYHUGE;
  283.             }
  284.     else                  dirinv[i] = 1.0 / direction[i];
  285.     }
  286.  
  287.     /* call the recursive part of this routine */
  288.     hit = intersection2(nodes,n1,&f1);
  289.  
  290.     /* return the result */
  291.     return hit;
  292.     }
  293.  
  294.