home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / libray / libobj / sweptsph.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  25.5 KB  |  1,074 lines
  1. /*
  2.  *
  3.  *
  4.  *
  5.  *    sweptsph.c
  6.  *
  7.  *
  8.  *    Copyright (C) 1992 by Lawrence K. Coffin, Craig Colb
  9.  *    All Rights reserved.
  10.  *
  11.  * This software may be freely copied, modified, and redistributed
  12.  * provided that this copyright notice is preserved on all copies.
  13.  *
  14.  * You may not distribute this software, in whole or in part, as part of
  15.  * any commercial product without the express consent of the authors.
  16.  *
  17.  * There is no warranty or other guarantee of fitness of this software
  18.  * for any purpose.  It is provided solely "as is".
  19.  *
  20.  * 06/20/92 Larry Coffin
  21.  * Initial version.
  22.  * 
  23.  */
  24.  
  25. /*
  26.  *    This primative is from "Ray Tracing Objects Defined by Sweeping a
  27.  * Sphere" by J.J. van Wijk in Eurographics '84, pg 73.  It is defined by
  28.  * a sphere of variying radius that is swept or extruded through space along a
  29.  * parametricly defined curve.  The equations for the x,y,z components of the
  30.  * curve and the radius of the sphere are defined by third order equations.
  31.  * This allows for the curve to be defined as the actual equation, as a bezier
  32.  * curve, and perhaps several other curves such as splines. 
  33.  *
  34.  *  The equation for the primitive is the parametric equation:
  35.  *
  36.  *  f(xyz,u) = |xyz - c(u)|^2 - r(u)^2 -- where xyz = {x,y,z} and
  37.  *    c(u) = {x(u),y(u),z(u)}  for umin <= u <= umax.
  38.  *    
  39.  *
  40.  *  Substituting the parametric form of a ray (xyz = Vt + P) where P is the
  41.  *  ray position (ray->pos), and V is the ray direction (ray->dir) we get an
  42.  *  equation
  43.  *
  44.  *  g(u,t) = f2 * t^2 + f1 * t^1 + f0
  45.  *
  46.  *  Where:
  47.  *    f2 = |V|^2 = dotp<ray->dir,ray->dir> = 1 since the ray direction is 
  48.  *        a unit vector.
  49.  *
  50.  *    f1 = 2*dotp<(P - c(u)),V)>
  51.  *
  52.  *    f0 = |P - c(u)|^2 - r(u)^2
  53.  *
  54.  *  Differentiating g(u,t) = 0 with respect to u gives us:
  55.  *
  56.  *    2*f2*t*t' + f1'*t + f1*t' + f0' = 0
  57.  *
  58.  *  Setting t' = 0 gives us:
  59.  *
  60.  *    f1'*t + f0' = 0
  61.  *
  62.  *  Which gives us t = -f1'/f0'
  63.  *
  64.  *  Substituting this back into g(u,t) gives us
  65.  *
  66.  *  g(u) = f2*f0'*f0' - f1*f1'*f0' + f1'*f1'*f0
  67.  *
  68.  *  working this out we get an equation in u of 2*(2*n -1) degree when starting
  69.  *  with center or radius equations of nth degree (10th degree for 3rd degree
  70.  *  center and radius equations).
  71.  *
  72.  *  After solving for u = 0, we can find t from t = -f1'/f0' and keep
  73.  *  the shortest t as our distance of closest intersection.  We must also check
  74.  *  the case where u = umin and u = umax to see if we intersect the sphere at
  75.  *  either end.
  76.  *
  77.  */
  78.  
  79.  
  80. #include "geom.h"
  81. #include "sweptsph.h"
  82.  
  83. #define BIGNUM 9.9e6  /* used as the min and max values for the root finding
  84.             intervals */
  85. #define SMALLNUM .1e-6 /* a number close enough to zero */
  86. #define    DBL_EPS    2.2204460492503131e-10 /* larger than the accuracy of
  87.                             the machine's Floats */ 
  88. #define COUNT 8    /* maximum number of iterations for the Newton's method
  89.                 root finder */ /* 8 best for linear start */
  90. #define COUNT2 200
  91.  
  92. #define sign(x) ((x) < 0 ? (-1) : (1))
  93.  
  94.  
  95. static Methods *iSweptSphMethods = NULL;
  96. static char sweptsphName[] = "sweptsph";
  97.  
  98. unsigned long SweptSphTests, SweptSphHits, maxcount;
  99. Float SweptSphRoot;  /* This holds the value for u for the most recent hit
  100.             This is neccecary inorder to calculate the normal */
  101. int steps[20];
  102. unsigned long int bichecks, bitests, newtons, newchecks;
  103.  
  104. FILE *sweptspherror;
  105.  
  106. /*
  107.  *  SweptSphCreate expects two prameters: the equations for the center and the
  108.  *  equation for the radius.  The center equations are passed as an array of
  109.  *  four vectors.  The center is defined by parametric equations such that
  110.  *
  111.  *  x(u) = cent[0].x + cent[1].x * u + cent[2].x * u^2 + cent[3].x * u^3
  112.  *  y(u) = cent[0].y + cent[1].y * u + cent[2].y * u^2 + cent[3].y * u^3
  113.  *  z(u) = cent[0].z + cent[1].z * u + cent[2].z * u^2 + cent[3].z * u^3
  114.  *
  115.  *  Likewise the radius is defined as:
  116.  *
  117.  *  r(u) = rad[0] + rad[1] * u + rad[2] * u^2 + rad[3] * u^3
  118.  *
  119.  */
  120.  
  121.  
  122. SweptSph *
  123. SweptSphCreate(cent, rad)
  124. Vector cent[4];
  125. Float rad[4];
  126. {
  127.     SweptSph *sweptsph;
  128.     int i;
  129.  
  130.     sweptsph = (SweptSph *)share_malloc(sizeof(SweptSph));
  131.  
  132.     for (i = 0; i < 4; i++){
  133.         if (fabs(cent[i].x) < DBL_EPS) cent[i].x = 0.0;
  134.         if (fabs(cent[i].y) < DBL_EPS) cent[i].y = 0.0;
  135.         if (fabs(cent[i].z) < DBL_EPS) cent[i].z = 0.0;
  136.         if (fabs(rad[i]) < DBL_EPS) rad[i] = 0.0;
  137.     }
  138.  
  139.     sweptsph->a0 = cent[0];
  140.     sweptsph->a1 = cent[1];
  141.     sweptsph->a2 = cent[2];
  142.     sweptsph->a3 = cent[3];
  143.  
  144.     sweptsph->rad[0] = rad[0];
  145.     sweptsph->rad[1] = rad[1];
  146.     sweptsph->rad[2] = rad[2];
  147.     sweptsph->rad[3] = rad[3];
  148.  
  149.     sweptsph->pb[0] = dotp(&(sweptsph->a0),&(sweptsph->a0)) - rad[0]*rad[0];
  150.     sweptsph->pb[1] = 2*dotp(&(sweptsph->a0),&(sweptsph->a1)) - 2*rad[0]*rad[1];
  151.     sweptsph->pb[2] = 2*dotp(&(sweptsph->a0),&(sweptsph->a2))  +   dotp(&(sweptsph->a1),&(sweptsph->a1)) - 2*rad[0]*rad[2] - rad[1]*rad[1];
  152.     sweptsph->pb[3] = 2*dotp(&(sweptsph->a0),&(sweptsph->a3))  + 2*dotp(&(sweptsph->a1),&(sweptsph->a2)) - 2*rad[0]*rad[3] - 2*rad[1]*rad[2];
  153.     sweptsph->pb[4] = 2*dotp(&(sweptsph->a1),&(sweptsph->a3)) +   dotp(&(sweptsph->a2),&(sweptsph->a2)) - 2*rad[1]*rad[3] - rad[2]*rad[2];
  154.     sweptsph->pb[5] = 2*dotp(&(sweptsph->a2),&(sweptsph->a3)) - 2*rad[2]*rad[3];
  155.     sweptsph->pb[6] = dotp(&(sweptsph->a3),&(sweptsph->a3)) - rad[3]*rad[3];
  156.  
  157.     sweptsph->pd[0] = 16*(sweptsph->pb[4])*(sweptsph->pb[4]);
  158.     sweptsph->pd[1] = 40*(sweptsph->pb[4])*(sweptsph->pb[5]);
  159.     sweptsph->pd[2] = 48*(sweptsph->pb[4])*(sweptsph->pb[6]) + 25*(sweptsph->pb[5])*(sweptsph->pb[5]);
  160.     sweptsph->pd[3] = 60*(sweptsph->pb[5])*(sweptsph->pb[6]);
  161.     sweptsph->pd[4] = 36*(sweptsph->pb[6])*(sweptsph->pb[6]);
  162.  
  163. /*    if ((sweptspherror = fopen("/usr/people/c/lcoffin/sweptsph.stat","w")) == NULL){
  164.         fprintf(stderr,"Unable to open stat file for sweptsph\n");
  165.         exit(1);
  166.     }
  167. */
  168.     return sweptsph;
  169. }
  170.  
  171. Methods *
  172. SweptSphMethods()
  173. {
  174.     if (iSweptSphMethods == (Methods *)NULL) {
  175.         iSweptSphMethods = MethodsCreate();
  176.         iSweptSphMethods->create = (GeomCreateFunc *)SweptSphCreate;
  177.         iSweptSphMethods->methods = SweptSphMethods;
  178.         iSweptSphMethods->name = SweptSphName;
  179.         iSweptSphMethods->intersect = SweptSphIntersect;
  180.         iSweptSphMethods->normal = SweptSphNormal;
  181.         iSweptSphMethods->bounds = SweptSphBounds;
  182.         iSweptSphMethods->stats = SweptSphStats;
  183.         iSweptSphMethods->checkbounds = TRUE;
  184.         iSweptSphMethods->closed = TRUE;
  185.     }
  186.     return iSweptSphMethods;
  187. }
  188.  
  189. int
  190. SweptSphIntersect(sweptsph, ray, mindist, maxdist)
  191. SweptSph *sweptsph;
  192. Ray *ray;
  193. Float mindist, *maxdist;
  194. {
  195.     Vector a0, a1, a2, a3, dir, pos;
  196.     Float b[7],c[4], d[11], e[6], h[5], i[11], j[11], g[11], roots[10], t;
  197.     Float rad[4], div, radic, root;
  198.     int num, hit, x, k;
  199.  
  200.     SweptSphTests++;
  201.  
  202.     dir = ray->dir;
  203.     pos = ray->pos;
  204.  
  205.     rad[0] = sweptsph->rad[0];
  206.     rad[1] = sweptsph->rad[1];
  207.     rad[2] = sweptsph->rad[2];
  208.     rad[3] = sweptsph->rad[3];
  209.  
  210.     a0 = sweptsph->a0;
  211.     a1 = sweptsph->a1;
  212.     a2 = sweptsph->a2;
  213.     a3 = sweptsph->a3;
  214.  
  215.     /* b[] = f0 */
  216.     b[0] =    dotp(&pos,&pos) - 2*dotp(&pos,&a0) + sweptsph->pb[0];
  217.     b[1] = -2*dotp(&pos,&a1 ) + sweptsph->pb[1];
  218.     b[2] = -2*dotp(&pos,&a2 ) + sweptsph->pb[2];
  219.     b[3] = -2*dotp(&pos,&a3 ) + sweptsph->pb[3];
  220.     b[4] = sweptsph->pb[4];
  221.     b[5] = sweptsph->pb[5];
  222.     b[6] = sweptsph->pb[6];
  223.  
  224.     /* c[] = f1 */
  225.     c[0] =  2*dotp(&pos,&dir) - 2*dotp(&a0,&dir);
  226.     c[1] = -2*dotp(&a1 ,&dir);
  227.     c[2] = -2*dotp(&a2 ,&dir);
  228.     c[3] = -2*dotp(&a3 ,&dir);
  229.  
  230.     /* d[] = f0'*f0' which also equals f2*f0'*f0' since our ray->dir is
  231.         a unit vector */
  232.     d[0] =     b[1]*b[1];
  233.     d[1] =   4*b[1]*b[2];
  234.     d[2] =   6*b[1]*b[3] +  4*b[2]*b[2];
  235.     d[3] =   8*b[1]*b[4] + 12*b[2]*b[3];
  236.     d[4] =  10*b[1]*b[5] + 16*b[2]*b[4] +  9*b[3]*b[3];
  237.     d[5] =  12*b[1]*b[6] + 20*b[2]*b[5] + 24*b[3]*b[4];
  238.     d[6] =  24*b[2]*b[6] + 30*b[3]*b[5] + sweptsph->pd[0];
  239.     d[7] =  36*b[3]*b[6] + sweptsph->pd[1];
  240.     d[8] =  sweptsph->pd[2];
  241.     d[9] =  sweptsph->pd[3];
  242.     d[10] = sweptsph->pd[4];
  243.  
  244.     /* e = f1*f1' */
  245.     e[0] =   c[0]*c[1];
  246.     e[1] = 2*c[0]*c[2] +   c[1]*c[1];
  247.     e[2] = 3*c[0]*c[3] + 3*c[1]*c[2];
  248.     e[3] = 4*c[1]*c[3] + 2*c[2]*c[2];
  249.     e[4] = 5*c[2]*c[3];
  250.     e[5] = 3*c[3]*c[3];
  251.  
  252.     /* h[] = f1'*f1' */
  253.     h[0] =    c[1]*c[1];
  254.     h[1] =  4*c[1]*c[2];
  255.     h[2] =  6*c[1]*c[3] + 4*c[2]*c[2];
  256.     h[3] = 12*c[2]*c[3];
  257.     h[4] =  9*c[3]*c[3];
  258.  
  259.     /* i[] = f1*f1'*f0' */
  260.     i[0] =    e[0]*b[1];
  261.     i[1] =  2*e[0]*b[2] +   e[1]*b[1];
  262.     i[2] =  3*e[0]*b[3] + 2*e[1]*b[2] +   e[2]*b[1];
  263.     i[3] =  4*e[0]*b[4] + 3*e[1]*b[3] + 2*e[2]*b[2] +   e[3]*b[1];
  264.     i[4] =  5*e[0]*b[5] + 4*e[1]*b[4] + 3*e[2]*b[3] + 2*e[3]*b[2] +   e[4]*b[1];
  265.     i[5] =  6*e[0]*b[6] + 5*e[1]*b[5] + 4*e[2]*b[4] + 3*e[3]*b[3] + 2*e[4]*b[2] + e[5]*b[1];
  266.     i[6] =  6*e[1]*b[6] + 5*e[2]*b[5] + 4*e[3]*b[4] + 3*e[4]*b[3] + 2*e[5]*b[2];
  267.     i[7] =  6*e[2]*b[6] + 5*e[3]*b[5] + 4*e[4]*b[4] + 3*e[5]*b[3];
  268.     i[8] =  6*e[3]*b[6] + 5*e[4]*b[5] + 4*e[5]*b[4];
  269.     i[9] =  6*e[4]*b[6] + 5*e[5]*b[5];
  270.     i[10] = 6*e[5]*b[6];
  271.  
  272.     /* j[] = f1'*f1'*f0 */
  273.     j[0] =  h[0]*b[0];
  274.     j[1] =  h[0]*b[1] + h[1]*b[0];
  275.     j[2] =  h[0]*b[2] + h[1]*b[1] + h[2]*b[0];
  276.     j[3] =  h[0]*b[3] + h[1]*b[2] + h[2]*b[1] + h[3]*b[0];
  277.     j[4] =  h[0]*b[4] + h[1]*b[3] + h[2]*b[2] + h[3]*b[1] + h[4]*b[0];
  278.     j[5] =  h[0]*b[5] + h[1]*b[4] + h[2]*b[3] + h[3]*b[2] + h[4]*b[1];
  279.     j[6] =  h[0]*b[6] + h[1]*b[5] + h[2]*b[4] + h[3]*b[3] + h[4]*b[2];
  280.     j[7] =  h[1]*b[6] + h[2]*b[5] + h[3]*b[4] + h[4]*b[3];
  281.     j[8] =  h[2]*b[6] + h[3]*b[5] + h[4]*b[4];
  282.     j[9] =  h[3]*b[6] + h[4]*b[5];
  283.     j[10] = h[4]*b[6];
  284.  
  285.     /* g[] = f2*f0'*f0' - f1*f1'*f0' + f1'*f1'*f0 */
  286.     for (x = 0; x <= 10; x++){
  287.         g[x] = d[x] - i[x] + j[x];
  288.     }
  289.  
  290.     /* find the roots */
  291.     num = FindRoots(g,10,0.0,1.0,roots);
  292.  
  293.     hit = FALSE;
  294.  
  295.     /* check each root */
  296.     for (k = 0; k < num; k++){
  297.         /* check to see if f1' = 0 */
  298.         div = (c[1] + (2*c[2] + 3*c[3]*roots[k])*roots[k]);
  299.  
  300.         if (div){
  301.             /* find t = -f0'/f1' */
  302.             t = -(b[1] + (2*b[2] + (3*b[3] + (4*b[4] + (5*b[5]
  303.                 + 6*b[6]*roots[k])*roots[k])*roots[k])*roots[k])*roots[k])/div;
  304.             /* set maxdist = t if t is smaller than our 
  305.                 maximum distance */
  306.             if (t > mindist + .0001){
  307.                 if (t < *maxdist){
  308.                     SweptSphRoot = roots[k];
  309.                     *maxdist = t;
  310.                     hit = TRUE;
  311.                 }
  312.             }
  313.         }
  314.     }
  315.  
  316.     /* check intersection with the end spheres */
  317.     g[2] = 1;
  318.     g[1] = c[0];
  319.     g[0] = b[0];
  320.  
  321.     radic = g[1]*g[1] - 4*g[0];
  322.  
  323.     if (radic > 0){
  324.         radic = sqrt(radic);
  325.         root = (-g[1] - radic)/2.0;
  326.  
  327.         if (root > mindist + .0001 && root < *maxdist){
  328.             SweptSphRoot = 0.0;
  329.             *maxdist = root;
  330.             hit = TRUE;
  331.         }
  332.  
  333.         root = (-g[1] + radic)/2.0;
  334.  
  335.         if (root > mindist + .0001 && root < *maxdist){
  336.             SweptSphRoot = 0.0;
  337.             *maxdist = root;
  338.             hit = TRUE;
  339.         }
  340.     }
  341.     else if (radic == 0){
  342.         root = -g[1]/2.0;
  343.         if (root > mindist + .0001 && root < *maxdist){
  344.             SweptSphRoot = 0.0;
  345.             *maxdist = root;
  346.             hit = TRUE;
  347.         }
  348.     }        
  349.  
  350.  
  351.     g[2] = 1;
  352.     g[1] = c[0]+c[1]+c[2]+c[3];
  353.     g[0] = b[0]+b[1]+b[2]+b[3]+b[4]+b[5]+b[6];
  354.  
  355.     radic = g[1]*g[1] - 4*g[0];
  356.  
  357.     if (radic > 0){
  358.         radic = sqrt(radic);
  359.         root = (-g[1] - radic)/2.0;
  360.  
  361.         if (root > mindist + .0001 && root < *maxdist){
  362.             SweptSphRoot = 1.0;
  363.             *maxdist = root;
  364.             hit = TRUE;
  365.         }
  366.  
  367.         root = (-g[1] + radic)/2.0;
  368.  
  369.         if (root > mindist + .0001 && root < *maxdist){
  370.             SweptSphRoot = 1.0;
  371.             *maxdist = root;
  372.             hit = TRUE;
  373.         }
  374.     }
  375.     else if (radic == 0){
  376.         root = -g[1]/2.0;
  377.         if (root > mindist + .0001 && root < *maxdist){
  378.             SweptSphRoot = 1.0;
  379.             *maxdist = root;
  380.             hit = TRUE;
  381.         }
  382.     }        
  383.  
  384.     if (hit == TRUE){
  385.         SweptSphHits++;
  386.     }
  387.     return hit;
  388. }
  389.  
  390.  
  391. int
  392. SweptSphNormal(sweptsph, pos, nrm, gnrm)
  393. SweptSph *sweptsph;
  394. Vector *pos, *nrm, *gnrm;
  395. {
  396.     Vector a0, a1, a2, a3;
  397.  
  398.     a0 = sweptsph->a0;
  399.     a1 = sweptsph->a1;
  400.     a2 = sweptsph->a2;
  401.     a3 = sweptsph->a3;
  402.  
  403.     /* normal is simply the normal to the sphere at c(u): u = SweptSphRoot */
  404.     nrm->x = pos->x - a0.x - (a1.x + (a2.x + a3.x*SweptSphRoot)*SweptSphRoot)*SweptSphRoot;
  405.     nrm->y = pos->y - a0.y - (a1.y + (a2.y + a3.y*SweptSphRoot)*SweptSphRoot)*SweptSphRoot;
  406.     nrm->z = pos->z - a0.z - (a1.z + (a2.z + a3.z*SweptSphRoot)*SweptSphRoot)*SweptSphRoot;
  407.  
  408.     VecNormalize(nrm);
  409.     *gnrm = *nrm;
  410.  
  411.     return FALSE;
  412. }
  413.  
  414. void
  415. SweptSphBounds(sweptsph, bounds)
  416. SweptSph *sweptsph;
  417. Float bounds[2][3];
  418. {
  419.     Float a[4], aprime[3],root, val, rad, radimax;
  420.  
  421.     /* for the bounding boxes, calculate the maximum radius of the swept
  422.         sphere then set the bounding box to the max and min of the
  423.         center curve plus or minus the maximum radius of the sphere
  424.     */
  425.  
  426.     /* find the maximum radius */
  427.  
  428.     radimax = 0;
  429.  
  430.     a[0] = sweptsph->rad[0];
  431.     a[1] = sweptsph->rad[1];
  432.     a[2] = sweptsph->rad[2];
  433.     a[3] = sweptsph->rad[3];
  434.  
  435.     aprime[0] = a[1];
  436.     aprime[1] = 2.0*a[2];
  437.     aprime[2] = 3.0*a[3];
  438.  
  439.     if (aprime[2] == 0){
  440.         if (aprime[1] != 0){
  441.             root = -aprime[0]/aprime[1];
  442.             if (root > 0.0 && root < 1.0){
  443.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  444.  
  445.                 if (val > radimax){
  446.                     radimax = val;
  447.                 }
  448.             }
  449.         }
  450.     }
  451.     else {
  452.         rad = aprime[1]*aprime[1] - 4.0 * aprime[2] *  aprime[0];
  453.         if (rad == 0) {
  454.             root = -aprime[1]/(2.0*aprime[2]);
  455.             if (root > 0.0 && root < 1.0){
  456.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  457.  
  458.                 if (val > radimax){
  459.                     radimax = val;
  460.                 }
  461.             }
  462.         }
  463.         else if (rad > 0){
  464.             rad = sqrt(rad);
  465.  
  466.             root = (-aprime[1] - rad)/(2.0*aprime[2]);
  467.  
  468.             if (root > 0.0 && root < 1.0){
  469.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  470.  
  471.                 if (val > radimax){
  472.                     radimax = val;
  473.                 }
  474.             }
  475.  
  476.             root = (-aprime[1] + rad)/(2.0*aprime[2]);
  477.  
  478.             if (root > 0.0 && root < 1.0){
  479.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  480.  
  481.                 if (val > radimax){
  482.                     radimax = val;
  483.                 }
  484.             }
  485.         }
  486.     }
  487.  
  488.     val = a[0];
  489.     if (val > radimax){
  490.         radimax = val;
  491.     }
  492.  
  493.     val = a[0] + a[1] + a[2] + a[3];
  494.     if (val > radimax){
  495.         radimax = val;
  496.     }
  497.  
  498.  
  499.     /* find the x bounds of the curve */
  500.     a[0] = sweptsph->a0.x;
  501.     a[1] = sweptsph->a1.x;
  502.     a[2] = sweptsph->a2.x;
  503.     a[3] = sweptsph->a3.x;
  504.  
  505.     aprime[0] = a[1];
  506.     aprime[1] = 2.0*a[2];
  507.     aprime[2] = 3.0*a[3];
  508.  
  509.  
  510.     if (aprime[2] == 0){
  511.         if (aprime[1] != 0){
  512.             root = -aprime[0]/aprime[1];
  513.             if (root > 0.0 && root < 1.0){
  514.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  515.  
  516.                 if (val - radimax < bounds[LOW][X]){
  517.                     bounds[LOW][X] = val - radimax;
  518.                 }
  519.                 if (val + radimax > bounds[HIGH][X]){
  520.                     bounds[HIGH][X] = val + radimax;
  521.                 }
  522.             }
  523.         }
  524.     }
  525.     else {
  526.         rad = aprime[1]*aprime[1] - 4.0 * aprime[2] *  aprime[0];
  527.  
  528.         if (rad == 0) {
  529.             root = -aprime[1]/(2.0*aprime[2]);
  530.             if (root > 0.0 && root < 1.0){
  531.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  532.     
  533.                 if (val - radimax < bounds[LOW][X]){
  534.                     bounds[LOW][X] = val - radimax;
  535.                 }
  536.                 if (val + radimax > bounds[HIGH][X]){
  537.                     bounds[HIGH][X] = val + radimax;
  538.                 }
  539.             }
  540.         }
  541.         else if (rad > 0){
  542.             rad = sqrt(rad);
  543.     
  544.             root = (-aprime[1] - rad)/(2.0*aprime[2]);
  545.     
  546.             if (root > 0.0 && root < 1.0){
  547.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  548.     
  549.                 if (val - radimax < bounds[LOW][X]){
  550.                     bounds[LOW][X] = val - radimax;
  551.                 }
  552.                 if (val + radimax > bounds[HIGH][X]){
  553.                     bounds[HIGH][X] = val + radimax;
  554.                 }
  555.             }
  556.     
  557.             root = (-aprime[1] + rad)/(2.0*aprime[2]);
  558.     
  559.             if (root > 0.0 && root < 1.0){
  560.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  561.     
  562.                 if (val - radimax < bounds[LOW][X]){
  563.                     bounds[LOW][X] = val - radimax;
  564.                 }
  565.                 if (val + radimax > bounds[HIGH][X]){
  566.                     bounds[HIGH][X] = val + radimax;
  567.                 }
  568.             }
  569.         }
  570.     }
  571.  
  572.     val = a[0];
  573.     if (val - radimax < bounds[LOW][X]){
  574.         bounds[LOW][X] = val - radimax;
  575.     }
  576.     if (val + radimax > bounds[HIGH][X]){
  577.         bounds[HIGH][X] = val + radimax;
  578.     }
  579.     val = a[0] + a[1] + a[2] + a[3];
  580.     if (val - radimax < bounds[LOW][X]){
  581.         bounds[LOW][X] = val - radimax;
  582.     }
  583.     if (val + radimax > bounds[HIGH][X]){
  584.         bounds[HIGH][X] = val + radimax;
  585.     }
  586.  
  587.  
  588.     /* find the y bounds of the curve */
  589.     a[0] = sweptsph->a0.y;
  590.     a[1] = sweptsph->a1.y;
  591.     a[2] = sweptsph->a2.y;
  592.     a[3] = sweptsph->a3.y;
  593.  
  594.     aprime[0] = a[1];
  595.     aprime[1] = 2.0*a[2];
  596.     aprime[2] = 3.0*a[3];
  597.  
  598.     if (aprime[2] == 0){
  599.         if (aprime[1] != 0){
  600.             root = -aprime[0]/aprime[1];
  601.             if (root > 0.0 && root < 1.0){
  602.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  603.  
  604.                 if (val - radimax < bounds[LOW][Y]){
  605.                     bounds[LOW][Y] = val - radimax;
  606.                 }
  607.                 if (val + radimax > bounds[HIGH][Y]){
  608.                     bounds[HIGH][Y] = val + radimax;
  609.                 }
  610.             }
  611.         }
  612.     }
  613.     else {
  614.         rad = aprime[1]*aprime[1] - 4.0 * aprime[2] *  aprime[0];
  615.         if (rad == 0) {
  616.             root = -aprime[1]/(2.0*aprime[2]);
  617.             if (root > 0.0 && root < 1.0){
  618.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  619.     
  620.                 if (val - radimax < bounds[LOW][Y]){
  621.                     bounds[LOW][Y] = val - radimax;
  622.                 }
  623.                 if (val + radimax > bounds[HIGH][Y]){
  624.                     bounds[HIGH][Y] = val + radimax;
  625.                 }
  626.             }
  627.         }
  628.         else if (rad > 0){
  629.             rad = sqrt(rad);
  630.             root = (-aprime[1] - rad)/(2.0*aprime[2]);
  631.             if (root > 0.0 && root < 1.0){
  632.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  633.     
  634.                 if (val - radimax < bounds[LOW][Y]){
  635.                     bounds[LOW][Y] = val - radimax;
  636.                 }
  637.                 if (val + radimax > bounds[HIGH][Y]){
  638.                     bounds[HIGH][Y] = val + radimax;
  639.                 }
  640.             }
  641.             root = (-aprime[1] + rad)/(2.0*aprime[2]);
  642.             if (root > 0.0 && root < 1.0){
  643.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;    
  644.     
  645.                 if (val - radimax < bounds[LOW][Y]){
  646.                     bounds[LOW][Y] = val - radimax;
  647.                 }
  648.                 if (val + radimax > bounds[HIGH][Y]){
  649.                     bounds[HIGH][Y] = val + radimax;
  650.                 }
  651.             }
  652.         }
  653.     }
  654.  
  655.     val = a[0];
  656.     if (val - radimax < bounds[LOW][Y]){
  657.         bounds[LOW][Y] = val - radimax;
  658.     }
  659.     if (val + radimax > bounds[HIGH][Y]){
  660.         bounds[HIGH][Y] = val + radimax;
  661.     }
  662.     val = a[0] + a[1] + a[2] + a[3];
  663.     if (val - radimax < bounds[LOW][Y]){
  664.         bounds[LOW][Y] = val - radimax;
  665.     }
  666.     if (val + radimax > bounds[HIGH][Y]){
  667.         bounds[HIGH][Y] = val + radimax;
  668.     }
  669.  
  670.  
  671.     /* find the z bounds of the curve */
  672.     a[0] = sweptsph->a0.z;
  673.     a[1] = sweptsph->a1.z;
  674.     a[2] = sweptsph->a2.z;
  675.     a[3] = sweptsph->a3.z;
  676.  
  677.     aprime[0] = a[1];
  678.     aprime[1] = 2.0*a[2];
  679.     aprime[2] = 3.0*a[3];
  680.  
  681.     if (aprime[2] == 0){
  682.         if (aprime[1] != 0){
  683.             root = -aprime[0]/aprime[1];
  684.             if (root > 0.0 && root < 1.0){
  685.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  686.  
  687.                 if (val - radimax < bounds[LOW][Z]){
  688.                     bounds[LOW][Z] = val - radimax;
  689.                 }
  690.                 if (val + radimax > bounds[HIGH][Z]){
  691.                     bounds[HIGH][Z] = val + radimax;
  692.                 }
  693.             }
  694.         }
  695.     }
  696.     else {
  697.         rad = aprime[1]*aprime[1] - 4.0 * aprime[2] *  aprime[0];
  698.  
  699.         if (rad == 0) {
  700.             root = -aprime[1]/(2.0*aprime[2]);
  701.             if (root > 0.0 && root < 1.0){
  702.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  703.     
  704.                 if (val - radimax < bounds[LOW][Z]){
  705.                     bounds[LOW][Z] = val - radimax;
  706.                 }
  707.                 if (val + radimax > bounds[HIGH][Z]){
  708.                     bounds[HIGH][Z] = val + radimax;
  709.                 }
  710.             }
  711.         }
  712.         else if (rad > 0){
  713.             rad = sqrt(rad);
  714.             root = (-aprime[1] - rad)/(2.0*aprime[2]);
  715.             if (root > 0.0 && root < 1.0){
  716.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  717.  
  718.                 if (val - radimax < bounds[LOW][Z]){
  719.                     bounds[LOW][Z] = val - radimax;
  720.                 }
  721.                 if (val + radimax > bounds[HIGH][Z]){
  722.                     bounds[HIGH][Z] = val + radimax;
  723.                 }
  724.             }
  725.             root = (-aprime[1] + rad)/(2.0*aprime[2]);
  726.             if (root > 0.0 && root < 1.0){
  727.                 val = a[0] + (a[1] + (a[2] + a[3]*root)*root)*root;
  728.     
  729.                 if (val - radimax < bounds[LOW][Z]){
  730.                     bounds[LOW][Z] = val - radimax;
  731.                 }
  732.                 if (val + radimax > bounds[HIGH][Z]){
  733.                     bounds[HIGH][Z] = val + radimax;
  734.                 }
  735.             }
  736.         }
  737.     }
  738.  
  739.     val = a[0];
  740.     if (val - radimax < bounds[LOW][Z]){
  741.         bounds[LOW][Z] = val - radimax;
  742.     }
  743.     if (val + radimax > bounds[HIGH][Z]){
  744.         bounds[HIGH][Z] = val + radimax;
  745.     }
  746.     val = a[0] + a[1] + a[2] + a[3];
  747.     if (val - radimax < bounds[LOW][Z]){
  748.         bounds[LOW][Z] = val - radimax;
  749.     }
  750.     if (val + radimax > bounds[HIGH][Z]){
  751.         bounds[HIGH][Z] = val + radimax;
  752.     }
  753.  
  754.  
  755. }
  756.  
  757. void
  758. SweptSphStats(tests, hits)
  759. unsigned long *tests, *hits;
  760. {
  761.     *tests = SweptSphTests;
  762.     *hits = SweptSphHits;
  763. }
  764.  
  765. char *SweptSphName()
  766. {
  767.     return sweptsphName;
  768. }
  769.  
  770. /* evaluate a polynomial at x */
  771.  
  772. Float evalpoly(p,order,x)
  773. Float *p,x;
  774. int order;
  775. {
  776.     Float val;
  777.     int i;
  778.  
  779.     val = p[order];
  780.  
  781.     for (i = order - 1; i >= 0; i--){
  782.         val = p[i] + val*x;
  783.     }
  784.     return val;
  785.  
  786. }
  787.  
  788. /* use bisection to find the root */
  789. int findit(p,order,minx,maxx,root)
  790. Float *p,minx,maxx,*root;
  791. int order;
  792. {
  793.     Float minval, val, midx;
  794.     int i;
  795.  
  796.     minval = evalpoly(p,order,minx);
  797.  
  798. /*    bichecks++;
  799.  
  800.     if((newchecks+bichecks)%10000 == 0){
  801.         fprintf(stderr,"Checks: %i(%i) Bisection: %.2f Newton: %.2f Both %.2f (%.2f) \n",bichecks+newchecks,newchecks, (Float)bitests/(Float)bichecks,(Float)newtons/(Float)newchecks,(Float)(bitests+newtons)/(Float)(bichecks+newchecks),(Float)(bitests+newtons)/(Float)newchecks);
  802.         fflush(stderr);
  803.     }
  804. */
  805.     i = 0;
  806.     while(fabs(minx-maxx) > DBL_EPS){
  807.  
  808.         midx = (maxx + minx)/2.0;
  809.  
  810.         val = evalpoly(p,order,midx);
  811.  
  812.         if (fabs(val) < SMALLNUM) break;
  813.  
  814.         if(val*minval < 0){
  815.             maxx = midx;
  816.         }
  817.         else {
  818.             minx = midx;
  819.             minval = val;
  820.         }
  821. /*        bitests++;
  822. */        i++;
  823.         if (i > COUNT2) {
  824.             maxcount++;
  825.             break;
  826.         }
  827.     }
  828.  
  829.     *root = midx;
  830.     return TRUE;
  831.  
  832. }
  833.  
  834. /* use Newton's method to find the root */
  835.  
  836. int newtonit(p,order,minx,maxx,root)
  837. Float *p,minx,maxx,*root;
  838. int order;
  839. {
  840.     Float val, x, *pprime, oldval, oldx, m;
  841.     int i;
  842.  
  843. /*
  844.     newchecks++;
  845.  
  846.     if((newchecks+bichecks)%10000 == 0){
  847.         fprintf(stderr,"Checks: %i(%i) Bisection: %.2f Newton: %.2f Both %.2f (%.2f) \n",bichecks+newchecks,newchecks, (Float)bitests/(Float)bichecks,(Float)newtons/(Float)newchecks,(Float)(bitests+newtons)/(Float)(bichecks+newchecks),(Float)(bitests+newtons)/(Float)newchecks);
  848.         fflush(stderr);
  849.     }
  850. */
  851.     pprime = (Float *)malloc(order*sizeof(Float));
  852.  
  853.     for (i = 0; i < order; i++){
  854.         pprime[i] = (i+1)*p[i+1];
  855.     }
  856.  
  857.     oldx = maxx;
  858.     oldval = evalpoly(p,order,maxx);
  859.  
  860.     /* first point is where the line through the interval end points
  861.         equals zero */
  862.     m = (oldval - evalpoly(p,order,minx))/(maxx - minx);
  863.     x = (m*maxx - oldval)/m;
  864.     val = evalpoly(p,order,x);
  865.  
  866.     
  867.     i = 0;
  868.  
  869.     while(fabs(val) > SMALLNUM || fabs(oldx-x) > DBL_EPS){
  870.         oldx = x;
  871.         oldval = val;
  872.         x = x - val/evalpoly(pprime,order-1,x);
  873.         val = evalpoly(p,order,x);
  874.         /* exit if the new x is outside [umin,umax] */
  875.         if (x > maxx || x < minx){
  876.             free(pprime);
  877.             return FALSE;
  878.         }
  879. /*        newtons++;
  880. */        i++;
  881.         /* exit if we haven't found the root after COUNT iterations
  882.             (in case we are at a point where we are bouncing
  883.             back and forth without getting closer to the root.
  884.             Perhaps a different test would be better)
  885.         */
  886.         if (i > COUNT){
  887.             free(pprime);
  888.             return FALSE;
  889.         }
  890.     }
  891.  
  892.     free(pprime);
  893.     *root = x;
  894.     return TRUE;
  895.  
  896. }
  897.  
  898.  
  899.  
  900. /* FindRoots finds the roots to a polynomial by reducing the polynomial
  901.     down to it's 2nd order derivative then using the roots of the derivative
  902.     as the max,min points between which there can be only one zero point, if
  903.     any, for the next higher derivative.
  904.  
  905.     Roots for greater than 2nd order polynomials are found by first trying
  906.     Newton's method and if that fails, then using bisection.
  907.  
  908.     This algorithm probably could use alot of work to make it faster than
  909.     it currently is.  I was not particularly sure what to use as determining
  910.     factors for the Newton's methods failure.  One was finding a new x value
  911.     outside the current search range, which would result in finding the
  912.     wrong root, the other was an upper limit on the number of iterations
  913.     performed.  This last test was to insure that we did not get caught in
  914.     a loop where each sucssesive iteration returned us to the previous x
  915.     value.
  916. */
  917.  
  918. int FindRoots(p,order,lower,upper,roots)
  919. Float *p,lower, upper;
  920. int order;
  921. Float *roots;
  922. {
  923.     Float *pprime, *droots, rad, root, lowerval,upperval, rootval, oldrootval;
  924.     int num, i, j, minindex, found;
  925.  
  926.     if (p[order] < SMALLNUM){
  927.         return FindRoots(p,order-1,lower,upper,roots);
  928.     }
  929.     if (order == 0) return 0;
  930.     if (order == 1){
  931.         root = -p[0]/p[1];
  932.         if (root > lower && root < upper){
  933.             roots[0] = root;
  934.             return 1;
  935.         }
  936.         return 0;
  937.     }
  938.  
  939.     found = 0;
  940.     if (order == 2) {
  941.         rad = p[1]*p[1] - 4*p[0]*p[2];
  942.         if (rad < 0) return 0;
  943.         if (rad == 0) {
  944.             root = -p[1]/(2.0*p[2]);
  945.             if (root > lower && root < upper){
  946.                 roots[0] = root;
  947.                 return 1;
  948.             }
  949.             return 0;
  950.         }
  951.  
  952.         rad = sqrt(rad);
  953.         root = (-p[1]-rad)/(2.0*p[2]);
  954.         if (root > lower && root < upper){
  955.             roots[found] = root;
  956.             found++;
  957.         }
  958.         root = (-p[1]+rad)/(2.0*p[2]);
  959.         if (root > lower && root < upper){
  960.             roots[found] = root;
  961.             found++;
  962.         }
  963.         if (found == 2){
  964.             if (roots[0] > roots[1]){
  965.                 root = roots[0];
  966.                 roots[0] = roots[1];
  967.                 roots[1] = root;
  968.             }
  969.         }
  970.         return found;
  971.     }
  972.  
  973.     droots = (Float*)malloc((order-1)*sizeof(Float));
  974.     pprime = (Float*)malloc((order)*sizeof(Float));
  975.  
  976.     if (droots == NULL || pprime == NULL){
  977.         fprintf(stderr,"\n\nNot enough memory\n\n");
  978.         exit (1);
  979.     }
  980.  
  981.     pprime[0] = p[1];
  982.  
  983.     /* calculate the derivative of p */
  984.     for (i = 1; i < order; i++){
  985.         pprime[i] = (i+1)*p[i+1];
  986.     }
  987.  
  988.     /* find the roots of the derivative which will be p's max and min
  989.         points */
  990.     num = FindRoots(pprime,order-1, lower, upper,droots);
  991.  
  992.     free(pprime);
  993.  
  994.     lowerval = evalpoly(p,order,lower);
  995.     upperval = evalpoly(p,order,upper);
  996.  
  997.     /* no max,min found check the interval (lower,upper) */
  998.     if (num == 0){
  999.         free(droots);
  1000.         if (lowerval*upperval < 0){
  1001.             if (newtonit(p,order,lower,upper,&roots[0]) == FALSE){
  1002.                 if (findit(p,order,lower,upper,&roots[0]) == FALSE){
  1003.                     fprintf(stderr,"Neither found the root!!");
  1004.                     return 0;
  1005.                 }
  1006.                 else {
  1007.                     return 1;
  1008.                 }
  1009.             }
  1010.             return 1;
  1011.         }
  1012.         return 0;
  1013.     }
  1014.  
  1015.     found = 0;
  1016.  
  1017.     /* check the interval (lower,1st_root] */
  1018.     rootval = evalpoly(p,order,droots[0]);
  1019.  
  1020.     if (fabs(rootval) < SMALLNUM){
  1021.         roots[found] = droots[0];
  1022.         found++;
  1023.     }
  1024.     else if (lowerval*rootval < 0){
  1025.         if (newtonit(p,order,lower,droots[0],&root) == TRUE) {
  1026.                 roots[found] = root;
  1027.                 found++;
  1028.         }
  1029.         else if (findit(p,order,lower,droots[0],&root) == TRUE) {
  1030.                 roots[found] = root;
  1031.                 found++;
  1032.         }
  1033.         
  1034.     }
  1035.  
  1036.     /*  check the intervals (nth_root,(n+1)th_root] for n = 1 to (num - 1) */
  1037.     for(i = 1 ; i < num; i++){
  1038.         oldrootval = rootval;
  1039.         rootval = evalpoly(p,order,droots[i]);
  1040.         if (fabs(rootval) < SMALLNUM){
  1041.             roots[found] = droots[i];
  1042.             found++;
  1043.         }
  1044.         else if (rootval*oldrootval < 0){
  1045.             if (newtonit(p,order,droots[i-1],droots[i],&root) == TRUE) {
  1046.                 roots[found] = root;
  1047.                 found++;
  1048.             }
  1049.             else if (findit(p,order,droots[i-1],droots[i],&root) == TRUE) {
  1050.                 roots[found] = root;
  1051.                 found++;
  1052.             }
  1053.         }
  1054.     }
  1055.  
  1056.     /* check the interval ((num)th_root,upper) */
  1057.     if (fabs(rootval) > SMALLNUM){
  1058.         if (rootval*upperval < 0){
  1059.             if (newtonit(p,order,droots[num-1],upper,&root) == TRUE) {
  1060.                 roots[found] = root;
  1061.                 found++;
  1062.             }
  1063.             else if (findit(p,order,droots[num-1],upper,&root) == TRUE) {
  1064.                 roots[found] = root;
  1065.                 found++;
  1066.             }
  1067.         }
  1068.     }
  1069.  
  1070.     free(droots);
  1071.     return found;
  1072.             
  1073. }
  1074.