home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / sphsweep.c < prev    next >
C/C++ Source or Header  |  2000-11-26  |  48KB  |  1,943 lines
  1. /****************************************************************************
  2. *                sphsweep.c
  3. *
  4. *  This module implements the sphere sweep primitive. Idea: J.J. Van Wijk,
  5. *  "Raytracing Objects Defined by Sweeping a Sphere", Eurographics 1984.
  6. *
  7. *  Copyright 1997-98 Jochen Lippert
  8. *
  9. *
  10. *
  11. *  Version history:
  12. *
  13. *  24. Feb. 1998: Fixed a bug in Compute_Sphere_Sweep_BBox that could
  14. *                 result in much too large bounding boxes for
  15. *                 Catmull-Rom-Spline Sphere Sweeps. Added statistics
  16. *                 output for Sphere Sweep primitive.
  17. *
  18. *  19. Jan. 1998: Corrected a problem with the calculation of single
  19. *                 spheres in Compute_Sphere_Sweep which I introduced
  20. *                 in the previous version, and cleaned up the code in
  21. *                 this function a bit. All_Sphere_Sweep_Intersections
  22. *                 now uses QSORT to sort the intersection list.
  23. *
  24. *  6.  Jan. 1998: Fixed a bug that affected the geometry of
  25. *                 Catmull-Rom-Spline Sphere Sweeps, fixed a bug in
  26. *                 Inside_Sphere_Sweep which affected non-proportionally
  27. *                 scaled Sphere Sweeps.
  28. *
  29. *  22. Dec. 1997: Simplyfied code for finding invalid intersections,
  30. *                 fixed some bugs related to this.
  31. *
  32. *  8.  Dec. 1997: Some cleanup & bug fixes.
  33. *
  34. *  1.  Dec. 1997: Made depth tolerance user-selectable.
  35. *
  36. *  29. Nov. 1997: Added support for Catmull-Rom-Splines.
  37. *
  38. *  26. Nov. 1997: Fixed a bug where the wrong value was returned by
  39. *                 All_Sphere_Sweep_Intersections when the clipped_by
  40. *                 statement was used with a sphere sweep.
  41. *
  42. *  24. Nov. 1997: Fixed a pretty silly bug in Copy_Sphere_Sweep where
  43. *                 the sphere sweep wasn't copied completely.
  44. *
  45. *  21. Nov. 1997: Initial release.
  46. *
  47. *****************************************************************************/
  48.  
  49. #include "frame.h"
  50. #include "povray.h"
  51. #include "vector.h"
  52. #include "povproto.h"
  53. #include "bbox.h"
  54. #include "matrices.h"
  55. #include "objects.h"
  56. #include "polysolv.h"
  57. #include "sphsweep.h"
  58.  
  59.  
  60.  
  61. /*****************************************************************************
  62. * Local preprocessor defines
  63. ******************************************************************************/
  64.  
  65. #define DEPTH_TOLERANCE        1.0e-6
  66. #define ZERO_TOLERANCE        1.0e-1
  67. #define SPHSWEEP_MAX_ISECT    64
  68.  
  69.  
  70.  
  71. /*****************************************************************************
  72. * Local typedefs
  73. ******************************************************************************/
  74.  
  75. typedef struct Sphere_Sweep_Intersection_Structure SPHSWEEP_INT;
  76.  
  77. /* Temporary storage for intersection values */
  78. struct Sphere_Sweep_Intersection_Structure
  79. {
  80.     DBL        t;            /* Distance along ray           */
  81.     VECTOR    Point;        /* Intersection point           */
  82.     VECTOR    Normal;        /* Normal at intersection point */
  83. };
  84.  
  85.  
  86.  
  87. /*****************************************************************************
  88. * Static functions
  89. ******************************************************************************/
  90. static int All_Sphere_Sweep_Intersections (OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack);
  91. static int Inside_Sphere_Sweep (VECTOR IPoint, OBJECT *Object);
  92. static void Sphere_Sweep_Normal (VECTOR Result, OBJECT *Object, INTERSECTION *Inter);
  93. static void Translate_Sphere_Sweep (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
  94. static void Rotate_Sphere_Sweep (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
  95. static void Scale_Sphere_Sweep (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
  96. static void Invert_Sphere_Sweep (OBJECT *Object);
  97.  
  98.  
  99.  
  100. /*****************************************************************************
  101. * Local functions
  102. ******************************************************************************/
  103. int Intersect_Sphere_Sweep_Sphere (RAY *Ray, SPHSWEEP_SPH *Sphere, SPHSWEEP_INT *Isect);
  104. int Intersect_Sphere_Sweep_Segment (RAY *Ray, SPHSWEEP_SEG *Segment, SPHSWEEP_INT *Isect);
  105. int Find_Valid_Points (SPHSWEEP_INT *Inter, int Num_Inter, RAY *Ray);
  106. static int Comp_Isects (CONST void *Intersection_1, CONST void *Intersection_2);
  107.  
  108.  
  109.  
  110. /*****************************************************************************
  111. * Local variables
  112. ******************************************************************************/
  113.  
  114. static METHODS Sphere_Sweep_Methods =
  115. {
  116.   All_Sphere_Sweep_Intersections,
  117.   Inside_Sphere_Sweep, Sphere_Sweep_Normal, Default_UVCoord,
  118.   Copy_Sphere_Sweep,
  119.   Translate_Sphere_Sweep, Rotate_Sphere_Sweep,
  120.   Scale_Sphere_Sweep, Transform_Sphere_Sweep, Invert_Sphere_Sweep,
  121.   Destroy_Sphere_Sweep
  122. };
  123.  
  124. MATRIX Catmull_Rom_Matrix =
  125. {
  126.      0.0 / 2.0,  2.0 / 2.0,  0.0 / 2.0,  0.0 / 2.0,
  127.     -1.0 / 2.0,  0.0 / 2.0,  1.0 / 2.0,  0.0 / 2.0,
  128.      2.0 / 2.0, -5.0 / 2.0,  4.0 / 2.0, -1.0 / 2.0,
  129.     -1.0 / 2.0,  3.0 / 2.0, -3.0 / 2.0,  1.0 / 2.0
  130. };
  131.  
  132. MATRIX B_Matrix =
  133. {
  134.      1.0 / 6.0,  4.0 / 6.0,  1.0 / 6.0,  0.0 / 6.0,
  135.     -3.0 / 6.0,  0.0 / 6.0,  3.0 / 6.0,  0.0 / 6.0,
  136.      3.0 / 6.0, -6.0 / 6.0,  3.0 / 6.0,  0.0 / 6.0,
  137.     -1.0 / 6.0,  3.0 / 6.0, -3.0 / 6.0,  1.0 / 6.0
  138. };
  139.  
  140.  
  141.  
  142. /*****************************************************************************
  143. *
  144. * FUNCTION
  145. *
  146. *   All_Sphere_Sweep_Intersections
  147. *
  148. * INPUT
  149. *
  150. *   Object, Ray, Depth stack
  151. *
  152. * OUTPUT
  153. *
  154. *   Depth stack
  155. *
  156. * RETURNS
  157. *
  158. *   Boolean - found any intersections?
  159. *
  160. * AUTHOR
  161. *
  162. *   Jochen Lippert
  163. *
  164. * DESCRIPTION
  165. *
  166. *   Find all valid intersections of a ray with a sphere sweep.
  167. *
  168. * CHANGES
  169. *
  170. *   -
  171. *
  172. ******************************************************************************/
  173.  
  174. static int All_Sphere_Sweep_Intersections(OBJECT    *Object, RAY *Ray,ISTACK *Depth_Stack)
  175. {
  176.     SPHERE_SWEEP    *Sphere_Sweep = (SPHERE_SWEEP *)Object;
  177.     RAY                New_Ray;
  178.     DBL                len;
  179.     int                Intersection_Found;
  180.     SPHSWEEP_INT    Isect[SPHSWEEP_MAX_ISECT];
  181.     int                Num_Isect;
  182.     SPHSWEEP_INT    Sphere_Isect[2];
  183.     SPHSWEEP_INT    Segment_Isect[12];
  184.     int                Num_Seg_Isect;
  185.     int                i, j;
  186.     
  187.     Increase_Counter(stats[Ray_Sphere_Sweep_Tests]);
  188.     Intersection_Found = FALSE;
  189.     
  190.     if (Sphere_Sweep->Trans == NULL)
  191.     {
  192.         Assign_Vector(New_Ray.Initial, Ray->Initial);
  193.         Assign_Vector(New_Ray.Direction, Ray->Direction);
  194.     }
  195.     else
  196.     {
  197.         MInvTransPoint(New_Ray.Initial, Ray->Initial, Sphere_Sweep->Trans);
  198.         MInvTransDirection(New_Ray.Direction, Ray->Direction, Sphere_Sweep->Trans);
  199.         
  200.         VLength(len, New_Ray.Direction);
  201.         VInverseScaleEq(New_Ray.Direction, len);
  202.     }
  203.     
  204.     Num_Isect = 0;
  205.     
  206.     /* Intersections with single spheres */
  207.     for (i = 0; i < Sphere_Sweep->Num_Spheres; i++)
  208.     {
  209.         /* Are there intersections with this sphere? */
  210.         if (Intersect_Sphere_Sweep_Sphere(&New_Ray, &Sphere_Sweep->Sphere[i], Sphere_Isect))
  211.         {
  212.             /* Test for end of vector */
  213.             if (Num_Isect + 2 <= SPHSWEEP_MAX_ISECT)
  214.             {
  215.             /* Valid intersection at Sphere_Isect[0]? */
  216.             if ((Sphere_Isect[0].t > -Max_Distance) && (Sphere_Isect[0].t < Max_Distance))
  217.             {
  218.                 /* Add intersection */
  219.                 Isect[Num_Isect] = Sphere_Isect[0];
  220.                 Num_Isect++;
  221.             }
  222.             
  223.             /* Valid intersection at Sphere_Isect[1]? */
  224.             if ((Sphere_Isect[1].t > -Max_Distance) && (Sphere_Isect[1].t < Max_Distance))
  225.             {
  226.                 /* Add intersection */
  227.                 Isect[Num_Isect] = Sphere_Isect[1];
  228.                 Num_Isect++;
  229.             }
  230.         }
  231.     }
  232.     }
  233.     
  234.     /* Intersections with segments */
  235.     for (i = 0; i < Sphere_Sweep->Num_Segments; i++)
  236.     {
  237.         /* Are there intersections with this segment? */
  238.         Num_Seg_Isect = Intersect_Sphere_Sweep_Segment(&New_Ray, &Sphere_Sweep->Segment[i], Segment_Isect);
  239.         
  240.         /* Test for end of vector */
  241.         if (Num_Isect + Num_Seg_Isect <= SPHSWEEP_MAX_ISECT)
  242.         {
  243.             for (j = 0; j < Num_Seg_Isect; j++)
  244.             {
  245.                 /* Add intersection */
  246.             Isect[Num_Isect] = Segment_Isect[j];
  247.             Num_Isect++;
  248.         }
  249.     }
  250.     }
  251.     
  252.     /* Any intersections so far? */
  253.     if (Num_Isect == 0)
  254.     {
  255.         /* No, return */
  256.         return (FALSE);
  257.     }
  258.     
  259.     /* Sort intersections */
  260.     QSORT((void *)Isect, (size_t)Num_Isect, sizeof(SPHSWEEP_INT), Comp_Isects);
  261.     
  262.     /* Delete invalid intersections inside the sphere sweep */
  263.     Num_Isect = Find_Valid_Points(Isect, Num_Isect, &New_Ray);
  264.     
  265.     /* Push valid intersections */
  266.     for (i = 0; i < Num_Isect; i++)
  267.     {
  268.         /* Was the ray transformed? */
  269.         if (Sphere_Sweep->Trans != NULL)
  270.         {
  271.             /* Yes, invert the transformation */
  272.             Isect[i].t /= len;
  273.             MTransPoint(Isect[i].Point, Isect[i].Point, Sphere_Sweep->Trans);
  274.             MTransNormal(Isect[i].Normal, Isect[i].Normal, Sphere_Sweep->Trans);
  275.             VNormalizeEq(Isect[i].Normal);
  276.         }
  277.         
  278.         /* Check for positive values of t (it's a ray after all...) */
  279.         if (Isect[i].t > Sphere_Sweep->Depth_Tolerance)
  280.         {
  281.             /* Test for clipping volume */
  282.             if (Point_In_Clip(Isect[i].Point, Object->Clip))
  283.             {
  284.                 push_normal_entry(Isect[i].t,
  285.                                   Isect[i].Point,
  286.                                   Isect[i].Normal, Object, Depth_Stack);
  287.                 Intersection_Found = TRUE;
  288.             }
  289.         }
  290.     }
  291.     
  292.     if (Intersection_Found == TRUE)
  293.     {
  294.         Increase_Counter(stats[Ray_Sphere_Sweep_Tests_Succeeded]);
  295.     }
  296.     
  297.     return (Intersection_Found);
  298. }
  299.  
  300.  
  301.  
  302. /*****************************************************************************
  303. *
  304. * FUNCTION
  305. *
  306. *   Intersect_Sphere_Sweep_Sphere
  307. *
  308. * INPUT
  309. *
  310. *   Ray     - Ray to test intersection with
  311. *   Sphere  - Sphere
  312. *   Inter   - Intersections (two element vector)
  313. *   
  314. * OUTPUT
  315. *
  316. *   -
  317. *
  318. * RETURNS
  319. *
  320. *   Boolean - is there at least one intersection?
  321. *
  322. * AUTHOR
  323. *
  324. *   Jochen Lippert
  325. *   
  326. * DESCRIPTION
  327. *
  328. *   Find intersections of ray (line in fact) with a sphere.
  329. *
  330. * CHANGES
  331. *
  332. *   -
  333. *
  334. ******************************************************************************/
  335.  
  336. int Intersect_Sphere_Sweep_Sphere(RAY *Ray,SPHSWEEP_SPH    *Sphere,SPHSWEEP_INT    *Inter)
  337. {
  338.     DBL        Radius2;
  339.     DBL        OCSquared;
  340.     DBL        t_Closest_Approach;
  341.     DBL        Half_Chord;
  342.     DBL        t_Half_Chord_Squared;
  343.     VECTOR    Origin_To_Center;
  344.     
  345.     Radius2 = Sqr(Sphere->Radius);
  346.     
  347.     VSub(Origin_To_Center, Sphere->Center, Ray->Initial);
  348.     
  349.     VDot(OCSquared, Origin_To_Center, Origin_To_Center);
  350.     
  351.     VDot(t_Closest_Approach, Origin_To_Center, Ray->Direction);
  352.     
  353.     if ((OCSquared >= Radius2) && (t_Closest_Approach < EPSILON))
  354.     {
  355.         return (FALSE);
  356.     }
  357.     
  358.     t_Half_Chord_Squared = Radius2 - OCSquared + Sqr(t_Closest_Approach);
  359.     
  360.     if (t_Half_Chord_Squared > EPSILON)
  361.     {
  362.         Half_Chord = sqrt(t_Half_Chord_Squared);
  363.         
  364.         /* Calculate smaller depth */
  365.         Inter[0].t = t_Closest_Approach - Half_Chord;
  366.         
  367.         /* Calculate point */
  368.         VEvaluateRay(Inter[0].Point, Ray->Initial, Inter[0].t, Ray->Direction);
  369.         
  370.         /* Calculate normal */
  371.         VSub(Inter[0].Normal, Inter[0].Point, Sphere->Center);
  372.         VInverseScaleEq(Inter[0].Normal, Sphere->Radius);
  373.         
  374.         /* Calculate bigger depth */
  375.         Inter[1].t = t_Closest_Approach + Half_Chord;
  376.         
  377.         /* Calculate point */
  378.         VEvaluateRay(Inter[1].Point, Ray->Initial, Inter[1].t, Ray->Direction);
  379.         
  380.         /* Calculate normal */
  381.         VSub(Inter[1].Normal, Inter[1].Point, Sphere->Center);
  382.         VInverseScaleEq(Inter[1].Normal, Sphere->Radius);
  383.         
  384.         /* Increase_Counter(stats[Ray_Sphere_Tests_Succeeded]); */
  385.         
  386.         return (TRUE);
  387.     }
  388.     return (FALSE);
  389. }
  390.  
  391.  
  392.  
  393. /*****************************************************************************
  394. *
  395. * FUNCTION
  396. *
  397. *   Intersect_Sphere_Sweep_Segment
  398. *
  399. * INPUT
  400. *
  401. *   Ray     - Ray to test intersection with
  402. *   Segment - Segment of a sphere sweep
  403. *   Isect   - intersection list (depth, point, normal)
  404. *
  405. * OUTPUT
  406. *
  407. *   Isect   - intersection list (depth, point, normal)
  408. *
  409. * RETURNS
  410. *
  411. *   Number of intersections
  412. *
  413. * AUTHOR
  414. *
  415. *   Jochen Lippert
  416. *   
  417. * DESCRIPTION
  418. *
  419. *   Find all intersections of a ray (line in fact) with one segment
  420. *   of a sphere sweep.
  421. *
  422. * CHANGES
  423. *
  424. *   -
  425. *
  426. ******************************************************************************/
  427.  
  428. int Intersect_Sphere_Sweep_Segment(RAY *Ray,SPHSWEEP_SEG *Segment,SPHSWEEP_INT *Isect)
  429. {
  430.     int        Isect_Count;
  431.     DBL        Dot1, Dot2;
  432.     DBL        t1, t2;
  433.     VECTOR    Vector;
  434.     VECTOR    IPoint;
  435.     VECTOR    DCenter;
  436.     DBL        DCenterDot;
  437.     DBL        temp;
  438.     DBL        b, c, d, e, f, g, h, i, j, k, l;
  439.     DBL        Coef[11];
  440.     DBL        Root[10];
  441.     int        Num_Poly_Roots, m, n;
  442.     DBL        fp0, fp1;
  443.     DBL        t;
  444.     SPHSWEEP_SPH    Temp_Sphere;
  445.     SPHSWEEP_INT    Temp_Isect[2];
  446.     VECTOR    Center;
  447.     
  448.     Isect_Count = 0;
  449.     
  450.     /* Calculate intersections with closing surface for u = 0 */
  451.     VDot(Dot1, Ray->Direction, Segment->Center_Deriv[0]);
  452.     if (fabs(Dot1) > EPSILON)
  453.     {
  454.         VSub(Vector, Ray->Initial, Segment->Closing_Sphere[0].Center);
  455.         VDot(Dot2, Vector, Segment->Center_Deriv[0]);
  456.         t1 = -(Dot2 + Segment->Closing_Sphere[0].Radius * Segment->Radius_Deriv[0]) / Dot1;
  457.         if ((t1 > -Max_Distance) && (t1 < Max_Distance))
  458.         {
  459.             /* Calculate point */
  460.             VEvaluateRay(IPoint, Ray->Initial, t1, Ray->Direction);
  461.             
  462.             /* Is the point inside the closing sphere? */
  463.             VSub(DCenter, IPoint, Segment->Closing_Sphere[0].Center);
  464.             VDot(DCenterDot, DCenter, DCenter);
  465.             if (DCenterDot < Sqr(Segment->Closing_Sphere[0].Radius))
  466.             {
  467.                 /* Add intersection */
  468.                 Isect[Isect_Count].t = t1;
  469.                 
  470.                 /* Copy point */
  471.                 Assign_Vector(Isect[Isect_Count].Point, IPoint);
  472.                 
  473.                 /* Calculate normal */
  474.                 Assign_Vector(Isect[Isect_Count].Normal, Segment->Center_Deriv[0]);
  475.                 VScaleEq(Isect[Isect_Count].Normal, -1.0);
  476.                 VNormalizeEq(Isect[Isect_Count].Normal);
  477.                 
  478.                 Isect_Count++;
  479.             }
  480.         }
  481.     }
  482.     
  483.     /* Calculate intersections with closing surface for u = 1 */
  484.     VDot(Dot1, Ray->Direction, Segment->Center_Deriv[1]);
  485.     if (fabs(Dot1) > EPSILON)
  486.     {
  487.         VSub(Vector, Ray->Initial, Segment->Closing_Sphere[1].Center);
  488.         VDot(Dot2, Vector, Segment->Center_Deriv[1]);
  489.         t2 = -(Dot2 + Segment->Closing_Sphere[1].Radius * Segment->Radius_Deriv[1]) / Dot1;
  490.         if ((t2 > -Max_Distance) && (t2 < Max_Distance))
  491.         {
  492.             /* Calculate point */
  493.             VEvaluateRay(IPoint, Ray->Initial, t2, Ray->Direction);
  494.             
  495.             /* Is the point inside the closing sphere? */
  496.             VSub(DCenter, IPoint, Segment->Closing_Sphere[1].Center);
  497.             VDot(DCenterDot, DCenter, DCenter);
  498.             if (DCenterDot < Sqr(Segment->Closing_Sphere[1].Radius))
  499.             {
  500.                 /* Add intersection */
  501.                 Isect[Isect_Count].t = t2;
  502.                 
  503.                 /* Copy point */
  504.                 Assign_Vector(Isect[Isect_Count].Point, IPoint);
  505.                 
  506.                 /* Calculate normal */
  507.                 Assign_Vector(Isect[Isect_Count].Normal, Segment->Center_Deriv[1]);
  508.                 VNormalizeEq(Isect[Isect_Count].Normal);
  509.                 
  510.                 Isect_Count++;
  511.             }
  512.         }
  513.     }
  514.     
  515.     /* Calculate intersections with sides of the segment */
  516.     
  517.     switch (Segment->Num_Coefs)
  518.     {
  519.         case 2:        /* First order Polynomial */
  520.             
  521.             VSub(Vector, Ray->Initial, Segment->Center_Coef[0]);
  522.             
  523.             /* a is always 1.0 */
  524.             
  525.             VDot(b, Segment->Center_Coef[1], Ray->Direction);
  526.             b *= -2.0;
  527.             
  528.             VDot(c, Vector, Ray->Direction);
  529.             c *= 2.0;
  530.             
  531.             VDot(d, Segment->Center_Coef[1], Segment->Center_Coef[1]);
  532.             d -= Sqr(Segment->Radius_Coef[1]);
  533.             
  534.             VDot(e, Vector, Segment->Center_Coef[1]);
  535.             e += Segment->Radius_Coef[0] * Segment->Radius_Coef[1];
  536.             e *= -2.0;
  537.             
  538.             VDot(f, Vector, Vector);
  539.             f -= Sqr(Segment->Radius_Coef[0]);
  540.             
  541.             Coef[0] = 4.0 * Sqr(d) - Sqr(b) * d;
  542.             Coef[1] = 4.0 * d * e - 2.0 * b * c * d;
  543.             Coef[2] = Sqr(e) - b * c * e + Sqr(b) * f;
  544.             
  545.             Num_Poly_Roots = Solve_Polynomial(2, Coef, Root, TRUE, 1e-10);
  546.             break;
  547.         
  548.         case 4:        /* Third order polynomial */
  549.             
  550.             VSub(Vector, Ray->Initial, Segment->Center_Coef[0]);
  551.             
  552.             /* a is always 1.0 */
  553.             
  554.             VDot(b, Segment->Center_Coef[3], Ray->Direction);
  555.             b *= -2.0;
  556.             
  557.             VDot(c, Segment->Center_Coef[2], Ray->Direction);
  558.             c *= -2.0;
  559.             
  560.             VDot(d, Segment->Center_Coef[1], Ray->Direction);
  561.             d *= -2.0;
  562.             
  563.             VDot(e, Vector, Ray->Direction);
  564.             e *= 2.0;
  565.             
  566.             VDot(f, Segment->Center_Coef[3], Segment->Center_Coef[3]);
  567.             f -= Sqr(Segment->Radius_Coef[3]);
  568.             
  569.             VDot(g, Segment->Center_Coef[3], Segment->Center_Coef[2]);
  570.             g -= Segment->Radius_Coef[3] * Segment->Radius_Coef[2];
  571.             g *= 2.0;
  572.             
  573.             VDot(h, Segment->Center_Coef[3], Segment->Center_Coef[1]);
  574.             h *= 2.0;
  575.             VDot(temp, Segment->Center_Coef[2], Segment->Center_Coef[2]);
  576.             h += temp;
  577.             h -= 2.0 * Segment->Radius_Coef[3] * Segment->Radius_Coef[1];
  578.             h -= Sqr(Segment->Radius_Coef[2]);
  579.             
  580.             VDot(i, Segment->Center_Coef[3], Vector);
  581.             VDot(temp, Segment->Center_Coef[2], Segment->Center_Coef[1]);
  582.             i -= temp;
  583.             i += Segment->Radius_Coef[3] * Segment->Radius_Coef[0];
  584.             i += Segment->Radius_Coef[2] * Segment->Radius_Coef[1];
  585.             i *= -2.0;
  586.             
  587.             VDot(j, Segment->Center_Coef[2], Vector);
  588.             j += Segment->Radius_Coef[2] * Segment->Radius_Coef[0];
  589.             j *= -2.0;
  590.             VDot(temp, Segment->Center_Coef[1], Segment->Center_Coef[1]);
  591.             j += temp;
  592.             j -= Sqr(Segment->Radius_Coef[1]);
  593.             
  594.             VDot(k, Segment->Center_Coef[1], Vector);
  595.             k += Segment->Radius_Coef[1] * Segment->Radius_Coef[0];
  596.             k *= -2.0;
  597.             
  598.             VDot(l, Vector, Vector);
  599.             l -= Sqr(Segment->Radius_Coef[0]);
  600.             
  601.             Coef[0] = 36.0 * Sqr(f) - 9.0 * f * Sqr(b);
  602.             Coef[1] = 60.0 * f * g - 6.0 * g * Sqr(b) - 18.0 * b * c * f;
  603.             Coef[2] = 48.0 * f * h + 25.0 * Sqr(g) - 3.0 * h * Sqr(b)
  604.                     - 13.0 * b * c * g - 8.0 * f * Sqr(c) - 18.0 * b * d * f;
  605.             Coef[3] = 36.0 * f * i + 40.0 * g * h - 18.0 * b * f * e - 8.0 * b * c * h
  606.                     - 6.0 * g * Sqr(c) - 14.0 * b * d * g - 14.0 * c * d * f;
  607.             Coef[4] = 24.0 * f * j + 30.0 * g * i + 16.0 * Sqr(h) - 15.0 * b * g * e
  608.                     - 12.0 * c * f * e + 3.0 * j * Sqr(b) - 3.0 * b * c * i - 4.0 * h * Sqr(c)
  609.                     - 10.0 * b * d * h - 11.0 * c * d * g - 5.0 * f * Sqr(d);
  610.             Coef[5] = 12.0 * f * k + 20.0 * g * j + 24.0 * h * i - 12.0 * b * h * e
  611.                     - 10.0 * c * g * e - 6.0 * d * f * e + 6.0 * k * Sqr(b) + 2.0 * b * c * j
  612.                     - 2.0 * i * Sqr(c) - 6.0 * b * d * i - 8.0 * c * d * h - 4.0 * g * Sqr(d);
  613.             Coef[6] = 10.0 * g * k + 16.0 * h * j + 9.0 * Sqr(i) - 9.0 * b * i * e
  614.                     - 8.0 * c * h * e - 5.0 * d * g * e + 9.0 * l * Sqr(b) + 7.0 * b * c * k
  615.                     - 2.0 * b * d * j - 5.0 * c * d * i - 3.0 * h * Sqr(d);
  616.             Coef[7] = 8.0 * h * k + 12.0 * i * j - 6.0 * b * j * e - 6.0 * c * i * e
  617.                     - 4.0 * d * h * e + 12.0 * b * c * l + 2.0 * k * Sqr(c) + 2.0 * b * d * k
  618.                     - 2.0 * c * d * j - 2.0 * i * Sqr(d);
  619.             Coef[8] = 6.0 * i * k + 4.0 * Sqr(j) - 3.0 * b * k * e - 4.0 * c * j * e
  620.                     - 3.0 * d * i * e + 4.0 * l * Sqr(c) + 6.0 * b * d * l
  621.                     + c * d * k - j * Sqr(d);
  622.             Coef[9] = 4.0 * j * k - 2.0 * c * k * e - 2.0 * d * j * e + 4.0 * c * d * l;
  623.             Coef[10] = Sqr(k) - d * k * e + l * Sqr(d);
  624.             
  625.             Num_Poly_Roots = Solve_Polynomial(10, Coef, Root, TRUE, 1e-10);
  626.             break;
  627.     }
  628.     
  629.     /* Remove roots not in interval [0, 1] */
  630.     
  631.     m = 0;
  632.     while (m < Num_Poly_Roots)
  633.     {
  634.         if (Root[m] < 0.0 || Root[m] > 1.0)
  635.         {
  636.             for (n = m; n < Num_Poly_Roots - 1; n++)
  637.             {
  638.                 Root[n] = Root[n + 1];
  639.             }
  640.             Num_Poly_Roots--;
  641.         }
  642.         else
  643.         {
  644.             m++;
  645.         }
  646.     }
  647.     
  648.     switch (Segment->Num_Coefs)
  649.     {
  650.         case 2:
  651.             
  652.             for (m = 0; m < Num_Poly_Roots; m++)
  653.             {
  654.                 fp0 = 2.0 * d * Root[m]
  655.                     + e;
  656.                 fp1 = b;
  657.                 
  658.                 if (fabs(fp1) > ZERO_TOLERANCE)
  659.                 {
  660.                     t = -fp0 / fp1;
  661.                     
  662.                     if ((t > -Max_Distance) && (t < Max_Distance))
  663.                     {
  664.                         /* Add intersection */
  665.                         Isect[Isect_Count].t = t;
  666.                         
  667.                         /* Calculate point */
  668.                         VEvaluateRay(Isect[Isect_Count].Point, Ray->Initial, t, Ray->Direction);
  669.                         
  670.                         /* Calculate normal */
  671.                         VAddScaled(Center, Segment->Center_Coef[0], Root[m], Segment->Center_Coef[1]);
  672.                         VSub(Isect[Isect_Count].Normal, Isect[Isect_Count].Point, Center);
  673.                         VNormalizeEq(Isect[Isect_Count].Normal);
  674.                         
  675.                         Isect_Count++;
  676.                     }
  677.                 }
  678.                 else
  679.                 {
  680.                     /* Calculate center of single sphere */
  681.                     VAddScaled(Temp_Sphere.Center, Segment->Center_Coef[0], Root[m], Segment->Center_Coef[1]);
  682.                     
  683.                     /* Calculate radius of single sphere */
  684.                     Temp_Sphere.Radius = Segment->Radius_Coef[1] * Root[m] + Segment->Radius_Coef[0];
  685.                     
  686.                     /* Calculate intersections */
  687.                     if (Intersect_Sphere_Sweep_Sphere(Ray, &Temp_Sphere, Temp_Isect))
  688.                     {
  689.                         /* Add intersections */
  690.                         Isect[Isect_Count] = Temp_Isect[0];
  691.                         Isect_Count++;
  692.                         
  693.                         Isect[Isect_Count] = Temp_Isect[1];
  694.                         Isect_Count++;
  695.                     }
  696.                 }
  697.             }
  698.             break;
  699.         
  700.         case 4:
  701.             
  702.             for (m = 0; m < Num_Poly_Roots; m++)
  703.             {
  704.                 fp0 = 6.0 * f * Root[m] * Root[m] * Root[m] * Root[m] * Root[m]
  705.                     + 5.0 * g * Root[m] * Root[m] * Root[m] * Root[m]
  706.                     + 4.0 * h * Root[m] * Root[m] * Root[m]
  707.                     + 3.0 * i * Root[m] * Root[m]
  708.                     + 2.0 * j * Root[m]
  709.                     + k;
  710.                 fp1 = 3.0 * b * Root[m] * Root[m]
  711.                     + 2.0 * c * Root[m]
  712.                     + d;
  713.                 
  714.                 if (fabs(fp1) > ZERO_TOLERANCE)
  715.                 {
  716.                     t = -fp0 / fp1;
  717.                     
  718.                     if ((t > -Max_Distance) && (t < Max_Distance))
  719.                     {
  720.                         /* Add intersection */
  721.                         Isect[Isect_Count].t = t;
  722.                         
  723.                         /* Calculate point */
  724.                         VEvaluateRay(Isect[Isect_Count].Point, Ray->Initial, t, Ray->Direction);
  725.                         
  726.                         /* Calculate normal */
  727.                         VAddScaled(Center, Segment->Center_Coef[0], Root[m], Segment->Center_Coef[1]);
  728.                         VAddScaledEq(Center, Root[m] * Root[m], Segment->Center_Coef[2]);
  729.                         VAddScaledEq(Center, Root[m] * Root[m] * Root[m], Segment->Center_Coef[3]);
  730.                         VSub(Isect[Isect_Count].Normal, Isect[Isect_Count].Point, Center);
  731.                         VNormalizeEq(Isect[Isect_Count].Normal);
  732.                         
  733.                         Isect_Count++;
  734.                     }
  735.                 }
  736.                 else
  737.                 {
  738.                     /* Calculate center of single sphere */
  739.                     VAddScaled(Temp_Sphere.Center, Segment->Center_Coef[0], Root[m], Segment->Center_Coef[1]);
  740.                     VAddScaledEq(Temp_Sphere.Center, Root[m] * Root[m], Segment->Center_Coef[2]);
  741.                     VAddScaledEq(Temp_Sphere.Center, Root[m] * Root[m] * Root[m], Segment->Center_Coef[3]);
  742.                     
  743.                     /* Calculate radius of single sphere */
  744.                     Temp_Sphere.Radius = Segment->Radius_Coef[3] * Root[m] * Root[m] * Root[m]
  745.                                        + Segment->Radius_Coef[2] * Root[m] * Root[m]
  746.                                        + Segment->Radius_Coef[1] * Root[m]
  747.                                        + Segment->Radius_Coef[0];
  748.                     
  749.                     /* Calculate intersections */
  750.                     if (Intersect_Sphere_Sweep_Sphere(Ray, &Temp_Sphere, Temp_Isect))
  751.                     {
  752.                         /* Add intersections */
  753.                         Isect[Isect_Count] = Temp_Isect[0];
  754.                         Isect_Count++;
  755.                         
  756.                         Isect[Isect_Count] = Temp_Isect[1];
  757.                         Isect_Count++;
  758.                     }
  759.                 }
  760.             }
  761.             break;
  762.     }
  763.     
  764.     return (Isect_Count);
  765. }
  766.  
  767.  
  768.  
  769. /*****************************************************************************
  770. *
  771. * FUNCTION
  772. *
  773. *   Inside_Sphere_Sweep
  774. *
  775. * INPUT
  776. *
  777. *   Point, Object
  778. *
  779. * OUTPUT
  780. *
  781. *   -
  782. *
  783. * RETURNS
  784. *
  785. *   Boolean - is the point inside the sphere sweep?
  786. *
  787. * AUTHOR
  788. *
  789. *   Jochen Lippert
  790. *
  791. * DESCRIPTION
  792. *
  793. *   Test if point is inside sphere sweep.
  794. *
  795. * CHANGES
  796. *
  797. *   -
  798. *
  799. ******************************************************************************/
  800.  
  801. static int Inside_Sphere_Sweep(VECTOR IPoint,OBJECT *Object)
  802. {
  803.     SPHERE_SWEEP *Sphere_Sweep = (SPHERE_SWEEP *)Object;
  804.     int        inside;
  805.     VECTOR    New_Point;
  806.     int        i, j;
  807.     VECTOR    Vector;
  808.     DBL        temp;
  809.     DBL        Coef[7];
  810.     DBL        Root[6];
  811.     int        Num_Poly_Roots;
  812.     
  813.     inside = FALSE;
  814.     
  815.     if (Sphere_Sweep->Trans == NULL)
  816.     {
  817.         Assign_Vector(New_Point, IPoint);
  818.     }
  819.     else
  820.     {
  821.         MInvTransPoint(New_Point, IPoint, Sphere_Sweep->Trans);
  822.     }
  823.     
  824.     switch (Sphere_Sweep->Interpolation)
  825.     {
  826.         case LINEAR_SPHERE_SWEEP:
  827.             
  828.             /* For each segment... */
  829.             for (i = 0; i < Sphere_Sweep->Num_Segments; i++)
  830.             {
  831.                 /* Pre-calculate vector */
  832.                 VSub(Vector, New_Point, Sphere_Sweep->Segment[i].Center_Coef[0]);
  833.                 
  834.                 /* Coefficient for u^2 */
  835.                 VDot(Coef[0], Sphere_Sweep->Segment[i].Center_Coef[1],
  836.                               Sphere_Sweep->Segment[i].Center_Coef[1]);
  837.                 Coef[0] -= Sqr(Sphere_Sweep->Segment[i].Radius_Coef[1]);
  838.                 
  839.                 /* Coefficient for u^1 */
  840.                 VDot(Coef[1], Vector, Sphere_Sweep->Segment[i].Center_Coef[1]);
  841.                 Coef[1] += Sphere_Sweep->Segment[i].Radius_Coef[0]
  842.                          * Sphere_Sweep->Segment[i].Radius_Coef[1];
  843.                 Coef[1] *= -2.0;
  844.                 
  845.                 /* Coefficient for u^0 */
  846.                 VDot(Coef[2], Vector, Vector);
  847.                 Coef[2] -= Sqr(Sphere_Sweep->Segment[i].Radius_Coef[0]);
  848.                 
  849.                 /* Find roots */
  850.                 Num_Poly_Roots = Solve_Polynomial(2, Coef, Root, TRUE, 1e-10);
  851.                 
  852.                 /* Test for interval [0, 1] */
  853.                 for (j = 0; j < Num_Poly_Roots; j++)
  854.                 {
  855.                     if (Root[j] >= 0.0 && Root[j] <= 1.0)
  856.                     {
  857.                         /* At least one root inside interval, */
  858.                         /* so we are inside sphere sweep      */
  859.                         inside = TRUE;
  860.                         break;
  861.                     }
  862.                 }
  863.             }
  864.             break;
  865.         
  866.         case CATMULL_ROM_SPLINE_SPHERE_SWEEP:
  867.         case B_SPLINE_SPHERE_SWEEP:
  868.             
  869.             /* For each segment... */
  870.             for (i = 0; i < Sphere_Sweep->Num_Segments; i++)
  871.             {
  872.                 /* Pre-calculate vector */
  873.                 VSub(Vector, New_Point, Sphere_Sweep->Segment[i].Center_Coef[0]);
  874.                 
  875.                 /* Coefficient for u^6 */
  876.                 VDot(Coef[0], Sphere_Sweep->Segment[i].Center_Coef[3],
  877.                               Sphere_Sweep->Segment[i].Center_Coef[3]);
  878.                 Coef[0] -= Sqr(Sphere_Sweep->Segment[i].Radius_Coef[3]);
  879.                 
  880.                 /* Coefficient for u^5 */
  881.                 VDot(Coef[1], Sphere_Sweep->Segment[i].Center_Coef[3],
  882.                               Sphere_Sweep->Segment[i].Center_Coef[2]);
  883.                 Coef[1] -= Sphere_Sweep->Segment[i].Radius_Coef[3]
  884.                          * Sphere_Sweep->Segment[i].Radius_Coef[2];
  885.                 Coef[1] *= 2.0;
  886.                 
  887.                 /* Coefficient for u^4 */
  888.                 VDot(Coef[2], Sphere_Sweep->Segment[i].Center_Coef[3],
  889.                               Sphere_Sweep->Segment[i].Center_Coef[1]);
  890.                 Coef[2] *= 2.0;
  891.                 VDot(temp, Sphere_Sweep->Segment[i].Center_Coef[2],
  892.                            Sphere_Sweep->Segment[i].Center_Coef[2]);
  893.                 Coef[2] += temp;
  894.                 Coef[2] -= 2.0 * Sphere_Sweep->Segment[i].Radius_Coef[3]
  895.                                * Sphere_Sweep->Segment[i].Radius_Coef[1];
  896.                 Coef[2] -= Sqr(Sphere_Sweep->Segment[i].Radius_Coef[2]);
  897.                 
  898.                 /* Coefficient for u^3 */
  899.                 VDot(Coef[3], Sphere_Sweep->Segment[i].Center_Coef[3], Vector);
  900.                 VDot(temp, Sphere_Sweep->Segment[i].Center_Coef[2],
  901.                            Sphere_Sweep->Segment[i].Center_Coef[1]);
  902.                 Coef[3] -= temp;
  903.                 Coef[3] += Sphere_Sweep->Segment[i].Radius_Coef[3]
  904.                          * Sphere_Sweep->Segment[i].Radius_Coef[0];
  905.                 Coef[3] += Sphere_Sweep->Segment[i].Radius_Coef[2]
  906.                          * Sphere_Sweep->Segment[i].Radius_Coef[1];
  907.                 Coef[3] *= -2.0;
  908.                 
  909.                 /* Coefficient for u^2 */
  910.                 VDot(Coef[4], Sphere_Sweep->Segment[i].Center_Coef[2], Vector);
  911.                 Coef[4] += Sphere_Sweep->Segment[i].Radius_Coef[2]
  912.                          * Sphere_Sweep->Segment[i].Radius_Coef[0];
  913.                 Coef[4] *= -2.0;
  914.                 VDot(temp, Sphere_Sweep->Segment[i].Center_Coef[1],
  915.                            Sphere_Sweep->Segment[i].Center_Coef[1]);
  916.                 Coef[4] += temp;
  917.                 Coef[4] -= Sqr(Sphere_Sweep->Segment[i].Radius_Coef[1]);
  918.                 
  919.                 /* Coefficient for u^1 */
  920.                 VDot(Coef[5], Sphere_Sweep->Segment[i].Center_Coef[1], Vector);
  921.                 Coef[5] += Sphere_Sweep->Segment[i].Radius_Coef[1]
  922.                          * Sphere_Sweep->Segment[i].Radius_Coef[0];
  923.                 Coef[5] *= -2.0;
  924.                 
  925.                 /* Coefficient for u^0 */
  926.                 VDot(Coef[6], Vector, Vector);
  927.                 Coef[6] -= Sqr(Sphere_Sweep->Segment[i].Radius_Coef[0]);
  928.                 
  929.                 /* Find roots */
  930.                 Num_Poly_Roots = Solve_Polynomial(6, Coef, Root, TRUE, 1e-10);
  931.                 
  932.                 /* Test for interval [0, 1] */
  933.                 for (j = 0; j < Num_Poly_Roots; j++)
  934.                 {
  935.                     if (Root[j] >= 0.0 && Root[j] <= 1.0)
  936.                     {
  937.                         /* At least one root inside interval, */
  938.                         /* so we are inside the sphere sweep  */
  939.                         inside = TRUE;
  940.                         break;
  941.                     }
  942.                 }
  943.             }
  944.             break;
  945.     }
  946.     
  947.     if (Test_Flag(Object, INVERTED_FLAG))
  948.     {
  949.         inside = !inside;
  950.     }
  951.     
  952.     return (inside);
  953. }
  954.  
  955.  
  956.  
  957. /*****************************************************************************
  958. *
  959. * FUNCTION
  960. *
  961. *   Sphere_Sweep_Normal
  962. *
  963. * INPUT
  964. *
  965. *   Object, Intersection
  966. *
  967. * OUTPUT
  968. *
  969. *   Normal
  970. *
  971. * RETURNS
  972. *
  973. *    -
  974. *
  975. * AUTHOR
  976. *
  977. *   Jochen Lippert
  978. *
  979. * DESCRIPTION
  980. *
  981. *   Calculate the surface normal of a sphere sweep.
  982. *
  983. * CHANGES
  984. *
  985. *   -
  986. *
  987. ******************************************************************************/
  988.  
  989. static void Sphere_Sweep_Normal(VECTOR    Result,OBJECT* Object,INTERSECTION    *Inter)
  990. {
  991.     Assign_Vector(Result, Inter->INormal);
  992. }
  993.  
  994.  
  995.  
  996. /*****************************************************************************
  997. *
  998. * FUNCTION
  999. *
  1000. *   Copy_Sphere_Sweep
  1001. *
  1002. * INPUT
  1003. *
  1004. *   Object
  1005. *
  1006. * OUTPUT
  1007. *
  1008. *   -
  1009. *
  1010. * RETURNS
  1011. *
  1012. *   Copy of sphere sweep
  1013. *
  1014. * AUTHOR
  1015. *
  1016. *   Jochen Lippert
  1017. *   
  1018. * DESCRIPTION
  1019. *
  1020. *   Copy a sphere sweep.
  1021. *
  1022. * CHANGES
  1023. *
  1024. *   -
  1025. *
  1026. ******************************************************************************/
  1027.  
  1028. void *Copy_Sphere_Sweep(OBJECT    *Object)
  1029. {
  1030.     SPHERE_SWEEP    *New;
  1031.     int                i;
  1032.     
  1033.     New = Create_Sphere_Sweep();
  1034.  
  1035.     New->UV_Trans = Copy_Transform( Object->UV_Trans );
  1036.  
  1037.     New->Interpolation = ((SPHERE_SWEEP *)Object)->Interpolation;
  1038.     
  1039.     New->Num_Modeling_Spheres = ((SPHERE_SWEEP *)Object)->Num_Modeling_Spheres;
  1040.     New->Modeling_Sphere = (SPHSWEEP_SPH*)POV_MALLOC(New->Num_Modeling_Spheres * sizeof(SPHSWEEP_SPH), "modeling sphere");
  1041.     for (i = 0; i < New->Num_Modeling_Spheres; i++)
  1042.     {
  1043.         New->Modeling_Sphere[i] = ((SPHERE_SWEEP *)Object)->Modeling_Sphere[i];
  1044.     }
  1045.     
  1046.     New->Depth_Tolerance = ((SPHERE_SWEEP *)Object)->Depth_Tolerance;
  1047.     
  1048.     Compute_Sphere_Sweep(New);
  1049.     
  1050.     New->Trans = Copy_Transform(((SPHERE_SWEEP *)Object)->Trans);
  1051.     
  1052.     return (New);
  1053. }
  1054.  
  1055.  
  1056.  
  1057. /*****************************************************************************
  1058. *
  1059. * FUNCTION
  1060. *
  1061. *   Translate_Sphere_Sweep
  1062. *
  1063. * INPUT
  1064. *
  1065. *   Object, Vector, Transformation
  1066. *
  1067. * OUTPUT
  1068. *
  1069. *   Object
  1070. *
  1071. * RETURNS
  1072. *
  1073. *   -
  1074. *
  1075. * AUTHOR
  1076. *
  1077. *   Jochen Lippert
  1078. *
  1079. * DESCRIPTION
  1080. *
  1081. *   Translate a sphere sweep.
  1082. *
  1083. * CHANGES
  1084. *
  1085. *   -
  1086. *
  1087. ******************************************************************************/
  1088.  
  1089. static void Translate_Sphere_Sweep(OBJECT *Object,VECTOR Vector,TRANSFORM *Trans)
  1090. {
  1091.     SPHERE_SWEEP *Sphere_Sweep = (SPHERE_SWEEP *)Object;
  1092.     int        i;
  1093.     
  1094.     if (Sphere_Sweep->Trans == NULL)
  1095.     {
  1096.         for (i = 0; i < Sphere_Sweep->Num_Modeling_Spheres; i++)
  1097.         {
  1098.             VAddEq(Sphere_Sweep->Modeling_Sphere[i].Center, Vector);
  1099.         }
  1100.         Compute_Sphere_Sweep(Sphere_Sweep);
  1101.         Compute_Sphere_Sweep_BBox(Sphere_Sweep);
  1102.     }
  1103.     else
  1104.     {
  1105.         Transform_Sphere_Sweep(Object, Trans);
  1106.     }
  1107. }
  1108.  
  1109.  
  1110.  
  1111. /*****************************************************************************
  1112. *
  1113. * FUNCTION
  1114. *
  1115. *   Rotate_Sphere_Sweep
  1116. *
  1117. * INPUT
  1118. *
  1119. *   Object, Vector, Transformation
  1120. *
  1121. * OUTPUT
  1122. *
  1123. *   Object
  1124. *
  1125. * RETURNS
  1126. *
  1127. *   -
  1128. *
  1129. * AUTHOR
  1130. *
  1131. *   Jochen Lippert
  1132. *
  1133. * DESCRIPTION
  1134. *
  1135. *   Rotate a sphere sweep.
  1136. *
  1137. * CHANGES
  1138. *
  1139. *   -
  1140. *
  1141. ******************************************************************************/
  1142.  
  1143. static void Rotate_Sphere_Sweep(OBJECT *Object,VECTOR Vector,TRANSFORM *Trans)
  1144. {
  1145.     SPHERE_SWEEP *Sphere_Sweep = (SPHERE_SWEEP *)Object;
  1146.     int        i;
  1147.     
  1148.     if (Sphere_Sweep->Trans == NULL)
  1149.     {
  1150.         for (i = 0; i < Sphere_Sweep->Num_Modeling_Spheres; i++)
  1151.         {
  1152.             MTransPoint(Sphere_Sweep->Modeling_Sphere[i].Center, Sphere_Sweep->Modeling_Sphere[i].Center, Trans);
  1153.         }
  1154.         Compute_Sphere_Sweep(Sphere_Sweep);
  1155.         Compute_Sphere_Sweep_BBox(Sphere_Sweep);
  1156.     }
  1157.     else
  1158.     {
  1159.         Transform_Sphere_Sweep(Object, Trans);
  1160.     }
  1161. }
  1162.  
  1163.  
  1164.  
  1165. /*****************************************************************************
  1166. *
  1167. * FUNCTION
  1168. *
  1169. *   Scale_Sphere_Sweep
  1170. *
  1171. * INPUT
  1172. *
  1173. *   Object, Vector, Transformation
  1174. *
  1175. * OUTPUT
  1176. *
  1177. *   Object
  1178. *
  1179. * RETURNS
  1180. *
  1181. *   -
  1182. *
  1183. * AUTHOR
  1184. *
  1185. *   Jochen Lippert
  1186. *
  1187. * DESCRIPTION
  1188. *
  1189. *   Scale a sphere sweep.
  1190. *
  1191. * CHANGES
  1192. *
  1193. *   -
  1194. *
  1195. ******************************************************************************/
  1196.  
  1197. static void Scale_Sphere_Sweep(OBJECT *Object,VECTOR Vector,TRANSFORM *Trans)
  1198. {
  1199.   SPHERE_SWEEP *Sphere_Sweep = (SPHERE_SWEEP *) Object;
  1200.   int i;
  1201.   
  1202.   if ((Vector[X] != Vector[Y]) || (Vector[X] != Vector[Z]))
  1203.   {
  1204.     if (Sphere_Sweep->Trans == NULL)
  1205.     {
  1206.       Sphere_Sweep->Trans = Create_Transform();
  1207.     }
  1208.   }
  1209.  
  1210.   if (Sphere_Sweep->Trans == NULL)
  1211.   {
  1212.       for (i = 0; i < Sphere_Sweep->Num_Modeling_Spheres; i++)
  1213.       {
  1214.         VScaleEq(Sphere_Sweep->Modeling_Sphere[i].Center, Vector[X]);
  1215.         
  1216.         Sphere_Sweep->Modeling_Sphere[i].Radius *= Vector[X];
  1217.       }
  1218.       Compute_Sphere_Sweep(Sphere_Sweep);
  1219.       Compute_Sphere_Sweep_BBox(Sphere_Sweep);
  1220.   }
  1221.   else
  1222.   {
  1223.     Transform_Sphere_Sweep(Object, Trans);
  1224.   }
  1225. }
  1226.  
  1227.  
  1228.  
  1229. /*****************************************************************************
  1230. *
  1231. * FUNCTION
  1232. *
  1233. *   Invert_Sphere_Sweep
  1234. *
  1235. * INPUT
  1236. *
  1237. *   Object
  1238. *
  1239. * OUTPUT
  1240. *
  1241. *   Object
  1242. *
  1243. * RETURNS
  1244. *
  1245. *   -
  1246. *
  1247. * AUTHOR
  1248. *
  1249. *   Jochen Lippert
  1250. *
  1251. * DESCRIPTION
  1252. *
  1253. *   Invert a sphere sweep.
  1254. *
  1255. * CHANGES
  1256. *
  1257. *   -
  1258. *
  1259. ******************************************************************************/
  1260.  
  1261. static void Invert_Sphere_Sweep(OBJECT *Object)
  1262. {
  1263.     Invert_Flag(Object, INVERTED_FLAG);
  1264. }
  1265.  
  1266.  
  1267.  
  1268. /*****************************************************************************
  1269. *
  1270. * FUNCTION
  1271. *
  1272. *   Create_Sphere_Sweep
  1273. *
  1274. * INPUT
  1275. *
  1276. *    -
  1277. *
  1278. * OUTPUT
  1279. *
  1280. *    -
  1281. *
  1282. * RETURNS
  1283. *
  1284. *    Sphere_Sweep
  1285. *
  1286. * AUTHOR
  1287. *
  1288. *   Jochen Lippert
  1289. *
  1290. * DESCRIPTION
  1291. *
  1292. *   Create a new, "empty" sphere sweep (no modeling spheres).
  1293. *
  1294. * CHANGES
  1295. *
  1296. *   -
  1297. *
  1298. ******************************************************************************/
  1299.  
  1300. SPHERE_SWEEP *Create_Sphere_Sweep(void)
  1301. {
  1302.     SPHERE_SWEEP    *New;
  1303.     
  1304.     New = (SPHERE_SWEEP *)POV_MALLOC(sizeof(SPHERE_SWEEP), "sphere sweep");
  1305.     
  1306.     INIT_OBJECT_FIELDS(New, SPHERE_SWEEP_OBJECT, &Sphere_Sweep_Methods)
  1307.     
  1308.     New->Interpolation = -1;
  1309.     
  1310.     New->Num_Modeling_Spheres = 0;
  1311.     New->Modeling_Sphere = NULL;
  1312.     
  1313.     New->Num_Spheres = 0;
  1314.     New->Sphere = NULL;
  1315.     
  1316.     New->Num_Segments = 0;
  1317.     New->Segment = NULL;
  1318.     
  1319.     New->Depth_Tolerance = DEPTH_TOLERANCE;
  1320.     
  1321.     New->Trans = NULL;
  1322.     
  1323.     return (New);
  1324. }
  1325.  
  1326.  
  1327.  
  1328. /*****************************************************************************
  1329. *
  1330. * FUNCTION
  1331. *
  1332. *   Transform_Sphere_Sweep
  1333. *
  1334. * INPUT
  1335. *
  1336. *   Object, Transformation
  1337. *
  1338. * OUTPUT
  1339. *
  1340. *   Object
  1341. *
  1342. * RETURNS
  1343. *
  1344. *   -
  1345. *
  1346. * AUTHOR
  1347. *
  1348. *   Jochen Lippert
  1349. *
  1350. * DESCRIPTION
  1351. *
  1352. *   Transform a sphere sweep.
  1353. *
  1354. * CHANGES
  1355. *
  1356. *   -
  1357. *
  1358. ******************************************************************************/
  1359.  
  1360. void Transform_Sphere_Sweep(OBJECT *Object,TRANSFORM *Trans)
  1361. {
  1362.   SPHERE_SWEEP *Sphere_Sweep = (SPHERE_SWEEP *)Object;
  1363.  
  1364.   if (Sphere_Sweep->Trans == NULL)
  1365.   {
  1366.     Sphere_Sweep->Trans = Create_Transform();
  1367.   }
  1368.  
  1369.   Compose_Transforms(Sphere_Sweep->Trans, Trans);
  1370.   
  1371.   Compute_Sphere_Sweep_BBox(Sphere_Sweep);
  1372. }
  1373.  
  1374.  
  1375.  
  1376. /*****************************************************************************
  1377. *
  1378. * FUNCTION
  1379. *
  1380. *   Destroy_Sphere_Sweep
  1381. *
  1382. * INPUT
  1383. *
  1384. *   Object
  1385. *
  1386. * OUTPUT
  1387. *
  1388. *   -
  1389. *
  1390. * RETURNS
  1391. *
  1392. *   -
  1393. *
  1394. * AUTHOR
  1395. *
  1396. *   Jochen Lippert
  1397. *
  1398. * DESCRIPTION
  1399. *
  1400. *   Free memory allocated for a sphere sweep.
  1401. *
  1402. * CHANGES
  1403. *
  1404. *   -
  1405. *
  1406. ******************************************************************************/
  1407.  
  1408. void Destroy_Sphere_Sweep(OBJECT *Object)
  1409. {
  1410.   Destroy_Transform(((SPHERE_SWEEP *)Object)->Trans);
  1411.   
  1412.   POV_FREE(((SPHERE_SWEEP *)Object)->Modeling_Sphere);
  1413.   POV_FREE(((SPHERE_SWEEP *)Object)->Sphere);
  1414.   POV_FREE(((SPHERE_SWEEP *)Object)->Segment);
  1415.  
  1416.   POV_FREE (Object);
  1417. }
  1418.  
  1419.  
  1420.  
  1421. /*****************************************************************************
  1422. *
  1423. * FUNCTION
  1424. *
  1425. *   Compute_Sphere_Sweep_BBox
  1426. *
  1427. * INPUT
  1428. *
  1429. *   Sphere Sweep
  1430. *
  1431. * OUTPUT
  1432. *
  1433. *   Sphere Sweep
  1434. *
  1435. * RETURNS
  1436. *
  1437. *   -
  1438. *
  1439. * AUTHOR
  1440. *
  1441. *   Jochen Lippert
  1442. *
  1443. * DESCRIPTION
  1444. *
  1445. *   Calculate the bounding box of a sphere sweep.
  1446. *
  1447. * CHANGES
  1448. *
  1449. *   -
  1450. *
  1451. ******************************************************************************/
  1452.  
  1453. void Compute_Sphere_Sweep_BBox(SPHERE_SWEEP *Sphere_Sweep)
  1454. {
  1455.     VECTOR    mins;
  1456.     VECTOR    maxs;
  1457.     int        i;
  1458.     
  1459.     mins[X] = mins[Y] = mins[Z] = maxs[X] = maxs[Y] = maxs[Z] = 0.0;
  1460.     
  1461.     for (i = 0; i < Sphere_Sweep->Num_Modeling_Spheres; i++)
  1462.     {
  1463.         if (Sphere_Sweep->Interpolation == CATMULL_ROM_SPLINE_SPHERE_SWEEP)
  1464.         {
  1465.         /* Make box a bit larger for Catmull-Rom-Spline sphere sweeps */
  1466.             mins[X] = min(mins[X], Sphere_Sweep->Modeling_Sphere[i].Center[X]
  1467.                              - 2 * Sphere_Sweep->Modeling_Sphere[i].Radius);
  1468.             mins[Y] = min(mins[Y], Sphere_Sweep->Modeling_Sphere[i].Center[Y]
  1469.                              - 2 * Sphere_Sweep->Modeling_Sphere[i].Radius);
  1470.             mins[Z] = min(mins[Z], Sphere_Sweep->Modeling_Sphere[i].Center[Z]
  1471.                              - 2 * Sphere_Sweep->Modeling_Sphere[i].Radius);
  1472.             maxs[X] = max(maxs[X], Sphere_Sweep->Modeling_Sphere[i].Center[X]
  1473.                              + 2 * Sphere_Sweep->Modeling_Sphere[i].Radius);
  1474.             maxs[Y] = max(maxs[Y], Sphere_Sweep->Modeling_Sphere[i].Center[Y]
  1475.                              + 2 * Sphere_Sweep->Modeling_Sphere[i].Radius);
  1476.             maxs[Z] = max(maxs[Z], Sphere_Sweep->Modeling_Sphere[i].Center[Z]
  1477.                              + 2 * Sphere_Sweep->Modeling_Sphere[i].Radius);
  1478.         }
  1479.         else
  1480.         {
  1481.             mins[X] = min(mins[X], Sphere_Sweep->Modeling_Sphere[i].Center[X]
  1482.                                  - Sphere_Sweep->Modeling_Sphere[i].Radius);
  1483.             mins[Y] = min(mins[Y], Sphere_Sweep->Modeling_Sphere[i].Center[Y]
  1484.                                  - Sphere_Sweep->Modeling_Sphere[i].Radius);
  1485.             mins[Z] = min(mins[Z], Sphere_Sweep->Modeling_Sphere[i].Center[Z]
  1486.                                  - Sphere_Sweep->Modeling_Sphere[i].Radius);
  1487.             maxs[X] = max(maxs[X], Sphere_Sweep->Modeling_Sphere[i].Center[X]
  1488.                                  + Sphere_Sweep->Modeling_Sphere[i].Radius);
  1489.             maxs[Y] = max(maxs[Y], Sphere_Sweep->Modeling_Sphere[i].Center[Y]
  1490.                                  + Sphere_Sweep->Modeling_Sphere[i].Radius);
  1491.             maxs[Z] = max(maxs[Z], Sphere_Sweep->Modeling_Sphere[i].Center[Z]
  1492.                                  + Sphere_Sweep->Modeling_Sphere[i].Radius);
  1493.         }
  1494.     }
  1495.     
  1496.     Make_BBox_from_min_max(Sphere_Sweep->BBox, mins, maxs);
  1497.     
  1498.     if (Sphere_Sweep->Trans != NULL)
  1499.     {
  1500.         Recompute_BBox(&Sphere_Sweep->BBox, Sphere_Sweep->Trans);
  1501.     }
  1502. }
  1503.  
  1504.  
  1505.  
  1506. /*****************************************************************************
  1507. *
  1508. * FUNCTION
  1509. *
  1510. *   Compute_Sphere_Sweep
  1511. *
  1512. * INPUT
  1513. *
  1514. *   Sphere sweep
  1515. *
  1516. * OUTPUT
  1517. *
  1518. *   Sphere sweep
  1519. *
  1520. * RETURNS
  1521. *
  1522. *   -
  1523. *
  1524. * AUTHOR
  1525. *
  1526. *   Jochen Lippert
  1527. *
  1528. * DESCRIPTION
  1529. *
  1530. *   Calculate the internal representation of a sphere sweep.
  1531. *
  1532. * CHANGES
  1533. *
  1534. *   -
  1535. *
  1536. ******************************************************************************/
  1537.  
  1538. void Compute_Sphere_Sweep(SPHERE_SWEEP *Sphere_Sweep)
  1539. {
  1540.     size_t    size;
  1541.     int        i;
  1542.     int        coef;
  1543.     int        msph;
  1544.     int        last_sph;
  1545.     int        last_seg;
  1546.     
  1547.     switch (Sphere_Sweep->Interpolation)
  1548.     {
  1549.         case LINEAR_SPHERE_SWEEP:
  1550.             
  1551.             /* Allocate memory if neccessary */
  1552.             if (Sphere_Sweep->Segment == NULL)
  1553.             {
  1554.                 Sphere_Sweep->Num_Segments = Sphere_Sweep->Num_Modeling_Spheres - 1;
  1555.                 size = Sphere_Sweep->Num_Segments * sizeof(SPHSWEEP_SEG);
  1556.                 Sphere_Sweep->Segment = (SPHSWEEP_SEG *)POV_MALLOC(size, "sphere sweep segments");
  1557.             }
  1558.             
  1559.             /* Calculate polynomials for each segment */
  1560.             for (i = 0; i < Sphere_Sweep->Num_Segments; i++)
  1561.             {
  1562.                 /* Polynomial has two coefficients */
  1563.                 Sphere_Sweep->Segment[i].Num_Coefs = 2;
  1564.                 
  1565.                 /* Coefficients for u^1 */
  1566.                 VSub(Sphere_Sweep->Segment[i].Center_Coef[1],
  1567.                      Sphere_Sweep->Modeling_Sphere[i + 1].Center,
  1568.                      Sphere_Sweep->Modeling_Sphere[i].Center);
  1569.                 
  1570.                 Sphere_Sweep->Segment[i].Radius_Coef[1] =
  1571.                      Sphere_Sweep->Modeling_Sphere[i + 1].Radius -
  1572.                      Sphere_Sweep->Modeling_Sphere[i].Radius;
  1573.                 
  1574.                 /* Coefficients for u^0 */
  1575.                 Assign_Vector(Sphere_Sweep->Segment[i].Center_Coef[0],
  1576.                      Sphere_Sweep->Modeling_Sphere[i].Center);
  1577.                 
  1578.                 Sphere_Sweep->Segment[i].Radius_Coef[0] =
  1579.                      Sphere_Sweep->Modeling_Sphere[i].Radius;
  1580.             }    
  1581.             
  1582.             break;    /* case LINEAR_SPHERE_SWEEP */
  1583.         
  1584.         case CATMULL_ROM_SPLINE_SPHERE_SWEEP:
  1585.             
  1586.             /* Allocate memory if neccessary */
  1587.             if (Sphere_Sweep->Segment == NULL)
  1588.             {
  1589.                 Sphere_Sweep->Num_Segments = Sphere_Sweep->Num_Modeling_Spheres - 3;
  1590.                 size = Sphere_Sweep->Num_Segments * sizeof(SPHSWEEP_SEG);
  1591.                 Sphere_Sweep->Segment = (SPHSWEEP_SEG *)POV_MALLOC(size, "sphere sweep segments");
  1592.             }
  1593.             
  1594.             /* Calculate polynomials for each segment */
  1595.             for (i = 0; i < Sphere_Sweep->Num_Segments; i++)
  1596.             {
  1597.                 /* Polynomial has four coefficients */
  1598.                 Sphere_Sweep->Segment[i].Num_Coefs = 4;
  1599.                 
  1600.                 /* Calculate coefficients */
  1601.                 for (coef = 0; coef < 4; coef++)
  1602.                 {
  1603.                     /* Center */
  1604.                     VScale(Sphere_Sweep->Segment[i].Center_Coef[coef],
  1605.                            Sphere_Sweep->Modeling_Sphere[i].Center,
  1606.                            Catmull_Rom_Matrix[coef][0]);
  1607.                     
  1608.                     /* Radius */
  1609.                     Sphere_Sweep->Segment[i].Radius_Coef[coef] =
  1610.                            Sphere_Sweep->Modeling_Sphere[i].Radius *
  1611.                            Catmull_Rom_Matrix[coef][0];
  1612.                     
  1613.                     for (msph = 1; msph < 4; msph++)
  1614.                     {
  1615.                         /* Center */
  1616.                         VAddScaledEq(Sphere_Sweep->Segment[i].Center_Coef[coef],
  1617.                                      Catmull_Rom_Matrix[coef][msph],
  1618.                                      Sphere_Sweep->Modeling_Sphere[i + msph].Center);
  1619.                         
  1620.                         /* Radius */
  1621.                         Sphere_Sweep->Segment[i].Radius_Coef[coef] +=
  1622.                                      Catmull_Rom_Matrix[coef][msph] *
  1623.                                      Sphere_Sweep->Modeling_Sphere[i + msph].Radius;
  1624.                     }
  1625.                 }
  1626.             }
  1627.             
  1628.             break;    /* case CATMULL_ROM_SPLINE_SPHERE_SWEEP */
  1629.         
  1630.         case B_SPLINE_SPHERE_SWEEP:
  1631.             
  1632.             /* Allocate memory if neccessary */
  1633.             if (Sphere_Sweep->Segment == NULL)
  1634.             {
  1635.                 Sphere_Sweep->Num_Segments = Sphere_Sweep->Num_Modeling_Spheres - 3;
  1636.                 size = Sphere_Sweep->Num_Segments * sizeof(SPHSWEEP_SEG);
  1637.                 Sphere_Sweep->Segment = (SPHSWEEP_SEG *)POV_MALLOC(size, "sphere sweep segments");
  1638.             }
  1639.             
  1640.             /* Calculate polynomials for each segment */
  1641.             for (i = 0; i < Sphere_Sweep->Num_Segments; i++)
  1642.             {
  1643.                 /* Polynomial has four coefficients */
  1644.                 Sphere_Sweep->Segment[i].Num_Coefs = 4;
  1645.                 
  1646.                 /* Calculate coefficients */
  1647.                 for (coef = 0; coef < 4; coef++)
  1648.                 {
  1649.                     /* Center */
  1650.                     VScale(Sphere_Sweep->Segment[i].Center_Coef[coef],
  1651.                            Sphere_Sweep->Modeling_Sphere[i].Center,
  1652.                            B_Matrix[coef][0]);
  1653.                     
  1654.                     /* Radius */
  1655.                     Sphere_Sweep->Segment[i].Radius_Coef[coef] =
  1656.                            Sphere_Sweep->Modeling_Sphere[i].Radius *
  1657.                            B_Matrix[coef][0];
  1658.                     
  1659.                     for (msph = 1; msph < 4; msph++)
  1660.                     {
  1661.                         /* Center */
  1662.                         VAddScaledEq(Sphere_Sweep->Segment[i].Center_Coef[coef],
  1663.                                      B_Matrix[coef][msph],
  1664.                                      Sphere_Sweep->Modeling_Sphere[i + msph].Center);
  1665.                         
  1666.                         /* Radius */
  1667.                         Sphere_Sweep->Segment[i].Radius_Coef[coef] +=
  1668.                                      B_Matrix[coef][msph] *
  1669.                                      Sphere_Sweep->Modeling_Sphere[i + msph].Radius;
  1670.                     }
  1671.                 }
  1672.             }
  1673.             
  1674.             break;    /* case B_SPLINE_SPHERE_SWEEP */
  1675.     }
  1676.     
  1677.     /*
  1678.      * Pre-calculate several constants
  1679.      */
  1680.     
  1681.     for (i = 0; i < Sphere_Sweep->Num_Segments; i++)
  1682.     {
  1683.         /*
  1684.          * Calculate closing sphere for u = 0
  1685.          */
  1686.         
  1687.         /* Center */
  1688.         Assign_Vector(Sphere_Sweep->Segment[i].Closing_Sphere[0].Center,
  1689.                       Sphere_Sweep->Segment[i].Center_Coef[0]);
  1690.         
  1691.         /* Radius */
  1692.         Sphere_Sweep->Segment[i].Closing_Sphere[0].Radius =
  1693.                       Sphere_Sweep->Segment[i].Radius_Coef[0];
  1694.         
  1695.         /*
  1696.          * Calculate derivatives for u = 0
  1697.          */
  1698.         
  1699.         /* Center */
  1700.         Assign_Vector(Sphere_Sweep->Segment[i].Center_Deriv[0],
  1701.                       Sphere_Sweep->Segment[i].Center_Coef[1]);
  1702.         
  1703.         /* Radius */
  1704.         Sphere_Sweep->Segment[i].Radius_Deriv[0] =
  1705.                       Sphere_Sweep->Segment[i].Radius_Coef[1];
  1706.         
  1707.         /*
  1708.          * Calculate closing sphere for u = 1
  1709.          */
  1710.         
  1711.         /* Center */
  1712.         Assign_Vector(Sphere_Sweep->Segment[i].Closing_Sphere[1].Center,
  1713.                       Sphere_Sweep->Segment[i].Center_Coef[0]);
  1714.         
  1715.         /* Radius */
  1716.         Sphere_Sweep->Segment[i].Closing_Sphere[1].Radius =
  1717.                       Sphere_Sweep->Segment[i].Radius_Coef[0];
  1718.         
  1719.         for (coef = 1; coef < Sphere_Sweep->Segment[i].Num_Coefs; coef++)
  1720.         {
  1721.             /* Center */
  1722.             VAddEq(Sphere_Sweep->Segment[i].Closing_Sphere[1].Center,
  1723.                    Sphere_Sweep->Segment[i].Center_Coef[coef]);
  1724.             
  1725.             /* Radius */
  1726.             Sphere_Sweep->Segment[i].Closing_Sphere[1].Radius +=
  1727.                    Sphere_Sweep->Segment[i].Radius_Coef[coef];
  1728.         }
  1729.         
  1730.         /*
  1731.          * Calculate derivatives for u = 1
  1732.          */
  1733.         
  1734.         /* Center */
  1735.         Assign_Vector(Sphere_Sweep->Segment[i].Center_Deriv[1],
  1736.                       Sphere_Sweep->Segment[i].Center_Coef[1]);
  1737.         
  1738.         /* Radius */
  1739.         Sphere_Sweep->Segment[i].Radius_Deriv[1] =
  1740.                       Sphere_Sweep->Segment[i].Radius_Coef[1];
  1741.         
  1742.         for (coef = 2; coef < Sphere_Sweep->Segment[i].Num_Coefs; coef++)
  1743.         {
  1744.             /* Center */
  1745.             VAddScaledEq(Sphere_Sweep->Segment[i].Center_Deriv[1], coef,
  1746.                          Sphere_Sweep->Segment[i].Center_Coef[coef]);
  1747.             
  1748.             /* Radius */
  1749.             Sphere_Sweep->Segment[i].Radius_Deriv[1] += coef *
  1750.                          Sphere_Sweep->Segment[i].Radius_Coef[coef];
  1751.         }
  1752.     }
  1753.     
  1754.     /*
  1755.      * Calculate single spheres
  1756.      */
  1757.     
  1758.     /* Allocate memory if neccessary */
  1759.     if (Sphere_Sweep->Sphere == NULL)
  1760.     {
  1761.         Sphere_Sweep->Num_Spheres = Sphere_Sweep->Num_Segments + 1;
  1762.         size = Sphere_Sweep->Num_Spheres * sizeof(SPHSWEEP_SPH);
  1763.         Sphere_Sweep->Sphere = (SPHSWEEP_SPH *)POV_MALLOC(size, "sphere sweep spheres");
  1764.     }
  1765.     
  1766.     /*
  1767.      * Calculate first sphere of every segment
  1768.      */
  1769.     
  1770.     for (i = 0; i < Sphere_Sweep->Num_Segments; i++)
  1771.     {
  1772.         /* Center */
  1773.         Assign_Vector(Sphere_Sweep->Sphere[i].Center,
  1774.                       Sphere_Sweep->Segment[i].Center_Coef[0]);
  1775.         
  1776.         /* Radius */
  1777.         Sphere_Sweep->Sphere[i].Radius =
  1778.                       Sphere_Sweep->Segment[i].Radius_Coef[0];
  1779.     }
  1780.     
  1781.     /*
  1782.      * Calculate last sphere of last segment
  1783.      */
  1784.     
  1785.     last_sph = Sphere_Sweep->Num_Spheres - 1;
  1786.     last_seg = Sphere_Sweep->Num_Segments - 1;
  1787.     
  1788.     /* Center */
  1789.     Assign_Vector(Sphere_Sweep->Sphere[last_sph].Center,
  1790.                   Sphere_Sweep->Segment[last_seg].Center_Coef[0]);
  1791.     /* Radius */
  1792.     Sphere_Sweep->Sphere[last_sph].Radius =
  1793.                   Sphere_Sweep->Segment[last_seg].Radius_Coef[0];
  1794.     
  1795.     for (coef = 1; coef < Sphere_Sweep->Segment[last_seg].Num_Coefs; coef++)
  1796.     {
  1797.         /* Center */
  1798.         VAddEq(Sphere_Sweep->Sphere[last_sph].Center,
  1799.                Sphere_Sweep->Segment[last_seg].Center_Coef[coef]);
  1800.         
  1801.         /* Radius */
  1802.         Sphere_Sweep->Sphere[last_sph].Radius +=
  1803.                Sphere_Sweep->Segment[last_seg].Radius_Coef[coef];
  1804.     }
  1805. }
  1806.  
  1807.  
  1808.  
  1809. /*****************************************************************************
  1810. *
  1811. * FUNCTION
  1812. *
  1813. *   Find_Valid_Points
  1814. *
  1815. * INPUT
  1816. *
  1817. *   Intersection list, number of intersections, ray
  1818. *
  1819. * OUTPUT
  1820. *
  1821. *   Intersection list
  1822. *
  1823. * RETURNS
  1824. *
  1825. *   Number of valid intersections
  1826. *
  1827. * AUTHOR
  1828. *
  1829. *   Jochen Lippert
  1830. *
  1831. * DESCRIPTION
  1832. *
  1833. *   Delete invalid intersections.
  1834. *
  1835. * CHANGES
  1836. *
  1837. *   -
  1838. *
  1839. ******************************************************************************/
  1840.  
  1841. int Find_Valid_Points(SPHSWEEP_INT    *Inter, int    Num_Inter, RAY*Ray)
  1842. {
  1843.     int        i;
  1844.     int        j;
  1845.     int        Inside;
  1846.     int        Keep;
  1847.     DBL        NormalDotDirection;
  1848.     
  1849.     Inside = 1;
  1850.     i = 1;
  1851.     while (i < Num_Inter - 1)
  1852.     {
  1853.         /* Angle between normal and ray */
  1854.         VDot(NormalDotDirection, Inter[i].Normal, Ray->Direction);
  1855.         
  1856.         /* Does the ray enter the part? */
  1857.         if (NormalDotDirection < 0.0)
  1858.         {
  1859.             /* Ray enters part, keep intersection if ray was outside any part */
  1860.             Keep = (Inside == 0);
  1861.             
  1862.             /* increase inside counter */
  1863.             Inside++;
  1864.         }
  1865.         else
  1866.         {
  1867.             /* Ray exits part, keep intersection if ray was inside one part */
  1868.             Keep = (Inside == 1);
  1869.             
  1870.             /* decrease inside counter */
  1871.             Inside--;
  1872.         }
  1873.         
  1874.         /* Keep intersection? */
  1875.         if (Keep)
  1876.         {
  1877.             /* Yes, advance to next one */
  1878.             i++;
  1879.         }
  1880.         else
  1881.         {
  1882.             /* No, delete it */
  1883.             for (j = i + 1; j < Num_Inter; j++)
  1884.             {
  1885.                 Inter[j - 1] = Inter[j];
  1886.             }
  1887.             Num_Inter--;
  1888.         }
  1889.     }
  1890.     return (Num_Inter);
  1891. }
  1892.  
  1893.  
  1894.  
  1895. /*****************************************************************************
  1896. *
  1897. * FUNCTION
  1898. *
  1899. *   Comp_Isects
  1900. *
  1901. * INPUT
  1902. *   
  1903. * OUTPUT
  1904. *   
  1905. * RETURNS
  1906. *   
  1907. * AUTHOR
  1908. *
  1909. *   Jochen Lippert
  1910. *   
  1911. * DESCRIPTION
  1912. *
  1913. *   Compare the depths of two sphere sweep intersections.
  1914. *
  1915. * CHANGES
  1916. *
  1917. *   -
  1918. *
  1919. ******************************************************************************/
  1920.  
  1921. static int Comp_Isects(CONST void    *Intersection_1,CONST void    *Intersection_2)
  1922. {
  1923.     SPHSWEEP_INT    *Int_1;
  1924.     SPHSWEEP_INT    *Int_2;
  1925.     
  1926.     Int_1 = (SPHSWEEP_INT *)Intersection_1;
  1927.     Int_2 = (SPHSWEEP_INT *)Intersection_2;
  1928.     
  1929.     if (Int_1->t < Int_2->t)
  1930.     {
  1931.         return (-1);
  1932.     }
  1933.     
  1934.     if (Int_1->t == Int_2->t)
  1935.     {
  1936.         return (0);
  1937.     }
  1938.     else
  1939.     {
  1940.         return (1);
  1941.     }
  1942. }
  1943.