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