home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch5_2 / quad.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-22  |  4.7 KB  |  134 lines
  1. /* ------------------------------------------------------------------------- *\
  2.    QUAD.C :
  3.  
  4.    by Christophe Schlick and Gilles Subrenat (15 May 1994)
  5.  
  6.    "Ray Intersection of Tessellated Surfaces : Quadrangles versus Triangles"
  7.    in Graphics Gems V (edited by A. Paeth), Academic Press
  8. \* ------------------------------------------------------------------------- */
  9.  
  10. #include "quad.h"
  11.  
  12. /*
  13. ** Macro definitions
  14. */
  15. #define MY_TOL                ((real) 0.0001)
  16.  
  17. #define LARGEST_COMPONENT(A)  (ABS((A).x) > ABS((A).y) ? \
  18.                               (ABS((A).x) > ABS((A).z) ? 'x' : 'z') : \
  19.                               (ABS((A).y) > ABS((A).z) ? 'y' : 'z'))
  20.  
  21. /*
  22. ** Check if the point is in the quadrangle
  23. */
  24. static bool point_in_quad (QUAD *Quad, HIT *Hit)
  25. {
  26.     char       LargestComponent;             /* of the normal vector         */
  27.     realvec2   A, B, C, D;                   /* Projected vertices           */
  28.     realvec2   M;                            /* Projected intersection point */
  29.     realvec2   AB, BC, CD, AD, AM, AE;       /* Miscellanous 3D-vectors      */
  30.     real       u, v;                         /* Parametric coordinates       */
  31.     real       a, b, c, SqrtDelta;           /* Quadratic equation           */
  32.     bool       Intersection = FALSE;         /* Intersection flag            */
  33.     realvec2   Vector;                       /* Temporary 2D-vector          */
  34.  
  35.     /*
  36.     ** Projection on the plane that is most parallel to the facet
  37.     */
  38.     LargestComponent = LARGEST_COMPONENT(Quad->Normal);
  39.  
  40.     if (LargestComponent == 'x') {
  41.         A.x = Quad->A.y; B.x = Quad->B.y; C.x = Quad->C.y; D.x = Quad->D.y;
  42.         M.x = Hit->Point.y;
  43.     }
  44.     else {
  45.         A.x = Quad->A.x; B.x = Quad->B.x; C.x = Quad->C.x; D.x = Quad->D.x;
  46.         M.x = Hit->Point.x;
  47.     }
  48.  
  49.     if (LargestComponent == 'z') {
  50.         A.y = Quad->A.y; B.y = Quad->B.y; C.y = Quad->C.y; D.y = Quad->D.y;
  51.         M.y = Hit->Point.y;
  52.     }
  53.     else {
  54.         A.y = Quad->A.z; B.y = Quad->B.z; C.y = Quad->C.z; D.y = Quad->D.z;
  55.         M.y = Hit->Point.z;
  56.     }
  57.  
  58.     SUB_VEC2 (AB, B, A); SUB_VEC2 (BC, C, B);
  59.     SUB_VEC2 (CD, D, C); SUB_VEC2 (AD, D, A);
  60.     ADD_VEC2 (AE, CD, AB); NEG_VEC2 (AE, AE); SUB_VEC2 (AM, M, A);
  61.  
  62.     if (ZERO_TOL (DELTA_VEC2(AB, CD), MY_TOL))             /* case AB // CD */
  63.     {
  64.         SUB_VEC2 (Vector, AB, CD);
  65.         v = DELTA_VEC2(AM, Vector) / DELTA_VEC2(AD, Vector);
  66.         if ((v >= 0.0) && (v <= 1.0)) {
  67.             b = DELTA_VEC2(AB, AD) - DELTA_VEC2(AM, AE);
  68.             c = DELTA_VEC2 (AM, AD);
  69.             u = ZERO_TOL(b, MY_TOL) ? -1.0 : c/b;
  70.             Intersection = ((u >= 0.0) && (u <= 1.0));
  71.         }
  72.     }
  73.     else if (ZERO_TOL(DELTA_VEC2(BC, AD), MY_TOL))         /* case AD // BC */
  74.     {
  75.         ADD_VEC2 (Vector, AD, BC);
  76.         u = DELTA_VEC2(AM, Vector) / DELTA_VEC2(AB, Vector);
  77.         if ((u >= 0.0) && (u <= 1.0)) {
  78.             b = DELTA_VEC2(AD, AB) - DELTA_VEC2(AM, AE);
  79.             c = DELTA_VEC2 (AM, AB);
  80.             v = ZERO_TOL(b, MY_TOL) ? -1.0 : c/b;
  81.             Intersection = ((v >= 0.0) && (v <= 1.0));
  82.         }
  83.     }
  84.     else                                                    /* general case */
  85.     {
  86.         a = DELTA_VEC2(AB, AE); c = - DELTA_VEC2 (AM,AD);
  87.         b = DELTA_VEC2(AB, AD) - DELTA_VEC2(AM, AE);
  88.         a = -0.5/a; b *= a; c *= (a + a); SqrtDelta = b*b + c;
  89.         if (SqrtDelta >= 0.0) {
  90.             SqrtDelta = sqrt(SqrtDelta);
  91.             u = b - SqrtDelta;
  92.             if ((u < 0.0) || (u > 1.0))        /* we want u between 0 and 1 */
  93.                 u = b + SqrtDelta;
  94.             if ((u >= 0.0) && (u <= 1.0)) {
  95.                 v = AD.x + u * AE.x;
  96.                 if (ZERO_TOL(v, MY_TOL))
  97.                     v = (AM.y - u * AB.y) / (AD.y + u * AE.y);
  98.                 else
  99.                     v = (AM.x - u * AB.x) / v;
  100.                 Intersection = ((v >= 0.0) && (v <= 1.0));
  101.             }
  102.         }
  103.     }
  104.  
  105.     if (Intersection) {
  106.         Hit->u = u;
  107.         Hit->v = v;
  108.     }
  109.     return (Intersection);
  110. }
  111.  
  112. /*
  113. ** Search for an intersection between a quadrangle and a ray
  114. */
  115. bool hit_ray_quad (RAY *Ray, QUAD *Quad, HIT *Hit)
  116. {
  117.     realvec3     Point;
  118.  
  119.     /* if the ray is parallel to the quadrangle, there is no intersection */
  120.     Hit->Distance = DOT_VEC3 (Ray->Vector, Quad->Normal);
  121.     if (ZERO_TOL(Hit->Distance, MY_TOL)) return (FALSE);
  122.  
  123.     /* compute ray intersection with the plane of the quadrangle */
  124.     SUB_VEC3 (Point, Quad->A, Ray->Point);
  125.     Hit->Distance = DOT_VEC3 (Point, Quad->Normal) / Hit->Distance;
  126.     MULS_VEC3 (Hit->Point, Ray->Vector, Hit->Distance);
  127.     INC_VEC3 (Hit->Point, Ray->Point);
  128.  
  129.     /* is the point in the quadrangle ? */
  130.     return (point_in_quad(Quad, Hit));
  131. }
  132.  
  133. /* ------------------------------------------------------------------------- */
  134.