ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/outil/src/vct_outils.cpp
Revision: 1019
Committed: Tue Jun 4 21:16:50 2019 UTC (6 years ago) by francois
File size: 16369 byte(s)
Log Message: 
restructuration de magic
outil est sorti de lib pour pouvoir etre utiliser en dehors de lib
template est merge avec outil
poly_occ et un sous projet de magic qui utilise le nouveau outil

File Contents

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