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