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

File Contents

# User Rev Content
1 francois 283 #include "gestionversion.h"
2     //---------------------------------------------------------------------------
3    
4     #pragma hdrstop
5    
6     #include "vct_outils.h"
7     //---------------------------------------------------------------------------
8     #pragma package(smart_init)
9    
10     VCT_OUTILS::VCT_OUTILS(double ord):ordre(ord)
11 souaissa 151 {
12 francois 283
13 souaissa 151 }
14 francois 283
15     double2 VCT_OUTILS::get_moment_dans_plan (double2 theta,vector<OT_VECTEUR_4DD>& ctrpts,int plan)
16 souaissa 151 {
17    
18 francois 283 //double li[]={78.10249676,78.10249676,78.10249676,78.10249676} ;
19     //double thi[]={39.80557109+10,140.19442891+10,219.80557109+10,320.19442891+10};
20     int nb_pts=ctrpts.size();
21 souaissa 151
22 francois 283 //Centrer les points autour du brycentre appartenant au plan;
23 souaissa 151
24 francois 283 //OT_VECTEUR_4DD bary;
25     // for(int i=0;i<nb_pts;i++)
26     // {
27     // bary=bary+ctrpts[i];
28     // }
29     //bary=1./nb_pts*bary;
30     //-----------------------------------
31 souaissa 151
32 francois 283
33     double2 zer=0.;
34     double2 sum=0.;
35     for (int i=0;i<nb_pts;i++)
36     {
37     // ctrpts[i]=ctrpts[i]-bary;
38     //thi[i]=thi[i]*PI/180;
39     //double2 nori=ctrpts[i].norme();
40     double2 x=ctrpts[i].x()^2;
41     double2 y=ctrpts[i].y()^2;
42     double2 z=ctrpts[i].z()^2;
43     double2 nori;
44    
45     // double2 nori=ctrpts[i].norme();
46     double2 cosi=0.;
47     double2 sini=0.;
48     if (plan==12)
49 souaissa 151 {
50 francois 283 nori=(x+y)^0.5;
51     if (nori!=zer) {
52     cosi=ctrpts[i].get_x()/nori;
53     sini=ctrpts[i].get_y()/nori;
54 souaissa 151 }
55 francois 283 }
56     if (plan==23)
57     {
58     nori=(y+z)^0.5;
59     if (nori!=zer) {
60     cosi=ctrpts[i].get_y()/nori;
61     sini=ctrpts[i].get_z()/nori;
62 souaissa 151 }
63    
64 francois 283 }
65     if (plan==13)
66     {
67     nori=(x+z)^0.5;
68     if (nori!=zer) {
69     cosi=ctrpts[i].get_x()/nori;
70     sini=ctrpts[i].get_z()/nori;
71 souaissa 151 }
72 francois 283 }
73     double2 costh=cos(theta);
74     double2 sinth=sin(theta);
75     double2 d_theta=sini*costh-cosi*sinth;
76     //sum=sum+ pow(li[i],4)*pow(d_theta,4);
77     double2 val1=nori^ordre;
78     double2 val2=d_theta^ordre;
79     double2 val=val1*val2;
80     sum=sum+ val;
81     }
82     return sum;
83     }
84 souaissa 151
85    
86 francois 283 double2 VCT_OUTILS:: get_max(double2 a,double2 b)
87     {
88     return((a<b)? b : a);
89     }
90 souaissa 151
91 francois 283 double2 VCT_OUTILS:: get_min(double2 a,double2 b)
92     {
93     return ((a>b)? b : a);
94     }
95 souaissa 151
96 francois 283 double2 VCT_OUTILS::get_signe(double2 a,double2 b)
97     {
98     double2 zer=0.;
99     double2 c=-1.;
100     if (b>=zer)
101     return a;
102     else
103     return c*a;
104     }
105    
106     void VCT_OUTILS::get_permutation(double2& a,double2& b)
107     {
108     double2 tmp;
109     tmp = a;
110     a = b;
111     b = tmp;
112     }
113     void VCT_OUTILS::get_shift(double2&a,double2&b,double2&c,double2 d)
114     {
115     a=b;
116     b=c;
117     c=d;
118     }
119    
120     void VCT_OUTILS::get_minimum_de_Brak(double2 &ax,double2 &bx,double2 &cx,double2 &fa,double2 &fb,double2 &fc,vector<OT_VECTEUR_4DD>& ctrpts,int plan)
121     {
122     static const double2 GOLD=1.618034;
123     static const double2 GLIMIT=100.0;
124     static const double2 TINY=1.0e-20;
125     double2 ulim,u,r,q,fu;
126     //double2(*func)(double2,vector<OT_VECTEUR_4DD>&,int)=&moment;
127     double2 dtheta=0.0001;
128    
129     //aptimisation de la valeur de theta
130    
131     int idx_pts_depart1=-1;
132     int idx_pts_depart2=-1;
133     int nb_pts=ctrpts.size();
134     double2 min_y;
135     double2 min_x;
136     int cfait1=-1;
137     int cfait2=-1;
138     int indx1,indx2;
139     if (plan==12) {
140     indx1=0;
141     indx2=1;
142     }
143     if (plan==23) {
144     indx1=1;
145     indx2=2;
146     }
147     if (plan==13) {
148     indx1=0;
149     indx2=2;
150     }
151    
152     OT_VECTEUR_4DD bary;
153     for (int i=0;i<nb_pts;i++)
154     {
155     bary=bary+ctrpts[i];
156     }
157     bary=1./nb_pts*bary;
158     //-----------------------------------
159     /*
160    
161     for (int i=0;i<nb_pts;i++)
162     {
163     ctrpts[i]=ctrpts[i]-bary;
164     OT_VECTEUR_4DD pt=ctrpts[i];
165     if (pt[indx1]>dtheta&&pt[indx2]>dtheta)
166     {
167     if (cfait1<0)
168     {
169     min_y=pt[indx2];
170     idx_pts_depart1=i;
171     cfait1=1;
172     }
173     if(pt[indx2]<min_y){
174     min_y=pt[indx2];
175     idx_pts_depart1=i;
176     }
177     }
178    
179     if (pt[indx1]>dtheta&&pt[indx2]>dtheta)
180     {
181     if(cfait2<0)
182     {
183     min_x=pt[indx1];
184     idx_pts_depart2=i;
185     cfait2=1;
186     }
187     if(pt[indx1]<min_x){
188     min_x=pt[indx1];
189     idx_pts_depart2=i;
190     }
191     }
192    
193     }
194     int idx_pts_depart;
195     if(idx_pts_depart1>=0) idx_pts_depart= idx_pts_depart1;
196     else idx_pts_depart= idx_pts_depart2;
197     OT_VECTEUR_4DD dept=ctrpts[idx_pts_depart];
198    
199     double2 theta_depart=atan(dept[indx2]/dept[indx1]);
200    
201     ax=theta_depart;//0.; */
202     ax=0.;
203     bx=ax+dtheta;
204     cx=bx+dtheta;
205     double2 fva=get_moment_dans_plan(ax,ctrpts,plan);//(*func)(ax,ctrpts,plan);
206     double2 fvb=get_moment_dans_plan(bx,ctrpts,plan);//(*func)(bx,ctrpts,plan);
207     double2 fvc=get_moment_dans_plan(cx,ctrpts,plan);//(*func)(cx,ctrpts,plan);
208     cx=cx+dtheta;
209     double2 fvc_plus_dth=get_moment_dans_plan(cx,ctrpts,plan);
210    
211     while (fvb> fva)
212     {
213 souaissa 151 ax=bx;
214     bx=ax+dtheta;
215     fva=get_moment_dans_plan(ax,ctrpts,plan);//(*func)(ax,ctrpts,plan);
216     fvb=get_moment_dans_plan(bx,ctrpts,plan);//(*func)(bx,ctrpts,plan);
217 francois 283 }
218 souaissa 151
219 francois 283 if (fvb<fva)
220     {
221     cx=bx+dtheta;
222     while (fvb> fvc)
223 souaissa 151 {
224     bx=cx;
225     cx=bx+dtheta;
226     fvb=get_moment_dans_plan(bx,ctrpts,plan);//(*func)(bx,ctrpts,plan);
227     fvc=get_moment_dans_plan(cx,ctrpts,plan);//(*func)(cx,ctrpts,plan);
228 francois 283 }
229     }
230    
231     fa=get_moment_dans_plan(ax,ctrpts,plan);//(*func)(ax,ctrpts,plan);
232     fb=get_moment_dans_plan(bx,ctrpts,plan);//(*func)(bx,ctrpts,plan);
233     if (fb>fa)
234     {
235     swap(ax,bx);
236     swap(fb,fa);
237     }
238     cx=bx+GOLD*(bx-ax);
239     fc=get_moment_dans_plan(cx,ctrpts,plan);//(*func)(cx,ctrpts,plan);
240     double2 d=2.;
241     double2 zer=0.;
242     while (fb>fc)
243     {
244     r=(bx-ax)*(fb-fc);
245     q=(bx-cx)*(fb-fa);
246     u=bx-((bx-cx)*q-(bx-ax)*r)/(d*get_signe(get_max((q-r).get_fabs(),TINY),q-r));//(d*Sign(Max(abs(q-r),TINY),q-r));
247     ulim=bx+GLIMIT*(cx-bx);
248     if ((bx-u)*(u-cx)>zer)
249     {
250     fu=get_moment_dans_plan(u,ctrpts,plan);//(*func)(u,ctrpts,plan);
251     if (fu<fc)
252     {
253     ax=bx;
254     bx=u;
255     fa=fb;
256     fb=fu;
257     return;
258 souaissa 151 }
259 francois 283 else if (fu>fb)
260     {
261     cx=u;
262     fc=fu;
263     return;
264     }
265     u=cx+GOLD*(cx-bx);
266     fu=get_moment_dans_plan(u,ctrpts,plan);
267 souaissa 151 }
268 francois 283 else if ((cx-u)*(u-ulim)>zer)
269     {
270     fu=get_moment_dans_plan(u,ctrpts,plan);
271     if (fu<fc)
272     {
273     bx=cx;
274     cx=u;
275     u=cx+GOLD*(cx-bx);
276     get_shift(fb,fc,fu,get_moment_dans_plan(u,ctrpts,plan));
277     }
278     }
279     else if ((u-ulim)*(ulim-cx)>=zer)
280     {
281     u=ulim;
282     fu=get_moment_dans_plan(u,ctrpts,plan);
283     }
284     else
285     {
286     u=cx+GOLD*(cx-bx);
287     fu=get_moment_dans_plan(u,ctrpts,plan);
288     }
289     get_shift(ax,bx,cx,u);
290     get_shift(fa,fb,fc,fu);
291     }
292     }
293 souaissa 151
294 francois 283
295    
296    
297     vector<OT_VECTEUR_4DD> VCT_OUTILS::get_system_axes(vector<OT_VECTEUR_4DD>& ctrpts)
298     {
299     double e1_proj[3];
300     double e2_proj[3];
301     double normal_palan_project[3];
302     double normal_au_plan[3];
303     double theta,c,s;
304     double root[3]={0.,0.,0.};
305     OT_ALGORITHME_GEOMETRIQUE algo;
306     int pris=0;
307     OT_MATRICE_3D P, P_1;
308    
309     vector<OT_VECTEUR_4DD> syst_axe;
310     int nb_points=ctrpts.size();
311     vector<OT_VECTEUR_4DD> ctrpts_proj;
312     OT_VECTEUR_4DD axe1,axe2,axe3,axe4;
313     double2 ax,bx,cx,fa,fb,fc;
314     int plan=12;
315    
316     double2 z0=ctrpts[0].get_z() ;
317     double2 y0=ctrpts[0].get_y() ;
318     double2 x0=ctrpts[0].get_x() ;
319    
320     int cmpt1=0;
321     int cmpt2=0;
322     int cmpt3=0;
323     for (int i=0;i<nb_points;i++)
324     {
325     if (ctrpts[i].get_x()==x0) cmpt1++;
326     if (ctrpts[i].get_y()==y0) cmpt2++;
327     if (ctrpts[i].get_z()==z0) cmpt3++;
328     }
329    
330     if (cmpt1==nb_points) {
331     plan=23;
332     normal_au_plan[0]=1.;
333     normal_au_plan[1]=0.;
334     normal_au_plan[2]=0.;
335     normal_au_plan[3]=0.;
336     }
337     if (cmpt2==nb_points) {
338     plan=13;
339     normal_au_plan[0]=0.;
340     normal_au_plan[1]=1.;
341     normal_au_plan[2]=0.;
342     normal_au_plan[3]=0.;
343     }
344     if (cmpt3==nb_points) {
345     plan=12;
346     normal_au_plan[0]=0.;
347     normal_au_plan[1]=0.;
348     normal_au_plan[2]=1.;
349     normal_au_plan[3]=0.;
350     }
351    
352    
353     get_minimum_de_Brak(ax,bx,cx,fa,fb,fc,ctrpts,plan);
354    
355     theta= bx.get_x();
356     if (theta==1e-6) theta=0.;
357     c=cos(theta);
358     s=sin(theta);
359    
360     if (plan==12)
361     {
362     axe1.change_x(c);
363     axe1.change_y(s);
364     axe1.change_z(0.);
365     axe1.change_w(0.);
366    
367     axe2.change_x(-s);
368     axe2.change_y(c);
369     axe2.change_z(0);
370     axe2.change_w(0);
371     }
372     if (plan==13)
373     {
374     axe1.change_x(c);
375     axe1.change_y(0.);
376     axe1.change_z(s);
377     axe1.change_w(0.);
378    
379     axe2.change_x(-s);
380     axe2.change_y(0);
381     axe2.change_z(c);
382     axe2.change_w(0);
383     }
384     if (plan==23)
385     {
386     axe1.change_x(0);
387     axe1.change_y(c);
388     axe1.change_z(s);
389     axe1.change_w(0.);
390    
391     axe2.change_x(0);
392     axe2.change_y(-s);
393     axe2.change_z(c);
394     axe2.change_w(0);
395     }
396    
397     if (cmpt1==nb_points||cmpt2==nb_points||cmpt3==nb_points)
398     {
399     // axe2[0]=normal_au_plan[1]*axe1[2]-normal_au_plan[2]*axe1[1];
400     // axe2[1]=normal_au_plan[2]*axe1[0]-normal_au_plan[2]*axe1[2];
401     // axe2[2]=normal_au_plan[0]*axe1[1]-normal_au_plan[2]*axe1[0];
402    
403     axe3[0]=normal_au_plan[0];
404     axe3[1]=normal_au_plan[1];
405     axe3[2]=normal_au_plan[2];
406    
407     syst_axe.insert( syst_axe.end(),axe1);
408     syst_axe.insert( syst_axe.end(),axe2);
409     syst_axe.insert( syst_axe.end(),axe3);
410    
411     return syst_axe;
412     }
413    
414     normal_palan_project[0] =s;
415     normal_palan_project[1] =-c;
416     normal_palan_project[2] =0.;
417    
418    
419     for (int i=0;i<nb_points;i++)
420     {
421     OT_VECTEUR_4DD pt_proj;
422     double pnt[3];
423     double proj_pnt[3];
424    
425     pnt[0] =(ctrpts[i].get_x()).get_x();
426     pnt[1] =(ctrpts[i].get_y()).get_x();
427     pnt[2] =(ctrpts[i].get_z()).get_x();
428    
429     algo.Proj3D_Point_Plan (normal_palan_project, root, pnt, proj_pnt);
430    
431     if ((proj_pnt[0]!=0.||proj_pnt[1]!=0.||proj_pnt[2]!=0)&&!pris)
432     {
433     e1_proj[0]=proj_pnt[0];
434     e1_proj[1]=proj_pnt[1];
435     e1_proj[2]=proj_pnt[2];
436     double norm=sqrt(pow(proj_pnt[0],2)+pow(proj_pnt[1],2)+pow(proj_pnt[2],2));
437    
438     e1_proj[0]=e1_proj[0]/norm;
439     e1_proj[1]=e1_proj[1]/norm;
440     e1_proj[2]=e1_proj[2]/norm;
441    
442     e2_proj[0]=normal_palan_project[1]*e1_proj[2]-normal_palan_project[2]*e1_proj[1];
443     e2_proj[1]=normal_palan_project[2]*e1_proj[0]-normal_palan_project[0]*e1_proj[2];
444     e2_proj[2]=normal_palan_project[0]*e1_proj[1]-normal_palan_project[1]*e1_proj[0];
445    
446     OT_VECTEUR_3D V1(e1_proj[0],e2_proj[0],normal_palan_project[0]);
447     OT_VECTEUR_3D V2(e1_proj[1],e2_proj[1],normal_palan_project[1]);
448     OT_VECTEUR_3D V3(e1_proj[2],e2_proj[2],normal_palan_project[2]);
449    
450     P.change_vecteur1(V1);
451     P.change_vecteur2(V2);
452     P.change_vecteur3(V3);
453    
454     P_1=P.inverse();
455     pris=1;
456     }
457    
458     OT_VECTEUR_3D V_PT(proj_pnt[0],proj_pnt[1],proj_pnt[2]);
459    
460     OT_VECTEUR_3D V=P_1*V_PT;
461    
462    
463     pt_proj.change_x(V[0]);
464     pt_proj.change_y(V[1]);
465     pt_proj.change_z(V[2]);
466    
467     ctrpts_proj.insert(ctrpts_proj.end(),pt_proj);
468     }
469    
470     get_minimum_de_Brak(ax,bx,cx,fa,fb,fc,ctrpts_proj,plan);
471    
472     theta= bx.get_x();
473     c=cos(theta);
474     s=sin(theta);
475     OT_VECTEUR_3D Vp(c,s,0.);
476     double nnv_x=Vp*P.get_vecteur1();
477     double nnv_y=Vp*P.get_vecteur2();
478     double nnv_z=Vp*P.get_vecteur3();
479    
480     axe2.change_x(nnv_x);
481     axe2.change_y(nnv_y);
482     axe2.change_z(nnv_z);
483    
484     double2 val=axe1*axe2;
485    
486     axe3[0]=axe1[1]*axe2[2]-axe1[2]*axe2[1];
487     axe3[1]=axe1[2]*axe2[0]-axe1[2]*axe2[2];
488     axe3[2]=axe1[0]*axe2[1]-axe1[2]*axe2[0];
489    
490     axe3[3]=0.;
491    
492    
493     syst_axe.insert( syst_axe.end(),axe1);
494     syst_axe.insert( syst_axe.end(),axe2);
495     syst_axe.insert( syst_axe.end(),axe3);
496    
497     return syst_axe;
498     }
499    
500    
501    
502    
503     //vector<OT_VECTEUR_4DD> VCT_OUTILS::get_system_axes(OT_VECTEUR_4DD& root,vector<OT_VECTEUR_4DD>& axes,vector<OT_VECTEUR_4DD>& ctrpts)
504     OT_MATRICE_3D VCT_OUTILS::get_system_axes(OT_VECTEUR_4DD& root,vector<OT_VECTEUR_4DD>& axes,vector<OT_VECTEUR_4DD>& ctrpts)
505     {
506    
507     double normal_palan_project[3];
508     double origine[3];
509     double2 ax,bx,cx,fa,fb,fc;
510     double theta,cosin,sinus;
511     OT_MATRICE_3D p,qT, pT,q;
512     vector<OT_VECTEUR_4DD> ctrpts_proj;
513     vector<OT_VECTEUR_4DD> syst_axe;
514     int nb_points=ctrpts.size();
515     OT_ALGORITHME_GEOMETRIQUE algo;
516    
517     normal_palan_project[0] =(axes[2].get_x()).get_x();
518     normal_palan_project[1] =(axes[2].get_y()).get_x();
519     normal_palan_project[2] =(axes[2].get_z()).get_x();
520     origine[0] =root[0].get_x();
521     origine[1] =root[1].get_x();
522     origine[2] =root[2].get_x();
523    
524     OT_VECTEUR_3D V1((axes[0].get_x()).get_x(),(axes[0].get_y()).get_x(),(axes[0].get_z()).get_x());
525     OT_VECTEUR_3D V2((axes[1].get_x()).get_x(),(axes[1].get_y()).get_x(),(axes[1].get_z()).get_x());
526     OT_VECTEUR_3D V3((axes[2].get_x()).get_x(),(axes[2].get_y()).get_x(),(axes[2].get_z()).get_x());
527    
528     p.change_vecteur1(V1);
529     p.change_vecteur2(V2);
530     p.change_vecteur3(V3);
531     pT=p.transpose();
532    
533     for (int i=0;i<nb_points;i++)
534     {
535     OT_VECTEUR_4DD pt_proj;
536     double point[3];
537     double proj_pnt[3];
538    
539    
540     point[0] =(ctrpts[i].get_x()).get_x();
541     point[1] =(ctrpts[i].get_y()).get_x();
542     point[2] =(ctrpts[i].get_z()).get_x();
543    
544     algo.Proj3D_Point_Plan (normal_palan_project, origine, point, proj_pnt);
545     //centrer les points de controles sur le barycentre
546     for (int i=0;i<3;i++)
547     {
548     proj_pnt[i]-= origine[i];
549     }
550    
551     OT_VECTEUR_3D V_PT(proj_pnt[0],proj_pnt[1],proj_pnt[2]);
552     //chercher les coord des points dans le nouveau systeme attach/ au plan
553    
554     OT_VECTEUR_3D V=pT*V_PT;
555    
556    
557     pt_proj.change_x(V[0]);
558     pt_proj.change_y(V[1]);
559     pt_proj.change_z(V[2]);
560    
561     ctrpts_proj.insert(ctrpts_proj.end(),pt_proj);
562    
563     }
564    
565     get_minimum_de_Brak(ax,bx,cx,fa,fb,fc,ctrpts_proj,12);
566    
567     theta= bx.get_x();
568     // if (theta<=0.1) theta=0.;
569     cosin=cos(theta);
570     sinus=sin(theta);
571    
572     OT_VECTEUR_4DD e1,e2,e3;
573    
574     e1.change_x(cosin);
575     e1.change_y(sinus);
576     e1.change_z(0.);
577     e1.change_w(0.);
578    
579     e2.change_x(-sinus);
580     e2.change_y(cosin);
581     e2.change_z(0);
582     e2.change_w(0);
583     //e3=axes[2];
584     e3=e1^e2;
585    
586    
587     // e3=e1^e2;
588    
589    
590     V1.change_x(cosin);
591     V1.change_y(sinus);
592     V1.change_z(0);
593    
594     V2.change_x((e2[0].get_x()));
595     V2.change_y((e2[1].get_x()));
596     V2.change_z((e2[2].get_x()));
597    
598     V3.change_x((e3[0].get_x()));
599     V3.change_y((e3[1].get_x()));
600     V3.change_z((e3[2].get_x()));
601     q.change_vecteur1(V1);
602     q.change_vecteur2(V2);
603     q.change_vecteur3(V3);
604     qT=q.transpose();
605     syst_axe.insert( syst_axe.end(),e1);
606     syst_axe.insert( syst_axe.end(),e2);
607     syst_axe.insert( syst_axe.end(),e3);
608    
609     OT_MATRICE_3D w=p*q; //les vecteurs de la base est en colone
610     //w=w.inverse();
611     // OT_MATRICE_3D w(V1,V2,V3);
612    
613     return w;
614     // return syst_axe;
615     }
616    
617    
618    
619     OT_VECTEUR_4DD VCT_OUTILS::get_plan_medien(vector<OT_VECTEUR_4DD>& ctrpts)
620     {
621     int nb_pts=ctrpts.size();
622     OT_MATRICE_3D p;
623     OT_VECTEUR_3D v;
624     for (int i=0;i<nb_pts;i++)
625     {
626     OT_VECTEUR_4DD pt=ctrpts[i];
627     p(0,0)=p(0,0)+pt[0].get_x()*pt[0].get_x();
628     p(0,1)=p(0,0)+pt[0].get_x()*pt[1].get_x();
629     p(0,2)=p(0,0)+pt[0].get_x();
630     p(1,1)=p(0,0)+pt[1].get_x()*pt[1].get_x();
631     p(1,2)=p(0,0)+pt[1].get_x();
632     v[0]=v[0]+pt[0].get_x()*pt[2].get_x();
633     v[1]=v[1]+pt[1].get_x()*pt[2].get_x();
634     v[2]=v[1]+pt[2].get_x();
635     }
636     p(1,0)=p(0,1);
637     p(2,0)=p(0,2);
638     p(2,1)=p(2,1);
639    
640     p=p.inverse();
641     OT_VECTEUR_3D param=p*v;
642     OT_VECTEUR_4DD tparam(param[0],param[1],1,0.);
643     tparam=tparam.vecteur_norme();
644     double2 nor= tparam.norme();
645     double2 prec=1e-5 ;
646     for (int j=0;j<4;j++)
647     {
648     if (tparam[j].get_fabs()<prec)
649     tparam[j]=0.;
650     }
651     tparam[3]=tparam[3]/nor;
652     return tparam;
653    
654     }