ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/step/src/stbspline.~cpp
Revision: 283
Committed: Tue Sep 13 21:11:20 2011 UTC (13 years, 8 months ago) by francois
File size: 15750 byte(s)
Log Message:
structure de l'écriture

File Contents

# User Rev Content
1 francois 283 //------------------------------------------------------------
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     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 d) {
144     int n,k;
145     // Setup the first line
146     Bin[0] = 1.0 ;
147     for (k=d-1;k>0;--k)
148     Bin[k] = 0.0 ;
149     // Setup the other lines
150     for (n=0;n<d-1;n++) {
151     Bin[(n+1)*d+0] = 1.0 ;
152     for (k=1;k<d;k++)
153     if (n+1<k)
154     Bin[(n+1)*d+k] = 0.0 ;
155     else
156     Bin[(n+1)*d+k] = Bin[n*d+k] + Bin[n*d+k-1] ;
157     }
158     }
159    
160     void ST_B_SPLINE::get_valeur_fonction(int inter, double t, int degre, std::vector<double> &knots,double *grand_n)
161     {
162     //inter=get_intervalle(nb_point,degre,t,knots);
163     double saved;
164     grand_n[0]=1.0;
165     double gauche[16];
166     double droite[16];
167     for (int j=1;j<=degre;j++)
168     {
169     gauche[j-1]= t-knots[inter-j+1];
170     droite[j-1]=knots[inter+j]-t;
171     saved=0.0;
172     for (int r=0;r<j;r++)
173     {
174     double temp=grand_n[r]/(droite[r]+ gauche[j-r-1]);
175     grand_n[r]=saved+droite[r]* temp;
176     saved=gauche[j-r-1]*temp;
177     }
178     grand_n[j]=saved;
179     }
180     }
181    
182    
183     void ST_B_SPLINE::evaluer(double t,double *Cw)
184     {
185     int inter = get_intervalle(nb_point,degre,t, knots);
186     double grand_n[256];
187     get_valeur_fonction(inter,t,degre,knots,grand_n);
188    
189     // somme ( Ni,p(t) * Wi * Pi ) / somme ( Ni,p(t) * Wi )
190    
191     double R[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
192     double Sum_Ni_wi = 0;
193    
194     for (int j=0; j<=degre; j++)
195     {
196     Sum_Ni_wi += grand_n[j] * poids[inter-degre+j];
197     }
198    
199     for (int j=0; j<=degre; j++)
200     {
201     R[j] = grand_n[j] * poids[inter-degre+j] / Sum_Ni_wi;
202     }
203    
204     for (int i=0; i<3; i++)
205     {
206     Cw[i] = 0;
207     for (int j=0; j<=degre; j++)
208     {
209     Cw[i] += R[j] * ptsctr[3*(inter-degre+j)+i];
210     }
211    
212     }
213    
214     }
215    
216     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)
217     {
218     int r = r2-r1;
219     #define PK_x(i,j) (*(PK_x+(i)*(degre+1)+j))
220     for (int i=0; i<=r; i++)
221     {
222     PK_x(0,i)= ptsctr_x[r1+i];
223     }
224     for (int k=1; k<=d; k++)
225     {
226     int tmp = degre-k+1;
227     for (int i=0; i<=r-k; i++)
228     {
229     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]);
230     }
231     }
232     #undef PK_x // enlever (i,j)
233     }
234    
235     void ST_B_SPLINE::get_tout_fonction(int inter, double t, int degre, std::vector<double> &knots,double *grand_n)
236     {
237     //inter=get_intervalle(nb_point,degre,t,knots);
238     #define grand_n(i,j) (*(grand_n+(i)*(degre+1)+j))
239     double saved;
240     grand_n(0,0)=1.0;
241     double gauche[16];
242     double droite[16];
243     for (int j=1;j<=degre;j++)
244     {
245     gauche[j-1]= t-knots[inter-j+1];
246     droite[j-1]=knots[inter+j]-t;
247     saved=0.0;
248     for (int r=0;r<j;r++)
249     {
250     grand_n(j,r)=(droite[r]+ gauche[j-r-1]);
251     double temp = grand_n(r,j-1)/grand_n(j,r);
252    
253     grand_n(r,j)= saved+droite[r]*temp;
254     saved=gauche[j-r-1]*temp;
255     }
256     grand_n(j,j)=saved;
257     }
258     #undef grand_n // enlever (i,j)
259     }
260    
261     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)
262     {
263     int du = std::min(d,degre);
264     for (int k=degre+1;k<=d;k++)
265     {
266     CK[k]=0.0;
267     }
268     int inter = get_intervalle(nb_point,degre,t, knots);
269     double grand_n[256];
270     get_tout_fonction(inter,t,degre,knots,grand_n);
271    
272     #define grand_n(i,j) (*(grand_n+(i)*(degre+1)+j))
273     double PK[256];
274     deriver_pts(nb_point,degre,knots,ptsctr_x,du,inter-degre,inter,PK);
275     #define PK(i,j) (*(PK+(i)*(degre+1)+j))
276     for (int k=0; k<=du;k++)
277     {
278     CK[k] = 0.0;
279     for (int j=0; j<=degre-k;j++)
280     {
281     CK[k] = CK[k]+ grand_n(j,degre-k)*PK(k,j); //PK(k,j) PK[k*degre+j]
282     }
283     }
284     #undef PK // enlever (i,j)
285     #undef grand_n // enlever (i,j)
286     }
287    
288    
289    
290     void ST_B_SPLINE::deriver_kieme(double t,int d,double *CK_x,double *CK_y,double *CK_z)
291     {
292     double Aders_x[16];
293     double Aders_y[16];
294     double Aders_z[16];
295     std::vector<double> ptsctr_wx;
296     std::vector<double> ptsctr_wy;
297     std::vector<double> ptsctr_wz;
298     double xyz[3];
299     for (int i=0;i<nb_point;i++)
300     {
301     xyz[0]=ptsctr[3*i];
302     xyz[1]=ptsctr[3*i+1];
303     xyz[2]=ptsctr[3*i+2];
304     ptsctr_wx.insert(ptsctr_wx.end(),poids[i]*xyz[0]);
305     ptsctr_wy.insert(ptsctr_wy.end(),poids[i]*xyz[1]);
306     ptsctr_wz.insert(ptsctr_wz.end(),poids[i]*xyz[2]);
307     }
308     deriver_bs_kieme(nb_point,degre,knots,ptsctr_wx,t,d,Aders_x);
309     deriver_bs_kieme(nb_point,degre,knots,ptsctr_wy,t,d,Aders_y);
310     deriver_bs_kieme(nb_point,degre,knots,ptsctr_wz,t,d,Aders_z);
311     double Wders[16];
312     deriver_bs_kieme(nb_point,degre,knots,poids,t,d,Wders);
313     double Bin[256];
314     binomialCoef(Bin, d);
315     #define Bin(i,j) (*(Bin+(i)*(3)+j))
316     for (int k=0;k<=d;k++)
317     {
318     double v_x = Aders_x[k];
319     double v_y = Aders_y[k];
320     double v_z = Aders_z[k];
321     for (int i=1;i<=k;i++)
322     {
323     v_x = v_x - Bin(k-1,i-1)*Wders[i]*CK_x[k-i];
324     v_y = v_y - Bin(k-1,i-1)*Wders[i]*CK_y[k-i];
325     v_z = v_z - Bin(k-1,i-1)*Wders[i]*CK_z[k-i];
326     }
327     CK_x[k] = v_x/Wders[0];
328     CK_y[k] = v_y/Wders[0];
329     CK_z[k] = v_z/Wders[0];
330     }
331     #undef Bin // enlever (i,j)
332     }
333    
334     void ST_B_SPLINE::deriver(double t,double *dxyz)
335     {
336     double CK_x[2];
337     double CK_y[2];
338     double CK_z[2];
339     deriver_kieme(t,1,CK_x,CK_y,CK_z);
340     dxyz[0]=CK_x[1];
341     dxyz[1]=CK_y[1];
342     dxyz[2]=CK_z[1];
343     }
344    
345     void ST_B_SPLINE::deriver_seconde(double t,double *ddxyz,double* dxyz ,double* xyz )
346     {
347     double CK_x[3];
348     double CK_y[3];
349     double CK_z[3];
350     deriver_kieme(t,2,CK_x,CK_y,CK_z);
351     ddxyz[0]=CK_x[2];
352     ddxyz[1]=CK_y[2];
353     ddxyz[2]=CK_z[2];
354     if (dxyz!=NULL)
355     {
356     dxyz[0]=CK_x[1];
357     dxyz[1]=CK_y[1];
358     dxyz[2]=CK_z[1];
359     }
360     if (xyz!=NULL)
361     {
362     xyz[0]=CK_x[0];
363     xyz[1]=CK_y[0];
364     xyz[2]=CK_z[0];
365     }
366     }
367    
368    
369    
370     void ST_B_SPLINE::inverser(double& t,double *xyz,double precision)
371     {
372     int code;
373     int num_point=nb_point;
374     do
375     {
376     code=inverser2(t,xyz,num_point,precision);
377     num_point=num_point*2;
378     }
379     while (code==0);
380     }
381    
382     int ST_B_SPLINE::inverser2(double& t,double *xyz,int num_test,double precision)
383     {
384     double xyz1[3];
385     double dxyz1[3];
386     double ddxyz1[3];
387     double ti;
388     double eps;
389     double tmin=get_tmin();
390     double tmax=get_tmax();
391     OT_VECTEUR_3D Pt(xyz[0],xyz[1],xyz[2]);
392     double distance_ref=1e308;
393     int ref;
394    
395     for (int i=0;i<num_test+1;i++)
396     {
397     double t=tmin+i*1./num_test*(tmax-tmin);
398     evaluer(t,xyz1);
399     OT_VECTEUR_3D Ct(xyz1[0],xyz1[1],xyz1[2]);
400     OT_VECTEUR_3D Distance = Ct-Pt;
401     double longueur=Distance.get_longueur2();
402     if (longueur<distance_ref)
403     {
404     distance_ref=longueur;
405     ref=i;
406     }
407     }
408    
409     double tii=tmin+ref*1./num_test*(tmax-tmin);
410     int compteur=0;
411     do
412     {
413     compteur++;
414     ti=tii;
415     deriver_seconde(ti,ddxyz1,dxyz1,xyz1);
416     OT_VECTEUR_3D Ct(xyz1[0],xyz1[1],xyz1[2]);
417     OT_VECTEUR_3D Ct_deriver(dxyz1[0],dxyz1[1],dxyz1[2]);
418     OT_VECTEUR_3D Ct_deriver_seconde(ddxyz1[0],ddxyz1[1],ddxyz1[2]);
419     OT_VECTEUR_3D Distance = Ct-Pt;
420     tii=ti-Ct_deriver*Distance/(Ct_deriver_seconde*Distance+Ct_deriver.get_longueur2());
421     if (periodique==1)
422     {
423     if (tii<get_tmin()) tii=get_tmax()-(get_tmin()-tii);
424     if (tii>get_tmax()) tii=get_tmin()+(tii-get_tmax());
425     }
426     else
427     {
428     if (tii<get_tmin()) tii=get_tmin();
429     if (tii>get_tmax()) tii=get_tmax();
430     }
431     eps=fabs(tii-ti);
432     if (compteur>500) return 0;
433     }
434     while (eps>precision);
435     t=ti;
436     return 1;
437     }
438    
439     double ST_B_SPLINE::get_tmin()
440     {
441     return knots[0];
442     }
443     double ST_B_SPLINE::get_tmax()
444     {
445     int i=knots.size();
446     return knots[i-1];
447     }
448    
449     double equation_longueur(ST_B_SPLINE& bsp,double t)
450     {
451     double dxyz[3];
452     bsp.deriver(t,dxyz);
453     return sqrt(dxyz[0]*dxyz[0]+dxyz[1]*dxyz[1]+dxyz[2]*dxyz[2]);
454     }
455    
456    
457     double ST_B_SPLINE::get_longueur(double t1,double t2,double precis)
458     {
459     TPL_FONCTION1<double,ST_B_SPLINE,double> longueur_bsp(*this,equation_longueur);
460     return longueur_bsp.integrer_gauss_2(t1,t2);
461     }
462    
463    
464     int ST_B_SPLINE::est_periodique(void)
465     {
466     return periodique;
467     }
468     double ST_B_SPLINE::get_periode(void)
469     {
470     return periode;
471     }
472    
473     int ST_B_SPLINE::get_type_geometrique(TPL_LISTE_ENTITE<double> &param)
474     {
475     for (int i=0;i<nb_point-(degre+1);i++)
476     {
477     param.ajouter(knots[i]);
478     }
479     double xyz[3];
480     for (int i=0;i<nb_point;i++)
481     {
482     xyz[0]=ptsctr[3*i];
483     xyz[1]=ptsctr[3*i+1];
484     xyz[2]=ptsctr[3*i+2];
485     param.ajouter(xyz[0]);
486     param.ajouter(xyz[1]);
487     param.ajouter(xyz[2]);
488     }
489     for (int i=0;i<nb_point;i++)
490     {
491     param.ajouter(poids[i]);
492     }
493     param.ajouter(degre);
494     return MGCo_BSPLINE;
495     }
496    
497     void ST_B_SPLINE::initialiser(ST_GESTIONNAIRE *gest)
498     {
499     int i=indexptsctr.size();
500     ST_POINT* point1=gest->lst_point.getid(indexptsctr[0]);
501     ST_POINT* point2=gest->lst_point.getid(indexptsctr[i-1]);
502     double xyz1[3];
503     double xyz2[3];
504     point1->evaluer(xyz1);
505     point2->evaluer(xyz2);
506     periodique=0;
507     if (OPERATEUR::egal (xyz1[0],xyz2[0],1E-6))
508     {
509     if (OPERATEUR::egal (xyz1[1],xyz2[1],1E-6))
510     {
511     if (OPERATEUR::egal (xyz1[2],xyz2[2],1E-6))
512     periodique=1;
513     }
514     }
515    
516     if (periodique==1)
517     {
518     int i=knots.size();
519     periode=(knots[i-1]-knots[0]);
520     }
521     else periode=0;
522    
523     int nbptsctr=indexptsctr.size();
524     for (int i=0;i<nbptsctr;i++)
525     {
526     ST_POINT* stpoint=gest->lst_point.getid(indexptsctr[i]);
527     double xyz[3];
528     stpoint->evaluer(xyz);
529     ptsctr.insert(ptsctr.end(),xyz[0]);
530     ptsctr.insert(ptsctr.end(),xyz[1]);
531     ptsctr.insert(ptsctr.end(),xyz[2]);
532     }
533     }
534    
535    
536    
537    
538    
539     void ST_B_SPLINE::est_util(ST_GESTIONNAIRE* gest)
540     {
541     util=true;
542     for (int i=0;i<nb_point;i++)
543     gest->lst_point.getid(indexptsctr[i])->est_util(gest);
544     }
545    
546    
547    
548    
549     void ST_B_SPLINE::get_param_NURBS(int& nb_param,TPL_LISTE_ENTITE<double> &param)
550     {
551    
552     // the order stock of the parameters is:
553     // 1/ the variante type =1 (dimension one)
554     // 2/ the order of the spline in two dirction the second direction is zero
555     // 3/ the number of the control points
556     // 4/ the knots vector
557     // 5/ the control points cordinates in homogeneos forms
558    
559     // The first parameter indicate the code access
560     param.ajouter(1);
561     // The follewing two parameters of the list indicate the orders of the net points
562    
563     param.ajouter(degre+1);
564     param.ajouter(0);
565    
566     // The follewing two parameters indicate the number of rows and colons of the control points
567     // respectively to the two parameters directions
568    
569     param.ajouter(nb_point);
570     param.ajouter(0);
571    
572     // this present the knot vector in the u-direction
573    
574     for (unsigned int i=0;i<knots.size();i++)
575     {
576     param.ajouter(knots[i]);
577     }
578    
579     // in the following we construct the control point vector
580    
581     for (unsigned int pt=0;pt<ptsctr.size();pt++)
582     {
583     param.ajouter(ptsctr[pt]);
584     }
585     nb_param=5+knots.size();
586     }