outil/src/ot_fonctions.h

//---------------------------------------------------------------------------

#ifndef ot_fonctionsH
#define ot_fonctionsH
//---------------------------------------------------------------------------

#ifdef WINDOWS_VERSION
#ifdef BUILT_DLL_OUTIL
#define DLLPORTOUTIL __declspec(dllexport)
#else
#define DLLPORTOUTIL __declspec(dllimport)
#endif
#else
#define DLLPORTOUTIL
#endif

#include <math.h>
#include "ot_doubleprecision.h"
#include "ot_mathematique.h"
#include "ot_algorithme_geometrique.h"


#define PI  3.1415926535897932384626433832795

double2 moment (double2 theta,vector<OT_VECTEUR_4DD>& ctrpts,int plan)
{

//double li[]={78.10249676,78.10249676,78.10249676,78.10249676} ;
//double thi[]={39.80557109+10,140.19442891+10,219.80557109+10,320.19442891+10};

    int nb_pts=ctrpts.size();
    double ordre=4.;
    double2 zer=0.;
    double2 sum=0.;
    for (int i=0;i<nb_pts;i++)
    {
        //thi[i]=thi[i]*PI/180;
        //double2 nori=ctrpts[i].norme();
        double2 x=ctrpts[i].x()^2;
        double2 y=ctrpts[i].y()^2;
        double2 z=ctrpts[i].z()^2;
        double2 nori;

        // double2 nori=ctrpts[i].norme();
        double2 cosi=0.;
        double2 sini=0.;
        if (plan==12)
        {
            nori=(x+y)^0.5;
            if (nori!=zer) {
                cosi=ctrpts[i].get_x()/nori;
                sini=ctrpts[i].get_y()/nori;
            }
        }
        if (plan==23)
        {
            nori=(y+z)^0.5;
            if (nori!=zer) {
                cosi=ctrpts[i].get_y()/nori;
                sini=ctrpts[i].get_z()/nori;
            }

        }
        if (plan==13)
        {
            nori=(x+z)^0.5;
            if (nori!=zer) {
                cosi=ctrpts[i].get_x()/nori;
                sini=ctrpts[i].get_z()/nori;
            }
        }
        double2 costh=cos(theta);
        double2 sinth=sin(theta);
        double2 d_theta=sini*costh-cosi*sinth;
        //sum=sum+ pow(li[i],4)*pow(d_theta,4);
        double2 val1=nori^ordre;
        double2 val2=d_theta^ordre;
        double2 val=val1*val2;
        sum=sum+ val;
    }
    return sum;
}


double2  Max(double2 a,double2 b)
{
    return((a<b)? b : a);
}

double2  Min(double2 a,double2 b)
{
    return ((a>b)? b : a);
}

double2 Sign(double2 a,double2 b)
{
    double2 zer=0.;
    double2 c=-1.;
    if (b>=zer)
        return a;
    else
        return c*a;
}

void swap(double2& a,double2& b)
{
    double2 tmp;
    tmp = a;
    a = b;
    b = tmp;
}
void Shift(double2&a,double2&b,double2&c,double2 d)
{
    a=b;
    b=c;
    c=d;
}

