home *** CD-ROM | disk | FTP | other *** search
/ Altsys Virtuoso 2.0K / virtuoso_20k.iso / DemoApps / Graphics / Viewers / raytracers / rpi / Source / intersect.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-10-08  |  3.7 KB  |  184 lines
  1.  
  2. /*
  3.  * (c) 1988 by George Kyriazis
  4.  */
  5.  
  6. /*
  7.  * Intersection routines
  8.  */
  9.  
  10. #include    "vector.h"
  11. #include    "ray.h"
  12. #include    <math.h>
  13.  
  14. /* 
  15.  * Intersect ray with sphere
  16.  */
  17. struct    intersect    sphere(obj, r)
  18. struct    obj    *obj;
  19. struct    ray    r;
  20. {
  21.     struct    vector    v;
  22.     struct    vector    n;    /* normal vector */
  23.     struct    intersect    i;
  24.     double    b, c, d;
  25.     double    sol1, sol2;
  26.     struct    vector    center;
  27.  
  28.     i.obj = NULL;
  29.  
  30. /* find out what the center is at this time */
  31.     center = svproduct( Time, obj->time );
  32.     center = vadd( obj->data.sphere.center, center );
  33.  
  34.     v = vsub( r.pos, center );
  35.     b = 2 * vdot( r.dir, v );
  36.     c = vdot(v, v) - obj->data.sphere.radius;
  37.  
  38.     d = b * b - 4 * c;
  39.     if( d < 0 )
  40.         return i;
  41.     d = sqrt(d);
  42.     sol1 = ( -b + d ) / 2;
  43.     sol2 = ( -b - d ) / 2;
  44.     if( sol1 <= 0 )
  45.         sol1 = sol2;
  46.     if( sol2 <= 0 )
  47.         sol2 = sol1;
  48.     i.t = (sol1 < sol2) ? sol1 : sol2 ;
  49. /* if intersection is behind eye */
  50.     if(i.t <= 0)
  51.         return i;
  52.     i.obj = obj;    
  53.  
  54. /* calculate the normal.  It is just the direction of the radius */
  55.     n = vsub(vadd(r.pos, svproduct(i.t, r.dir)), center);
  56.     i.n = norm(n);
  57.  
  58.     return i;
  59. }
  60.  
  61. /*
  62.  * intersect ray with a quadrangle
  63.  */
  64. struct    intersect    quad(obj, r)
  65. struct    obj    *obj;
  66. struct    ray    r;
  67. {
  68.     struct    intersect    i;
  69.     double    x0, y0, z0;
  70.     double    dx1, dy1, dz1, dx2, dy2, dz2, qx, qy, qz;
  71.     double    alpha, a, b, t;
  72.     double    d, dalpha, da, db;
  73.     struct    vector    n;
  74.     double    pdx, pdy, pdz;
  75.     double    size;
  76.  
  77.     i.obj = NULL;
  78.  
  79.     x0 = obj->data.quad.p1.x;
  80.     y0 = obj->data.quad.p1.y;
  81.     z0 = obj->data.quad.p1.z;
  82.     dx1 = x0 - obj->data.quad.p2.x;
  83.     dy1 = y0 - obj->data.quad.p2.y;
  84.     dz1 = z0 - obj->data.quad.p2.z;
  85.     dx2 = x0 - obj->data.quad.p3.x;
  86.     dy2 = y0 - obj->data.quad.p3.y;
  87.     dz2 = z0 - obj->data.quad.p3.z;
  88.     qx = r.dir.x;
  89.     qy = r.dir.y;
  90.     qz = r.dir.z;
  91.  
  92.     d = qx * ( ( dy1 * dz2 ) - ( dy2 * dz1 ) )
  93.         - dx1 * ( ( qy * dz2 ) - ( qz * dy2 ) )
  94.         + dx2 * ( ( qy * dz1 ) - ( qz * dy1 ) ) ;
  95. /* if no intersection */
  96.     if( ABS(d) < MINT )
  97.         return i;
  98.  
  99. /* use the right time */
  100.     pdx = x0 - r.pos.x + Time * obj->time.x;
  101.     pdy = y0 - r.pos.y + Time * obj->time.y;
  102.     pdz = z0 - r.pos.z + Time * obj->time.z;
  103.  
  104.     dalpha = pdx * ( ( dy1 * dz2 ) - ( dz1 * dy2 ) )
  105.         - dx1 * ( ( pdy * dz2 ) - ( pdz * dy2 ) )
  106.         + dx2 * ( ( pdy * dz1 ) - ( pdz * dy1 ) ) ;
  107.     alpha = dalpha / d;
  108.  
  109. /* if intersection behind the eye */
  110.     if( alpha <= 0 )
  111.         return i;
  112.  
  113.     da = qx * ( ( pdy * dz2 ) - ( pdz * dy2 ) )
  114.         - pdx * ( ( qy * dz2 ) - ( qz * dy2 ) )
  115.         + dx2 * ( ( qy * pdz ) - ( qz * pdy ) ) ;
  116.  
  117.     db = qx * ( ( dy1 * pdz ) - ( dz1 * pdy ) )
  118.         - dx1 * ( ( qy * pdz ) - ( qz * pdy ) )
  119.         + pdx * ( ( qy * dz1 ) - ( qz * dy1 ) );
  120.  
  121.     a = da / d;
  122.     b = db / d;
  123.  
  124. /* check if intersection within quad */
  125.     if( ( a < 0 ) || ( a > 1 ) )
  126.         return i;
  127.     if( ( b < 0 ) || ( b > 1 ) )
  128.         return i;
  129.  
  130. /* assign the object */
  131.     i.obj = obj;
  132.     i.t = alpha;
  133. /* normal to the plane */
  134.     n.x = dy1 * dz2 - dy2 * dz1;
  135.     n.y = - ( dx1 * dz2 - dx2 * dz1 );
  136.     n.z = dx1 * dy2 - dx2 * dy1;
  137.     i.n = norm(n);
  138.  
  139. /* reverse the normal if desired */
  140.     if( vdot( i.n, r.dir ) > 0 )
  141.         i.n = vneg( i.n );
  142.  
  143.     return i;
  144. }
  145.  
  146.  
  147.  
  148. struct    intersect    intersect(r)
  149. struct    ray    r;
  150. {
  151.     int    i;
  152.     struct    intersect    inter, intermin;
  153.  
  154.     intermin.obj = NULL;
  155.  
  156.     for(i = 0; i < noo; i++) {
  157.         objtestline++;
  158. /* choose the appropriate routine for each object */
  159.         switch( obj[i].type ) {
  160.             case SPHERE:
  161.                 inter = sphere( &obj[i], r);
  162.                 break;
  163.             case SQUARE:
  164.                 inter = quad( &obj[i], r);
  165.                 break;
  166.             default:
  167.                 inter.obj = NULL;
  168.                 break;
  169.         }
  170. /* update the minimum intersection distance if the new intersection  */
  171. /* exists (ray intersects the object), and the intersection distance */
  172. /* is smaller that the one logged and also the object intersected is */
  173. /* not the object that the ray is originating from.                  */
  174.         if( inter.obj &&
  175.             ( !intermin.obj || 
  176.                 (inter.t < intermin.t) ) &&
  177.             ( inter.obj != r.obj ) )
  178.             intermin = inter;
  179.     }
  180.  
  181.     return intermin;
  182. }
  183.  
  184.