ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/outil/src/vct_outils.cpp
Revision: 160
Committed: Thu Sep 18 15:19:08 2008 UTC (16 years, 7 months ago) by souaissa
Original Path: magic/lib/outil/outil/src/vct_outils.cpp
File size: 15317 byte(s)
Log Message:
mise a jour de vct_outils ot_fonctions

File Contents

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