home *** CD-ROM | disk | FTP | other *** search
/ Chestnut's Multimedia Mania / MM_MANIA.ISO / graphics / povsrc20 / bezier.c < prev    next >
C/C++ Source or Header  |  1993-08-04  |  34KB  |  1,149 lines
  1. /****************************************************************************
  2. *                   bezier.c
  3. *
  4. *  This module implements the code for Bezier bicubic patch shapes
  5. *
  6. *  This file was written by Alexander Enzmann.    He wrote the code for
  7. *  bezier bicubic patches and generously provided us these enhancements.
  8. *
  9. *  from Persistence of Vision Raytracer
  10. *  Copyright 1993 Persistence of Vision Team
  11. *---------------------------------------------------------------------------
  12. *  NOTICE: This source code file is provided so that users may experiment
  13. *  with enhancements to POV-Ray and to port the software to platforms other 
  14. *  than those supported by the POV-Ray Team.  There are strict rules under
  15. *  which you are permitted to use this file.  The rules are in the file
  16. *  named POVLEGAL.DOC which should be distributed with this file. If 
  17. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  18. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  19. *  Forum.  The latest version of POV-Ray may be found there as well.
  20. *
  21. * This program is based on the popular DKB raytracer version 2.12.
  22. * DKBTrace was originally written by David K. Buck.
  23. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  24. *
  25. *****************************************************************************/
  26.  
  27. #include "frame.h"
  28. #include "vector.h"
  29. #include "povproto.h"
  30.  
  31. METHODS Bicubic_Patch_Methods =
  32.   { 
  33.   All_Bicubic_Patch_Intersections,
  34.   Inside_Bicubic_Patch, Bicubic_Patch_Normal, Copy_Bicubic_Patch,
  35.   Translate_Bicubic_Patch, Rotate_Bicubic_Patch,
  36.   Scale_Bicubic_Patch, Transform_Bicubic_Patch, Invert_Bicubic_Patch,
  37.   Destroy_Bicubic_Patch
  38.   };
  39.  
  40. extern long Ray_Bicubic_Tests, Ray_Bicubic_Tests_Succeeded;
  41. extern unsigned int Options;
  42. extern RAY *CM_Ray;
  43. extern int Shadow_Test_Flag;
  44.  
  45. int max_depth_reached;
  46.  
  47. static int InvertMatrix PARAMS((VECTOR in[3], VECTOR out[3]));
  48. static void bezier_value PARAMS((VECTOR (*Control_Points)[4][4], DBL u0, DBL v0,
  49. VECTOR *P, VECTOR *N));
  50. static int intersect_subpatch PARAMS((BICUBIC_PATCH *, RAY *, VECTOR [3],
  51. DBL [3], DBL [3], DBL *, VECTOR *, VECTOR *, DBL *, DBL *));
  52. static void find_average PARAMS((int, VECTOR *, VECTOR *, DBL *));
  53. static int spherical_bounds_check PARAMS((RAY *, VECTOR *, DBL));
  54. static int intersect_bicubic_patch0 PARAMS((RAY *, BICUBIC_PATCH *, DBL *));
  55. static int intersect_bicubic_patch1 PARAMS((RAY *, BICUBIC_PATCH *, DBL *));
  56. static DBL point_plane_distance PARAMS((VECTOR *, VECTOR *, DBL *));
  57. static DBL determine_subpatch_flatness PARAMS((VECTOR (*)[4][4]));
  58. static int flat_enough PARAMS((BICUBIC_PATCH *, VECTOR (*)[4][4]));
  59. static void bezier_bounding_sphere PARAMS((VECTOR (*)[4][4], VECTOR *,DBL *));
  60. static void bezier_subpatch_intersect PARAMS((RAY *, BICUBIC_PATCH *,
  61. VECTOR (*)[4][4], DBL, DBL, DBL, DBL,
  62. int *, DBL *));
  63. static void bezier_split_left_right PARAMS((VECTOR (*)[4][4],VECTOR (*)[4][4],
  64. VECTOR (*)[4][4]));
  65. static void bezier_split_up_down PARAMS((VECTOR (*)[4][4], VECTOR (*)[4][4],
  66. VECTOR (*)[4][4]));
  67. static void bezier_subdivider PARAMS((RAY *, BICUBIC_PATCH *,VECTOR (*)[4][4],
  68. DBL, DBL, DBL, DBL, int, int *, DBL *));
  69. static void bezier_tree_deleter PARAMS((BEZIER_NODE *Node));
  70. static BEZIER_NODE *bezier_tree_builder PARAMS((BICUBIC_PATCH *Object,
  71. VECTOR(*Patch)[4][4], DBL u0, DBL u1,
  72. DBL v0, DBL v1, int depth));
  73. static void bezier_tree_walker PARAMS((RAY *, BICUBIC_PATCH *, BEZIER_NODE *,
  74. int, int *, DBL *));
  75. static BEZIER_NODE *create_new_bezier_node PARAMS((void));
  76. static BEZIER_VERTICES *create_bezier_vertex_block PARAMS((void));
  77. static BEZIER_CHILDREN *create_bezier_child_block PARAMS((void));
  78. static int subpatch_normal PARAMS((VECTOR *v1, VECTOR *v2, VECTOR *v3,
  79. VECTOR *Result, DBL *d));
  80.  
  81.   static DBL C[4] = {
  82.   1.0, 3.0, 3.0, 1.0
  83.   };
  84.  
  85. #define BEZIER_EPSILON 1.0e-10
  86. #define BEZIER_TOLERANCE 1.0e-5
  87.  
  88. static BEZIER_NODE *create_new_bezier_node()
  89.   {
  90.   BEZIER_NODE *Node = (BEZIER_NODE *)malloc(sizeof(BEZIER_NODE));
  91.   if (Node == NULL) 
  92.     MAError("bezier node");
  93.   Node->Data_Ptr = NULL;
  94.   return Node;
  95.   }
  96.  
  97. static BEZIER_VERTICES *create_bezier_vertex_block()
  98.   {
  99.   BEZIER_VERTICES *Vertices = (BEZIER_VERTICES *)malloc(sizeof(BEZIER_VERTICES));
  100.   if (Vertices == NULL) 
  101.     MAError("bezier vertices");
  102.   return Vertices;
  103.   }
  104.  
  105. static BEZIER_CHILDREN *create_bezier_child_block()
  106.   {
  107.   BEZIER_CHILDREN *Children = (BEZIER_CHILDREN *)malloc(sizeof(BEZIER_CHILDREN));
  108.   if (Children == NULL) 
  109.     MAError("bezier children");
  110.   return Children;
  111.   }
  112.  
  113. static BEZIER_NODE *bezier_tree_builder(Object, Patch, u0, u1, v0, v1, depth)
  114. BICUBIC_PATCH *Object;
  115. VECTOR (*Patch)[4][4];
  116. DBL u0, u1, v0, v1;
  117. int depth;
  118.   {
  119.   VECTOR Lower_Left[4][4], Lower_Right[4][4];
  120.   VECTOR Upper_Left[4][4], Upper_Right[4][4];
  121.   BEZIER_CHILDREN *Children;
  122.   BEZIER_VERTICES *Vertices;
  123.   BEZIER_NODE *Node = create_new_bezier_node();
  124.  
  125.   if (depth > max_depth_reached) max_depth_reached = depth;
  126.  
  127.   /* Build the bounding sphere for this subpatch */
  128.   bezier_bounding_sphere(Patch, &(Node->Center), &(Node->Radius_Squared));
  129.  
  130.   /* If the patch is close to being flat, then just perform a ray-plane
  131.       intersection test. */
  132.   if (flat_enough(Object, Patch)) 
  133.     {
  134.     /* The patch is now flat enough to simply store the corners */
  135.     Node->Node_Type = BEZIER_LEAF_NODE;
  136.     Vertices = create_bezier_vertex_block();
  137.     Vertices->Vertices[0] = (*Patch)[0][0];
  138.     Vertices->Vertices[1] = (*Patch)[0][3];
  139.     Vertices->Vertices[2] = (*Patch)[3][3];
  140.     Vertices->Vertices[3] = (*Patch)[3][0];
  141.     Vertices->uvbnds[0] = u0;
  142.     Vertices->uvbnds[1] = u1;
  143.     Vertices->uvbnds[2] = v0;
  144.     Vertices->uvbnds[3] = v1;
  145.     Node->Data_Ptr = (void *)Vertices;
  146.     }
  147.   else if (depth >= Object->U_Steps)
  148.     if (depth >= Object->V_Steps) 
  149.       {
  150.       /* We are at the max recursion depth. Just store corners. */
  151.       Node->Node_Type = BEZIER_LEAF_NODE;
  152.       Vertices = create_bezier_vertex_block();
  153.       Vertices->Vertices[0] = (*Patch)[0][0];
  154.       Vertices->Vertices[1] = (*Patch)[0][3];
  155.       Vertices->Vertices[2] = (*Patch)[3][3];
  156.       Vertices->Vertices[3] = (*Patch)[3][0];
  157.       Vertices->uvbnds[0] = u0;
  158.       Vertices->uvbnds[1] = u1;
  159.       Vertices->uvbnds[2] = v0;
  160.       Vertices->uvbnds[3] = v1;
  161.       Node->Data_Ptr = (void *)Vertices;
  162.       }
  163.     else 
  164.       {
  165.       bezier_split_up_down(Patch, (VECTOR (*)[4][4])Lower_Left,
  166.         (VECTOR (*)[4][4])Upper_Left);
  167.       Node->Node_Type = BEZIER_INTERIOR_NODE;
  168.       Children = create_bezier_child_block();
  169.       Children->Children[0] =
  170.       bezier_tree_builder(Object, (VECTOR (*)[4][4])Lower_Left,
  171.         u0, u1, v0, (v0 + v1) / 2.0, depth+1);
  172.       Children->Children[1] =
  173.       bezier_tree_builder(Object, (VECTOR (*)[4][4])Upper_Left,
  174.         u0, u1, (v0 + v1) / 2.0, v1, depth+1);
  175.       Node->Count = 2;
  176.       Node->Data_Ptr = (void *)Children;
  177.       }
  178.   else if (depth >= Object->V_Steps) 
  179.     {
  180.     bezier_split_left_right(Patch, (VECTOR (*)[4][4])Lower_Left,
  181.       (VECTOR (*)[4][4])Lower_Right);
  182.     Node->Node_Type = BEZIER_INTERIOR_NODE;
  183.     Children = create_bezier_child_block();
  184.     Children->Children[0] =
  185.     bezier_tree_builder(Object, (VECTOR (*)[4][4])Lower_Left,
  186.       u0, (u0 + u1) / 2.0, v0, v1, depth+1);
  187.     Children->Children[1] =
  188.     bezier_tree_builder(Object, (VECTOR (*)[4][4])Lower_Right,
  189.       (u0 + u1) / 2.0, u1, v0, v1, depth+1);
  190.     Node->Count = 2;
  191.     Node->Data_Ptr = (void *)Children;
  192.     }
  193.   else 
  194.     {
  195.     bezier_split_left_right(Patch, (VECTOR (*)[4][4])Lower_Left,
  196.       (VECTOR (*)[4][4])Lower_Right);
  197.     bezier_split_up_down((VECTOR (*)[4][4])Lower_Left,
  198.       (VECTOR (*)[4][4])Lower_Left,
  199.       (VECTOR (*)[4][4])Upper_Left);
  200.     bezier_split_up_down((VECTOR (*)[4][4])Lower_Right,
  201.       (VECTOR (*)[4][4])Lower_Right,
  202.       (VECTOR (*)[4][4])Upper_Right);
  203.     Node->Node_Type = BEZIER_INTERIOR_NODE;
  204.     Children = create_bezier_child_block();
  205.     Children->Children[0] =
  206.     bezier_tree_builder(Object, (VECTOR (*)[4][4])Lower_Left,
  207.       u0, (u0 + u1) / 2.0, v0, (v0 + v1) / 2.0, depth+1);
  208.     Children->Children[1] =
  209.     bezier_tree_builder(Object, (VECTOR (*)[4][4])Upper_Left,
  210.       u0, (u0 + u1) / 2.0, (v0 + v1) / 2.0, v1, depth+1);
  211.     Children->Children[2] =
  212.     bezier_tree_builder(Object, (VECTOR (*)[4][4])Lower_Right,
  213.       (u0 + u1) / 2.0, u1, v0, (v0 + v1) / 2.0, depth+1);
  214.     Children->Children[3] =
  215.     bezier_tree_builder(Object, (VECTOR (*)[4][4])Upper_Right,
  216.       (u0 + u1) / 2.0, u1, (v0 + v1) / 2.0, v1, depth+1);
  217.     Node->Count = 4;
  218.     Node->Data_Ptr = (void *)Children;
  219.     }
  220.   return Node;
  221.   }
  222.  
  223. /* Determine the position and normal at a single coordinate point (u, v)
  224.    on a Bezier patch */
  225. static void bezier_value(Control_Points, u0, v0, P, N)
  226. VECTOR (*Control_Points)[4][4];
  227. DBL u0, v0;
  228. VECTOR *P, *N;
  229.   {
  230.   DBL c, t, ut, vt;
  231.   DBL u[4], uu[4], v[4], vv[4];
  232.   DBL du[4], duu[4], dv[4], dvv[4];
  233.   int i, j;
  234.   VECTOR U, V, T;
  235.  
  236.   /* Calculate binomial coefficients times coordinate positions */
  237.   u[0]  = 1.0; uu[0] = 1.0; du[0] = 0.0; duu[0] = 0.0;
  238.   v[0]  = 1.0; vv[0] = 1.0; dv[0] = 0.0; dvv[0] = 0.0;
  239.   for (i=1;i<4;i++) 
  240.     {
  241.     u[i]   = u[i-1] * u0;
  242.     uu[i]  = uu[i-1] * (1.0 - u0);
  243.     v[i]   = v[i-1] * v0;
  244.     vv[i]  = vv[i-1] * (1.0 - v0);
  245.     du[i]  =  i * u[i-1];
  246.     duu[i] = -i * uu[i-1];
  247.     dv[i]  =  i * v[i-1];
  248.     dvv[i] = -i * vv[i-1];
  249.     }
  250.  
  251.   /* Now evaluate position and tangents based on control points */
  252.   Make_Vector(P, 0, 0, 0);
  253.   Make_Vector(&U, 0, 0, 0);
  254.   Make_Vector(&V, 0, 0, 0);
  255.   for (i=0;i<4;i++)
  256.     for (j=0;j<4;j++) 
  257.     {
  258.     c = C[i] * C[j];
  259.     ut = u[i] * uu[3 - i];
  260.     vt = v[j] * vv[3 - j];
  261.     t = c * ut * vt;
  262.     VScale(T, (*Control_Points)[i][j], t);
  263.     VAddEq(*P, T);
  264.     t = c * vt * (du[i] * uu[3-i] + u[i] * duu[3-i]);
  265.     VScale(T, (*Control_Points)[i][j], t);
  266.     VAddEq(U, T);
  267.     t = c * ut * (dv[j] * vv[3-j] + v[j] * dvv[3-j]);
  268.     VScale(T, (*Control_Points)[i][j], t);
  269.     VAddEq(V, T);
  270.     }
  271.  
  272.   /* Make the normal from the cross product of the tangents */
  273.   VCross(*N, U, V);
  274.   VDot(t, *N, *N);
  275.   if (t > BEZIER_EPSILON) 
  276.     {
  277.     t = 1.0 / sqrt(t);
  278.     VScaleEq(*N, t);
  279.     }
  280.   else
  281.     Make_Vector(N, 1, 0, 0)
  282.       }
  283.  
  284. /* Calculate the normal to a subpatch (triangle) return the vector
  285. <1.0 0.0 0.0> if the triangle is degenerate. */
  286. static int subpatch_normal(v1, v2, v3, Result, d)
  287. VECTOR *v1, *v2, *v3, *Result;
  288. DBL *d;
  289.   {
  290.   VECTOR V1, V2;
  291.   DBL Length;
  292.  
  293.   VSub(V1, *v1, *v2);
  294.   VSub(V2, *v3, *v2);
  295.   VCross(*Result, V1, V2);
  296.   VLength(Length, *Result);
  297.   if (Length < BEZIER_EPSILON) 
  298.     {
  299.     Result->x = 1.0;
  300.     Result->y = 0.0;
  301.     Result->z = 0.0;
  302.     *d = -1.0 * v1->x;
  303.     return 0;
  304.     }
  305.   else 
  306.     {
  307.     VInverseScale(*Result, *Result, Length);
  308.     VDot(*d, *Result, *v1);
  309.     *d = 0.0 - *d;
  310.     return 1;
  311.     }
  312.   }
  313.  
  314. static int
  315. InvertMatrix(in, out)
  316. VECTOR in[3], out[3];
  317.   {
  318.   int i;
  319.   DBL det;
  320.  
  321.   out[0].x =  (in[1].y * in[2].z - in[1].z * in[2].y);
  322.   out[1].x = -(in[0].y * in[2].z - in[0].z * in[2].y);
  323.   out[2].x =  (in[0].y * in[1].z - in[0].z * in[1].y);
  324.  
  325.   out[0].y = -(in[1].x * in[2].z - in[1].z * in[2].x);
  326.   out[1].y =  (in[0].x * in[2].z - in[0].z * in[2].x);
  327.   out[2].y = -(in[0].x * in[1].z - in[0].z * in[1].x);
  328.  
  329.   out[0].z =  (in[1].x * in[2].y - in[1].y * in[2].x);
  330.   out[1].z = -(in[0].x * in[2].y - in[0].y * in[2].x);
  331.   out[2].z =  (in[0].x * in[1].y - in[0].y * in[1].x);
  332.  
  333.   det = 
  334.   in[0].x * in[1].y * in[2].z +
  335.   in[0].y * in[1].z * in[2].x +
  336.   in[0].z * in[1].x * in[2].y -
  337.   in[0].z * in[1].y * in[2].x -
  338.   in[0].x * in[1].z * in[2].y -
  339.   in[0].y * in[1].x * in[2].z;
  340.  
  341.   if (det > -BEZIER_EPSILON && det < BEZIER_EPSILON)
  342.     return 0;
  343.  
  344.   det = 1.0 / det;
  345.  
  346.   for (i=0;i<3;i++)
  347.     VScaleEq(out[i], det)
  348.  
  349.       return 1;
  350.   }
  351.  
  352. static int
  353. intersect_subpatch(Shape, ray, V, uu, vv, Depth, P, N, u, v)
  354. BICUBIC_PATCH *Shape;
  355. RAY *ray;
  356. VECTOR V[3];
  357. DBL uu[3], vv[3], *Depth;
  358. VECTOR *P, *N;
  359. DBL *u, *v;
  360.   {
  361.   VECTOR B[3], IB[3], NN[3], Q, T;
  362.   DBL d, n, a, b, r;
  363.  
  364.   VSub(B[0], V[1], V[0]);
  365.   VSub(B[1], V[2], V[0]);
  366.   VCross(B[2], B[0], B[1]);
  367.   VDot(d, B[2], B[2]);
  368.   if (d < BEZIER_EPSILON)
  369.     return 0;
  370.   d = 1.0 / sqrt(d);
  371.   VScaleEq(B[2], d);
  372.  
  373.   if (!InvertMatrix(B, IB))
  374.     /* Degenerate triangle */
  375.     return 0;
  376.  
  377.   VDot(d, ray->Direction, IB[2]);
  378.   if (d > -BEZIER_EPSILON && d < BEZIER_EPSILON)
  379.     return 0;
  380.  
  381.   VSub(Q, V[0], ray->Initial);
  382.   VDot(n, Q, IB[2]);
  383.   *Depth = n / d;
  384.   if (*Depth < BEZIER_TOLERANCE)
  385.     return 0;
  386.   VScale(T, ray->Direction, *Depth);
  387.   VAdd(*P, ray->Initial, T);
  388.   VSub(Q, *P, V[0]);
  389.  
  390.   VDot(a, Q, IB[0]);
  391.   VDot(b, Q, IB[1]);
  392.   if (a < 0.0 || b < 0.0 || a + b > 1.0) 
  393.     return 0;
  394.  
  395.   r = 1.0 - a - b;
  396.  
  397.   Make_Vector(N, 0.0, 0.0, 0.0);
  398.   bezier_value((VECTOR (*)[4][4])&Shape->Control_Points,
  399.     uu[0], vv[0], &T, &NN[0]);
  400.   bezier_value((VECTOR (*)[4][4])&Shape->Control_Points,
  401.     uu[1], vv[1], &T, &NN[1]);
  402.   bezier_value((VECTOR (*)[4][4])&Shape->Control_Points,
  403.     uu[2], vv[2], &T, &NN[2]);
  404.   VScale(T, NN[0], r); VAdd(*N, *N, T);
  405.   VScale(T, NN[1], a); VAdd(*N, *N, T);
  406.   VScale(T, NN[2], b); VAdd(*N, *N, T);
  407.   *u = r * uu[0] + a * uu[1] + b * uu[2];
  408.   *v = r * vv[0] + a * vv[1] + b * vv[2];
  409.   VDot(d, *N, *N);
  410.   if (d > BEZIER_EPSILON) 
  411.     {
  412.     d = 1.0 / sqrt(d);
  413.     VScaleEq(*N, d);
  414.     }
  415.   else
  416.     Make_Vector(N, 1, 0, 0);
  417.   return 1;
  418.   }
  419.  
  420. /* Find a sphere that contains all of the points in the list "vectors" */
  421. static void find_average(vector_count, vectors, center, radius)
  422. int vector_count;
  423. VECTOR *vectors, *center;
  424. DBL *radius;
  425.   {
  426.   DBL r0, r1, xc=0, yc=0, zc=0;
  427.   DBL x0, y0, z0;
  428.   int i;
  429.   for (i=0;i<vector_count;i++) 
  430.     {
  431.     xc += vectors[i].x;
  432.     yc += vectors[i].y;
  433.     zc += vectors[i].z;
  434.     }
  435.   xc /= (DBL)vector_count;
  436.   yc /= (DBL)vector_count;
  437.   zc /= (DBL)vector_count;
  438.   r0 = 0.0;
  439.   for (i=0;i<vector_count;i++) 
  440.     {
  441.     x0 = vectors[i].x - xc;
  442.     y0 = vectors[i].y - yc;
  443.     z0 = vectors[i].z - zc;
  444.     r1 = x0 * x0 + y0 * y0 + z0 * z0;
  445.     if (r1 > r0) r0 = r1;
  446.     }
  447.   center->x = xc; center->y = yc; center->z = zc;
  448.   *radius = r0;
  449.   }
  450.  
  451. static int spherical_bounds_check(Ray, center, radius)
  452. RAY *Ray;
  453. VECTOR *center;
  454. DBL radius;
  455.   {
  456.   DBL x, y, z, dist1, dist2;
  457.   x = center->x - Ray->Initial.x;
  458.   y = center->y - Ray->Initial.y;
  459.   z = center->z - Ray->Initial.z;
  460.   dist1 = x * x + y * y + z * z;
  461.   if (dist1 < radius)
  462.     /* ray starts inside sphere - assume it intersects. */
  463.     return 1;
  464.   else 
  465.     {
  466.     dist2 = x*Ray->Direction.x + y*Ray->Direction.y + z*Ray->Direction.z;
  467.     dist2 = dist2 * dist2;
  468.     if (dist2 > 0 && (dist1 - dist2 < radius))
  469.       return 1;
  470.     }
  471.   return 0;
  472.   }
  473.  
  474. /* Find a sphere that bounds all of the control points of a Bezier patch.
  475.    The values returned are: the center of the bounding sphere, and the
  476.    square of the radius of the bounding sphere. */
  477. static void bezier_bounding_sphere(Patch, center, radius)
  478. VECTOR (*Patch)[4][4], *center;
  479. DBL *radius;
  480.   {
  481.   DBL r0, r1, xc=0, yc=0, zc=0;
  482.   DBL x0, y0, z0;
  483.   int i, j;
  484.   for (i=0;i<4;i++) 
  485.     {
  486.     for (j=0;j<4;j++) 
  487.       {
  488.       xc += (*Patch)[i][j].x;
  489.       yc += (*Patch)[i][j].y;
  490.       zc += (*Patch)[i][j].z;
  491.       }
  492.     }
  493.   xc /= 16.0;
  494.   yc /= 16.0;
  495.   zc /= 16.0;
  496.   r0 = 0.0;
  497.   for (i=0;i<4;i++) 
  498.     {
  499.     for (j=0;j<4;j++) 
  500.       {
  501.       x0 = (*Patch)[i][j].x - xc;
  502.       y0 = (*Patch)[i][j].y - yc;
  503.       z0 = (*Patch)[i][j].z - zc;
  504.       r1 = x0 * x0 + y0 * y0 + z0 * z0;
  505.       if (r1 > r0) r0 = r1;
  506.       }
  507.     }
  508.   center->x = xc; center->y = yc; center->z = zc;
  509.   *radius = r0;
  510.   }
  511.  
  512. /* Precompute grid points and normals for a bezier patch */
  513. void Precompute_Patch_Values(Shape)
  514. BICUBIC_PATCH *Shape;
  515.   {
  516.   int i, j;
  517.   VECTOR Control_Points[16];
  518.   VECTOR (*Patch_Ptr)[4][4] = (VECTOR (*)[4][4]) Shape->Control_Points;
  519.  
  520.   /* Calculate the bounding sphere for the entire patch. */
  521.   for (i=0;i<4;i++) for (j=0;j<4;j++)
  522.     Control_Points[4*i+j] = Shape->Control_Points[i][j];
  523.   find_average(16, &Control_Points[0], &Shape->Bounding_Sphere_Center,
  524.     &Shape->Bounding_Sphere_Radius);
  525.  
  526.   if (Shape->Patch_Type == 0)
  527.     return;
  528.   else 
  529.     {
  530.     if (Shape->Node_Tree != NULL)
  531.       bezier_tree_deleter(Shape->Node_Tree);
  532.     Shape->Node_Tree = bezier_tree_builder(Shape, Patch_Ptr,
  533.       0.0, 1.0, 0.0, 1.0, 0);
  534.     return;
  535.     }
  536.   }
  537.  
  538. /* Determine the distance from a point to a plane. */
  539. static DBL point_plane_distance(p, n, d)
  540. VECTOR *p, *n;
  541. DBL *d;
  542.   {
  543.   DBL temp1, temp2;
  544.  
  545.   VDot(temp1, *p, *n);
  546.   temp1 += *d;
  547.   VLength(temp2, *n);
  548.   if (fabs(temp2) < EPSILON) return 0;
  549.   temp1 /= temp2;
  550.   return temp1;
  551.   }
  552.  
  553. static void
  554. bezier_subpatch_intersect(ray, Shape, Patch, u0, u1, v0, v1,
  555. depth_count, Depths)
  556. RAY *ray;
  557. BICUBIC_PATCH *Shape;
  558. VECTOR (*Patch)[4][4];
  559. DBL u0, u1, v0, v1;
  560. int *depth_count;
  561. DBL *Depths;
  562.   {
  563.   int tcnt = Shape->Intersection_Count;
  564.   VECTOR V[3];
  565.   DBL u, v, Depth, uu[3], vv[3];
  566.   VECTOR P, N;
  567.  
  568.   if (tcnt + *depth_count >= MAX_BICUBIC_INTERSECTIONS) return;
  569.   V[0] = (*Patch)[0][0];
  570.   V[1] = (*Patch)[0][3];
  571.   V[2] = (*Patch)[3][0];
  572.  
  573.   uu[0] = u0; uu[1] = u0; uu[2] = u1;
  574.   vv[0] = v0; vv[1] = v1; vv[2] = v1;
  575.  
  576.   if (intersect_subpatch(Shape, ray, V, uu, vv, &Depth, &P, &N, &u, &v)) 
  577.     {
  578.     Shape->IPoint[tcnt + *depth_count] = P;
  579.     Shape->Normal_Vector[tcnt + *depth_count] = N;
  580.     Depths[*depth_count] = Depth;
  581.     *depth_count += 1;
  582.     }
  583.  
  584.   if (tcnt + *depth_count >= MAX_BICUBIC_INTERSECTIONS) return;
  585.  
  586.   V[1] = V[2];
  587.   V[2] = (*Patch)[3][0];
  588.   uu[1] = uu[2]; uu[2] = u1;
  589.   vv[1] = vv[2]; vv[2] = v0;
  590.  
  591.   if (intersect_subpatch(Shape, ray, V, uu, vv, &Depth, &P, &N, &u, &v)) 
  592.     {
  593.     Shape->IPoint[tcnt + *depth_count] = P;
  594.     Shape->Normal_Vector[tcnt + *depth_count] = N;
  595.     Depths[*depth_count] = Depth;
  596.     *depth_count += 1;
  597.     }
  598.   }
  599.  
  600.  
  601. static void bezier_split_left_right(Patch, Left_Patch, Right_Patch)
  602. VECTOR (*Patch)[4][4], (*Left_Patch)[4][4], (*Right_Patch)[4][4];
  603.   {
  604.   VECTOR Temp1[4], Temp2[4], Half;
  605.   int i, j;
  606.   for (i=0;i<4;i++) 
  607.     {
  608.     Temp1[0] = (*Patch)[0][i];
  609.     VHalf(Temp1[1], (*Patch)[0][i], (*Patch)[1][i]);
  610.     VHalf(Half, (*Patch)[1][i], (*Patch)[2][i]);
  611.     VHalf(Temp1[2], Temp1[1], Half);
  612.     VHalf(Temp2[2], (*Patch)[2][i], (*Patch)[3][i]);
  613.     VHalf(Temp2[1], Half, Temp2[2]);
  614.     VHalf(Temp1[3], Temp1[2], Temp2[1]);
  615.     Temp2[0] = Temp1[3];
  616.     Temp2[3] = (*Patch)[3][i];
  617.     for (j=0;j<4;j++) 
  618.       {
  619.       (*Left_Patch)[j][i]  = Temp1[j];
  620.       (*Right_Patch)[j][i] = Temp2[j];
  621.       }
  622.     }
  623.   }
  624.  
  625. static void bezier_split_up_down(Patch, Bottom_Patch, Top_Patch)
  626. VECTOR (*Patch)[4][4], (*Top_Patch)[4][4], (*Bottom_Patch)[4][4];
  627.   {
  628.   VECTOR Temp1[4], Temp2[4], Half;
  629.   int i, j;
  630.  
  631.   for (i=0;i<4;i++) 
  632.     {
  633.     Temp1[0] = (*Patch)[i][0];
  634.     VHalf(Temp1[1], (*Patch)[i][0], (*Patch)[i][1]);
  635.     VHalf(Half, (*Patch)[i][1], (*Patch)[i][2]);
  636.     VHalf(Temp1[2], Temp1[1], Half);
  637.     VHalf(Temp2[2], (*Patch)[i][2], (*Patch)[i][3]);
  638.     VHalf(Temp2[1], Half, Temp2[2]);
  639.     VHalf(Temp1[3], Temp1[2], Temp2[1]);
  640.     Temp2[0] = Temp1[3];
  641.     Temp2[3] = (*Patch)[i][3];
  642.     for (j=0;j<4;j++) 
  643.       {
  644.       (*Bottom_Patch)[i][j] = Temp1[j];
  645.       (*Top_Patch)[i][j]    = Temp2[j];
  646.       }
  647.     }
  648.   }
  649.  
  650. /* See how close to a plane a subpatch is, the patch must have at least
  651.    three distinct vertices. A negative result from this function indicates
  652.    that a degenerate value of some sort was encountered. */
  653. static DBL determine_subpatch_flatness(Patch)
  654. VECTOR (*Patch)[4][4];
  655.   {
  656.   VECTOR vertices[4], n, TempV;
  657.   DBL d, dist, temp1;
  658.   int i, j;
  659.  
  660.   vertices[0] = (*Patch)[0][0];
  661.   vertices[1] = (*Patch)[0][3];
  662.   VSub(TempV, vertices[0], vertices[1]);
  663.   VLength(temp1, TempV);
  664.   if (fabs(temp1) < EPSILON) 
  665.     {
  666.     /* Degenerate in the V direction for U = 0. This is ok if the other
  667.          two corners are distinct from the lower left corner - I'm sure there
  668.          are cases where the corners coincide and the middle has good values,
  669.          but that is somewhat pathalogical and won't be considered. */
  670.     vertices[1] = (*Patch)[3][3];
  671.     VSub(TempV, vertices[0], vertices[1]);
  672.     VLength(temp1, TempV);
  673.     if (fabs(temp1) < EPSILON) return -1.0;
  674.     vertices[2] = (*Patch)[3][0];
  675.     VSub(TempV, vertices[0], vertices[1]);
  676.     VLength(temp1, TempV);
  677.     if (fabs(temp1) < EPSILON) return -1.0;
  678.     VSub(TempV, vertices[1], vertices[2]);
  679.     VLength(temp1, TempV);
  680.     if (fabs(temp1) < EPSILON) return -1.0;
  681.     }
  682.   else 
  683.     {
  684.     vertices[2] = (*Patch)[3][0];
  685.     VSub(TempV, vertices[0], vertices[1]);
  686.     VLength(temp1, TempV);
  687.     if (fabs(temp1) < EPSILON) 
  688.       {
  689.       vertices[2] = (*Patch)[3][3];
  690.       VSub(TempV, vertices[0], vertices[2]);
  691.       VLength(temp1, TempV);
  692.       if (fabs(temp1) < EPSILON)
  693.         return -1.0;
  694.       VSub(TempV, vertices[1], vertices[2]);
  695.       VLength(temp1, TempV);
  696.       if (fabs(temp1) < EPSILON)
  697.         return -1.0;
  698.       }
  699.     else 
  700.       {
  701.       VSub(TempV, vertices[1], vertices[2]);
  702.       VLength(temp1, TempV);
  703.       if (fabs(temp1) < EPSILON)
  704.         return -1.0;
  705.       }
  706.     }
  707.   /* Now that a good set of candidate points has been found, find the
  708.       plane equations for the patch */
  709.   if (subpatch_normal(&vertices[0], &vertices[1], &vertices[2], &n, &d)) 
  710.     {
  711.     /* Step through all vertices and see what the maximum distance from the
  712.              plane happens to be. */
  713.     dist = 0.0;
  714.     for (i=0;i<4;i++) 
  715.       {
  716.       for (j=0;j<4;j++) 
  717.         {
  718.         temp1 = fabs(point_plane_distance(&((*Patch)[i][j]), &n, &d));
  719.         if (temp1 > dist)
  720.           dist = temp1;
  721.         }
  722.       }
  723.     return dist;
  724.     }
  725.   else 
  726.     {
  727.     if (Options & DEBUGGING) 
  728.       {
  729.       printf("Subpatch normal failed in determine_subpatch_flatness\n");
  730.       fflush(stdout);
  731.       }
  732.     return -1.0;
  733.     }
  734.   }
  735.  
  736. static int flat_enough(Object, Patch)
  737. BICUBIC_PATCH *Object;
  738. VECTOR (*Patch)[4][4];
  739.   {
  740.   DBL Dist;
  741.  
  742.   Dist = determine_subpatch_flatness(Patch);
  743.   if (Dist < 0.0)
  744.     return 0;
  745.   else if (Dist < Object->Flatness_Value)
  746.     return 1;
  747.   else
  748.     return 0;
  749.   }
  750.  
  751. static void bezier_subdivider(Ray, Object, Patch, u0, u1, v0, v1,
  752. recursion_depth, depth_count, Depths)
  753. RAY *Ray;
  754. BICUBIC_PATCH *Object;
  755. VECTOR (*Patch)[4][4];
  756. DBL u0, u1, v0, v1;
  757. int recursion_depth, *depth_count;
  758. DBL *Depths;
  759.   {
  760.   VECTOR Lower_Left[4][4], Lower_Right[4][4];
  761.   VECTOR Upper_Left[4][4], Upper_Right[4][4];
  762.   VECTOR center;
  763.   DBL ut, vt, radius;
  764.   int tcnt = Object->Intersection_Count;
  765.  
  766.   /* Don't waste time if there are already too many intersections */
  767.   if (tcnt >= MAX_BICUBIC_INTERSECTIONS) return;
  768.  
  769.   /* Make sure the ray passes through a sphere bounding the control points of
  770.       the patch */
  771.   bezier_bounding_sphere(Patch, ¢er, &radius);
  772.   if (!spherical_bounds_check(Ray, ¢er, radius))
  773.     return;
  774.  
  775.   /* If the patch is close to being flat, then just perform a ray-plane
  776.       intersection test. */
  777.   if (flat_enough(Object, Patch))
  778.     bezier_subpatch_intersect(Ray, Object, Patch, u0, u1, v0, v1,
  779.       depth_count, Depths);
  780.  
  781.   if (recursion_depth >= Object->U_Steps)
  782.     if (recursion_depth >= Object->V_Steps)
  783.       bezier_subpatch_intersect(Ray, Object, Patch, u0, u1, v0, v1,
  784.         depth_count, Depths);
  785.     else 
  786.       {
  787.       bezier_split_up_down(Patch, (VECTOR (*)[4][4])Lower_Left,
  788.         (VECTOR (*)[4][4])Upper_Left);
  789.       vt = (v1 - v0) / 2.0;
  790.       bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Lower_Left,
  791.         u0, u1, v0, vt,
  792.         recursion_depth+1, depth_count, Depths);
  793.       bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Upper_Left,
  794.         u0, u1, vt, v1,
  795.         recursion_depth+1, depth_count, Depths);
  796.       }
  797.   else if (recursion_depth >= Object->V_Steps) 
  798.     {
  799.     bezier_split_left_right(Patch, (VECTOR (*)[4][4])Lower_Left,
  800.       (VECTOR (*)[4][4])Lower_Right);
  801.     ut = (u1 - u0) / 2.0;
  802.     bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Lower_Left,
  803.       u0, ut, v0, v1,
  804.       recursion_depth+1, depth_count, Depths);
  805.     bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Lower_Right,
  806.       ut, u1, v0, v1,
  807.       recursion_depth+1, depth_count, Depths);
  808.     }
  809.   else 
  810.     {
  811.     ut = (u1 - u0) / 2.0;
  812.     vt = (v1 - v0) / 2.0;
  813.     bezier_split_left_right(Patch, (VECTOR (*)[4][4])Lower_Left,
  814.       (VECTOR (*)[4][4])Lower_Right);
  815.     bezier_split_up_down((VECTOR (*)[4][4])Lower_Left,
  816.       (VECTOR (*)[4][4])Lower_Left,
  817.       (VECTOR (*)[4][4])Upper_Left);
  818.     bezier_split_up_down((VECTOR (*)[4][4])Lower_Right,
  819.       (VECTOR (*)[4][4])Lower_Right,
  820.       (VECTOR (*)[4][4])Upper_Right);
  821.     bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Lower_Left,
  822.       u0, ut, v0, vt,
  823.       recursion_depth+1, depth_count, Depths);
  824.     bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Upper_Left,
  825.       u0, ut, vt, v1,
  826.       recursion_depth+1, depth_count, Depths);
  827.     bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Lower_Right,
  828.       ut, u1, v0, vt,
  829.       recursion_depth+1, depth_count, Depths);
  830.     bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Upper_Right,
  831.       ut, u1, vt, v1,
  832.       recursion_depth+1, depth_count, Depths);
  833.     }
  834.   }
  835.  
  836. static void bezier_tree_deleter(Node)
  837. BEZIER_NODE *Node;
  838.   {
  839.   BEZIER_CHILDREN *Children;
  840.   int i;
  841.  
  842.   /* If this is an interior node then continue the descent */
  843.   if (Node->Node_Type == BEZIER_INTERIOR_NODE) 
  844.     {
  845.     Children = (BEZIER_CHILDREN *)Node->Data_Ptr;
  846.     for (i=0;i<Node->Count;i++)
  847.       bezier_tree_deleter(Children->Children[i]);
  848.     free((void *)Children);
  849.     }
  850.   else if (Node->Node_Type == BEZIER_LEAF_NODE) 
  851.     {
  852.     /* Free the memory used for the vertices. */
  853.     free((void *)Node->Data_Ptr);
  854.     }
  855.   /* Free the memory used for the node. */
  856.   free((void *)Node);
  857.   }
  858.  
  859. static void bezier_tree_walker(Ray, Shape, Node, depth, depth_count, Depths)
  860. RAY *Ray;
  861. BICUBIC_PATCH *Shape;
  862. BEZIER_NODE *Node;
  863. int depth, *depth_count;
  864. DBL *Depths;
  865.   {
  866.   BEZIER_CHILDREN *Children;
  867.   BEZIER_VERTICES *Vertices;
  868.   VECTOR N, P, V[3];
  869.   DBL Depth, u, v, uu[3], vv[3];
  870.   int i, tcnt = Shape->Intersection_Count;
  871.  
  872.   /* Don't waste time if there are already too many intersections */
  873.   if (tcnt >= MAX_BICUBIC_INTERSECTIONS) return;
  874.  
  875.   /* Make sure the ray passes through a sphere bounding the control points of
  876.       the patch */
  877.   if (!spherical_bounds_check(Ray, &(Node->Center), Node->Radius_Squared))
  878.     return;
  879.  
  880.   /* If this is an interior node then continue the descent, else
  881.       do a check against the vertices. */
  882.   if (Node->Node_Type == BEZIER_INTERIOR_NODE) 
  883.     {
  884.     Children = (BEZIER_CHILDREN *)Node->Data_Ptr;
  885.     for (i=0;i<Node->Count;i++)
  886.       bezier_tree_walker(Ray, Shape, Children->Children[i],
  887.         depth+1, depth_count, Depths);
  888.     }
  889.   else if (Node->Node_Type == BEZIER_LEAF_NODE) 
  890.     {
  891.     Vertices = (BEZIER_VERTICES *)Node->Data_Ptr;
  892.     V[0] = Vertices->Vertices[0];
  893.     V[1] = Vertices->Vertices[1];
  894.     V[2] = Vertices->Vertices[2];
  895.  
  896.     uu[0] = Vertices->uvbnds[0];
  897.     uu[1] = Vertices->uvbnds[0];
  898.     uu[2] = Vertices->uvbnds[1];
  899.     vv[0] = Vertices->uvbnds[2];
  900.     vv[1] = Vertices->uvbnds[3];
  901.     vv[2] = Vertices->uvbnds[3];
  902.  
  903.     /* Triangulate this subpatch, then check for intersections in
  904.          the triangles. */
  905.     if (intersect_subpatch(Shape, Ray, V, uu, vv, &Depth, &P, &N, &u, &v)) 
  906.       {
  907.       Shape->IPoint[tcnt + *depth_count] = P;
  908.       Shape->Normal_Vector[tcnt + *depth_count] = N;
  909.       Depths[*depth_count] = Depth;
  910.       *depth_count += 1;
  911.       }
  912.     if (*depth_count + tcnt >= MAX_BICUBIC_INTERSECTIONS) return;
  913.  
  914.     V[1] = V[2];
  915.     V[2] = Vertices->Vertices[3];
  916.     uu[1] = uu[2]; uu[2] = Vertices->uvbnds[1];
  917.     vv[1] = vv[2]; vv[2] = Vertices->uvbnds[2];
  918.  
  919.     if (intersect_subpatch(Shape, Ray, V, uu, vv, &Depth, &P, &N, &u, &v)) 
  920.       {
  921.       Shape->IPoint[tcnt + *depth_count] = P;
  922.       Shape->Normal_Vector[tcnt + *depth_count] = N;
  923.       Depths[*depth_count] = Depth;
  924.       *depth_count += 1;
  925.       }
  926.     }
  927.   else 
  928.     {
  929.     printf("Bad Node type at depth %d\n", depth);
  930.     }
  931.   }
  932.  
  933. static int intersect_bicubic_patch0(Ray, Shape, Depths)
  934. RAY *Ray;
  935. BICUBIC_PATCH *Shape;
  936. DBL *Depths;
  937.   {
  938.   int cnt = 0;
  939.   VECTOR (*Patch)[4][4] = (VECTOR (*)[4][4]) Shape->Control_Points;
  940.  
  941.   bezier_subdivider(Ray, Shape, Patch, 0.0, 1.0, 0.0, 1.0,
  942.     0, &cnt, Depths);
  943.   return cnt;
  944.   }
  945.  
  946. static int intersect_bicubic_patch1(Ray, Shape, Depths)
  947. RAY *Ray;
  948. BICUBIC_PATCH *Shape;
  949. DBL *Depths;
  950.   {
  951.   int cnt = 0;
  952.   bezier_tree_walker(Ray, Shape, Shape->Node_Tree, 0, &cnt, Depths);
  953.   return cnt;
  954.   }
  955.  
  956. int All_Bicubic_Patch_Intersections(Object, Ray, Depth_Stack)
  957. OBJECT *Object;
  958. RAY *Ray;
  959. ISTACK *Depth_Stack;
  960.   {
  961.   DBL Depths[MAX_BICUBIC_INTERSECTIONS];
  962.   VECTOR IPoint;
  963.   int cnt, tcnt, i, Found;
  964.  
  965.   Found = FALSE;
  966.   Ray_Bicubic_Tests++;
  967.  
  968.   if (Ray == CM_Ray)
  969.     ((BICUBIC_PATCH *)Object)->Intersection_Count = 0;
  970.  
  971.   tcnt = ((BICUBIC_PATCH *)Object)->Intersection_Count;
  972.  
  973.   switch (((BICUBIC_PATCH *)Object)->Patch_Type)
  974.   {
  975.   case 0: 
  976.     cnt = intersect_bicubic_patch0(Ray, ((BICUBIC_PATCH *)Object), &Depths[0]);
  977.     break;
  978.   case 1: 
  979.     cnt = intersect_bicubic_patch1(Ray, ((BICUBIC_PATCH *)Object), &Depths[0]);
  980.     break;
  981.   default: 
  982.     Error("Bad patch type\n");
  983.   }
  984.  
  985.   if (cnt > 0) Ray_Bicubic_Tests_Succeeded++;
  986.   for (i=0;i<cnt;i++) 
  987.     {
  988.     if (!Shadow_Test_Flag)
  989.       ((BICUBIC_PATCH *)Object)->Intersection_Count++;
  990.     IPoint = ((BICUBIC_PATCH *)Object)->IPoint[tcnt + i];
  991.     if (Point_In_Clip(&IPoint,Object->Clip))
  992.       {
  993.       push_entry(Depths[i], IPoint, Object, Depth_Stack);
  994.       Found = TRUE;
  995.       }
  996.     }
  997.   return (Found);
  998.   }
  999.  
  1000. /* A patch is not a solid, so an inside test doesn't make sense. */
  1001. int Inside_Bicubic_Patch (IPoint, Object)
  1002. VECTOR *IPoint;
  1003. OBJECT *Object;
  1004.   {
  1005.   return 0;
  1006.   }
  1007.  
  1008. void Bicubic_Patch_Normal (Result, Object, IPoint)
  1009. OBJECT *Object;
  1010. VECTOR *Result, *IPoint;
  1011.   {
  1012.   BICUBIC_PATCH *Patch = (BICUBIC_PATCH *)Object;
  1013.   int i;
  1014.  
  1015.   /* If all is going well, the normal was computed at the time the intersection
  1016.       was computed.  Look on the list of associated intersection points and normals */
  1017.   for (i=0;i<Patch->Intersection_Count;i++)
  1018.     if (IPoint->x == Patch->IPoint[i].x &&
  1019.       IPoint->y == Patch->IPoint[i].y &&
  1020.       IPoint->z == Patch->IPoint[i].z) 
  1021.       {
  1022.       Result->x = Patch->Normal_Vector[i].x;
  1023.       Result->y = Patch->Normal_Vector[i].y;
  1024.       Result->z = Patch->Normal_Vector[i].z;
  1025.       return;
  1026.       }
  1027.   if (Options & DEBUGGING) 
  1028.     {
  1029.     printf("Bicubic patch normal for unknown intersection point\n");
  1030.     fflush(stdout);
  1031.     }
  1032.   Result->x = 1.0;
  1033.   Result->y = 0.0;
  1034.   Result->z = 0.0;
  1035.   }
  1036.  
  1037. void Translate_Bicubic_Patch (Object, Vector)
  1038. OBJECT *Object;
  1039. VECTOR *Vector;
  1040.   {
  1041.   BICUBIC_PATCH *Patch = (BICUBIC_PATCH *) Object;
  1042.   int i, j;
  1043.   for (i=0;i<4;i++) for (j=0;j<4;j++)
  1044.     VAdd (Patch->Control_Points[i][j], Patch->Control_Points[i][j], *Vector)
  1045.       Precompute_Patch_Values(Patch);
  1046.   }
  1047.  
  1048. void Rotate_Bicubic_Patch (Object, Vector)
  1049. OBJECT *Object;
  1050. VECTOR *Vector;
  1051.   {
  1052.   TRANSFORM Trans;
  1053.  
  1054.   Compute_Rotation_Transform (&Trans, Vector);
  1055.   Transform_Bicubic_Patch (Object,&Trans);
  1056.   }
  1057.  
  1058. void Scale_Bicubic_Patch (Object, Vector)
  1059. OBJECT *Object;
  1060. VECTOR *Vector;
  1061.   {
  1062.   BICUBIC_PATCH *Patch = (BICUBIC_PATCH *) Object;
  1063.   int i, j;
  1064.   for (i=0;i<4;i++) for (j=0;j<4;j++)
  1065.     VEvaluate(Patch->Control_Points[i][j],
  1066.       Patch->Control_Points[i][j], *Vector);
  1067.   Precompute_Patch_Values(Patch);
  1068.   }
  1069.  
  1070.  
  1071. /* Inversion of a patch really doesn't make sense. */
  1072. void Invert_Bicubic_Patch (Object)
  1073. OBJECT *Object;
  1074.   {
  1075.   ;
  1076.   }
  1077.  
  1078. BICUBIC_PATCH *Create_Bicubic_Patch()
  1079.   {
  1080.   BICUBIC_PATCH *New;
  1081.  
  1082.   if ((New = (BICUBIC_PATCH *) malloc (sizeof (BICUBIC_PATCH))) == NULL)
  1083.     MAError ("bicubic patch");
  1084.  
  1085.   INIT_OBJECT_FIELDS(New,BICUBIC_PATCH_OBJECT,&Bicubic_Patch_Methods)
  1086.  
  1087.     New->Patch_Type = -1;
  1088.   New->U_Steps = 0;
  1089.   New->V_Steps = 0;
  1090.   New->Intersection_Count = 0;
  1091.   New->Flatness_Value = 0.0;
  1092.   New->Node_Tree = NULL;
  1093.  
  1094.   /* NOTE: Control_Points[4][4] is initialized in Parse_Bicubic_Patch.
  1095.    Bounding_Sphere_Center,Bounding_Sphere_Radius, Normal_Vector[], and 
  1096.    IPoint[] are initialized in Precompute_Patch_Values. 
  1097.  */
  1098.  
  1099.   return (New);
  1100.   }
  1101.  
  1102. void *Copy_Bicubic_Patch (Object)
  1103. OBJECT *Object;
  1104.   {
  1105.   BICUBIC_PATCH *New;
  1106.   int i,j;
  1107.  
  1108.   New = Create_Bicubic_Patch ();
  1109.  
  1110.   /*  Do not do *New = *Old so that Precompute works right */
  1111.  
  1112.   New->Patch_Type = ((BICUBIC_PATCH *)Object)->Patch_Type;
  1113.   New->U_Steps    = ((BICUBIC_PATCH *)Object)->U_Steps;
  1114.   New->V_Steps    = ((BICUBIC_PATCH *)Object)->V_Steps;
  1115.   for (i=0;i<4;i++) for (j=0;j<4;j++)
  1116.     New->Control_Points[i][j]
  1117.     = ((BICUBIC_PATCH *)Object)->Control_Points[i][j];
  1118.  
  1119.   New->Flatness_Value     = ((BICUBIC_PATCH *)Object)->Flatness_Value;
  1120.   New->Intersection_Count = ((BICUBIC_PATCH *)Object)->Intersection_Count;
  1121.  
  1122.   Precompute_Patch_Values(New);
  1123.  
  1124.   return (New);
  1125.   }
  1126.  
  1127. void Transform_Bicubic_Patch (Object, Trans)
  1128. OBJECT *Object;
  1129. TRANSFORM *Trans;
  1130.   {
  1131.   BICUBIC_PATCH *Patch = (BICUBIC_PATCH *) Object;
  1132.   int i, j;
  1133.  
  1134.   for (i=0;i<4;i++) for (j=0;j<4;j++)
  1135.     MTransPoint(&(Patch->Control_Points[i][j]),
  1136.       &(Patch->Control_Points[i][j]),
  1137.       Trans);
  1138.   Precompute_Patch_Values(Patch);
  1139.   }
  1140.  
  1141. void Destroy_Bicubic_Patch (Object)
  1142. OBJECT *Object;
  1143.   {
  1144.  
  1145.   /* Need to add code to free() all malloc'ed memory created for this
  1146.      object. */
  1147.   free (Object);     
  1148.   }
  1149.