ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/step/src/stbspline.~cpp
Revision: 5
Committed: Tue Jun 12 20:26:34 2007 UTC (17 years, 11 months ago)
Original Path: magic/lib/step/step/src/stbspline.~cpp
File size: 15084 byte(s)
Log Message:

File Contents

# User Rev Content
1 5 //------------------------------------------------------------
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     #include "ot_mathematique.h"
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     ptsctr.insert(ptsctr.end(),vec_point[3*i]);
63     ptsctr.insert(ptsctr.end(),vec_point[3*i+1]);
64     ptsctr.insert(ptsctr.end(),vec_point[3*i+2]);
65     poids.insert(poids.end(),vec_poids[i]);
66     }
67     double xyz1[3];
68     double xyz2[3];
69     xyz1[0]=vec_point[0];
70     xyz1[1]=vec_point[1];
71     xyz1[2]=vec_point[2];
72     xyz2[0]=vec_point[3*nb_point-3];
73     xyz2[1]=vec_point[3*nb_point-2];
74     xyz2[2]=vec_point[3*nb_point-1];
75     periodique=0;
76     if (OPERATEUR::egal (xyz1[0],xyz2[0],1E-6))
77     {
78     if (OPERATEUR::egal (xyz1[1],xyz2[1],1E-6))
79     {
80     if (OPERATEUR::egal (xyz1[2],xyz2[2],1E-6))
81     periodique=1;
82     }
83     }
84    
85     if (periodique==1)
86     {
87     int i=knots.size();
88     periode=(knots[i-1]-knots[0]);
89     }
90     else periode=0;
91     }
92    
93     ST_B_SPLINE::~ST_B_SPLINE()
94     {
95     }
96    
97     int ST_B_SPLINE::get_intervalle(int nb_point, int degre, double t, std::vector<double> &knots)
98     {
99     int inter;
100     if (OPERATEUR::egal(t,knots[nb_point],1E-12)==1) inter=nb_point-1;
101     else
102     {
103     int low=degre;
104     int high=nb_point+1;
105     int mid=((low+high)/2);
106     while ((t<knots[mid-1]) || (t>=knots[mid]))
107     {
108     if (t<knots[mid-1]) high=mid;
109     else low=mid;
110     mid=(low+high)/2;
111     }
112     inter=mid-1;
113     }
114     return inter;
115     }
116    
117     // Setup the binomial coefficients into th matrix Bin
118     // Bin(i,j) = (i j)
119     // The binomical coefficients are defined as follow
120     // (n) n!
121     // (k) = k!(n-k)! 0<=k<=n
122     // and the following relationship applies
123     // (n+1) (n) ( n )
124     // ( k ) = (k) + (k-1)
125     /*!
126     \brief Setup a matrix containing binomial coefficients
127    
128     Setup the binomial coefficients into th matrix Bin
129     \htmlonly
130     \[ Bin(i,j) = \left( \begin{array}{c}i \\ j\end{array} \right)\]
131     The binomical coefficients are defined as follow
132     \[ \left(\begin{array}{c} n \\ k \end{array} \right)= \frac{ n!}{k!(n-k)!} \mbox{for $0\leq k \leq n$} \]
133     and the following relationship applies
134     \[ \left(\begin{array}{c} n+1 \\ k \end{array} \right) =
135     \left(\begin{array}{c} n \\ k \end{array} \right) +
136     \left(\begin{array}{c} n \\ k-1 \end{array} \right) \]
137     \endhtmlonly
138    
139     \param Bin the binomial matrix
140     \author Philippe Lavoie
141     \date 24 January, 1997
142     */
143     void ST_B_SPLINE::binomialCoef(double * Bin, int degre){
144     #define Bin(i,j) (*(Bin+(i)*(3)+j))
145     int n,k;
146     // Setup the first line
147     Bin(0,0) = 1.0 ;
148     for(k=degre-1;k>0;--k)
149     Bin(0,k) = 0.0 ;
150     // Setup the other lines
151     for(n=0;n<degre-1;n++){
152     Bin(n+1,0) = 1.0 ;
153     for(k=1;k<degre;k++)
154     if(n+1<k)
155     Bin(n,k) = 0.0 ;
156     else
157     Bin(n+1,k) = Bin(n,k) + Bin(n,k-1) ;
158     }
159     #undef Bin
160     }
161    
162     void ST_B_SPLINE::get_valeur_fonction(int inter, double t, int degre, std::vector<double> &knots,double *grand_n)
163     {
164     //inter=get_intervalle(nb_point,degre,t,knots);
165     double saved;
166     grand_n[0]=1.0;
167     double gauche[16];
168     double droite[16];
169     for (int j=1;j<=degre;j++)
170     {
171     gauche[j-1]= t-knots[inter-j+1];
172     droite[j-1]=knots[inter+j]-t;
173     saved=0.0;
174     for (int r=0;r<j;r++)
175     {
176     double temp=grand_n[r]/(droite[r]+ gauche[j-r-1]);
177     grand_n[r]=saved+droite[r]* temp;
178     saved=gauche[j-r-1]*temp;
179     }
180     grand_n[j]=saved;
181     }
182     }
183    
184    
185     void ST_B_SPLINE::evaluer(double t,double *Cw)
186     {
187     int inter = get_intervalle(nb_point,degre,t, knots);
188     double grand_n[256];
189     get_valeur_fonction(inter,t,degre,knots,grand_n);
190    
191     // somme ( Ni,p(t) * Wi * Pi ) / somme ( Ni,p(t) * Wi )
192    
193     double R[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
194     double Sum_Ni_wi = 0;
195    
196     for (int j=0; j<=degre; j++)
197     {
198     Sum_Ni_wi += grand_n[j] * poids[inter-degre+j];
199     }
200    
201     for (int j=0; j<=degre; j++)
202     {
203     R[j] = grand_n[j] * poids[inter-degre+j] / Sum_Ni_wi;
204     }
205    
206     for (int i=0; i<3; i++)
207     {
208     Cw[i] = 0;
209     for (int j=0; j<=degre; j++)
210     {
211     Cw[i] += R[j] * ptsctr[3*(inter-degre+j)+i];
212     }
213    
214     }
215    
216     }
217    
218     void ST_B_SPLINE::deriver_pts(int nb_point,int degre,std::vector<double> &knots,std::vector<double> &ptsctr_x,int d,int r1,int r2,double *PK_x)
219     {
220     int r = r2-r1;
221     #define PK_x(i,j) (*(PK_x+(i)*(degre+1)+j))
222     for (int i=0; i<=r; i++)
223     {
224     PK_x(0,i)= ptsctr_x[r1+i];
225     }
226     for (int k=1; k<=d; k++)
227     {
228     int tmp = degre-k+1;
229     for (int i=0; i<=r-k; i++)
230     {
231     PK_x(k,i)= tmp*(PK_x(k-1,i+1)-PK_x(k-1,i))/(knots[r1+i+degre+1]-knots[r1+i+k]);
232     }
233     }
234     #undef PK_x // enlever (i,j)
235     }
236    
237     void ST_B_SPLINE::get_tout_fonction(int inter, double t, int degre, std::vector<double> &knots,double *grand_n)
238     {
239     //inter=get_intervalle(nb_point,degre,t,knots);
240     #define grand_n(i,j) (*(grand_n+(i)*(degre+1)+j))
241     double saved;
242     grand_n(0,0)=1.0;
243     double gauche[16];
244     double droite[16];
245     for (int j=1;j<=degre;j++)
246     {
247     gauche[j-1]= t-knots[inter-j+1];
248     droite[j-1]=knots[inter+j]-t;
249     saved=0.0;
250     for (int r=0;r<j;r++)
251     {
252     grand_n(j,r)=(droite[r]+ gauche[j-r-1]);
253     double temp = grand_n(r,j-1)/grand_n(j,r);
254    
255     grand_n(r,j)= saved+droite[r]*temp;
256     saved=gauche[j-r-1]*temp;
257     }
258     grand_n(j,j)=saved;
259     }
260     #undef grand_n // enlever (i,j)
261     }
262    
263     void ST_B_SPLINE::deriver_bs_kieme(int nb_point,int degre,std::vector<double> &knots,std::vector<double> &ptsctr_x,double t,int d,double *CK)
264     {
265     int du = std::min(d,degre);
266     for (int k=degre+1;k<=d;k++)
267     {
268     CK[k]=0.0;
269     }
270     int inter = get_intervalle(nb_point,degre,t, knots);
271     double grand_n[256];
272     get_tout_fonction(inter,t,degre,knots,grand_n);
273    
274     #define grand_n(i,j) (*(grand_n+(i)*(degre+1)+j))
275     double PK[256];
276     deriver_pts(nb_point,degre,knots,ptsctr_x,du,inter-degre,inter,PK);
277     #define PK(i,j) (*(PK+(i)*(degre+1)+j))
278     for(int k=0; k<=du;k++)
279     {
280     CK[k] = 0.0;
281     for(int j=0; j<=degre-k;j++)
282     {
283     CK[k] = CK[k]+ grand_n(j,degre-k)*PK(k,j); //PK(k,j) PK[k*degre+j]
284     }
285     }
286     #undef PK // enlever (i,j)
287     #undef grand_n // enlever (i,j)
288     }
289    
290    
291    
292     void ST_B_SPLINE::deriver_kieme(double t,int d,double *CK_x,double *CK_y,double *CK_z)
293     {
294     double Aders_x[16];
295     double Aders_y[16];
296     double Aders_z[16];
297     std::vector<double> ptsctr_wx;
298     std::vector<double> ptsctr_wy;
299     std::vector<double> ptsctr_wz;
300     double xyz[3];
301     for (int i=0;i<nb_point;i++)
302     {
303     xyz[0]=ptsctr[3*i];
304     xyz[1]=ptsctr[3*i+1];
305     xyz[2]=ptsctr[3*i+2];
306     ptsctr_wx.insert(ptsctr_wx.end(),poids[i]*xyz[0]);
307     ptsctr_wy.insert(ptsctr_wy.end(),poids[i]*xyz[1]);
308     ptsctr_wz.insert(ptsctr_wz.end(),poids[i]*xyz[2]);
309     }
310     deriver_bs_kieme(nb_point,degre,knots,ptsctr_wx,t,d,Aders_x);
311     deriver_bs_kieme(nb_point,degre,knots,ptsctr_wy,t,d,Aders_y);
312     deriver_bs_kieme(nb_point,degre,knots,ptsctr_wz,t,d,Aders_z);
313     double Wders[16];
314     deriver_bs_kieme(nb_point,degre,knots,poids,t,d,Wders);
315     double Bin[256];
316     binomialCoef(Bin, d);
317     #define Bin(i,j) (*(Bin+(i)*(3)+j))
318     for(int k=0;k<=d;k++)
319     {
320     double v_x = Aders_x[k];
321     double v_y = Aders_y[k];
322     double v_z = Aders_z[k];
323     for(int i=1;i<=k;i++)
324     {
325     v_x = v_x - Bin(k-1,i-1)*Wders[i]*CK_x[k-i];
326     v_y = v_y - Bin(k-1,i-1)*Wders[i]*CK_y[k-i];
327     v_z = v_z - Bin(k-1,i-1)*Wders[i]*CK_z[k-i];
328     }
329     CK_x[k] = v_x/Wders[0];
330     CK_y[k] = v_y/Wders[0];
331     CK_z[k] = v_z/Wders[0];
332     }
333     #undef Bin // enlever (i,j)
334     }
335    
336     void ST_B_SPLINE::deriver(double t,double *dxyz)
337     {
338     double CK_x[2];
339     double CK_y[2];
340     double CK_z[2];
341     deriver_kieme(t,1,CK_x,CK_y,CK_z);
342     dxyz[0]=CK_x[1];
343     dxyz[1]=CK_y[1];
344     dxyz[2]=CK_z[1];
345     }
346    
347     void ST_B_SPLINE::deriver_seconde(double t,double *ddxyz,double* dxyz ,double* xyz )
348     {
349     double CK_x[3];
350     double CK_y[3];
351     double CK_z[3];
352     deriver_kieme(t,2,CK_x,CK_y,CK_z);
353     ddxyz[0]=CK_x[2];
354     ddxyz[1]=CK_y[2];
355     ddxyz[2]=CK_z[2];
356     if (dxyz!=NULL)
357     {
358     dxyz[0]=CK_x[1];
359     dxyz[1]=CK_y[1];
360     dxyz[2]=CK_z[1];
361     }
362     if (xyz!=NULL)
363     {
364     xyz[0]=CK_x[0];
365     xyz[1]=CK_y[0];
366     xyz[2]=CK_z[0];
367     }
368     }
369    
370    
371    
372     void ST_B_SPLINE::inverser(double& t,double *xyz,double precision)
373     {
374     int code;
375     int num_point=nb_point;
376     do
377     {
378     code=inverser2(t,xyz,num_point,precision);
379     num_point=num_point*2;
380     }
381     while (code==0);
382     }
383    
384     int ST_B_SPLINE::inverser2(double& t,double *xyz,int num_test,double precision)
385     {
386     double xyz1[3];
387     double dxyz1[3];
388     double ddxyz1[3];
389     double ti;
390     double eps;
391     double tmin=get_tmin();
392     double tmax=get_tmax();
393     OT_VECTEUR_3D Pt(xyz[0],xyz[1],xyz[2]);
394     double distance_ref=1e308;
395     int ref;
396    
397     for (int i=0;i<num_test+1;i++)
398     {
399     double t=tmin+i*1./num_test*(tmax-tmin);
400     evaluer(t,xyz1);
401     OT_VECTEUR_3D Ct(xyz1[0],xyz1[1],xyz1[2]);
402     OT_VECTEUR_3D Distance = Ct-Pt;
403     double longueur=Distance.get_longueur2();
404     if (longueur<distance_ref)
405     {
406     distance_ref=longueur;
407     ref=i;
408     }
409     }
410    
411     double tii=tmin+ref*1./num_test*(tmax-tmin);
412     int compteur=0;
413     do
414     {
415     compteur++;
416     ti=tii;
417     deriver_seconde(ti,ddxyz1,dxyz1,xyz1);
418     OT_VECTEUR_3D Ct(xyz1[0],xyz1[1],xyz1[2]);
419     OT_VECTEUR_3D Ct_deriver(dxyz1[0],dxyz1[1],dxyz1[2]);
420     OT_VECTEUR_3D Ct_deriver_seconde(ddxyz1[0],ddxyz1[1],ddxyz1[2]);
421     OT_VECTEUR_3D Distance = Ct-Pt;
422     tii=ti-Ct_deriver*Distance/(Ct_deriver_seconde*Distance+Ct_deriver.get_longueur2());
423     if (periodique==1)
424     {
425     if (tii<get_tmin()) tii=get_tmax()-(get_tmin()-tii);
426     if (tii>get_tmax()) tii=get_tmin()+(tii-get_tmax());
427     }
428     else
429     {
430     if (tii<get_tmin()) tii=get_tmin();
431     if (tii>get_tmax()) tii=get_tmax();
432     }
433     eps=fabs(tii-ti);
434     if (compteur>500) return 0;
435     }
436     while (eps>precision);
437     t=ti;
438     if ((periodique==1) && (OPERATEUR::egal(t,get_tmax(),1e-12)) t=get_tmin();
439     return 1;
440     }
441    
442     double ST_B_SPLINE::get_tmin()
443     {
444     return knots[0];
445     }
446     double ST_B_SPLINE::get_tmax()
447     {
448     int i=knots.size();
449     return knots[i-1];
450     }
451    
452     double equation_longueur(ST_B_SPLINE& bsp,double t)
453     {
454     double dxyz[3];
455     bsp.deriver(t,dxyz);
456     return sqrt(dxyz[0]*dxyz[0]+dxyz[1]*dxyz[1]+dxyz[2]*dxyz[2]);
457     }
458    
459    
460     double ST_B_SPLINE::get_longueur(double t1,double t2,double precis)
461     {
462     TPL_FONCTION1<double,ST_B_SPLINE,double> longueur_bsp(*this,equation_longueur);
463     return longueur_bsp.integrer_gauss_2(t1,t2);
464     }
465    
466    
467     int ST_B_SPLINE::est_periodique(void)
468     {
469     return periodique;
470     }
471     double ST_B_SPLINE::get_periode(void)
472     {
473     return periode;
474     }
475    
476     int ST_B_SPLINE::get_type_geometrique(TPL_LISTE_ENTITE<double> &param)
477     {
478     for(int i=0;i<nb_point-(degre+1);i++)
479     {
480     param.ajouter(knots[i]);
481     }
482     double xyz[3];
483     for(int i=0;i<nb_point;i++)
484     {
485     xyz[0]=ptsctr[3*i];
486     xyz[1]=ptsctr[3*i+1];
487     xyz[2]=ptsctr[3*i+2];
488     param.ajouter(xyz[0]);
489     param.ajouter(xyz[1]);
490     param.ajouter(xyz[2]);
491     }
492     for(int i=0;i<nb_point;i++)
493     {
494     param.ajouter(poids[i]);
495     }
496     param.ajouter(degre);
497     return MGCo_BSPLINE;
498     }
499    
500     void ST_B_SPLINE::initialiser(ST_GESTIONNAIRE *gest)
501     {
502     int i=indexptsctr.size();
503     ST_POINT* point1=gest->lst_point.getid(indexptsctr[0]);
504     ST_POINT* point2=gest->lst_point.getid(indexptsctr[i-1]);
505     double xyz1[3];
506     double xyz2[3];
507     point1->evaluer(xyz1);
508     point2->evaluer(xyz2);
509     periodique=0;
510     if (OPERATEUR::egal (xyz1[0],xyz2[0],1E-6))
511     {
512     if (OPERATEUR::egal (xyz1[1],xyz2[1],1E-6))
513     {
514     if (OPERATEUR::egal (xyz1[2],xyz2[2],1E-6))
515     periodique=1;
516     }
517     }
518    
519     if (periodique==1)
520     {
521     int i=knots.size();
522     periode=(knots[i-1]-knots[0]);
523     }
524     else periode=0;
525    
526     int nbptsctr=indexptsctr.size();
527     for (int i=0;i<nbptsctr;i++)
528     {
529     ST_POINT* stpoint=gest->lst_point.getid(indexptsctr[i]);
530     double xyz[3];
531     stpoint->evaluer(xyz);
532     ptsctr.insert(ptsctr.end(),xyz[0]);
533     ptsctr.insert(ptsctr.end(),xyz[1]);
534     ptsctr.insert(ptsctr.end(),xyz[2]);
535     }
536     }
537    
538    
539    
540    
541    
542     void ST_B_SPLINE::est_util(ST_GESTIONNAIRE* gest)
543     {
544     util=true;
545     for (int i=0;i<nb_point;i++)
546     gest->lst_point.getid(indexptsctr[i])->est_util(gest);
547     }
548