void min_Brak(double2 &ax,double2 &bx,double2 &cx,double2 &fa,double2 &fb,double2 &fc,vector<OT_VECTEUR_4DD>& ctrpts,int plan)
{
    static const double2 GOLD=1.618034;
    static const double2 GLIMIT=100.0;
    static const double2 TINY=1.0e-20;
    double2 ulim,u,r,q,fu;
    double2(*func)(double2,vector<OT_VECTEUR_4DD>&,int)=&moment;
    double2 dtheta=0.0001;
    ax=0.;
    bx=ax+dtheta;
    cx=bx+dtheta;
    double2 fva=(*func)(ax,ctrpts,plan);
    double2 fvb=(*func)(bx,ctrpts,plan);
    double2 fvc=(*func)(cx,ctrpts,plan);

    while (fvb> fva)
    {
        ax=bx;
        bx=ax+dtheta;
        fva=(*func)(ax,ctrpts,plan);
        fvb=(*func)(bx,ctrpts,plan);
    }

    if (fvb<fva)
    {
        cx=bx+dtheta;
        while (fvb> fvc)
        {
            bx=cx;
            cx=bx+dtheta;
            fvb=(*func)(bx,ctrpts,plan);
            fvc=(*func)(cx,ctrpts,plan);
        }
    }

    fa=(*func)(ax,ctrpts,plan);
    fb=(*func)(bx,ctrpts,plan);
    if (fb>fa)
    {
        swap(ax,bx);
        swap(fb,fa);
    }
    cx=bx+GOLD*(bx-ax);
    fc=(*func)(cx,ctrpts,plan);
    double2 d=2.;
    double2 zer=0.;
    while (fb>fc)
    {
        r=(bx-ax)*(fb-fc);
        q=(bx-cx)*(fb-fa);
        u=bx-((bx-cx)*q-(bx-ax)*r)/(d*Sign(Max((q-r).get_fabs(),TINY),q-r));//(d*Sign(Max(abs(q-r),TINY),q-r));
        ulim=bx+GLIMIT*(cx-bx);
        if ((bx-u)*(u-cx)>zer)
        {
            fu=(*func)(u,ctrpts,plan);
            if (fu<fc)
            {
                ax=bx;
                bx=u;
                fa=fb;
                fb=fu;
                return;
            }
            else if (fu>fb)
            {
                cx=u;
                fc=fu;
                return;
            }
            u=cx+GOLD*(cx-bx);
            fu=(*func)(u,ctrpts,plan);
        }
        else if ((cx-u)*(u-ulim)>zer)
        {
            fu=(*func)(u,ctrpts,plan);
            if (fu<fc)
            {
                bx=cx;
                cx=u;
                u=cx+GOLD*(cx-bx);
                Shift(fb,fc,fu,(*func)(u,ctrpts,plan));
            }
        }
        else if ((u-ulim)*(ulim-cx)>=zer)
        {
            u=ulim;
            fu=(*func)(u,ctrpts,plan);
        }
        else
        {
            u=cx+GOLD*(cx-bx);
            fu=(*func)(u,ctrpts,plan);
        }
        Shift(ax,bx,cx,u);
        Shift(fa,fb,fc,fu);
    }
}


vector<OT_VECTEUR_4DD> system_axes(vector<OT_VECTEUR_4DD>& ctrpts)
{
    vector<OT_VECTEUR_4DD> syst_axe;
    int nb_points=ctrpts.size();
    vector<OT_VECTEUR_4DD> ctrpts_proj;
    OT_VECTEUR_4DD axe1,axe2,axe3,axe4;
    double2 ax,bx,cx,fa,fb,fc;
    int plan=12;
    double2 z0=ctrpts[0].get_z() ;
    double2 y0=ctrpts[0].get_z() ;
    double2 x0=ctrpts[0].get_z() ;
    int cmpt1=0;
    int cmpt2=0;
    int cmpt3=0;
    for (int i=0;i<nb_points;i++)
    {
        if (ctrpts[i].get_x()==x0) cmpt1++;
        if (ctrpts[i].get_y()==y0) cmpt2++;
        if (ctrpts[i].get_z()==z0) cmpt3++;
    }

    if (cmpt1==nb_points) plan=23;
    if (cmpt2==nb_points) plan=13;
    if (cmpt3==nb_points) plan=12;


    min_Brak(ax,bx,cx,fa,fb,fc,ctrpts,plan);
    double theta= bx.get_x();
    double c=cos(theta);
    double s=sin(theta);
    axe1.change_x(c);
    axe1.change_y(s);
    axe1.change_z(0.);
    axe1.change_w(0.);
    double e1_proj[3];
    double e2_proj[3];
    double normal[3];
    normal[0] =(axe1.get_x()).get_x();
    normal[1] =(axe1.get_y()).get_x();
    normal[2] =(axe1.get_z()).get_x();
    double root[3]={0.,0.,0.};
    OT_ALGORITHME_GEOMETRIQUE algo;
    int pris=0;
    OT_MATRICE_3D P, P_1;

    for (int i=0;i<nb_points;i++)
    {
        double pnt[3];
        double proj_pnt[3];
        pnt[0] =(ctrpts[i].get_x()).get_x();
        pnt[1] =(ctrpts[i].get_y()).get_x();
        pnt[2] =(ctrpts[i].get_z()).get_x();


        algo.Proj3D_Point_Plan (normal, root, pnt, proj_pnt);

        OT_VECTEUR_4DD pt_proj;


        if ((proj_pnt[0]!=0.||proj_pnt[1]!=0.||proj_pnt[2]!=0)&&!pris)
        {
            e1_proj[0]=proj_pnt[0];
            e1_proj[1]=proj_pnt[1];
            e1_proj[2]=proj_pnt[2];
            double norm=sqrt(pow(proj_pnt[0],2)+pow(proj_pnt[1],2)+pow(proj_pnt[2],2));

            e1_proj[0]=e1_proj[0]/norm;
            e1_proj[1]=e1_proj[1]/norm;
            e1_proj[2]=e1_proj[2]/norm;

            e2_proj[0]=normal[1]*e1_proj[2]-normal[2]*e1_proj[1];
            e2_proj[1]=normal[2]*e1_proj[0]-normal[2]*e1_proj[2];
            e2_proj[2]=normal[0]*e1_proj[1]-normal[2]*e1_proj[0];

            OT_VECTEUR_3D V1(e1_proj[0],e2_proj[0],normal[0]);
            OT_VECTEUR_3D V2(e1_proj[1],e2_proj[1],normal[1]);
            OT_VECTEUR_3D V3(e1_proj[2],e2_proj[2],normal[2]);


            P.change_vecteur1(V1);
            P.change_vecteur2(V2);
            P.change_vecteur3(V3);

            P_1=P.inverse();
            pris=1;
        }

        OT_VECTEUR_3D V_PT(proj_pnt[0],proj_pnt[1],proj_pnt[2]);

        OT_VECTEUR_3D V=P_1*V_PT;


        pt_proj.change_x(V[0]);
        pt_proj.change_y(V[1]);
        pt_proj.change_z(V[2]);

        ctrpts_proj.insert(ctrpts_proj.end(),pt_proj);
    }

    min_Brak(ax,bx,cx,fa,fb,fc,ctrpts_proj,plan);

    theta= bx.get_x();
    c=cos(theta);
    s=sin(theta);
    OT_VECTEUR_3D Vp(c,s,0.);
    double   nnv_x=Vp*P.get_vecteur1();
    double   nnv_y=Vp*P.get_vecteur2();
    double   nnv_z=Vp*P.get_vecteur3();

    axe2.change_x(nnv_x);
    axe2.change_y(nnv_y);
    axe2.change_z(nnv_z);

    double2 val=axe1*axe2;

    axe3[0]=axe1[1]*axe2[2]-axe1[2]*axe2[1];
    axe3[1]=axe1[2]*axe2[0]-axe1[2]*axe2[2];
    axe3[2]=axe1[0]*axe2[1]-axe1[2]*axe2[0];

    axe3[3]=0.;


    syst_axe.insert( syst_axe.end(),axe1);
    syst_axe.insert( syst_axe.end(),axe2);
    syst_axe.insert( syst_axe.end(),axe3);

    return syst_axe;
}


