1 /* -------------------------------------------------------------------- */
2 /* CALCULIX */
3 /* - GRAPHICAL INTERFACE - */
4 /* */
5 /* A 3-dimensional pre- and post-processor for finite elements */
6 /* Copyright (C) 1996 Klaus Wittig */
7 /* */
8 /* This program is free software; you can reistribute it and/or */
9 /* modify it under the terms of the GNU General Public License as */
10 /* published by the Free Software Foundation; version 2 of */
11 /* the License. */
12 /* */
13 /* This program is distributed in the hope that it will be useful, */
14 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
15 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
16 /* GNU General Public License for more details. */
17 /* */
18 /* You should have received a copy of the GNU General Public License */
19 /* along with this program; if not, write to the Free Software */
20 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
21 /* -------------------------------------------------------------------- */
22
23 /* TO DO
24 */
25
26 #define TEST 0 /* debugging */
27 #define TEST2 0 /* debugging */
28 #define PRINT_SNL 0 /* debugging */
29 #define FAST_SHADING 0 /* scip the line division oriented resolution of filled surfaces and use the coarsest possibility */
30
31 #include "uselibSNL.h"
32
33 #include <iostream>
34
35 #include <snlSurface.h>
36 #include <snlUtil.h>
37 #include <snlNurbsCommon.h>
38
39 #define ITER_TOL 1.e-6
40 #define COS_TOL 1.e-2
41 #define MAX_PASS 500
42 #define PROJ_SENSITIVITY 10
43 #define MAX_RETURNED_PNTS_PER_PNT 10
44 #define MIN_DIST 1.e-5
45 #define MAXVALUE 1.e99
46
47 #ifdef WIN32
48 extern "C" char printFlag; /* printf 1:on 0:off */
49 extern "C" Nurbs *nurbs;
50 extern "C" Points *point;
51 extern "C" Lines *line;
52 extern "C" Lcmb *lcmb;
53 extern "C" Gsur *surf;
54 extern "C" sem_t sem_g;
55 extern "C" sem_t sem_n;
56 #else
57 extern char printFlag; /* printf 1:on 0:off */
58 extern Nurbs *nurbs;
59 extern Points *point;
60 extern Lines *line;
61 extern Lcmb *lcmb;
62 extern Gsur *surf;
63 extern sem_t sem_g;
64 extern sem_t sem_n;
65 #endif
66 #undef max
67 #undef min
68
69
70
71 typedef struct {
72 int n;
73 double paramU[MAX_RETURNED_PNTS_PER_PNT];
74 double paramV[MAX_RETURNED_PNTS_PER_PNT];
75 double dist[MAX_RETURNED_PNTS_PER_PNT];
76 }PntProj;
77
78
vl_result(double * A,double * B,double * C)79 void vl_result( double *A, double *B, double *C )
80 /**********************************************************/
81 /* Vektorbetrag: C = Vektor(B)-Vektor(A) == Vector(AB)*/
82 /**********************************************************/
83 {
84 register int i;
85
86 for (i=0; i<3; i++)
87 C[i]=B[i]-A[i];
88 }
vl_prod(double * A,double * B,double * C)89 void vl_prod( double *A, double *B, double *C )
90 /*********************************************************/
91 /* Vektormultiplikation: C = A x B */
92 /*********************************************************/
93 {
94 C[0]=A[1]*B[2]-A[2]*B[1];
95 C[1]=A[2]*B[0]-A[0]*B[2];
96 C[2]=A[0]*B[1]-A[1]*B[0];
97 }
98 /*********************************************************/
99 /* Scalar Product: C = A o B */
100 /*********************************************************/
vl_sprod(double * A,double * B,double * C)101 void vl_sprod( double *A, double *B, double *C )
102 {
103 *C=A[0]*B[0]+A[1]*B[1]+A[2]*B[2];
104 }
vl_norm(double * A,double * C)105 double vl_norm( double *A, double *C )
106 /*********************************************************/
107 /* Vector(B)=Vektor(A)/|A| */
108 /*********************************************************/
109 {
110 int i;
111 double B;
112
113 B=0.;
114 for ( i=0; i<3; i++)
115 B=B+A[i]*A[i];
116
117 B = sqrt(B);
118 if(B==0.) {C[0]=C[1]=C[2]=0.; return(B);}
119 for ( i=0; i<3; i++)
120 C[i]=A[i]/B;
121 return(B);
122 }
vl_betrag(double * a)123 double vl_betrag(double *a)
124 /* ********************************************************* */
125 /* laenge von Vektor a */
126 /* ********************************************************* */
127 {
128 register int i;
129 double b;
130 b=0.;
131
132 for (i=0; i<3; i++)
133 {
134 b=b+a[i]*a[i];
135 }
136 b=sqrt(b);
137 return(b);
138 }
vl_scal(double * A,double * B,double * C)139 void vl_scal( double *A, double *B, double *C )
140 /**********************************************************/
141 /* Vektormultiplikation: vektor(C) = scalar(A)*Vektor(B) */
142 /**********************************************************/
143 {
144 register int i;
145
146 for (i=0; i<3; i++){
147 C[i]= *A * B[i];
148 }
149 }
vl_add(double * A,double * B,double * C)150 void vl_add( double *A, double *B, double *C )
151 /**********************************************************/
152 /* Vektoraddition: C = Vektor(B)+Vektor(A) */
153 /**********************************************************/
154 {
155 register int i;
156
157 for (i=0; i<3; i++){
158 C[i]=B[i]+A[i];
159 }
160 }
161
162 /**********************************************************/
163 /* Coordinate transformation : A = M__ * P + t_ */
164 /**********************************************************/
vl_trans(double M[][3],double * t,double * P,double * A)165 void vl_trans(double M [][3],double * t, double * P, double * A)
166 {
167 int i,j;
168 for(i=0;i<3;i++){
169 A[i]=0;
170 for(j=0;j<3;j++){
171 A[i]=A[i]+M[i][j]*P[j];
172 }
173 }
174 vl_add(A,t,A);
175 }
176
177 //==============================================================//
178 // a function to generate a torodial NURBS surface //
179 // input: p1 ^= origin, p2 ^= second point on symmetry axis //
180 // r1 ^= major radius, r2 ^= minur radius, mySurf ^= result //
181 //==============================================================//
makeTorus(double * p1,double * p2,double r1,double r2,BSplineSurface * mySurf)182 void makeTorus(double *p1, double *p2, double r1, double r2, BSplineSurface * mySurf)
183 {
184 // construct a B-Spline cirlce that can later be rotated to form a torus
185 int l,m,i,j;
186 double *resKnot = NULL;
187 double *resWeight = NULL;
188 double *rescX=NULL, *rescY = NULL, *rescZ = NULL;
189
190 BSplineCurve myCurve;
191 myCurve.nPol = 9;
192 myCurve.nKnt = 12;
193 myCurve.deg = 2;
194
195 if((rescX = (double *)realloc((double *)rescX, (myCurve.nPol)*sizeof(double))) == NULL )
196 { printf(" ERROR: realloc failure rescX\n\n"); }
197
198 if((rescY = (double *)realloc((double *)rescY, (myCurve.nPol)*sizeof(double))) == NULL )
199 { printf(" ERROR: realloc failure rescY\n\n"); }
200
201 if((rescZ = (double *)realloc((double *)rescZ, (myCurve.nPol)*sizeof(double))) == NULL )
202 { printf(" ERROR: realloc failure rescZ\n\n"); }
203
204 double circlePoint[9][3]= {{1,0,0},{1,0,-1},{0,0,-1},
205 {-1,0,-1},{-1,0,0},{-1,0,1},
206 {0,0,1},{1,0,1},{1,0,0}};
207
208 double wCircle[9]={1,0.707107,1,0.707107,1,0.707107,1,0.707107,1};
209
210 double kCircle[12]={0,0,0,0.25,0.25,0.5,0.5,0.75,0.75,1,1,1};
211
212 double circlePointScaled[9][3];
213 double resPol[9][3];
214
215 double xAxis[3] ,zAxis[3];//local coordinate system
216 double p1p2[3];
217 double vIndep[3] = {1,0,0}, scalProd;
218
219 //calculate the local coordinate system
220 vl_result(p1,p2,p1p2);
221 vl_norm(p1p2,zAxis);//scale to length 1
222 vl_sprod(zAxis,vIndep,&scalProd);
223 if(abs(scalProd)>0.99){// not linear independent
224 vIndep[0] = 0;
225 vIndep[1] = 1;
226 vIndep[2] = 0;
227 }
228 vl_prod(zAxis,vIndep,xAxis);
229 vl_norm(xAxis,xAxis);
230
231 //Matrix M: (x1, y1, z1)
232 // (x2, y2, z2)
233 // (x3, y3, z3)
234 // [column,row] bzw. [zeile,spalte]
235 double M [3][3];//matrix for new coordinate system
236 double t[3];//translation for new coordinate system
237
238 //calculate the transformation Matrix M and translation vector t
239 for(i=0;i<3;i++)
240 {
241 M[i][0]=xAxis[i];//wirte x-Axis to fisrt column
242 M[i][1]=0;//wirte y-Axis to second column
243 M[i][2]=zAxis[i];//wirte z-Axis to third column
244 t[i] = p1[i]+r1*xAxis[i];//set translation to origin
245 }
246
247 //scale the circular curve to the minor radius r2
248 for(l=0;l<9;l++){
249 for(m=0;m<3;m++){
250 circlePointScaled[l][m] = r2 * circlePoint[l][m];
251 }
252 }
253
254 //calcutale position of control points
255 for(j=0;j<9;j++){//calculate control points of circle
256 vl_trans(M,t,circlePointScaled[j], resPol[j]);//calculate controll points of the rotated surface
257 rescX[j]=resPol[j][0];
258 rescY[j]=resPol[j][1];
259 rescZ[j]=resPol[j][2];
260 }
261
262 //write result to BSplineCurve struct
263 myCurve.k = kCircle;
264 myCurve.w = wCircle;
265 myCurve.cX =rescX;
266 myCurve.cY =rescY;
267 myCurve.cZ =rescZ;
268
269 //rotate circular profile curve to get a Torodial NURBS surface
270 rotateBSpline(p1, p2, &myCurve, mySurf);
271 }
272
273 //==============================================================//
274 // a function to generate a torodial NURBS surface //
275 // input: p1 ^= origin, p2 ^= second point on symmetry axis //
276 // r1 ^= major radius, r2 ^= minur radius, mySurf ^= result //
277 //==============================================================//
translateBSpline(double * p1,double * p2,BSplineCurve * myCurve,BSplineSurface * mySurf)278 void translateBSpline(double *p1, double *p2, BSplineCurve * myCurve ,BSplineSurface * mySurf)
279 {
280 double p1p2[3];
281 int i,j;
282
283 double *resVKnot = NULL;
284 double *resUKnot = NULL;
285 double *resWeight = NULL;
286 double *rescX=NULL, *rescY = NULL, *rescZ = NULL;
287
288 if((resVKnot = (double *)realloc((double *)resVKnot, 4*sizeof(double))) == NULL )
289 { printf(" ERROR: realloc failure resVKnot\n\n"); }
290
291 if((resUKnot = (double *)realloc((double *)resUKnot, (myCurve->nKnt)*sizeof(double))) == NULL )
292 { printf(" ERROR: realloc failure resUKnot\n\n"); }
293
294 if((resWeight = (double *)realloc((double *)resWeight, (myCurve->nPol)*2*sizeof(double))) == NULL )
295 { printf(" ERROR: realloc failure resWeights\n\n"); }
296
297 if((rescX = (double *)realloc((double *)rescX, (myCurve->nPol)*2*sizeof(double))) == NULL )
298 { printf(" ERROR: realloc failure rescX\n\n"); }
299
300 if((rescY = (double *)realloc((double *)rescY, (myCurve->nPol)*2*sizeof(double))) == NULL )
301 { printf(" ERROR: realloc failure rescY\n\n"); }
302
303 if((rescZ = (double *)realloc((double *)rescZ, (myCurve->nPol)*2*sizeof(double))) == NULL )
304 { printf(" ERROR: realloc failure rescZ\n\n"); }
305
306 //calcutale control points
307 vl_result(p1,p2,p1p2);
308 i=0;
309 for(j=0;j<(myCurve->nPol);j++)
310 {
311 rescX[i]=(myCurve->cX[j]);
312 rescY[i]=(myCurve->cY[j]);
313 rescZ[i]=(myCurve->cZ[j]);
314 resWeight[i] = myCurve->w[j];//calculate weights
315 i++;
316 rescX[i]=(myCurve->cX[j])+p1p2[0];
317 rescY[i]=(myCurve->cY[j])+p1p2[1];
318 rescZ[i]=(myCurve->cZ[j])+p1p2[2];
319 resWeight[i] = myCurve->w[j];//calculate weights
320 i++;
321 }
322
323 //calculate u-Knots
324 for(j = 0;j<(myCurve->nKnt);j++){
325 resUKnot[j] = myCurve->k[j];
326 }
327 //calculate v-Knots
328 resVKnot[0] = 0;
329 resVKnot[1] = 0;
330 resVKnot[2] = 1;
331 resVKnot[3] = 1;
332
333 //write result to struct
334 mySurf->uDeg = myCurve->deg;
335 mySurf->vDeg = 1;
336 mySurf->nUPol = myCurve->nPol;
337 mySurf->nVPol = 2;
338 mySurf->nUKnt = myCurve->nKnt;
339 mySurf->nVKnt = 4;
340 mySurf->uKnt = resUKnot;
341 mySurf->vKnt = resVKnot;
342
343 mySurf->weights = resWeight;
344 mySurf->cX =rescX;
345 mySurf->cY =rescY;
346 mySurf->cZ =rescZ;
347 }
348
349
350
rotateBSpline(double * p1,double * p2,BSplineCurve * myCurve,BSplineSurface * mySurf)351 void rotateBSpline(double *p1, double *p2, BSplineCurve * myCurve ,BSplineSurface * mySurf)
352 {
353 double xAxis[3], yAxis[3] ,zAxis[3];//local coordinate system
354 double origin[3];
355 double p1cPoint[3], cPointTmp[3];
356 double tmpVec[3],p1p2[3];
357 double scalProd, length;
358 double radius;
359
360 int i,j,l,m;
361
362 //Matrix M: (x1, y1, z1)
363 // (x2, y2, z2)
364 // (x3, y3, z3)
365 // [column,row] bzw. [zeile,spalte]
366 double M [3][3];//matrix for new coordinate system
367 double t[3];//translation for new coordinate system
368
369 //Define BSpline-circle
370 double circlePoint[9][3]= {{1,0,0},{1,1,0},{0,1,0},
371 {-1,1,0},{-1,0,0},{-1,-1,0},
372 {0,-1,0},{1,-1,0},{1,0,0}};
373 double circlePointScaled[9][3];
374
375 double wCircle[9]={1,0.707107,1,0.707107,1,0.707107,1,0.707107,1};
376
377 double kCircle[12]={0,0,0,0.25,0.25,0.5,0.5,0.75,0.75,1,1,1};
378
379 double *resVKnot = NULL;
380 double *resUKnot = NULL;
381 double *resWeight = NULL;
382 double *rescX=NULL, *rescY = NULL, *rescZ = NULL;
383
384 if((resVKnot = (double *)realloc((double *)resVKnot, 12*sizeof(double))) == NULL )
385 { printf(" ERROR: realloc failure resVKnot\n\n"); }
386
387 if((resUKnot = (double *)realloc((double *)resUKnot, (myCurve->nKnt)*sizeof(double))) == NULL )
388 { printf(" ERROR: realloc failure resUKnot\n\n"); }
389
390 if((resWeight = (double *)realloc((double *)resWeight, (myCurve->nPol)*9*sizeof(double))) == NULL )
391 { printf(" ERROR: realloc failure resWeights\n\n"); }
392
393 if((rescX = (double *)realloc((double *)rescX, (myCurve->nPol)*9*sizeof(double))) == NULL )
394 { printf(" ERROR: realloc failure rescX\n\n"); }
395
396 if((rescY = (double *)realloc((double *)rescY, (myCurve->nPol)*9*sizeof(double))) == NULL )
397 { printf(" ERROR: realloc failure rescY\n\n"); }
398
399 if((rescZ = (double *)realloc((double *)rescZ, (myCurve->nPol)*9*sizeof(double))) == NULL )
400 { printf(" ERROR: realloc failure rescZ\n\n"); }
401
402 double resPol[(myCurve->nPol)*9][3];
403 //calculate coordinate axis
404
405 vl_result(p1,p2,zAxis);//z-Axis of local coordinate system
406 int count=0;
407 do{
408 cPointTmp[0]=myCurve->cX[count];//Fist control point
409 cPointTmp[1]=myCurve->cY[count];
410 cPointTmp[2]=myCurve->cZ[count];
411 vl_result(p1,cPointTmp,p1cPoint);
412 count++;
413 vl_sprod(p1cPoint,zAxis,&scalProd);
414 length=vl_betrag(zAxis);
415 length = scalProd/(length*length);
416 vl_scal(&length,zAxis,tmpVec);
417 vl_add(p1,tmpVec,origin);//origin of local coordinate system
418 //cout<<"cPointTmp: "<<cPointTmp[0]<<" "<<cPointTmp[1]<<" "<<cPointTmp[2]<<endl;
419 //cout<<"go:"<<vl_betrag(p1cPoint)<<endl;
420 //cout<<"Origin: "<<origin[0]<<" "<<origin[1]<<" "<<origin[2]<<endl;
421 vl_result(origin,p2,zAxis);//z-Axis of local coordinate system
422 if(vl_betrag(zAxis) == 0) vl_result(p2,p1,zAxis);
423 vl_norm(zAxis,zAxis);//scale to length 1
424 vl_result(origin,cPointTmp,xAxis);//x-Axis of local coordinate system
425 }
426 while(vl_betrag(xAxis)<0.0000000001);//constant is used to avoid rounding errors
427 vl_norm(xAxis,xAxis);//scale to length 1
428 vl_prod(zAxis,xAxis,yAxis);//y-Axis of local coordinate system
429 //calculate rotaion matrix
430 for(i=0;i<3;i++){
431 M[i][0]=xAxis[i];//wirte x-Axis to fisrt column
432 M[i][1]=yAxis[i];//wirte y-Axis to second column
433 M[i][2]=zAxis[i];//wirte z-Axis to third column
434 //cout<<"( "<<xAxis[i]<<" "<<yAxis[i]<<" "<<zAxis[i]<<" )"<<endl;
435 }
436
437 for(i=0;i<(myCurve->nPol);i++)
438 {
439 vl_result(p1,p2,p1p2);
440 cPointTmp[0]=myCurve->cX[i];//control point Ci
441 cPointTmp[1]=myCurve->cY[i];
442 cPointTmp[2]=myCurve->cZ[i];
443 vl_result(p1,cPointTmp,p1cPoint);
444 vl_sprod(p1cPoint,p1p2,&scalProd);
445 length=vl_betrag(p1p2);
446 length = scalProd/(length*length);
447 vl_scal(&length,p1p2,tmpVec);
448 vl_add(p1,tmpVec,t);
449 vl_result(t,cPointTmp,tmpVec);
450 radius = vl_betrag(tmpVec);
451
452 //"scale circle points
453 for(l=0;l<9;l++){
454 for(m=0;m<3;m++){
455 circlePointScaled[l][m] = radius * circlePoint[l][m];
456 }
457 }
458
459 //calcutale control points
460 for(j=0;j<9;j++){//calculate control points of circle
461 vl_trans(M,t,circlePointScaled[j], resPol[i*9+j]);//calculate controll points of the rotated surface
462 rescX[i*9+j]=resPol[i*9+j][0];
463 rescY[i*9+j]=resPol[i*9+j][1];
464 rescZ[i*9+j]=resPol[i*9+j][2];
465 resWeight[i*9+j] = wCircle[j]*(myCurve->w[i]);//calculate weights
466 }
467 }
468 //calculate u-Knots
469 for(j=0;j<(myCurve->nKnt);j++){
470 resUKnot[j] = myCurve->k[j];
471 }
472 //calculate v-Knots
473 for(j=0;j<12;j++){
474 resVKnot[j] = kCircle[j];
475 }
476 //write result to struct
477 mySurf->uDeg = myCurve->deg;
478 mySurf->vDeg = 2;
479 mySurf->nUPol = myCurve->nPol;
480 mySurf->nVPol = 9;
481 mySurf->nUKnt = myCurve->nKnt;
482 mySurf->nVKnt = 12;
483 mySurf->uKnt = resUKnot;
484 mySurf->vKnt = resVKnot;
485
486 mySurf->weights = resWeight;
487 mySurf->cX =rescX;
488 mySurf->cY =rescY;
489 mySurf->cZ =rescZ;
490 }
491
wik(int k,int i,double * u,double x)492 inline double wik(int k,int i,double *u,double x)
493 {
494 if((u[i+k-1]-u[i])!=0){
495 return (x-u[i])/(u[i+k-1]-u[i]);
496 }
497 else
498 return 0;
499 }
500
deBoor(int k,int i,double * u,double x)501 inline double deBoor(int k, int i, double * u ,double x)
502 {
503 if(x>1)x=1;
504 if(k==1){
505 if(x >= u[i] && x <= u[i+1]) return 1;
506 else return 0;
507 }
508 else {
509 return wik(k,i,u,x)*deBoor(k-1,i,u,x)+(1-wik(k,i+1,u,x))*deBoor(k-1,i+1,u,x);
510 }
511 }
512
calculateBSpline(double * pnt,BSplineCurve * myCurve,double u)513 inline void calculateBSpline(double * pnt, BSplineCurve * myCurve,double u)
514 {
515 double sumX=0,sumY=0,sumZ=0,N,tmp=0;
516 int i,j;
517 for(j=0;j<(myCurve->nPol);j++){
518 tmp+=(myCurve->w[j])*deBoor(myCurve->deg+1,j,(myCurve->k),u);
519 }
520 if(tmp!=0){
521 for(i=0;i<(myCurve->nPol);i++)
522 {
523 N=(myCurve->w[i])*deBoor(myCurve->deg+1,i,(myCurve->k),u)/tmp;
524 sumX+= myCurve->cX[i]*(N);
525 sumY+= myCurve->cY[i]*(N);
526 sumZ+= myCurve->cZ[i]*(N);
527 }
528 }
529 else{
530 sumX = sumY = sumZ = 0;
531 }
532 pnt[0]=sumX;
533 pnt[1]=sumY;
534 pnt[2]=sumZ;
535 }
536
537 /**************************************************/
538 /*** fitting algorithm for B-Splines ***/
539 /*** progressive-iterative approximation (PIA) ***/
540 /**************************************************/
piaFitting(double pCloud[][3],int nPnt,BSplineCurve * fitCurve,int deg,double tolerance)541 void piaFitting(double pCloud [][3],int nPnt,BSplineCurve * fitCurve, int deg, double tolerance)
542 {
543 if(nPnt<2) cout<<"WARNING: too few Points"<<endl;
544 fitCurve->k = NULL;fitCurve->cX = NULL;fitCurve->cY = NULL;fitCurve->cZ = NULL;fitCurve->w = NULL;
545 int nCPnt = nPnt + 2, i;
546 if((fitCurve->k = (double *)realloc((double *)fitCurve->k,(nPnt+2*deg)*sizeof(double))) == NULL )
547 { printf(" ERROR: realloc failure\n\n"); }
548 if((fitCurve->cX = (double *)realloc((double *)fitCurve->cX,nCPnt*sizeof(double))) == NULL )
549 { printf(" ERROR: realloc failure\n\n"); }
550 if((fitCurve->cY = (double *)realloc((double *)fitCurve->cY,nCPnt*sizeof(double))) == NULL )
551 { printf(" ERROR: realloc failure\n\n"); }
552 if((fitCurve->cZ = (double *)realloc((double *)fitCurve->cZ,nCPnt*sizeof(double))) == NULL )
553 { printf(" ERROR: realloc failure\n\n"); }
554 if((fitCurve->w = (double *)realloc((double *)fitCurve->w,nCPnt*sizeof(double))) == NULL )
555 { printf(" ERROR: realloc failure\n\n"); }
556
557
558 fitCurve->deg = deg;//cubic Bspline
559 fitCurve->nPol = nCPnt;
560 fitCurve->nKnt = nPnt+2*deg;
561 fitCurve->deg = deg;
562 //Define curve as clamped:
563 //Start knodes have multiplicity deg + 1
564 //End knodes have multiplicity deg + 1
565 for(i=0;i<=deg;i++)
566 {
567 fitCurve->k[i]=0;
568 fitCurve->k[fitCurve->nKnt-(i+1)]=1;
569 }
570
571 double dist=0, p1p2[3], p1[3],p2[3];
572 double du,dControll[3];
573
574 //sum up distances between points
575 for(i=1;i<nPnt;i++){
576 vl_result(pCloud[i],pCloud[i-1],p1p2);
577 dist+=vl_betrag(p1p2);
578 }
579
580 //set pCloud points as control points
581 fitCurve->cX[0]=pCloud[0][0];
582 fitCurve->cY[0]=pCloud[0][1];
583 fitCurve->cZ[0]=pCloud[0][2];
584 fitCurve->w[0]=1;
585 fitCurve->cX[nCPnt-1]=pCloud[nPnt-1][0];
586 fitCurve->cY[nCPnt-1]=pCloud[nPnt-1][1];
587 fitCurve->cZ[nCPnt-1]=pCloud[nPnt-1][2];
588 fitCurve->w[nCPnt-1]=1;
589
590 for(i=1;i<=nPnt;i++){
591 fitCurve->cX[i]=pCloud[i-1][0];
592 fitCurve->cY[i]=pCloud[i-1][1];
593 fitCurve->cZ[i]=pCloud[i-1][2];
594 fitCurve->w[i]=1;
595 }
596
597 //calculate initial knot vector via chord-length parametrization
598 for(i=1;i<nPnt;i++){
599 p1[0]=pCloud[i-1][0];
600 p1[1]=pCloud[i-1][1];
601 p1[2]=pCloud[i-1][2];
602 p2[0]=pCloud[i][0];
603 p2[1]=pCloud[i][1];
604 p2[2]=pCloud[i][2];
605 vl_result(p1,p2,p1p2);
606 du=(vl_betrag(p1p2))/dist;
607 fitCurve->k[i+deg]=fitCurve->k[i+deg-1]+du;
608 }
609
610 double error=10, errorOld=0;
611 int c, iter = 1;
612 //interation loop for PIA fitting
613 do{
614 //iterate over all inner knots
615 c=0;
616 errorOld = error;
617 error=0;
618 for(i=(deg);i<((fitCurve->nKnt)-(deg));i++)
619 {
620 calculateBSpline(dControll,fitCurve,fitCurve->k[i]);//calculate point for parameter with current control points
621 p1[0]=pCloud[c][0];
622 p1[1]=pCloud[c][1];
623 p1[2]=pCloud[c][2];
624 vl_result(dControll,p1,p2);//get the difference between calculated spline point and control point
625 error+=vl_betrag(p2);//get the difference to current control point
626 p1[0]=fitCurve->cX[c+1];
627 p1[1]=fitCurve->cY[c+1];
628 p1[2]=fitCurve->cZ[c+1];
629 vl_add(p1,p2,dControll);
630 fitCurve->cX[c+1]=dControll[0];
631 fitCurve->cY[c+1]=dControll[1];
632 fitCurve->cZ[c+1]=dControll[2];
633 c++;
634 }
635 #if PLOT == 1
636 cout<<"Error: "<<error<<" after "<<iter<<"-iterations"<<endl;
637 #endif
638 iter++;
639 }while(sqrt((error-errorOld)*(error-errorOld))>dist*tolerance);
640 cout<<"BSpline fitting ... done"<<endl;
641 cout<<"Error: "<<error<<" after "<<iter<<"-iterations"<<endl;
642 }
643
644
645 #ifdef __cplusplus
646 extern "C" {
647 #endif
648 extern SumGeo anzGeo[1];
649 #ifdef __cplusplus
650 }
651 #endif
652
653 // free cpp array (to be called from c code to free cpp objects)
cppfreearray(double * array)654 void cppfreearray(double *array)
655 {
656 delete[] array;
657 }
658
659 //int createBlendedNurbs(int nr, Points **pntpntr, Lines *line, Lcmb *lcmb, Gsur *surf )
createBlendedNurbs(int nr)660 int createBlendedNurbs(int nr)
661 {
662 int i,j,k,l,cl,n,p, nip, flag, np[4], index;
663 static double *px[4], *py[4], *pz[4];
664 char name[MAX_LINE_LENGTH], buffer[MAX_LINE_LENGTH];
665 char sname[MAX_LINE_LENGTH];
666
667 int e;
668 double p0[3],p1[3],p0p1[3],ep0p1[3],vl,dvl,p0pn[3],pn[3];
669 extern SumAsci sumAsci[1];
670
671 // printf("surf:%s\n", surf[nr].name);
672 if(surf[nr].nl!=4)
673 {
674 printf("ERROR: nr of edges must be 4 but is %d for surf:%s to create a blended nurbs\n\n",surf[nr].nl,surf[nr].name);
675 return(-1);
676 }
677 /* create an array of points for each edge */
678 for(j=0; j<surf[nr].nl; j++)
679 {
680 p=0;
681 nip=0;
682 if(surf[nr].typ[j]=='l')
683 {
684 // printf("line:%s\n", line[surf[nr].l[j]].name);
685 l=surf[nr].l[j];
686 nip+=line[l].nip;
687
688 if( (px[j] = (double *)realloc( (double *)px[j], (nip/3+1)*sizeof(double ) ))==NULL )
689 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
690 return(-1); }
691 if( (py[j] = (double *)realloc( (double *)py[j], (nip/3+1)*sizeof(double ) ))==NULL )
692 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
693 return(-1); }
694 if( (pz[j] = (double *)realloc( (double *)pz[j], (nip/3+1)*sizeof(double ) ))==NULL )
695 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
696 return(-1); }
697
698 if(surf[nr].o[j]=='+')
699 {
700 n=0; flag=line[l].nip;
701 do
702 {
703 px[j][p]=line[l].ip[n++];
704 py[j][p]=line[l].ip[n++];
705 pz[j][p]=line[l].ip[n++];
706 p++;
707 }while(n<flag);
708 }
709 else
710 {
711 n=line[l].nip; flag=0;
712 while(n>flag)
713 {
714 pz[j][p]=line[l].ip[--n];
715 py[j][p]=line[l].ip[--n];
716 px[j][p]=line[l].ip[--n];
717 p++;
718 }
719 }
720 }
721 else
722 {
723 // printf("lcmb:%s\n", lcmb[surf[nr].l[j]].name);
724 if(surf[nr].o[j]=='+')
725 {
726 for(cl=0; cl<lcmb[surf[nr].l[j]].nl; cl++)
727 {
728 l=lcmb[surf[nr].l[j]].l[cl];
729
730 nip+=line[l].nip;
731 if( (px[j] = (double *)realloc( (double *)px[j], (nip/3+1)*sizeof(double ) ))==NULL )
732 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
733 return(-1); }
734 if( (py[j] = (double *)realloc( (double *)py[j], (nip/3+1)*sizeof(double ) ))==NULL )
735 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
736 return(-1); }
737 if( (pz[j] = (double *)realloc( (double *)pz[j], (nip/3+1)*sizeof(double ) ))==NULL )
738 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
739 return(-1); }
740
741 if(lcmb[surf[nr].l[j]].o[cl]=='-') flag=-1;
742 else flag=1;
743
744 if(flag==1)
745 {
746
747 if(cl==0) { n=0; flag=line[l].nip-3; }
748 else if(cl==lcmb[surf[nr].l[j]].nl-1) { n=0; flag=line[l].nip; }
749 else { n=0; flag=line[l].nip-3; }
750
751 do
752 {
753 px[j][p]=line[l].ip[n++];
754 py[j][p]=line[l].ip[n++];
755 pz[j][p]=line[l].ip[n++];
756 p++;
757 }while(n<flag);
758 }
759 else
760 {
761
762 if(cl==0) { n=line[l].nip; flag=3; }
763 else if(cl==lcmb[surf[nr].l[j]].nl-1) { n=line[l].nip; flag=0; }
764 else { n=line[l].nip; flag=3; }
765
766 while(n>flag)
767 {
768 pz[j][p]=line[l].ip[--n];
769 py[j][p]=line[l].ip[--n];
770 px[j][p]=line[l].ip[--n];
771 p++;
772 }
773 }
774 }
775 }
776 else
777 {
778 for(cl=lcmb[surf[nr].l[j]].nl-1; cl>=0; cl--)
779 {
780 l=lcmb[surf[nr].l[j]].l[cl];
781
782 nip+=line[l].nip;
783 if( (px[j] = (double *)realloc( (double *)px[j], (nip/3+1)*sizeof(double ) ))==NULL )
784 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
785 return(-1); }
786 if( (py[j] = (double *)realloc( (double *)py[j], (nip/3+1)*sizeof(double ) ))==NULL )
787 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
788 return(-1); }
789 if( (pz[j] = (double *)realloc( (double *)pz[j], (nip/3+1)*sizeof(double ) ))==NULL )
790 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
791 return(-1); }
792
793 if(lcmb[surf[nr].l[j]].o[cl]=='-') flag=1;
794 else flag=-1;
795
796 if(flag==1)
797 {
798
799 if(cl==0) { n=0; flag=line[l].nip; }
800 else if(cl==lcmb[surf[nr].l[j]].nl-1) { n=0; flag=line[l].nip-3; }
801 else { n=0; flag=line[l].nip-3; }
802
803 do
804 {
805 px[j][p]=line[l].ip[n++];
806 py[j][p]=line[l].ip[n++];
807 pz[j][p]=line[l].ip[n++];
808 p++;
809 }while(n<flag);
810 }
811 else
812 {
813
814 if(cl==0) { n=line[l].nip; flag=0; }
815 else if(cl==lcmb[surf[nr].l[j]].nl-1) { n=line[l].nip; flag=3; }
816 else { n=line[l].nip; flag=3; }
817
818 while(n>flag)
819 {
820 pz[j][p]=line[l].ip[--n];
821 py[j][p]=line[l].ip[--n];
822 px[j][p]=line[l].ip[--n];
823 p++;
824 }
825 }
826 }
827 }
828 }
829 np[j]=p;
830 }
831
832 /* each edge must have at least 4 points */
833 for (e=0; e<4; e++) if(np[e]<4)
834 {
835 if( (px[e] = (double *)realloc( (double *)px[e], (5)*sizeof(double ) ))==NULL )
836 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
837 return(-1); }
838 if( (py[e] = (double *)realloc( (double *)py[e], (5)*sizeof(double ) ))==NULL )
839 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
840 return(-1); }
841 if( (pz[e] = (double *)realloc( (double *)pz[e], (5)*sizeof(double ) ))==NULL )
842 { printf(" ERROR: realloc failure in createBlendedNurbs()\n\n");
843 return(-1); }
844 p0[0]=px[e][np[e]-2];
845 p0[1]=py[e][np[e]-2];
846 p0[2]=pz[e][np[e]-2];
847 p1[0]=px[e][np[e]-1];
848 p1[1]=py[e][np[e]-1];
849 p1[2]=pz[e][np[e]-1];
850 vl_result(p0,p1,p0p1);
851 vl_norm(p0p1,ep0p1);
852 vl=vl_betrag(p0p1);
853 nip=4-np[e];
854 dvl=vl/(nip+1);
855 vl=0.;
856 for(i=0; i<=nip; i++)
857 {
858 vl+=dvl;
859 vl_scal(&vl,ep0p1, p0pn);
860 vl_add(p0,p0pn,pn);
861 // printf("PNT ! %f %f %f\n", pn[0],pn[1],pn[2]);
862 px[e][np[e]-1+i]=pn[0];
863 py[e][np[e]-1+i]=pn[1];
864 pz[e][np[e]-1+i]=pn[2];
865 }
866 np[e]+=nip;
867 }
868 /*
869 for(j=0; j<4; j++)
870 {
871 printf("# pnts:%d\n", np[j]);
872 for ( i= 0; i < np[j] ; i ++ ) printf(" PNT ! %f %f %f\n", px[j][i], py[j][i], pz[j][i]);
873 }
874 */
875
876 // Create array of points to be passed to snlCurve constructor.
877 // printf(" Create array of points\n");
878 snlPoint* curvePointsEdge1 = new snlPoint [ np[0] ];
879 snlPoint* curvePointsEdge2 = new snlPoint [ np[1] ];
880 snlPoint* curvePointsEdge3 = new snlPoint [ np[2] ];
881 snlPoint* curvePointsEdge4 = new snlPoint [ np[3] ];
882
883 j=0; for ( i= 0; i < np[j] ; i ++ ) curvePointsEdge1 [ i ].components ( px[j][i], py[j][i], pz[j][i] );
884 j=1; for ( i= 0; i < np[j] ; i ++ ) curvePointsEdge2 [ i ].components ( px[j][i], py[j][i], pz[j][i] );
885 j=2; for ( i= np[j]-1; i >= 0 ; i -- ) curvePointsEdge3 [ np[j]-1-i ].components ( px[j][i], py[j][i], pz[j][i] );
886 j=3; for ( i= np[j]-1; i >= 0 ; i -- ) curvePointsEdge4 [ np[j]-1-i ].components ( px[j][i], py[j][i], pz[j][i] );
887
888 // Create curves. All curves must have the same degree. "degree" is an integer.
889 int degree=3;
890 snlCurve* curveEdgeU1 = new snlCurve ( curvePointsEdge1,
891 np[0],
892 snlCurve::SNL_GLOBAL_INTERP_CENTRIFUGAL,
893 degree );
894 snlCurve* curveEdgeV2 = new snlCurve ( curvePointsEdge2,
895 np[1],
896 snlCurve::SNL_GLOBAL_INTERP_CENTRIFUGAL,
897 degree );
898 snlCurve* curveEdgeU2 = new snlCurve ( curvePointsEdge3,
899 np[2],
900 snlCurve::SNL_GLOBAL_INTERP_CENTRIFUGAL,
901 degree );
902 snlCurve* curveEdgeV1 = new snlCurve ( curvePointsEdge4,
903 np[3],
904 snlCurve::SNL_GLOBAL_INTERP_CENTRIFUGAL,
905 degree );
906
907 // Create bilinear Coons patch. The orientation of the curves given to this
908 // function is very important. There are two curves in the U direction
909 // and two in the V direction. They should be oriented as follows:
910 //
911 // V1
912 // -------->
913 // | |
914 // U1 | | U2
915 // | |
916 // V V
917 // -------->
918 // V2
919 //
920 // The arrow heads are the end of the curve. If you want to reverse
921 // the direction of the curve call it's reverseEvalDirection() function.
922
923 // You will have to correspond curveEdge1 etc to one of the directions in
924 // the following constructor.
925
926 // printf(" Create bilinear Coons patch\n");
927 snlSurface * surface = new snlSurface ( curveEdgeU1,curveEdgeU2,curveEdgeV1,curveEdgeV2 );
928
929 // You now have a new NURBS Coons Patch.
930
931 // Clean up allocated memory.
932
933 delete[] curvePointsEdge1;
934 delete[] curvePointsEdge2;
935 delete[] curvePointsEdge3;
936 delete[] curvePointsEdge4;
937
938 delete curveEdgeU1;
939 delete curveEdgeU2;
940 delete curveEdgeV1;
941 delete curveEdgeV2;
942
943
944 // (5) create the cgx nurbs
945 // printf(" create blended nurbs\n");
946 buffer[0]='S';
947 buffer[1]='\0';
948
949 sem_wait(&sem_g);
950
951 getNewName( name, buffer );
952 for (i=0; i<MAX_LINE_LENGTH; i++) sname[i]=name[i];
953
954 nr=anzGeo->nurs++;
955 if(printFlag) printf ("store NURS Nr:%d Name:%s\n", anzGeo->nurs, name);
956
957 if ((nurbs = (Nurbs *)realloc( (Nurbs *)nurbs, (anzGeo->nurs)*sizeof(Nurbs)) ) == NULL )
958 { printf("\n\nERROR: realloc failure in Nurs, nurbs:%s not installed\n\n", name); anzGeo->nurs--; }
959
960 // hashNurs was deactivated but its needed. Could be moved into the calling function
961 hashNurs( sumAsci, name, nr );
962
963 k=strlen(sname);
964 if((nurbs[nr].name= (char *)malloc((k+1)*sizeof(char))) == NULL )
965 { printf("ERROR: malloc failed\n\n" ); }
966 for (i=0; i<=k; i++) nurbs[nr].name[i]=name[i];
967
968 /* the nurbs fits exactly, no trimming is necessary */
969 nurbs[nr].trimFlag=0;
970 nurbs[nr].u_exp=surface->degreeU();
971 nurbs[nr].v_exp=surface->degreeV();
972 nurbs[nr].u_npnt=surface->sizeU();
973 nurbs[nr].v_npnt=surface->sizeV();
974 nurbs[nr].u_nknt=1+nurbs[nr].u_npnt+ nurbs[nr].u_exp;
975 nurbs[nr].v_nknt=1+nurbs[nr].v_npnt+ nurbs[nr].v_exp;
976 nurbs[nr].v_stride=4;
977 nurbs[nr].u_stride=4* nurbs[nr].v_npnt;
978
979 nurbs[nr].uv=NULL;
980 nurbs[nr].xyz=NULL;
981 nurbs[nr].np=NULL;
982 nurbs[nr].nc=NULL;
983 nurbs[nr].vmax=NULL;
984 nurbs[nr].umax=NULL;
985 nurbs[nr].vstep=NULL;
986 nurbs[nr].ustep=NULL;
987 nurbs[nr].sum_ambiguousPnts=NULL;
988 nurbs[nr].uvflipped=NULL;
989
990 nurbs[nr].ctlarray=NULL;
991
992 nurbs[nr].patches=0;
993 nurbs[nr].endFlag=1;
994 nurbs[nr].type=GL_MAP2_VERTEX_4;
995 nurbs[nr].Nurb = (GLUnurbsObj *)gluNewNurbsRenderer();
996
997 // knots
998 if ( (nurbs[nr].uknt = (float *)malloc((nurbs[nr].u_nknt+1) * sizeof(float))) == NULL )
999 printf("\n\n ERROR: malloc failed uknt\n\n");
1000 if ( (nurbs[nr].vknt = (float *)malloc((nurbs[nr].v_nknt+1) * sizeof(float))) == NULL )
1001 printf("\n\n ERROR: malloc failed vknt\n\n");
1002 for(i=0; i<nurbs[nr].u_nknt; i++) { nurbs[nr].uknt[i]=surface->knotsU()[i]; }
1003 for(i=0; i<nurbs[nr].v_nknt; i++) { nurbs[nr].vknt[i]=surface->knotsV()[i]; }
1004
1005 // discrete controll points and weights
1006 if ( (nurbs[nr].ctlpnt = (int **)malloc( (nurbs[nr].u_npnt+1) * sizeof(int *))) == NULL )
1007 printf("\n\n ERROR: malloc failed ctlpnt\n\n");
1008 for (i=0; i<nurbs[nr].u_npnt; i++)
1009 {
1010 if ( (nurbs[nr].ctlpnt[i] = (int *)malloc( (nurbs[nr].v_npnt+1) * sizeof( int ))) == NULL )
1011 printf("\n\n ERROR: malloc failed ctlpnt[i]\n\n");
1012 }
1013 if ( (nurbs[nr].weight = (GLfloat **)malloc( (nurbs[nr].u_npnt+1) * sizeof(GLfloat *))) == NULL )
1014 printf("\n\n ERROR: malloc failed weight\n\n");
1015 for (i=0; i<nurbs[nr].u_npnt; i++)
1016 {
1017 if ( (nurbs[nr].weight[i] = (GLfloat *)malloc( (nurbs[nr].v_npnt+1) * sizeof(GLfloat))) == NULL )
1018 printf("\n\n ERROR: malloc failed weight[i]\n\n");
1019 }
1020
1021 index=0;
1022 for (i=0; i<nurbs[nr].u_npnt; i++)
1023 {
1024 for (j=0; j<nurbs[nr].v_npnt; j++)
1025 {
1026 buffer[0]='p';
1027 buffer[1]='\0';
1028 getNewName( name, buffer );
1029 nurbs[nr].ctlpnt[i][j] = pnt( name, surface->controlPoints()[ index ].x(), surface->controlPoints()[ index ].y(), surface->controlPoints()[ index ].z(), 0);
1030 nurbs[nr].weight[i][j] = surface->controlPoints()[ index ].w();
1031 index++;
1032 }
1033 }
1034
1035 #if TEST
1036 printf("out: degUV %d %d sizeUV %u %u\n", surface->degreeU(), surface->degreeV(), surface->sizeU(), surface->sizeV());
1037 printf("knotu:%d\n", nurbs[nr].u_nknt);
1038 printf("knotv:%d\n", nurbs[nr].v_nknt);
1039 for(i=0; i<nurbs[nr].u_nknt; i++) { printf("knu:%f\n", nurbs[nr].uknt[i]); }
1040 for(i=0; i<nurbs[nr].v_nknt; i++) { printf("knv:%f\n", nurbs[nr].vknt[i]); }
1041 index=0;
1042 for (i=0; i<nurbs[nr].u_npnt; i++)
1043 {
1044 for (j=0; j<nurbs[nr].v_npnt; j++)
1045 {
1046 printf("pnt: %f %f %f %f \n", surface->controlPoints()[ index ].x(),surface->controlPoints()[ index ].y(),surface->controlPoints()[ index ].z(),surface->controlPoints()[ index ].w());
1047 index++;
1048 }
1049 }
1050 #endif
1051 sem_post(&sem_g);
1052
1053 delete surface; // Release surface object.
1054 return(nr);
1055 }
1056
1057
1058
1059 // uses double reduceDegree( int dir, unsigned numDeg, double tolerance )
1060 // returns the error during reduction. "dir" is the direction you wish to degree
1061 // reduce, it is one of the constants:
1062 //
1063 // snlSurface::SNL_U_DIR
1064 // snlSurface::SNL_V_DIR
1065 //
1066 // for the U and V directions respectively. "numDeg" is the number of degrees you w
1067 // ish the surface to be reduced by and "tolerance" is the maximum error allowed du
1068 // ring processing. If tolerance is exceeded then the surface is restored to origin
1069 // al. If you don't want a tolerance and will take any error then use a value of 0.
1070 // This function will increase the size of the data describing the surface consider
1071 // ably. I am yet to implement a function that remedies that, but in the mean time
1072 // it will work just the same.
1073 //
1074 // WARNING: The size will change in certain situations (180deg sector of a cylinder)
1075 // in this case cgx is not able to render the new nurbs because the master surf
1076 // might be bigger than the remaining nurbs. The program will crash.
1077 //
repairNurbs(int nr,int deg,int dir)1078 double repairNurbs( int nr, int deg, int dir )
1079 {
1080 int i,j,index=0;
1081
1082 // warning sem_wait in this function not allowed since it is called inside a sem block
1083 //sem_wait(&sem_g);
1084
1085 int numberOfControlPoints =nurbs[nr].u_npnt*nurbs[nr].v_npnt;
1086 snlCtrlPoint* controlPoints = new snlCtrlPoint [ numberOfControlPoints ];
1087 double error;
1088
1089
1090 // (1) Create an array of control points from your cgx control points:
1091 for (i=0; i<nurbs[nr].u_npnt; i++)
1092 {
1093 for (j=0; j<nurbs[nr].v_npnt; j++)
1094 {
1095 controlPoints [ index++ ].components (
1096 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0],
1097 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1],
1098 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2],
1099 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3] );
1100 }
1101 }
1102
1103 // (2) Create an array of knots for each knot vector in your cgx surface:
1104 knot* knotsU = new knot [ nurbs[nr].u_nknt ];
1105 knot* knotsV = new knot [ nurbs[nr].v_nknt ];
1106
1107 for(i=0; i<nurbs[nr].u_nknt; i++) knotsU[i]=nurbs[nr].uknt[i];
1108 for(i=0; i<nurbs[nr].v_nknt; i++) knotsV[i]=nurbs[nr].vknt[i];
1109
1110 // (3) Create an snlSurface:
1111 // You need to supply degreeU, degreeV, sizeU and sizeV.
1112 snlSurface* surface = new snlSurface ( nurbs[nr].u_exp, nurbs[nr].v_exp, nurbs[nr].u_npnt, nurbs[nr].v_npnt, controlPoints, knotsU, knotsV );
1113
1114 #if TEST
1115 // print the nurbs-parameter before reduction
1116 printf("in: degUV %d %d sizeUV %u %u\n", surface->degreeU(), surface->degreeV(), surface->sizeU(), surface->sizeV());
1117 for(i=0; i<nurbs[nr].u_nknt; i++) { printf("knu:%f\n", nurbs[nr].uknt[i]); }
1118 for(i=0; i<nurbs[nr].v_nknt; i++) { printf("knv:%f\n", nurbs[nr].vknt[i]); }
1119 index=0;
1120 for (i=0; i<nurbs[nr].u_npnt; i++)
1121 {
1122 for (j=0; j<nurbs[nr].v_npnt; j++)
1123 {
1124 printf("pnt: %f %f %f %f \n", controlPoints [ index ].x(),controlPoints [ index ].y(),controlPoints [ index ].z(),controlPoints [ index ].w());
1125 index++;
1126 }
1127 }
1128 #endif
1129
1130 // (4) reduce degree
1131 if(dir==0) error=surface -> reduceDegree ( snlSurface::SNL_U_DIR, deg, 0.0 );
1132 else error=surface -> reduceDegree ( snlSurface::SNL_V_DIR, deg, 0.0 );
1133
1134 // (5) change the cgx nurbs
1135 nurbs[nr].u_exp=surface->degreeU();
1136 nurbs[nr].v_exp=surface->degreeV();
1137 nurbs[nr].u_npnt=surface->sizeU();
1138 nurbs[nr].v_npnt=surface->sizeV();
1139 nurbs[nr].u_nknt=1+nurbs[nr].u_npnt+ nurbs[nr].u_exp;
1140 nurbs[nr].v_nknt=1+nurbs[nr].v_npnt+ nurbs[nr].v_exp;
1141 nurbs[nr].u_stride=4* nurbs[nr].v_npnt;
1142
1143 // knots
1144 if ( (nurbs[nr].uknt = (float *)realloc( (float *)nurbs[nr].uknt, (nurbs[nr].u_nknt+1) * sizeof(float))) == NULL )
1145 printf("\n\n ERROR: realloc failed uknt\n\n");
1146 if ( (nurbs[nr].vknt = (float *)realloc( (float *)nurbs[nr].vknt, (nurbs[nr].v_nknt+1) * sizeof(float))) == NULL )
1147 printf("\n\n ERROR: realloc failed vknt\n\n");
1148 for(i=0; i<nurbs[nr].u_nknt; i++) { nurbs[nr].uknt[i]=surface->knotsU()[i]; }
1149 for(i=0; i<nurbs[nr].v_nknt; i++) { nurbs[nr].vknt[i]=surface->knotsV()[i]; }
1150
1151 // controll-array, discrete controll points and weights
1152 if( (nurbs[nr].ctlarray = (GLfloat *)realloc( (GLfloat *)nurbs[nr].ctlarray, (nurbs[nr].u_npnt*nurbs[nr].v_npnt*nurbs[nr].v_stride+5)*sizeof(GLfloat) )) == NULL )
1153 { printf(" ERROR: realloc failure in repairNurbs(), nurbs:%s can not be shaped\n\n", nurbs[nr].name);
1154 sem_post(&sem_g); return(-1); }
1155 index=0;
1156 for (i=0; i<nurbs[nr].u_npnt; i++)
1157 {
1158 for (j=0; j<nurbs[nr].v_npnt; j++)
1159 {
1160 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]=surface->controlPoints()[ index ].x();
1161 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]=surface->controlPoints()[ index ].y();
1162 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]=surface->controlPoints()[ index ].z();
1163 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]=surface->controlPoints()[ index ].w();
1164 index++;
1165 }
1166 }
1167
1168 #if TEST
1169 // print
1170 printf("out: degUV %d %d sizeUV %u %u\n", surface->degreeU(), surface->degreeV(), surface->sizeU(), surface->sizeV());
1171 printf("knotu:%d\n", nurbs[nr].u_nknt);
1172 printf("knotv:%d\n", nurbs[nr].v_nknt);
1173 for(i=0; i<nurbs[nr].u_nknt; i++) { printf("knu:%f\n", nurbs[nr].uknt[i]); }
1174 for(i=0; i<nurbs[nr].v_nknt; i++) { printf("knv:%f\n", nurbs[nr].vknt[i]); }
1175 index=0;
1176 for (i=0; i<nurbs[nr].u_npnt; i++)
1177 {
1178 for (j=0; j<nurbs[nr].v_npnt; j++)
1179 {
1180 printf("pnt: %f %f %f %f \n", surface->controlPoints()[ index ].x(),surface->controlPoints()[ index ].y(),surface->controlPoints()[ index ].z(),surface->controlPoints()[ index ].w());
1181 index++;
1182 }
1183 }
1184 #endif
1185 //sem_post(&sem_g);
1186
1187 delete surface; // Release surface object.
1188
1189 // The rest of the objects are deleted by libSNL.
1190
1191 return(error);
1192 }
1193
1194
evalNurbs(int nr,int sum_p,double * pnt_u,double * pnt_v,Points * pnt)1195 int evalNurbs( int nr, int sum_p, double *pnt_u, double *pnt_v, Points *pnt)
1196 {
1197 int i,j,index=0;
1198 snlPoint snlPnt;
1199
1200 sem_wait(&sem_g);
1201
1202 int numberOfControlPoints =nurbs[nr].u_npnt*nurbs[nr].v_npnt;
1203 snlCtrlPoint* controlPoints = new snlCtrlPoint [ numberOfControlPoints ];
1204
1205
1206 // (1) Create an array of control points from your cgx control points:
1207 for (i=0; i<nurbs[nr].u_npnt; i++)
1208 {
1209 for (j=0; j<nurbs[nr].v_npnt; j++)
1210 {
1211 controlPoints [ index++ ].components (
1212 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0],
1213 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1],
1214 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2],
1215 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3] );
1216 }
1217 }
1218
1219 // (2) Create an array of knots for each knot vector in your cgx surface:
1220 knot* knotsU = new knot [ nurbs[nr].u_nknt ];
1221 knot* knotsV = new knot [ nurbs[nr].v_nknt ];
1222
1223 for(i=0; i<nurbs[nr].u_nknt; i++) knotsU[i]=(double)nurbs[nr].uknt[i];
1224 for(i=0; i<nurbs[nr].v_nknt; i++) knotsV[i]=(double)nurbs[nr].vknt[i];
1225
1226 // (3) Create an snlSurface:
1227 // You need to supply degreeU, degreeV, sizeU and sizeV.
1228 snlSurface* surface = new snlSurface ( nurbs[nr].u_exp, nurbs[nr].v_exp, nurbs[nr].u_npnt, nurbs[nr].v_npnt, controlPoints, knotsU, knotsV );
1229
1230 // (4) get coordinates
1231 for(i=0; i<sum_p; i++)
1232 {
1233 // negative values lead to a segfault
1234 if(pnt_u[i]<0.) pnt_u[i]=0.;
1235 if(pnt_v[i]<0.) pnt_v[i]=0.;
1236 snlPnt= surface -> eval ( pnt_u[i], pnt_v[i] );
1237 pnt[i].px=snlPnt.x();
1238 pnt[i].py=snlPnt.y();
1239 pnt[i].pz=snlPnt.z();
1240 }
1241 sem_post(&sem_g);
1242
1243 delete surface; // Release surface object.
1244
1245 // The rest of the objects are deleted by libSNL.
1246 return(1);
1247 }
1248
1249
evalNurbsWithNormalVector(int nr,int sum_p,double * pnt_u,double * pnt_v,Points * pnt,Points * nv)1250 int evalNurbsWithNormalVector( int nr, int sum_p, double *pnt_u, double *pnt_v, Points *pnt, Points *nv)
1251 {
1252 int i,j,index=0;
1253 snlPoint snlPnt;
1254
1255 sem_wait(&sem_g);
1256
1257 int numberOfControlPoints =nurbs[nr].u_npnt*nurbs[nr].v_npnt;
1258 snlCtrlPoint* controlPoints = new snlCtrlPoint [ numberOfControlPoints ];
1259
1260
1261 // (1) Create an array of control points from your cgx control points:
1262 for (i=0; i<nurbs[nr].u_npnt; i++)
1263 {
1264 for (j=0; j<nurbs[nr].v_npnt; j++)
1265 {
1266 controlPoints [ index++ ].components (
1267 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0],
1268 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1],
1269 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2],
1270 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3] );
1271 }
1272 }
1273
1274 // (2) Create an array of knots for each knot vector in your cgx surface:
1275 knot* knotsU = new knot [ nurbs[nr].u_nknt ];
1276 knot* knotsV = new knot [ nurbs[nr].v_nknt ];
1277
1278 for(i=0; i<nurbs[nr].u_nknt; i++) knotsU[i]=nurbs[nr].uknt[i];
1279 for(i=0; i<nurbs[nr].v_nknt; i++) knotsV[i]=nurbs[nr].vknt[i];
1280
1281 // (3) Create an snlSurface:
1282 // You need to supply degreeU, degreeV, sizeU and sizeV.
1283 snlSurface* surface = new snlSurface ( nurbs[nr].u_exp, nurbs[nr].v_exp, nurbs[nr].u_npnt, nurbs[nr].v_npnt, controlPoints, knotsU, knotsV );
1284
1285 snlVector refNorm;
1286 // (4) get coordinates
1287 for(i=0; i<sum_p; i++)
1288 {
1289 // negative values lead to a segfault
1290 if(pnt_u[i]<0.) pnt_u[i]=0.;
1291 if(pnt_v[i]<0.) pnt_v[i]=0.;
1292 snlPnt= surface -> eval ( pnt_u[i], pnt_v[i] );
1293 refNorm = surface -> calcNormal ( pnt_u[i], pnt_v[i] );
1294 nv[i].px=refNorm.x();
1295 nv[i].py=refNorm.y();
1296 nv[i].pz=refNorm.z();
1297 pnt[i].px=snlPnt.x();
1298 pnt[i].py=snlPnt.y();
1299 pnt[i].pz=snlPnt.z();
1300 }
1301 sem_post(&sem_g);
1302
1303 delete surface; // Release surface object.
1304
1305 // The rest of the objects are deleted by libSNL.
1306 return(1);
1307 }
1308
1309
1310
projSurfToNurbs(int nr,Gsur * surf,int snr,Nodes ** node)1311 void projSurfToNurbs( int nr, Gsur *surf, int snr, Nodes **node )
1312 {
1313 int i,j,n,k,index=0;
1314 int sum_p=0, sum_inverted=0, returnedPntsPerPnt;
1315
1316 double convergTol, normTol;
1317 int maxPass;
1318 snlSurfLocn* inverted;
1319 int *isort=NULL;
1320
1321 sem_wait(&sem_g);
1322
1323 int numberOfControlPoints =nurbs[nr].u_npnt*nurbs[nr].v_npnt;
1324 snlCtrlPoint* controlPoints = new snlCtrlPoint [ numberOfControlPoints ];
1325
1326 normTol=COS_TOL;
1327 convergTol=ITER_TOL;
1328 maxPass=MAX_PASS;
1329
1330 if(node==0) sum_p=surf[snr].npgn;
1331 else sum_p=surf[snr].nn;
1332
1333 if(sum_p<=0) { sem_post(&sem_g); return; }
1334
1335 // (1) Create an array of control points from your cgx control points:
1336 for (i=0; i<nurbs[nr].u_npnt; i++)
1337 {
1338 for (j=0; j<nurbs[nr].v_npnt; j++)
1339 {
1340 controlPoints [ index++ ].components (
1341 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0],
1342 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1],
1343 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2],
1344 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3] );
1345 }
1346 }
1347
1348 // (2) Create an array of knots for each knot vector in your cgx surface:
1349 knot* knotsU = new knot [ nurbs[nr].u_nknt ];
1350 knot* knotsV = new knot [ nurbs[nr].v_nknt ];
1351
1352 for(i=0; i<nurbs[nr].u_nknt; i++) knotsU[i]=nurbs[nr].uknt[i];
1353 for(i=0; i<nurbs[nr].v_nknt; i++) knotsV[i]=nurbs[nr].vknt[i];
1354
1355 // (3) Create an snlSurface:
1356 // You need to supply degreeU, degreeV, sizeU and sizeV.
1357 snlSurface* surface = new snlSurface ( nurbs[nr].u_exp, nurbs[nr].v_exp, nurbs[nr].u_npnt, nurbs[nr].v_npnt,
1358 controlPoints, knotsU, knotsV );
1359
1360 sem_post(&sem_g);
1361
1362 // (4) Project some points to the surface:
1363 // Create an array of snlPoints to project onto the surface.
1364 snlPoint* toProject = new snlPoint [ sum_p ]; // You supply numPointsToProject.
1365
1366 // Fill the snlPoint array with data from cgx
1367 if(node==0)
1368 {
1369 sum_p=n=0;
1370 while((surf[snr].npgn-n))
1371 {
1372 /* create temporary nodes */
1373 n++;
1374 j=(int)surf[snr].pgn[n++];
1375 n+=3;
1376 for(k=0; k<j; k++)
1377 {
1378 toProject[sum_p].components (surf[snr].pgn[n], surf[snr].pgn[n+1], surf[snr].pgn[n+2]);
1379 sum_p++;
1380 n+=3;
1381 }
1382 }
1383 }
1384 else
1385 {
1386 if(sem_wait(&sem_n)) printf("Error in:sem_wait\n");
1387 for(i=0; i<sum_p; i++)
1388 {
1389 toProject[i].components ((*node)[surf[snr].nod[i]].nx, (*node)[surf[snr].nod[i]].ny, (*node)[surf[snr].nod[i]].nz);
1390 }
1391 if(sem_post(&sem_n)) printf("Error in:sem_post\n");
1392 }
1393
1394 // Project points to surface.
1395 returnedPntsPerPnt =MAX_RETURNED_PNTS_PER_PNT/5;
1396 isort = new int [ sum_p ];
1397
1398 finerProjection:;
1399 inverted =surface -> fastProject( toProject, sum_p, &sum_inverted, convergTol, normTol, maxPass, PROJ_SENSITIVITY, returnedPntsPerPnt );
1400
1401 //printf("sum_inverted:%d\n", sum_inverted);
1402 //for(i=0; i<sum_inverted; i++) printf("%d node:%d dist:%f\n", inverted[i].origPtIndex, surf[snr].nod[inverted[i].origPtIndex], inverted[i].dist);
1403
1404 /* store the index of the clostest entity */
1405 for(i=0; i<sum_p; i++) isort[i]=-1;
1406 for(i=0; i<sum_inverted; i++)
1407 {
1408 n=inverted[i].origPtIndex;
1409 if(isort[n]==-1) isort[n]=i;
1410 else if(inverted[i].dist<inverted[isort[n]].dist) isort[n]=i;
1411 }
1412 //printf("sum:%d\n", sum_p);
1413 //for(i=0; i<sum_p; i++) printf("%d node:%d dist:%f\n", isort[i], surf[snr].nod[inverted[isort[i]].origPtIndex], inverted[isort[i]].dist);
1414
1415 if(node==0)
1416 {
1417 i=n=0;
1418 while((surf[snr].npgn-n))
1419 {
1420 /* create temporary nodes */
1421 n++;
1422 j=(int)surf[snr].pgn[n++];
1423 n+=3;
1424 for(k=0; k<j; k++)
1425 {
1426 // check if the distance is "close"
1427 if((inverted[isort[i]].dist>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
1428 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; goto finerProjection; }
1429
1430 surf[snr].pgn[n] = inverted [ isort[i]].pt.x();
1431 surf[snr].pgn[n+1] = inverted [ isort[i]].pt.y();
1432 surf[snr].pgn[n+2] = inverted [ isort[i]].pt.z();
1433 i++;
1434 n+=3;
1435 }
1436 }
1437 }
1438 else
1439 {
1440 for(i=0; i<sum_p; i++)
1441 {
1442 // check if the distance is "close"
1443 if((inverted[isort[i]].dist>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
1444 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; goto finerProjection; }
1445 if(sem_wait(&sem_n)) printf("Error in:sem_wait\n");
1446 (*node)[surf[snr].nod[i]].nx = inverted[isort[i]].pt.x();
1447 (*node)[surf[snr].nod[i]].ny = inverted[isort[i]].pt.y();
1448 (*node)[surf[snr].nod[i]].nz = inverted[isort[i]].pt.z();
1449 if(sem_post(&sem_n)) printf("Error in:sem_post\n");
1450 }
1451 }
1452 delete[] isort;
1453 delete surface; // Release surface object.
1454 delete[] inverted; // Release snlVertex array returned from projection function.
1455 delete[] toProject; // Release points that were projected onto surface.
1456
1457 // The rest of the objects are deleted by libSNL.
1458 }
1459
1460
projSetToNurbs(int nr,Sets * set,int setNr,Points * pnt)1461 void projSetToNurbs( int nr, Sets *set, int setNr, Points *pnt)
1462 {
1463 int i,j,n,index=0;
1464
1465 double convergTol, normTol;
1466 int maxPass,sum_p, sum_inverted=0, returnedPntsPerPnt;
1467 snlSurfLocn* inverted;
1468
1469 sem_wait(&sem_g);
1470
1471 int numberOfControlPoints =nurbs[nr].u_npnt*nurbs[nr].v_npnt;
1472 snlCtrlPoint* controlPoints = new snlCtrlPoint [ numberOfControlPoints ];
1473
1474 normTol=COS_TOL;
1475 convergTol=ITER_TOL;
1476 maxPass=MAX_PASS;
1477
1478 // (1) Create an array of control points from your cgx control points:
1479 for (i=0; i<nurbs[nr].u_npnt; i++)
1480 {
1481 for (j=0; j<nurbs[nr].v_npnt; j++)
1482 {
1483 controlPoints [ index++ ].components (
1484 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0],
1485 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1],
1486 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2],
1487 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3] );
1488 }
1489 }
1490
1491 // (2) Create an array of knots for each knot vector in your cgx surface:
1492
1493 knot* knotsU = new knot [ nurbs[nr].u_nknt ];
1494 knot* knotsV = new knot [ nurbs[nr].v_nknt ];
1495
1496 for(i=0; i<nurbs[nr].u_nknt; i++) knotsU[i]=nurbs[nr].uknt[i];
1497 for(i=0; i<nurbs[nr].v_nknt; i++) knotsV[i]=nurbs[nr].vknt[i];
1498
1499 // (3) Create an snlSurface:
1500 // You need to supply degreeU, degreeV, sizeU and sizeV.
1501 snlSurface* surface = new snlSurface ( nurbs[nr].u_exp, nurbs[nr].v_exp, nurbs[nr].u_npnt, nurbs[nr].v_npnt,
1502 controlPoints, knotsU, knotsV );
1503 sem_post(&sem_g);
1504
1505 // (4) Project some points to the surface:
1506 // Create an array of snlPoints to project onto the surface.
1507 snlPoint* toProject = new snlPoint [ set[setNr].anz_p ]; // You supply numPointsToProject.
1508
1509 // Fill the snlPoint array with data from cgx
1510 for(i=0; i<set[setNr].anz_p; i++)
1511 {
1512 toProject[i].components (pnt[set[setNr].pnt[i]].px, pnt[set[setNr].pnt[i]].py, pnt[set[setNr].pnt[i]].pz);
1513 #if TEST
1514 printf("1xyz: %f %f %f\n",pnt[set[setNr].pnt[i]].px, pnt[set[setNr].pnt[i]].py, pnt[set[setNr].pnt[i]].pz);
1515 #endif
1516 }
1517
1518 // Project points to surface.
1519 sum_p=set[setNr].anz_p;
1520 returnedPntsPerPnt =MAX_RETURNED_PNTS_PER_PNT/5;
1521 int* isort = new int [ sum_p ];
1522
1523 finerProjection:;
1524 inverted =surface -> fastProject( toProject, sum_p, &sum_inverted, convergTol, normTol, maxPass, PROJ_SENSITIVITY, returnedPntsPerPnt );
1525
1526 //printf("sum_inverted:%d\n", sum_inverted);
1527 //for(i=0; i<sum_inverted; i++) printf("%d node:%d dist:%f\n", inverted[i].origPtIndex, surf[snr].nod[inverted[i].origPtIndex], inverted[i].dist);
1528
1529 /* store the index of the clostest entity */
1530 for(i=0; i<sum_p; i++) isort[i]=-1;
1531 for(i=0; i<sum_inverted; i++)
1532 {
1533 n=inverted[i].origPtIndex;
1534 if(isort[n]==-1) isort[n]=i;
1535 else if(inverted[i].dist<inverted[isort[n]].dist) isort[n]=i;
1536 }
1537 //printf("sum:%d\n", sum_p);
1538 //for(i=0; i<sum_p; i++) printf("%d node:%d dist:%f\n", isort[i], surf[snr].nod[inverted[isort[i]].origPtIndex], inverted[isort[i]].dist);
1539
1540 for(i=0; i<set[setNr].anz_p; i++)
1541 {
1542 // check if the distance is "close"
1543 if((inverted[isort[i]].dist>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
1544 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; goto finerProjection; }
1545
1546 pnt[set[setNr].pnt[i]].px = inverted [ isort[i]].pt.x();
1547 pnt[set[setNr].pnt[i]].py = inverted [ isort[i]].pt.y();
1548 pnt[set[setNr].pnt[i]].pz = inverted [ isort[i]].pt.z();
1549
1550 #if TEST
1551 printf("pnt P%d %f %f %f\n", i, pnt[set[setNr].pnt[i]].px, pnt[set[setNr].pnt[i]].py, pnt[set[setNr].pnt[i]].pz);
1552 #endif
1553 }
1554
1555 delete[] isort;
1556 delete surface; // Release surface object.
1557 delete[] inverted; // Release snlVertex array returned from projection function.
1558 delete[] toProject; // Release points that were projected onto surface.
1559
1560 // The rest of the objects are deleted by libSNL.
1561 }
1562
1563
1564
projPntsToNurbs(int nr,int anz_p,Points * pnt)1565 double *projPntsToNurbs( int nr, int anz_p, Points *pnt)
1566 {
1567 int i,j,n,index=0;
1568
1569 double convergTol, normTol;
1570 int maxPass,sum_p, sum_inverted=0, returnedPntsPerPnt;
1571 snlSurfLocn* inverted;
1572
1573 sem_wait(&sem_g);
1574
1575 int numberOfControlPoints =nurbs[nr].u_npnt*nurbs[nr].v_npnt;
1576 snlCtrlPoint* controlPoints = new snlCtrlPoint [ numberOfControlPoints ];
1577
1578 normTol=COS_TOL;
1579 convergTol=ITER_TOL;
1580 maxPass=MAX_PASS;
1581
1582
1583 // (1) Create an array of control points from your cgx control points:
1584 for (i=0; i<nurbs[nr].u_npnt; i++)
1585 {
1586 for (j=0; j<nurbs[nr].v_npnt; j++)
1587 {
1588 controlPoints [ index++ ].components (
1589 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0],
1590 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1],
1591 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2],
1592 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3] );
1593 }
1594 }
1595
1596 // (2) Create an array of knots for each knot vector in your cgx surface:
1597
1598 knot* knotsU = new knot [ nurbs[nr].u_nknt ];
1599 knot* knotsV = new knot [ nurbs[nr].v_nknt ];
1600
1601 for(i=0; i<nurbs[nr].u_nknt; i++) knotsU[i]=nurbs[nr].uknt[i];
1602 for(i=0; i<nurbs[nr].v_nknt; i++) knotsV[i]=nurbs[nr].vknt[i];
1603
1604 // (3) Create an snlSurface:
1605 // You need to supply degreeU, degreeV, sizeU and sizeV.
1606 snlSurface* surface = new snlSurface ( nurbs[nr].u_exp, nurbs[nr].v_exp, nurbs[nr].u_npnt, nurbs[nr].v_npnt,
1607 controlPoints, knotsU, knotsV );
1608
1609 sem_post(&sem_g);
1610
1611 // (4) Project some points to the surface:
1612 // Create an array of snlPoints to project onto the surface.
1613 snlPoint* toProject = new snlPoint [ anz_p ]; // You supply numPointsToProject.
1614
1615 // Fill the snlPoint array with data from cgx
1616 for(i=0; i<anz_p; i++)
1617 {
1618 toProject[i].components (pnt[i].px, pnt[i].py, pnt[i].pz);
1619 #if TEST
1620 printf("1xyz: %f %f %f\n",pnt[i].px, pnt[i].py, pnt[i].pz);
1621 #endif
1622 }
1623
1624 // Project points to surface.
1625 sum_p=anz_p;
1626 returnedPntsPerPnt =MAX_RETURNED_PNTS_PER_PNT/5;
1627 int* isort = new int [ sum_p ];
1628 double* dist = new double [ sum_p ];
1629
1630 finerProjection:;
1631 //printf("surface -> fastProject\n");
1632 inverted =surface -> fastProject( toProject, sum_p, &sum_inverted, convergTol, normTol, maxPass, PROJ_SENSITIVITY, returnedPntsPerPnt );
1633
1634 //printf("sum_inverted:%d\n", sum_inverted);
1635 //for(i=0; i<sum_inverted; i++) printf("%d node:%d dist:%f\n", inverted[i].origPtIndex, surf[snr].nod[inverted[i].origPtIndex], inverted[i].dist);
1636
1637 /* store the index of the clostest entity */
1638 for(i=0; i<sum_p; i++) isort[i]=-1;
1639 for(i=0; i<sum_inverted; i++)
1640 {
1641 n=inverted[i].origPtIndex;
1642 if(isort[n]==-1) isort[n]=i;
1643 else if(inverted[i].dist<inverted[isort[n]].dist) isort[n]=i;
1644 }
1645 //printf("sum:%d\n", sum_p);
1646 //for(i=0; i<sum_p; i++) printf("%d node:%d dist:%f\n", isort[i], surf[snr].nod[inverted[isort[i]].origPtIndex], inverted[isort[i]].dist);
1647
1648 for(i=0; i<anz_p; i++)
1649 {
1650 // check if the distance is "close"
1651 if((inverted[isort[i]].dist>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
1652 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; goto finerProjection; }
1653
1654 pnt[i].px = inverted [ isort[i]].pt.x();
1655 pnt[i].py = inverted [ isort[i]].pt.y();
1656 pnt[i].pz = inverted [ isort[i]].pt.z();
1657 dist[i] = inverted [ isort[i]].dist;
1658
1659 #if TEST
1660 printf("2xyz: %f %f %f\n",pnt[i].px, pnt[i].py, pnt[i].pz);
1661 #endif
1662 }
1663
1664 delete[] isort;
1665 delete surface; // Release surface object.
1666 delete[] inverted; // Release snlVertex array returned from projection function.
1667 delete[] toProject; // Release points that were projected onto surface.
1668
1669 // The rest of the objects are deleted by libSNL.
1670 return(dist);
1671 }
1672
1673
1674 // Attention: return pnt[5] (xyz,uv)
proj1PntToNurbs(int nr,double * pnt)1675 double proj1PntToNurbs( int nr, double *pnt)
1676 {
1677 int i,j,n,index=0;
1678
1679 double convergTol, normTol, result=0.;
1680 int maxPass,sum_p, sum_inverted=0, returnedPntsPerPnt;
1681 snlSurfLocn* inverted;
1682
1683 sem_wait(&sem_g);
1684
1685 int numberOfControlPoints =nurbs[nr].u_npnt*nurbs[nr].v_npnt;
1686 snlCtrlPoint* controlPoints = new snlCtrlPoint [ numberOfControlPoints ];
1687
1688 normTol=COS_TOL;
1689 convergTol=ITER_TOL;
1690 maxPass=MAX_PASS;
1691
1692 // (1) Create an array of control points from your cgx control points:
1693 for (i=0; i<nurbs[nr].u_npnt; i++)
1694 {
1695 for (j=0; j<nurbs[nr].v_npnt; j++)
1696 {
1697 controlPoints [ index++ ].components (
1698 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0],
1699 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1],
1700 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2],
1701 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3] );
1702 }
1703 }
1704
1705 // (2) Create an array of knots for each knot vector in your cgx surface:
1706
1707 knot* knotsU = new knot [ nurbs[nr].u_nknt ];
1708 knot* knotsV = new knot [ nurbs[nr].v_nknt ];
1709
1710 for(i=0; i<nurbs[nr].u_nknt; i++) knotsU[i]=nurbs[nr].uknt[i];
1711 for(i=0; i<nurbs[nr].v_nknt; i++) knotsV[i]=nurbs[nr].vknt[i];
1712
1713 // (3) Create an snlSurface:
1714 // You need to supply degreeU, degreeV, sizeU and sizeV.
1715 snlSurface* surface = new snlSurface ( nurbs[nr].u_exp, nurbs[nr].v_exp, nurbs[nr].u_npnt, nurbs[nr].v_npnt,
1716 controlPoints, knotsU, knotsV );
1717
1718 sem_post(&sem_g);
1719
1720 // (4) Project some points to the surface:
1721 // Create an array of snlPoints to project onto the surface.
1722 snlPoint* toProject = new snlPoint [ 1 ]; // You supply numPointsToProject.
1723
1724 // Fill the snlPoint array with data from cgx
1725 toProject[0].components (pnt[0], pnt[1], pnt[2]);
1726
1727 // Project points to surface.
1728 sum_p=1;
1729 returnedPntsPerPnt =MAX_RETURNED_PNTS_PER_PNT/5;
1730 int* isort = new int [ sum_p ];
1731
1732 finerProjection:;
1733 inverted =surface -> fastProject( toProject, sum_p, &sum_inverted, convergTol, normTol, maxPass, PROJ_SENSITIVITY, returnedPntsPerPnt );
1734 //
1735 //printf("sum_inverted:%d\n", sum_inverted);
1736 //for(i=0; i<sum_inverted; i++) printf("%d dist:%f\n", inverted[i].origPtIndex, inverted[i].dist);
1737
1738 /* store the index of the clostest entity */
1739 isort[0]=-1;
1740 for(i=0; i<sum_inverted; i++)
1741 {
1742 n=inverted[i].origPtIndex;
1743 if(isort[n]==-1) isort[n]=i;
1744 else if(inverted[i].dist<inverted[isort[n]].dist) isort[n]=i;
1745 }
1746 //
1747 //printf("sum:%d\n", sum_p);
1748 //for(i=0; i<sum_p; i++) printf("%d node:%d dist:%f\n", isort[i],inverted[isort[i]].origPtIndex, inverted[isort[i]].dist);
1749
1750 i=0;
1751 {
1752 // check if the distance is "close"
1753 if((inverted[isort[i]].dist>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
1754 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; goto finerProjection; }
1755
1756 pnt[0] = inverted [ isort[i]].pt.x();
1757 pnt[1] = inverted [ isort[i]].pt.y();
1758 pnt[2] = inverted [ isort[i]].pt.z();
1759 pnt[3] = inverted [ isort[i]].paramU;
1760 pnt[4] = inverted [ isort[i]].paramV;
1761 }
1762
1763 result=inverted[isort[i]].dist;
1764 delete[] isort;
1765 delete surface; // Release surface object.
1766 delete[] inverted; // Release snlVertex array returned from projection function.
1767 delete[] toProject; // Release points that were projected onto surface.
1768
1769 return(result);
1770 // The rest of the objects are deleted by libSNL.
1771 }
1772
1773
1774
1775 // The "ball" type is difficult if the trimming loop extend in the ambiguous zone. To avoid this the nurbs is rotated in space.
1776 // return 1 if the nurbs was rotated in space
rotateBall(int nr,int axis,int patch,double utol_ambig,double vtol_ambig)1777 int rotateBall( int nr, int axis, int patch, double utol_ambig, double vtol_ambig)
1778 {
1779 int i=0,j=0,k=0,n=0, p;
1780 int ipax1=0, ipax2=0, pntset=0, curve=0, atPole=0;
1781 double u=0.5, angle[2]={0.,0.};
1782 double pax1[3]={0.,0.,0.}, pax2[3]={0.,0.,0.}, p1p2[3]={0.,0.,0.}, cg_sphere[3]={0.,0.,0.}, cg_patch[3]={0.,0.,0.}, vax1[3]={0.,0.,0.}, vax2[3]={0.,0.,0.}, vcgp[3]={0.,0.,0.};
1783 char buffer[MAX_LINE_LENGTH];
1784 int foundPntAtPole[2]={0,0}, foundPntAtMeridian[2]={0,0};
1785 double umin, umax,vmin,vmax, du,dv;
1786 double pumin, pumax,pvmin,pvmax;
1787 int *ptmp=NULL;
1788 extern Points *point;
1789
1790
1791 // determine the axis cross
1792
1793 // 1.axis defined by the points at the poles
1794 if(axis==1) //u is axis, pax1 at top
1795 {
1796 sem_wait(&sem_g);
1797 i=nurbs[nr].u_npnt-1;j=0;
1798 pax1[0]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0];
1799 pax1[1]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1];
1800 pax1[2]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2];
1801 i=0;
1802 pax2[0]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0];
1803 pax2[1]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1];
1804 pax2[2]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2];
1805 sem_post(&sem_g);
1806 }
1807 else if(axis==2) //v is axis, pax1 at top
1808 {
1809 sem_wait(&sem_g);
1810 j=nurbs[nr].v_npnt-1;i=0;
1811 pax1[0]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0];
1812 pax1[1]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1];
1813 pax1[2]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2];
1814 j=0;
1815 pax2[0]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0];
1816 pax2[1]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1];
1817 pax2[2]= nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2];
1818 sem_post(&sem_g);
1819 }
1820 else return(0);
1821
1822 v_result( pax1, pax2, p1p2 );
1823 v_scal( &u, p1p2, vax1);
1824 v_add( pax1, vax1, cg_sphere);
1825
1826 // 2.axis defined by the cross-product of v_1.axis and v_CGsurf
1827
1828 // CG-surf (trimming loop)
1829 // over all points of the outer curve
1830 // scip the last point
1831 j=0;
1832 sem_wait(&sem_g);
1833 for(k=0; k<nurbs[nr].np[patch][0]-1; k++)
1834 {
1835 for(n=0;n<3; n++)
1836 {
1837 cg_patch[n]+=nurbs[nr].xyz[patch][0][j++];
1838 }
1839 }
1840 sem_post(&sem_g);
1841 for(n=0;n<3; n++) cg_patch[n]/=j/3;
1842
1843 v_result( cg_sphere, cg_patch, vcgp );
1844 v_prod( vcgp, vax1, vax2);
1845 v_add(cg_sphere, vax2, pax2);
1846
1847 // determine the nurbs-u,v range by looking into the knots (to detect ambiguous points)
1848 sem_wait(&sem_g);
1849 umin=nurbs[nr].uknt[0];
1850 umax=nurbs[nr].uknt[nurbs[nr].u_nknt-1];
1851 vmin=nurbs[nr].vknt[0];
1852 vmax=nurbs[nr].vknt[nurbs[nr].v_nknt-1];
1853 sem_post(&sem_g);
1854
1855 // determine the rotations. "0" around axis through the poles (1.axis)
1856 i=0;
1857 pumax=pvmax=-MAXVALUE;
1858 pumin=pvmin=MAXVALUE;
1859 sem_wait(&sem_g);
1860 if(axis==1)
1861 {
1862 for(curve=0; curve<nurbs[nr].nc[patch]; curve++)
1863 for(j=0; j<nurbs[nr].np[patch][curve]; j++)
1864 {
1865 if( nurbs[nr].uv[patch][curve][i] < pumin ) pumin=nurbs[nr].uv[patch][curve][i];
1866 if( nurbs[nr].uv[patch][curve][i] > pumax ) pumax=nurbs[nr].uv[patch][curve][i];
1867 atPole=0;
1868 if( (nurbs[nr].uv[patch][curve][i]- utol_ambig ) < umin ) { foundPntAtPole[0]=1; atPole=1; }
1869 if( (nurbs[nr].uv[patch][curve][i]+ utol_ambig ) > umax ) { foundPntAtPole[1]=1; atPole=1; }
1870 i++;
1871
1872 if( nurbs[nr].uv[patch][curve][i] < pvmin ) pvmin=nurbs[nr].uv[patch][curve][i];
1873 if( nurbs[nr].uv[patch][curve][i] > pvmax ) pvmax=nurbs[nr].uv[patch][curve][i];
1874 if(!atPole)
1875 {
1876 if( (nurbs[nr].uv[patch][curve][i]- vtol_ambig ) < vmin ) foundPntAtMeridian[0]=1;
1877 if( (nurbs[nr].uv[patch][curve][i]+ vtol_ambig ) > vmax ) foundPntAtMeridian[1]=1;
1878 }
1879 i++;
1880 }
1881 }
1882 if(axis==2)
1883 {
1884 for(curve=0; curve<nurbs[nr].nc[patch]; curve++)
1885 for(j=0; j<nurbs[nr].np[patch][curve]; j++)
1886 {
1887 if( nurbs[nr].uv[patch][curve][i] < pumin ) pumin=nurbs[nr].uv[patch][curve][i];
1888 if( nurbs[nr].uv[patch][curve][i] > pumax ) pumax=nurbs[nr].uv[patch][curve][i];
1889 if( nurbs[nr].uv[patch][curve][i+1] < pvmin ) pvmin=nurbs[nr].uv[patch][curve][i+1];
1890 if( nurbs[nr].uv[patch][curve][i+1] > pvmax ) pvmax=nurbs[nr].uv[patch][curve][i+1];
1891 atPole=0;
1892 if( (nurbs[nr].uv[patch][curve][i+1]- vtol_ambig ) < vmin ) { foundPntAtPole[0]=1; atPole=1; }
1893 if( (nurbs[nr].uv[patch][curve][i+1]+ vtol_ambig ) > vmax ) { foundPntAtPole[1]=1; atPole=1; }
1894 if(!atPole)
1895 {
1896 if( (nurbs[nr].uv[patch][curve][i]- utol_ambig ) < umin ) foundPntAtMeridian[0]=1;
1897 if( (nurbs[nr].uv[patch][curve][i]+ utol_ambig ) > umax ) foundPntAtMeridian[1]=1;
1898 }
1899 i+=2;
1900 }
1901 }
1902 sem_post(&sem_g);
1903
1904 // determine the uv displacement necessary to move the patch in the middle of the nurbs
1905 du=0.5*((umax-umin) - (pumax-pumin)) + (umin-pumin);
1906 dv=0.5*((vmax-vmin) - (pvmax-pvmin)) + (vmin-pvmin);
1907 //du=(0.5*(pumin-umin) + (umax-pumax)) - (pumin-umin);
1908 //dv=(0.5*(pvmin-vmin) + (vmax-pvmax)) - (pvmin-vmin);
1909
1910
1911 // convert in rotation (deg):
1912 if(axis==1) //axis is u
1913 {
1914 angle[0]=dv*360./(vmax-vmin);
1915 // check if it is likely that v=0 is in the surf (the surf extends nearly around the aequator)
1916 if( ((pvmax-pvmin)/(vmax-vmin)) > 0.9) angle[0]+=180.;
1917 angle[1]=du*180./(umax-umin);
1918 }
1919 else //axis is v
1920 {
1921 angle[1]=dv*180./(vmax-vmin);
1922 angle[0]=du*360./(umax-umin);
1923 // check if it is likely that u=0 is in the surf (the surf extends nearly around the aequator)
1924 if( ((pumax-pumin)/(umax-umin)) > 0.9) angle[0]+=180.;
1925 }
1926
1927 #if TEST
1928 printf(" axis:%d uvtol_ambig:%f %f pu: %f %f pv: %f %f ball rotated by angle: %f %f duv:%f %f\n", axis, utol_ambig, vtol_ambig, pumin,pumax, pvmin,pvmax, angle[0],angle[1], du,dv);
1929 #endif
1930
1931 if((abs(angle[0])<1.)&&(abs(angle[1])<1.)) return(0);
1932
1933 // rotate the control-points. Do not rotate the original points because a substitute nurbs could be referenced and they uses modified coordinates
1934 if ( (ptmp = (int *)realloc((int *)ptmp, nurbs[nr].u_npnt*nurbs[nr].v_npnt* sizeof(int))) == NULL )
1935 printf("\n\n ERROR: realloc failed\n\n") ;
1936
1937 sem_wait(&sem_g);
1938 delSet("-tmp");
1939 pntset=pre_seta("-tmp1","i",0);
1940 k=0;
1941 for (i=0; i<nurbs[nr].u_npnt; i++)
1942 {
1943 for (j=0; j<nurbs[nr].v_npnt; j++)
1944 {
1945 getNewName( buffer, "p" );
1946 p=pnt( buffer, nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0],
1947 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1],
1948 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2], 0 );
1949 seta(pntset,"p",p);
1950 ptmp[k++]=p;
1951 }
1952 }
1953 sem_post(&sem_g);
1954
1955 if(angle[0]!=0.)
1956 {
1957 /* generate points on rotation-axis */
1958 sem_wait(&sem_g);
1959 ipax1=pnt( "+pax11",cg_sphere[0],cg_sphere[1],cg_sphere[2], 0 );
1960 ipax2=pnt( "+pax12",pax1[0],pax1[1],pax1[2], 0 );
1961 sem_post(&sem_g);
1962 sprintf(buffer,"-tmp1 rot +pax11 +pax12 %f", angle[0] );
1963 pre_move_pfast(buffer);
1964 sem_wait(&sem_g);
1965 delPnt( 1, &ipax1 );
1966 delPnt( 1, &ipax2 );
1967 sem_post(&sem_g);
1968 }
1969 if(angle[1]!=0.)
1970 {
1971 /* generate points on rotation-axis */
1972 sem_wait(&sem_g);
1973 ipax1=pnt( "+pax21",cg_sphere[0],cg_sphere[1],cg_sphere[2], 0 );
1974 ipax2=pnt( "+pax22",pax2[0],pax2[1],pax2[2], 0 );
1975 sem_post(&sem_g);
1976 sprintf(buffer,"-tmp1 rot +pax21 +pax22 %f", angle[1] );
1977 pre_move_pfast(buffer);
1978 sem_wait(&sem_g);
1979 delPnt( 1, &ipax1 );
1980 delPnt( 1, &ipax2 );
1981 sem_post(&sem_g);
1982 }
1983
1984 // adapt the nurbs definition
1985 k=0;
1986 sem_wait(&sem_g);
1987 for (i=0; i<nurbs[nr].u_npnt; i++)
1988 {
1989 for (j=0; j<nurbs[nr].v_npnt; j++)
1990 {
1991 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]=point[ptmp[k]].px;
1992 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]=point[ptmp[k]].py;
1993 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]=point[ptmp[k++]].pz;
1994 }
1995 }
1996
1997 delPnt( k, ptmp );
1998 delSet("-tmp1");
1999 sem_post(&sem_g);
2000 return(1);
2001 }
2002
2003
2004
2005 /* returns the type ( -1error, 0plate, 1cyl, 2torus, 3ball, 4half-ball-bottom, 5half-ball-top */
2006 /* ( half-ball-top: closed at top) */
2007 /* axis=1 ->u (ambig point has either vmax or vmin) */
2008 /* axis=2 ->v (ambig point has either umax or umin) */
2009 /* axis=3 ->uv (ambig point has either umax or umin and either vmax or vmin) */
2010
2011 /* ambig=1 if the end-points of one axis are closer than ITER_TOL in the same xyz location (ie a cylinder) */
2012 /* ambig[0]==1 u-ambig at vmin axis=+2 */
2013 /* ambig[1]==1 u-ambig at vmax axis=+2 */
2014 /* ambig[2]==1 v-ambig at umin axis=1 */
2015 /* ambig[3]==1 v-ambig at umax axis=1 */
2016
2017 /* collapsed=1 if one ambiguous edge is collapsed (a point, top of a ball or cone) */
2018 /* this is the case if all control points on that edge are closer than ITER_TOL in the same xyz location */
2019 /* collapsed[0]==1 u-collapsed at vmin axis=2 */
2020 /* collapsed[1]==1 u-collapsed at vmax axis=2 */
2021 /* collapsed[2]==1 v-collapsed at umin axis=1 */
2022 /* collapsed[3]==1 v-collapsed at umax axis=1 */
getNurbsType(int nr,int * axis,int * sector)2023 int getNurbsType( int nr, int *axis, int *sector)
2024 {
2025 int i,j,k,n;
2026 int collapsed[6], ambig[6];
2027 #if TEST2
2028 extern Points *point;
2029 #endif
2030
2031 double diff1;
2032
2033 *axis=0;
2034 *sector=0;
2035 for(k=0; k<6; k++) collapsed[k]=ambig[k]=0;
2036
2037
2038 /* check if one edge is ambiguous */
2039 /* this is the case if control points meet each other */
2040
2041 sem_wait(&sem_g);
2042
2043 /* check v at umin and umax */
2044 for(k=2; k<4; k++)
2045 {
2046 if(k==2) i=0; else i=nurbs[nr].u_npnt-1;
2047 #if TEST2
2048 printf("check at umin and umax\n");
2049 j=0;
2050 printf("x:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]);
2051 printf("y:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]);
2052 printf("z:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]);
2053 printf("w:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2054 printf(" i(u):%d j(v):%d %s\n", i, j, point[nurbs[nr].ctlpnt[i][j]].name);
2055 #endif
2056
2057 /* check the first and last only */
2058 /* if they are equal its an ambiguous edge */
2059 /* difference to the first point */
2060 ambig[k]=1;
2061 j=nurbs[nr].v_npnt-1;
2062 #if TEST2
2063 printf("x:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]);
2064 printf("y:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]);
2065 printf("z:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]);
2066 printf("w:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2067 printf(" i:%d j:%d %s\n", i, j, point[nurbs[nr].ctlpnt[i][j]].name);
2068 #endif
2069 for(n=0; n<3; n++)
2070 {
2071 diff1=( nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+n] / nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+3] -
2072 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n] / nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2073 if( diff1*diff1 >ITER_TOL*ITER_TOL ) ambig[k]=0;
2074 #if TEST2
2075 printf("diff1:%f ambig:%d\n",diff1,ambig[k]);
2076 #endif
2077 }
2078
2079 /* check if one ambiguous edge is collapsed (a point!) */
2080 /* this is the case if all control points on that edge are on the same spot */
2081 if(ambig[k])
2082 {
2083 collapsed[k]=1;
2084 for(j=1; j<nurbs[nr].v_npnt; j++)
2085 {
2086 #if TEST2
2087 printf("x:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]);
2088 printf("y:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]);
2089 printf("z:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]);
2090 printf("w:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2091 printf(" %s\n", point[nurbs[nr].ctlpnt[i][j]].name);
2092 #endif
2093
2094 /* difference to the first point */
2095 for(n=0; n<3; n++)
2096 {
2097 //diff1=( nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+n] -
2098 // nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n]);
2099 diff1=( nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+n] / nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+3] -
2100 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n] / nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2101 if( diff1*diff1 >ITER_TOL*ITER_TOL ) collapsed[k]=0;
2102 }
2103 }
2104 }
2105 }
2106
2107 /* check v at u in the middle */
2108 if(nurbs[nr].u_npnt>2)
2109 {
2110 k=4;
2111 i=nurbs[nr].u_npnt/2;
2112 #if TEST2
2113 printf("check at umean\n");
2114 j=0;
2115 printf("x:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]);
2116 printf("y:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]);
2117 printf("z:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]);
2118 printf("w:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2119 printf(" i:%d j:%d %s\n", i, j, point[nurbs[nr].ctlpnt[i][j]].name);
2120 #endif
2121
2122 /* check the first and last only */
2123 /* if they are equal its an ambiguous edge */
2124 /* difference to the first point */
2125 ambig[k]=1;
2126 j=nurbs[nr].v_npnt-1;
2127 #if TEST2
2128 printf("x:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]);
2129 printf("y:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]);
2130 printf("z:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]);
2131 printf("w:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2132 printf(" i:%d j:%d %s\n", i, j, point[nurbs[nr].ctlpnt[i][j]].name);
2133 #endif
2134 for(n=0; n<3; n++)
2135 {
2136 //diff1=( nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+n] -
2137 // nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n]);
2138 diff1=( nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+n] / nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+3] -
2139 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n] / nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2140 if( diff1*diff1 >ITER_TOL*ITER_TOL ) ambig[k]=0;
2141 #if TEST2
2142 printf("diff1:%f ambig:%d\n",diff1,ambig[k]);
2143 #endif
2144 }
2145 }
2146
2147 /* check u at vmin and vmax */
2148 for(k=0; k<2; k++)
2149 {
2150 if(k==0) j=0; else j=nurbs[nr].v_npnt-1;
2151
2152 #if TEST2
2153 printf("check at vmin and vmax\n");
2154 i=0;
2155 printf("x:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]);
2156 printf("y:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]);
2157 printf("z:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]);
2158 printf("w:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2159 printf(" i:%d j:%d %s\n", i, j, point[nurbs[nr].ctlpnt[i][j]].name);
2160 #endif
2161 /* check the first and last only */
2162 /* if they are equal its an ambiguous edge */
2163 /* difference to the first point */
2164 ambig[k]=1;
2165 i=nurbs[nr].u_npnt-1;
2166 #if TEST2
2167 printf("x:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]);
2168 printf("y:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]);
2169 printf("z:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]);
2170 printf("w:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2171 printf(" i:%d j:%d %s\n", i, j, point[nurbs[nr].ctlpnt[i][j]].name);
2172 #endif
2173 for(n=0; n<3; n++)
2174 {
2175 //diff1=( nurbs[nr].ctlarray[j*(nurbs[nr].v_stride)+n] -
2176 // nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n]);
2177 diff1=( nurbs[nr].ctlarray[j*(nurbs[nr].v_stride)+n] / nurbs[nr].ctlarray[j*(nurbs[nr].v_stride)+3] -
2178 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n] / nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2179 if( diff1*diff1 >ITER_TOL*ITER_TOL ) ambig[k]=0;
2180 #if TEST2
2181 printf("diff1:%f ambig:%d\n",diff1,ambig[k]);
2182 #endif
2183 }
2184
2185 /* check if one ambiguous edge is collapsed (a point!) */
2186 /* this is the case if all control points on that edge are on the same spot */
2187 if(ambig[k])
2188 {
2189 collapsed[k]=1;
2190 for(i=1; i<nurbs[nr].u_npnt; i++)
2191 {
2192 #if TEST2
2193 printf("x:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]);
2194 printf("y:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]);
2195 printf("z:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]);
2196 printf("w:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2197 printf(" i:%d j:%d %s\n", i, j, point[nurbs[nr].ctlpnt[i][j]].name);
2198 #endif
2199
2200 /* difference to the first point */
2201 for(n=0; n<3; n++)
2202 {
2203 diff1=(nurbs[nr].ctlarray[j*(nurbs[nr].v_stride)+n] / nurbs[nr].ctlarray[j*(nurbs[nr].v_stride)+3] -
2204 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n] / nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2205 if( diff1*diff1 >ITER_TOL*ITER_TOL ) collapsed[k]=0;
2206 }
2207 }
2208 }
2209 }
2210
2211 /* check u at v in the middle */
2212 if(nurbs[nr].v_npnt>2)
2213 {
2214 k=5;
2215 i=nurbs[nr].v_npnt/2;
2216 #if TEST2
2217 printf("check at vmean\n");
2218 i=0;
2219 printf("x:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]);
2220 printf("y:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]);
2221 printf("z:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]);
2222 printf("w:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2223 printf(" i:%d j:%d %s\n", i, j, point[nurbs[nr].ctlpnt[i][j]].name);
2224 #endif
2225 /* check the first and last only */
2226 /* if they are equal its an ambiguous edge */
2227 /* difference to the first point */
2228 ambig[k]=1;
2229 i=nurbs[nr].u_npnt-1;
2230 #if TEST2
2231 printf("x:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0]);
2232 printf("y:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1]);
2233 printf("z:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2]);
2234 printf("w:%f ", nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2235 printf(" i:%d j:%d %s\n", i, j, point[nurbs[nr].ctlpnt[i][j]].name);
2236 #endif
2237 for(n=0; n<4; n++)
2238 {
2239 //diff1=(nurbs[nr].ctlarray[j*(nurbs[nr].v_stride)+n] -
2240 // nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n]);
2241 diff1=(nurbs[nr].ctlarray[j*(nurbs[nr].v_stride)+n] / nurbs[nr].ctlarray[j*(nurbs[nr].v_stride)+3] -
2242 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n] / nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3]);
2243 if( diff1*diff1 >ITER_TOL*ITER_TOL ) ambig[k]=0;
2244 #if TEST2
2245 printf("diff1:%f ambig:%d\n",diff1,ambig[k]);
2246 #endif
2247 }
2248 }
2249
2250 sem_post(&sem_g);
2251
2252 if((ambig[2]==1)||(ambig[3]==1)) *axis=1;
2253 if((ambig[0]==1)||(ambig[1]==1)) *axis+=2;
2254 #if TEST2
2255 i=0; for(k=0; k<4; k++)
2256 {
2257 printf("%d collapsed=%d (1:yes)\n",k, collapsed[k]);
2258 printf(" ambig=%d (1:yes)\n", ambig[k]);
2259 printf(" axis=%d (1:2)\n", *axis);
2260 }
2261 #endif
2262 /* determine the case */
2263 /* -1error, 0plate, 1cyl, 2torus, 3ball, 4half-ball-bot, 5half-ball-top */
2264 /* ( half-ball-top: closed at top) */
2265 i=0; for(k=0; k<4; k++)
2266 {
2267 i+=collapsed[k];
2268 }
2269 if(i==0)
2270 {
2271 if((*axis==1)||(*axis==2)) return(1);
2272 else if(*axis==3) return(2);
2273 else return(0);
2274 }
2275 else
2276 {
2277 if((collapsed[0]==1)&&(collapsed[1]==1)) { *axis=2; if(ambig[5]==0) *sector=1; return(3); }
2278 else if((collapsed[2]==1)&&(collapsed[3]==1)) { *axis=1; if(ambig[4]==0) *sector=1; return(3); }
2279 else if((collapsed[0]==1)&&(collapsed[1]==0)) { *axis=2; if(ambig[1]==0) *sector=1; return(4); }
2280 else if((collapsed[2]==1)&&(collapsed[3]==0)) { *axis=1; if(ambig[3]==0) *sector=1; return(4); }
2281 else if((collapsed[0]==0)&&(collapsed[1]==1)) { *axis=2; if(ambig[0]==0) *sector=1; return(5); }
2282 else if((collapsed[2]==0)&&(collapsed[3]==1)) { *axis=1; if(ambig[2]==0) *sector=1; return(5); }
2283 else return(-1);
2284 }
2285 }
2286
2287
2288
trimNurbs(int nr,int patch,double ini_tol_ambig)2289 int trimNurbs( int nr, int patch, double ini_tol_ambig)
2290 {
2291 int i,j,k,n,nn,p;
2292 int nbuf[4];
2293 int curve, sum_p=0, sum_inverted=0, nurbsType=0,index=0, axis=0, sectorFlag=0, buf=0, returnedPntsPerPnt ;
2294 int orig_np=0, pcollapsed, collapsedFlag=0, moveFlag=1 ;
2295
2296 double convergTol, normTol;
2297 int maxPass;
2298
2299 sem_wait(&sem_g);
2300 int numberOfControlPoints =nurbs[nr].u_npnt*nurbs[nr].v_npnt;
2301 sem_post(&sem_g);
2302
2303 snlCtrlPoint* controlPoints=NULL;
2304 knot* knotsU=NULL;
2305 knot* knotsV=NULL;
2306 snlSurface* surface=NULL;
2307
2308 double v[3],maxv[3], minv[3];
2309
2310 int *nppc;
2311 int disti=0, startFlag, pskip, nextUflag=0, nextVflag=0, lastUflag=-1, lastVflag=-1;
2312 double umin, umax,vmin,vmax, du, dv, distp;
2313 double dist1=0., dist2=0., tol_ambig=0., utol_ambig, vtol_ambig;
2314 double av_local_elength,p0[3],p1[3],p0p1[3];
2315
2316 PntProj *pntproj=NULL;
2317 snlSurfLocn* inverted=NULL;
2318
2319 normTol=COS_TOL;
2320 convergTol=ITER_TOL;
2321 maxPass=MAX_PASS;
2322
2323 sem_wait(&sem_g);
2324 if(nurbs[nr].endFlag==0)
2325 {
2326 printf(" ERROR: nurbs not valid\n",nurbs[nr].name);
2327 sem_post(&sem_g);
2328 return(-1);
2329 }
2330 if(nurbs[nr].ctlarray==NULL)
2331 {
2332 printf(" ERROR: programmers error: nurbs ctlarray was not initialized\n");
2333 sem_post(&sem_g);
2334 exit(0);
2335 }
2336 for(i=0; i<nurbs[nr].nc[patch]; i++) sum_p+=nurbs[nr].np[patch][i];
2337 sem_post(&sem_g);
2338 if(sum_p==0) return(-1);
2339
2340 // (4) Create an array of snlPoints to project onto the surface.
2341 snlPoint* toProject = new snlPoint [ sum_p ]; // You supply numPointsToProject.
2342
2343 // determine the nurbs-u,v range by looking into the knots (to detect ambiguous points)
2344 sem_wait(&sem_g);
2345 umin=nurbs[nr].uknt[0];
2346 umax=nurbs[nr].uknt[nurbs[nr].u_nknt-1];
2347 vmin=nurbs[nr].vknt[0];
2348 vmax=nurbs[nr].vknt[nurbs[nr].v_nknt-1];
2349 sem_post(&sem_g);
2350 du=umax-umin;
2351 dv=vmax-vmin;
2352 #if TEST
2353 printf(" u_nknt:%d v_nknt:%d\n",nurbs[nr].u_nknt,nurbs[nr].v_nknt );
2354 printf(" umax:%f umin:%f vmax:%f vmin:%f\n", umax,umin,vmax,vmin);
2355 #endif
2356
2357 // Fill the snlPoint array with data from cgx
2358 // and determine the patch-x,y,z range by looking into the points (tolerance of ambig)
2359 for(n=0; n<3; n++) { maxv[n]=-MAXVALUE; minv[n]=MAXVALUE; }
2360 p=0;
2361
2362 sem_wait(&sem_g);
2363 nbuf[0]=nurbs[nr].nc[patch];
2364 sem_post(&sem_g);
2365 for(curve=0; curve<nbuf[0]; curve++)
2366 {
2367 j=0;
2368 sem_wait(&sem_g);
2369 nbuf[1]=nurbs[nr].np[patch][curve];
2370 sem_post(&sem_g);
2371 for(i=0; i<nbuf[1]; i++)
2372 {
2373 for(n=0;n<3; n++)
2374 {
2375 sem_wait(&sem_g);
2376 v[n]=nurbs[nr].xyz[patch][curve][j++];
2377 sem_post(&sem_g);
2378 if(maxv[n]<v[n]) maxv[n]=v[n];
2379 if(minv[n]>v[n]) minv[n]=v[n];
2380 }
2381 toProject[p++].components (v[0],v[1],v[2]);
2382 #if TEST
2383 printf("PNT P%d %f %f %f\n",p,v[0],v[1],v[2]);
2384 #endif
2385 }
2386 }
2387 // else
2388 //for(n=0; n<3; n++) dist1+=maxv[n]-minv[n];
2389 // or
2390 dist1=-MAXVALUE;
2391 for(n=0; n<3; n++)
2392 {
2393 maxv[n]-=minv[n];
2394 if(dist1<maxv[n]) dist1=maxv[n];
2395 }
2396
2397 #if PRINT_SNL
2398 surface ->print_cpp();
2399 exit(-1);
2400 #endif
2401
2402 // determine the nurbs-x,y,z range by looking into the points (tolerance of ambig)
2403 for(n=0; n<3; n++) { maxv[n]=-MAXVALUE; minv[n]=MAXVALUE; }
2404 sem_wait(&sem_g);
2405 for(i=0; i<nurbs[nr].u_npnt; i++)
2406 {
2407 for(j=1; j<nurbs[nr].v_npnt; j++)
2408 {
2409 for(n=0; n<3; n++)
2410 {
2411 if(maxv[n]<nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n]) maxv[n]=nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n];
2412 if(minv[n]>nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n]) minv[n]=nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+n];
2413 }
2414 }
2415 }
2416 sem_post(&sem_g);
2417 // else
2418 //for(n=0; n<3; n++) dist2+=maxv[n]-minv[n];
2419 // or
2420 dist2=-MAXVALUE;
2421 for(n=0; n<3; n++)
2422 {
2423 maxv[n]-=minv[n];
2424 if(dist2<maxv[n]) dist2=maxv[n];
2425 }
2426 tol_ambig=dist1/dist2*ini_tol_ambig;
2427 utol_ambig=tol_ambig*du;
2428 vtol_ambig=tol_ambig*dv;
2429
2430 #if TEST
2431 printf("extension of region (xyz) dist1:%f\n", dist1);
2432 printf("extension of nurbs (xyz) dist2:%f\n", dist2);
2433 printf("dist1/dist2*ini_tol_ambig = tol_ambig:%f\n", tol_ambig);
2434 #endif
2435
2436 // temporary buffer for the projected points
2437 pntproj= new PntProj [sum_p];
2438
2439 // additional number of points per curve object
2440 sem_wait(&sem_g);
2441 nppc= new int [nurbs[nr].nc[patch]+1];
2442 sem_post(&sem_g);
2443 for(n=0; n<=nurbs[nr].nc[patch]; n++) nppc[n]=0;
2444
2445
2446 // (5) Project points to surface.
2447 if(printFlag) printf("inversion of %d points to nurbs:%s\n", sum_p, nurbs[nr].name);
2448 returnedPntsPerPnt =MAX_RETURNED_PNTS_PER_PNT/5;
2449
2450
2451 // determine the type of the nurbs
2452 // (0plate, 1cyl, 2torus, 3ball, 4half-ball-bot, 5half-ball-top)
2453 // ( half-ball-top: closed at top)
2454 // axis=1 v is ambiguous (point has either vmax or vmin)
2455 // axis=2 u is ambiguous (point has either umax or umin)
2456 // axis=3 u and v are ambiguous (torus)
2457 if((nurbsType=getNurbsType( nr, &axis, §orFlag ))<0)
2458 { printf("ERROR: type of nurbs:%d not known\n",nurbsType); exit(-1); }
2459 #if TEST
2460 printf("nurbsType:%d (0plate, 1cyl, 2torus, 3ball, 4half-ball-bot, 5half-ball-top) axis:%d (1:u v:2 uv:3) tol_ambig:%f\n",nurbsType,axis,tol_ambig);
2461 #endif
2462
2463 sem_wait(&sem_g);
2464 nurbs[nr].nurbsType=nurbsType;
2465 sem_post(&sem_g);
2466
2467 finerProjection:;
2468 pcollapsed=-1;
2469 if(surface) delete surface; // Release surface object.
2470
2471 // (1) Create an array of control points from your cgx control points:
2472 controlPoints = new snlCtrlPoint [ numberOfControlPoints ];
2473 index=0;
2474
2475 sem_wait(&sem_g);
2476
2477 for (i=0; i<nurbs[nr].u_npnt; i++)
2478 {
2479 for (j=0; j<nurbs[nr].v_npnt; j++)
2480 {
2481 controlPoints [ index++ ].components (
2482 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+0],
2483 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+1],
2484 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+2],
2485 nurbs[nr].ctlarray[i*(nurbs[nr].v_stride*nurbs[nr].v_npnt)+j*(nurbs[nr].v_stride)+3] );
2486 }
2487 }
2488
2489 // (2) Create an array of knots for each knot vector in your cgx surface:
2490 knotsU = new knot [ nurbs[nr].u_nknt ];
2491 knotsV = new knot [ nurbs[nr].v_nknt ];
2492
2493 for(i=0; i<nurbs[nr].u_nknt; i++) knotsU[i]=nurbs[nr].uknt[i];
2494 for(i=0; i<nurbs[nr].v_nknt; i++) knotsV[i]=nurbs[nr].vknt[i];
2495
2496 // (3) Create an snlSurface:
2497 // You need to supply degreeU, degreeV, sizeU and sizeV.
2498 #if TEST
2499 printf("proj attempt\n");
2500 #endif
2501
2502 surface = new snlSurface ( nurbs[nr].u_exp, nurbs[nr].v_exp, nurbs[nr].u_npnt, nurbs[nr].v_npnt, controlPoints, knotsU, knotsV );
2503
2504 sem_post(&sem_g);
2505
2506 for(i=0; i<sum_p; i++) pntproj[i].n=0;
2507 inverted = surface -> fastProject( toProject, sum_p, &sum_inverted, convergTol, normTol, maxPass, PROJ_SENSITIVITY, returnedPntsPerPnt );
2508
2509 // sort all projections according to the point indexes in "toProject"
2510 for(j=0; j<sum_inverted; j++)
2511 {
2512 if(pntproj[inverted[j].origPtIndex].n>=MAX_RETURNED_PNTS_PER_PNT) break;
2513 pntproj[inverted[j].origPtIndex].paramU[pntproj[inverted[j].origPtIndex].n]
2514 =inverted[j].paramU;
2515 pntproj[inverted[j].origPtIndex].paramV[pntproj[inverted[j].origPtIndex].n]
2516 =inverted[j].paramV;
2517 pntproj[inverted[j].origPtIndex].dist[pntproj[inverted[j].origPtIndex].n]
2518 =inverted[j].dist;
2519 pntproj[inverted[j].origPtIndex].n++;
2520 }
2521
2522 tryThePlate:;
2523 // plate
2524 if(nurbsType==0)
2525 {
2526 n=0;
2527 sem_wait(&sem_g);
2528 nbuf[0]=nurbs[nr].nc[patch];
2529 sem_post(&sem_g);
2530 for(curve=0; curve<nbuf[0]; curve++)
2531 {
2532 j=0;
2533 // over all points of that curve
2534 // scip the last point to avoid differences between 1st and last
2535 sem_wait(&sem_g);
2536 nbuf[1]=nurbs[nr].np[patch][curve]-1;
2537 sem_post(&sem_g);
2538 for(k=0; k<nbuf[1]; k++)
2539 {
2540 distp=MAXVALUE;
2541 for(i=0; i<pntproj[n].n; i++)
2542 {
2543 if(pntproj[n].dist[i]<distp)
2544 { distp=pntproj[n].dist[i]; disti=i; }
2545 }
2546 // check if the distance is "close"
2547 if((distp>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
2548 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; goto finerProjection; }
2549
2550 sem_wait(&sem_g);
2551 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramU[disti];
2552 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramV[disti];
2553 sem_post(&sem_g);
2554
2555 #if TEST
2556 printf("plate pnt:%d dist:%f disti:%d uv:%e %e orig-xyz: %f %f %f \n",n, distp,disti, pntproj[n].paramU[disti],pntproj[n].paramV[disti], toProject[n].x(), toProject[n].y(), toProject[n].z() );
2557 #endif
2558 n++;
2559 }
2560
2561 /* last point == 1st point */
2562 sem_wait(&sem_g);
2563 nurbs[nr].uv[patch][curve][j++]=nurbs[nr].uv[patch][curve][0];
2564 nurbs[nr].uv[patch][curve][j++]=nurbs[nr].uv[patch][curve][1];
2565 sem_post(&sem_g);
2566 n++;
2567 }
2568 }
2569 else
2570 {
2571 // go to the first not ambiguous point and store the skipped ones
2572 // go over all points and chose the closest one if ambiguous
2573 // finally care about the skipped leading points
2574
2575 // if ball look if the max or min direction of u or v is collapsed
2576 if(nurbsType>2)
2577 {
2578 if(axis==1) // u is axis
2579 {
2580 if(nurbsType==3) collapsedFlag=3; // collapsed at umin and umax
2581 if(nurbsType==4) collapsedFlag=1; // collapsed at umin
2582 if(nurbsType==5) collapsedFlag=2; // collapsed at umax
2583 }
2584 if(axis==2) // v is axis
2585 {
2586 if(nurbsType==3) collapsedFlag=3; // collapsed at vmin and vmax
2587 if(nurbsType==4) collapsedFlag=1; // collapsed at vmin
2588 if(nurbsType==5) collapsedFlag=2; // collapsed at vmax
2589 }
2590 if(axis==3)
2591 {
2592 collapsedFlag=3;
2593 }
2594 #if TEST
2595 printf(" axis:%d collapsedFlag:%d\n", axis, collapsedFlag);
2596 #endif
2597 }
2598
2599 n=0;
2600 sem_wait(&sem_g);
2601 nbuf[0]=nurbs[nr].nc[patch];
2602 sem_post(&sem_g);
2603 for(curve=0; curve<nbuf[0]; curve++)
2604 {
2605 nppc[curve]=0;
2606 startFlag=1;
2607 pskip=0;
2608 j=0;
2609
2610 // over all points of that curve
2611 // scip the last point to avoid differences between 1st and last
2612 sem_wait(&sem_g);
2613 nbuf[1]=nurbs[nr].np[patch][curve]-1;
2614 sem_post(&sem_g);
2615 for(k=0; k<nbuf[1]; k++)
2616 {
2617 // search the closest one to the nurbs (we have several possible points per xyz-location)
2618 distp=MAXVALUE;
2619 for(i=0; i<pntproj[n].n; i++) if(pntproj[n].dist[i]<distp)
2620 {
2621 distp=pntproj[n].dist[i]; disti=i;
2622 }
2623
2624 // check if the distance is "close"
2625 if((distp>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
2626 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; goto finerProjection; }
2627 #if TEST
2628 printf("xx pnt:%d dist:%f disti:%d uv:%f %f orig-xyz: %f %f %f \n",n, distp,disti, pntproj[n].paramU[disti],pntproj[n].paramV[disti], toProject[n].x(), toProject[n].y(), toProject[n].z() );
2629 #endif
2630
2631 // look if one edge is collapsed and set coords of close points to ideal values
2632 if((nurbsType>2)&&(collapsedFlag))
2633 {
2634 if(axis==1)
2635 {
2636 if(((collapsedFlag==1)||(collapsedFlag==3))&&((pntproj[n].paramU[disti]- utol_ambig )<umin)) { pntproj[n].paramV[disti]=vmin; }
2637 if((collapsedFlag>1)&&((pntproj[n].paramU[disti]+ utol_ambig )>umax)) { pntproj[n].paramV[disti]=vmin; }
2638 }
2639 else if(axis==2)
2640 {
2641 if(((collapsedFlag==1)||(collapsedFlag==3))&&((pntproj[n].paramV[disti]- vtol_ambig )<vmin)) { pntproj[n].paramU[disti]=umin; }
2642 if((collapsedFlag>1)&&((pntproj[n].paramV[disti]+ vtol_ambig )>vmax)) { pntproj[n].paramU[disti]=umin; }
2643 }
2644 else if(axis==3)
2645 {
2646 if((pntproj[n].paramV[disti]- vtol_ambig )<vmin) { pntproj[n].paramU[disti]=umin; }
2647 if((pntproj[n].paramV[disti]+ vtol_ambig )>vmax) { pntproj[n].paramU[disti]=umin; }
2648 if((pntproj[n].paramU[disti]- utol_ambig )<umin) { pntproj[n].paramV[disti]=vmin; }
2649 if((pntproj[n].paramU[disti]+ utol_ambig )>umax) { pntproj[n].paramV[disti]=vmin; }
2650 }
2651 }
2652
2653 // look for leading ambiguous points (at the poles) and skip the leading ones but only until the first valid was found
2654 if(startFlag)
2655 {
2656 if((axis==2)||(axis==3)) // in udir
2657 {
2658 if(((pntproj[n].paramU[disti]- utol_ambig )<umin)||((pntproj[n].paramU[disti]+ utol_ambig )>umax))
2659 {
2660 pskip++;
2661 #if TEST
2662 printf("pnt:%d scipped\n",n);
2663 #endif
2664 goto skip;
2665 }
2666 }
2667 if((axis==1)||(axis==3)) // in vdir
2668 {
2669 if(((pntproj[n].paramV[disti]- vtol_ambig )<vmin)||((pntproj[n].paramV[disti]+ vtol_ambig )>vmax))
2670 {
2671 pskip++;
2672 #if TEST
2673 printf("pnt:%d scipped\n",n);
2674 #endif
2675 goto skip;
2676 }
2677 }
2678 }
2679 startFlag=0;
2680
2681 if(!sectorFlag)
2682 {
2683 // if the point is ambig look which coord would give the shortest dist to the prev point
2684 if((axis==2)||(axis==3)) // in udir
2685 {
2686 if(((pntproj[n].paramU[disti]-utol_ambig)<umin)||((pntproj[n].paramU[disti]+utol_ambig)>umax))
2687 {
2688 sem_wait(&sem_g);
2689 dist1=nurbs[nr].uv[patch][curve][j-2]-umin;
2690 dist2=nurbs[nr].uv[patch][curve][j-2]-umax;
2691 sem_post(&sem_g);
2692 if((dist1*dist1)<(dist2*dist2))
2693 {
2694 /* pntproj is at umin. Take the value of umin if pntproj is close to umax, else do nothing */
2695 if((pntproj[n].paramU[disti]+utol_ambig)>umax) pntproj[n].paramU[disti]=umin;
2696 }
2697 else
2698 {
2699 /* pntproj is at umax. Take the value of umax if pntproj is close to umin, else do nothing */
2700 if((pntproj[n].paramU[disti]-utol_ambig)<umin) pntproj[n].paramU[disti]=umax;
2701 }
2702 #if TEST
2703 printf("dist1:%f dist2:%f corr.u:%f\n", dist1, dist2, pntproj[n].paramU[disti]);
2704 #endif
2705 }
2706 }
2707 if((axis==1)||(axis==3)) // in vdir
2708 {
2709 if(((pntproj[n].paramV[disti]-vtol_ambig)<vmin)||((pntproj[n].paramV[disti]+vtol_ambig)>vmax))
2710 {
2711 sem_wait(&sem_g);
2712 dist1=nurbs[nr].uv[patch][curve][j-1]-vmin;
2713 dist2=nurbs[nr].uv[patch][curve][j-1]-vmax;
2714 sem_post(&sem_g);
2715 if((dist1*dist1)<(dist2*dist2))
2716 {
2717 /* pntproj is at vmin. Take the value of vmin if pntproj is close to vmax, else do nothing */
2718 if((pntproj[n].paramV[disti]+vtol_ambig)>vmax) pntproj[n].paramV[disti]=vmin;
2719 }
2720 else
2721 {
2722 /* pntproj is at vmax. Take the value of vmax if pntproj is close to vmin, else do nothing */
2723 if((pntproj[n].paramV[disti]-vtol_ambig)<vmin) pntproj[n].paramV[disti]=vmax;
2724 }
2725 #if TEST
2726 printf("dist1:%f dist2:%f corr.v:%f\n", dist1, dist2, pntproj[n].paramV[disti]);
2727 #endif
2728 }
2729 }
2730 }
2731
2732 // if the coordinates of this point have to be used for the previously duplicated one
2733 #if TEST
2734 if((nextUflag>0)||(nextVflag>0))
2735 printf(" nextu:%d nextv:%d prev pnt:%d changed to uv:%f %f \n",nextUflag,nextVflag,n-1, pntproj[n].paramU[disti],pntproj[n].paramV[disti] );
2736 #endif
2737 sem_wait(&sem_g);
2738 if(nextUflag>0)
2739 { nurbs[nr].uv[patch][curve][nextUflag]=pntproj[n].paramU[disti]; nextUflag=0; }
2740 if(nextVflag>0)
2741 { nurbs[nr].uv[patch][curve][nextVflag]=pntproj[n].paramV[disti]; nextVflag=0; }
2742 sem_post(&sem_g);
2743
2744 // if collapsed edges exist add an additional point
2745 if(nurbsType>2)
2746 {
2747 if(axis==1) // in udir
2748 {
2749 if(nurbsType==3)
2750 { if(((pntproj[n].paramU[disti]-utol_ambig)<umin)||((pntproj[n].paramU[disti]+utol_ambig)>umax)) buf=1; else buf=0; }
2751 if(nurbsType==4)
2752 { if((pntproj[n].paramU[disti]-utol_ambig)<umin) buf=1; else buf=0; }
2753 if(nurbsType==5)
2754 { if((pntproj[n].paramU[disti]+utol_ambig)>umax) buf=1; else buf=0; }
2755 if(buf) // on the collapsed edge
2756 {
2757 // additional point required
2758 // use the v value from the previous point
2759 sem_wait(&sem_g);
2760 #if TEST
2761 if(j>2) printf("use v-val from prev point:%f\n", nurbs[nr].uv[patch][curve][j-1] );
2762 else printf("use v-val from prev point: (wait till end)\n" );
2763 #endif
2764 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramU[disti];
2765 if(j>2) { nurbs[nr].uv[patch][curve][j]=nurbs[nr].uv[patch][curve][j-2];j++; }
2766 else lastVflag=j++;
2767
2768 // include a new point
2769 nppc[curve]++;
2770 if(nppc[curve]>NURS_ADD_AMBIG_PNTS)
2771 {
2772 if( (nurbs[nr].uv[patch][curve] = (GLfloat *)realloc( (GLfloat *)nurbs[nr].uv[patch][curve], ((nurbs[nr].np[patch][curve]+1+nppc[curve])*2)*sizeof(GLfloat) ))==NULL )
2773 { printf(" ERROR: realloc failure in trimNurbs(), nurbs:%s can not be shaped\n\n", nurbs[nr].name); }
2774 }
2775 // use the v value from the next point
2776 // the value is not really known now, save a flag for later correction
2777 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramU[disti];
2778 distp=MAXVALUE; for(i=0; i<pntproj[n+1].n; i++) if(pntproj[n+1].dist[i]<distp)
2779 { distp=pntproj[n+1].dist[i]; nurbs[nr].uv[patch][curve][j]=pntproj[n+1].paramV[i]; }
2780 sem_post(&sem_g);
2781
2782 // check if the distance is "close"
2783 if((distp>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
2784 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; goto finerProjection; }
2785 nextVflag=j++;
2786 pcollapsed=j/2;
2787 #if TEST
2788 printf("collapsed at u:%f add pnt:%d\n", pntproj[n].paramU[disti],pcollapsed );
2789 #endif
2790 }
2791 else // not on the collapsed edge
2792 {
2793 sem_wait(&sem_g);
2794 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramU[disti];
2795 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramV[disti];
2796 #if TEST
2797 printf(" ori:%d uv:%f %f\n", n, nurbs[nr].uv[patch][curve][j-2], nurbs[nr].uv[patch][curve][j-1] );
2798 #endif
2799 sem_post(&sem_g);
2800 }
2801 }
2802 else // in vdir
2803 {
2804 if(nurbsType==3)
2805 { if(((pntproj[n].paramV[disti]-vtol_ambig)<vmin)||((pntproj[n].paramV[disti]+vtol_ambig)>vmax)) buf=1; else buf=0; }
2806 if(nurbsType==4)
2807 { if((pntproj[n].paramV[disti]-vtol_ambig)<vmin) buf=1; else buf=0; }
2808 if(nurbsType==5)
2809 { if((pntproj[n].paramV[disti]+vtol_ambig)>vmax) buf=1; else buf=0; }
2810 if(buf) // on the collapsed edge
2811 {
2812 // additional point required
2813 // use the u value from the previous point
2814 sem_wait(&sem_g);
2815 #if TEST
2816 if(j>2) printf("use u-val from prev point:%f\n", nurbs[nr].uv[patch][curve][j-2] );
2817 else printf("use u-val from prev point: (wait till end)\n" );
2818 #endif
2819 if(j>2) { nurbs[nr].uv[patch][curve][j]=nurbs[nr].uv[patch][curve][j-2];j++; }
2820 else lastUflag=j++;
2821 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramV[disti];
2822
2823 // include a new point
2824 nppc[curve]++;
2825 if(nppc[curve]>NURS_ADD_AMBIG_PNTS)
2826 {
2827 if( (nurbs[nr].uv[patch][curve] = (GLfloat *)realloc( (GLfloat *)nurbs[nr].uv[patch][curve], ((nurbs[nr].np[patch][curve]+1+nppc[curve])*2)*sizeof(GLfloat) ))==NULL )
2828 { printf(" ERROR: realloc failure in trimNurbs(), nurbs:%s can not be shaped\n\n", nurbs[nr].name); }
2829 }
2830 // use the u value from the next point
2831 // the value is not really known now, save a flag for later correction
2832 distp=MAXVALUE; for(i=0; i<pntproj[n+1].n; i++) if(pntproj[n+1].dist[i]<distp)
2833 { distp=pntproj[n+1].dist[i]; nurbs[nr].uv[patch][curve][j]=pntproj[n+1].paramU[i]; }
2834
2835 // check if the distance is "close"
2836 if((distp>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
2837 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; sem_post(&sem_g); goto finerProjection; }
2838 nextUflag=j++;
2839 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramV[disti];
2840 sem_post(&sem_g);
2841 pcollapsed=j/2;
2842 #if TEST
2843 printf("collapsed at v:%f add pnt:%d\n", pntproj[n].paramV[disti], pcollapsed );
2844 #endif
2845 }
2846 else // not on the collapsed edge
2847 {
2848 sem_wait(&sem_g);
2849 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramU[disti];
2850 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramV[disti];
2851 #if TEST
2852 printf(" ori:%d uv:%f %f\n", n, nurbs[nr].uv[patch][curve][j-2], nurbs[nr].uv[patch][curve][j-1] );
2853 #endif
2854 sem_post(&sem_g);
2855 }
2856 }
2857 }
2858 // no collapsed edges
2859 else
2860 {
2861 sem_wait(&sem_g);
2862 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramU[disti];
2863 nurbs[nr].uv[patch][curve][j++]=pntproj[n].paramV[disti];
2864 sem_post(&sem_g);
2865 }
2866
2867 skip:;
2868 n++;
2869 if(n>=sum_p) { printf("ERROR in trimNurbs, nr of points do not match. Talk to the programmer\n"); break; }
2870 }
2871 if(j<2) { nurbsType=0; goto tryThePlate; } //return(1); // no regular point was found, all points are on a collapsed edge
2872
2873 #if TEST
2874 printf("\nLEADING ambiguous points:%d, not ambiguous:%d\n\n",pskip, j/2);
2875 #endif
2876
2877 // store leading ambiguous points
2878 sem_wait(&sem_g);
2879 nbuf[1]=n-nurbs[nr].np[patch][curve]+1;
2880 nbuf[2]=n-nurbs[nr].np[patch][curve]+1+pskip;
2881 sem_post(&sem_g);
2882 for(k=nbuf[1]; k<nbuf[2]; k++)
2883 {
2884 // search the closest one to the nurbs
2885 distp=MAXVALUE;
2886 for(i=0; i<pntproj[k].n; i++) if(pntproj[k].dist[i]<distp)
2887 { distp=pntproj[k].dist[i]; disti=i; }
2888
2889 // check if the distance is "close"
2890 if((distp>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
2891 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; goto finerProjection; }
2892 #if TEST
2893 printf("3 pnt:%d dist:%f disti:%d uv:%f %f orig-xyz: %f %f %f \n", k, distp,disti, pntproj[k].paramU[disti],pntproj[k].paramV[disti], toProject[k].x(), toProject[k].y(), toProject[k].z() );
2894 #endif
2895
2896
2897 // look which coord would give the shortest dist to the last-1 point of the curve
2898 if((axis==2)||(axis==3)) // in udir
2899 {
2900 if(((pntproj[k].paramU[disti]-utol_ambig)<umin)||((pntproj[k].paramU[disti]+utol_ambig)>umax))
2901 {
2902 sem_wait(&sem_g);
2903 dist1=nurbs[nr].uv[patch][curve][j-2]-umin;
2904 dist2=nurbs[nr].uv[patch][curve][j-2]-umax;
2905 sem_post(&sem_g);
2906 if((dist1*dist1)<(dist2*dist2))
2907 {
2908 /* pntproj is at umin. Take the value of umin if pntproj is close to umax, else do nothing */
2909 if((pntproj[k].paramU[disti]+utol_ambig)>umax) pntproj[k].paramU[disti]=umin;
2910 }
2911 else
2912 {
2913 /* pntproj is at umax. Take the value of umax if pntproj is close to umin, else do nothing */
2914 if((pntproj[k].paramU[disti]-utol_ambig)<umin) pntproj[k].paramU[disti]=umax;
2915 }
2916 #if TEST
2917 printf("dist1:%f dist2:%f corr.u:%f\n", dist1, dist2, pntproj[k].paramU[disti]);
2918 #endif
2919 }
2920 }
2921 if((axis==1)||(axis==3)) // in vdir
2922 {
2923 if(((pntproj[k].paramV[disti]-vtol_ambig)<vmin)||((pntproj[k].paramV[disti]+vtol_ambig)>vmax))
2924 {
2925 sem_wait(&sem_g);
2926 dist1=nurbs[nr].uv[patch][curve][j-1]-vmin;
2927 dist2=nurbs[nr].uv[patch][curve][j-1]-vmax;
2928 sem_post(&sem_g);
2929 if((dist1*dist1)<(dist2*dist2))
2930 {
2931 /* pntproj is at vmin. Take the value of vmin if pntproj is close to vmax, else do nothing */
2932 if((pntproj[k].paramV[disti]+vtol_ambig)>vmax) pntproj[k].paramV[disti]=vmin;
2933 }
2934 else
2935 {
2936 /* pntproj is at vmax. Take the value of vmax if pntproj is close to vmin, else do nothing */
2937 if((pntproj[k].paramV[disti]-vtol_ambig)<vmin) pntproj[k].paramV[disti]=vmax;
2938 }
2939 #if TEST
2940 printf("dist1:%f dist2:%f corr.v:%f\n", dist1, dist2, pntproj[k].paramV[disti]);
2941 #endif
2942 }
2943 }
2944
2945
2946 // if the coordinates of this point have to be used for the previously duplicated one
2947 sem_wait(&sem_g);
2948 if(nextUflag>0)
2949 { nurbs[nr].uv[patch][curve][nextUflag]=pntproj[k].paramU[disti]; nextUflag=0; }
2950 if(nextVflag>0)
2951 { nurbs[nr].uv[patch][curve][nextVflag]=pntproj[k].paramV[disti]; nextVflag=0; }
2952 sem_post(&sem_g);
2953
2954 // if collapsed edges exist add an additional point
2955 if(nurbsType>2)
2956 {
2957 if(axis==1) // in udir
2958 {
2959 if(nurbsType==3)
2960 { if(((pntproj[k].paramU[disti]-utol_ambig)<umin)||((pntproj[k].paramU[disti]+utol_ambig)>umax)) buf=1; else buf=0; }
2961 if(nurbsType==4)
2962 { if((pntproj[k].paramU[disti]-utol_ambig)<umin) buf=1; else buf=0; }
2963 if(nurbsType==5)
2964 { if((pntproj[k].paramU[disti]+utol_ambig)>umax) buf=1; else buf=0; }
2965 if(buf) // on the collapsed edge
2966 {
2967 // additional point required
2968 // use the v value from the previous point
2969 sem_wait(&sem_g);
2970 #if TEST
2971 if(j>2) printf("use v-val from prev point:%f\n", nurbs[nr].uv[patch][curve][j-1] );
2972 else printf("use v-val from prev point: (wait till end)\n" );
2973 #endif
2974 nurbs[nr].uv[patch][curve][j++]=pntproj[k].paramU[disti];
2975 //nurbs[nr].uv[patch][curve][j]=nurbs[nr].uv[patch][curve][j-2]; j++;
2976 if(j>2) { nurbs[nr].uv[patch][curve][j]=nurbs[nr].uv[patch][curve][j-2];j++; }
2977 else lastVflag=j++;
2978
2979 // include a new point
2980 nppc[curve]++;
2981 if(nppc[curve]>NURS_ADD_AMBIG_PNTS)
2982 {
2983 if( (nurbs[nr].uv[patch][curve] = (GLfloat *)realloc( (GLfloat *)nurbs[nr].uv[patch][curve], ((nurbs[nr].np[patch][curve]+1+nppc[curve])*2)*sizeof(GLfloat) ))==NULL )
2984 { printf(" ERROR: realloc failure in trimNurbs(), nurbs:%s can not be shaped\n\n", nurbs[nr].name); }
2985 }
2986 // use the v value from the next point
2987 // the value is not really known now, save a flag for later correction
2988 nurbs[nr].uv[patch][curve][j++]=pntproj[k].paramU[disti];
2989 distp=MAXVALUE; for(i=0; i<pntproj[k+1].n; i++) if(pntproj[k+1].dist[i]<distp)
2990 { distp=pntproj[k+1].dist[i]; nurbs[nr].uv[patch][curve][j]=pntproj[k+1].paramV[i]; }
2991 sem_post(&sem_g);
2992
2993 // check if the distance is "close"
2994 if((distp>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
2995 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; goto finerProjection; }
2996 nextVflag=j++;
2997 pcollapsed=j/2;
2998 #if TEST
2999 printf("collapsed at u:%f add pnt:%d\n", pntproj[k].paramU[disti],pcollapsed );
3000 #endif
3001 }
3002 else // not on the collapsed edge
3003 {
3004 sem_wait(&sem_g);
3005 nurbs[nr].uv[patch][curve][j++]=pntproj[k].paramU[disti];
3006 nurbs[nr].uv[patch][curve][j++]=pntproj[k].paramV[disti];
3007 #if TEST
3008 printf(" ori:%d uv:%f %f\n", k, nurbs[nr].uv[patch][curve][j-2], nurbs[nr].uv[patch][curve][j-1] );
3009 #endif
3010 sem_post(&sem_g);
3011 }
3012 }
3013 else // in vdir
3014 {
3015 if(nurbsType==3)
3016 /* this is the original line but its suspicious. changed from u to v:
3017 { if(((pntproj[k].paramU[disti]-utol_ambig)<umin)||((pntproj[k].paramU[disti]+utol_ambig)>umax)) buf=1; else buf=0; }
3018 */
3019 { if(((pntproj[k].paramV[disti]-vtol_ambig)<vmin)||((pntproj[k].paramV[disti]+vtol_ambig)>vmax)) buf=1; else buf=0; }
3020 if(nurbsType==4)
3021 { if((pntproj[k].paramV[disti]-vtol_ambig)<vmin) buf=1; else buf=0; }
3022 if(nurbsType==5)
3023 { if((pntproj[k].paramV[disti]+vtol_ambig)>vmax) buf=1; else buf=0; }
3024 if(buf) // on the collapsed edge
3025 {
3026 // additional point required
3027 // use the u value from the previous point
3028 sem_wait(&sem_g);
3029 #if TEST
3030 if(j>2) printf("1 use u-val from prev point:%f\n", nurbs[nr].uv[patch][curve][j-2] );
3031 else printf("use u-val from prev point: (wait till end)\n" );
3032 #endif
3033 if(j>2) { nurbs[nr].uv[patch][curve][j]=nurbs[nr].uv[patch][curve][j-2];j++; }
3034 else lastUflag=j++;
3035 nurbs[nr].uv[patch][curve][j++]=pntproj[k].paramV[disti];
3036
3037 // include a new point
3038 nppc[curve]++;
3039 if(nppc[curve]>NURS_ADD_AMBIG_PNTS)
3040 {
3041 if( (nurbs[nr].uv[patch][curve] = (GLfloat *)realloc( (GLfloat *)nurbs[nr].uv[patch][curve], ((nurbs[nr].np[patch][curve]+1+nppc[curve])*2)*sizeof(GLfloat) ))==NULL )
3042 { printf(" ERROR: realloc failure in trimNurbs(), nurbs:%s can not be shaped\n\n", nurbs[nr].name); }
3043 }
3044 // use the u value from the next point
3045 // the value is not really known now, save a flag for later correction
3046 distp=MAXVALUE; for(i=0; i<pntproj[k+1].n; i++) if(pntproj[k+1].dist[i]<distp)
3047 { distp=pntproj[k+1].dist[i]; nurbs[nr].uv[patch][curve][j]=pntproj[k+1].paramU[i]; }
3048
3049 // check if the distance is "close"
3050 if((distp>MIN_DIST)&&(returnedPntsPerPnt<MAX_RETURNED_PNTS_PER_PNT))
3051 { returnedPntsPerPnt=MAX_RETURNED_PNTS_PER_PNT; delete[] inverted; sem_post(&sem_g); goto finerProjection; }
3052 nextUflag=j++;
3053 nurbs[nr].uv[patch][curve][j++]=pntproj[k].paramV[disti];
3054 sem_post(&sem_g);
3055 pcollapsed=j/2;
3056 #if TEST
3057 printf("collapsed at v:%f add pnt:%d\n", pntproj[k].paramV[disti], j/2 );
3058 #endif
3059 }
3060 else // not on the collapsed edge
3061 {
3062 sem_wait(&sem_g);
3063 nurbs[nr].uv[patch][curve][j++]=pntproj[k].paramU[disti];
3064 nurbs[nr].uv[patch][curve][j++]=pntproj[k].paramV[disti];
3065 #if TEST
3066 printf(" ori:%d uv:%f %f\n", k, nurbs[nr].uv[patch][curve][j-2], nurbs[nr].uv[patch][curve][j-1] );
3067 #endif
3068 sem_post(&sem_g);
3069 }
3070 }
3071 }
3072 else
3073 {
3074 sem_wait(&sem_g);
3075 nurbs[nr].uv[patch][curve][j++]=pntproj[k].paramU[disti];
3076 nurbs[nr].uv[patch][curve][j++]=pntproj[k].paramV[disti];
3077 #if TEST
3078 printf("add pnt:%d uv:%f %f %d\n", j/2, pntproj[k].paramU[disti],pntproj[k].paramV[disti], nurbs[nr].np[patch][curve]);
3079 #endif
3080 sem_post(&sem_g);
3081 }
3082 }
3083
3084 // if the coordinates of this point have to be used
3085 sem_wait(&sem_g);
3086 if(lastUflag>-1)
3087 { nurbs[nr].uv[patch][curve][lastUflag]=nurbs[nr].uv[patch][curve][j-2]; lastUflag=-1; }
3088 if(lastVflag>-1)
3089 { nurbs[nr].uv[patch][curve][lastVflag]=nurbs[nr].uv[patch][curve][j-1]; lastVflag=-1; }
3090 sem_post(&sem_g);
3091
3092 // if a collapsed point exists, start the loop at that point
3093 if(pcollapsed!=-1)
3094 {
3095 sem_wait(&sem_g);
3096 nurbs[nr].sum_ambiguousPnts[patch][curve]=nppc[curve];
3097 nn=0;
3098 double* uvbuf = new double [ (nurbs[nr].np[patch][curve]+nppc[curve])*2 +1 ];
3099
3100 #if TEST
3101 i=0;
3102 for(k=0; k<nurbs[nr].np[patch][curve]+nppc[curve]-1; k++)
3103 {
3104 printf(" pnt P%d %f",k, nurbs[nr].uv[patch][curve][i++] );
3105 printf(" %f\n", nurbs[nr].uv[patch][curve][i++] );
3106 }
3107 #endif
3108
3109 i=pcollapsed*2-2;
3110 for(k=0; k<nurbs[nr].np[patch][curve]+nppc[curve]-pcollapsed; k++)
3111 {
3112 //printf(" pnt P%d %f",k, nurbs[nr].uv[patch][curve][i++] );
3113 //printf(" %f\n", nurbs[nr].uv[patch][curve][i++] ); i--; i--;
3114 uvbuf[nn++]=nurbs[nr].uv[patch][curve][i++];
3115 uvbuf[nn++]=nurbs[nr].uv[patch][curve][i++];
3116 }
3117 i=0;
3118 for(k=0; k<pcollapsed-1; k++)
3119 {
3120 //printf(" pnt ! %f", nurbs[nr].uv[patch][curve][i++] );
3121 //printf(" %f\n", nurbs[nr].uv[patch][curve][i++] );
3122 uvbuf[nn++]=nurbs[nr].uv[patch][curve][i++];
3123 uvbuf[nn++]=nurbs[nr].uv[patch][curve][i++];
3124 }
3125 i=0;
3126 for(k=0; k<nurbs[nr].np[patch][curve]+nppc[curve]-1; k++)
3127 {
3128 //printf(" pnt ! %f", nurbs[nr].uv[patch][curve][i++] );
3129 //printf(" %f\n", nurbs[nr].uv[patch][curve][i++] );
3130 nurbs[nr].uv[patch][curve][i]=uvbuf[i]; i++;
3131 nurbs[nr].uv[patch][curve][i]=uvbuf[i]; i++;
3132 }
3133 sem_post(&sem_g);
3134 delete[] uvbuf;
3135 }
3136
3137 // else if points were skipped, start the loop at the fist skipped point
3138 else if(pskip)
3139 {
3140 sem_wait(&sem_g);
3141 double* uvbuf = new double [ (nurbs[nr].np[patch][curve])*2 +1 ];
3142 i=0;
3143 for(k=0; k<nurbs[nr].np[patch][curve]-1; k++)
3144 {
3145 //printf(" pnt P%d %f",k, nurbs[nr].uv[patch][curve][i++] );
3146 //printf(" %f\n", nurbs[nr].uv[patch][curve][i++] );
3147 uvbuf[i]=nurbs[nr].uv[patch][curve][i]; i++;
3148 uvbuf[i]=nurbs[nr].uv[patch][curve][i]; i++;
3149 }
3150
3151 i=0;
3152 nn=(nurbs[nr].np[patch][curve]-1-pskip)*2;
3153 for(k=0; k<pskip; k++)
3154 {
3155 //printf(" pnt ! %f", nurbs[nr].uv[patch][curve][i++] );
3156 //printf(" %f\n", nurbs[nr].uv[patch][curve][i++] );
3157 nurbs[nr].uv[patch][curve][i++]=uvbuf[nn++];
3158 nurbs[nr].uv[patch][curve][i++]=uvbuf[nn++];
3159 }
3160 nn=0;
3161 for(k=0; k<nurbs[nr].np[patch][curve]-pskip-1; k++)
3162 {
3163 //printf(" pnt ! %f", nurbs[nr].uv[patch][curve][i++] );
3164 //printf(" %f\n", nurbs[nr].uv[patch][curve][i++] );
3165 nurbs[nr].uv[patch][curve][i++]=uvbuf[nn++];
3166 nurbs[nr].uv[patch][curve][i++]=uvbuf[nn++];
3167 }
3168 sem_post(&sem_g);
3169
3170 delete[] uvbuf;
3171 }
3172
3173 // last point == 1st point
3174 sem_wait(&sem_g);
3175 nurbs[nr].uv[patch][curve][j++]=nurbs[nr].uv[patch][curve][0];
3176 nurbs[nr].uv[patch][curve][j++]=nurbs[nr].uv[patch][curve][1];
3177
3178 // if the coordinates of this point have to be used for the previously duplicated one
3179 if(nextUflag>0)
3180 { nurbs[nr].uv[patch][curve][nextUflag]=nurbs[nr].uv[patch][curve][j-2]; }
3181 if(nextVflag>0)
3182 { nurbs[nr].uv[patch][curve][nextVflag]=nurbs[nr].uv[patch][curve][j-1]; }
3183
3184 #if TEST
3185 // check if we have too many points now
3186 if(((nurbs[nr].np[patch][curve]+nppc[curve])*2)!=j)
3187 {
3188 printf(" ERROR too many points:%d/2 (!=%d) in curve)\n", j,nurbs[nr].np[patch][curve]);
3189 exit(-1);
3190 }
3191 #endif
3192 sem_post(&sem_g);
3193 n++;
3194 }
3195
3196 // The "ball" type is difficult if the trimming loop extend in the ambiguous zone. To avoid this the nurbs is rotated in space.
3197 if((nurbsType==3)&&(moveFlag))
3198 {
3199 moveFlag=0;
3200 if( rotateBall( nr, axis, patch, utol_ambig, vtol_ambig) ) goto finerProjection;
3201 }
3202
3203 }
3204
3205 sem_wait(&sem_g);
3206 #if FAST_SHADING
3207 // for the moment use the coarsest resolution:
3208 nurbs[nr].ustep[patch]=nurbs[nr].u_exp*4./(umax-umin);
3209 nurbs[nr].vstep[patch]=nurbs[nr].v_exp*4./(vmax-vmin);
3210 #else
3211 /* calc the average local spacing between the trimming points (only outer loop) */
3212 av_local_elength=0.;
3213 i=0;
3214 p0[0]=nurbs[nr].xyz[patch][0][i++];
3215 p0[1]=nurbs[nr].xyz[patch][0][i++];
3216 p0[2]=nurbs[nr].xyz[patch][0][i++];
3217 for(j=0; j<nurbs[nr].np[patch][0]-1; j++)
3218 {
3219 p1[0]=nurbs[nr].xyz[patch][0][i++];
3220 p1[1]=nurbs[nr].xyz[patch][0][i++];
3221 p1[2]=nurbs[nr].xyz[patch][0][i++];
3222 v_result(p0,p1,p0p1);
3223 av_local_elength+=v_betrag(p0p1);
3224 //printf("pnt ! %f %f %f av_local_elength %f\n", p1[0],p1[1],p1[2], av_local_elength);
3225 p0[0]=p1[0];
3226 p0[1]=p1[1];
3227 p0[2]=p1[2];
3228 }
3229 av_local_elength/=nurbs[nr].np[patch][0]-1;
3230
3231 // adapt the nurbs-parameter; number of points (np)
3232 for(curve=0; curve<nurbs[nr].nc[patch]; curve++) nurbs[nr].np[patch][curve]+=nppc[curve];
3233 #if TEST
3234 i=0;
3235 for(curve=0; curve<nurbs[nr].nc[patch]; curve++)
3236 for(j=0; j<nurbs[nr].np[patch][curve]; j++)
3237 {
3238 printf(" pnt ! %e", nurbs[nr].uv[patch][curve][i++] );
3239 printf(" %e\n", nurbs[nr].uv[patch][curve][i++] );
3240 }
3241 #endif
3242
3243 // average resolution based on gtol which is related to the width of the model and was used in calcNurbsResolution()
3244 nurbs[nr].ustep[patch]=nurbs[nr].ures/av_local_elength/4;
3245 nurbs[nr].vstep[patch]=nurbs[nr].vres/av_local_elength/4;
3246 //printf("i(pnts*3):%d uvres %f %f elength %f uvstep %f %f\n",i, nurbs[nr].ures,nurbs[nr].vres,av_local_elength,nurbs[nr].ustep[patch],nurbs[nr].vstep[patch]);
3247
3248 // avoid bad values
3249 if(nurbs[nr].ustep[patch]<(nurbs[nr].u_exp*2./(umax-umin))) nurbs[nr].ustep[patch]=nurbs[nr].u_exp*2./(umax-umin);
3250 if(nurbs[nr].vstep[patch]<(nurbs[nr].u_exp*2./(vmax-vmin))) nurbs[nr].vstep[patch]=nurbs[nr].v_exp*2./(vmax-vmin);
3251 if(nurbs[nr].ustep[patch]>CGX_GLU_MAX_STEPS) nurbs[nr].ustep[patch]=CGX_GLU_MAX_STEPS;
3252 if(nurbs[nr].vstep[patch]>CGX_GLU_MAX_STEPS) nurbs[nr].vstep[patch]=CGX_GLU_MAX_STEPS;
3253
3254 if(nurbs[nr].ustep[patch]<nurbs[nr].u_exp*4./(umax-umin)) nurbs[nr].ustep[patch]=nurbs[nr].u_exp*4./(umax-umin);
3255 if(nurbs[nr].vstep[patch]<nurbs[nr].v_exp*4./(vmax-vmin)) nurbs[nr].vstep[patch]=nurbs[nr].v_exp*4./(vmax-vmin);
3256
3257 #endif
3258 sem_post(&sem_g);
3259
3260 delete[] pntproj;
3261 delete[] nppc; // Release number of points per curve object
3262 delete surface; // Release surface object.
3263 delete[] inverted; // Release snlVertex array returned from projection function.
3264 delete[] toProject; // Release points that were inverted onto surface.
3265
3266 // The rest of the objects are deleted by libSNL.
3267 return(0);
3268 }
3269