src/glu/sgi/libnurbs/interface/insurfeval.cc

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