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