home *** CD-ROM | disk | FTP | other *** search
/ gondwana.ecr.mu.oz.au/pub/ / Graphics.tar / Graphics / spline-patch.tar.gz / spline-patch.tar / patch / patch.c < prev    next >
C/C++ Source or Header  |  1991-11-18  |  18KB  |  583 lines
  1. /*
  2.  *   File     : patch.c
  3.  *   Author   : Sean Graves
  4.  *
  5.  *   Computer : Any 
  6.  *   Date     : 6/21/90 
  7.  *
  8.  *   A set of routines for creating and raytracing bicubic surface
  9.  *   patches.  The file patch.h contains function prototypes for
  10.  *   the routines in this file.
  11.  *
  12.  *   User Routines:
  13.  *     NewPatch - Allocates the patch data structures and returns a pointer.
  14.  *     DefPatch - Defines a patch, given a patch pointer.
  15.  *     IsectPatch - Intersects a ray and a patch. Returns u, v, t, pnt. 
  16.  *     NormalPatch - Given a patch and a u, v position, returns the normal.
  17.  *     FreePatch - Deallocates a patch.
  18.  *     
  19.  *   Additionally, there is one useful routine not specific to patches:
  20.  *     BoxIntersect - Given a ray and a box, returns the intersection point.
  21.  *
  22.  * Copyright (c) 1990, by Sean Graves and Texas A&M University 
  23.  *
  24.  * Permission is hereby granted for non-commercial reproduction and use of
  25.  * this program, provided that this notice is included in any material copied
  26.  * from it. The author assumes no responsibility for damages resulting from
  27.  * the use of this software, however caused.
  28.  *
  29.  */
  30.  
  31. #include "patch.h"
  32. #include "vector.h"
  33. #include "math.h"
  34.  
  35. #define NEWTON_TOL 0.01
  36. #define MAX_ITER 3
  37. #define OVERLAP 0.1
  38. #define INF (1e99)
  39. #define MAXDIV 16
  40. #define TRUE 1
  41. #define FALSE 0
  42. #define NULL 0 
  43.  
  44. /*  Function Prototypes for local functions */
  45.  
  46. static void samplepatch();
  47. static compute_overlap();
  48. static TreePtr buildtree();
  49. static ListPtr insert_node(); 
  50. static void DumpList();          /* Diagnostic */
  51. static void free_list();
  52. static float eval_patch();
  53. static void newton();   
  54.  
  55. static POINT samples[MAXDIV +1][MAXDIV +1]; 
  56.  
  57. /*************************************************************************/
  58.  
  59. PatchPtr NewPatch()
  60. {
  61.    return((PatchPtr) malloc(sizeof(Patch)));
  62. }
  63.  
  64. /*************************************************************************/
  65.  
  66. void FreePatch(p)
  67. PatchPtr p;
  68. {
  69.    free(p);
  70. }
  71.  
  72. /*************************************************************************/
  73.  
  74. void DefPatch(p, M, geomx, geomy, geomz, div)
  75. PatchPtr p;                      /* Pointer to patch struct */
  76. MATRIX M, geomx, geomy, geomz;   /* Spline basis, geometry matrices */
  77. int div;                         /* Number of divisions for tree */
  78. {
  79.    int i, j;
  80.  
  81.  
  82.    /* Copy basic information into patch struct */
  83.  
  84.    for (i = 0; i < 4; i++) 
  85.       for (j = 0; j < 4; j++) { 
  86.          p->M[i][j] = M[i][j];
  87.          p->geomx[i][j] = geomx[i][j];
  88.          p->geomy[i][j] = geomy[i][j];
  89.          p->geomz[i][j] = geomz[i][j];
  90.       }
  91.  
  92.    /* Compute derivative of basis matrix */
  93.  
  94.    p->Md[0][0] = p->Md[0][1] = p->Md[0][2] = p->Md[0][3] = 0.0;
  95.    for (i = 1; i < 4; i++) 
  96.       for (j = 0; j < 4; j++) 
  97.          p->Md[i][j] = p->M[i-1][j] * (4-i); 
  98.        
  99.    /* Fill in transpose fields of patch struct */
  100.  
  101.    transpose(p->M , p->MT );
  102.    transpose(p->Md, p->MdT);
  103.  
  104.    /* Compute intermediate blending matrices */
  105.  
  106.    mmult4x4_4x4(p->Mx, p->M, p->geomx); /* Calculate blending matrix for x,y,z */
  107.    mmult4x4_4x4(p->My, p->M, p->geomy);
  108.    mmult4x4_4x4(p->Mz, p->M, p->geomz);
  109.    mmult4x4_4x4(p->Mx, p->Mx, p->MT);
  110.    mmult4x4_4x4(p->My, p->My, p->MT);
  111.    mmult4x4_4x4(p->Mz, p->Mz, p->MT);
  112.  
  113.    mmult4x4_4x4(p->Mxdu, p->Md, p->geomx);  /* Calculate d/du matrix for x,y,z */
  114.    mmult4x4_4x4(p->Mydu, p->Md, p->geomy);
  115.    mmult4x4_4x4(p->Mzdu, p->Md, p->geomz);
  116.    mmult4x4_4x4(p->Mxdu, p->Mxdu, p->MT);
  117.    mmult4x4_4x4(p->Mydu, p->Mydu, p->MT);
  118.    mmult4x4_4x4(p->Mzdu, p->Mzdu, p->MT);
  119.  
  120.    mmult4x4_4x4(p->Mxdv, p->M, p->geomx);  /* Calculate d/dv matrix for x,y,z */
  121.    mmult4x4_4x4(p->Mydv, p->M, p->geomy);
  122.    mmult4x4_4x4(p->Mzdv, p->M, p->geomz);
  123.    mmult4x4_4x4(p->Mxdv, p->Mxdv, p->MdT);
  124.    mmult4x4_4x4(p->Mydv, p->Mydv, p->MdT);
  125.    mmult4x4_4x4(p->Mzdv, p->Mzdv, p->MdT);
  126.  
  127.    samplepatch(p, div);
  128.    p->tree = buildtree(0, div, 0, div, div);
  129.  
  130.  }
  131.     
  132. /*******************************************************************/
  133.  
  134. static void samplepatch(p, div)
  135. PatchPtr p;
  136. int div;
  137. {
  138.    float U[1][4], V[4][1];         /* u, v vectors */
  139.    MATRIX C;                       /* sol matrix (only [0][0] used) */
  140.    MATRIX M2x, M2y, M2z;           /* Intermediate results */
  141.    float u, v;                     /* Parameters */
  142.    int i,j;
  143.  
  144.    for (i = 0; i <= div; i++) {
  145.       u = (float) i / div;
  146.       U[0][3] = 1;
  147.       U[0][2] = u;
  148.       U[0][1] = u * u;
  149.       U[0][0] = U[0][1] * u;
  150.       mmult1x4_4x4 (M2x, U, p->Mx);
  151.       mmult1x4_4x4 (M2y, U, p->My);
  152.       mmult1x4_4x4 (M2z, U, p->Mz);
  153.       for (j = 0; j <= div; j++) {
  154.         v = (float) j / div;
  155.         V[3][0] = 1;
  156.         V[2][0] = v;
  157.         V[1][0] = v * v;
  158.         V[0][0] = V[1][0] * v;
  159.         mmult1x4_4x1 (C, M2x, V);
  160.         samples[i][j].x = C[0][0];
  161.         mmult1x4_4x1 (C, M2y, V);
  162.         samples[i][j].y = C[0][0];
  163.         mmult1x4_4x1 (C, M2z, V);
  164.         samples[i][j].z = C[0][0];
  165.       }
  166.    }
  167. }
  168.  
  169. /***************************************************************************/
  170.  
  171. static compute_overlap(curr_node)
  172. TreePtr curr_node;
  173. {
  174.    int i;
  175.    float incrs;
  176. /*
  177.    for (i=0; i<3; i++) {
  178.       incrs = fabs((curr_node->box_max[i] - curr_node->box_min[i]) * OVERLAP);
  179.       curr_node->box_min[i] -= incrs;
  180.       curr_node->box_max[i] += incrs;
  181.    } 
  182. */
  183. }
  184.  
  185. /***************************************************************************/
  186.  
  187. static TreePtr buildtree(u_low, u_hi, v_low, v_hi, div)
  188. int u_low, u_hi, v_low, v_hi, div;
  189. {
  190.    int u_mid, v_mid, i, j;
  191.    TreePtr curr_node;
  192.    
  193.    curr_node = (TreePtr) malloc(sizeof(TreeNode));
  194.    if (!curr_node) { printf("Error - malloc failed\n"); exit(0); }
  195.    if ((u_low+1 == u_hi) || (v_low+1 == v_hi)) {         /* Found leaf */
  196.                                                        /* Compute midpoint */
  197.       curr_node->u_mid = ((((float) (u_hi - u_low)) / 2.0) + u_low) / div ; 
  198.       curr_node->v_mid = ((((float) (v_hi - v_low)) / 2.0) + v_low) / div ; 
  199.                                                        /* Compute box */
  200.       if (length(samples[u_low][v_low], samples[u_hi][v_hi]) >
  201.           length(samples[u_hi][v_low], samples[u_low][v_hi])) {
  202.          curr_node->box_min.x = samples[u_low][v_low].x;
  203.          curr_node->box_max.x = samples[u_hi ][v_hi ].x;
  204.          curr_node->box_min.y = samples[u_low][v_low].y;
  205.          curr_node->box_max.y = samples[u_hi ][v_hi ].y;
  206.          curr_node->box_min.z = samples[u_low][v_low].z;
  207.          curr_node->box_max.z = samples[u_hi ][v_hi ].z;
  208.       } else {
  209.          curr_node->box_min.x = samples[u_hi ][v_low].x;
  210.          curr_node->box_max.x = samples[u_low][v_hi ].x;
  211.          curr_node->box_min.y = samples[u_hi ][v_low].y;
  212.          curr_node->box_max.y = samples[u_low][v_hi ].y;
  213.          curr_node->box_min.z = samples[u_hi ][v_low].z;
  214.          curr_node->box_max.z = samples[u_low][v_hi ].z;
  215.       }
  216.       for (j = 0; j < 4; j++) curr_node->child[j] = NULL;
  217.  
  218.    } else {                                           /* Not leaf node */ 
  219.       u_mid = (u_hi - u_low) / 2 + u_low; 
  220.       v_mid = (v_hi - v_low) / 2 + v_low;
  221.  
  222.       curr_node->child[0] = buildtree(u_low, u_mid, v_low, v_mid, div); 
  223.       curr_node->child[1] = buildtree(u_mid, u_hi, v_low, v_mid, div); 
  224.       curr_node->child[2] = buildtree(u_low, u_mid, v_mid, v_hi, div);
  225.       curr_node->child[3] = buildtree(u_mid, u_hi, v_mid, v_hi, div);
  226.                                                       /* Compute box */
  227.       curr_node->box_min.x = 1e99;
  228.       curr_node->box_max.x = -1e99;
  229.       for (j = 0; j < 4; j++) {;
  230.          if (curr_node->child[j]->box_min.x < curr_node->box_min.x)
  231.             curr_node->box_min.x = curr_node->child[j]->box_min.x;
  232.          if (curr_node->child[j]->box_min.x > curr_node->box_max.x)
  233.             curr_node->box_max.x = curr_node->child[j]->box_min.x;
  234.          if (curr_node->child[j]->box_max.x < curr_node->box_min.x)
  235.             curr_node->box_min.x = curr_node->child[j]->box_max.x;
  236.          if (curr_node->child[j]->box_max.x > curr_node->box_max.x)
  237.             curr_node->box_max.x = curr_node->child[j]->box_max.x;
  238.       }
  239.       curr_node->box_min.y = 1e99;
  240.       curr_node->box_max.y = -1e99;
  241.       for (j = 0; j < 4; j++) {;
  242.          if (curr_node->child[j]->box_min.y < curr_node->box_min.y)
  243.             curr_node->box_min.y = curr_node->child[j]->box_min.y;
  244.          if (curr_node->child[j]->box_min.y > curr_node->box_max.y)
  245.             curr_node->box_max.y = curr_node->child[j]->box_min.y;
  246.          if (curr_node->child[j]->box_max.y < curr_node->box_min.y)
  247.             curr_node->box_min.y = curr_node->child[j]->box_max.y;
  248.          if (curr_node->child[j]->box_max.y > curr_node->box_max.y)
  249.             curr_node->box_max.y = curr_node->child[j]->box_max.y;
  250.       }
  251.       curr_node->box_min.z = 1e99;
  252.       curr_node->box_max.z = -1e99;
  253.       for (j = 0; j < 4; j++) {;
  254.          if (curr_node->child[j]->box_min.z < curr_node->box_min.z)
  255.             curr_node->box_min.z = curr_node->child[j]->box_min.z;
  256.          if (curr_node->child[j]->box_min.z > curr_node->box_max.z)
  257.             curr_node->box_max.z = curr_node->child[j]->box_min.z;
  258.          if (curr_node->child[j]->box_max.z < curr_node->box_min.z)
  259.             curr_node->box_min.z = curr_node->child[j]->box_max.z;
  260.          if (curr_node->child[j]->box_max.z > curr_node->box_max.z)
  261.             curr_node->box_max.z = curr_node->child[j]->box_max.z;
  262.       }
  263.    }
  264. /*   compute_overlap(curr_node); */
  265.    return curr_node;
  266. }
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  
  273.  
  274. /*************************************************************************/
  275.  
  276. float BoxIntersect(r, b1, b2) 
  277. RAY r;
  278. POINT b1, b2;
  279. {
  280.    float t1, t2, tswap, tnear, tfar, dinv;
  281.  
  282.    /*  Check X slabs */
  283.    if (fcmp (r.dir.x, 0.0)) {  /* Ray is parallel to X planes of box */
  284.       if (!(((b1.x <= r.org.x ) && (r.org.x <= b2.x)) ||
  285.             ((b2.x <= r.org.x ) && (r.org.x <= b1.x))))/* But not between */
  286.          return -1.0;
  287.       tfar = INF;
  288.       tnear = -INF;
  289.    } else {
  290.       dinv = 1.0/r.dir.x;         /* Calculate inverse for speed */
  291.       t1 = (b1.x - r.org.x) * dinv;
  292.       t2 = (b2.x - r.org.x) * dinv;
  293.  
  294.       if (t1 > t2) { tfar = t1; tnear = t2; } 
  295.       else         { tfar = t2; tnear = t1; } 
  296.       if (tfar < 0.0)              /* Box is behind ray */
  297.          return -1.0;
  298.    }
  299.  
  300.    /* Check Y slabs */
  301.  
  302.    if (fcmp (r.dir.y, 0.0)) {   /* Ray is parallel to Y planes of box */
  303.       if (!(((b1.y <= r.org.y ) && (r.org.y <= b2.y)) ||
  304.             ((b2.y <= r.org.y ) && (r.org.y <= b1.y))))/* But not between */
  305.          return -1.0;
  306.    } else {
  307.       dinv = 1.0/r.dir.y;         /* Calculate inverse for speed */
  308.       t1 = (b1.y - r.org.y) * dinv;
  309.       t2 = (b2.y - r.org.y) * dinv;
  310.  
  311.       if (t1 > t2) { tswap = t1; t1 = t2; t2 = tswap; }
  312.  
  313.       if (t1 > tnear) tnear = t1;
  314.       if (t2 < tfar)  tfar = t2;
  315.  
  316.       if (tnear > tfar)            /* Box is missed */
  317.          return -1.0;
  318.  
  319.       if (tfar < 0.0)              /* Box is behind ray */
  320.       return -1.0;     
  321.    }
  322.  
  323.    /* Check Z slabs */
  324.  
  325.    if (fcmp (r.dir.z, 0.0)) {   /* Ray is parallel to Z planes of box */
  326.       if (!(((b1.z <= r.org.z ) && (r.org.z <= b2.z)) ||
  327.             ((b2.z <= r.org.z ) && (r.org.z <= b1.z))))/* But not between */
  328.          return -1.0;
  329.    } else {
  330.       dinv = 1.0/r.dir.z;         /* Calculate inverse for speed */
  331.       t1 = (b1.z - r.org.z) * dinv;
  332.       t2 = (b2.z - r.org.z) * dinv;
  333.  
  334.       if (t1 > t2) { tswap = t1; t1 = t2; t2 = tswap; }
  335.  
  336.       if (t1 > tnear) tnear = t1;
  337.       if (t2 < tfar)  tfar = t2;
  338.  
  339.       if (tnear > tfar)            /* Box is missed */
  340.          return -1.0;
  341.  
  342.       if (tfar < 0.0)              /* Box is behind ray */
  343.          return -1.0;     
  344.    } 
  345.    return tnear;
  346. } 
  347.   
  348. /*************************************************************************/
  349.  
  350. /*************************************************************************/
  351. /* Create a intersection list node.  Put it in the list based on the
  352. /* parameter t.  Return the new head of the list
  353. /*  Pass: list ptr, node to insert, t
  354. /************************************************************************/
  355.  
  356. static ListPtr insert_node(list, node, t)
  357. ListPtr list;
  358. TreePtr node;
  359. float t;
  360. {
  361.    ListPtr new, this, back;
  362.  
  363.    new = (ListPtr) malloc(sizeof(ListNode));
  364.    if (!new) { printf("Error - malloc failed\n"); exit(0); }
  365.    new->t = t;
  366.    new->node = node;
  367.    new->next = NULL;
  368.  
  369.    if (!list) 
  370.       return new;
  371.  
  372.    back = NULL;
  373.    this = list;
  374.    while (this && (this->t < new->t)) {
  375.       back = this;
  376.       this = this->next;
  377.    }
  378.    if (!back) {                       /* First item */
  379.      new->next = list;
  380.      return new;
  381.    } else {                             /* Second through last item */
  382.      back->next = new;
  383.      new->next  = this;
  384.      return list;
  385.    }
  386. }   
  387.    
  388. /************************************************************************/
  389.  
  390. static void dump_list(list) 
  391. ListPtr list;
  392. {
  393.    ListPtr this;
  394.  
  395.    this = list;
  396.    while (this) {
  397.       printf("%4.2f  ",this->t);
  398.       this = this->next;
  399.    }
  400.    printf("\n");
  401. }
  402.  
  403. /************************************************************************/
  404.  
  405. static void free_list(list) 
  406. ListPtr list;
  407. {
  408.    ListPtr this;
  409.  
  410.    this = list;
  411.    while (this) {
  412.       list = this->next;
  413.       free(this);
  414.       this = list;
  415.    }
  416. }
  417.  
  418. /************************************************************************/
  419.  
  420. static float eval_patch(u, M, G, MT, v)
  421. float u, v;
  422. MATRIX M, G, MT;
  423. {
  424.    float U[1][4], V[4][1];              /* u, v vectors */
  425.    MATRIX M2;                           /* Intermediate results */
  426.    int i,j;
  427.  
  428.    U[0][3] = 1;
  429.    U[0][2] = u;
  430.    U[0][1] = u * u;
  431.    U[0][0] = U[0][1] * u;
  432.    V[3][0] = 1;
  433.    V[2][0] = v;
  434.    V[1][0] = v * v;
  435.    V[0][0] = V[1][0] * v;
  436.  
  437.    mmult1x4_4x4 (M2, U, M);
  438.    mmult1x4_4x4 (M2, M2, G);
  439.    mmult1x4_4x4 (M2, M2, MT);
  440.    mmult1x4_4x1 (M2, M2, V); 
  441.    return M2[0][0];
  442. }
  443.  
  444. /************************************************************************/
  445.  
  446. static void newton(surface, node, r, irec)
  447. PatchPtr surface;
  448. TreePtr node;
  449. RAY r;
  450. Isectrec *irec;
  451. {
  452.    int iterations, success, failure, i, j;
  453.    float u, v, error, last_error;
  454.    float plane1d, plane2d, E1, E2;
  455.    VECTOR plane1abc, plane2abc, temp; 
  456.    MATRIX H1, H2;
  457.    float dfdu, dfdv, dgdu, dgdv;
  458.  
  459.    /* get initial values for u and v from tree node */
  460.    u = node->u_mid;
  461.    v = node->v_mid;
  462.  
  463.    /* calculate plane1, plane2 */
  464.    cross_product(&plane1abc, r.org, r.dir);   
  465.    normalize(&plane1abc);
  466.    cross_product(&plane2abc, plane1abc, r.dir);
  467.    normalize(&plane2abc);
  468.  
  469.    plane1d = dot_product(plane1abc, r.org);
  470.    plane2d = dot_product(plane2abc, r.org); 
  471.  
  472.    /* calculate the H matrices */
  473.    for (i=0; i<4; i++) {
  474.       for (j=0; j<4; j++) {
  475.          temp.x = surface->geomx[i][j];
  476.          temp.y = surface->geomy[i][j];
  477.          temp.z = surface->geomz[i][j];
  478.          H1[i][j] = dot_product(plane1abc, temp);
  479.          H2[i][j] = dot_product(plane2abc, temp); 
  480.       }
  481.    }
  482.  
  483.    /* begin iteration */
  484.  
  485.    error = 0;
  486.    iterations = 0;
  487.    success = failure = FALSE;
  488.    while (!success && !failure) {
  489.       E1 = eval_patch(u, surface->M, H1, surface->MT, v) - plane1d; 
  490.       E2 = eval_patch(u, surface->M, H2, surface->MT, v) - plane2d; 
  491.       last_error = error;
  492.       error = fabs(E1) + fabs(E2);
  493.       if (error < NEWTON_TOL) 
  494.          success = TRUE;
  495.       else if ((iterations >= MAX_ITER) && (error >= last_error)) 
  496.          failure = TRUE;
  497.       else {                /****** calc newton step ***  f = E1, g = E2 **/ 
  498.  
  499.          dfdu = eval_patch(u, surface->Md, H1, surface->MT, v);
  500.          dfdv = eval_patch(u, surface->M, H1, surface->MdT, v);
  501.          dgdu = eval_patch(u, surface->Md, H2, surface->MT, v);
  502.          dgdv = eval_patch(u, surface->M, H2, surface->MdT, v);
  503.  
  504.          u += ((E2*dfdv-E1*dgdv)/(dfdu*dgdv-dfdv*dgdu));  /* Delta u */
  505.          v += ((E2*dfdu-E1*dgdu)/(dfdv*dgdu-dfdu*dgdv));  /* Delta v */
  506.  
  507.          iterations++; 
  508.       }
  509.    }
  510.   
  511.    if (success) {
  512.       irec->isect.x = eval_patch(u, surface->M, surface->geomx, surface->MT, v);
  513.       irec->isect.y = eval_patch(u, surface->M, surface->geomy, surface->MT, v);
  514.       irec->isect.z = eval_patch(u, surface->M, surface->geomz, surface->MT, v);
  515.       irec->t = length (r.org, irec->isect);
  516.       irec->u = u;
  517.       irec->v = v;
  518.    }
  519.    else irec->t = -1.0;
  520. }
  521.  
  522. /************************************************************************/
  523.  
  524. Isectrec IsectPatch (surface, r)
  525. PatchPtr surface;
  526. RAY r;
  527. {
  528.    int i;
  529.    float tmin, t;
  530.    ListPtr this, head; 
  531.    ListNode thisnode;
  532.    Isectrec isect, isectmin;
  533.  
  534.    head = insert_node(NULL, surface->tree, 0); 
  535.    tmin = INF;
  536.    isectmin.t = -1.0;              /* Use -1.0 to indicate no intersection */
  537.  
  538.    while (head) {
  539.       this = head;
  540.       head = this->next;
  541.       thisnode = *this;
  542.       free(this);
  543.       if (!thisnode.node->child[0]) {     /* if node is a leaf node */
  544.          newton(surface, thisnode.node, r, &isect); 
  545.          if ((isect.t != -1.0) && (isect.t < tmin)) {
  546.             tmin = isect.t;    
  547.             isectmin = isect;
  548.          }
  549.       } else {
  550.          for (i = 0; i < 4; i++) {                  /* for each child of node */
  551.             t = BoxIntersect(r, thisnode.node->child[i]->box_min,
  552.                                 thisnode.node->child[i]->box_max);
  553.             if ( t != -1.0)                         /* if ray intersects box */
  554.                head = insert_node(head, thisnode.node->child[i], t);
  555.          }
  556.       }
  557.       if ((head) && (tmin < head->t)) { free_list(head); head = NULL; }
  558.    }
  559.    return isectmin;
  560. }   
  561.       
  562. /********************************************************************/
  563.  
  564.  
  565. VECTOR NormalPatch (p, u, v)
  566. PatchPtr p;
  567. float u, v;
  568. {
  569.    VECTOR ddu, ddv, norm;
  570.  
  571.    ddu.x = eval_patch(u, p->Md, p->geomx, p->MT, v);
  572.    ddu.y = eval_patch(u, p->Md, p->geomy, p->MT, v);
  573.    ddu.z = eval_patch(u, p->Md, p->geomz, p->MT, v);
  574.  
  575.    ddv.x = eval_patch(u, p->M, p->geomx, p->MdT, v);
  576.    ddv.y = eval_patch(u, p->M, p->geomy, p->MdT, v);
  577.    ddv.z = eval_patch(u, p->M, p->geomz, p->MdT, v);
  578.  
  579.    cross_product(&norm, ddu, ddv);
  580.    normalize(&norm); 
  581.    return norm;
  582. }
  583.