ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/addin/step/src/stbsplines.cpp
Revision: 1156
Committed: Thu Jun 13 22:02:48 2024 UTC (14 months ago) by francois
File size: 21654 byte(s)
Log Message:
compatibilité Ubuntu 22.04
Suppression des refeences à Windows
Ajout d'une banière

File Contents

# User Rev Content
1 francois 1156 //####//------------------------------------------------------------
2     //####//------------------------------------------------------------
3     //####// MAGiC
4     //####// Jean Christophe Cuilliere et Vincent FRANCOIS
5     //####// Departement de Genie Mecanique - UQTR
6     //####//------------------------------------------------------------
7     //####// MAGIC est un projet de recherche de l equipe ERICCA
8     //####// du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
9     //####// http://www.uqtr.ca/ericca
10     //####// http://www.uqtr.ca/
11     //####//------------------------------------------------------------
12     //####//------------------------------------------------------------
13     //####//
14     //####// stbsplines.cpp
15     //####//
16     //####//------------------------------------------------------------
17     //####//------------------------------------------------------------
18     //####// COPYRIGHT 2000-2024
19     //####// jeu 13 jun 2024 11:53:59 EDT
20     //####//------------------------------------------------------------
21     //####//------------------------------------------------------------
22 foucault 27
23     #pragma hdrstop
24    
25    
26     #include "stbsplines.h"
27     #include <vector>
28     #include "st_gestionnaire.h"
29     #include "ot_systeme.h"
30     #include "constantegeo.h"
31    
32     #include <math.h>
33    
34     #pragma package(smart_init)
35    
36    
37     ST_B_SPLINE_SURF::ST_B_SPLINE_SURF(long LigneCourante,std::string idori,int bs_degre_u,int bs_degre_v,std::vector<int> bs_indexptsctr,std::vector<int> bs_knots_multiplicities_u,std::vector<int> bs_knots_multiplicities_v,std::vector<double> bs_knots_u,std::vector<double> bs_knots_v):ST_SURFACE(LigneCourante,idori),degre_u(bs_degre_u),degre_v(bs_degre_v),sens(1)
38     {
39 francois 283 int r_u=bs_knots_multiplicities_u.size();
40     for (int k=0;k<r_u;k++)
41     {
42 foucault 27 for (int j=0;j<bs_knots_multiplicities_u[k];j++)
43 francois 283 knots_u.insert(knots_u.end(),bs_knots_u[k]);
44     }
45     int r_v=bs_knots_multiplicities_v.size();
46     for (int k=0;k<r_v;k++)
47     {
48 foucault 27 for (int j=0;j<bs_knots_multiplicities_v[k];j++)
49 francois 283 knots_v.insert(knots_v.end(),bs_knots_v[k]);
50     }
51     nb_point=bs_indexptsctr.size();
52     nb_ptsctr_u=knots_u.size()-degre_u-1;
53     nb_ptsctr_v=knots_v.size()-degre_v-1;
54     for (int i=0;i<nb_point;i++)
55     {
56 foucault 27 indexptsctr.insert(indexptsctr.end(),bs_indexptsctr[i]);
57 francois 283 }
58     umin=knots_u[0];
59     umax=knots_u[knots_u.size()-1];
60     vmin=knots_v[0];
61     vmax=knots_v[knots_v.size()-1];
62 foucault 27 }
63    
64     ST_B_SPLINE_SURF::ST_B_SPLINE_SURF(int bs_degre_u,int bs_degre_v,std::vector<double> &bs_knots_u,std::vector<double> &bs_knots_v,std::vector<double> &bs_ptsctr,std::vector<double> &bs_poids,int sense):ST_SURFACE(),degre_u(bs_degre_u),degre_v(bs_degre_v),sens(sense)
65     {
66 francois 283 int r_u=bs_knots_u.size();
67     for (int k=0;k<r_u;k++)
68     knots_u.insert(knots_u.end(),bs_knots_u[k]);
69     int r_v=bs_knots_v.size();
70     for (int k=0;k<r_v;k++)
71 foucault 27 knots_v.insert(knots_v.end(),bs_knots_v[k]);
72 francois 283 nb_point=bs_ptsctr.size()/3;
73     nb_ptsctr_u=knots_u.size()-degre_u-1;
74     nb_ptsctr_v=knots_v.size()-degre_v-1;
75     for (int i=0;i<nb_point;i++)
76     {
77     double x=bs_ptsctr[3*i];
78     double y=bs_ptsctr[3*i+1];
79     double z=bs_ptsctr[3*i+2];
80     double w=bs_poids[i];
81     ptsctr.push_back(OT_VECTEUR_4D(w*x,w*y,w*z,w));
82     }
83     periodique_u=1;
84     for (int i=0; i<nb_ptsctr_v; i++)
85     {
86 foucault 27 double *xyz1=ptsctr[i];
87     double *xyz2=ptsctr[nb_ptsctr_v*(nb_ptsctr_u-1)+i];
88     if ((!(OPERATEUR::egal (xyz1[0],xyz2[0],1E-10))) || (!(OPERATEUR::egal (xyz1[1],xyz2[1],1E-10))) || (!(OPERATEUR::egal (xyz1[2],xyz2[2],1E-10))))
89 francois 283 {
90     periodique_u=0;
91     break;
92 foucault 27 }
93 francois 283 }
94     if (periodique_u==1)
95     {
96 foucault 27 int i=knots_u.size();
97     periode_u=(knots_u[i-1]-knots_u[0]);
98 francois 283 }
99     else periode_u=0;
100     periodique_v=1;
101     for (int j=0; j<nb_ptsctr_u; j++)
102     {
103 foucault 27 double *xyz3=ptsctr[j];
104     double *xyz4=ptsctr[(j*nb_ptsctr_v+nb_ptsctr_v-1)];
105     if ((!(OPERATEUR::egal (xyz3[0],xyz4[0],1E-10))) || (!(OPERATEUR::egal (xyz3[1],xyz4[1],1E-10))) || (!(OPERATEUR::egal (xyz3[2],xyz4[2],1E-10))))
106 francois 283 {
107     periodique_v=0;
108     break;
109 foucault 27 }
110 francois 283 }
111     if (periodique_v==1)
112     {
113 foucault 27 int j=knots_v.size();
114     periode_v=(knots_v[j-1]-knots_v[0]);
115 francois 283 }
116     else periode_v=0;
117     umin=knots_u[0];
118     umax=knots_u[knots_u.size()-1];
119     vmin=knots_v[0];
120     vmax=knots_v[knots_v.size()-1];
121 foucault 27 }
122    
123    
124     ST_B_SPLINE_SURF::~ST_B_SPLINE_SURF()
125     {
126     }
127    
128     void ST_B_SPLINE_SURF::set_periodique_u(int val)
129 francois 283 {
130     periodique_u=val;
131     if (periodique_u==1)
132     {
133 foucault 27 int i=knots_u.size();
134     periode_u=(knots_u[i-1]-knots_u[0]);
135 francois 283 }
136 foucault 27 }
137    
138     void ST_B_SPLINE_SURF::set_periodique_v(int val)
139     {
140 francois 283 periodique_v=val;
141     if (periodique_v==1)
142     {
143 foucault 27 int j=knots_v.size();
144     periode_v=(knots_v[j-1]-knots_v[0]);
145 francois 283 }
146 foucault 27 }
147    
148     int ST_B_SPLINE_SURF::get_intervalle(int nb_pt, int degre, double t, std::vector<double> &knots)
149     {
150 francois 283 int compteur = 0;
151     int inter;
152     if (OPERATEUR::egal(t,knots[nb_pt],1E-10)==1) return nb_pt-1;
153     if (OPERATEUR::egal(t,knots[degre],1E-10)==1) return degre;
154     else
155     {
156 foucault 27 int low=degre;
157     int high=nb_pt+1;
158     int mid=((low+high)/2);
159     while ((t<knots[mid-1]) || (t>=knots[mid]))
160 francois 283 {
161     if (t<knots[mid-1]) high=mid;
162     else low=mid;
163     mid=(low+high)/2;
164     compteur++;
165     }
166 foucault 27 inter=mid-1;
167 francois 283 }
168     return inter;
169 foucault 27 }
170    
171     void ST_B_SPLINE_SURF::get_valeur_fonction(int inter, double t, int degre, std::vector<double> &knots,double *grand_n)
172     {
173 francois 283 double saved;
174     grand_n[0]=1.0;
175     double *gauche=new double[degre];
176     double *droite=new double[degre];
177     for (int j=1;j<=degre;j++)
178     {
179 foucault 27 gauche[j-1]= t-knots[inter-j+1];
180     droite[j-1]=knots[inter+j]-t;
181     saved=0.0;
182     for (int r=0;r<j;r++)
183 francois 283 {
184     double temp=grand_n[r]/(droite[r]+ gauche[j-r-1]);
185     grand_n[r]=saved+droite[r]* temp;
186     saved=gauche[j-r-1]*temp;
187     }
188 foucault 27 grand_n[j]=saved;
189 francois 283 }
190     delete [] gauche;
191     delete [] droite;
192 foucault 27 }
193    
194    
195     void ST_B_SPLINE_SURF::evaluer(double *uv,double *xyz)
196     {
197    
198 francois 283 double u = uv[0];
199     double v = uv[1];
200 foucault 27
201 francois 283 if (sens==-1) u=umin+umax-u;
202 foucault 27
203     #define P(i,j) ptsctr [ (uspan - degre_u + i) * nb_ptsctr_v + (vspan - degre_v + j) ]
204    
205 francois 283 int uspan = get_intervalle(nb_ptsctr_u,degre_u,u, knots_u);
206 foucault 27
207 francois 283 int vspan = get_intervalle(nb_ptsctr_v,degre_v,v, knots_v);
208 foucault 27
209 francois 283 double Nu[20];
210 foucault 27 get_valeur_fonction(uspan,u,degre_u,knots_u,Nu);
211    
212     double Nv[20];
213     get_valeur_fonction(vspan,v,degre_v,knots_v,Nv);
214    
215 francois 283 OT_VECTEUR_4D temp[20];
216 foucault 27
217 francois 283 int l;
218     for (l=0;l<=degre_v;l++) {
219     temp[l] =0.0 ;
220     for (int k=0;k<=degre_u;k++) {
221     temp[l] += Nu[k]*P(k,l) ;
222     }
223 foucault 27 }
224 francois 283 OT_VECTEUR_4D sp(0,0,0,0) ;
225     for (l=0;l<=degre_v;l++) {
226     sp += Nv[l]*temp[l];
227     }
228 foucault 27
229 francois 283 // transform homogeneous coordinates to 3D coordinates
230     for (int i=0; i<3; i++)
231     xyz[i] = sp[i]/sp.w();
232 foucault 27
233     #undef P
234     }
235    
236     void ST_B_SPLINE_SURF::deriver_fonction(int inter,double t,int degre,int dd,std::vector<double> &knots,double *f_deriver)
237     {
238     #define f_deriver(i,j) (*(f_deriver+(i)*(degre+1)+j))
239     #define grand_n(i,j) (*(grand_n+(i)*(degre+1)+j))
240     #define a(i,j) (*(a+(i)*(dd+1)+j))
241 francois 283 double grand_n[256];
242     double saved;
243     grand_n(0,0)=1.0;
244     double gauche[16];
245     double droite[16];
246     for (int j=1;j<=degre;j++)
247     {
248 foucault 27 gauche[j]= t-knots[inter-j+1];
249     droite[j]=knots[inter+j]-t;
250     saved=0.0;
251     for (int r=0;r<j;r++)
252 francois 283 {
253     grand_n(j,r)=(droite[r+1]+ gauche[j-r]);
254     double temp = grand_n(r,j-1)/grand_n(j,r);
255 foucault 27
256 francois 283 grand_n(r,j)= saved+droite[r+1]*temp;
257     saved=gauche[j-r]*temp;
258     }
259 foucault 27 grand_n(j,j)=saved;
260 francois 283 }
261     for (int j=0;j<=degre;j++)
262     {
263 foucault 27 f_deriver(0,j)= grand_n(j,degre);
264 francois 283 }
265     double a[256];
266     for (int r=0;r<=degre;r++)
267     {
268     int s1=0;
269     int s2=1;
270 foucault 27 a(0,0)=1.0;
271     for (int k=1;k<=dd; k++)
272 francois 283 {
273     double d=0.0;
274     int rk=r-k;
275     int pk=degre-k;
276     if (r>=k)
277     {
278     a(s2,0)=a(s1,0)/grand_n(pk+1,rk);
279     d= a(s2,0)* grand_n(rk,pk);
280     }
281     int j1;
282     int j2;
283     if (rk>=-1) j1=1;
284     else j1= -rk;
285     if (r-1<=pk) j2=k-1;
286     else j2=degre-r;
287     for (int j=j1;j<=j2;j++)
288     {
289     a(s2,j) = (a(s1,j)-a(s1,j-1))/grand_n(pk+1,rk+j);
290     d+=a(s2,j)*grand_n(rk+j,pk);
291     }
292     if (r<=pk)
293     {
294     a(s2,k) = -a(s1,k-1)/grand_n(pk+1,r);
295     d+=a(s2,k)*grand_n(r,pk);
296     }
297     f_deriver(k,r)=d;
298     int j=s1;
299     s1=s2;
300     s2=j;
301 foucault 27 }
302 francois 283 }
303     int r=degre;
304     for (int k=1;k<=dd;k++)
305     {
306     for (int j=0;j<=degre;j++)
307 foucault 27 {
308 francois 283 f_deriver(k,j)*=r;
309     }
310 foucault 27 r*=(degre-k);
311 francois 283 }
312 foucault 27 #undef f_deriver
313     #undef grand_n
314 francois 283 #undef a
315 foucault 27 }
316    
317 francois 283 void ST_B_SPLINE_SURF::binomialCoef(double * Bin, int d) {
318     int n,k;
319     // Setup the first line
320     Bin[0] = 1.0 ;
321     for (k=d-1;k>0;--k)
322     Bin[k] = 0.0 ;
323     // Setup the other lines
324     for (n=0;n<d-1;n++) {
325     Bin[(n+1)*d+0] = 1.0 ;
326     for (k=1;k<d;k++)
327     if (n+1<k)
328     Bin[(n+1)*d+k] = 0.0 ;
329     else
330     Bin[(n+1)*d+k] = Bin[n*d+k] + Bin[n*d+k-1] ;
331     }
332 foucault 27 }
333    
334     void ST_B_SPLINE_SURF::deriver_bs_kieme(int nb_ptsctr_u,int degre_u,std::vector<double> &knots_u,int nb_ptsctr_v,int degre_v,std::vector<double> &knots_v,double u,double v,int d,OT_VECTEUR_4D * skl)
335     {
336     #define skl(i,j) skl[(i)*(d+1)+j]
337     #define Nu(i,j) Nu[(i)*(degre_u+1)+j]
338     #define Nv(i,j) Nv[(i)*(degre_v+1)+j]
339     #define P(i,j) ptsctr[(i)*(nb_ptsctr_v)+j]
340    
341 francois 283 int i,j,k,l,s,r,dd;
342 foucault 27
343 francois 283 int du = std::min(d,degre_u);
344     for ( k=degre_u+1;k<=d;k++)
345     {
346     for (int l=0;l<=d-k;k++)
347 foucault 27 {
348 francois 283 skl(k,l)=0.0;
349 foucault 27 }
350 francois 283 }
351     int dv = std::min(d,degre_v);
352     for ( l=degre_v+1;l<=d;l++)
353     {
354     for ( k=0;k<=d-l;k++)
355 foucault 27 {
356 francois 283 skl(k,l)=0.0;
357 foucault 27 }
358 francois 283 }
359 foucault 27
360 francois 283 int inter_u = get_intervalle(nb_ptsctr_u,degre_u,u, knots_u);
361     double Nu [ 256 ];
362     deriver_fonction(inter_u,u,degre_u,du,knots_u,Nu);
363 foucault 27
364 francois 283 int inter_v = get_intervalle(nb_ptsctr_v,degre_v,v, knots_v);
365     double Nv [ 256 ];
366     deriver_fonction(inter_v,v,degre_v,dv,knots_v,Nv);
367    
368     OT_VECTEUR_4D temp [ 16 ];
369     for ( k=0; k<=du;k++)
370     {
371     for ( s=0; s<=degre_v;s++)
372 foucault 27 {
373 francois 283 temp[s]=0.0;
374     for ( r=0; r<=degre_u;r++)
375     temp[s]+=P((inter_u-degre_u+r),(inter_v-degre_v+s))*Nu(k,r);
376 foucault 27 }
377     dd = std::min(d-k,dv);
378 francois 283 for ( l=0; l<=dd;l++)
379 foucault 27 {
380 francois 283 skl(k,l) = 0.0;
381     for ( s=0; s<=degre_v;s++)
382     skl(k,l) += temp[s]*Nv(l,s);
383 foucault 27 }
384 francois 283 }
385 foucault 27
386     #undef skl
387     #undef f_deriver_u
388     #undef f_deriver_v
389     #undef P
390     }
391    
392     void ST_B_SPLINE_SURF::deriver_kieme(double *uv,int d,OT_VECTEUR_4D *skl)
393     {
394     #define skl(i,j) (*(skl+(i)*(d+1)+j))
395     #define ders(i,j) (*(ders+(i)*(d+1)+j))
396 francois 283 double u = uv[0];
397     double v = uv[1];
398     int i,j,k,l;
399     OT_VECTEUR_4D pv,pv2;
400 foucault 27
401 francois 283 OT_VECTEUR_4D ders [256];
402 foucault 27
403 francois 283 double Bin[256];
404     int dbin=d+1;
405     binomialCoef(Bin, dbin);
406 foucault 27 #define Bin(i,j) (*(Bin+(i)*(dbin)+j))
407    
408 francois 283 deriver_bs_kieme(nb_ptsctr_u,degre_u,knots_u,nb_ptsctr_v,degre_v,knots_v,u,v,d,ders);
409 foucault 27
410 francois 283 for (k=0;k<=d;++k)
411 foucault 27 {
412 francois 283 for (l=0;l<=d-k;++l)
413     {
414     pv = ders(k,l);
415     for (j=1;j<=l;j++)
416     pv -= Bin(l,j)*ders(0,j).w()*skl(k,l-j) ;
417     for (i=1;i<=k;i++)
418     {
419     pv -= Bin(k,i)*ders(i,0).w()*skl(k-i,l) ;
420     pv2 = 0.0 ;
421     for (j=1;j<=l;j++)
422     pv2 += Bin(l,j)*ders(i,j).w()*skl(k-i,l-j) ;
423     pv -= Bin(k,i)*pv2 ;
424     }
425     skl(k,l) = pv;
426     skl(k,l) /= ders[0].w();
427     }
428 foucault 27 }
429     #undef skl
430     #undef ders
431 francois 283 #undef Bin
432 foucault 27 }
433    
434    
435     void ST_B_SPLINE_SURF::deriver(double *uv,double *xyzdu, double *xyzdv)
436     {
437 francois 283 OT_VECTEUR_4D skl[4];
438     double uvl[3];
439     uvl[0]=uv[0];
440     uvl[1]=uv[1];
441     if (sens==-1) uvl[0]=umin+umax-uvl[0];
442     deriver_kieme(uvl,1,skl);
443 foucault 27 #define skl(i,j) skl[(i)*(2)+j]
444    
445 francois 283 for (int i=0; i<3; i++)
446     xyzdu[i] = sens * skl(1,0)[i];
447    
448     for (int i=0; i<3; i++)
449     xyzdv[i]=skl(0,1)[i];
450    
451     #undef skl
452 foucault 27 }
453    
454     void ST_B_SPLINE_SURF::deriver_seconde(double *uv,double* xyzduu,double* xyzduv,double* xyzdvv,double *xyz, double *xyzdu, double *xyzdv)
455 francois 283 {
456 foucault 27 #define skl(i,j) skl[(i)*(3)+j]
457 francois 283 OT_VECTEUR_4D skl[9];
458 foucault 27
459 francois 283 double uvl[3];
460     uvl[0]=uv[0];
461     uvl[1]=uv[1];
462     if (sens==-1) uvl[0]=umin+umax-uvl[0];
463     deriver_kieme(uvl,2,skl);
464 foucault 27
465 francois 283 for (int i=0; i<3; i++)
466     xyzduu[i]=skl(2,0)[i];
467 foucault 27
468 francois 283 for (int i=0; i<3; i++)
469     xyzduv[i] = sens * skl(1,1)[i];
470 foucault 27
471 francois 283 for (int i=0; i<3; i++)
472     xyzdvv[i]=skl(0,2)[i];
473 foucault 27
474 francois 283 if ((xyzdu!=NULL) && (xyzdv!=NULL ) )
475     {
476     for (int i=0; i<3; i++)
477     xyzdu[i] = sens * skl(1,0)[i];
478 foucault 27
479 francois 283 for (int i=0; i<3; i++)
480     xyzdv[i]=skl(0,1)[i];
481     }
482     if (xyz!=NULL)
483     {
484     for (int i=0; i<3; i++)
485     xyz[i]=skl(0,0)[i];
486     }
487 foucault 27
488 francois 283 unsigned kk, kkk;
489     double LIMITE = 1E99;
490     for (kk = 0; kk<3; kk++)
491 foucault 27 if ( (xyzduu[kk] > LIMITE ) || (xyzduv[kk] > LIMITE ) || (xyzdvv[kk] > LIMITE ) )
492 francois 283 break;
493     if (kk < 3)
494     {
495 foucault 27 deriver_seconde(uv,xyzduu,xyzduv,xyzdvv,xyz, xyzdu, xyzdv);
496     printf("pb");
497 francois 283 }
498 foucault 27
499     #undef skl
500     }
501    
502     void ST_B_SPLINE_SURF::inverser(double *uv,double *xyz,double precision)
503     {
504 francois 283 int code=0;
505     int num_point=(nb_ptsctr_u+nb_ptsctr_v)/2;
506     do
507     {
508 foucault 27 code=inverser2(uv,xyz,num_point,precision);
509     num_point=num_point*2;
510 francois 283 }
511     while (code==0&&num_point<1500);
512     }
513 foucault 27
514    
515     int ST_B_SPLINE_SURF::inverser2(double *uv,double *xyz,int num_test, double precision)
516     {
517 francois 283 double uvi[2];
518     double xyz1[3];
519     double xyzdu[3];
520     double xyzdv[3];
521     double xyzduu[3];
522     double xyzduv[3];
523     double xyzdvv[3];
524     OT_VECTEUR_3D Pt(xyz[0],xyz[1],xyz[2]);
525     double a[4];
526     double delta[2];
527     double b[9];
528     double delta_u;
529     double delta_v;
530     double ui;
531     double vi;
532 foucault 27
533 francois 283 double distance_ref=1e308;
534     int refu;
535     int refv;
536 foucault 27
537 francois 283 for (int ii=0;ii<num_test+1;ii++)
538     {
539 foucault 27 double u=umin+ii*1./num_test*(umax-umin);
540     uvi[0]=u;
541 francois 283 for (int jj=0;jj<num_test+1;jj++)
542     {
543     double v=vmin+jj*1./num_test*(vmax-vmin);
544     uvi[1]=v;
545     evaluer(uvi,xyz1);
546     OT_VECTEUR_3D S(xyz1[0],xyz1[1],xyz1[2]);
547     OT_VECTEUR_3D Distance = S-Pt;
548     double longueur=Distance.get_longueur2();
549     if (longueur<distance_ref)
550     {
551     distance_ref=longueur;
552     refu=ii;
553     refv=jj;
554     }
555 foucault 27 }
556 francois 283 }
557 foucault 27
558 francois 283 double uii=umin+refu*1./num_test*(umax-umin);
559     double vii=vmin+refv*1./num_test*(vmax-vmin);
560     int compteur = 0;
561     do
562     {
563 foucault 27 compteur++;
564     ui=uii;
565     vi=vii;
566     uvi[0]=ui;
567     uvi[1]=vi;
568     deriver_seconde(uvi,xyzduu,xyzduv,xyzdvv,xyz1,xyzdu,xyzdv);
569     OT_VECTEUR_3D S(xyz1[0],xyz1[1],xyz1[2]);
570     OT_VECTEUR_3D Su(xyzdu[0],xyzdu[1],xyzdu[2]);
571     OT_VECTEUR_3D Sv(xyzdv[0],xyzdv[1],xyzdv[2]);
572     OT_VECTEUR_3D Suu(xyzduu[0],xyzduu[1],xyzduu[2]);
573     OT_VECTEUR_3D Suv(xyzduv[0],xyzduv[1],xyzduv[2]);
574     OT_VECTEUR_3D Svv(xyzdvv[0],xyzdvv[1],xyzdvv[2]);
575     OT_VECTEUR_3D Distance = S-Pt;
576     a[0]=Su.get_longueur2()+Distance*Suu;
577     a[1]=Su*Sv+Distance*Suv;
578     a[2]=Su*Sv+Distance*Suv;
579     a[3]=Sv.get_longueur2()+Distance*Svv;
580 francois 283 b[0]=Distance*Su;
581     b[0]=-b[0];
582     b[1]=Distance*Sv;
583     b[1]=-b[1];
584 foucault 27 double denominateur_delta=(a[0]*a[3]-a[2]*a[1]);
585 francois 35 if (fabs(denominateur_delta) < ( (fabs(a[0])+fabs(a[1])+fabs(a[2])+fabs(a[3]))*1e-12 ) )
586 foucault 27 return 0;
587     delta_u=(b[0]*a[3]-b[1]*a[1])/denominateur_delta;
588     delta_v=(a[0]*b[1]-a[2]*b[0])/denominateur_delta;
589 francois 35 if (fabs(delta_u)>umax-umin) return 0;
590     if (fabs(delta_v)>vmax-vmin) return 0;
591 foucault 27 if (Su.get_longueur2()>1E-14)
592     uii=ui+delta_u;
593     else
594     uii=ui;
595     if (Sv.get_longueur2()>1E-14)
596     vii=vi+delta_v;
597     else
598     vii=vi;
599     if (periodique_u==0)
600 francois 283 {
601     if (uii<umin) uii=umin;
602     if (uii>umax) uii=umax;
603     }
604 foucault 27 if (periodique_v==0)
605 francois 283 {
606     if (vii<vmin) vii=vmin;
607     if (vii>vmax) vii=vmax;
608     }
609 foucault 27 if (periodique_u==1)
610 francois 283 {
611     if (uii<umin) uii=umax - (umin-uii);
612     if (uii>umax) uii=umin + (uii-umax);
613     }
614 foucault 27 if (periodique_v==1)
615 francois 283 {
616     if (vii<vmin) vii=vmax - (vmin-vii);
617     if (vii>vmax) vii=vmin + (vii-vmax);
618     }
619 foucault 27 delta_u=uii-ui;
620     delta_v=vii-vi;
621     if (compteur>500) return 0;
622 francois 283 }
623 foucault 27
624 francois 283 while ((fabs(delta_u)>precision)||(fabs(delta_v)>precision));
625     uv[0]=uii;
626     uv[1]=vii;
627     return 1;
628 foucault 27 }
629    
630     int ST_B_SPLINE_SURF::est_periodique_u(void)
631     {
632 francois 283 return periodique_u;
633 foucault 27 }
634     int ST_B_SPLINE_SURF::est_periodique_v(void)
635     {
636 francois 283 return periodique_v;
637 foucault 27 }
638     double ST_B_SPLINE_SURF::get_periode_u(void)
639     {
640 francois 283 return periode_u;
641 foucault 27 }
642     double ST_B_SPLINE_SURF::get_periode_v(void)
643     {
644 francois 283 return periode_v;
645 foucault 27 }
646     double ST_B_SPLINE_SURF::get_umin(void)
647     {
648 francois 283 return umin;
649 foucault 27 }
650     double ST_B_SPLINE_SURF::get_umax(void)
651     {
652 francois 283 return umax;
653 foucault 27 }
654     double ST_B_SPLINE_SURF::get_vmin(void)
655     {
656 francois 283 return vmin;
657 foucault 27 }
658     double ST_B_SPLINE_SURF::get_vmax(void)
659     {
660 francois 283 return vmax;
661 foucault 27 }
662    
663     int ST_B_SPLINE_SURF::get_type_geometrique(TPL_LISTE_ENTITE<double> &param)
664     {
665 francois 283 for (int i=0;i<nb_ptsctr_u-(degre_u+1);i++)
666     {
667     param.ajouter(knots_u[i]);
668     }
669     for (int i=0;i<nb_ptsctr_v-(degre_v+1);i++)
670     {
671     param.ajouter(knots_v[i]);
672     }
673     double xyz[3];
674     for (int i=0;i<nb_point;i++)
675     {
676     param.ajouter(ptsctr[i].x()/ptsctr[i].w());
677     param.ajouter(ptsctr[i].y()/ptsctr[i].w());
678     param.ajouter(ptsctr[i].z()/ptsctr[i].w());
679     }
680     for (int i=0;i<nb_point;i++)
681     {
682     param.ajouter(ptsctr[i].w());
683     }
684     param.ajouter(degre_u);
685     param.ajouter(degre_v);
686 francois 1149 return GEOMETRIE::CONST::Co_BSPLINES;
687 foucault 27 }
688    
689     void ST_B_SPLINE_SURF::initialiser(class ST_GESTIONNAIRE* gest)
690     {
691 francois 283 for (int i=0;i<nb_point;i++)
692     {
693 foucault 27 ST_POINT* stpoint=gest->lst_point.getid(indexptsctr[i]);
694     double xyz[3];
695     stpoint->evaluer(xyz);
696     double w=1.0;
697     ptsctr.push_back(OT_VECTEUR_4D(w*xyz[0],w*xyz[1],w*xyz[2],w));
698 francois 283 }
699 foucault 27
700 francois 283 periodique_u=1;
701     for (int i=0; i<nb_ptsctr_v; i++)
702     {
703 foucault 27 ST_POINT* point1_u=gest->lst_point.getid(indexptsctr[i]);
704     ST_POINT* point2_u=gest->lst_point.getid(indexptsctr[nb_ptsctr_v*(nb_ptsctr_u-1)+i]);
705     double xyz1[3];
706     double xyz2[3];
707     point1_u->evaluer(xyz1);
708     point2_u->evaluer(xyz2);
709     if ((!(OPERATEUR::egal (xyz1[0],xyz2[0],1E-10))) || (!(OPERATEUR::egal (xyz1[1],xyz2[1],1E-10))) || (!(OPERATEUR::egal (xyz1[2],xyz2[2],1E-10))))
710 francois 283 {
711     periodique_u=0;
712 foucault 27 }
713 francois 283 }
714     if (periodique_u==1)
715     {
716 foucault 27 int i=knots_u.size();
717     periode_u=(knots_u[i-1]-knots_u[0]);
718 francois 283 }
719     else periode_u=0;
720     periodique_v=1;
721     for (int j=0; j<nb_ptsctr_u; j++)
722     {
723 foucault 27 ST_POINT* point1_v=gest->lst_point.getid(indexptsctr[j*nb_ptsctr_v]);
724     ST_POINT* point2_v=gest->lst_point.getid(indexptsctr[j*nb_ptsctr_v+nb_ptsctr_v-1]);
725     double xyz3[3];
726     double xyz4[3];
727     point1_v->evaluer(xyz3);
728     point2_v->evaluer(xyz4);
729     if ((!(OPERATEUR::egal (xyz3[0],xyz4[0],1E-10))) || (!(OPERATEUR::egal (xyz3[1],xyz4[1],1E-10))) || (!(OPERATEUR::egal (xyz3[2],xyz4[2],1E-10))))
730 francois 283 {
731     periodique_v=0;
732 foucault 27 }
733 francois 283 }
734     if (periodique_v==1)
735     {
736 foucault 27 int j=knots_v.size();
737     periode_v=(knots_v[j-1]-knots_v[0]);
738 francois 283 }
739     else periode_v=0;
740 foucault 27 }
741    
742    
743    
744     void ST_B_SPLINE_SURF::est_util(ST_GESTIONNAIRE* gest)
745     {
746 francois 283 util=true;
747     for (int i=0;i<nb_point;i++)
748 foucault 27 gest->lst_point.getid(indexptsctr[i])->est_util(gest);
749     }
750    
751    
752    
753     void ST_B_SPLINE_SURF::get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
754     {
755    
756    
757 francois 283 param.ajouter(2);
758 foucault 27
759 francois 283 param.ajouter(degre_u+1);
760     param.ajouter(degre_v+1);
761 foucault 27
762    
763 francois 283 param.ajouter(nb_ptsctr_u);
764     param.ajouter(nb_ptsctr_v);
765 foucault 27
766    
767 francois 283 for (unsigned int i=0;i<knots_u.size();i++)
768     {
769     param.ajouter(knots_u[i]);
770     }
771 foucault 27
772    
773 francois 283 for (unsigned int j=0;j<knots_v.size();j++)
774     {
775     param.ajouter(knots_v[j]);
776     }
777 foucault 27
778 francois 283 for (int v=0;v<nb_ptsctr_v;v++)
779     {
780     for (int u=0;u<nb_ptsctr_u;u++)
781     {
782     int pt=u*nb_ptsctr_v+v;
783     double w=ptsctr[pt].w();
784     double inv_w=1/w;
785     double x=ptsctr[pt].x()*inv_w;
786     double y=ptsctr[pt].y()*inv_w;
787     double z=ptsctr[pt].z()*inv_w;
788     param.ajouter(x);
789     param.ajouter(y);
790     param.ajouter(z);
791     param.ajouter(w);
792     }
793 foucault 27
794 francois 283 }
795 souaissa 57
796 francois 283 /*
797     for(unsigned int pt=0;pt<ptsctr.size();pt++)
798     {
799 foucault 27
800 francois 283 double w=ptsctr[pt].w();
801     double inv_w=1/w;
802     double x=ptsctr[pt].x()*inv_w;
803     double y=ptsctr[pt].y()*inv_w;
804     double z=ptsctr[pt].z()*inv_w;
805     param.ajouter(x);
806     param.ajouter(y);
807     param.ajouter(z);
808     param.ajouter(w);
809 foucault 27
810 francois 283 // OT_VECTEUR_4D ctrpoint=ptsctr[pt];
811     // for(int i=0;i<4;i++)
812     // {
813     // param.ajouter(ctrpoint[i]);
814     // }
815 foucault 27
816    
817    
818 francois 283 } */
819 foucault 27
820 francois 283 indx_premier_ptctr=5+knots_u.size()+knots_v.size();
821 foucault 27
822    
823    
824    
825 souaissa 57
826 francois 283
827    
828    
829 foucault 27 }