home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / mac / developm / source / povsrc.sit / SOURCE / BEZIER.C next >
Encoding:
C/C++ Source or Header  |  1992-07-03  |  49.6 KB  |  1,481 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 1992 Persistence of Vision Team
  11. *---------------------------------------------------------------------------
  12. *  Copying, distribution and legal info is in the file povlegal.doc which
  13. *  should be distributed with this file. If povlegal.doc is not available
  14. *  or for more info please contact:
  15. *
  16. *       Drew Wells [POV-Team Leader] 
  17. *       CIS: 73767,1244  Internet: 73767.1244@compuserve.com
  18. *       Phone: (213) 254-4041
  19. * 
  20. * This program is based on the popular DKB raytracer version 2.12.
  21. * DKBTrace was originally written by David K. Buck.
  22. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  23. *
  24. *****************************************************************************/
  25.  
  26. #include "frame.h"
  27. #include "vector.h"
  28. #include "povproto.h"
  29.  
  30. #undef EPSILON
  31. #ifndef EPSILON
  32. #define EPSILON 1.0e-10
  33. #endif
  34.  
  35. METHODS Bicubic_Patch_Methods =
  36. { Object_Intersect, All_Bicubic_Patch_Intersections,
  37.    Inside_Bicubic_Patch, Bicubic_Patch_Normal, Copy_Bicubic_Patch,
  38.    Translate_Bicubic_Patch, Rotate_Bicubic_Patch,
  39.    Scale_Bicubic_Patch, Invert_Bicubic_Patch};
  40.  
  41. extern long Ray_Bicubic_Tests, Ray_Bicubic_Tests_Succeeded;
  42. extern unsigned int Options;
  43. extern RAY *VP_Ray;
  44. extern int Shadow_Test_Flag;
  45.  
  46. int max_depth_reached;
  47.  
  48. static void bezier_value PARAMS((VECTOR *Result, DBL u, DBL v,
  49. VECTOR (*Control_Points)[4][4]));
  50. static void bezier_partial PARAMS((VECTOR *Result, DBL u, DBL v,
  51. BICUBIC_PATCH *Shape));
  52. static int subpatch_normal PARAMS((VECTOR *v1, VECTOR *v2, VECTOR *v3,
  53. VECTOR *Result, DBL *d));
  54. static int intersect_subpatch PARAMS((int Patch_Type, RAY *Ray,
  55. VECTOR *v1, VECTOR *v2, VECTOR *v3,
  56. VECTOR *n, DBL d,
  57. VECTOR *n1, VECTOR *n2, VECTOR *n3,
  58. DBL *Depth, VECTOR *IP, VECTOR *IP_Norm));
  59. static void find_average PARAMS((int vector_count, VECTOR *vectors,
  60. VECTOR *center, DBL *radius));
  61. static int spherical_bounds_check PARAMS((RAY *Ray, VECTOR *center,
  62. DBL radius));
  63. static int intersect_bicubic_patch0 PARAMS((RAY *Ray, BICUBIC_PATCH *Shape,
  64. DBL *Depths));
  65. static int intersect_bicubic_patch1 PARAMS((RAY *Ray, BICUBIC_PATCH *Shape,
  66. DBL *Depths));
  67. static int intersect_bicubic_patch2 PARAMS((RAY *Ray, BICUBIC_PATCH *Shape,
  68. DBL *Depths));
  69. static int intersect_bicubic_patch3 PARAMS((RAY *Ray, BICUBIC_PATCH *Shape,
  70. DBL *Depths));
  71. static int intersect_bicubic_patch4 PARAMS((RAY *Ray, BICUBIC_PATCH *Shape,
  72. DBL *Depths));
  73. static DBL point_plane_distance PARAMS((VECTOR *, VECTOR *, DBL *));
  74. static DBL determine_subpatch_flatness PARAMS((VECTOR (*)[4][4]));
  75. static int flat_enough PARAMS((BICUBIC_PATCH *, VECTOR (*)[4][4]));
  76. static void bezier_bounding_sphere PARAMS((VECTOR (*)[4][4], VECTOR *,DBL *));
  77. static void bezier_subpatch_intersect PARAMS((RAY *, BICUBIC_PATCH *,
  78. VECTOR (*)[4][4], DBL, DBL, DBL,
  79. int, int *, DBL *, DBL *, DBL *));
  80. static void bezier_split_left_right PARAMS((VECTOR (*)[4][4],VECTOR (*)[4][4],
  81. VECTOR (*)[4][4]));
  82. static void bezier_split_up_down PARAMS((VECTOR (*)[4][4], VECTOR (*)[4][4],
  83. VECTOR (*)[4][4]));
  84. static void bezier_subdivider PARAMS((RAY *, BICUBIC_PATCH *,VECTOR (*)[4][4],
  85. DBL, DBL, DBL, DBL, int, int *, DBL *,
  86. DBL *, DBL *));
  87. static void bezier_tree_deleter PARAMS((BEZIER_NODE *Node));
  88. static BEZIER_NODE *bezier_tree_builder PARAMS((BICUBIC_PATCH *,
  89. VECTOR(*)[4][4], int));
  90. static void bezier_tree_walker PARAMS((RAY *, BICUBIC_PATCH *, BEZIER_NODE *,
  91. int, int *, DBL *));
  92. static BEZIER_NODE *create_new_bezier_node PARAMS((void));
  93. static BEZIER_VERTICES *create_bezier_vertex_block PARAMS((void));
  94. static BEZIER_CHILDREN *create_bezier_child_block PARAMS((void));
  95.  
  96.  
  97. static BEZIER_NODE *
  98. create_new_bezier_node()
  99. {
  100.    BEZIER_NODE *Node = (BEZIER_NODE *)malloc(sizeof(BEZIER_NODE));
  101.    if (Node == NULL) {
  102.       printf("Failed to allocate Bezier node\n");
  103.       exit(0);
  104.    }
  105.    Node->Data_Ptr = NULL;
  106.    return Node;
  107. }
  108.  
  109. static BEZIER_VERTICES *
  110. create_bezier_vertex_block()
  111. {
  112.    BEZIER_VERTICES *Vertices = (BEZIER_VERTICES *)malloc(sizeof(BEZIER_VERTICES));
  113.    if (Vertices == NULL) {
  114.       printf("Failed to allocate Bezier vertices\n");
  115.       exit(0);
  116.    }
  117.    return Vertices;
  118. }
  119.  
  120. static BEZIER_CHILDREN *
  121. create_bezier_child_block()
  122. {
  123.    BEZIER_CHILDREN *Children = (BEZIER_CHILDREN *)malloc(sizeof(BEZIER_CHILDREN));
  124.    if (Children == NULL) {
  125.       printf("Failed to allocate Bezier children\n");
  126.       exit(0);
  127.    }
  128.    return Children;
  129. }
  130.  
  131. static BEZIER_NODE *bezier_tree_builder(Object, Patch, depth)
  132. BICUBIC_PATCH *Object;
  133. VECTOR (*Patch)[4][4];
  134. int depth;
  135. {
  136.    VECTOR Lower_Left[4][4], Lower_Right[4][4];
  137.    VECTOR Upper_Left[4][4], Upper_Right[4][4];
  138.    BEZIER_CHILDREN *Children;
  139.    BEZIER_VERTICES *Vertices;
  140.    BEZIER_NODE *Node = create_new_bezier_node();
  141.  
  142.    if (depth > max_depth_reached) max_depth_reached = depth;
  143.  
  144.    /* Build the bounding sphere for this subpatch */
  145.    bezier_bounding_sphere(Patch, &(Node->Center), &(Node->Radius_Squared));
  146.  
  147.    /* If the patch is close to being flat, then just perform a ray-plane
  148.       intersection test. */
  149.    if (flat_enough(Object, Patch)) {
  150.       /* The patch is now flat enough to simply store the 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.       Node->Data_Ptr = (void *)Vertices;
  158.    }
  159.    else if (depth >= Object->U_Steps)
  160.       if (depth >= Object->V_Steps) {
  161.          /* We are at the max recursion depth. Just store corners. */
  162.          Node->Node_Type = BEZIER_LEAF_NODE;
  163.          Vertices = create_bezier_vertex_block();
  164.          Vertices->Vertices[0] = (*Patch)[0][0];
  165.          Vertices->Vertices[1] = (*Patch)[0][3];
  166.          Vertices->Vertices[2] = (*Patch)[3][3];
  167.          Vertices->Vertices[3] = (*Patch)[3][0];
  168.          Node->Data_Ptr = (void *)Vertices;
  169.       }
  170.       else {
  171.          bezier_split_up_down(Patch, (VECTOR (*)[4][4])Lower_Left, (VECTOR (*)[4][4])Upper_Left);
  172.          Node->Node_Type = BEZIER_INTERIOR_NODE;
  173.          Children = create_bezier_child_block();
  174.          Children->Children[0]=bezier_tree_builder(Object,(VECTOR (*)[4][4])Lower_Left,depth+1);
  175.          Children->Children[1]=bezier_tree_builder(Object,(VECTOR (*)[4][4])Upper_Left,depth+1);
  176.          Node->Count = 2;
  177.          Node->Data_Ptr = (void *)Children;
  178.       }
  179.    else if (depth >= Object->V_Steps) {
  180.       bezier_split_left_right(Patch, (VECTOR (*)[4][4])Lower_Left, (VECTOR (*)[4][4])Lower_Right);
  181.       Node->Node_Type = BEZIER_INTERIOR_NODE;
  182.       Children = create_bezier_child_block();
  183.       Children->Children[0] = bezier_tree_builder(Object,(VECTOR (*)[4][4])Lower_Left,depth+1);
  184.       Children->Children[1] = bezier_tree_builder(Object,(VECTOR (*)[4][4])Lower_Right,depth+1);
  185.       Node->Count = 2;
  186.       Node->Data_Ptr = (void *)Children;
  187.    }
  188.    else {
  189.       bezier_split_left_right(Patch, (VECTOR (*)[4][4])Lower_Left, (VECTOR (*)[4][4])Lower_Right);
  190.       bezier_split_up_down((VECTOR (*)[4][4])Lower_Left, (VECTOR (*)[4][4])Lower_Left, (VECTOR (*)[4][4])Upper_Left);
  191.       bezier_split_up_down((VECTOR (*)[4][4])Lower_Right, (VECTOR (*)[4][4])Lower_Right, (VECTOR (*)[4][4])Upper_Right);
  192.       Node->Node_Type = BEZIER_INTERIOR_NODE;
  193.       Children = create_bezier_child_block();
  194.       Children->Children[0] = bezier_tree_builder(Object,(VECTOR (*)[4][4])Lower_Left,depth+1);
  195.       Children->Children[1] = bezier_tree_builder(Object,(VECTOR (*)[4][4])Upper_Left,depth+1);
  196.       Children->Children[2] = bezier_tree_builder(Object,(VECTOR (*)[4][4])Lower_Right,depth+1);
  197.       Children->Children[3] = bezier_tree_builder(Object,(VECTOR (*)[4][4])Upper_Right,depth+1);
  198.       Node->Count = 4;
  199.       Node->Data_Ptr = (void *)Children;
  200.    }
  201.    return Node;
  202. }
  203.  
  204.    /* Evaluate a single coordinate point (u, v) on a bezier patch. */
  205.    static void bezier_value(Result, u, v, Control_Points)
  206.       VECTOR *Result;
  207. DBL u, v;
  208. VECTOR (*Control_Points)[4][4];
  209. {
  210.    DBL u2, u3, v2, v3, uu1, uu2, uu3, vv1, vv2, vv3;
  211.    DBL t[4][4];
  212.    int i, j;
  213.  
  214.    u2 = u * u; u3 = u * u2; v2 = v * v; v3 = v * v2;
  215.    uu1 = 1.0 - u; uu2 = uu1 * uu1; uu3 = uu1 * uu2;
  216.    vv1 = 1.0 - v; vv2 = vv1 * vv1; vv3 = vv1 * vv2;
  217.    t[0][0] =         uu3 *        vv3;
  218.    t[0][1] = 3.0 *     uu3 * v  * vv2;
  219.    t[0][2] = 3.0 *     uu3 * v2 * vv1;
  220.    t[0][3] =         uu3 * v3;
  221.    t[1][0] = 3.0 *  u  * uu2 *        vv3;
  222.    t[1][1] = 9.0 *  u  * uu2 * v  * vv2;
  223.    t[1][2] = 9.0 *  u  * uu2 * v2 * vv1;
  224.    t[1][3] = 3.0 *  u  * uu2 * v3;
  225.    t[2][0] = 3.0 *  u2 * uu1 *        vv3;
  226.    t[2][1] = 9.0 *  u2 * uu1 * v  * vv2;
  227.    t[2][2] = 9.0 *  u2 * uu1 * v2 * vv1;
  228.    t[2][3] = 3.0 *  u2 * uu1 * v3;
  229.    t[3][0] =        u3 *        vv3;
  230.    t[3][1] = 3.0 *  u3 *       v  * vv2;
  231.    t[3][2] = 3.0 *  u3 *       v2 * vv1;
  232.    t[3][3] =        u3 *       v3;
  233.  
  234.    Result->x = 0.0;
  235.    Result->y = 0.0;
  236.    Result->z = 0.0;
  237.    for (i=0;i<4;i++) for (j=0;j<4;j++) {
  238.       Result->x += t[i][j] * (*Control_Points)[i][j].x;
  239.       Result->y += t[i][j] * (*Control_Points)[i][j].y;
  240.       Result->z += t[i][j] * (*Control_Points)[i][j].z;
  241.    }
  242. }
  243.  
  244. /* Calculate the normal to a bezier patch for a particular axis, at
  245.    a particular point on the patch.
  246.  
  247.    The normal at a point of a parametric surface z = f(u, v) is:
  248.  
  249.       (|[[dy/du, dy/dv],[dz/du, dz/dv]]|,
  250.        |[[dz/du, dz/dv],[dx/du, dx/dv]]|,
  251.        |[[dx/du, dx/dv],[dy/du, dy/dv]]|)
  252.  
  253.    The normal is undefined where the determinants vanish.
  254. */
  255. static void bezier_partial(Result, u, v, Shape)
  256. VECTOR *Result;
  257. DBL u, v;
  258. BICUBIC_PATCH *Shape;
  259. {
  260.    VECTOR U, V; /* Partial derivatives with respect to u, and v. */
  261.    DBL u2, u3, v2, v3;
  262.    DBL t[4][4], Temp;
  263.    int i, j;
  264.  
  265.    u2 = u * u; u3 = u * u2; v2 = v * v; v3 = v * v2;
  266.  
  267.    /* Calculate the derivative with respect to u */
  268.    t[0][0] = 3.0 *    (v3-3.0*v2+3.0*v-1.0) * (u2-2.0*u+1.0);
  269.    t[0][1] = 9.0 * v  * (v2-2.0*v+1.0) * (u2-2.0*u+1.0);
  270.    t[0][2] = 9.0 * v2 * (v-1.0) * (u2-2.0*u+1.0);
  271.    t[0][3] = 3.0 * v3 * (u2-2.0*u+1.0);
  272.    t[1][0] = 3.0 *    (v3-3.0*v2+3.0*v-1) * (3.0*u2-4.0*u+1);
  273.    t[1][1] = 9.0 * v  * (v2-2.0*v+1.0) * (3.0*u2-4.0*u+1.0);
  274.    t[1][2] = 9.0 * v2 * (v-1.0) * (3.0*u2-4.0*u+1.0);
  275.    t[1][3] = 3.0 * v3 * (3.0*u2-4.0*u+1.0);
  276.    t[2][0] = 3.0 * u  * (v3-3.0*v2+3.0*v-1.0)*(3.0*u-2.0);
  277.    t[2][1] = 9.0 * u  * v  * (v2-2.0*v+1.0)*(3.0*u-2.0);
  278.    t[2][2] = 9.0 * u  * v2 * (v-1.0)*(3.0*u-2.0);
  279.    t[2][3] = 3.0 * u  * v3 * (3.0*u-2.0);
  280.    t[3][0] = 3.0 * u2 * (v3-3.0*v2+3.0*v-1.0);
  281.    t[3][1] = 9.0 * u2 * v  * (v2-2.0*v+1.0);
  282.    t[3][2] = 9.0 * u2 * v2 * (v-1.0);
  283.    t[3][3] = 3.0 * u2 * v3;
  284.    U.x = 0.0; U.y = 0.0; U.z = 0.0;
  285.    for (i=0;i<4;i++) for (j=0;j<4;j++) {
  286.       U.x += t[i][j] * Shape->Control_Points[i][j].x;
  287.       U.y += t[i][j] * Shape->Control_Points[i][j].y;
  288.       U.z += t[i][j] * Shape->Control_Points[i][j].z;
  289.    }
  290.    Temp = U.x * U.x + U.y * U.y + U.z * U.z;
  291.    if (Temp < EPSILON) {
  292.       /* Partial with respect to u is undefined. */
  293.       Result->x = 1.0; Result->y = 0.0; Result->z = 0.0;
  294.       /* *Result = *n; */
  295.       return;
  296.    }
  297.    else {
  298.       Temp = sqrt(Temp);
  299.       VInverseScale(U, U, Temp);
  300.    }
  301.    /* Calculate the derivative with respect to v */
  302.       t[0][0] = 3.0 * (v2-2.0*v+1.0) * (u3-3.0*u2+3.0*u-1.0);
  303.    t[0][1] = 3.0 * (3.0*v2-4.0*v+1.0) * (u3-3.0*u2+3.0*u-1.0);
  304.    t[0][2] = 3.0 * v  * (3.0*v-2.0) * (u3-3.0*u2+3.0*u-1.0);
  305.    t[0][3] = 3.0 * v2 * (u3-3.0*u2+3.0*u-1.0);
  306.    t[1][0] = 9.0 * u  * (v2-2.0*v+1.0) * (u2-2.0*u+1.0);
  307.    t[1][1] = 9.0 * u  * (3.0*v2-4.0*v+1.0) * (u2-2.0*u+1.0);
  308.    t[1][2] = 9.0 * u  * v * (3.0*v-2.0) * (u2-2.0*u+1.0);
  309.    t[1][3] = 9.0 * u  * v2 * (u2-2.0*u+1.0);
  310.    t[2][0] = 9.0 * u2 * (v2-2.0*v+1.0) * (u-1.0);
  311.    t[2][1] = 9.0 * u2 * (3.0*v2-4.0*v+1.0) * (u-1.0);
  312.    t[2][2] = 9.0 * u2 * v * (3.0*v-2.0) * (u-1.0);
  313.    t[2][3] = 9.0 * u2 * v2 * (u-1.0);
  314.    t[3][0] = 3.0 * u3 * (v2-2.0*v+1.0);
  315.    t[3][1] = 3.0 * u3 * (3.0*v2-4.0*v+1.0);
  316.    t[3][2] = 3.0 * u3 * v  * (3.0*v-2.0);
  317.    t[3][3] = 3.0 * u3 * v2;
  318.    V.x = 0.0; V.y = 0.0; V.z = 0.0;
  319.    for (i=0;i<4;i++) for (j=0;j<4;j++) {
  320.       V.x += t[i][j] * Shape->Control_Points[i][j].x;
  321.       V.y += t[i][j] * Shape->Control_Points[i][j].y;
  322.       V.z += t[i][j] * Shape->Control_Points[i][j].z;
  323.    }
  324.    Temp = V.x * V.x + V.y * V.y + V.z * V.z;
  325.    if (Temp < EPSILON) {
  326.       /* Partial with respect to u is undefined. */
  327.       Result->x = 1.0; Result->y = 0.0; Result->z = 0.0;
  328.       return;
  329.    }
  330.    else {
  331.       Temp = sqrt(Temp);
  332.       VInverseScale(V, V, Temp);
  333.    }
  334.    VCross(*Result, U, V);
  335. }
  336.  
  337.    /* Calculate the normal to a subpatch (triangle) return the vector
  338.    <1.0 0.0 0.0> if the triangle is degenerate. */
  339.    static int subpatch_normal(v1, v2, v3, Result, d)
  340.       VECTOR *v1, *v2, *v3, *Result;
  341. DBL *d;
  342. {
  343.    VECTOR V1, V2;
  344.    DBL Length;
  345.  
  346.    VSub(V1, *v1, *v2);
  347.    VSub(V2, *v3, *v2);
  348.    VCross(*Result, V1, V2);
  349.    VLength(Length, *Result);
  350.    if (Length < EPSILON) {
  351.       Result->x = 1.0;
  352.       Result->y = 0.0;
  353.       Result->z = 0.0;
  354.       *d = -1.0 * v1->x;
  355.       return 0;
  356.    }
  357.    else {
  358.       VInverseScale(*Result, *Result, Length);
  359.       VDot(*d, *Result, *v1);
  360.       *d = 0.0 - *d;
  361.       return 1;
  362.    }
  363. }
  364.  
  365.    /* See if Ray intersects the triangular subpatch defined by the three
  366.    points v1, v2, v3 and the normal to the triangle n. If so, Depth will
  367.    be set to the distance along Ray at which the intersection occurs. */
  368.    static int intersect_subpatch(Patch_Type, Ray, v1, v2, v3, n, d, n1, n2, n3,
  369.       Depth, IP, IP_Norm)
  370.       int Patch_Type;
  371. RAY *Ray;
  372. VECTOR *v1, *v2, *v3, *n, *n1, *n2, *n3, *IP, *IP_Norm;
  373. DBL d, *Depth;
  374. {
  375.    VECTOR TempV1, TempV2, Perp;
  376.    DBL s, t, Proj, mu;
  377.    DBL x1, y1, x2, y2, x3, y3;
  378.    DBL x, y, z;
  379.    int Sign_Holder, Next_Sign, Crossings;
  380.  
  381.    /* Calculate the point of intersection and the depth. */
  382.    VDot(s, Ray->Direction, *n);
  383.    if (s == 0.0)
  384.       return 0;
  385.    VDot(t, Ray->Initial, *n);
  386.    *Depth = 0.0 - (d + t) / s;
  387.    if (*Depth < Small_Tolerance)
  388.       return 0;
  389.    VScale(*IP, Ray->Direction, *Depth);
  390.    VAdd(*IP, *IP, Ray->Initial);
  391.  
  392.    /* Map the intersection point and the triangle onto a plane. */
  393.    x = fabs(n->x); y = fabs(n->y); z = fabs(n->z);
  394.    if (x>y)
  395.       if (x>z) {
  396.          x1 = v1->y - IP->y; y1 = v1->z - IP->z;
  397.          x2 = v2->y - IP->y; y2 = v2->z - IP->z;
  398.          x3 = v3->y - IP->y; y3 = v3->z - IP->z;
  399.       }
  400.       else {
  401.          x1 = v1->x - IP->x; y1 = v1->y - IP->y;
  402.          x2 = v2->x - IP->x; y2 = v2->y - IP->y;
  403.          x3 = v3->x - IP->x; y3 = v3->y - IP->y;
  404.       }
  405.    else if (y>z) {
  406.       x1 = v1->x - IP->x; y1 = v1->z - IP->z;
  407.       x2 = v2->x - IP->x; y2 = v2->z - IP->z;
  408.       x3 = v3->x - IP->x; y3 = v3->z - IP->z;
  409.    }
  410.    else {
  411.       x1 = v1->x - IP->x; y1 = v1->y - IP->y;
  412.       x2 = v2->x - IP->x; y2 = v2->y - IP->y;
  413.       x3 = v3->x - IP->x; y3 = v3->y - IP->y;
  414.    }
  415.  
  416.       /* Determine crossing count */
  417.       Crossings = 0;
  418.    if (y1 < 0) Sign_Holder = -1;
  419.    else Sign_Holder = 1;
  420.  
  421.       /* Start of winding tests, test the segment from v1 to v2 */
  422.       if (y2 < 0) Next_Sign = -1;
  423.       else Next_Sign = 1;
  424.    if (Sign_Holder != Next_Sign) {
  425.       if (x1 > 0.0) {
  426.          if (x2 > 0.0)
  427.             Crossings++;
  428.          else if ((x1 - y1 * (x2 - x1) / (y2 - y1)) > 0.0)
  429.             Crossings++;
  430.       }
  431.       else if (x2 > 0.0) {
  432.          if ((x1 - y1 * (x2 - x1) / (y2 - y1)) > 0.0)
  433.             Crossings++;
  434.       }
  435.       Sign_Holder = Next_Sign;
  436.    }
  437.  
  438.    /* Test the segment from v2 to v3 */
  439.    if (y3 < 0) Next_Sign = -1;
  440.    else Next_Sign = 1;
  441.    if (Sign_Holder != Next_Sign) {
  442.       if (x2 > 0.0) {
  443.          if (x3 > 0.0)
  444.             Crossings++;
  445.          else if ((x2 - y2 * (x3 - x2) / (y3 - y2)) > 0.0)
  446.             Crossings++;
  447.       }
  448.       else if (x3 > 0.0) {
  449.          if ((x2 - y2 * (x3 - x2) / (y3 - y2)) > 0.0)
  450.             Crossings++;
  451.       }
  452.       Sign_Holder = Next_Sign;
  453.    }
  454.  
  455.    /* Test the segment from v3 to v1 */
  456.    if (y1 < 0) Next_Sign = -1;
  457.    else Next_Sign = 1;
  458.    if (Sign_Holder != Next_Sign) {
  459.       if (x3 > 0.0) {
  460.          if (x1 > 0.0)
  461.             Crossings++;
  462.          else if ((x3 - y3 * (x1 - x3) / (y1 - y3)) > 0.0)
  463.             Crossings++;
  464.       }
  465.       else if (x1 > 0.0) {
  466.          if ((x3 - y3 * (x1 - x3) / (y1 - y3)) > 0.0)
  467.             Crossings++;
  468.       }
  469.       Sign_Holder = Next_Sign;
  470.    }
  471.    if (Crossings != 1)
  472.       return 0;
  473.  
  474.    if (Patch_Type == 0 || Patch_Type == 1 || Patch_Type == 2 ||
  475.       Patch_Type == 3) {
  476.       *IP_Norm = *n;
  477.       return 1;
  478.    }
  479.    else if (Patch_Type == 4) {
  480.       VSub(TempV1, *v2, *v3);
  481.       VNormalize(TempV1, TempV1);
  482.       VSub(TempV2, *v1, *v3);
  483.       VDot(Proj, TempV2, TempV1);
  484.       VScale(TempV1, TempV1, Proj);
  485.       VSub(Perp, TempV1, TempV2);
  486.       VNormalize(Perp, Perp);
  487.       VDot(mu, TempV2, Perp);
  488.       mu = -1.0 / mu;
  489.       VScale(Perp, Perp, mu);
  490.       VSub(TempV1, *IP, *v1);
  491.       VDot(s, TempV1, Perp);
  492.       if (s < EPSILON) {
  493.          *IP_Norm = *n1;
  494.          return 1;
  495.       }
  496.       VSub(TempV1, *v3, *v2);
  497.       x = fabs(TempV1.x); y = fabs(TempV1.y); z = fabs(TempV1.z);
  498.       if (x>y)
  499.          if (x>z)
  500.             t = (TempV1.x / s + v1->x - v2->x) / TempV1.x;
  501.          else
  502.             t = (TempV1.z / s + v1->z - v2->z) / TempV1.z;
  503.       else if (y > z)
  504.          t = (TempV1.y / s + v1->y - v2->y) / TempV1.y;
  505.       else
  506.          t = (TempV1.z / s + v1->z - v2->z) / TempV1.z;
  507.       VSub(TempV1, *n2, *n1);
  508.       VScale(TempV1, TempV1, s);
  509.       VAdd(TempV1, TempV1, *n1);
  510.       VSub(TempV2, *n3, *n1);
  511.       VScale(TempV2, TempV2, s);
  512.       VAdd(TempV2, TempV2, *n1);
  513.       VSub(*IP_Norm, TempV2, TempV1);
  514.       VScale(*IP_Norm, *IP_Norm, t);
  515.       VAdd(*IP_Norm, *IP_Norm, TempV1);
  516.       VNormalize(*IP_Norm, *IP_Norm);
  517.       return 1;
  518.    }
  519.    else
  520.       return 0;
  521. }
  522.  
  523. /* Find a sphere that contains all of the points in the list "vectors" */
  524. static void find_average(vector_count, vectors, center, radius)
  525. int vector_count;
  526. VECTOR *vectors, *center;
  527. DBL *radius;
  528. {
  529.    DBL r0, r1, xc=0, yc=0, zc=0;
  530.    DBL x0, y0, z0;
  531.    int i;
  532.    for (i=0;i<vector_count;i++) {
  533.       xc += vectors[i].x;
  534.       yc += vectors[i].y;
  535.       zc += vectors[i].z;
  536.    }
  537.    xc /= (DBL)vector_count;
  538.    yc /= (DBL)vector_count;
  539.    zc /= (DBL)vector_count;
  540.    r0 = 0.0;
  541.    for (i=0;i<vector_count;i++) {
  542.       x0 = vectors[i].x - xc;
  543.       y0 = vectors[i].y - yc;
  544.       z0 = vectors[i].z - zc;
  545.       r1 = x0 * x0 + y0 * y0 + z0 * z0;
  546.       if (r1 > r0) r0 = r1;
  547.    }
  548.    center->x = xc; center->y = yc; center->z = zc;
  549.    *radius = r0;
  550. }
  551.  
  552. #define SUBDIVISION_EPSILON 0.001
  553. #define MAX_RECURSION_DEPTH 20
  554.  
  555. static int spherical_bounds_check(Ray, center, radius)
  556. RAY *Ray;
  557. VECTOR *center;
  558. DBL radius;
  559. {
  560.    DBL x, y, z, dist1, dist2;
  561.    x = center->x - Ray->Initial.x;
  562.    y = center->y - Ray->Initial.y;
  563.    z = center->z - Ray->Initial.z;
  564.    dist1 = x * x + y * y + z * z;
  565.    if (dist1 < radius)
  566.       /* ray starts inside sphere - assume it intersects. */
  567.       return 1;
  568.    else {
  569.       dist2 = x*Ray->Direction.x + y*Ray->Direction.y + z*Ray->Direction.z;
  570.       dist2 = dist2 * dist2;
  571.       if (dist2 > 0 && (dist1 - dist2 < radius))
  572.          return 1;
  573.    }
  574.    return 0;
  575. }
  576.  
  577. /* Find a sphere that bounds all of the control points of a Bezier patch.
  578.    The values returned are: the center of the bounding sphere, and the
  579.    square of the radius of the bounding sphere. */
  580. static void bezier_bounding_sphere(Patch, center, radius)
  581. VECTOR (*Patch)[4][4], *center;
  582. DBL *radius;
  583. {
  584.    DBL r0, r1, xc=0, yc=0, zc=0;
  585.    DBL x0, y0, z0;
  586.    int i, j;
  587.    for (i=0;i<4;i++) {
  588.       for (j=0;j<4;j++) {
  589.          xc += (*Patch)[i][j].x;
  590.          yc += (*Patch)[i][j].y;
  591.          zc += (*Patch)[i][j].z;
  592.       }
  593.    }
  594.    xc /= 16.0;
  595.    yc /= 16.0;
  596.    zc /= 16.0;
  597.    r0 = 0.0;
  598.    for (i=0;i<4;i++) {
  599.       for (j=0;j<4;j++) {
  600.          x0 = (*Patch)[i][j].x - xc;
  601.          y0 = (*Patch)[i][j].y - yc;
  602.          z0 = (*Patch)[i][j].z - zc;
  603.          r1 = x0 * x0 + y0 * y0 + z0 * z0;
  604.          if (r1 > r0) r0 = r1;
  605.       }
  606.    }
  607.    center->x = xc; center->y = yc; center->z = zc;
  608.    *radius = r0;
  609. }
  610.  
  611. /* Precompute grid points and normals for a bezier patch */
  612. void
  613. Precompute_Patch_Values(Shape)
  614. BICUBIC_PATCH *Shape;
  615. {
  616.    int i, j;
  617.    DBL d, u, v, delta_u, delta_v;
  618.    VECTOR v0, v1, v2, v3, n;
  619.    VECTOR Control_Points[16];
  620.    VECTOR (*Patch_Ptr)[4][4] = (VECTOR (*)[4][4]) Shape->Control_Points;
  621.  
  622.    /* Calculate the bounding sphere for the entire patch. */
  623.    for (i=0;i<4;i++) for (j=0;j<4;j++)
  624.       Control_Points[4*i+j] = Shape->Control_Points[i][j];
  625.    find_average(16, &Control_Points[0], &Shape->Bounding_Sphere_Center,
  626.       &Shape->Bounding_Sphere_Radius);
  627.    /* Shape->Node_Tree = NULL; */
  628.    if (Shape->Patch_Type == 0 || Shape->Patch_Type == 2)
  629.       return;
  630.    else if (Shape->Patch_Type == 3) {
  631.       if (Shape->Node_Tree != NULL)
  632.          bezier_tree_deleter(Shape->Node_Tree);
  633.       Shape->Node_Tree = bezier_tree_builder(Shape, Patch_Ptr, 0);
  634.       return;
  635.    }
  636.    delta_u = 1.0 / (DBL)Shape->U_Steps;
  637.    delta_v = 1.0 / (DBL)Shape->V_Steps;
  638.    if (Shape->Interpolated_Grid == NULL) {
  639.       Shape->Interpolated_Grid = (VECTOR **)
  640.          malloc((Shape->U_Steps+1) * sizeof(VECTOR *));
  641.       if (Shape->Interpolated_Grid == NULL)
  642.          Error("Failed to allocate Interpolated_Grid");
  643.       for (i=0;i<=Shape->U_Steps;i++) {
  644.          Shape->Interpolated_Grid[i] = (VECTOR *)
  645.             malloc((Shape->V_Steps+1) * sizeof(VECTOR));
  646.          if (Shape->Interpolated_Grid == NULL)
  647.             Error("Failed to allocate component of Interpolated_Grid");
  648.       }
  649.       Shape->Interpolated_Normals = (VECTOR **)
  650.          malloc((Shape->U_Steps+1) * sizeof(VECTOR *));
  651.       if (Shape->Interpolated_Normals == NULL)
  652.          Error("Failed to allocate Interpolated_Normals");
  653.       for (i=0;i<=Shape->U_Steps;i++) {
  654.          Shape->Interpolated_Normals[i] = (VECTOR *)
  655.             malloc(2*(Shape->V_Steps+1)*sizeof(VECTOR));
  656.          if (Shape->Interpolated_Normals == NULL)
  657.             Error("Failed to allocate component of Interpolated_Normals");
  658.       }
  659.  
  660.       if (Shape->Patch_Type == 4) {
  661.          Shape->Smooth_Normals = (VECTOR **)
  662.             malloc((Shape->U_Steps+1) * sizeof(VECTOR *));
  663.          if (Shape->Smooth_Normals == NULL)
  664.             Error("Failed to allocate Smooth_Normals");
  665.          for (i=0;i<=Shape->U_Steps;i++) {
  666.             Shape->Smooth_Normals[i] = (VECTOR *)
  667.                malloc((Shape->V_Steps+1)*sizeof(VECTOR));
  668.             if (Shape->Smooth_Normals == NULL)
  669.                Error("Failed to allocate component of Smooth_Normals");
  670.          }
  671.       }
  672.  
  673.       Shape->Interpolated_D = (DBL **)malloc((Shape->U_Steps+1) * sizeof(DBL *));
  674.       if (Shape->Interpolated_D == NULL)
  675.          Error("Failed to allocate Interpolated_D");
  676.       for (i=0;i<=Shape->U_Steps;i++) {
  677.          Shape->Interpolated_D[i] = (DBL *)
  678.             malloc(2*(Shape->V_Steps+1)*sizeof(DBL));
  679.          if (Shape->Interpolated_D == NULL)
  680.             Error("Failed to allocate component of Interpolated_D");
  681.       }
  682.    }
  683.  
  684.    /* Calculate the grid values for the given subdivision values. */
  685.    for (i=0;i<=Shape->U_Steps;i++) {
  686.       u = (DBL)i / (DBL)Shape->U_Steps;
  687.       for (j=0;j<Shape->V_Steps;j++) {
  688.          v = (DBL) j / (DBL)Shape->V_Steps;
  689.          bezier_value(&Shape->Interpolated_Grid[i][j], u, v, Patch_Ptr);
  690.       }
  691.    }
  692.  
  693.    for (i=0;i<Shape->U_Steps;i++) {
  694.       u = (DBL)i / (DBL)Shape->U_Steps;
  695.       for (j=0;j<Shape->V_Steps;j++) {
  696.          v = (DBL)j / (DBL)Shape->V_Steps;
  697.  
  698.          /* Calculate surface values for the current patch. */
  699.          bezier_value(&v0, u, v, Patch_Ptr);
  700.          bezier_value(&v1, u+delta_u, v, Patch_Ptr);
  701.          bezier_value(&v2, u, v+delta_v, Patch_Ptr);
  702.          bezier_value(&v3, u+delta_u, v+delta_v, Patch_Ptr);
  703.  
  704.          Shape->Interpolated_Grid[i][j] = v0;
  705.          Shape->Interpolated_Grid[i+1][j] = v1;
  706.          Shape->Interpolated_Grid[i][j+1] = v2;
  707.          Shape->Interpolated_Grid[i+1][j+1] = v3;
  708.          if (Shape->Patch_Type == 1 || Shape->Patch_Type == 4) {
  709.             /* Calculate the normals */
  710.             if (subpatch_normal(&v0, &v2, &v1, &n, &d)) {
  711.                Shape->Interpolated_Normals[i][2*j] = n;
  712.                Shape->Interpolated_D[i][2*j] = d;
  713.             }
  714.             else {
  715.                Shape->Interpolated_Normals[i][2*j].x = 0.0;
  716.                Shape->Interpolated_Normals[i][2*j].y = 0.0;
  717.                Shape->Interpolated_Normals[i][2*j].z = 0.0;
  718.                Shape->Interpolated_D[i][2*j] = 0.0;
  719.             }
  720.  
  721.                if (subpatch_normal(&v1, &v2, &v3, &n, &d)) {
  722.                   Shape->Interpolated_Normals[i][2*j+1] = n;
  723.                   Shape->Interpolated_D[i][2*j+1] = d;
  724.                }
  725.                else {
  726.                   Shape->Interpolated_Normals[i][2*j+1].x = 0.0;
  727.                   Shape->Interpolated_Normals[i][2*j+1].y = 0.0;
  728.                   Shape->Interpolated_Normals[i][2*j+1].z = 0.0;
  729.                   Shape->Interpolated_D[i][2*j+1] = 0.0;
  730.                }
  731.          }
  732.       }
  733.    }
  734.  
  735.       if (Shape->Patch_Type == 4) {
  736.          /* Calculate normals at the corners of the subpatches */
  737.          for (i=0;i<=Shape->U_Steps;i++) {
  738.             u = (DBL)i / (DBL)Shape->U_Steps;
  739.             for (j=0;j<=Shape->V_Steps;j++) {
  740.                v = (DBL)j / (DBL)Shape->V_Steps;
  741.                bezier_partial(&Shape->Smooth_Normals[i][j], u, v, Shape);
  742.             }
  743.          }
  744.       }
  745. }
  746.  
  747. /* Determine the distance from a point to a plane. */
  748. static DBL point_plane_distance(p, n, d)
  749. VECTOR *p, *n;
  750. DBL *d;
  751. {
  752.    DBL temp1, temp2;
  753.  
  754.    VDot(temp1, *p, *n);
  755.    temp1 += *d;
  756.    VLength(temp2, *n);
  757.    if (fabs(temp2) < EPSILON) return 0;
  758.    temp1 /= temp2;
  759.    return temp1;
  760. }
  761.  
  762. static void bezier_subpatch_intersect(Ray, Shape, Patch, u0, u1, v0,
  763. recursion_depth, depth_count, Depths, U_Values, V_Values)
  764. RAY *Ray;
  765. BICUBIC_PATCH *Shape;
  766. VECTOR (*Patch)[4][4];
  767. DBL u0, u1, v0;
  768. int recursion_depth, *depth_count;
  769. DBL *Depths, *U_Values, *V_Values;
  770. {
  771.    int tcnt = Shape->Intersection_Count;
  772.    VECTOR vv0, vv1, vv2, vv3, n, ip;
  773.    DBL d, Depth;
  774.  
  775.    if (tcnt + *depth_count >= MAX_BICUBIC_INTERSECTIONS) return;
  776.  
  777.    vv0 = (*Patch)[0][0]; vv1 = (*Patch)[0][3];
  778.    vv2 = (*Patch)[3][3]; vv3 = (*Patch)[3][0];
  779.  
  780.    /* Triangulate this subpatch, then check for intersections in
  781.       the triangles. */
  782.    if (subpatch_normal(&vv0, &vv1, &vv2, &n, &d)) {
  783.       if (intersect_subpatch(Shape->Patch_Type, Ray, &vv0, &vv1, &vv2, &n, d,
  784.          NULL, NULL, NULL, &Depth, &ip, &n)) {
  785.          Shape->Intersection_Point[tcnt + *depth_count] = ip;
  786.          Shape->Normal_Vector[tcnt + *depth_count] = n;
  787.          Depths[*depth_count] = Depth;
  788.          *depth_count += 1;
  789.       }
  790.    }
  791.  
  792.    if (*depth_count + tcnt >= MAX_BICUBIC_INTERSECTIONS) return;
  793.  
  794.    if (subpatch_normal(&vv0, &vv2, &vv3, &n, &d)) {
  795.       if (intersect_subpatch(Shape->Patch_Type, Ray, &vv0, &vv2, &vv3, &n, d,
  796.          NULL, NULL, NULL, &Depth, &ip, &n)) {
  797.          Shape->Intersection_Point[tcnt + *depth_count] = ip;
  798.          Shape->Normal_Vector[tcnt + *depth_count] = n;
  799.          Depths[*depth_count] = Depth;
  800.          *depth_count += 1;
  801.       }
  802.    }
  803. }
  804.  
  805. static void bezier_split_left_right(Patch, Left_Patch, Right_Patch)
  806. VECTOR (*Patch)[4][4], (*Left_Patch)[4][4], (*Right_Patch)[4][4];
  807. {
  808.    VECTOR Temp1[4], Temp2[4], Half;
  809.    int i, j;
  810.    for (i=0;i<4;i++) {
  811.       Temp1[0] = (*Patch)[i][0];
  812.       VHalf(Temp1[1], (*Patch)[i][0], (*Patch)[i][1]);
  813.       VHalf(Half, (*Patch)[i][1], (*Patch)[i][2]);
  814.       VHalf(Temp1[2], Temp1[1], Half);
  815.       VHalf(Temp2[2], (*Patch)[i][2], (*Patch)[i][3]);
  816.       VHalf(Temp2[1], Half, Temp2[2]);
  817.       VHalf(Temp1[3], Temp1[2], Temp2[1]);
  818.       Temp2[0] = Temp1[3];
  819.       Temp2[3] = (*Patch)[i][3];
  820.       for (j=0;j<4;j++) {
  821.          (*Left_Patch)[i][j]  = Temp1[j];
  822.          (*Right_Patch)[i][j] = Temp2[j];
  823.       }
  824.    }
  825. }
  826.  
  827. static void bezier_split_up_down(Patch, Top_Patch, Bottom_Patch)
  828. VECTOR (*Patch)[4][4], (*Top_Patch)[4][4], (*Bottom_Patch)[4][4];
  829. {
  830.    VECTOR Temp1[4], Temp2[4], Half;
  831.    int i, j;
  832.  
  833.    for (i=0;i<4;i++) {
  834.       /* Split Left */
  835.       Temp1[0] = (*Patch)[0][i];
  836.       VHalf(Temp1[1], (*Patch)[0][i], (*Patch)[1][i]);
  837.       VHalf(Half, (*Patch)[1][i], (*Patch)[2][i]);
  838.       VHalf(Temp1[2], Temp1[1], Half);
  839.       VHalf(Temp2[2], (*Patch)[2][i], (*Patch)[3][i]);
  840.       VHalf(Temp2[1], Half, Temp2[2]);
  841.       VHalf(Temp1[3], Temp1[2], Temp2[1]);
  842.       Temp2[0] = Temp1[3];
  843.       Temp2[3] = (*Patch)[3][i];
  844.       for (j=0;j<4;j++) {
  845.          (*Bottom_Patch)[j][i] = Temp1[j];
  846.          (*Top_Patch)[j][i]       = Temp2[j];
  847.       }
  848.    }
  849. }
  850.  
  851. /* See how close to a plane a subpatch is, the patch must have at least
  852.    three distinct vertices. A negative result from this function indicates
  853.    that a degenerate value of some sort was encountered. */
  854. static DBL determine_subpatch_flatness(Patch)
  855. VECTOR (*Patch)[4][4];
  856. {
  857.    VECTOR vertices[4], n, TempV;
  858.    DBL d, dist, temp1;
  859.    int i, j;
  860.  
  861.    vertices[0] = (*Patch)[0][0];
  862.    vertices[1] = (*Patch)[0][3];
  863.    VSub(TempV, vertices[0], vertices[1]);
  864.    VLength(temp1, TempV);
  865.    if (fabs(temp1) < EPSILON) {
  866.       /* Degenerate in the V direction for U = 0. This is ok if the other
  867.      two corners are distinct from the lower left corner - I'm sure there
  868.      are cases where the corners coincide and the middle has good values,
  869.      but that is somewhat pathalogical and won't be considered. */
  870.       vertices[1] = (*Patch)[3][3];
  871.       VSub(TempV, vertices[0], vertices[1]);
  872.       VLength(temp1, TempV);
  873.       if (fabs(temp1) < EPSILON) return -1.0;
  874.       vertices[2] = (*Patch)[3][0];
  875.       VSub(TempV, vertices[0], vertices[1]);
  876.       VLength(temp1, TempV);
  877.       if (fabs(temp1) < EPSILON) return -1.0;
  878.       VSub(TempV, vertices[1], vertices[2]);
  879.       VLength(temp1, TempV);
  880.       if (fabs(temp1) < EPSILON) return -1.0;
  881.    }
  882.    else {
  883.       vertices[2] = (*Patch)[3][0];
  884.       VSub(TempV, vertices[0], vertices[1]);
  885.       VLength(temp1, TempV);
  886.       if (fabs(temp1) < EPSILON) {
  887.          vertices[2] = (*Patch)[3][3];
  888.          VSub(TempV, vertices[0], vertices[2]);
  889.          VLength(temp1, TempV);
  890.          if (fabs(temp1) < EPSILON)
  891.             return -1.0;
  892.          VSub(TempV, vertices[1], vertices[2]);
  893.          VLength(temp1, TempV);
  894.          if (fabs(temp1) < EPSILON)
  895.             return -1.0;
  896.       }
  897.       else {
  898.          VSub(TempV, vertices[1], vertices[2]);
  899.          VLength(temp1, TempV);
  900.          if (fabs(temp1) < EPSILON)
  901.             return -1.0;
  902.       }
  903.    }
  904.    /* Now that a good set of candidate points has been found, find the
  905.       plane equations for the patch */
  906.    if (subpatch_normal(&vertices[0], &vertices[1], &vertices[2], &n, &d)) {
  907.       /* Step through all vertices and see what the maximum distance from the
  908.          plane happens to be. */
  909.       dist = 0.0;
  910.       for (i=0;i<4;i++) {
  911.          for (j=0;j<4;j++) {
  912.             temp1 = fabs(point_plane_distance(&((*Patch)[i][j]), &n, &d));
  913.             if (temp1 > dist)
  914.                dist = temp1;
  915.          }
  916.       }
  917.       return dist;
  918.    }
  919.    else {
  920.       /* printf("Subpatch normal failed in determine_subpatch_flatness\n"); */
  921.          return -1.0;
  922.    }
  923. }
  924.  
  925. static int flat_enough(Object, Patch)
  926. BICUBIC_PATCH *Object;
  927. VECTOR (*Patch)[4][4];
  928. {
  929.    DBL Dist;
  930.  
  931.    Dist = determine_subpatch_flatness(Patch);
  932.    if (Dist < 0.0)
  933.       return 0;
  934.    else if (Dist < Object->Flatness_Value)
  935.       return 1;
  936.    else
  937.       return 0;
  938. }
  939.  
  940. static void bezier_subdivider(Ray, Object, Patch, u0, u1, v0, v1,
  941. recursion_depth, depth_count, Depths, U_Values, V_Values)
  942. RAY *Ray;
  943. BICUBIC_PATCH *Object;
  944. VECTOR (*Patch)[4][4];
  945. DBL u0, u1, v0, v1;
  946. int recursion_depth, *depth_count;
  947. DBL *Depths, *U_Values, *V_Values;
  948. {
  949.    VECTOR Lower_Left[4][4], Lower_Right[4][4];
  950.    VECTOR Upper_Left[4][4], Upper_Right[4][4];
  951.    VECTOR center;
  952.    DBL ut, vt, radius;
  953.    int tcnt = Object->Intersection_Count;
  954.  
  955.    /* Don't waste time if there are already too many intersections */
  956.    if (tcnt >= MAX_BICUBIC_INTERSECTIONS) return;
  957.  
  958.    /* Make sure the ray passes through a sphere bounding the control points of
  959.       the patch */
  960.    bezier_bounding_sphere(Patch, ¢er, &radius);
  961.    if (!spherical_bounds_check(Ray, ¢er, radius))
  962.       return;
  963.  
  964.    /* If the patch is close to being flat, then just perform a ray-plane
  965.       intersection test. */
  966.    if (flat_enough(Object, Patch))
  967.       bezier_subpatch_intersect(Ray, Object, Patch, u0, u1, v0,
  968.          recursion_depth+1, depth_count, Depths,
  969.          U_Values, V_Values);
  970.  
  971.    if (recursion_depth >= Object->U_Steps)
  972.       if (recursion_depth >= Object->V_Steps)
  973.          bezier_subpatch_intersect(Ray, Object, Patch, u0, u1, v0,
  974.             recursion_depth+1, depth_count, Depths,
  975.             U_Values, V_Values);
  976.       else {
  977.          bezier_split_up_down(Patch, (VECTOR (*)[4][4])Lower_Left, (VECTOR (*)[4][4])Upper_Left);
  978.          vt = (v1 - v0) / 2.0;
  979.          bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Lower_Left, u0, u1, v0, vt,
  980.             recursion_depth+1, depth_count, Depths,
  981.                U_Values, V_Values);
  982.          bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Upper_Left, u0, u1, vt, v1,
  983.             recursion_depth+1, depth_count, Depths,
  984.             U_Values, V_Values);
  985.       }
  986.    else if (recursion_depth >= Object->V_Steps) {
  987.       bezier_split_left_right(Patch, (VECTOR (*)[4][4])Lower_Left, (VECTOR (*)[4][4])Lower_Right);
  988.       ut = (u1 - u0) / 2.0;
  989.       bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Lower_Left, u0, ut, v0, v1,
  990.          recursion_depth+1, depth_count, Depths,
  991.          U_Values, V_Values);
  992.       bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Lower_Right, ut, u1, v0, v1,
  993.          recursion_depth+1, depth_count, Depths,
  994.          U_Values, V_Values);
  995.    }
  996.    else {
  997.       ut = (u1 - u0) / 2.0;
  998.       vt = (v1 - v0) / 2.0;
  999.       bezier_split_left_right(Patch, (VECTOR (*)[4][4])Lower_Left, (VECTOR (*)[4][4])Lower_Right);
  1000.       bezier_split_up_down((VECTOR (*)[4][4])Lower_Left, (VECTOR (*)[4][4])Lower_Left, (VECTOR (*)[4][4])Upper_Left);
  1001.       bezier_split_up_down((VECTOR (*)[4][4])Lower_Right, (VECTOR (*)[4][4])Lower_Right, (VECTOR (*)[4][4])Upper_Right);
  1002.       bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Lower_Left, u0, ut, v0, vt,
  1003.          recursion_depth+1, depth_count, Depths,
  1004.             U_Values, V_Values);
  1005.       bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Upper_Left, u0, ut, vt, v1,
  1006.          recursion_depth+1, depth_count, Depths,
  1007.          U_Values, V_Values);
  1008.       bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Lower_Right, ut, u1, v0, vt,
  1009.          recursion_depth+1, depth_count, Depths,
  1010.          U_Values, V_Values);
  1011.       bezier_subdivider(Ray, Object, (VECTOR (*)[4][4])Upper_Right, ut, u1, vt, v1,
  1012.          recursion_depth+1, depth_count, Depths,
  1013.          U_Values, V_Values);
  1014.    }
  1015. }
  1016.  
  1017. static void bezier_tree_deleter(Node)
  1018. BEZIER_NODE *Node;
  1019. {
  1020.    BEZIER_CHILDREN *Children;
  1021.    int i;
  1022.  
  1023.    /* If this is an interior node then continue the descent */
  1024.    if (Node->Node_Type == BEZIER_INTERIOR_NODE) {
  1025.       Children = (BEZIER_CHILDREN *)Node->Data_Ptr;
  1026.       for (i=0;i<Node->Count;i++)
  1027.          bezier_tree_deleter(Children->Children[i]);
  1028.       free((void *)Children);
  1029.    }
  1030.    else if (Node->Node_Type == BEZIER_LEAF_NODE) {
  1031.       /* Free the memory used for the vertices. */
  1032.       free((void *)Node->Data_Ptr);
  1033.    }
  1034.    /* Free the memory used for the node. */
  1035.    free((void *)Node);
  1036. }
  1037.  
  1038. static void bezier_tree_walker(Ray, Shape, Node, depth, depth_count, Depths)
  1039. RAY *Ray;
  1040. BICUBIC_PATCH *Shape;
  1041. BEZIER_NODE *Node;
  1042. int depth, *depth_count;
  1043. DBL *Depths;
  1044. {
  1045.    BEZIER_CHILDREN *Children;
  1046.    BEZIER_VERTICES *Vertices;
  1047.    VECTOR n, ip, vv0, vv1, vv2, vv3;
  1048.    DBL d, Depth;
  1049.    int i, tcnt = Shape->Intersection_Count;
  1050.  
  1051.    /* Don't waste time if there are already too many intersections */
  1052.    if (tcnt >= MAX_BICUBIC_INTERSECTIONS) return;
  1053.  
  1054.    /* Make sure the ray passes through a sphere bounding the control points of
  1055.       the patch */
  1056.    if (!spherical_bounds_check(Ray, &(Node->Center), Node->Radius_Squared))
  1057.       return;
  1058.  
  1059.    /* If this is an interior node then continue the descent, else
  1060.       do a check against the vertices. */
  1061.    if (Node->Node_Type == BEZIER_INTERIOR_NODE) {
  1062.       Children = (BEZIER_CHILDREN *)Node->Data_Ptr;
  1063.       for (i=0;i<Node->Count;i++)
  1064.          bezier_tree_walker(Ray, Shape, Children->Children[i],
  1065.             depth+1, depth_count, Depths);
  1066.    }
  1067.    else if (Node->Node_Type == BEZIER_LEAF_NODE) {
  1068.       Vertices = (BEZIER_VERTICES *)Node->Data_Ptr;
  1069.       vv0 = Vertices->Vertices[0]; vv1 = Vertices->Vertices[1];
  1070.       vv2 = Vertices->Vertices[2]; vv3 = Vertices->Vertices[3];
  1071.  
  1072.       /* Triangulate this subpatch, then check for intersections in
  1073.      the triangles. */
  1074.       if (subpatch_normal(&vv0, &vv1, &vv2, &n, &d)) {
  1075.          if (intersect_subpatch(Shape->Patch_Type, Ray, &vv0, &vv1, &vv2, &n, d,
  1076.             NULL, NULL, NULL, &Depth, &ip, &n)) {
  1077.             Shape->Intersection_Point[tcnt + *depth_count] = ip;
  1078.             Shape->Normal_Vector[tcnt + *depth_count] = n;
  1079.             Depths[*depth_count] = Depth;
  1080.             *depth_count += 1;
  1081.          }
  1082.       }
  1083.  
  1084.       if (*depth_count + tcnt >= MAX_BICUBIC_INTERSECTIONS) return;
  1085.  
  1086.       if (subpatch_normal(&vv0, &vv2, &vv3, &n, &d)) {
  1087.          if (intersect_subpatch(Shape->Patch_Type, Ray, &vv0, &vv2, &vv3, &n, d,
  1088.             NULL, NULL, NULL, &Depth, &ip, &n)) {
  1089.             Shape->Intersection_Point[tcnt + *depth_count] = ip;
  1090.             Shape->Normal_Vector[tcnt + *depth_count] = n;
  1091.             Depths[*depth_count] = Depth;
  1092.             *depth_count += 1;
  1093.          }
  1094.       }
  1095.    }
  1096.    else {
  1097.       printf("Bad Node type at depth %d\n", depth);
  1098.    }
  1099. }
  1100.  
  1101.    /* Intersection of a ray and a bezier patch. */
  1102.    /* Note: There is MUCH repeated work being done here. During the subdivision
  1103.    process, the values of surface points are not retained from one subpatch
  1104.    to the next, even though there are two points in common between one subpatch
  1105.    and the one adjacent to it.    An obvious optimization is to retain this
  1106.    surface information during the subdivision process.
  1107.  
  1108.    A second optimization is to make use of the fact that the surface never
  1109.    goes outside the convex hull generated by the control points.  By testing
  1110.    the ray against that hull first many unnecessary tests can be avoided.  The
  1111.    hull really should be generated at the time the patch is first defined -
  1112.    not at run time.
  1113. */
  1114.    static int intersect_bicubic_patch0(Ray, Shape, Depths)
  1115.       RAY *Ray;
  1116. BICUBIC_PATCH *Shape;
  1117. DBL *Depths;
  1118. {
  1119.    int cnt = 0;
  1120.    int tcnt = Shape->Intersection_Count;
  1121.    int i, j;
  1122.    DBL Depth, d, u, v, delta_u, delta_v;
  1123.    VECTOR v0, v1, v2, v3, n, ip;
  1124.    VECTOR (*Patch_Ptr)[4][4] = (VECTOR (*)[4][4]) Shape->Control_Points;
  1125.  
  1126.    if (!spherical_bounds_check(Ray, &(Shape->Bounding_Sphere_Center),
  1127.       Shape->Bounding_Sphere_Radius))
  1128.       return 0;
  1129.  
  1130.    delta_u = 1.0 / (DBL)Shape->U_Steps;
  1131.    delta_v = 1.0 / (DBL)Shape->V_Steps;
  1132.  
  1133.    /* Calculate the initial point */
  1134.    for (i=0;i<Shape->U_Steps;i++) {
  1135.       u = (DBL)i / (DBL)Shape->U_Steps;
  1136.       for (j=0;j<Shape->V_Steps;j++) {
  1137.          v = (DBL)j / (DBL)Shape->V_Steps;
  1138.  
  1139.          /* Calculate surface values for the current patch. */
  1140.          bezier_value(&v0, u, v, Patch_Ptr);
  1141.          bezier_value(&v1, u+delta_u, v, Patch_Ptr);
  1142.          bezier_value(&v2, u, v+delta_v, Patch_Ptr);
  1143.          bezier_value(&v3, u+delta_u, v+delta_v, Patch_Ptr);
  1144.  
  1145.          /* Triangulate this subpatch, then check for intersections in
  1146.         the triangles. */
  1147.          if (subpatch_normal(&v0, &v2, &v1, &n, &d)) {
  1148.             if (intersect_subpatch(Shape->Patch_Type, Ray, &v0, &v2, &v1, &n, d,
  1149.                NULL, NULL, NULL, &Depth, &ip, &n)) {
  1150.                Shape->Intersection_Point[tcnt+cnt] = ip;
  1151.                Shape->Normal_Vector[tcnt+cnt] = n;
  1152.                Depths[cnt] = Depth;
  1153.                if (tcnt + ++cnt >= MAX_BICUBIC_INTERSECTIONS)
  1154.                   /* Too many intersections. Stop looking for more. */
  1155.                   return cnt;
  1156.             }
  1157.          }
  1158.          if (subpatch_normal(&v1, &v2, &v3, &n, &d)) {
  1159.             if (intersect_subpatch(Shape->Patch_Type, Ray, &v1, &v2, &v3, &n, d,
  1160.                NULL, NULL, NULL, &Depth, &ip, &n)) {
  1161.                Shape->Intersection_Point[tcnt+cnt] = ip;
  1162.                Shape->Normal_Vector[tcnt+cnt] = n;
  1163.                Depths[cnt] = Depth;
  1164.                if (tcnt + ++cnt >= MAX_BICUBIC_INTERSECTIONS)
  1165.                   /* Too many intersections. Stop looking for more. */
  1166.                   return cnt;
  1167.             }
  1168.          }
  1169.       }
  1170.    }
  1171.    return cnt;
  1172. }
  1173.  
  1174. static int intersect_bicubic_patch1(Ray, Shape, Depths)
  1175. RAY *Ray;
  1176. BICUBIC_PATCH *Shape;
  1177. DBL *Depths;
  1178. {
  1179.    int cnt = 0;
  1180.    int tcnt = Shape->Intersection_Count;
  1181.    int i, j;
  1182.    DBL Depth, d, radius;
  1183.    VECTOR v[4], n, ip, center;
  1184.  
  1185.    if (!spherical_bounds_check(Ray, &(Shape->Bounding_Sphere_Center),
  1186.       Shape->Bounding_Sphere_Radius))
  1187.       return 0;
  1188.  
  1189.    /* Calculate the initial point */
  1190.    for (i=0;i<Shape->U_Steps;i++) {
  1191.       for (j=0;j<Shape->V_Steps;j++) {
  1192.  
  1193.          /* Grab precomputed surface values for the current patch. */
  1194.          v[0] = Shape->Interpolated_Grid[i][j];
  1195.          v[1] = Shape->Interpolated_Grid[i+1][j];
  1196.          v[2] = Shape->Interpolated_Grid[i][j+1];
  1197.          v[3] = Shape->Interpolated_Grid[i+1][j+1];
  1198.  
  1199.          /* Check the ray against the bounding sphere for this subpatch */
  1200.          find_average(4, &v[0], ¢er, &radius);
  1201.          if (!spherical_bounds_check(Ray, ¢er, radius)) continue;
  1202.  
  1203.          n = Shape->Interpolated_Normals[i][2*j];
  1204.          if (n.x == 0.0 && n.y == 0.0 && n.z == 0.0) {
  1205.             goto l0;
  1206.          }
  1207.          d = Shape->Interpolated_D[i][2*j];
  1208.  
  1209.          /* Check for intersections in this subpatch. */
  1210.          if (intersect_subpatch(Shape->Patch_Type, Ray, &v[0], &v[2], &v[1], &n, d,
  1211.             NULL, NULL, NULL, &Depth, &ip, &n)) {
  1212.             Shape->Intersection_Point[tcnt+cnt] = ip;
  1213.             Shape->Normal_Vector[tcnt+cnt] = n;
  1214.             Depths[cnt] = Depth;
  1215.             if (tcnt + ++cnt >= MAX_BICUBIC_INTERSECTIONS)
  1216.                /* Too many intersections. Stop looking for more. */
  1217.                return cnt;
  1218.          }
  1219. l0:
  1220.          n = Shape->Interpolated_Normals[i][2*j+1];
  1221.          if (n.x == 0.0 && n.y == 0.0 && n.z == 0.0) {
  1222.             continue;
  1223.          }
  1224.          d = Shape->Interpolated_D[i][2*j+1];
  1225.          if (intersect_subpatch(Shape->Patch_Type, Ray, &v[1], &v[2], &v[3], &n, d,
  1226.             NULL, NULL, NULL, &Depth, &ip, &n)) {
  1227.             Shape->Intersection_Point[tcnt+cnt] = ip;
  1228.             Shape->Normal_Vector[tcnt+cnt] = n;
  1229.             Depths[cnt] = Depth;
  1230.             if (tcnt + ++cnt >= MAX_BICUBIC_INTERSECTIONS)
  1231.                /* Too many intersections. Stop looking for more. */
  1232.                return cnt;
  1233.          }
  1234.       }
  1235.    }
  1236.    return cnt;
  1237. }
  1238.  
  1239. static int intersect_bicubic_patch2(Ray, Shape, Depths)
  1240. RAY *Ray;
  1241. BICUBIC_PATCH *Shape;
  1242. DBL *Depths;
  1243. {
  1244.    int cnt = 0;
  1245.    DBL U_Values[MAX_BICUBIC_INTERSECTIONS];
  1246.    DBL V_Values[MAX_BICUBIC_INTERSECTIONS];
  1247.    VECTOR (*Patch)[4][4] = (VECTOR (*)[4][4]) Shape->Control_Points;
  1248.  
  1249.    bezier_subdivider(Ray, Shape, Patch, 0.0, 1.0, 0.0, 1.0,
  1250.       0, &cnt, Depths, &U_Values[0], &V_Values[0]);
  1251.    return cnt;
  1252. }
  1253.  
  1254. static int intersect_bicubic_patch3(Ray, Shape, Depths)
  1255. RAY *Ray;
  1256. BICUBIC_PATCH *Shape;
  1257. DBL *Depths;
  1258. {
  1259.    int cnt = 0;
  1260.    bezier_tree_walker(Ray, Shape, Shape->Node_Tree, 0, &cnt, Depths);
  1261.    return cnt;
  1262. }
  1263.  
  1264. static int intersect_bicubic_patch4(Ray, Shape, Depths)
  1265. RAY *Ray;
  1266. BICUBIC_PATCH *Shape;
  1267. DBL *Depths;
  1268. {
  1269.    int cnt = 0;
  1270.    int tcnt = Shape->Intersection_Count;
  1271.    int i, j;
  1272.    DBL Depth, d, t;
  1273.    VECTOR v0, v1, v2, v3, n;
  1274.    VECTOR n0, n1, n2, n3, ip, ip_norm;
  1275.  
  1276.    if (!spherical_bounds_check(Ray, &(Shape->Bounding_Sphere_Center),
  1277.       Shape->Bounding_Sphere_Radius))
  1278.       return 0;
  1279.  
  1280.    /* Calculate the initial point */
  1281.    for (i=0;i<Shape->U_Steps;i++) {
  1282.       for (j=0;j<Shape->V_Steps;j++) {
  1283.  
  1284.          /* Grab precomputed surface values for the current patch. */
  1285.          v0 = Shape->Interpolated_Grid[i][j];
  1286.          v1 = Shape->Interpolated_Grid[i+1][j];
  1287.          v2 = Shape->Interpolated_Grid[i][j+1];
  1288.          v3 = Shape->Interpolated_Grid[i+1][j+1];
  1289.  
  1290.          n0 = Shape->Smooth_Normals[i][j];
  1291.          n1 = Shape->Smooth_Normals[i+1][j];
  1292.          n2 = Shape->Smooth_Normals[i][j+1];
  1293.          n3 = Shape->Smooth_Normals[i+1][j+1];
  1294.  
  1295.          n = Shape->Interpolated_Normals[i][2*j];
  1296.          if (n.x == 0.0 && n.y == 0.0 && n.z == 0.0) goto l0;
  1297.          d = Shape->Interpolated_D[i][2*j];
  1298.  
  1299.          /* Make sure the smooth normals point in the same direction as the normal */
  1300.          VDot(t, n0, n); if (t < 0) VScale(n0, n0, -1.0);
  1301.          VDot(t, n1, n); if (t < 0) VScale(n1, n1, -1.0);
  1302.          VDot(t, n2, n); if (t < 0) VScale(n2, n2, -1.0);
  1303.  
  1304.          /* Check for intersections in this subpatch. */
  1305.          if (intersect_subpatch(Shape->Patch_Type, Ray, &v0, &v2, &v1, &n,
  1306.             d, &n0, &n2, &n1, &Depth, &ip, &ip_norm)) {
  1307.             Shape->Intersection_Point[tcnt + cnt] = ip;
  1308.             Shape->Normal_Vector[tcnt + cnt] = ip_norm;
  1309.             Depths[cnt] = Depth;
  1310.             if (tcnt + ++cnt >= MAX_BICUBIC_INTERSECTIONS)
  1311.                /* Too many intersections. Stop looking for more. */
  1312.                return cnt;
  1313.          }
  1314. l0:
  1315.          n = Shape->Interpolated_Normals[i][2*j+1];
  1316.          if (n.x == 0.0 && n.y == 0.0 && n.z == 0.0) continue;
  1317.          d = Shape->Interpolated_D[i][2*j+1];
  1318.  
  1319.          /* Make sure the smooth normals point in the same direction as the normal */
  1320.          VDot(t, n1, n); if (t > 0) VScale(n1, n0, -1.0);
  1321.          VDot(t, n2, n); if (t > 0) VScale(n2, n1, -1.0);
  1322.          VDot(t, n3, n); if (t > 0) VScale(n3, n2, -1.0);
  1323.  
  1324.          if (intersect_subpatch(Shape->Patch_Type, Ray, &v1, &v2, &v3, &n,
  1325.             d, &n1, &n2, &n3, &Depth, &ip, &ip_norm)) {
  1326.             Shape->Intersection_Point[tcnt + cnt] = ip;
  1327.             Shape->Normal_Vector[tcnt + cnt] = ip_norm;
  1328.             Depths[cnt] = Depth;
  1329.             if (tcnt + ++cnt >= MAX_BICUBIC_INTERSECTIONS)
  1330.                /* Too many intersections. Stop looking for more. */
  1331.                return cnt;
  1332.          }
  1333.       }
  1334.    }
  1335.    return cnt;
  1336. }
  1337.  
  1338. int All_Bicubic_Patch_Intersections(Object, Ray, Depth_Queue)
  1339. OBJECT *Object;
  1340. RAY *Ray;
  1341. PRIOQ *Depth_Queue;
  1342. {
  1343.    BICUBIC_PATCH *Shape = (BICUBIC_PATCH *) Object;
  1344.    DBL Depths[MAX_BICUBIC_INTERSECTIONS];
  1345.    INTERSECTION Local_Element;
  1346.    int cnt, tcnt, i, Intersection_Found;
  1347.  
  1348.    Intersection_Found = 0;
  1349.    Ray_Bicubic_Tests++;
  1350.    if (Ray == VP_Ray)
  1351.       Shape->Intersection_Count = 0;
  1352.    tcnt = Shape->Intersection_Count;
  1353.    if (Shape->Patch_Type == 0)
  1354.       cnt = intersect_bicubic_patch0(Ray, Shape, &Depths[0]);
  1355.    else if (Shape->Patch_Type == 1)
  1356.       cnt = intersect_bicubic_patch1(Ray, Shape, &Depths[0]);
  1357.    else if (Shape->Patch_Type == 2)
  1358.       cnt = intersect_bicubic_patch2(Ray, Shape, &Depths[0]);
  1359.    else if (Shape->Patch_Type == 3)
  1360.       cnt = intersect_bicubic_patch3(Ray, Shape, &Depths[0]);
  1361.    else if (Shape->Patch_Type == 4)
  1362.       cnt = intersect_bicubic_patch4(Ray, Shape, &Depths[0]);
  1363.    else {
  1364.       Error("Bad patch type\n");
  1365.    }
  1366.    if (cnt > 0) Ray_Bicubic_Tests_Succeeded++;
  1367.    for (i=0;i<cnt;i++) {
  1368.       if (!Shadow_Test_Flag)
  1369.          Shape->Intersection_Count++;
  1370.       Local_Element.Depth = Depths[i];
  1371.       Local_Element.Object = Shape->Parent_Object;
  1372.       Local_Element.Point = Shape->Intersection_Point[tcnt + i];
  1373.       Local_Element.Shape = (SHAPE *)Shape;
  1374.       pq_add (Depth_Queue, &Local_Element);
  1375.       Intersection_Found = 1;
  1376.    }
  1377.    return (Intersection_Found);
  1378. }
  1379.  
  1380. /* A patch is not a solid, so an inside test doesn't make sense. */
  1381. int Inside_Bicubic_Patch (Test_Point, Object)
  1382. VECTOR *Test_Point;
  1383. OBJECT *Object;
  1384. {
  1385.    return 0;
  1386. }
  1387.  
  1388. void Bicubic_Patch_Normal (Result, Object, Intersection_Point)
  1389. OBJECT *Object;
  1390. VECTOR *Result, *Intersection_Point;
  1391. {
  1392.    BICUBIC_PATCH *Patch = (BICUBIC_PATCH *)Object;
  1393.    int i;
  1394.  
  1395.    /* If all is going well, the normal was computed at the time the intersection
  1396.       was computed.  Look on the list of associated intersection points and normals */
  1397.    for (i=0;i<Patch->Intersection_Count;i++)
  1398.       if (Intersection_Point->x == Patch->Intersection_Point[i].x &&
  1399.          Intersection_Point->y == Patch->Intersection_Point[i].y &&
  1400.          Intersection_Point->z == Patch->Intersection_Point[i].z) {
  1401.          Result->x = Patch->Normal_Vector[i].x;
  1402.          Result->y = Patch->Normal_Vector[i].y;
  1403.          Result->z = Patch->Normal_Vector[i].z;
  1404.          return;
  1405.       }
  1406.    if (Options & DEBUGGING) {
  1407.       printf("Bicubic patch normal for unknown intersection point\n");
  1408.       fflush(stdout);
  1409.    }
  1410.    Result->x = 1.0;
  1411.    Result->y = 0.0;
  1412.    Result->z = 0.0;
  1413. }
  1414.  
  1415. void *Copy_Bicubic_Patch (Object)
  1416. OBJECT *Object;
  1417. {
  1418.    BICUBIC_PATCH *New_Shape;
  1419.  
  1420.    New_Shape = Get_Bicubic_Patch_Shape ();
  1421.    *New_Shape = * ((BICUBIC_PATCH *)Object);
  1422.    New_Shape -> Next_Object = NULL;
  1423.  
  1424.    New_Shape->Interpolated_Grid = NULL;
  1425.    Precompute_Patch_Values(New_Shape);
  1426.    if (New_Shape->Shape_Texture != NULL)
  1427.       New_Shape->Shape_Texture = Copy_Texture (New_Shape->Shape_Texture);
  1428.  
  1429.    return (void *)(New_Shape);
  1430. }
  1431.  
  1432. void Translate_Bicubic_Patch (Object, Vector)
  1433. OBJECT *Object;
  1434. VECTOR *Vector;
  1435. {
  1436.    BICUBIC_PATCH *Patch = (BICUBIC_PATCH *) Object;
  1437.    int i, j;
  1438.    for (i=0;i<4;i++) for (j=0;j<4;j++)
  1439.       VAdd (Patch->Control_Points[i][j], Patch->Control_Points[i][j], *Vector)
  1440.          Precompute_Patch_Values(Patch);
  1441.    Translate_Texture (&((BICUBIC_PATCH *) Object)->Shape_Texture, Vector);
  1442. }
  1443.  
  1444. void Rotate_Bicubic_Patch (Object, Vector)
  1445. OBJECT *Object;
  1446. VECTOR *Vector;
  1447. {
  1448.    TRANSFORMATION Transformation;
  1449.    BICUBIC_PATCH *Patch = (BICUBIC_PATCH *) Object;
  1450.    int i, j;
  1451.  
  1452.    Get_Rotation_Transformation (&Transformation, Vector);
  1453.    for (i=0;i<4;i++) for (j=0;j<4;j++)
  1454.       MTransformVector(&(Patch->Control_Points[i][j]),
  1455.          &(Patch->Control_Points[i][j]),
  1456.          &Transformation);
  1457.    Precompute_Patch_Values(Patch);
  1458.    Rotate_Texture (&((BICUBIC_PATCH *) Object)->Shape_Texture, Vector);
  1459. }
  1460.  
  1461. void Scale_Bicubic_Patch (Object, Vector)
  1462. OBJECT *Object;
  1463. VECTOR *Vector;
  1464. {
  1465.    BICUBIC_PATCH *Patch = (BICUBIC_PATCH *) Object;
  1466.    int i, j;
  1467.    for (i=0;i<4;i++) for (j=0;j<4;j++)
  1468.       VEvaluate(Patch->Control_Points[i][j],
  1469.          Patch->Control_Points[i][j], *Vector);
  1470.    Precompute_Patch_Values(Patch);
  1471.    Scale_Texture (&((BICUBIC_PATCH *) Object)->Shape_Texture, Vector);
  1472. }
  1473.  
  1474.  
  1475. /* Inversion of a patch really doesn't make sense. */
  1476. void Invert_Bicubic_Patch (Object)
  1477. OBJECT *Object;
  1478. {
  1479.    ;
  1480. }
  1481.