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
64static 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
92void 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);*/
107glBegin(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();
136glEnd();
137    }
138}
139
140void 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
183extern float xmin, xmax, ymin, ymax, zmin, zmax; /*bounding box*/
184extern int temp_signal;
185
186static 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
210static 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
238void OpenGLSurfaceEvaluator::inBPMListEval(bezierPatchMesh* list)
239{
240  bezierPatchMesh* temp;
241  for(temp = list; temp != NULL; temp = temp->next)
242    {
243      inBPMEval(temp);
244    }
245}
246
247void 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
271if(  global_ev_u1 ==2 &&   global_ev_u2 == 3
272  && global_ev_v1 ==2 &&   global_ev_v2 == 3)
273{
274REAL vertex[4];
275REAL normal[4];
276#ifdef DEBUG
277printf("***number 1\n");
278#endif
279
280beginCallBack(GL_QUAD_STRIP, NULL);
281inEvalCoord2f(3.0, 3.0);
282inEvalCoord2f(2.0, 3.0);
283inEvalCoord2f(3.0, 2.7);
284inEvalCoord2f(2.0, 2.7);
285inEvalCoord2f(3.0, 2.0);
286inEvalCoord2f(2.0, 2.0);
287endCallBack(NULL);
288
289
290beginCallBack(GL_TRIANGLE_STRIP, NULL);
291inEvalCoord2f(2.0, 3.0);
292inEvalCoord2f(2.0, 2.0);
293inEvalCoord2f(2.0, 2.7);
294endCallBack(NULL);
295
296}
297
298/*
299if(  global_ev_u1 ==2 &&   global_ev_u2 == 3
300  && global_ev_v1 ==1 &&   global_ev_v2 == 2)
301{
302#ifdef DEBUG
303printf("***number 2\n");
304#endif
305beginCallBack(GL_QUAD_STRIP);
306inEvalCoord2f(2.0, 2.0);
307inEvalCoord2f(2.0, 1.0);
308inEvalCoord2f(3.0, 2.0);
309inEvalCoord2f(3.0, 1.0);
310endCallBack();
311}
312*/
313if(  global_ev_u1 ==1 &&   global_ev_u2 == 2
314  && global_ev_v1 ==2 &&   global_ev_v2 == 3)
315{
316#ifdef DEBUG
317printf("***number 3\n");
318#endif
319beginCallBack(GL_QUAD_STRIP, NULL);
320inEvalCoord2f(2.0, 3.0);
321inEvalCoord2f(1.0, 3.0);
322inEvalCoord2f(2.0, 2.3);
323inEvalCoord2f(1.0, 2.3);
324inEvalCoord2f(2.0, 2.0);
325inEvalCoord2f(1.0, 2.0);
326endCallBack(NULL);
327
328beginCallBack(GL_TRIANGLE_STRIP, NULL);
329inEvalCoord2f(2.0, 2.3);
330inEvalCoord2f(2.0, 2.0);
331inEvalCoord2f(2.0, 3.0);
332endCallBack(NULL);
333
334}
335return;
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
362void 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
375void 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 */
389void 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
400void 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
444void 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 */
501void 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 */
518void 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 */
541void 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 */
621void 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 */
684void 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 */
749void 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
809void 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
838void 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
868void 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
903void 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 */
949void 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 */
1022void 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 */
1066void 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
1125void 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];
1130inPreEvaluateBV_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
1144void 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];
1149inPreEvaluateBU_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 */
1166void 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 */
1355void 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-------------------*/
1540void 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
1600void 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
1653void 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
1702void 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
1826void 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
1963if(  bpm->bpatch->umin == 2 &&   bpm->bpatch->umax == 3
1964  && bpm->bpatch->vmin ==2 &&    bpm->bpatch->vmax == 3)
1965{
1966REAL vertex[4];
1967REAL normal[4];
1968#ifdef DEBUG
1969printf("***number ****1\n");
1970#endif
1971
1972beginCallBack(GL_QUAD_STRIP, NULL);
1973inDoEvalCoord2EM(3.0, 3.0);
1974inDoEvalCoord2EM(2.0, 3.0);
1975inDoEvalCoord2EM(3.0, 2.7);
1976inDoEvalCoord2EM(2.0, 2.7);
1977inDoEvalCoord2EM(3.0, 2.0);
1978inDoEvalCoord2EM(2.0, 2.0);
1979endCallBack(NULL);
1980
1981beginCallBack(GL_TRIANGLE_STRIP, NULL);
1982inDoEvalCoord2EM(2.0, 3.0);
1983inDoEvalCoord2EM(2.0, 2.0);
1984inDoEvalCoord2EM(2.0, 2.7);
1985endCallBack(NULL);
1986
1987}
1988if(  bpm->bpatch->umin == 1 &&   bpm->bpatch->umax == 2
1989  && bpm->bpatch->vmin ==2 &&    bpm->bpatch->vmax == 3)
1990{
1991#ifdef DEBUG
1992printf("***number 3\n");
1993#endif
1994beginCallBack(GL_QUAD_STRIP, NULL);
1995inDoEvalCoord2EM(2.0, 3.0);
1996inDoEvalCoord2EM(1.0, 3.0);
1997inDoEvalCoord2EM(2.0, 2.3);
1998inDoEvalCoord2EM(1.0, 2.3);
1999inDoEvalCoord2EM(2.0, 2.0);
2000inDoEvalCoord2EM(1.0, 2.0);
2001endCallBack(NULL);
2002
2003beginCallBack(GL_TRIANGLE_STRIP, NULL);
2004inDoEvalCoord2EM(2.0, 2.3);
2005inDoEvalCoord2EM(2.0, 2.0);
2006inDoEvalCoord2EM(2.0, 3.0);
2007endCallBack(NULL);
2008
2009}
2010return;
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
2056void OpenGLSurfaceEvaluator::inBPMListEvalEM(bezierPatchMesh* list)
2057{
2058  bezierPatchMesh* temp;
2059  for(temp = list; temp != NULL; temp = temp->next)
2060    {
2061      inBPMEvalEM(temp);
2062    }
2063}
2064
2065