home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 231.lha / URay_v1.0 / intersect.c < prev    next >
C/C++ Source or Header  |  1989-04-04  |  8KB  |  293 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 COPYING. 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.  *                                    *
  33.  ************************************************************************/
  34.  
  35. #include "uray.h"
  36.  
  37. /*************************************************************************/
  38. /******** This module computes all ray object intersections **************/
  39. /*************************************************************************/
  40.  
  41. /* see if we hit a bounding box (n1 = near side, f1 = far side if good) */
  42. /* near point = from + n1 * direction                    */
  43. /* far  point = from + f1 * direction                    */
  44.  
  45. NODE    *hitextent(s,n,f)
  46. EXTENT    *s;
  47. FLTDBL    *n,*f;
  48.     {
  49.  
  50.     /* Use registers to be as efficient as possible */
  51.     register FLTDBL  n1 = -MYHUGE,f1 = MYHUGE,n2,f2,tval,sval;
  52.     register int     i;
  53.  
  54.     /* check each dimension of the bounding box */
  55.     for (i=0; i<3; i++) {
  56.  
  57.     /* see which side of the box that we're on */
  58.     if ((tval=dirinv[i]) < 0.0) {
  59.         n2  = (s->max[i] - (sval=from[i]))    * tval;
  60.         f2  = (s->min[i] -  sval)        * tval;
  61.         }
  62.     else {
  63.         n2  = (s->min[i] - (sval=from[i]))    * tval;
  64.         f2  = (s->max[i] -  sval)        * tval;
  65.         }
  66.  
  67.     /* if new near is worse than the old near (and old far) then bad */
  68.     if (n2 > n1 && (n1=n2) > f1) return NULL;
  69.  
  70.  
  71.     /* if new far is worse than old far and (and old near) or if we're
  72.         too close to 0 then return bad */
  73.     if (f2 < f1 && ((f1=f2) < n1 || f1 < TOL)) return NULL;
  74.     }
  75.  
  76.     /* are we worse than the far we ENTERED the routine with? */
  77.     if (n1 > *f) return NULL;
  78.  
  79.     /* set the return values */
  80.     *n = n1;
  81.     *f = f1;
  82.  
  83.     /* return the object hit */
  84.     return (NODE *)s;
  85.     }
  86.  
  87. /* see if we hit the sphere in question */
  88. NODE    *hitsphere(s,n1)
  89. SPHERE    *s;
  90. FLTDBL    *n1;
  91.     {
  92.     register FLTDBL    r_r,d_r,radical,t1,t2;
  93.     VEC            r;
  94.  
  95.     vcomb(-1.0, from, s->cen, r);
  96.     r_r = vdot(r,r);
  97.     d_r = vdot(direction,r);
  98.  
  99.     if ((radical = (d_r*d_r) + s->rad - r_r) < 0.0) return NULL;
  100.     radical = sqrt(radical);
  101.  
  102.     if (d_r < radical)  { t1 = d_r + radical; t2 = d_r - radical; }
  103.     else        { t1 = d_r - radical; t2 = d_r + radical; }
  104.  
  105.     /* since there are 2 possible intersection points, pick the best */
  106.     if (t1 <= TOL)    t1 = t2;
  107.  
  108.     /* are we too close to the intersection point */
  109.     if (t1 <= TOL)    return NULL;
  110.  
  111.     /* are we further than the bounding box intersection */
  112.     if (t1  > *n1)    return NULL;
  113.  
  114.     /* return the intersection point */
  115.     *n1 = t1;
  116.  
  117.     /* return the object that we hit */
  118.     return (NODE *)s;
  119.     }
  120.  
  121. /* take care of all planar intersections (quad, triangle, ring) */
  122. NODE    *hitplane(q,n1)
  123. RING        *q;
  124. FLTDBL        *n1;
  125.     {
  126.     VEC        sfs,h;
  127.     FLTDBL  det,r;
  128.  
  129.     /* First get the determinant of the interseection plane */
  130.     det  = q->v1[0] *      ((q->v2[1] * direction[2])-(q->v2[2] * direction[1]));
  131.     det -= q->v2[0] *      ((q->v1[1] * direction[2])-(q->v1[2] * direction[1]));
  132.     det += direction[0] * ((q->v1[1] * q->v2[2])    -(q->v1[2] * q->v2[1]));
  133.  
  134.     /* 0 determinant says that we missed the plane of intersection */
  135.     if (det == 0.0) return NULL;
  136.  
  137.     vcomb(-1.0,q->p0,from,h);
  138.  
  139.     /* now build the exact intersection point */
  140.     sfs[2]  = h[0] * ((q->v1[1] * q->v2[2]) - (q->v1[2] * q->v2[1]));
  141.     sfs[2] -= h[1] * ((q->v1[0] * q->v2[2]) - (q->v1[2] * q->v2[0]));
  142.     sfs[2] += h[2] * ((q->v1[0] * q->v2[1]) - (q->v1[1] * q->v2[0]));
  143.     sfs[2] /= det;
  144.     sfs[2]  = -sfs[2];
  145.  
  146.     /* Too close to 0? */
  147.     if (sfs[2] <= TOL)    return NULL;
  148.  
  149.     /* Further than bounding box? */
  150.     if (sfs[2] >  *n1)    return NULL;
  151.  
  152.     /* Now build comparison points for each object type */
  153.     sfs[0]  = h[0] * ((q->v2[1] * direction[2]) - (q->v2[2] * direction[1]));
  154.     sfs[0] -= h[1] * ((q->v2[0] * direction[2]) - (q->v2[2] * direction[0]));
  155.     sfs[0] += h[2] * ((q->v2[0] * direction[1]) - (q->v2[1] * direction[0]));
  156.     sfs[0] /= det;
  157.  
  158.     sfs[1]  = h[0] * ((q->v1[2] * direction[1]) - (q->v1[1] * direction[2]));
  159.     sfs[1] += h[1] * ((q->v1[0] * direction[2]) - (q->v1[2] * direction[0]));
  160.     sfs[1] -= h[2] * ((q->v1[0] * direction[1]) - (q->v1[1] * direction[0]));
  161.     sfs[1] /= det;
  162.  
  163.     switch (q->typ) {
  164.  
  165.     /* Make sure that we aren't off the triangles surface */
  166.     case TYP_T:
  167.     if (sfs[0] < 0.0 || sfs[1] < 0.0 || sfs[0]+sfs[1] > 1.0) return NULL;
  168.     break;
  169.  
  170.     /* make sure that we're within the bounds of the quad */
  171.     case TYP_Q:
  172.     if (sfs[0] < 0.0 || sfs[1] < 0.0 || sfs[0] > 1.0 || sfs[1] > 1.0)
  173.         return NULL;
  174.     break;
  175.  
  176.     /* make sure that we're between the ring's radii */
  177.     case TYP_R:
  178.     r = sfs[0] * sfs[0] + sfs[1] * sfs[1];
  179.     if (r > q->rad2 || r < q->rad1) return NULL;
  180.     break;
  181.     }
  182.  
  183.     /* return the intersection point */
  184.     *n1 = sfs[2];
  185.     return (NODE *)q;
  186.     }
  187.  
  188. /* recursive part of the intersection() routine */
  189. NODE *intersection2(cur,n1,f1)
  190. NODE    *cur;
  191. FLTDBL    *n1,*f1;
  192.     {
  193.  
  194.     register NODE    *best = NULL,*hit;
  195.     register int    i;
  196.     FLTDBL        ln1,lf1,ln2,lf2;
  197.  
  198.     /* save copies of the near and far intersection points (for extents) */
  199.     ln1 = *n1;
  200.     lf1 = *f1;
  201.  
  202.     /* check each object in the hierarchy */
  203.     while (cur) {
  204.  
  205.     /* re-init the intersection points */
  206.     hit = NULL;
  207.     ln2 = ln1;
  208.     lf2 = lf1;
  209.  
  210.     /* do the specific object */
  211.     switch (cur->typ) {
  212.  
  213.         /* if we hit an extent, call ourselves recursively */
  214.         case TYP_E:    
  215.         if (hitextent(cur,&ln2,&lf2)) {
  216.         gehits++;
  217.         hit = intersection2(((EXTENT *)cur)->child,&ln2,&lf2);
  218.         }
  219.         else fehits++;
  220.         break;
  221.  
  222.         /* if we hit a sphere, save the intersection info */
  223.         case TYP_S:
  224.         hit = hitsphere(cur,&lf2);
  225.         ln2 = lf2;
  226.         break;
  227.  
  228.         /* if we hit a planar, save the intersection info */
  229.         case TYP_T:
  230.         case TYP_Q:
  231.         case TYP_R:
  232.         hit = hitplane(cur,&lf2);
  233.         ln2 = lf2;
  234.         break;
  235.         }
  236.  
  237.     /* if NOT an extent, update the node statistics */
  238.     if (cur->typ != TYP_E) {
  239.         if (hit)    gnhits++;
  240.         else    fnhits++;
  241.         }
  242.  
  243.     /* if hit a non extent, then save this is the best so far */
  244.     if (hit && hit->typ != TYP_E) {
  245.         best = hit;
  246.         ln1  = ln2;
  247.         lf1  = lf2;
  248.         }
  249.  
  250.     /* go get next object */
  251.     cur  = cur->next;
  252.     }
  253.  
  254.     /* if we got something... return in */
  255.     if (best) {
  256.     *n1 = ln1;
  257.     *f1 = lf1;
  258.     }
  259.     return best;
  260.     }
  261.  
  262. /* return closest node (frm + n1 * dir) or NULL */
  263. NODE *intersection(frm, dir, n1)
  264. VEC    frm, dir;
  265. FLTDBL    *n1;
  266.     {
  267.     register int    i;
  268.     NODE        *hit = NULL;
  269.     FLTDBL        f1;
  270.  
  271.     /* initialize near and far */
  272.     *n1    = -MYHUGE;
  273.     f1  =  MYHUGE;
  274.  
  275.     /* set up global vectors FROM, DIRECTION and DIRINV (inverse of dir) */
  276.     for (i=0; i<3; i++) {
  277.     from[i]        = frm[i];
  278.     direction[i]    = dir[i];
  279.     if (FABS(direction[i]) < TOL)  {
  280.         if (direction[i] < 0.0)   dirinv[i] = -MYHUGE;
  281.         else              dirinv[i] =  MYHUGE;
  282.             }
  283.     else                  dirinv[i] = 1.0 / direction[i];
  284.     }
  285.  
  286.     /* call the recursive part of this routine */
  287.     hit = intersection2(nodes,n1,&f1);
  288.  
  289.     /* return the result */
  290.     return hit;
  291.     }
  292.  
  293.