home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / gems / g_gemsii.lha / GemsII / intersect / inttor.c < prev   
Encoding:
C/C++ Source or Header  |  1991-09-22  |  3.7 KB  |  121 lines
  1.  
  2. #include    <math.h>
  3. #include    "GraphicsGems.h"
  4.  
  5. /* ----    inttor - Intersect a ray with a torus. ------------------------    */
  6. /*                                    */
  7. /*                                    */
  8. /*    Description:                            */
  9. /*        Inttor determines the intersection of a ray with a torus.    */
  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. /*        center  = The center location of the torus.            */
  16. /*        radius  = The major radius of the torus.            */
  17. /*        rplane  = The minor radius in the plane of the torus.    */
  18. /*        rnorm   = The minor radius normal to the plane of the torus.*/
  19. /*        tran    = A 4x4 transformation matrix that will position    */
  20. /*              the torus at the origin and orient it such that    */
  21. /*              the plane of the torus lyes in the x-z plane.    */
  22. /*                                    */
  23. /*    On return:                            */
  24. /*        nhits   = The number of intersections the ray makes with    */
  25. /*              the torus.                    */
  26. /*        rhits   = The entering/leaving distances of the        */
  27. /*              intersections.                    */
  28. /*                                    */
  29. /*    Returns:  True if the ray intersects the torus.            */
  30. /*                                    */
  31. /* --------------------------------------------------------------------    */
  32.  
  33.  
  34. int    inttor    (raybase,raycos,center,radius,rplane,rnorm,tran,nhits,rhits)
  35.  
  36.     Point3    raybase;        /* Base of the intersection ray    */
  37.     Vector3    raycos;            /* Direction cosines of the ray    */
  38.     Point3    center;            /* Center of the torus        */
  39.     double    radius;            /* Major radius of the torus    */
  40.     double    rplane;            /* Minor planer radius        */
  41.     double    rnorm;            /* Minor normal radius        */
  42.     Matrix4    tran;            /* Transformation matrix    */
  43.     int *    nhits;            /* Number of intersections    */
  44.     double    rhits[4];        /* Intersection distances    */
  45.  
  46. {
  47.     int    hit;            /* True if ray intersects torus    */
  48.     double    rsphere;        /* Bounding sphere radius    */
  49.     Vector3    Base, Dcos;        /* Transformed intersection ray    */
  50.     double    rmin, rmax;        /* Root bounds            */
  51.     double    yin, yout;
  52.     double    rho, a0, b0;        /* Related constants        */
  53.     double    f, l, t, g, q, m, u;    /* Ray dependent terms        */
  54.     double    C[5];            /* Quartic coefficents        */
  55.  
  56. extern    int    intsph ();        /* Intersect ray with sphere    */
  57. extern    int    SolveQuartic ();    /* Solve quartic equation    */
  58.  
  59.  
  60.     *nhits  = 0;
  61.  
  62. /*    Compute the intersection of the ray with a bounding sphere.    */
  63.  
  64.     rsphere = radius + MAX (rplane,rnorm);
  65.     hit     = intsph (raybase,raycos,center,rsphere,&rmin,&rmax);
  66.  
  67.     if  (!hit) return (hit);    /* If ray misses bounding sphere*/
  68.  
  69. /*    Transform the intersection ray                    */
  70.  
  71.     Base = raybase;
  72.     Dcos = raycos;
  73.     V3MulPointByMatrix  (&Base,&tran);
  74.     V3MulVectorByMatrix (&Dcos,&tran);
  75.  
  76. /*    Bound the torus by two parallel planes rnorm from the x-z plane.*/
  77.  
  78.     yin  = Base.y + rmin * Dcos.y;
  79.     yout = Base.y + rmax * Dcos.y;
  80.     hit  = !( (yin >  rnorm && yout >  rnorm) ||
  81.           (yin < -rnorm && yout < -rnorm) );
  82.  
  83.     if  (!hit) return (hit);    /* If ray is above/below torus.    */
  84.  
  85. /*    Compute constants related to the torus.                */
  86.  
  87.     rho = rplane*rplane / (rnorm*rnorm);
  88.     a0  = 4. * radius*radius;
  89.     b0  = radius*radius - rplane*rplane;
  90.  
  91. /*    Compute ray dependent terms.                    */
  92.  
  93.     f = 1. - Dcos.y*Dcos.y;
  94.     l = 2. * (Base.x*Dcos.x + Base.z*Dcos.z);
  95.     t = Base.x*Base.x + Base.z*Base.z;
  96.     g = f + rho * Dcos.y*Dcos.y;
  97.     q = a0 / (g*g);
  98.     m = (l + 2.*rho*Dcos.y*Base.y) / g;
  99.     u = (t +    rho*Base.y*Base.y + b0) / g;
  100.  
  101. /*    Compute the coefficients of the quartic.            */
  102.  
  103.     C[4] = 1.0;
  104.     C[3] = 2. * m;
  105.     C[2] = m*m + 2.*u - q*f;
  106.     C[1] = 2.*m*u - q*l;
  107.     C[0] = u*u - q*t;
  108.     
  109. /*    Use quartic root solver found in "Graphics Gems" by Jochen    */
  110. /*    Schwarze.                            */
  111.  
  112.     *nhits = SolveQuartic (C,rhits);
  113.  
  114. /*    SolveQuartic returns root pairs in reversed order.        */
  115.     m = rhits[0]; u = rhits[1]; rhits[0] = u; rhits[1] = m;
  116.     m = rhits[2]; u = rhits[3]; rhits[2] = u; rhits[3] = m;
  117.  
  118.     return (*nhits != 0);
  119. }
  120.  
  121.