~ubuntu-branches/ubuntu/quantal/mesa-glw/quantal

« back to all changes in this revision

Viewing changes to src/glu/sgi/libnurbs/interface/insurfeval.cc

  • Committer: Bazaar Package Importer
  • Author(s): Morten Kjeldgaard
  • Date: 2008-05-06 16:19:15 UTC
  • Revision ID: james.westby@ubuntu.com-20080506161915-uynz7nftmfixu6bq
Tags: upstream-7.0.3
ImportĀ upstreamĀ versionĀ 7.0.3

Show diffs side-by-side

added added

removed removed

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