src/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 */
  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