home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / msdos / raytrace / rayshade / src / torus.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-07-18  |  9.1 KB  |  371 lines
  1. /*
  2.  * torus.c
  3.  *
  4.  * Copyright (C) 1990, 1991, Mark Podlipec, Craig E. Kolb
  5.  * All rights reserved.
  6.  *
  7.  * This software may be freely copied, modified, and redistributed
  8.  * provided that this copyright notice is preserved on all copies.
  9.  *
  10.  * You may not distribute this software, in whole or in part, as part of
  11.  * any commercial product without the express consent of the authors.
  12.  *
  13.  * There is no warranty or other guarantee of fitness of this software
  14.  * for any purpose.  It is provided solely "as is".
  15.  *
  16.  * $Id: torus.c,v 4.0 91/07/17 14:39:28 kolb Exp Locker: kolb $
  17.  *
  18.  * $Log:    torus.c,v $
  19.  * Revision 4.0  91/07/17  14:39:28  kolb
  20.  * Initial version.
  21.  * 
  22.  */
  23. #include "geom.h"
  24. #include "torus.h"
  25.  
  26. static Methods *iTorusMethods = NULL;
  27. static char torusName[] = "torus";
  28. unsigned long TorusTests, TorusHits;
  29.  
  30. /*
  31.  * Create & return reference to a torus.
  32.  */
  33. Torus *
  34. TorusCreate(a, b, pos, norm)
  35. Float a, b;
  36. Vector *pos, *norm;
  37. {
  38.     Torus  *torus;
  39.     Vector tmpnrm;
  40.  
  41.     if ((a < EPSILON) || (b < EPSILON)) {
  42.         RLerror(RL_WARN, "Degenerate torus.\n");
  43.         return (Torus *)NULL;
  44.     }
  45.  
  46.     tmpnrm = *norm;
  47.     if (VecNormalize(&tmpnrm) == 0.) {
  48.         RLerror(RL_WARN, "Degenerate torus normal.\n");
  49.         return (Torus *)NULL;
  50.     }
  51.  
  52.     torus = (Torus *)share_malloc(sizeof(Torus));
  53.  
  54.     /*
  55.      * torus->aa holds the square of the swept radius.
  56.      * torus->bb holds the square of the tube radius.
  57.      */
  58.     torus->a = a;
  59.     torus->b = b;
  60.     torus->aa = a*a;
  61.     torus->bb = b*b;
  62.     CoordSysTransform(pos, &tmpnrm, 1., 1., &torus->trans);
  63.  
  64.     return torus;
  65. }
  66.  
  67. /*
  68.  * Ray/torus intersection test.
  69.  */
  70. int
  71. TorusIntersect(torus, inray, mindist, maxdist)
  72. Torus *torus;
  73. Ray *inray;
  74. Float mindist, *maxdist;
  75. {
  76.     Vector pos,ray;
  77.     double c[5],s[4], dist, nmin;
  78.     Float distfactor;
  79.     register int num,i;
  80.  
  81.     TorusTests++;
  82.  
  83.     /* Transform ray into toroid space */
  84.     {
  85.         Ray tmpray;
  86.         tmpray = *inray;
  87.         distfactor = RayTransform(&tmpray, &torus->trans.itrans);
  88.         ray = tmpray.dir;
  89.         pos = tmpray.pos;
  90.         nmin = mindist * distfactor;
  91.     }
  92.  
  93.     /*
  94.  * Original Equations for Toroid with position of (0,0,0) and axis (0,0,1)
  95.  *
  96.  * Equation for two circles of radius b centered at (-a,0,0) and (a,0,0) 
  97.  *
  98.  *      ((R-a)^2 + z*2 - b*b) * ((R+a)^2 + z*z - b*b) = 0 
  99.  *
  100.  *       a         is swept radius
  101.  *       b         is tube  radius
  102.  *
  103.  * subsitute R*R = x*x + y*y  to rotate about z-axis
  104.  *
  105.  * and substitute the parametric ray equations:
  106.  *
  107.  *       x = x0 + t * x1;
  108.  *       y = y0 + t * y1;
  109.  *       z = z0 + t * z1;
  110.  *
  111.  * to get a Quartic in t.
  112.  *
  113.  *       c4*t^4 + c3*t^3 + c2*t^2 + c1*t + c0 = 0
  114.  *
  115.  * where the coefficients are:
  116.  *
  117.  *       c4 =   (x1s + y1s + z1s) * (x1s + y1s + z1s); 
  118.  *       c3 =   4.0 * (tx + ty + tz) * (x1s + y1s + z1s);
  119.  *       c2 =   2.0 * (x1s + y1s + z1s) * (x0s + y0s + z0s - as - bs)
  120.  *            + 4.0 * (tx + ty + tz)    * (tx + ty + tz)
  121.  *            + 4.0 * as * z1s;
  122.  *       c1 =   4.0 * (tx + ty + tz) * (x0s + y0s + z0s - as - bs)
  123.  *            + 8.0 * as * tz;
  124.  *       c0 =   (x0s + y0s + z0s - as - bs) * (x0s + y0s + z0s - as - bs)
  125.  *            + 4.0 * as * (z0s - bs);
  126.  *
  127.  *       as        is swept radius squared
  128.  *       bs        is tube  radius squared
  129.  *      (x0,y0,z0) is origin of ray to be tested
  130.  *      (x1,y1,z1) is direction vector of ray to be tested
  131.  *       tx        is x0 * x1
  132.  *       ty        is y0 * y1
  133.  *       tz        is z0 * z1
  134.  *
  135.  *   Since the direction vector (x1,y1,z1) is normalized:
  136.  *              (x1s + y1s + z1s) = 1.0
  137.  *
  138.  *   Also let     g2s = (x1 * x0) + (y1 * y0) + (z1 * z0)
  139.  *    and let     g0s = (x0 * x0) * (y0 * y0) + (z0 * z0) - as - bs 
  140.  *    since these terms are used fairly often
  141.  */
  142.     {
  143.         register Float g0s,g2s;
  144.         register Float as,bs;
  145.         register Float z0s,z1s,tz;
  146.  
  147.         as  = torus->aa;
  148.         bs  = torus->bb;
  149.         z0s = pos.z * pos.z;
  150.         z1s = ray.z * ray.z;
  151.         tz  = pos.z * ray.z;
  152.         g0s = pos.x * pos.x  +  pos.y * pos.y  +  z0s  -  as  -  bs;
  153.         g2s = pos.x * ray.x  +  pos.y * ray.y  +  tz;
  154.  
  155.         c[4] =   1.0;
  156.         c[3] =   4.0 * g2s;
  157.         c[2] =   2.0 * (g0s  +  2.0 * g2s * g2s  +  2.0 * as * z1s);
  158.         c[1] =   4.0 * (g2s*g0s  +  2.0*as*tz);
  159.         c[0] =   g0s * g0s  +  4.0 * as * (z0s - bs);
  160.     }
  161.  
  162.     /* use GraphGem's Solve Quartic to find roots */
  163.     num = SolveQuartic(c,s);
  164.  
  165.     /* no roots - return 0. */
  166.     if (num==0) return FALSE;
  167.  
  168.     /* of roots return the smallest root > EPSILON */
  169.     dist = 0.0;
  170.     for(i=0;i<num;i++)
  171.     {
  172.         /* if root is in front of ray origin */
  173.         if (s[i] > nmin) {
  174.             /* first valid root */
  175.             if (dist == 0.0) dist = s[i];
  176.             /* else update only if it's closer to ray origin */
  177.             else if (s[i] < dist) dist = s[i];
  178.         }
  179.     }
  180.     dist /= distfactor;
  181.     if (dist > mindist && dist < *maxdist) {
  182.         *maxdist = dist;
  183.         TorusHits++;
  184.         return TRUE;
  185.     }
  186.     return FALSE;
  187. }
  188.  
  189. /*
  190.  * Compute the normal to a torus at a given location on its surface
  191.  */
  192. int
  193. TorusNormal(torus, rawpos, nrm, gnrm)
  194. Torus *torus;
  195. Vector *rawpos, *nrm, *gnrm;
  196. {
  197.     Vector pos;
  198.     register Float dist,posx,posy,xm,ym;
  199.  
  200.     /* Transform intersection point to torus space. */
  201.     pos = *rawpos;
  202.     PointTransform(&pos, &torus->trans.itrans);
  203.  
  204. /* 
  205.  *  The code for the toroid is simpified by always having the axis
  206.  *  be the z-axis and then transforming information to and from
  207.  *  toroid space.
  208.  *
  209.  *  Flatten toroid by ignoring z. Now imagine a knife cutting from
  210.  *  center of toroid to the ray intersection point(x,y). The point
  211.  *  on the tube axis(a circle about the origin with radius 'a') 
  212.  *  where the knife cuts is (xm,ym,zm=0). Unflattening the toroid,
  213.  *  the normal at the point [x,y,z] is (x-xm,y-ym,z). Of course, we
  214.  *  must transform the normal back into world coordinates.
  215.  *  Instead of messing with tan-1,sin and cos, we can find (xm,ym)
  216.  *  by using the proportions:
  217.  *
  218.  *     xm     x           ym     y
  219.  *    ---- = ----   and  ---- = ----
  220.  *     a     dist         a     dist
  221.  *
  222.  *       a         is the swept radius
  223.  *    [x,y,z]      is the point on the toroids surface
  224.  *      dist       is the distance from the z-axis (x*x + y*y).
  225.  *    [xm,ym,zm=0] is the point on the tube's axis 
  226.  *
  227.  */
  228.  
  229.     /* find distance from axis */
  230.     posx = pos.x;
  231.     posy = pos.y;
  232.     dist = sqrt(posx * posx + posy * posy);
  233.  
  234.     if (dist > EPSILON)
  235.     {
  236.         xm = torus->a * posx / dist;
  237.         ym = torus->a * posy / dist;
  238.     }
  239.     else /* ERROR - dist should not be < EPSILON (should never happen)*/
  240.     {
  241.         xm = 0.0;
  242.         ym = 0.0;
  243.     }
  244.  
  245.     /* normal is vector from [xm,ym,zm] to [x,y,z] */
  246.     nrm->x = posx - xm;
  247.     nrm->y = posy - ym;
  248.     nrm->z = pos.z;   /* note by default zm is 0 */
  249.  
  250.     /* Transform normal back to world space. */
  251.     NormalTransform(nrm, &torus->trans.itrans);
  252.     *gnrm = *nrm;
  253.     return FALSE;
  254. }
  255.  
  256. void
  257. TorusUV(torus, pos, norm, uv, dpdu, dpdv)
  258. Torus *torus;
  259. Vector *pos, *norm, *dpdu, *dpdv;
  260. Vec2d *uv;
  261. {
  262.     Vector npos;
  263.     Float costheta, sintheta, rad, cosphi;
  264.  
  265.     npos = *pos;
  266.     PointTransform(&npos, &torus->trans.itrans);
  267.     /*
  268.      * u = theta / 2PI
  269.      */
  270.     rad = sqrt(npos.x*npos.x + npos.y*npos.y);
  271.     costheta = npos.x / rad;
  272.     sintheta = npos.y / rad;
  273.     if (costheta > 1.)    /* roundoff */
  274.         uv->u = 0.;
  275.     else if (costheta < -1.)
  276.         uv->u = 0.5;
  277.     else
  278.         uv->u = acos(costheta) / TWOPI;
  279.     if (sintheta < 0.)
  280.         uv->u = 1. - uv->u;
  281.     if (dpdu) {
  282.         dpdu->x = -npos.y;
  283.         dpdu->y = npos.x;
  284.         dpdu->z = 0.;
  285.         VecTransform(dpdu, &torus->trans.trans);
  286.         (void)VecNormalize(dpdu);
  287.     }
  288.     /*
  289.      * sinphi = npos.z / tor->b;
  290.      * cosphi = rad - tor->a;
  291.      * cosphi is negated in order to make texture 'seam'
  292.      * occur on the interior of the torus.
  293.      */
  294.     cosphi = -(rad - torus->a) / torus->b;
  295.     if (cosphi > 1.)
  296.         uv->v = 0.;
  297.     else if (cosphi < -1.)
  298.         uv->v = 0.5;
  299.     else
  300.         uv->v = acos(cosphi) / TWOPI;
  301.     if (npos.z > 0.)    /* if sinphi > 0... */
  302.         uv->v = 1. - uv->v;
  303.     /*
  304.      * dpdv = norm X dpdu
  305.      */
  306.     if (dpdv) {
  307.         VecCross(norm, dpdu, dpdv);
  308.         VecTransform(dpdv, &torus->trans.trans);
  309.         (void)VecNormalize(dpdv);
  310.     }
  311. }
  312.  
  313. /*
  314.  * Return the extent of a torus.
  315.  */
  316. void
  317. TorusBounds(torus, bounds)
  318. Torus *torus;
  319. Float bounds[2][3];
  320. {
  321.     bounds[LOW][X]  = bounds[LOW][Y]  = -(torus->a+torus->b);
  322.     bounds[HIGH][X] = bounds[HIGH][Y] =  torus->a+torus->b;
  323.     bounds[LOW][Z]  = -torus->b;
  324.     bounds[HIGH][Z] =  torus->b;
  325.     /*
  326.          * Transform bounding box to world space.
  327.          */
  328.     BoundsTransform(&torus->trans.trans, bounds);
  329. }
  330.  
  331. char *
  332. TorusName()
  333. {
  334.     return torusName;
  335. }
  336.  
  337. void
  338. TorusStats(tests, hits)
  339. unsigned long *tests, *hits;
  340. {
  341.     *tests = TorusTests;
  342.     *hits = TorusHits;
  343. }
  344.  
  345. Methods *
  346. TorusMethods()
  347. {
  348.     if (iTorusMethods == NULL) {
  349.         iTorusMethods = MethodsCreate();
  350.         iTorusMethods->create = (GeomCreateFunc *)TorusCreate;
  351.         iTorusMethods->methods  = TorusMethods;
  352.         iTorusMethods->name = TorusName;
  353.         iTorusMethods->intersect = TorusIntersect;
  354.         iTorusMethods->bounds = TorusBounds;
  355.         iTorusMethods->normal = TorusNormal;
  356.         iTorusMethods->uv = TorusUV;
  357.         iTorusMethods->stats = TorusStats;
  358.         iTorusMethods->checkbounds = TRUE;
  359.         iTorusMethods->closed = TRUE;
  360.     }
  361.     return iTorusMethods;
  362. }
  363.  
  364. void
  365. TorusMethodRegister(meth)
  366. UserMethodType meth;
  367. {
  368.     if (iTorusMethods)
  369.         iTorusMethods->user = meth;
  370. }
  371.