#endif
Revision:	283
Committed:	Tue Sep 13 21:11:20 2011 UTC (13 years, 8 months ago) by francois
Content type:	text/plain
File size:	8391 byte(s)
Log Message:	structure de l'écriture
#	User	Rev	Content
1	francois	283	//---------------------------------------------------------------------------
2
3			#ifndef ot_fonctionsH
4			#define ot_fonctionsH
5			//---------------------------------------------------------------------------
6
7			#ifdef WINDOWS_VERSION
8			#ifdef BUILT_DLL_OUTIL
9			#define DLLPORTOUTIL __declspec(dllexport)
10			#else
11			#define DLLPORTOUTIL __declspec(dllimport)
12			#endif
13			#else
14			#define DLLPORTOUTIL
15			#endif
16
17			#include <math.h>
18			#include "ot_doubleprecision.h"
19			#include "ot_mathematique.h"
20			#include "ot_algorithme_geometrique.h"
21
22
23			#define PI 3.1415926535897932384626433832795
24
25			double2 moment (double2 theta,vector<OT_VECTEUR_4DD>& ctrpts,int plan)
26	souaissa	107	{
27	francois	283
28			//double li[]={78.10249676,78.10249676,78.10249676,78.10249676} ;
29			//double thi[]={39.80557109+10,140.19442891+10,219.80557109+10,320.19442891+10};
30
31			int nb_pts=ctrpts.size();
32			double ordre=4.;
33			double2 zer=0.;
34			double2 sum=0.;
35			for (int i=0;i<nb_pts;i++)
36			{
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			{
49			nori=(x+y)^0.5;
50			if (nori!=zer) {
51			cosi=ctrpts[i].get_x()/nori;
52			sini=ctrpts[i].get_y()/nori;
53			}
54			}
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			}
62
63			}
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			}
71			}
72			double2 costh=cos(theta);
73			double2 sinth=sin(theta);
74			double2 d_theta=sinicosth-cosisinth;
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	souaissa	107	}
83	francois	283
84
85			double2 Max(double2 a,double2 b)
86	souaissa	107	{
87	francois	283	return((a<b)? b : a);
88			}
89
90			double2 Min(double2 a,double2 b)
91	souaissa	107	{
92	francois	283	return ((a>b)? b : a);
93			}
94	souaissa	107
95	francois	283	double2 Sign(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 swap(double2& a,double2& b)
106			{
107			double2 tmp;
108			tmp = a;
109			a = b;
110			b = tmp;
111			}
112			void Shift(double2&a,double2&b,double2&c,double2 d)
113			{
114			a=b;
115			b=c;
116			c=d;
117			}
118
119			void min_Brak(double2 &ax,double2 &bx,double2 &cx,double2 &fa,double2 &fb,double2 &fc,vector<OT_VECTEUR_4DD>& ctrpts,int plan)
120			{
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			ax=0.;
128			bx=ax+dtheta;
129			cx=bx+dtheta;
130			double2 fva=(*func)(ax,ctrpts,plan);
131			double2 fvb=(*func)(bx,ctrpts,plan);
132			double2 fvc=(*func)(cx,ctrpts,plan);
133
134			while (fvb> fva)
135			{
136	souaissa	107	ax=bx;
137			bx=ax+dtheta;
138			fva=(*func)(ax,ctrpts,plan);
139			fvb=(*func)(bx,ctrpts,plan);
140	francois	283	}
141	souaissa	107
142	francois	283	if (fvb<fva)
143			{
144			cx=bx+dtheta;
145			while (fvb> fvc)
146	souaissa	107	{
147			bx=cx;
148			cx=bx+dtheta;
149			fvb=(*func)(bx,ctrpts,plan);
150			fvc=(*func)(cx,ctrpts,plan);
151	francois	283	}
152			}
153
154			fa=(*func)(ax,ctrpts,plan);
155			fb=(*func)(bx,ctrpts,plan);
156			if (fb>fa)
157			{
158			swap(ax,bx);
159			swap(fb,fa);
160			}
161			cx=bx+GOLD*(bx-ax);
162			fc=(*func)(cx,ctrpts,plan);
163			double2 d=2.;
164			double2 zer=0.;
165			while (fb>fc)
166			{
167			r=(bx-ax)*(fb-fc);
168			q=(bx-cx)*(fb-fa);
169			u=bx-((bx-cx)q-(bx-ax)r)/(dSign(Max((q-r).get_fabs(),TINY),q-r));//(dSign(Max(abs(q-r),TINY),q-r));
170			ulim=bx+GLIMIT*(cx-bx);
171			if ((bx-u)*(u-cx)>zer)
172			{
173			fu=(*func)(u,ctrpts,plan);
174			if (fu<fc)
175			{
176			ax=bx;
177			bx=u;
178			fa=fb;
179			fb=fu;
180			return;
181	souaissa	107	}
182	francois	283	else if (fu>fb)
183			{
184			cx=u;
185			fc=fu;
186			return;
187			}
188			u=cx+GOLD*(cx-bx);
189			fu=(*func)(u,ctrpts,plan);
190	souaissa	107	}
191	francois	283	else if ((cx-u)*(u-ulim)>zer)
192			{
193			fu=(*func)(u,ctrpts,plan);
194			if (fu<fc)
195			{
196			bx=cx;
197			cx=u;
198			u=cx+GOLD*(cx-bx);
199			Shift(fb,fc,fu,(*func)(u,ctrpts,plan));
200			}
201			}
202			else if ((u-ulim)*(ulim-cx)>=zer)
203			{
204			u=ulim;
205			fu=(*func)(u,ctrpts,plan);
206			}
207			else
208			{
209			u=cx+GOLD*(cx-bx);
210			fu=(*func)(u,ctrpts,plan);
211			}
212			Shift(ax,bx,cx,u);
213			Shift(fa,fb,fc,fu);
214			}
215			}
216	souaissa	107
217	francois	283
218
219
220			vector<OT_VECTEUR_4DD> system_axes(vector<OT_VECTEUR_4DD>& ctrpts)
221			{
222			vector<OT_VECTEUR_4DD> syst_axe;
223			int nb_points=ctrpts.size();
224			vector<OT_VECTEUR_4DD> ctrpts_proj;
225			OT_VECTEUR_4DD axe1,axe2,axe3,axe4;
226			double2 ax,bx,cx,fa,fb,fc;
227			int plan=12;
228			double2 z0=ctrpts[0].get_z() ;
229			double2 y0=ctrpts[0].get_z() ;
230			double2 x0=ctrpts[0].get_z() ;
231			int cmpt1=0;
232			int cmpt2=0;
233			int cmpt3=0;
234			for (int i=0;i<nb_points;i++)
235			{
236			if (ctrpts[i].get_x()==x0) cmpt1++;
237			if (ctrpts[i].get_y()==y0) cmpt2++;
238			if (ctrpts[i].get_z()==z0) cmpt3++;
239			}
240
241			if (cmpt1==nb_points) plan=23;
242			if (cmpt2==nb_points) plan=13;
243			if (cmpt3==nb_points) plan=12;
244
245
246
247			min_Brak(ax,bx,cx,fa,fb,fc,ctrpts,plan);
248			double theta= bx.get_x();
249			double c=cos(theta);
250			double s=sin(theta);
251			axe1.change_x(c);
252			axe1.change_y(s);
253			axe1.change_z(0.);
254			axe1.change_w(0.);
255			double e1_proj[3];
256			double e2_proj[3];
257			double normal[3];
258			normal[0] =(axe1.get_x()).get_x();
259			normal[1] =(axe1.get_y()).get_x();
260			normal[2] =(axe1.get_z()).get_x();
261			double root[3]={0.,0.,0.};
262			OT_ALGORITHME_GEOMETRIQUE algo;
263			int pris=0;
264			OT_MATRICE_3D P, P_1;
265
266			for (int i=0;i<nb_points;i++)
267			{
268			double pnt[3];
269			double proj_pnt[3];
270			pnt[0] =(ctrpts[i].get_x()).get_x();
271			pnt[1] =(ctrpts[i].get_y()).get_x();
272			pnt[2] =(ctrpts[i].get_z()).get_x();
273
274
275			algo.Proj3D_Point_Plan (normal, root, pnt, proj_pnt);
276
277			OT_VECTEUR_4DD pt_proj;
278
279
280			if ((proj_pnt[0]!=0.\|\|proj_pnt[1]!=0.\|\|proj_pnt[2]!=0)&&!pris)
281			{
282			e1_proj[0]=proj_pnt[0];
283			e1_proj[1]=proj_pnt[1];
284			e1_proj[2]=proj_pnt[2];
285			double norm=sqrt(pow(proj_pnt[0],2)+pow(proj_pnt[1],2)+pow(proj_pnt[2],2));
286
287			e1_proj[0]=e1_proj[0]/norm;
288			e1_proj[1]=e1_proj[1]/norm;
289			e1_proj[2]=e1_proj[2]/norm;
290
291			e2_proj[0]=normal[1]e1_proj[2]-normal[2]e1_proj[1];
292			e2_proj[1]=normal[2]e1_proj[0]-normal[2]e1_proj[2];
293			e2_proj[2]=normal[0]e1_proj[1]-normal[2]e1_proj[0];
294
295			OT_VECTEUR_3D V1(e1_proj[0],e2_proj[0],normal[0]);
296			OT_VECTEUR_3D V2(e1_proj[1],e2_proj[1],normal[1]);
297			OT_VECTEUR_3D V3(e1_proj[2],e2_proj[2],normal[2]);
298
299
300
301			P.change_vecteur1(V1);
302			P.change_vecteur2(V2);
303			P.change_vecteur3(V3);
304
305			P_1=P.inverse();
306			pris=1;
307			}
308
309			OT_VECTEUR_3D V_PT(proj_pnt[0],proj_pnt[1],proj_pnt[2]);
310
311			OT_VECTEUR_3D V=P_1*V_PT;
312
313
314			pt_proj.change_x(V[0]);
315			pt_proj.change_y(V[1]);
316			pt_proj.change_z(V[2]);
317
318			ctrpts_proj.insert(ctrpts_proj.end(),pt_proj);
319			}
320
321			min_Brak(ax,bx,cx,fa,fb,fc,ctrpts_proj,plan);
322
323			theta= bx.get_x();
324			c=cos(theta);
325			s=sin(theta);
326			OT_VECTEUR_3D Vp(c,s,0.);
327			double nnv_x=Vp*P.get_vecteur1();
328			double nnv_y=Vp*P.get_vecteur2();
329			double nnv_z=Vp*P.get_vecteur3();
330
331			axe2.change_x(nnv_x);
332			axe2.change_y(nnv_y);
333			axe2.change_z(nnv_z);
334
335			double2 val=axe1*axe2;
336
337			axe3[0]=axe1[1]axe2[2]-axe1[2]axe2[1];
338			axe3[1]=axe1[2]axe2[0]-axe1[2]axe2[2];
339			axe3[2]=axe1[0]axe2[1]-axe1[2]axe2[0];
340
341			axe3[3]=0.;
342
343
344			syst_axe.insert( syst_axe.end(),axe1);
345			syst_axe.insert( syst_axe.end(),axe2);
346			syst_axe.insert( syst_axe.end(),axe3);
347
348			return syst_axe;
349			}
350
351
352
353
354			#endif