Home | History | Annotate | Line # | Download | only in interface
      1 /*
      2 ** License Applicability. Except to the extent portions of this file are
      3 ** made subject to an alternative license as permitted in the SGI Free
      4 ** Software License B, Version 1.1 (the "License"), the contents of this
      5 ** file are subject only to the provisions of the License. You may not use
      6 ** this file except in compliance with the License. You may obtain a copy
      7 ** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
      8 ** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
      9 **
     10 ** http://oss.sgi.com/projects/FreeB
     11 **
     12 ** Note that, as provided in the License, the Software is distributed on an
     13 ** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
     14 ** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
     15 ** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
     16 ** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
     17 **
     18 ** Original Code. The Original Code is: OpenGL Sample Implementation,
     19 ** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
     20 ** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
     21 ** Copyright in any portions created by third parties is as indicated
     22 ** elsewhere herein. All Rights Reserved.
     23 **
     24 ** Additional Notice Provisions: The application programming interfaces
     25 ** established by SGI in conjunction with the Original Code are The
     26 ** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
     27 ** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
     28 ** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
     29 ** Window System(R) (Version 1.3), released October 19, 1998. This software
     30 ** was created using the OpenGL(R) version 1.2.1 Sample Implementation
     31 ** published by SGI, but has not been independently verified as being
     32 ** compliant with the OpenGL(R) version 1.2.1 Specification.
     33 **
     34 */
     35 /*
     36 */
     37 
     38 #include "gluos.h"
     39 #include <stdlib.h>
     40 #include <stdio.h>
     41 #include <GL/gl.h>
     42 #include <math.h>
     43 #include <assert.h>
     44 
     45 #include "glsurfeval.h"
     46 
     47 //extern int surfcount;
     48 
     49 //#define CRACK_TEST
     50 
     51 #define AVOID_ZERO_NORMAL
     52 
     53 #ifdef AVOID_ZERO_NORMAL
     54 #define myabs(x)  ((x>0)? x: (-x))
     55 #define MYZERO 0.000001
     56 #define MYDELTA 0.001
     57 #endif
     58 
     59 //#define USE_LOD
     60 #ifdef USE_LOD
     61 //#define LOD_EVAL_COORD(u,v) inDoEvalCoord2EM(u,v)
     62 #define LOD_EVAL_COORD(u,v) glEvalCoord2f(u,v)
     63 
     64 static void LOD_interpolate(REAL A[2], REAL B[2], REAL C[2], int j, int k, int pow2_level,
     65 			    REAL& u, REAL& v)
     66 {
     67   REAL a,a1,b,b1;
     68 
     69   a = ((REAL) j) / ((REAL) pow2_level);
     70   a1 = 1-a;
     71 
     72   if(j != 0)
     73     {
     74       b = ((REAL) k) / ((REAL)j);
     75       b1 = 1-b;
     76     }
     77   REAL x,y,z;
     78   x = a1;
     79   if(j==0)
     80     {
     81       y=0; z=0;
     82     }
     83   else{
     84     y = b1*a;
     85     z = b *a;
     86   }
     87 
     88   u = x*A[0] + y*B[0] + z*C[0];
     89   v = x*A[1] + y*B[1] + z*C[1];
     90 }
     91 
     92 void OpenGLSurfaceEvaluator::LOD_triangle(REAL A[2], REAL B[2], REAL C[2],
     93 			 int level)
     94 {
     95   int k,j;
     96   int pow2_level;
     97   /*compute 2^level*/
     98   pow2_level = 1;
     99 
    100   for(j=0; j<level; j++)
    101     pow2_level *= 2;
    102   for(j=0; j<=pow2_level-1; j++)
    103     {
    104       REAL u,v;
    105 
    106 /*      beginCallBack(GL_TRIANGLE_STRIP);*/
    107 glBegin(GL_TRIANGLE_STRIP);
    108       LOD_interpolate(A,B,C, j+1, j+1, pow2_level, u,v);
    109 #ifdef USE_LOD
    110       LOD_EVAL_COORD(u,v);
    111 //      glEvalCoord2f(u,v);
    112 #else
    113       inDoEvalCoord2EM(u,v);
    114 #endif
    115 
    116       for(k=0; k<=j; k++)
    117 	{
    118 	  LOD_interpolate(A,B,C,j,j-k,pow2_level, u,v);
    119 #ifdef USE_LOD
    120           LOD_EVAL_COORD(u,v);
    121 //	  glEvalCoord2f(u,v);
    122 #else
    123 	  inDoEvalCoord2EM(u,v);
    124 #endif
    125 
    126 	  LOD_interpolate(A,B,C,j+1,j-k,pow2_level, u,v);
    127 
    128 #ifdef USE_LOD
    129 	  LOD_EVAL_COORD(u,v);
    130 //	  glEvalCoord2f(u,v);
    131 #else
    132 	  inDoEvalCoord2EM(u,v);
    133 #endif
    134 	}
    135 //      endCallBack();
    136 glEnd();
    137     }
    138 }
    139 
    140 void OpenGLSurfaceEvaluator::LOD_eval(int num_vert, REAL* verts, int type,
    141 		     int level
    142 		     )
    143 {
    144   int i,k;
    145   switch(type){
    146   case GL_TRIANGLE_STRIP:
    147   case GL_QUAD_STRIP:
    148     for(i=2, k=4; i<=num_vert-2; i+=2, k+=4)
    149       {
    150 	LOD_triangle(verts+k-4, verts+k-2, verts+k,
    151 		     level
    152 		     );
    153 	LOD_triangle(verts+k-2, verts+k+2, verts+k,
    154 		     level
    155 		     );
    156       }
    157     if(num_vert % 2 ==1)
    158       {
    159 	LOD_triangle(verts+2*(num_vert-3), verts+2*(num_vert-2), verts+2*(num_vert-1),
    160 		     level
    161 		     );
    162       }
    163     break;
    164   case GL_TRIANGLE_FAN:
    165     for(i=1, k=2; i<=num_vert-2; i++, k+=2)
    166       {
    167 	LOD_triangle(verts,verts+k, verts+k+2,
    168 		     level
    169 		     );
    170       }
    171     break;
    172 
    173   default:
    174     fprintf(stderr, "typy not supported in LOD_\n");
    175   }
    176 }
    177 
    178 
    179 #endif //USE_LOD
    180 
    181 //#define  GENERIC_TEST
    182 #ifdef GENERIC_TEST
    183 extern float xmin, xmax, ymin, ymax, zmin, zmax; /*bounding box*/
    184 extern int temp_signal;
    185 
    186 static void gTessVertexSphere(float u, float v, float temp_normal[3], float temp_vertex[3])
    187 {
    188   float r=2.0;
    189   float Ox = 0.5*(xmin+xmax);
    190   float Oy = 0.5*(ymin+ymax);
    191   float Oz = 0.5*(zmin+zmax);
    192   float nx = cos(v) * sin(u);
    193   float ny = sin(v) * sin(u);
    194   float nz = cos(u);
    195   float x= Ox+r * nx;
    196   float y= Oy+r * ny;
    197   float z= Oz+r * nz;
    198 
    199   temp_normal[0] = nx;
    200   temp_normal[1] = ny;
    201   temp_normal[2] =  nz;
    202   temp_vertex[0] = x;
    203   temp_vertex[1] = y;
    204   temp_vertex[2] = z;
    205 
    206 //  glNormal3f(nx,ny,nz);
    207 //  glVertex3f(x,y,z);
    208 }
    209 
    210 static void gTessVertexCyl(float u, float v, float temp_normal[3], float temp_vertex[3])
    211 {
    212    float r=2.0;
    213   float Ox = 0.5*(xmin+xmax);
    214   float Oy = 0.5*(ymin+ymax);
    215   float Oz = 0.5*(zmin+zmax);
    216   float nx = cos(v);
    217   float ny = sin(v);
    218   float nz = 0;
    219   float x= Ox+r * nx;
    220   float y= Oy+r * ny;
    221   float z= Oz - 2*u;
    222 
    223   temp_normal[0] = nx;
    224   temp_normal[1] = ny;
    225   temp_normal[2] =  nz;
    226   temp_vertex[0] = x;
    227   temp_vertex[1] = y;
    228   temp_vertex[2] = z;
    229 
    230 /*
    231   glNormal3f(nx,ny,nz);
    232   glVertex3f(x,y,z);
    233 */
    234 }
    235 
    236 #endif //GENERIC_TEST
    237 
    238 void OpenGLSurfaceEvaluator::inBPMListEval(bezierPatchMesh* list)
    239 {
    240   bezierPatchMesh* temp;
    241   for(temp = list; temp != NULL; temp = temp->next)
    242     {
    243       inBPMEval(temp);
    244     }
    245 }
    246 
    247 void OpenGLSurfaceEvaluator::inBPMEval(bezierPatchMesh* bpm)
    248 {
    249   int i,j,k,l;
    250   float u,v;
    251 
    252   int ustride = bpm->bpatch->dimension * bpm->bpatch->vorder;
    253   int vstride = bpm->bpatch->dimension;
    254   inMap2f(
    255 	  (bpm->bpatch->dimension == 3)? GL_MAP2_VERTEX_3 : GL_MAP2_VERTEX_4,
    256 	  bpm->bpatch->umin,
    257 	  bpm->bpatch->umax,
    258 	  ustride,
    259 	  bpm->bpatch->uorder,
    260 	  bpm->bpatch->vmin,
    261 	  bpm->bpatch->vmax,
    262 	  vstride,
    263 	  bpm->bpatch->vorder,
    264 	  bpm->bpatch->ctlpoints);
    265 
    266   bpm->vertex_array = (float*) malloc(sizeof(float)* (bpm->index_UVarray/2) * 3+1); /*in case the origional dimenion is 4, then we need 4 space to pass to evaluator.*/
    267   assert(bpm->vertex_array);
    268   bpm->normal_array = (float*) malloc(sizeof(float)* (bpm->index_UVarray/2) * 3);
    269   assert(bpm->normal_array);
    270 #ifdef CRACK_TEST
    271 if(  global_ev_u1 ==2 &&   global_ev_u2 == 3
    272   && global_ev_v1 ==2 &&   global_ev_v2 == 3)
    273 {
    274 REAL vertex[4];
    275 REAL normal[4];
    276 #ifdef DEBUG
    277 printf("***number 1\n");
    278 #endif
    279 
    280 beginCallBack(GL_QUAD_STRIP, NULL);
    281 inEvalCoord2f(3.0, 3.0);
    282 inEvalCoord2f(2.0, 3.0);
    283 inEvalCoord2f(3.0, 2.7);
    284 inEvalCoord2f(2.0, 2.7);
    285 inEvalCoord2f(3.0, 2.0);
    286 inEvalCoord2f(2.0, 2.0);
    287 endCallBack(NULL);
    288 
    289 
    290 beginCallBack(GL_TRIANGLE_STRIP, NULL);
    291 inEvalCoord2f(2.0, 3.0);
    292 inEvalCoord2f(2.0, 2.0);
    293 inEvalCoord2f(2.0, 2.7);
    294 endCallBack(NULL);
    295 
    296 }
    297 
    298 /*
    299 if(  global_ev_u1 ==2 &&   global_ev_u2 == 3
    300   && global_ev_v1 ==1 &&   global_ev_v2 == 2)
    301 {
    302 #ifdef DEBUG
    303 printf("***number 2\n");
    304 #endif
    305 beginCallBack(GL_QUAD_STRIP);
    306 inEvalCoord2f(2.0, 2.0);
    307 inEvalCoord2f(2.0, 1.0);
    308 inEvalCoord2f(3.0, 2.0);
    309 inEvalCoord2f(3.0, 1.0);
    310 endCallBack();
    311 }
    312 */
    313 if(  global_ev_u1 ==1 &&   global_ev_u2 == 2
    314   && global_ev_v1 ==2 &&   global_ev_v2 == 3)
    315 {
    316 #ifdef DEBUG
    317 printf("***number 3\n");
    318 #endif
    319 beginCallBack(GL_QUAD_STRIP, NULL);
    320 inEvalCoord2f(2.0, 3.0);
    321 inEvalCoord2f(1.0, 3.0);
    322 inEvalCoord2f(2.0, 2.3);
    323 inEvalCoord2f(1.0, 2.3);
    324 inEvalCoord2f(2.0, 2.0);
    325 inEvalCoord2f(1.0, 2.0);
    326 endCallBack(NULL);
    327 
    328 beginCallBack(GL_TRIANGLE_STRIP, NULL);
    329 inEvalCoord2f(2.0, 2.3);
    330 inEvalCoord2f(2.0, 2.0);
    331 inEvalCoord2f(2.0, 3.0);
    332 endCallBack(NULL);
    333 
    334 }
    335 return;
    336 #endif
    337 
    338   k=0;
    339   l=0;
    340 
    341   for(i=0; i<bpm->index_length_array; i++)
    342     {
    343       beginCallBack(bpm->type_array[i], userData);
    344       for(j=0; j<bpm->length_array[i]; j++)
    345 	{
    346 	  u = bpm->UVarray[k];
    347 	  v = bpm->UVarray[k+1];
    348 	  inDoEvalCoord2NOGE(u,v,
    349 			     bpm->vertex_array+l,
    350 			     bpm->normal_array+l);
    351 
    352 	  normalCallBack(bpm->normal_array+l, userData);
    353 	  vertexCallBack(bpm->vertex_array+l, userData);
    354 
    355 	  k += 2;
    356 	  l += 3;
    357 	}
    358       endCallBack(userData);
    359     }
    360 }
    361 
    362 void OpenGLSurfaceEvaluator::inEvalPoint2(int i, int j)
    363 {
    364   REAL du, dv;
    365   REAL point[4];
    366   REAL normal[3];
    367   REAL u,v;
    368   du = (global_grid_u1 - global_grid_u0) / (REAL)global_grid_nu;
    369   dv = (global_grid_v1 - global_grid_v0) / (REAL)global_grid_nv;
    370   u = (i==global_grid_nu)? global_grid_u1:(global_grid_u0 + i*du);
    371   v = (j == global_grid_nv)? global_grid_v1: (global_grid_v0 +j*dv);
    372   inDoEvalCoord2(u,v,point,normal);
    373 }
    374 
    375 void OpenGLSurfaceEvaluator::inEvalCoord2f(REAL u, REAL v)
    376 {
    377 
    378   REAL point[4];
    379   REAL normal[3];
    380   inDoEvalCoord2(u,v,point, normal);
    381 }
    382 
    383 
    384 
    385 /*define a grid. store the values into the global variabls:
    386  * global_grid_*
    387  *These values will be used later by evaluating functions
    388  */
    389 void OpenGLSurfaceEvaluator::inMapGrid2f(int nu, REAL u0, REAL u1,
    390 		 int nv, REAL v0, REAL v1)
    391 {
    392  global_grid_u0 = u0;
    393  global_grid_u1 = u1;
    394  global_grid_nu = nu;
    395  global_grid_v0 = v0;
    396  global_grid_v1 = v1;
    397  global_grid_nv = nv;
    398 }
    399 
    400 void OpenGLSurfaceEvaluator::inEvalMesh2(int lowU, int lowV, int highU, int highV)
    401 {
    402   REAL du, dv;
    403   int i,j;
    404   REAL point[4];
    405   REAL normal[3];
    406   if(global_grid_nu == 0 || global_grid_nv == 0)
    407     return; /*no points need to be output*/
    408   du = (global_grid_u1 - global_grid_u0) / (REAL)global_grid_nu;
    409   dv = (global_grid_v1 - global_grid_v0) / (REAL)global_grid_nv;
    410 
    411   if(global_grid_nu >= global_grid_nv){
    412     for(i=lowU; i<highU; i++){
    413       REAL u1 = (i==global_grid_nu)? global_grid_u1:(global_grid_u0 + i*du);
    414       REAL u2 = ((i+1) == global_grid_nu)? global_grid_u1: (global_grid_u0+(i+1)*du);
    415 
    416       bgnqstrip();
    417       for(j=highV; j>=lowV; j--){
    418 	REAL v1 = (j == global_grid_nv)? global_grid_v1: (global_grid_v0 +j*dv);
    419 
    420 	inDoEvalCoord2(u1, v1, point, normal);
    421 	inDoEvalCoord2(u2, v1, point, normal);
    422       }
    423       endqstrip();
    424     }
    425   }
    426 
    427   else{
    428     for(i=lowV; i<highV; i++){
    429       REAL v1 = (i==global_grid_nv)? global_grid_v1:(global_grid_v0 + i*dv);
    430       REAL v2 = ((i+1) == global_grid_nv)? global_grid_v1: (global_grid_v0+(i+1)*dv);
    431 
    432       bgnqstrip();
    433       for(j=highU; j>=lowU; j--){
    434 	REAL u1 = (j == global_grid_nu)? global_grid_u1: (global_grid_u0 +j*du);
    435 	inDoEvalCoord2(u1, v2, point, normal);
    436 	inDoEvalCoord2(u1, v1, point, normal);
    437       }
    438       endqstrip();
    439     }
    440   }
    441 
    442 }
    443 
    444 void OpenGLSurfaceEvaluator::inMap2f(int k,
    445 	     REAL ulower,
    446 	     REAL uupper,
    447 	     int ustride,
    448 	     int uorder,
    449 	     REAL vlower,
    450 	     REAL vupper,
    451 	     int vstride,
    452 	     int vorder,
    453 	     REAL *ctlPoints)
    454 {
    455   int i,j,x;
    456   REAL *data = global_ev_ctlPoints;
    457 
    458 
    459 
    460   if(k == GL_MAP2_VERTEX_3) k=3;
    461   else if (k==GL_MAP2_VERTEX_4) k =4;
    462   else {
    463     printf("error in inMap2f, maptype=%i is wrong, k,map is not updated\n", k);
    464     return;
    465   }
    466 
    467   global_ev_k = k;
    468   global_ev_u1 = ulower;
    469   global_ev_u2 = uupper;
    470   global_ev_ustride = ustride;
    471   global_ev_uorder = uorder;
    472   global_ev_v1 = vlower;
    473   global_ev_v2 = vupper;
    474   global_ev_vstride = vstride;
    475   global_ev_vorder = vorder;
    476 
    477   /*copy the contrl points from ctlPoints to global_ev_ctlPoints*/
    478   for (i=0; i<uorder; i++) {
    479     for (j=0; j<vorder; j++) {
    480       for (x=0; x<k; x++) {
    481 	data[x] = ctlPoints[x];
    482       }
    483       ctlPoints += vstride;
    484       data += k;
    485     }
    486     ctlPoints += ustride - vstride * vorder;
    487   }
    488 
    489 }
    490 
    491 
    492 /*
    493  *given a point p with homegeneous coordiante (x,y,z,w),
    494  *let pu(x,y,z,w) be its partial derivative vector with
    495  *respect to u
    496  *and pv(x,y,z,w) be its partial derivative vector with repect to v.
    497  *This function returns the partial derivative vectors of the
    498  *inhomegensous coordinates, i.e.,
    499  * (x/w, y/w, z/w) with respect to u and v.
    500  */
    501 void OpenGLSurfaceEvaluator::inComputeFirstPartials(REAL *p, REAL *pu, REAL *pv)
    502 {
    503     pu[0] = pu[0]*p[3] - pu[3]*p[0];
    504     pu[1] = pu[1]*p[3] - pu[3]*p[1];
    505     pu[2] = pu[2]*p[3] - pu[3]*p[2];
    506 
    507     pv[0] = pv[0]*p[3] - pv[3]*p[0];
    508     pv[1] = pv[1]*p[3] - pv[3]*p[1];
    509     pv[2] = pv[2]*p[3] - pv[3]*p[2];
    510 }
    511 
    512 /*compute the cross product of pu and pv and normalize.
    513  *the normal is returned in retNormal
    514  * pu: dimension 3
    515  * pv: dimension 3
    516  * n: return normal, of dimension 3
    517  */
    518 void OpenGLSurfaceEvaluator::inComputeNormal2(REAL *pu, REAL *pv, REAL *n)
    519 {
    520   REAL mag;
    521 
    522   n[0] = pu[1]*pv[2] - pu[2]*pv[1];
    523   n[1] = pu[2]*pv[0] - pu[0]*pv[2];
    524   n[2] = pu[0]*pv[1] - pu[1]*pv[0];
    525 
    526   mag = sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
    527 
    528   if (mag > 0.0) {
    529      n[0] /= mag;
    530      n[1] /= mag;
    531      n[2] /= mag;
    532   }
    533 }
    534 
    535 
    536 
    537 /*Compute point and normal
    538  *see the head of inDoDomain2WithDerivs
    539  *for the meaning of the arguments
    540  */
    541 void OpenGLSurfaceEvaluator::inDoEvalCoord2(REAL u, REAL v,
    542 			   REAL *retPoint, REAL *retNormal)
    543 {
    544 
    545   REAL du[4];
    546   REAL dv[4];
    547 
    548 
    549   assert(global_ev_k>=3 && global_ev_k <= 4);
    550   /*compute homegeneous point and partial derivatives*/
    551   inDoDomain2WithDerivs(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);
    552 
    553 #ifdef AVOID_ZERO_NORMAL
    554 
    555   if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
    556     {
    557 
    558       REAL tempdu[4];
    559       REAL tempdata[4];
    560       REAL u1 = global_ev_u1;
    561       REAL u2 = global_ev_u2;
    562       if(u-MYDELTA*(u2-u1) < u1)
    563 	u = u+ MYDELTA*(u2-u1);
    564       else
    565 	u = u-MYDELTA*(u2-u1);
    566       inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
    567     }
    568   if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
    569     {
    570       REAL tempdv[4];
    571       REAL tempdata[4];
    572       REAL v1 = global_ev_v1;
    573       REAL v2 = global_ev_v2;
    574       if(v-MYDELTA*(v2-v1) < v1)
    575 	v = v+ MYDELTA*(v2-v1);
    576       else
    577 	v = v-MYDELTA*(v2-v1);
    578       inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
    579     }
    580 #endif
    581 
    582 
    583   /*compute normal*/
    584   switch(global_ev_k){
    585   case 3:
    586     inComputeNormal2(du, dv, retNormal);
    587 
    588     break;
    589   case 4:
    590     inComputeFirstPartials(retPoint, du, dv);
    591     inComputeNormal2(du, dv, retNormal);
    592     /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
    593     retPoint[0] /= retPoint[3];
    594     retPoint[1] /= retPoint[3];
    595     retPoint[2] /= retPoint[3];
    596     break;
    597   }
    598   /*output this vertex*/
    599 /*  inMeshStreamInsert(global_ms, retPoint, retNormal);*/
    600 
    601 
    602 
    603   glNormal3fv(retNormal);
    604   glVertex3fv(retPoint);
    605 
    606 
    607 
    608 
    609   #ifdef DEBUG
    610   printf("vertex(%f,%f,%f)\n", retPoint[0],retPoint[1],retPoint[2]);
    611   #endif
    612 
    613 
    614 
    615 }
    616 
    617 /*Compute point and normal
    618  *see the head of inDoDomain2WithDerivs
    619  *for the meaning of the arguments
    620  */
    621 void OpenGLSurfaceEvaluator::inDoEvalCoord2NOGE_BU(REAL u, REAL v,
    622 			   REAL *retPoint, REAL *retNormal)
    623 {
    624 
    625   REAL du[4];
    626   REAL dv[4];
    627 
    628 
    629   assert(global_ev_k>=3 && global_ev_k <= 4);
    630   /*compute homegeneous point and partial derivatives*/
    631 //   inPreEvaluateBU(global_ev_k, global_ev_uorder, global_ev_vorder, (u-global_ev_u1)/(global_ev_u2-global_ev_u1), global_ev_ctlPoints);
    632   inDoDomain2WithDerivsBU(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);
    633 
    634 
    635 #ifdef AVOID_ZERO_NORMAL
    636 
    637   if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
    638     {
    639 
    640       REAL tempdu[4];
    641       REAL tempdata[4];
    642       REAL u1 = global_ev_u1;
    643       REAL u2 = global_ev_u2;
    644       if(u-MYDELTA*(u2-u1) < u1)
    645 	u = u+ MYDELTA*(u2-u1);
    646       else
    647 	u = u-MYDELTA*(u2-u1);
    648       inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
    649     }
    650   if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
    651     {
    652       REAL tempdv[4];
    653       REAL tempdata[4];
    654       REAL v1 = global_ev_v1;
    655       REAL v2 = global_ev_v2;
    656       if(v-MYDELTA*(v2-v1) < v1)
    657 	v = v+ MYDELTA*(v2-v1);
    658       else
    659 	v = v-MYDELTA*(v2-v1);
    660       inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
    661     }
    662 #endif
    663 
    664   /*compute normal*/
    665   switch(global_ev_k){
    666   case 3:
    667     inComputeNormal2(du, dv, retNormal);
    668     break;
    669   case 4:
    670     inComputeFirstPartials(retPoint, du, dv);
    671     inComputeNormal2(du, dv, retNormal);
    672     /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
    673     retPoint[0] /= retPoint[3];
    674     retPoint[1] /= retPoint[3];
    675     retPoint[2] /= retPoint[3];
    676     break;
    677   }
    678 }
    679 
    680 /*Compute point and normal
    681  *see the head of inDoDomain2WithDerivs
    682  *for the meaning of the arguments
    683  */
    684 void OpenGLSurfaceEvaluator::inDoEvalCoord2NOGE_BV(REAL u, REAL v,
    685 			   REAL *retPoint, REAL *retNormal)
    686 {
    687 
    688   REAL du[4];
    689   REAL dv[4];
    690 
    691 
    692   assert(global_ev_k>=3 && global_ev_k <= 4);
    693   /*compute homegeneous point and partial derivatives*/
    694 //   inPreEvaluateBV(global_ev_k, global_ev_uorder, global_ev_vorder, (v-global_ev_v1)/(global_ev_v2-global_ev_v1), global_ev_ctlPoints);
    695 
    696   inDoDomain2WithDerivsBV(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);
    697 
    698 
    699 #ifdef AVOID_ZERO_NORMAL
    700 
    701   if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
    702     {
    703 
    704       REAL tempdu[4];
    705       REAL tempdata[4];
    706       REAL u1 = global_ev_u1;
    707       REAL u2 = global_ev_u2;
    708       if(u-MYDELTA*(u2-u1) < u1)
    709 	u = u+ MYDELTA*(u2-u1);
    710       else
    711 	u = u-MYDELTA*(u2-u1);
    712       inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
    713     }
    714   if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
    715     {
    716       REAL tempdv[4];
    717       REAL tempdata[4];
    718       REAL v1 = global_ev_v1;
    719       REAL v2 = global_ev_v2;
    720       if(v-MYDELTA*(v2-v1) < v1)
    721 	v = v+ MYDELTA*(v2-v1);
    722       else
    723 	v = v-MYDELTA*(v2-v1);
    724       inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
    725     }
    726 #endif
    727 
    728   /*compute normal*/
    729   switch(global_ev_k){
    730   case 3:
    731     inComputeNormal2(du, dv, retNormal);
    732     break;
    733   case 4:
    734     inComputeFirstPartials(retPoint, du, dv);
    735     inComputeNormal2(du, dv, retNormal);
    736     /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
    737     retPoint[0] /= retPoint[3];
    738     retPoint[1] /= retPoint[3];
    739     retPoint[2] /= retPoint[3];
    740     break;
    741   }
    742 }
    743 
    744 
    745 /*Compute point and normal
    746  *see the head of inDoDomain2WithDerivs
    747  *for the meaning of the arguments
    748  */
    749 void OpenGLSurfaceEvaluator::inDoEvalCoord2NOGE(REAL u, REAL v,
    750 			   REAL *retPoint, REAL *retNormal)
    751 {
    752 
    753   REAL du[4];
    754   REAL dv[4];
    755 
    756 
    757   assert(global_ev_k>=3 && global_ev_k <= 4);
    758   /*compute homegeneous point and partial derivatives*/
    759   inDoDomain2WithDerivs(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);
    760 
    761 
    762 #ifdef AVOID_ZERO_NORMAL
    763 
    764   if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
    765     {
    766 
    767       REAL tempdu[4];
    768       REAL tempdata[4];
    769       REAL u1 = global_ev_u1;
    770       REAL u2 = global_ev_u2;
    771       if(u-MYDELTA*(u2-u1) < u1)
    772 	u = u+ MYDELTA*(u2-u1);
    773       else
    774 	u = u-MYDELTA*(u2-u1);
    775       inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
    776     }
    777   if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
    778     {
    779       REAL tempdv[4];
    780       REAL tempdata[4];
    781       REAL v1 = global_ev_v1;
    782       REAL v2 = global_ev_v2;
    783       if(v-MYDELTA*(v2-v1) < v1)
    784 	v = v+ MYDELTA*(v2-v1);
    785       else
    786 	v = v-MYDELTA*(v2-v1);
    787       inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
    788     }
    789 #endif
    790 
    791   /*compute normal*/
    792   switch(global_ev_k){
    793   case 3:
    794     inComputeNormal2(du, dv, retNormal);
    795     break;
    796   case 4:
    797     inComputeFirstPartials(retPoint, du, dv);
    798     inComputeNormal2(du, dv, retNormal);
    799     /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
    800     retPoint[0] /= retPoint[3];
    801     retPoint[1] /= retPoint[3];
    802     retPoint[2] /= retPoint[3];
    803     break;
    804   }
    805 //  glNormal3fv(retNormal);
    806 //  glVertex3fv(retPoint);
    807 }
    808 
    809 void OpenGLSurfaceEvaluator::inPreEvaluateBV(int k, int uorder, int vorder, REAL vprime, REAL *baseData)
    810 {
    811   int j,row,col;
    812   REAL p, pdv;
    813   REAL *data;
    814 
    815   if(global_vprime != vprime || global_vorder != vorder) {
    816     inPreEvaluateWithDeriv(vorder, vprime, global_vcoeff, global_vcoeffDeriv);
    817     global_vprime = vprime;
    818     global_vorder = vorder;
    819   }
    820 
    821   for(j=0; j<k; j++){
    822     data = baseData+j;
    823     for(row=0; row<uorder; row++){
    824       p = global_vcoeff[0] * (*data);
    825       pdv = global_vcoeffDeriv[0] * (*data);
    826       data += k;
    827       for(col = 1; col < vorder; col++){
    828 	p += global_vcoeff[col] *  (*data);
    829 	pdv += global_vcoeffDeriv[col] * (*data);
    830 	data += k;
    831       }
    832       global_BV[row][j]  = p;
    833       global_PBV[row][j]  = pdv;
    834     }
    835   }
    836 }
    837 
    838 void OpenGLSurfaceEvaluator::inPreEvaluateBU(int k, int uorder, int vorder, REAL uprime, REAL *baseData)
    839 {
    840   int j,row,col;
    841   REAL p, pdu;
    842   REAL *data;
    843 
    844   if(global_uprime != uprime || global_uorder != uorder) {
    845     inPreEvaluateWithDeriv(uorder, uprime, global_ucoeff, global_ucoeffDeriv);
    846     global_uprime = uprime;
    847     global_uorder = uorder;
    848   }
    849 
    850   for(j=0; j<k; j++){
    851     data = baseData+j;
    852     for(col=0; col<vorder; col++){
    853       data = baseData+j + k*col;
    854       p = global_ucoeff[0] * (*data);
    855       pdu = global_ucoeffDeriv[0] * (*data);
    856       data += k*uorder;
    857       for(row = 1; row < uorder; row++){
    858 	p += global_ucoeff[row] *  (*data);
    859 	pdu += global_ucoeffDeriv[row] * (*data);
    860 	data += k * uorder;
    861       }
    862       global_BU[col][j]  = p;
    863       global_PBU[col][j]  = pdu;
    864     }
    865   }
    866 }
    867 
    868 void OpenGLSurfaceEvaluator::inDoDomain2WithDerivsBU(int k, REAL u, REAL v,
    869 						      REAL u1, REAL u2, int uorder,
    870 						      REAL v1, REAL v2, int vorder,
    871 						      REAL *baseData,
    872 						      REAL *retPoint, REAL* retdu, REAL *retdv)
    873 {
    874   int j, col;
    875 
    876   REAL vprime;
    877 
    878 
    879   if((u2 == u1) || (v2 == v1))
    880     return;
    881 
    882   vprime = (v - v1) / (v2 - v1);
    883 
    884 
    885   if(global_vprime != vprime || global_vorder != vorder) {
    886     inPreEvaluateWithDeriv(vorder, vprime, global_vcoeff, global_vcoeffDeriv);
    887     global_vprime = vprime;
    888     global_vorder = vorder;
    889   }
    890 
    891 
    892   for(j=0; j<k; j++)
    893     {
    894       retPoint[j] = retdu[j] = retdv[j] = 0.0;
    895       for (col = 0; col < vorder; col++)  {
    896 	retPoint[j] += global_BU[col][j] * global_vcoeff[col];
    897 	retdu[j] += global_PBU[col][j] * global_vcoeff[col];
    898 	retdv[j] += global_BU[col][j] * global_vcoeffDeriv[col];
    899       }
    900     }
    901 }
    902 
    903 void OpenGLSurfaceEvaluator::inDoDomain2WithDerivsBV(int k, REAL u, REAL v,
    904 						      REAL u1, REAL u2, int uorder,
    905 						      REAL v1, REAL v2, int vorder,
    906 						      REAL *baseData,
    907 						      REAL *retPoint, REAL* retdu, REAL *retdv)
    908 {
    909   int j, row;
    910   REAL uprime;
    911 
    912 
    913   if((u2 == u1) || (v2 == v1))
    914     return;
    915   uprime = (u - u1) / (u2 - u1);
    916 
    917 
    918   if(global_uprime != uprime || global_uorder != uorder) {
    919     inPreEvaluateWithDeriv(uorder, uprime, global_ucoeff, global_ucoeffDeriv);
    920     global_uprime = uprime;
    921     global_uorder = uorder;
    922   }
    923 
    924 
    925   for(j=0; j<k; j++)
    926     {
    927       retPoint[j] = retdu[j] = retdv[j] = 0.0;
    928       for (row = 0; row < uorder; row++)  {
    929 	retPoint[j] += global_BV[row][j] * global_ucoeff[row];
    930 	retdu[j] += global_BV[row][j] * global_ucoeffDeriv[row];
    931 	retdv[j] += global_PBV[row][j] * global_ucoeff[row];
    932       }
    933     }
    934 }
    935 
    936 
    937 /*
    938  *given a Bezier surface, and parameter (u,v), compute the point in the object space,
    939  *and the normal
    940  *k: the dimension of the object space: usually 2,3,or 4.
    941  *u,v: the paramter pair.
    942  *u1,u2,uorder: the Bezier polynomial of u coord is defined on [u1,u2] with order uorder.
    943  *v1,v2,vorder: the Bezier polynomial of v coord is defined on [v1,v2] with order vorder.
    944  *baseData: contrl points. arranged as: (u,v,k).
    945  *retPoint:  the computed point (one point) with dimension k.
    946  *retdu: the computed partial derivative with respect to u.
    947  *retdv: the computed partial derivative with respect to v.
    948  */
    949 void OpenGLSurfaceEvaluator::inDoDomain2WithDerivs(int k, REAL u, REAL v,
    950 				REAL u1, REAL u2, int uorder,
    951 				REAL v1,  REAL v2, int vorder,
    952 				REAL *baseData,
    953 				REAL *retPoint, REAL *retdu, REAL *retdv)
    954 {
    955     int j, row, col;
    956     REAL uprime;
    957     REAL vprime;
    958     REAL p;
    959     REAL pdv;
    960     REAL *data;
    961 
    962     if((u2 == u1) || (v2 == v1))
    963 	return;
    964     uprime = (u - u1) / (u2 - u1);
    965     vprime = (v - v1) / (v2 - v1);
    966 
    967     /* Compute coefficients for values and derivs */
    968 
    969     /* Use already cached values if possible */
    970     if(global_uprime != uprime || global_uorder != uorder) {
    971         inPreEvaluateWithDeriv(uorder, uprime, global_ucoeff, global_ucoeffDeriv);
    972 	global_uorder = uorder;
    973 	global_uprime = uprime;
    974     }
    975     if (global_vprime != vprime ||
    976 	  global_vorder != vorder) {
    977 	inPreEvaluateWithDeriv(vorder, vprime, global_vcoeff, global_vcoeffDeriv);
    978 	global_vorder = vorder;
    979 	global_vprime = vprime;
    980     }
    981 
    982     for (j = 0; j < k; j++) {
    983 	data=baseData+j;
    984 	retPoint[j] = retdu[j] = retdv[j] = 0.0;
    985 	for (row = 0; row < uorder; row++)  {
    986 	    /*
    987 	    ** Minor optimization.
    988 	    ** The col == 0 part of the loop is extracted so we don't
    989 	    ** have to initialize p and pdv to 0.
    990 	    */
    991 	    p = global_vcoeff[0] * (*data);
    992 	    pdv = global_vcoeffDeriv[0] * (*data);
    993 	    data += k;
    994 	    for (col = 1; col < vorder; col++) {
    995 		/* Incrementally build up p, pdv value */
    996 		p += global_vcoeff[col] * (*data);
    997 		pdv += global_vcoeffDeriv[col] * (*data);
    998 		data += k;
    999 	    }
   1000 	    /* Use p, pdv value to incrementally add up r, du, dv */
   1001 	    retPoint[j] += global_ucoeff[row] * p;
   1002 	    retdu[j] += global_ucoeffDeriv[row] * p;
   1003 	    retdv[j] += global_ucoeff[row] * pdv;
   1004 	}
   1005     }
   1006 }
   1007 
   1008 
   1009 /*
   1010  *compute the Bezier polynomials C[n,j](v) for all j at v with
   1011  *return values stored in coeff[], where
   1012  *  C[n,j](v) = (n,j) * v^j * (1-v)^(n-j),
   1013  *  j=0,1,2,...,n.
   1014  *order : n+1
   1015  *vprime: v
   1016  *coeff : coeff[j]=C[n,j](v), this array store the returned values.
   1017  *The algorithm is a recursive scheme:
   1018  *   C[0,0]=1;
   1019  *   C[n,j](v) = (1-v)*C[n-1,j](v) + v*C[n-1,j-1](v), n>=1
   1020  *This code is copied from opengl/soft/so_eval.c:PreEvaluate
   1021  */
   1022 void OpenGLSurfaceEvaluator::inPreEvaluate(int order, REAL vprime, REAL *coeff)
   1023 {
   1024   int i, j;
   1025   REAL oldval, temp;
   1026   REAL oneMinusvprime;
   1027 
   1028   /*
   1029    * Minor optimization
   1030    * Compute orders 1 and 2 outright, and set coeff[0], coeff[1] to
   1031      * their i==1 loop values to avoid the initialization and the i==1 loop.
   1032      */
   1033   if (order == 1) {
   1034     coeff[0] = 1.0;
   1035     return;
   1036   }
   1037 
   1038   oneMinusvprime = 1-vprime;
   1039   coeff[0] = oneMinusvprime;
   1040   coeff[1] = vprime;
   1041   if (order == 2) return;
   1042 
   1043   for (i = 2; i < order; i++) {
   1044     oldval = coeff[0] * vprime;
   1045     coeff[0] = oneMinusvprime * coeff[0];
   1046     for (j = 1; j < i; j++) {
   1047       temp = oldval;
   1048       oldval = coeff[j] * vprime;
   1049 	    coeff[j] = temp + oneMinusvprime * coeff[j];
   1050     }
   1051     coeff[j] = oldval;
   1052   }
   1053 }
   1054 
   1055 /*
   1056  *compute the Bezier polynomials C[n,j](v) and derivatives for all j at v with
   1057  *return values stored in coeff[] and coeffDeriv[].
   1058  *see the head of function inPreEvaluate for the definition of C[n,j](v)
   1059  *and how to compute the values.
   1060  *The algorithm to compute the derivative is:
   1061  *   dC[0,0](v) = 0.
   1062  *   dC[n,j](v) = n*(dC[n-1,j-1](v) - dC[n-1,j](v)).
   1063  *
   1064  *This code is copied from opengl/soft/so_eval.c:PreEvaluateWidthDeriv
   1065  */
   1066 void OpenGLSurfaceEvaluator::inPreEvaluateWithDeriv(int order, REAL vprime,
   1067     REAL *coeff, REAL *coeffDeriv)
   1068 {
   1069   int i, j;
   1070   REAL oldval, temp;
   1071   REAL oneMinusvprime;
   1072 
   1073   oneMinusvprime = 1-vprime;
   1074   /*
   1075    * Minor optimization
   1076    * Compute orders 1 and 2 outright, and set coeff[0], coeff[1] to
   1077    * their i==1 loop values to avoid the initialization and the i==1 loop.
   1078    */
   1079   if (order == 1) {
   1080     coeff[0] = 1.0;
   1081     coeffDeriv[0] = 0.0;
   1082     return;
   1083   } else if (order == 2) {
   1084     coeffDeriv[0] = -1.0;
   1085     coeffDeriv[1] = 1.0;
   1086     coeff[0] = oneMinusvprime;
   1087     coeff[1] = vprime;
   1088     return;
   1089   }
   1090   coeff[0] = oneMinusvprime;
   1091   coeff[1] = vprime;
   1092   for (i = 2; i < order - 1; i++) {
   1093     oldval = coeff[0] * vprime;
   1094     coeff[0] = oneMinusvprime * coeff[0];
   1095     for (j = 1; j < i; j++) {
   1096       temp = oldval;
   1097       oldval = coeff[j] * vprime;
   1098       coeff[j] = temp + oneMinusvprime * coeff[j];
   1099     }
   1100     coeff[j] = oldval;
   1101   }
   1102   coeffDeriv[0] = -coeff[0];
   1103   /*
   1104    ** Minor optimization:
   1105    ** Would make this a "for (j=1; j<order-1; j++)" loop, but it is always
   1106    ** executed at least once, so this is more efficient.
   1107    */
   1108   j=1;
   1109   do {
   1110     coeffDeriv[j] = coeff[j-1] - coeff[j];
   1111     j++;
   1112   } while (j < order - 1);
   1113   coeffDeriv[j] = coeff[j-1];
   1114 
   1115   oldval = coeff[0] * vprime;
   1116   coeff[0] = oneMinusvprime * coeff[0];
   1117   for (j = 1; j < i; j++) {
   1118     temp = oldval;
   1119     oldval = coeff[j] * vprime;
   1120     coeff[j] = temp + oneMinusvprime * coeff[j];
   1121   }
   1122   coeff[j] = oldval;
   1123 }
   1124 
   1125 void OpenGLSurfaceEvaluator::inEvalULine(int n_points, REAL v, REAL* u_vals,
   1126 	int stride, REAL ret_points[][3], REAL ret_normals[][3])
   1127 {
   1128   int i,k;
   1129   REAL temp[4];
   1130 inPreEvaluateBV_intfac(v);
   1131 
   1132   for(i=0,k=0; i<n_points; i++, k += stride)
   1133     {
   1134       inDoEvalCoord2NOGE_BV(u_vals[k],v,temp, ret_normals[i]);
   1135 
   1136       ret_points[i][0] = temp[0];
   1137       ret_points[i][1] = temp[1];
   1138       ret_points[i][2] = temp[2];
   1139 
   1140     }
   1141 
   1142 }
   1143 
   1144 void OpenGLSurfaceEvaluator::inEvalVLine(int n_points, REAL u, REAL* v_vals,
   1145 	int stride, REAL ret_points[][3], REAL ret_normals[][3])
   1146 {
   1147   int i,k;
   1148   REAL temp[4];
   1149 inPreEvaluateBU_intfac(u);
   1150   for(i=0,k=0; i<n_points; i++, k += stride)
   1151     {
   1152       inDoEvalCoord2NOGE_BU(u, v_vals[k], temp, ret_normals[i]);
   1153       ret_points[i][0] = temp[0];
   1154       ret_points[i][1] = temp[1];
   1155       ret_points[i][2] = temp[2];
   1156     }
   1157 }
   1158 
   1159 
   1160 /*triangulate a strip bounded by two lines which are parallel  to U-axis
   1161  *upperVerts: the verteces on the upper line
   1162  *lowerVertx: the verteces on the lower line
   1163  *n_upper >=1
   1164  *n_lower >=1
   1165  */
   1166 void OpenGLSurfaceEvaluator::inEvalUStrip(int n_upper, REAL v_upper, REAL* upper_val, int n_lower, REAL v_lower, REAL* lower_val)
   1167 {
   1168   int i,j,k,l;
   1169   REAL leftMostV[2];
   1170  typedef REAL REAL3[3];
   1171 
   1172   REAL3* upperXYZ = (REAL3*) malloc(sizeof(REAL3)*n_upper);
   1173   assert(upperXYZ);
   1174   REAL3* upperNormal = (REAL3*) malloc(sizeof(REAL3) * n_upper);
   1175   assert(upperNormal);
   1176   REAL3* lowerXYZ = (REAL3*) malloc(sizeof(REAL3)*n_lower);
   1177   assert(lowerXYZ);
   1178   REAL3* lowerNormal = (REAL3*) malloc(sizeof(REAL3) * n_lower);
   1179   assert(lowerNormal);
   1180 
   1181   inEvalULine(n_upper, v_upper, upper_val,  1, upperXYZ, upperNormal);
   1182   inEvalULine(n_lower, v_lower, lower_val,  1, lowerXYZ, lowerNormal);
   1183 
   1184 
   1185 
   1186   REAL* leftMostXYZ;
   1187   REAL* leftMostNormal;
   1188 
   1189   /*
   1190    *the algorithm works by scanning from left to right.
   1191    *leftMostV: the left most of the remaining verteces (on both upper and lower).
   1192    *           it could an element of upperVerts or lowerVerts.
   1193    *i: upperVerts[i] is the first vertex to the right of leftMostV on upper line   *j: lowerVerts[j] is the first vertex to the right of leftMostV on lower line   */
   1194 
   1195   /*initialize i,j,and leftMostV
   1196    */
   1197   if(upper_val[0] <= lower_val[0])
   1198     {
   1199       i=1;
   1200       j=0;
   1201 
   1202       leftMostV[0] = upper_val[0];
   1203       leftMostV[1] = v_upper;
   1204       leftMostXYZ = upperXYZ[0];
   1205       leftMostNormal = upperNormal[0];
   1206     }
   1207   else
   1208     {
   1209       i=0;
   1210       j=1;
   1211 
   1212       leftMostV[0] = lower_val[0];
   1213       leftMostV[1] = v_lower;
   1214 
   1215       leftMostXYZ = lowerXYZ[0];
   1216       leftMostNormal = lowerNormal[0];
   1217     }
   1218 
   1219   /*the main loop.
   1220    *the invariance is that:
   1221    *at the beginning of each loop, the meaning of i,j,and leftMostV are
   1222    *maintained
   1223    */
   1224   while(1)
   1225     {
   1226       if(i >= n_upper) /*case1: no more in upper*/
   1227         {
   1228           if(j<n_lower-1) /*at least two vertices in lower*/
   1229             {
   1230               bgntfan();
   1231 	      glNormal3fv(leftMostNormal);
   1232               glVertex3fv(leftMostXYZ);
   1233 
   1234               while(j<n_lower){
   1235 		glNormal3fv(lowerNormal[j]);
   1236 		glVertex3fv(lowerXYZ[j]);
   1237 		j++;
   1238 
   1239               }
   1240               endtfan();
   1241             }
   1242           break; /*exit the main loop*/
   1243         }
   1244       else if(j>= n_lower) /*case2: no more in lower*/
   1245         {
   1246           if(i<n_upper-1) /*at least two vertices in upper*/
   1247             {
   1248               bgntfan();
   1249 	      glNormal3fv(leftMostNormal);
   1250 	      glVertex3fv(leftMostXYZ);
   1251 
   1252               for(k=n_upper-1; k>=i; k--) /*reverse order for two-side lighting*/
   1253 		{
   1254 		  glNormal3fv(upperNormal[k]);
   1255 		  glVertex3fv(upperXYZ[k]);
   1256 		}
   1257 
   1258               endtfan();
   1259             }
   1260           break; /*exit the main loop*/
   1261         }
   1262       else /* case3: neither is empty, plus the leftMostV, there is at least one triangle to output*/
   1263         {
   1264           if(upper_val[i] <= lower_val[j])
   1265             {
   1266 	      bgntfan();
   1267 
   1268 	      glNormal3fv(lowerNormal[j]);
   1269 	      glVertex3fv(lowerXYZ[j]);
   1270 
   1271               /*find the last k>=i such that
   1272                *upperverts[k][0] <= lowerverts[j][0]
   1273                */
   1274               k=i;
   1275 
   1276               while(k<n_upper)
   1277                 {
   1278                   if(upper_val[k] > lower_val[j])
   1279                     break;
   1280                   k++;
   1281 
   1282                 }
   1283               k--;
   1284 
   1285 
   1286               for(l=k; l>=i; l--)/*the reverse is for two-side lighting*/
   1287                 {
   1288 		  glNormal3fv(upperNormal[l]);
   1289 		  glVertex3fv(upperXYZ[l]);
   1290 
   1291                 }
   1292 	      glNormal3fv(leftMostNormal);
   1293 	      glVertex3fv(leftMostXYZ);
   1294 
   1295               endtfan();
   1296 
   1297               /*update i and leftMostV for next loop
   1298                */
   1299               i = k+1;
   1300 
   1301 	      leftMostV[0] = upper_val[k];
   1302 	      leftMostV[1] = v_upper;
   1303 	      leftMostNormal = upperNormal[k];
   1304 	      leftMostXYZ = upperXYZ[k];
   1305             }
   1306           else /*upperVerts[i][0] > lowerVerts[j][0]*/
   1307             {
   1308 	      bgntfan();
   1309 	      glNormal3fv(upperNormal[i]);
   1310 	      glVertex3fv(upperXYZ[i]);
   1311 
   1312               glNormal3fv(leftMostNormal);
   1313 	      glVertex3fv(leftMostXYZ);
   1314 
   1315 
   1316               /*find the last k>=j such that
   1317                *lowerverts[k][0] < upperverts[i][0]
   1318                */
   1319               k=j;
   1320               while(k< n_lower)
   1321                 {
   1322                   if(lower_val[k] >= upper_val[i])
   1323                     break;
   1324 		  glNormal3fv(lowerNormal[k]);
   1325 		  glVertex3fv(lowerXYZ[k]);
   1326 
   1327                   k++;
   1328                 }
   1329               endtfan();
   1330 
   1331               /*update j and leftMostV for next loop
   1332                */
   1333               j=k;
   1334 	      leftMostV[0] = lower_val[j-1];
   1335 	      leftMostV[1] = v_lower;
   1336 
   1337 	      leftMostNormal = lowerNormal[j-1];
   1338 	      leftMostXYZ = lowerXYZ[j-1];
   1339             }
   1340         }
   1341     }
   1342   //clean up
   1343   free(upperXYZ);
   1344   free(lowerXYZ);
   1345   free(upperNormal);
   1346   free(lowerNormal);
   1347 }
   1348 
   1349 /*triangulate a strip bounded by two lines which are parallel  to V-axis
   1350  *leftVerts: the verteces on the left line
   1351  *rightVertx: the verteces on the right line
   1352  *n_left >=1
   1353  *n_right >=1
   1354  */
   1355 void OpenGLSurfaceEvaluator::inEvalVStrip(int n_left, REAL u_left, REAL* left_val, int n_right, REAL u_right, REAL* right_val)
   1356 {
   1357   int i,j,k,l;
   1358   REAL botMostV[2];
   1359   typedef REAL REAL3[3];
   1360 
   1361   REAL3* leftXYZ = (REAL3*) malloc(sizeof(REAL3)*n_left);
   1362   assert(leftXYZ);
   1363   REAL3* leftNormal = (REAL3*) malloc(sizeof(REAL3) * n_left);
   1364   assert(leftNormal);
   1365   REAL3* rightXYZ = (REAL3*) malloc(sizeof(REAL3)*n_right);
   1366   assert(rightXYZ);
   1367   REAL3* rightNormal = (REAL3*) malloc(sizeof(REAL3) * n_right);
   1368   assert(rightNormal);
   1369 
   1370   inEvalVLine(n_left, u_left, left_val,  1, leftXYZ, leftNormal);
   1371   inEvalVLine(n_right, u_right, right_val,  1, rightXYZ, rightNormal);
   1372 
   1373 
   1374 
   1375   REAL* botMostXYZ;
   1376   REAL* botMostNormal;
   1377 
   1378   /*
   1379    *the algorithm works by scanning from bot to top.
   1380    *botMostV: the bot most of the remaining verteces (on both left and right).
   1381    *           it could an element of leftVerts or rightVerts.
   1382    *i: leftVerts[i] is the first vertex to the top of botMostV on left line
   1383    *j: rightVerts[j] is the first vertex to the top of botMostV on rightline   */
   1384 
   1385   /*initialize i,j,and botMostV
   1386    */
   1387   if(left_val[0] <= right_val[0])
   1388     {
   1389       i=1;
   1390       j=0;
   1391 
   1392       botMostV[0] = u_left;
   1393       botMostV[1] = left_val[0];
   1394       botMostXYZ = leftXYZ[0];
   1395       botMostNormal = leftNormal[0];
   1396     }
   1397   else
   1398     {
   1399       i=0;
   1400       j=1;
   1401 
   1402       botMostV[0] = u_right;
   1403       botMostV[1] = right_val[0];
   1404 
   1405       botMostXYZ = rightXYZ[0];
   1406       botMostNormal = rightNormal[0];
   1407     }
   1408 
   1409   /*the main loop.
   1410    *the invariance is that:
   1411    *at the beginning of each loop, the meaning of i,j,and botMostV are
   1412    *maintained
   1413    */
   1414   while(1)
   1415     {
   1416       if(i >= n_left) /*case1: no more in left*/
   1417         {
   1418           if(j<n_right-1) /*at least two vertices in right*/
   1419             {
   1420               bgntfan();
   1421 	      glNormal3fv(botMostNormal);
   1422               glVertex3fv(botMostXYZ);
   1423 
   1424               while(j<n_right){
   1425 		glNormal3fv(rightNormal[j]);
   1426 		glVertex3fv(rightXYZ[j]);
   1427 		j++;
   1428 
   1429               }
   1430               endtfan();
   1431             }
   1432           break; /*exit the main loop*/
   1433         }
   1434       else if(j>= n_right) /*case2: no more in right*/
   1435         {
   1436           if(i<n_left-1) /*at least two vertices in left*/
   1437             {
   1438               bgntfan();
   1439 	      glNormal3fv(botMostNormal);
   1440 	      glVertex3fv(botMostXYZ);
   1441 
   1442               for(k=n_left-1; k>=i; k--) /*reverse order for two-side lighting*/
   1443 		{
   1444 		  glNormal3fv(leftNormal[k]);
   1445 		  glVertex3fv(leftXYZ[k]);
   1446 		}
   1447 
   1448               endtfan();
   1449             }
   1450           break; /*exit the main loop*/
   1451         }
   1452       else /* case3: neither is empty, plus the botMostV, there is at least one triangle to output*/
   1453         {
   1454           if(left_val[i] <= right_val[j])
   1455             {
   1456 	      bgntfan();
   1457 
   1458 	      glNormal3fv(rightNormal[j]);
   1459 	      glVertex3fv(rightXYZ[j]);
   1460 
   1461               /*find the last k>=i such that
   1462                *leftverts[k][0] <= rightverts[j][0]
   1463                */
   1464               k=i;
   1465 
   1466               while(k<n_left)
   1467                 {
   1468                   if(left_val[k] > right_val[j])
   1469                     break;
   1470                   k++;
   1471 
   1472                 }
   1473               k--;
   1474 
   1475 
   1476               for(l=k; l>=i; l--)/*the reverse is for two-side lighting*/
   1477                 {
   1478 		  glNormal3fv(leftNormal[l]);
   1479 		  glVertex3fv(leftXYZ[l]);
   1480 
   1481                 }
   1482 	      glNormal3fv(botMostNormal);
   1483 	      glVertex3fv(botMostXYZ);
   1484 
   1485               endtfan();
   1486 
   1487               /*update i and botMostV for next loop
   1488                */
   1489               i = k+1;
   1490 
   1491 	      botMostV[0] = u_left;
   1492 	      botMostV[1] = left_val[k];
   1493 	      botMostNormal = leftNormal[k];
   1494 	      botMostXYZ = leftXYZ[k];
   1495             }
   1496           else /*left_val[i] > right_val[j])*/
   1497             {
   1498 	      bgntfan();
   1499 	      glNormal3fv(leftNormal[i]);
   1500 	      glVertex3fv(leftXYZ[i]);
   1501 
   1502               glNormal3fv(botMostNormal);
   1503 	      glVertex3fv(botMostXYZ);
   1504 
   1505 
   1506               /*find the last k>=j such that
   1507                *rightverts[k][0] < leftverts[i][0]
   1508                */
   1509               k=j;
   1510               while(k< n_right)
   1511                 {
   1512                   if(right_val[k] >= left_val[i])
   1513                     break;
   1514 		  glNormal3fv(rightNormal[k]);
   1515 		  glVertex3fv(rightXYZ[k]);
   1516 
   1517                   k++;
   1518                 }
   1519               endtfan();
   1520 
   1521               /*update j and botMostV for next loop
   1522                */
   1523               j=k;
   1524 	      botMostV[0] = u_right;
   1525 	      botMostV[1] = right_val[j-1];
   1526 
   1527 	      botMostNormal = rightNormal[j-1];
   1528 	      botMostXYZ = rightXYZ[j-1];
   1529             }
   1530         }
   1531     }
   1532   //clean up
   1533   free(leftXYZ);
   1534   free(rightXYZ);
   1535   free(leftNormal);
   1536   free(rightNormal);
   1537 }
   1538 
   1539 /*-----------------------begin evalMachine-------------------*/
   1540 void OpenGLSurfaceEvaluator::inMap2fEM(int which, int k,
   1541 	     REAL ulower,
   1542 	     REAL uupper,
   1543 	     int ustride,
   1544 	     int uorder,
   1545 	     REAL vlower,
   1546 	     REAL vupper,
   1547 	     int vstride,
   1548 	     int vorder,
   1549 	     REAL *ctlPoints)
   1550 {
   1551   int i,j,x;
   1552   surfEvalMachine *temp_em;
   1553   switch(which){
   1554   case 0: //vertex
   1555     vertex_flag = 1;
   1556     temp_em = &em_vertex;
   1557     break;
   1558   case 1: //normal
   1559     normal_flag = 1;
   1560     temp_em = &em_normal;
   1561     break;
   1562   case 2: //color
   1563     color_flag = 1;
   1564     temp_em = &em_color;
   1565     break;
   1566   default:
   1567     texcoord_flag = 1;
   1568     temp_em = &em_texcoord;
   1569     break;
   1570   }
   1571 
   1572   REAL *data = temp_em->ctlPoints;
   1573 
   1574   temp_em->uprime = -1;//initilized
   1575   temp_em->vprime = -1;
   1576 
   1577   temp_em->k = k;
   1578   temp_em->u1 = ulower;
   1579   temp_em->u2 = uupper;
   1580   temp_em->ustride = ustride;
   1581   temp_em->uorder = uorder;
   1582   temp_em->v1 = vlower;
   1583   temp_em->v2 = vupper;
   1584   temp_em->vstride = vstride;
   1585   temp_em->vorder = vorder;
   1586 
   1587   /*copy the contrl points from ctlPoints to global_ev_ctlPoints*/
   1588   for (i=0; i<uorder; i++) {
   1589     for (j=0; j<vorder; j++) {
   1590       for (x=0; x<k; x++) {
   1591 	data[x] = ctlPoints[x];
   1592       }
   1593       ctlPoints += vstride;
   1594       data += k;
   1595     }
   1596     ctlPoints += ustride - vstride * vorder;
   1597   }
   1598 }
   1599 
   1600 void OpenGLSurfaceEvaluator::inDoDomain2WithDerivsEM(surfEvalMachine *em, REAL u, REAL v,
   1601 				REAL *retPoint, REAL *retdu, REAL *retdv)
   1602 {
   1603     int j, row, col;
   1604     REAL the_uprime;
   1605     REAL the_vprime;
   1606     REAL p;
   1607     REAL pdv;
   1608     REAL *data;
   1609 
   1610     if((em->u2 == em->u1) || (em->v2 == em->v1))
   1611 	return;
   1612     the_uprime = (u - em->u1) / (em->u2 - em->u1);
   1613     the_vprime = (v - em->v1) / (em->v2 - em->v1);
   1614 
   1615     /* Compute coefficients for values and derivs */
   1616 
   1617     /* Use already cached values if possible */
   1618     if(em->uprime != the_uprime) {
   1619         inPreEvaluateWithDeriv(em->uorder, the_uprime, em->ucoeff, em->ucoeffDeriv);
   1620 	em->uprime = the_uprime;
   1621     }
   1622     if (em->vprime != the_vprime) {
   1623 	inPreEvaluateWithDeriv(em->vorder, the_vprime, em->vcoeff, em->vcoeffDeriv);
   1624 	em->vprime = the_vprime;
   1625     }
   1626 
   1627     for (j = 0; j < em->k; j++) {
   1628 	data=em->ctlPoints+j;
   1629 	retPoint[j] = retdu[j] = retdv[j] = 0.0;
   1630 	for (row = 0; row < em->uorder; row++)  {
   1631 	    /*
   1632 	    ** Minor optimization.
   1633 	    ** The col == 0 part of the loop is extracted so we don't
   1634 	    ** have to initialize p and pdv to 0.
   1635 	    */
   1636 	    p = em->vcoeff[0] * (*data);
   1637 	    pdv = em->vcoeffDeriv[0] * (*data);
   1638 	    data += em->k;
   1639 	    for (col = 1; col < em->vorder; col++) {
   1640 		/* Incrementally build up p, pdv value */
   1641 		p += em->vcoeff[col] * (*data);
   1642 		pdv += em->vcoeffDeriv[col] * (*data);
   1643 		data += em->k;
   1644 	    }
   1645 	    /* Use p, pdv value to incrementally add up r, du, dv */
   1646 	    retPoint[j] += em->ucoeff[row] * p;
   1647 	    retdu[j] += em->ucoeffDeriv[row] * p;
   1648 	    retdv[j] += em->ucoeff[row] * pdv;
   1649 	}
   1650     }
   1651 }
   1652 
   1653 void OpenGLSurfaceEvaluator::inDoDomain2EM(surfEvalMachine *em, REAL u, REAL v,
   1654 				REAL *retPoint)
   1655 {
   1656     int j, row, col;
   1657     REAL the_uprime;
   1658     REAL the_vprime;
   1659     REAL p;
   1660     REAL *data;
   1661 
   1662     if((em->u2 == em->u1) || (em->v2 == em->v1))
   1663 	return;
   1664     the_uprime = (u - em->u1) / (em->u2 - em->u1);
   1665     the_vprime = (v - em->v1) / (em->v2 - em->v1);
   1666 
   1667     /* Compute coefficients for values and derivs */
   1668 
   1669     /* Use already cached values if possible */
   1670     if(em->uprime != the_uprime) {
   1671         inPreEvaluate(em->uorder, the_uprime, em->ucoeff);
   1672 	em->uprime = the_uprime;
   1673     }
   1674     if (em->vprime != the_vprime) {
   1675 	inPreEvaluate(em->vorder, the_vprime, em->vcoeff);
   1676 	em->vprime = the_vprime;
   1677     }
   1678 
   1679     for (j = 0; j < em->k; j++) {
   1680 	data=em->ctlPoints+j;
   1681 	retPoint[j] = 0.0;
   1682 	for (row = 0; row < em->uorder; row++)  {
   1683 	    /*
   1684 	    ** Minor optimization.
   1685 	    ** The col == 0 part of the loop is extracted so we don't
   1686 	    ** have to initialize p and pdv to 0.
   1687 	    */
   1688 	    p = em->vcoeff[0] * (*data);
   1689 	    data += em->k;
   1690 	    for (col = 1; col < em->vorder; col++) {
   1691 		/* Incrementally build up p, pdv value */
   1692 		p += em->vcoeff[col] * (*data);
   1693 		data += em->k;
   1694 	    }
   1695 	    /* Use p, pdv value to incrementally add up r, du, dv */
   1696 	    retPoint[j] += em->ucoeff[row] * p;
   1697 	}
   1698     }
   1699 }
   1700 
   1701 
   1702 void OpenGLSurfaceEvaluator::inDoEvalCoord2EM(REAL u, REAL v)
   1703 {
   1704   REAL temp_vertex[5];
   1705   REAL temp_normal[3];
   1706   REAL temp_color[4];
   1707   REAL temp_texcoord[4];
   1708 
   1709   if(texcoord_flag)
   1710     {
   1711       inDoDomain2EM(&em_texcoord, u,v, temp_texcoord);
   1712       texcoordCallBack(temp_texcoord, userData);
   1713     }
   1714   if(color_flag)
   1715     {
   1716       inDoDomain2EM(&em_color, u,v, temp_color);
   1717       colorCallBack(temp_color, userData);
   1718     }
   1719 
   1720   if(normal_flag) //there is a normla map
   1721     {
   1722       inDoDomain2EM(&em_normal, u,v, temp_normal);
   1723       normalCallBack(temp_normal, userData);
   1724 
   1725       if(vertex_flag)
   1726 	{
   1727 	  inDoDomain2EM(&em_vertex, u,v,temp_vertex);
   1728 	  if(em_vertex.k == 4)
   1729 	    {
   1730 	      temp_vertex[0] /= temp_vertex[3];
   1731 	      temp_vertex[1] /= temp_vertex[3];
   1732 	      temp_vertex[2] /= temp_vertex[3];
   1733 	    }
   1734           temp_vertex[3]=u;
   1735           temp_vertex[4]=v;
   1736 	  vertexCallBack(temp_vertex, userData);
   1737 	}
   1738     }
   1739   else if(auto_normal_flag) //no normal map but there is a normal callbackfunctin
   1740     {
   1741       REAL du[4];
   1742       REAL dv[4];
   1743 
   1744       /*compute homegeneous point and partial derivatives*/
   1745       inDoDomain2WithDerivsEM(&em_vertex, u,v,temp_vertex,du,dv);
   1746 
   1747       if(em_vertex.k ==4)
   1748 	inComputeFirstPartials(temp_vertex, du, dv);
   1749 
   1750 #ifdef AVOID_ZERO_NORMAL
   1751       if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
   1752 	{
   1753 
   1754 	  REAL tempdu[4];
   1755 	  REAL tempdata[4];
   1756 	  REAL u1 = em_vertex.u1;
   1757 	  REAL u2 = em_vertex.u2;
   1758 	  if(u-MYDELTA*(u2-u1) < u1)
   1759 	    u = u+ MYDELTA*(u2-u1);
   1760 	  else
   1761 	    u = u-MYDELTA*(u2-u1);
   1762 	  inDoDomain2WithDerivsEM(&em_vertex,u,v, tempdata, tempdu, dv);
   1763 
   1764 	  if(em_vertex.k ==4)
   1765 	    inComputeFirstPartials(temp_vertex, du, dv);
   1766 	}
   1767       else if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
   1768 	{
   1769 	  REAL tempdv[4];
   1770 	  REAL tempdata[4];
   1771 	  REAL v1 = em_vertex.v1;
   1772 	  REAL v2 = em_vertex.v2;
   1773 	  if(v-MYDELTA*(v2-v1) < v1)
   1774 	    v = v+ MYDELTA*(v2-v1);
   1775 	  else
   1776 	    v = v-MYDELTA*(v2-v1);
   1777 	  inDoDomain2WithDerivsEM(&em_vertex,u,v, tempdata, du, tempdv);
   1778 
   1779 	  if(em_vertex.k ==4)
   1780 	    inComputeFirstPartials(temp_vertex, du, dv);
   1781 	}
   1782 #endif
   1783 
   1784       /*compute normal*/
   1785       switch(em_vertex.k){
   1786       case 3:
   1787 
   1788 	inComputeNormal2(du, dv, temp_normal);
   1789 	break;
   1790       case 4:
   1791 
   1792 //	inComputeFirstPartials(temp_vertex, du, dv);
   1793 	inComputeNormal2(du, dv, temp_normal);
   1794 
   1795 	/*transform the homegeneous coordinate of retPoint into inhomogenous one*/
   1796 	temp_vertex[0] /= temp_vertex[3];
   1797 	temp_vertex[1] /= temp_vertex[3];
   1798 	temp_vertex[2] /= temp_vertex[3];
   1799 	break;
   1800       }
   1801       normalCallBack(temp_normal, userData);
   1802       temp_vertex[3] = u;
   1803       temp_vertex[4] = v;
   1804       vertexCallBack(temp_vertex, userData);
   1805 
   1806     }/*end if auto_normal*/
   1807   else //no normal map, and no normal callback function
   1808     {
   1809       if(vertex_flag)
   1810 	{
   1811 	  inDoDomain2EM(&em_vertex, u,v,temp_vertex);
   1812 	  if(em_vertex.k == 4)
   1813 	    {
   1814 	      temp_vertex[0] /= temp_vertex[3];
   1815 	      temp_vertex[1] /= temp_vertex[3];
   1816 	      temp_vertex[2] /= temp_vertex[3];
   1817 	    }
   1818           temp_vertex[3] = u;
   1819           temp_vertex[4] = v;
   1820 	  vertexCallBack(temp_vertex, userData);
   1821 	}
   1822     }
   1823 }
   1824 
   1825 
   1826 void OpenGLSurfaceEvaluator::inBPMEvalEM(bezierPatchMesh* bpm)
   1827 {
   1828   int i,j,k;
   1829   float u,v;
   1830 
   1831   int ustride;
   1832   int vstride;
   1833 
   1834 #ifdef USE_LOD
   1835   if(bpm->bpatch != NULL)
   1836     {
   1837       bezierPatch* p=bpm->bpatch;
   1838       ustride = p->dimension * p->vorder;
   1839       vstride = p->dimension;
   1840 
   1841       glMap2f( (p->dimension == 3)? GL_MAP2_VERTEX_3 : GL_MAP2_VERTEX_4,
   1842 	      p->umin,
   1843 	      p->umax,
   1844 	      ustride,
   1845 	      p->uorder,
   1846 	      p->vmin,
   1847 	      p->vmax,
   1848 	      vstride,
   1849 	      p->vorder,
   1850 	      p->ctlpoints);
   1851 
   1852 
   1853 /*
   1854     inMap2fEM(0, p->dimension,
   1855 	  p->umin,
   1856 	  p->umax,
   1857 	  ustride,
   1858 	  p->uorder,
   1859 	  p->vmin,
   1860 	  p->vmax,
   1861 	  vstride,
   1862 	  p->vorder,
   1863 	  p->ctlpoints);
   1864 */
   1865     }
   1866 #else
   1867 
   1868   if(bpm->bpatch != NULL){
   1869     bezierPatch* p = bpm->bpatch;
   1870     ustride = p->dimension * p->vorder;
   1871     vstride = p->dimension;
   1872     inMap2fEM(0, p->dimension,
   1873 	  p->umin,
   1874 	  p->umax,
   1875 	  ustride,
   1876 	  p->uorder,
   1877 	  p->vmin,
   1878 	  p->vmax,
   1879 	  vstride,
   1880 	  p->vorder,
   1881 	  p->ctlpoints);
   1882   }
   1883   if(bpm->bpatch_normal != NULL){
   1884     bezierPatch* p = bpm->bpatch_normal;
   1885     ustride = p->dimension * p->vorder;
   1886     vstride = p->dimension;
   1887     inMap2fEM(1, p->dimension,
   1888 	  p->umin,
   1889 	  p->umax,
   1890 	  ustride,
   1891 	  p->uorder,
   1892 	  p->vmin,
   1893 	  p->vmax,
   1894 	  vstride,
   1895 	  p->vorder,
   1896 	  p->ctlpoints);
   1897   }
   1898   if(bpm->bpatch_color != NULL){
   1899     bezierPatch* p = bpm->bpatch_color;
   1900     ustride = p->dimension * p->vorder;
   1901     vstride = p->dimension;
   1902     inMap2fEM(2, p->dimension,
   1903 	  p->umin,
   1904 	  p->umax,
   1905 	  ustride,
   1906 	  p->uorder,
   1907 	  p->vmin,
   1908 	  p->vmax,
   1909 	  vstride,
   1910 	  p->vorder,
   1911 	  p->ctlpoints);
   1912   }
   1913   if(bpm->bpatch_texcoord != NULL){
   1914     bezierPatch* p = bpm->bpatch_texcoord;
   1915     ustride = p->dimension * p->vorder;
   1916     vstride = p->dimension;
   1917     inMap2fEM(3, p->dimension,
   1918 	  p->umin,
   1919 	  p->umax,
   1920 	  ustride,
   1921 	  p->uorder,
   1922 	  p->vmin,
   1923 	  p->vmax,
   1924 	  vstride,
   1925 	  p->vorder,
   1926 	  p->ctlpoints);
   1927   }
   1928 #endif
   1929 
   1930 
   1931   k=0;
   1932   for(i=0; i<bpm->index_length_array; i++)
   1933     {
   1934 #ifdef USE_LOD
   1935       if(bpm->type_array[i] == GL_POLYGON) //a mesh
   1936 	{
   1937 	  GLfloat *temp = bpm->UVarray+k;
   1938 	  GLfloat u0 = temp[0];
   1939 	  GLfloat v0 = temp[1];
   1940 	  GLfloat u1 = temp[2];
   1941 	  GLfloat v1 = temp[3];
   1942 	  GLint nu = (GLint) ( temp[4]);
   1943 	  GLint nv = (GLint) ( temp[5]);
   1944 	  GLint umin = (GLint) ( temp[6]);
   1945 	  GLint vmin = (GLint) ( temp[7]);
   1946 	  GLint umax = (GLint) ( temp[8]);
   1947 	  GLint vmax = (GLint) ( temp[9]);
   1948 
   1949 	  glMapGrid2f(LOD_eval_level*nu, u0, u1, LOD_eval_level*nv, v0, v1);
   1950 	  glEvalMesh2(GL_FILL, LOD_eval_level*umin, LOD_eval_level*umax, LOD_eval_level*vmin, LOD_eval_level*vmax);
   1951 	}
   1952       else
   1953 	{
   1954 	  LOD_eval(bpm->length_array[i], bpm->UVarray+k, bpm->type_array[i],
   1955 		   0
   1956 		   );
   1957 	}
   1958 	  k+= 2*bpm->length_array[i];
   1959 
   1960 #else //undef  USE_LOD
   1961 
   1962 #ifdef CRACK_TEST
   1963 if(  bpm->bpatch->umin == 2 &&   bpm->bpatch->umax == 3
   1964   && bpm->bpatch->vmin ==2 &&    bpm->bpatch->vmax == 3)
   1965 {
   1966 REAL vertex[4];
   1967 REAL normal[4];
   1968 #ifdef DEBUG
   1969 printf("***number ****1\n");
   1970 #endif
   1971 
   1972 beginCallBack(GL_QUAD_STRIP, NULL);
   1973 inDoEvalCoord2EM(3.0, 3.0);
   1974 inDoEvalCoord2EM(2.0, 3.0);
   1975 inDoEvalCoord2EM(3.0, 2.7);
   1976 inDoEvalCoord2EM(2.0, 2.7);
   1977 inDoEvalCoord2EM(3.0, 2.0);
   1978 inDoEvalCoord2EM(2.0, 2.0);
   1979 endCallBack(NULL);
   1980 
   1981 beginCallBack(GL_TRIANGLE_STRIP, NULL);
   1982 inDoEvalCoord2EM(2.0, 3.0);
   1983 inDoEvalCoord2EM(2.0, 2.0);
   1984 inDoEvalCoord2EM(2.0, 2.7);
   1985 endCallBack(NULL);
   1986 
   1987 }
   1988 if(  bpm->bpatch->umin == 1 &&   bpm->bpatch->umax == 2
   1989   && bpm->bpatch->vmin ==2 &&    bpm->bpatch->vmax == 3)
   1990 {
   1991 #ifdef DEBUG
   1992 printf("***number 3\n");
   1993 #endif
   1994 beginCallBack(GL_QUAD_STRIP, NULL);
   1995 inDoEvalCoord2EM(2.0, 3.0);
   1996 inDoEvalCoord2EM(1.0, 3.0);
   1997 inDoEvalCoord2EM(2.0, 2.3);
   1998 inDoEvalCoord2EM(1.0, 2.3);
   1999 inDoEvalCoord2EM(2.0, 2.0);
   2000 inDoEvalCoord2EM(1.0, 2.0);
   2001 endCallBack(NULL);
   2002 
   2003 beginCallBack(GL_TRIANGLE_STRIP, NULL);
   2004 inDoEvalCoord2EM(2.0, 2.3);
   2005 inDoEvalCoord2EM(2.0, 2.0);
   2006 inDoEvalCoord2EM(2.0, 3.0);
   2007 endCallBack(NULL);
   2008 
   2009 }
   2010 return;
   2011 #endif //CRACK_TEST
   2012 
   2013       beginCallBack(bpm->type_array[i], userData);
   2014 
   2015       for(j=0; j<bpm->length_array[i]; j++)
   2016 	{
   2017 	  u = bpm->UVarray[k];
   2018 	  v = bpm->UVarray[k+1];
   2019 #ifdef USE_LOD
   2020           LOD_EVAL_COORD(u,v);
   2021 //	  glEvalCoord2f(u,v);
   2022 #else
   2023 
   2024 #ifdef  GENERIC_TEST
   2025           float temp_normal[3];
   2026           float temp_vertex[3];
   2027           if(temp_signal == 0)
   2028 	    {
   2029 	      gTessVertexSphere(u,v, temp_normal, temp_vertex);
   2030 //printf("normal=(%f,%f,%f)\n", temp_normal[0], temp_normal[1], temp_normal[2])//printf("veretx=(%f,%f,%f)\n", temp_vertex[0], temp_vertex[1], temp_vertex[2]);
   2031               normalCallBack(temp_normal, userData);
   2032 	      vertexCallBack(temp_vertex, userData);
   2033 	    }
   2034           else if(temp_signal == 1)
   2035 	    {
   2036 	      gTessVertexCyl(u,v, temp_normal, temp_vertex);
   2037 //printf("normal=(%f,%f,%f)\n", temp_normal[0], temp_normal[1], temp_normal[2])//printf("veretx=(%f,%f,%f)\n", temp_vertex[0], temp_vertex[1], temp_vertex[2]);
   2038               normalCallBack(temp_normal, userData);
   2039 	      vertexCallBack(temp_vertex, userData);
   2040 	    }
   2041 	  else
   2042 #endif //GENERIC_TEST
   2043 
   2044 	    inDoEvalCoord2EM(u,v);
   2045 
   2046 #endif //USE_LOD
   2047 
   2048 	  k += 2;
   2049 	}
   2050       endCallBack(userData);
   2051 
   2052 #endif //USE_LOD
   2053     }
   2054 }
   2055 
   2056 void OpenGLSurfaceEvaluator::inBPMListEvalEM(bezierPatchMesh* list)
   2057 {
   2058   bezierPatchMesh* temp;
   2059   for(temp = list; temp != NULL; temp = temp->next)
   2060     {
   2061       inBPMEvalEM(temp);
   2062     }
   2063 }
   2064 
   2065