home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / gems / gemsiii / intqdr.c < prev    next >
C/C++ Source or Header  |  1992-03-16  |  3KB  |  104 lines
  1. #include    "GraphicsGems.h"
  2.  
  3. /* ----    intqdr - Intersect a ray with a quadric surface. --------------    */
  4. /*                                    */
  5. /*                                    */
  6. /*    Description:                            */
  7. /*        Intqdr determines the intersection of a ray with a quadric    */
  8. /*        surface in matrix form.                    */
  9. /*                                    */
  10. /*    On entry:                            */
  11. /*        raybase = The coordinate defining the base of the        */
  12. /*              intersecting ray.                    */
  13. /*        raycos  = The direction cosines of the above ray.        */
  14. /*        base    = The location of the base of the quadric.        */
  15. /*        axis    = The axis of symmetry for the quadric.        */
  16. /*        radius  = The bounding radius from the axis of symmetry.    */
  17. /*        Q       = A 4x4 coefficient matrix representing the quadric    */
  18. /*              surface.                        */
  19. /*                                    */
  20. /*    On return:                            */
  21. /*        rin     = The entering distance of the intersection.    */
  22. /*        rout    = The leaving  distance of the intersection.    */
  23. /*                                    */
  24. /*                                    */
  25. /*    Returns:  True if the ray intersects the quadric surface.    */
  26. /*                                    */
  27. /* --------------------------------------------------------------------    */
  28.  
  29.  
  30. boolean    intqdr    (raybase,raycos,base,axis,radius,Q,rin,rout)
  31.  
  32.     Point3    raybase;        /* Base of the intersection ray    */
  33.     Vector3    raycos;            /* Direction cosines of the ray    */
  34.     Vector3    base;            /* Base point for the quadric    */
  35.     Vector3    axis;            /* Axis of symmetry        */
  36.     double    radius;            /* Bounding radius        */
  37.     Matrix4    Q;            /* Quadric coefficient matrix    */
  38.     double    *rin;            /* Entering distance        */
  39.     double    *rout;            /* Leaving  distance        */
  40.  
  41. {
  42.     boolean    hit;            /* True if ray intersects quad.    */
  43.     double    qa, qb, qc, qd, qe,    /* Coefficient matrix terms    */
  44.         qf, qg, qh, qi, qj;
  45.     double    k2, k1, k0;        /* Quadric coefficients        */
  46.     double    disc;
  47.     Vector3    D;            /* Ray base to quadric base    */
  48.     double    d;            /* Shortest distance between    */
  49.                     /* ray and axis of symmetry    */
  50.     Vector3    N;
  51.  
  52.  
  53. /*    Compute the perpendicular distance between the symmetry axis    */
  54. /*    and the ray.                            */
  55.  
  56.     V3Sub   (&base,&raybase,&D);    /* compute shortest distance    */
  57.     V3Cross (&raycos,&axis,&N);
  58.     if  ((d=V3Length(&N)) == 0) {            /* if parrallel    */
  59.         V3Cross (&raycos,&D,&N);
  60.         d = fabs (magntd(N));
  61.     } else {
  62.         V3Normalize (&N);
  63.         d = fabs (V3Dot(&D,&N)/d);
  64.     }
  65.  
  66.     hit = (d <= radius);
  67.     if  (!hit) return (hit);    /* If outside of cylinder    */
  68.  
  69.  
  70.     qa = Q.element[0][0];
  71.     qb = Q.element[0][1] + Q.element[1][0];
  72.     qc = Q.element[0][2] + Q.element[2][0];
  73.     qd = Q.element[0][3] + Q.element[3][0];
  74.     qe = Q.element[1][1];
  75.     qf = Q.element[1][2] + Q.element[2][1];
  76.     qg = Q.element[1][3] + Q.element[3][1];
  77.     qh = Q.element[2][2];
  78.     qi = Q.element[2][3] + Q.element[3][2];
  79.     qj = Q.element[3][3];
  80.  
  81.     k2 = raycos.x * (raycos.x*qa + raycos.y*qb) +
  82.          raycos.y * (raycos.y*qe + raycos.z*qf) +
  83.          raycos.z * (raycos.z*qh + raycos.z*qc);
  84.     k1 = raycos.x * ((raybase.x*2.*qa + raybase.y*qb + raybase.z*qc) + qd) +
  85.          raycos.y * ((raybase.x*qb + raybase.y*2.*qe + raybase.z*qf) + qg) +
  86.          raycos.z * ((raybase.x*qc + raybase.y*qg + raybase.z*2.*qh) + qi);
  87.     k0 = raybase.x * (raybase.x*qa + raybase.y*qb + raybase.z*qc + qd) +
  88.          raybase.y * (               raybase.y*qe + raybase.z*qf + qg) +
  89.          raybase.z * (                              raybase.z*qh + qi) +
  90.          qj;
  91.  
  92.     disc = k1*k1 - 4.*k2*k0;
  93.   
  94.     hit  = (disc >= 0.0);
  95.  
  96.     if  (hit) {                 /* If ray hits quadric    */
  97.         disc  =  sqrt(disc);
  98.         *rin  = (-k1 - disc) / (2.*k2);    /*    entering distance    */
  99.         *rout = (-k1 + disc) / (2.*k2);    /*    leaving distance    */
  100.     }
  101.  
  102.     return (hit);
  103. }
  104.