ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/step/src/stbspline.cpp
Revision: 27
Committed: Thu Jul 5 15:26:40 2007 UTC (17 years, 10 months ago) by foucault
Original Path: magic/lib/step/step/src/stbspline.cpp
File size: 15814 byte(s)
Log Message:

File Contents

# User Rev Content
1 foucault 27 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4     // Jean Christophe Cuillière et Vincent FRANCOIS
5     // Département de Génie Mécanique - UQTR
6     //------------------------------------------------------------
7     // Le projet MAGIC est un projet de recherche du département
8     // de génie mécanique de l'Université du Québec à
9     // Trois Rivières
10     // Les librairies ne peuvent être utilisées sans l'accord
11     // des auteurs (contact : francois@uqtr.ca)
12     //------------------------------------------------------------
13     //------------------------------------------------------------
14     //
15     // stbspline.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20     // Version du 02/03/2006 à 11H24
21     //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24    
25     #include "gestionversion.h"
26    
27     #include "stbspline.h"
28     #include <vector>
29     #include "st_gestionnaire.h"
30     #include "tpl_fonction.h"
31     #include "constantegeo.h"
32    
33     #include <math.h>
34    
35    
36    
37     ST_B_SPLINE::ST_B_SPLINE(long LigneCourante,std::string idori,int bs_degre,std::vector<int> bs_indexptsctr,std::vector<int> bs_knots_multiplicities,std::vector<double> bs_knots):ST_COURBE(LigneCourante,idori),degre(bs_degre)
38     {
39     int nb=bs_knots_multiplicities.size();
40     for (int i=0;i<nb;i++)
41     {
42     for (int j=0;j<bs_knots_multiplicities[i];j++)
43     knots.insert(knots.end(),bs_knots[i]);
44     }
45     nb_point=bs_indexptsctr.size();
46     for (int i=0;i<nb_point;i++)
47     {
48     indexptsctr.insert(indexptsctr.end(),bs_indexptsctr[i]);
49     //poids.insert(poids.end(),1.);
50     }
51    
52     }
53    
54     ST_B_SPLINE::ST_B_SPLINE(int bs_degre,std::vector<double> &vec_knots,std::vector<double> &vec_point,std::vector<double> &vec_poids):ST_COURBE(),degre(bs_degre)
55     {
56     int nb=vec_knots.size();
57     for (int i=0;i<nb;i++)
58     knots.insert(knots.end(),vec_knots[i]);
59     nb_point=vec_point.size()/3;
60     for (int i=0;i<nb_point;i++)
61     {
62     double w=vec_poids[i];
63     double x=vec_point[3*i];
64     double y=vec_point[3*i+1];
65     double z=vec_point[3*i+2];
66     OT_VECTEUR_4D pt(w*x,w*y,w*z,w);
67     ptsctr.push_back(pt);
68     }
69     double xyz1[3];
70     double xyz2[3];
71     xyz1[0]=vec_point[0];
72     xyz1[1]=vec_point[1];
73     xyz1[2]=vec_point[2];
74     xyz2[0]=vec_point[3*nb_point-3];
75     xyz2[1]=vec_point[3*nb_point-2];
76     xyz2[2]=vec_point[3*nb_point-1];
77     periodique=0;
78     if (OPERATEUR::egal (xyz1[0],xyz2[0],1E-6))
79     {
80     if (OPERATEUR::egal (xyz1[1],xyz2[1],1E-6))
81     {
82     if (OPERATEUR::egal (xyz1[2],xyz2[2],1E-6))
83     periodique=1;
84     }
85     }
86    
87     if (periodique==1)
88     {
89     int i=knots.size();
90     periode=(knots[i-1]-knots[0]);
91     }
92     else periode=0;
93     }
94    
95     ST_B_SPLINE::~ST_B_SPLINE()
96     {
97     }
98    
99     int ST_B_SPLINE::get_intervalle(int nb_point, int degre, double t, std::vector<double> &knots)
100     {
101     int inter;
102     if (OPERATEUR::egal(t,knots[nb_point],1E-10)==1) inter=nb_point-1;
103     else if (OPERATEUR::egal(t,knots[degre],1E-10)==1) inter=degre;
104     else
105     {
106     int low=degre;
107     int high=nb_point+1;
108     int mid=((low+high)/2);
109     while ((t<knots[mid-1]) || (t>=knots[mid]))
110     {
111     if (t<knots[mid-1]) high=mid;
112     else low=mid;
113     mid=(low+high)/2;
114     }
115     inter=mid-1;
116     }
117     return inter;
118     }
119    
120     void ST_B_SPLINE::binomialCoef(double * Bin, int d){
121     int n,k;
122     // Setup the first line
123     Bin[0] = 1.0 ;
124     for(k=d-1;k>0;--k)
125     Bin[k] = 0.0 ;
126     // Setup the other lines
127     for(n=0;n<d-1;n++){
128     Bin[(n+1)*d+0] = 1.0 ;
129     for(k=1;k<d;k++)
130     if(n+1<k)
131     Bin[(n+1)*d+k] = 0.0 ;
132     else
133     Bin[(n+1)*d+k] = Bin[n*d+k] + Bin[n*d+k-1] ;
134     }
135     }
136    
137     void ST_B_SPLINE::get_valeur_fonction(int inter, double t, int degre, std::vector<double> &knots,double *grand_n)
138     {
139     double saved;
140     grand_n[0]=1.0;
141     double gauche[16];
142     double droite[16];
143     for (int j=1;j<=degre;j++)
144     {
145     gauche[j-1]= t-knots[inter-j+1];
146     droite[j-1]=knots[inter+j]-t;
147     saved=0.0;
148     for (int r=0;r<j;r++)
149     {
150     double temp=grand_n[r]/(droite[r]+ gauche[j-r-1]);
151     grand_n[r]=saved+droite[r]* temp;
152     saved=gauche[j-r-1]*temp;
153     }
154     grand_n[j]=saved;
155     }
156     }
157    
158    
159     void ST_B_SPLINE::evaluer(double t,double *xyz)
160     {
161     if (est_periodique())
162     {
163     double tmin=get_tmin();
164     double tmax=get_tmax();
165     while (t>tmax) t -= periode;
166     while (t<tmin) t += periode;
167     }
168     int inter = get_intervalle(nb_point,degre,t, knots);
169     double grand_n[256];
170     get_valeur_fonction(inter,t,degre,knots,grand_n);
171    
172     OT_VECTEUR_4D sp(0,0,0,0) ;
173     for (int j=0; j<=degre; j++)
174     {
175     sp += grand_n[j] * ptsctr[inter-degre+j];
176     }
177    
178     // transform homogeneous coordinates to 3D coordinates
179     for (int i=0; i<3; i++)
180     xyz[i] = sp[i]/sp.w();
181     }
182     /*
183     void ST_B_SPLINE::get_tout_fonction(int inter, double t, int degre, std::vector<double> &knots,double *grand_n)
184     {
185     #define grand_n(i,j) (*(grand_n+(i)*(degre+1)+j))
186     double saved;
187     grand_n(0,0)=1.0;
188     double gauche[16];
189     double droite[16];
190     for (int j=1;j<=degre;j++)
191     {
192     gauche[j-1]= t-knots[inter-j+1];
193     droite[j-1]=knots[inter+j]-t;
194     saved=0.0;
195     for (int r=0;r<j;r++)
196     {
197     grand_n(j,r)=(droite[r]+ gauche[j-r-1]);
198     double temp = grand_n(r,j-1)/grand_n(j,r);
199    
200     grand_n(r,j)= saved+droite[r]*temp;
201     saved=gauche[j-r-1]*temp;
202     }
203     grand_n(j,j)=saved;
204     }
205     #undef grand_n // enlever (i,j)
206     } */
207    
208     void ST_B_SPLINE::deriver_fonction(int inter,double t,int degre,int dd,std::vector<double> &knots,double *f_deriver)
209     {
210     #define f_deriver(i,j) (*(f_deriver+(i)*(degre+1)+j))
211     #define grand_n(i,j) (*(grand_n+(i)*(degre+1)+j))
212     #define a(i,j) (*(a+(i)*(dd+1)+j))
213     double *grand_n=new double[(degre+1)*(degre+1)];
214     double saved;
215     grand_n(0,0)=1.0;
216     double *gauche=new double[degre+1];
217     double *droite=new double[degre+1];
218     for (int j=1;j<=degre;j++)
219     {
220     gauche[j]= t-knots[inter-j+1];
221     droite[j]=knots[inter+j]-t;
222     saved=0.0;
223     for (int r=0;r<j;r++)
224     {
225     grand_n(j,r)=(droite[r+1]+ gauche[j-r]);
226     double temp = grand_n(r,j-1)/grand_n(j,r);
227    
228     grand_n(r,j)= saved+droite[r+1]*temp;
229     saved=gauche[j-r]*temp;
230     }
231     grand_n(j,j)=saved;
232     }
233     for (int j=0;j<=degre;j++)
234     {
235     f_deriver(0,j)= grand_n(j,degre);
236     }
237     double *a=new double[(degre+1)*(degre+1)];
238     for (int r=0;r<=degre;r++)
239     {
240     int s1=0; int s2=1;
241     a(0,0)=1.0;
242     for (int k=1;k<=dd; k++)
243     {
244     double d=0.0;
245     int rk=r-k; int pk=degre-k;
246     if(r>=k)
247     {
248     a(s2,0)=a(s1,0)/grand_n(pk+1,rk);
249     d= a(s2,0)* grand_n(rk,pk);
250     }
251     int j1; int j2;
252     if(rk>=-1) j1=1;
253     else j1= -rk;
254     if(r-1<=pk) j2=k-1;
255     else j2=degre-r;
256     for (int j=j1;j<=j2;j++)
257     {
258     a(s2,j) = (a(s1,j)-a(s1,j-1))/grand_n(pk+1,rk+j);
259     d+=a(s2,j)*grand_n(rk+j,pk);
260     }
261     if(r<=pk)
262     {
263     a(s2,k) = -a(s1,k-1)/grand_n(pk+1,r);
264     d+=a(s2,k)*grand_n(r,pk);
265     }
266     f_deriver(k,r)=d;
267     int j=s1; s1=s2; s2=j;
268     }
269     }
270     int r=degre;
271     for (int k=1;k<=dd;k++)
272     {
273     for (int j=0;j<=degre;j++)
274     {
275     f_deriver(k,j)*=r;
276     }
277     r*=(degre-k);
278     }
279     delete [] a;
280     delete [] grand_n;
281     delete [] gauche;
282     delete [] droite;
283     #undef f_deriver
284     #undef grand_n
285     #undef a
286     }
287    
288     void ST_B_SPLINE::deriver_bs_kieme(int nb_point,int degre,std::vector<double> &knots,std::vector<OT_VECTEUR_4D> &ptsctr_x,double t,int d, OT_VECTEUR_4D * CK)
289     {
290     int du = std::min(d,degre);
291     int inter = get_intervalle(nb_point,degre,t, knots);
292     double derF[256];
293     deriver_fonction(inter, t,degre,du,knots,derF);
294    
295     #define derF(i,j) (*(derF+(i)*(degre+1)+j))
296     for(int k=du;k>=0;--k)
297     {
298     CK[k] = OT_VECTEUR_4D(0,0,0,0);
299     for(int j=degre;j>=0;--j)
300     {
301     CK[k] += derF(k,j)*ptsctr_x[inter-degre+j];
302     }
303     }
304     #undef derF
305     }
306    
307    
308    
309     void ST_B_SPLINE::deriver_kieme(double t,int d, double *CK_x,double *CK_y,double *CK_z)
310     {
311     int i,k;
312     OT_VECTEUR_4D dersW[16], ders[16];
313     OT_VECTEUR_4D v;
314     deriver_bs_kieme(nb_point, degre, knots, ptsctr, t, d, dersW);
315    
316     double Bin[256];
317     int dbin=d+1;
318     binomialCoef(Bin, dbin);
319     #define Bin(i,j) (*(Bin+(i)*(dbin)+j))
320    
321     for(k=0;k<=d;k++)
322     {
323     ders[k] = dersW[k];
324     for(i=k ;i>0 ;--i)
325     ders[k] -= (Bin(k,i)*dersW[i].w())*ders[k-i];
326     ders[k]/=dersW[0].w();
327     }
328    
329     for (int i=0; i<=d; i++)
330     {
331     CK_x[i]=ders[i][0];
332     CK_y[i]=ders[i][1];
333     CK_z[i]=ders[i][2];
334     }
335     }
336    
337     void ST_B_SPLINE::deriver(double t,double *dxyz)
338     {
339     if (est_periodique())
340     {
341     double tmin=get_tmin();
342     double tmax=get_tmax();
343     while (t>tmax) t -= periode;
344     while (t<tmin) t += periode;
345     }
346     double CK_x[2];
347     double CK_y[2];
348     double CK_z[2];
349     deriver_kieme(t,1,CK_x,CK_y,CK_z);
350     dxyz[0]=CK_x[1];
351     dxyz[1]=CK_y[1];
352     dxyz[2]=CK_z[1];
353     }
354    
355     void ST_B_SPLINE::deriver_seconde(double t,double *ddxyz,double* dxyz ,double* xyz )
356     {
357     if (est_periodique())
358     {
359     double tmin=get_tmin();
360     double tmax=get_tmax();
361     while (t>tmax) t -= periode;
362     while (t<tmin) t += periode;
363     }
364     double CK_x[3];
365     double CK_y[3];
366     double CK_z[3];
367     deriver_kieme(t,2,CK_x,CK_y,CK_z);
368     ddxyz[0]=CK_x[2];
369     ddxyz[1]=CK_y[2];
370     ddxyz[2]=CK_z[2];
371     if (dxyz!=NULL)
372     {
373     dxyz[0]=CK_x[1];
374     dxyz[1]=CK_y[1];
375     dxyz[2]=CK_z[1];
376     }
377     if (xyz!=NULL)
378     {
379     xyz[0]=CK_x[0];
380     xyz[1]=CK_y[0];
381     xyz[2]=CK_z[0];
382     }
383     }
384    
385    
386    
387     void ST_B_SPLINE::inverser(double& t,double *xyz,double precision)
388     {
389     int code;
390     int num_point=nb_point;
391     do
392     {
393     code=inverser2(t,xyz,num_point,precision);
394     num_point=num_point*2;
395     }
396     while (code==0 && num_point < 100000);
397     }
398    
399     int ST_B_SPLINE::inverser2(double& t,double *xyz,int num_test,double precision)
400     {
401     double xyz1[3];
402     double dxyz1[3];
403     double ddxyz1[3];
404     double ti;
405     double eps;
406     double tmin=get_tmin();
407     double tmax=get_tmax();
408     OT_VECTEUR_3D Pt(xyz[0],xyz[1],xyz[2]);
409     double distance_ref=1e308;
410     int ref;
411    
412     for (int i=0;i<num_test+1;i++)
413     {
414     double t=tmin+i*1./num_test*(tmax-tmin);
415     evaluer(t,xyz1);
416     OT_VECTEUR_3D Ct(xyz1[0],xyz1[1],xyz1[2]);
417     OT_VECTEUR_3D Distance = Ct-Pt;
418     double longueur=Distance.get_longueur2();
419     if (longueur<distance_ref)
420     {
421     distance_ref=longueur;
422     ref=i;
423     }
424     }
425    
426     double tii=tmin+ref*1./num_test*(tmax-tmin);
427     int compteur=0;
428     do
429     {
430     compteur++;
431     ti=tii;
432     deriver_seconde(ti,ddxyz1,dxyz1,xyz1);
433     OT_VECTEUR_3D Ct(xyz1[0],xyz1[1],xyz1[2]);
434     OT_VECTEUR_3D Ct_deriver(dxyz1[0],dxyz1[1],dxyz1[2]);
435     OT_VECTEUR_3D Ct_deriver_seconde(ddxyz1[0],ddxyz1[1],ddxyz1[2]);
436     OT_VECTEUR_3D Distance = Ct-Pt;
437     tii=ti-Ct_deriver*Distance/(Ct_deriver_seconde*Distance+Ct_deriver.get_longueur2());
438     if (periodique==1)
439     {
440     if (tii<get_tmin()) tii=get_tmax()-(get_tmin()-tii);
441     if (tii>get_tmax()) tii=get_tmin()+(tii-get_tmax());
442     }
443     else
444     {
445     if (tii<get_tmin()) tii=get_tmin();
446     if (tii>get_tmax()) tii=get_tmax();
447     }
448     eps=fabs(tii-ti);
449     if (compteur>500) return 0;
450     }
451     while (eps>precision);
452     t=ti;
453     return 1;
454     }
455    
456     double ST_B_SPLINE::get_tmin()
457     {
458     return knots[0];
459     }
460     double ST_B_SPLINE::get_tmax()
461     {
462     int i=knots.size();
463     return knots[i-1];
464     }
465    
466     double equation_longueur(ST_B_SPLINE& bsp,double t)
467     {
468     double dxyz[3];
469     bsp.deriver(t,dxyz);
470     return sqrt(dxyz[0]*dxyz[0]+dxyz[1]*dxyz[1]+dxyz[2]*dxyz[2]);
471     }
472    
473    
474     double ST_B_SPLINE::get_longueur(double t1,double t2,double precis)
475     {
476     TPL_FONCTION1<double,ST_B_SPLINE,double> longueur_bsp(*this,equation_longueur);
477     return longueur_bsp.integrer_gauss_2(t1,t2);
478     }
479    
480    
481     int ST_B_SPLINE::est_periodique(void)
482     {
483     return periodique;
484     }
485     double ST_B_SPLINE::get_periode(void)
486     {
487     return periode;
488     }
489    
490     int ST_B_SPLINE::get_type_geometrique(TPL_LISTE_ENTITE<double> &param)
491     {
492     for(int i=0;i<nb_point-(degre+1);i++)
493     {
494     param.ajouter(knots[i]);
495     }
496     double xyz[3];
497     for(int i=0;i<nb_point;i++)
498     {
499     xyz[0]=ptsctr[i].x()/ptsctr[i].w();
500     xyz[1]=ptsctr[i].y()/ptsctr[i].w();
501     xyz[2]=ptsctr[i].z()/ptsctr[i].w();
502     param.ajouter(xyz[0]);
503     param.ajouter(xyz[1]);
504     param.ajouter(xyz[2]);
505     }
506     for(int i=0;i<nb_point;i++)
507     {
508     param.ajouter(ptsctr[i].w());
509     }
510     param.ajouter(degre);
511     return MGCo_BSPLINE;
512     }
513    
514     void ST_B_SPLINE::initialiser(ST_GESTIONNAIRE *gest)
515     {
516     int i=indexptsctr.size();
517     ST_POINT* point1=gest->lst_point.getid(indexptsctr[0]);
518     ST_POINT* point2=gest->lst_point.getid(indexptsctr[i-1]);
519     double xyz1[3];
520     double xyz2[3];
521     point1->evaluer(xyz1);
522     point2->evaluer(xyz2);
523     periodique=0;
524     if (OPERATEUR::egal (xyz1[0],xyz2[0],1E-6))
525     {
526     if (OPERATEUR::egal (xyz1[1],xyz2[1],1E-6))
527     {
528     if (OPERATEUR::egal (xyz1[2],xyz2[2],1E-6))
529     periodique=1;
530     }
531     }
532    
533     if (periodique==1)
534     {
535     int i=knots.size();
536     periode=(knots[i-1]-knots[0]);
537     }
538     else periode=0;
539    
540     int nbptsctr=indexptsctr.size();
541     for (int i=0;i<nbptsctr;i++)
542     {
543     ST_POINT* stpoint=gest->lst_point.getid(indexptsctr[i]);
544     double xyz[3];
545     stpoint->evaluer(xyz);
546     OT_VECTEUR_4D pt(xyz[0],xyz[1],xyz[2],1);
547     ptsctr.insert(ptsctr.end(),pt);
548     }
549     }
550    
551    
552    
553    
554    
555     void ST_B_SPLINE::est_util(ST_GESTIONNAIRE* gest)
556     {
557     util=true;
558     for (int i=0;i<nb_point;i++)
559     gest->lst_point.getid(indexptsctr[i])->est_util(gest);
560     }
561    
562    
563    
564    
565     void ST_B_SPLINE::get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
566     {
567    
568     // the order stock of the parameters is:
569     // 1/ the variante type =1 (dimension one)
570     // 2/ the order of the spline in two dirction the second direction is zero
571     // 3/ the number of the control points
572     // 4/ the knots vector
573     // 5/ the control points cordinates in homogeneos forms
574    
575     // The first parameter indicate the code access
576     param.ajouter(1);
577     // The follewing two parameters of the list indicate the orders of the net points
578    
579     param.ajouter(degre+1);
580     param.ajouter(0);
581    
582     // The follewing two parameters indicate the number of rows and colons of the control points
583     // respectively to the two parameters directions
584    
585     param.ajouter(nb_point);
586     param.ajouter(0);
587    
588     // this present the knot vector in the u-direction
589    
590     for(unsigned int i=0;i<knots.size();i++)
591     {
592     param.ajouter(knots[i]);
593     }
594    
595     // in the following we construct the control point vector
596    
597     for(unsigned int pt=0;pt<ptsctr.size();pt++)
598     {
599     double w=ptsctr[pt].w();
600     double inv_w=1/w;
601     double x=ptsctr[pt].x()*inv_w;
602     double y=ptsctr[pt].y()*inv_w;
603     double z=ptsctr[pt].z()*inv_w;
604     param.ajouter(x);
605     param.ajouter(y);
606     param.ajouter(z);
607     param.ajouter(w);
608     }
609     indx_premier_ptctr=5+knots.size();
610     }