ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mg_lissage.cpp
Revision: 253
Committed: Tue Jul 13 19:40:46 2010 UTC (14 years, 10 months ago) by francois
File size: 55836 byte(s)
Log Message:
changement de hiearchie et utilisation de ccmake + mise a jour

File Contents

# User Rev Content
1 francois 222 #include "gestionversion.h"
2     #include "mg_lissage.h"
3     #include "fem_solution.h"
4     #include "fem_maillage.h"
5     #include "mg_maillage.h"
6     #include "mg_gestionnaire.h"
7     #include "mg_triangle_peau.h"
8     #include <stdio.h>
9 picher 230 #include <stdlib.h>
10     #include <time.h>
11 francois 222 #include <string.h>
12 francois 224 #include <math.h>
13 francois 222 //---------------------------------------------------------------------------
14    
15    
16 francois 224
17 francois 232 MG_LISSAGE::MG_LISSAGE():affichageactif(0)
18 francois 222 {
19 francois 224 for (int i=0;i<9;i++) etat[i]=0;
20 francois 222 }
21    
22    
23     MG_LISSAGE::~MG_LISSAGE()
24     {
25 francois 224 for (int i=0;i<lst_peau.size();i++)
26     delete lst_peau[i];
27 francois 222 }
28 picher 233
29 francois 232 void MG_LISSAGE::active_affichage(void (*fonc)(char*))
30     {
31     affiche=fonc;
32     affichageactif=1;
33     }
34 francois 222
35 picher 233
36 francois 224 void MG_LISSAGE::conserve(int origine)
37     {
38     etat[(origine-1000)/10]=1;
39     }
40    
41 picher 233
42     void MG_LISSAGE::copieorigine(FEM_MAILLAGE* mai,MG_GESTIONNAIRE& gest2)
43 picher 231 {
44 picher 233 LISTE_FEM_TETRA::iterator it_tetra;
45     for (FEM_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
46     {
47     MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
48     mgtet->change_origine(mgtet->get_origine());
49     }
50     }
51     void MG_LISSAGE::gain_poids(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
52     {
53     LISTE_MG_TETRA::iterator it_tetra;
54     double volume_initial = 0.;
55     double volume_final = 0.;
56     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
57     {
58     int originetet=mgtet->get_origine();
59     MG_NOEUD* noeud1 = mgtet->get_noeud1();
60     double x1 = noeud1->get_x();
61     double y1 = noeud1->get_y();
62     double z1 = noeud1->get_z();
63     MG_NOEUD* noeud2 = mgtet->get_noeud2();
64     double x2 = noeud2->get_x();
65     double y2 = noeud2->get_y();
66     double z2 = noeud2->get_z();
67     MG_NOEUD* noeud3 = mgtet->get_noeud3();
68     double x3 = noeud3->get_x();
69     double y3 = noeud3->get_y();
70     double z3 = noeud3->get_z();
71     MG_NOEUD* noeud4 = mgtet->get_noeud4();
72     double x4 = noeud4->get_x();
73     double y4 = noeud4->get_y();
74     double z4 = noeud4->get_z();
75     double volume_tet = (x2-x1)*((y3-y1)*(z4-z1)-(z3-z1)*(y4-y1))-(y2-y1)*((x3-x1)*(z4-z1)-(z3-z1)*(x4-x1)) + (z2-z1)*((x3-x1)*(y4-y1)-(y3-y1)*(x4-x1));
76     volume_tet = volume_tet/6.;
77     volume_initial = volume_initial + volume_tet;
78     if (originetet != MAILLEUR_AUTO)
79     {
80     volume_final = volume_final + volume_tet;
81     }
82     }
83     cout << "volume_initial = " << volume_initial << endl;
84     cout << "volume_final = " << volume_final << endl;
85     }
86     void MG_LISSAGE::reactivation(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
87     {
88     cout << "Reactivation d'elements manquants sur aretes vives et faces planes" << endl;
89     LISTE_MG_TETRA::iterator it_tetra;
90     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
91     {
92     int originetet=mgtet->get_origine();
93     int compteur = 0;
94     if (originetet == MAILLEUR_AUTO)
95     {
96     MG_TRIANGLE* tri1 = mgtet->get_triangle1();
97     MG_TRIANGLE* tri2 = mgtet->get_triangle2();
98     MG_TRIANGLE* tri3 = mgtet->get_triangle3();
99     MG_TRIANGLE* tri4 = mgtet->get_triangle4();
100     int nb_tet1 = tri1->get_lien_tetra()->get_nb();
101     for (int k = 0;k<nb_tet1;k++)
102     {
103     MG_TETRA* tet1 = tri1->get_lien_tetra()->get(k);
104     int ori1 = tet1->get_origine();
105     if ((tet1 != mgtet) && (ori1 == MAILLEUR_AUTO)) compteur++;
106     }
107     int nb_tet2 = tri2->get_lien_tetra()->get_nb();
108     for (int k = 0;k<nb_tet2;k++)
109     {
110     MG_TETRA* tet2 = tri2->get_lien_tetra()->get(k);
111     int ori2 = tet2->get_origine();
112     if ((tet2 != mgtet) && (ori2 == MAILLEUR_AUTO)) compteur++;
113     }
114     int nb_tet3 = tri3->get_lien_tetra()->get_nb();
115     for (int k = 0;k<nb_tet3;k++)
116     {
117     MG_TETRA* tet3 = tri3->get_lien_tetra()->get(k);
118     int ori3 = tet3->get_origine();
119     if ((tet3 != mgtet) && (ori3 == MAILLEUR_AUTO)) compteur++;
120     }
121     int nb_tet4 = tri4->get_lien_tetra()->get_nb();
122     for (int k = 0;k<nb_tet4;k++)
123     {
124     MG_TETRA* tet4 = tri4->get_lien_tetra()->get(k);
125     int ori4 = tet4->get_origine();
126     if ((tet4 != mgtet) && (ori4 == MAILLEUR_AUTO)) compteur++;
127     }
128 picher 248 if (compteur == 0) mgtet->change_origine(OPTIMISE);
129 picher 233 }
130     }
131     }
132     void MG_LISSAGE::visualisation2(MG_MAILLAGE* mai2,MG_GESTIONNAIRE& gest2,char *nomfichier,int id)
133     {
134     MG_MAILLAGE* mai3 = new MG_MAILLAGE(NULL);
135     MG_GESTIONNAIRE gest3;
136     gest3.ajouter_mg_maillage(mai3);
137     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
138     MG_NOEUD* noeud1 = mai2->get_mg_noeudid(id);
139     int nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
140     for (int j = 0;j<nb_voisins1;j++)
141     {
142     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
143     liste_voisins.ajouter(mgtri_1);
144     MG_NOEUD* no1 = mgtri_1->get_noeud1();
145     int nb_vois1 = no1->get_lien_triangle()->get_nb();
146     for (int k = 0;k<nb_vois1;k++)
147     {
148     MG_TRIANGLE_PEAU* mgtri1 = (MG_TRIANGLE_PEAU*) no1->get_lien_triangle()->get(k);
149     liste_voisins.ajouter(mgtri1);
150     }
151     MG_NOEUD* no2 = mgtri_1->get_noeud2();
152     int nb_vois2 = no2->get_lien_triangle()->get_nb();
153     for (int k = 0;k<nb_vois2;k++)
154     {
155     MG_TRIANGLE_PEAU* mgtri2 = (MG_TRIANGLE_PEAU*) no2->get_lien_triangle()->get(k);
156     liste_voisins.ajouter(mgtri2);
157     }
158     MG_NOEUD* no3 = mgtri_1->get_noeud3();
159     int nb_vois3 = no3->get_lien_triangle()->get_nb();
160     for (int k = 0;k<nb_vois3;k++)
161     {
162     MG_TRIANGLE_PEAU* mgtri3 = (MG_TRIANGLE_PEAU*) no3->get_lien_triangle()->get(k);
163     liste_voisins.ajouter(mgtri3);
164     }
165     }
166     int nb_voisins = liste_voisins.get_nb();
167     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
168     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
169     {
170     MG_NOEUD* noeud1 = mgtri_j->get_noeud1();
171     MG_NOEUD* noeud2 = mgtri_j->get_noeud2();
172     MG_NOEUD* noeud3 = mgtri_j->get_noeud3();
173     mai3->ajouter_mg_noeud(noeud1);
174     mai3->ajouter_mg_noeud(noeud2);
175     mai3->ajouter_mg_noeud(noeud3);
176     mai3->ajouter_mg_triangle(mgtri_j);
177     }
178     gest3.enregistrer(nomfichier);
179     }
180    
181     void MG_LISSAGE::visualisation(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2,char *nomfichier)
182     {
183     //Routine pour visualiser le traitement des cas de non-manifold
184     MG_MAILLAGE* mai3 = new MG_MAILLAGE(NULL);
185     MG_GESTIONNAIRE gest3;
186     gest3.ajouter_mg_maillage(mai3);
187     MG_NOEUD* noeud_nm1 = mg_mai->get_mg_noeudid(41526);//placer ici le ID d'un noeud non manifold
188 francois 234 //MG_NOEUD* noeud_nm2 = mg_mai->get_mg_noeudid(12437);//placer ici le ID d'un second noeud non manifold (pour une ar�te)
189 picher 233 TPL_MAP_ENTITE<MG_TETRA*> liste_voisins;
190     int nb_voisins1 = noeud_nm1->get_lien_tetra()->get_nb();
191     for (int j = 0;j<nb_voisins1;j++)
192     {
193     MG_TETRA* mgtet_1 = (MG_TETRA*) noeud_nm1->get_lien_tetra()->get(j);
194     int origine = mgtet_1->get_origine();
195     if ((origine == OPTIMISE) || (origine == IMPOSE))
196     {
197     liste_voisins.ajouter(mgtet_1);
198     MG_NOEUD* no1 = mgtet_1->get_noeud1();
199     int nb_vois1 = no1->get_lien_tetra()->get_nb();
200     for (int k = 0;k<nb_vois1;k++)
201     {
202     MG_TETRA* mgtet1 = (MG_TETRA*) no1->get_lien_tetra()->get(k);
203     int origine1 = mgtet1->get_origine();
204     if ((origine1 == OPTIMISE) || (origine1 == IMPOSE)) liste_voisins.ajouter(mgtet1);
205     }
206     MG_NOEUD* no2 = mgtet_1->get_noeud2();
207     int nb_vois2 = no2->get_lien_tetra()->get_nb();
208     for (int k = 0;k<nb_vois2;k++)
209     {
210     MG_TETRA* mgtet2 = (MG_TETRA*) no2->get_lien_tetra()->get(k);
211     int origine2 = mgtet2->get_origine();
212     if ((origine2 == OPTIMISE) || (origine2 == IMPOSE)) liste_voisins.ajouter(mgtet2);
213     }
214     MG_NOEUD* no3 = mgtet_1->get_noeud3();
215     int nb_vois3 = no3->get_lien_tetra()->get_nb();
216     for (int k = 0;k<nb_vois3;k++)
217     {
218     MG_TETRA* mgtet3 = (MG_TETRA*) no3->get_lien_tetra()->get(k);
219     int origine3 = mgtet3->get_origine();
220     if ((origine3 == OPTIMISE) || (origine3 == IMPOSE)) liste_voisins.ajouter(mgtet3);
221     }
222     MG_NOEUD* no4 = mgtet_1->get_noeud4();
223     int nb_vois4 = no4->get_lien_tetra()->get_nb();
224     for (int k = 0;k<nb_vois4;k++)
225     {
226     MG_TETRA* mgtet4 = (MG_TETRA*) no4->get_lien_tetra()->get(k);
227     int origine4 = mgtet4->get_origine();
228     if ((origine4 == OPTIMISE) || (origine4 == IMPOSE)) liste_voisins.ajouter(mgtet4);
229     }
230     }
231     }
232     // int nb_voisins2 = noeud_nm2->get_lien_tetra()->get_nb();
233     // for (int j = 0;j<nb_voisins2;j++)
234     // {
235     // MG_TETRA* mgtet_1 = (MG_TETRA*) noeud_nm2->get_lien_tetra()->get(j);
236     // int origine = mgtet_1->get_origine();
237     // if ((origine == OPTIMISE) || (origine == IMPOSE))
238     // {
239     // liste_voisins.ajouter(mgtet_1);
240     // // MG_NOEUD* no1 = mgtet_1->get_noeud1();
241     // // int nb_vois1 = no1->get_lien_tetra()->get_nb();
242     // // for (int k = 0;k<nb_vois1;k++)
243     // // {
244     // // MG_TETRA* mgtet1 = (MG_TETRA*) no1->get_lien_tetra()->get(k);
245     // // int origine1 = mgtet1->get_origine();
246     // // if ((origine1 == OPTIMISE) || (origine1 == IMPOSE)) liste_voisins.ajouter(mgtet1);
247     // // }
248     // // MG_NOEUD* no2 = mgtet_1->get_noeud2();
249     // // int nb_vois2 = no2->get_lien_tetra()->get_nb();
250     // // for (int k = 0;k<nb_vois2;k++)
251     // // {
252     // // MG_TETRA* mgtet2 = (MG_TETRA*) no2->get_lien_tetra()->get(k);
253     // // int origine2 = mgtet2->get_origine();
254     // // if ((origine2 == OPTIMISE) || (origine2 == IMPOSE)) liste_voisins.ajouter(mgtet2);
255     // // }
256     // // MG_NOEUD* no3 = mgtet_1->get_noeud3();
257     // // int nb_vois3 = no3->get_lien_tetra()->get_nb();
258     // // for (int k = 0;k<nb_vois3;k++)
259     // // {
260     // // MG_TETRA* mgtet3 = (MG_TETRA*) no3->get_lien_tetra()->get(k);
261     // // int origine3 = mgtet3->get_origine();
262     // // if ((origine3 == OPTIMISE) || (origine3 == IMPOSE)) liste_voisins.ajouter(mgtet3);
263     // // }
264     // // MG_NOEUD* no4 = mgtet_1->get_noeud4();
265     // // int nb_vois4 = no4->get_lien_tetra()->get_nb();
266     // // for (int k = 0;k<nb_vois4;k++)
267     // // {
268     // // MG_TETRA* mgtet4 = (MG_TETRA*) no4->get_lien_tetra()->get(k);
269     // // int origine4 = mgtet4->get_origine();
270     // // if ((origine4 == OPTIMISE) || (origine4 == IMPOSE)) liste_voisins.ajouter(mgtet4);
271     // // }
272     // }
273     // }
274     int nb_voisins = liste_voisins.get_nb();
275     TPL_MAP_ENTITE<MG_TETRA*>::ITERATEUR it;
276     for(MG_TETRA* mgtet_j=liste_voisins.get_premier(it);mgtet_j!=NULL;mgtet_j=liste_voisins.get_suivant(it))
277     {
278     MG_NOEUD* noeud1 = mgtet_j->get_noeud1();
279     double xx1 = noeud1->get_x();
280     double yy1 = noeud1->get_y();
281     double zz1 = noeud1->get_z();
282     int ori_no1 = noeud1->get_origine();
283     MG_NOEUD* noeud11 = new MG_NOEUD(NULL,xx1,yy1,zz1,ori_no1);
284     MG_NOEUD* noeud2 = mgtet_j->get_noeud2();
285     double xx2 = noeud2->get_x();
286     double yy2 = noeud2->get_y();
287     double zz2 = noeud2->get_z();
288     int ori_no2 = noeud2->get_origine();
289     MG_NOEUD* noeud22 = new MG_NOEUD(NULL,xx2,yy2,zz2,ori_no2);
290     MG_NOEUD* noeud3 = mgtet_j->get_noeud3();
291     double xx3 = noeud3->get_x();
292     double yy3 = noeud3->get_y();
293     double zz3 = noeud3->get_z();
294     int ori_no3 = noeud3->get_origine();
295     MG_NOEUD* noeud33 = new MG_NOEUD(NULL,xx3,yy3,zz3,ori_no3);
296     MG_NOEUD* noeud4 = mgtet_j->get_noeud4();
297     double xx4 = noeud4->get_x();
298     double yy4 = noeud4->get_y();
299     double zz4 = noeud4->get_z();
300     int ori_no4 = noeud4->get_origine();
301     MG_NOEUD* noeud44 = new MG_NOEUD(NULL,xx4,yy4,zz4,ori_no4);
302     mai3->ajouter_mg_noeud(noeud11);
303     mai3->ajouter_mg_noeud(noeud22);
304     mai3->ajouter_mg_noeud(noeud33);
305     mai3->ajouter_mg_noeud(noeud44);
306     MG_TRIANGLE* tri1 = mgtet_j->get_triangle1();
307     MG_NOEUD* no11 = tri1->get_noeud1();
308     MG_NOEUD* no12 = tri1->get_noeud2();
309     MG_NOEUD* no13 = tri1->get_noeud3();
310     MG_SEGMENT* seg11 = tri1->get_segment1();
311     MG_SEGMENT* seg12 = tri1->get_segment2();
312     MG_SEGMENT* seg13 = tri1->get_segment3();
313     int ori_tri1 = tri1->get_origine();
314     MG_TRIANGLE* tri11 = new MG_TRIANGLE(NULL,no11,no12,no13,seg11,seg12,seg13,ori_tri1);
315     MG_TRIANGLE* tri2 = mgtet_j->get_triangle2();
316     MG_NOEUD* no21 = tri2->get_noeud1();
317     MG_NOEUD* no22 = tri2->get_noeud2();
318     MG_NOEUD* no23 = tri2->get_noeud3();
319     MG_SEGMENT* seg21 = tri2->get_segment1();
320     MG_SEGMENT* seg22 = tri2->get_segment2();
321     MG_SEGMENT* seg23 = tri2->get_segment3();
322     int ori_tri2 = tri2->get_origine();
323     MG_TRIANGLE* tri22 = new MG_TRIANGLE(NULL,no21,no22,no23,seg21,seg22,seg23,ori_tri2);
324     MG_TRIANGLE* tri3 = mgtet_j->get_triangle3();
325     MG_NOEUD* no31 = tri3->get_noeud1();
326     MG_NOEUD* no32 = tri3->get_noeud2();
327     MG_NOEUD* no33 = tri3->get_noeud3();
328     MG_SEGMENT* seg31 = tri3->get_segment1();
329     MG_SEGMENT* seg32 = tri3->get_segment2();
330     MG_SEGMENT* seg33 = tri3->get_segment3();
331     int ori_tri3 = tri3->get_origine();
332     MG_TRIANGLE* tri33 = new MG_TRIANGLE(NULL,no31,no32,no33,seg31,seg32,seg33,ori_tri3);
333     MG_TRIANGLE* tri4 = mgtet_j->get_triangle4();
334     MG_NOEUD* no41 = tri4->get_noeud1();
335     MG_NOEUD* no42 = tri4->get_noeud2();
336     MG_NOEUD* no43 = tri4->get_noeud3();
337     MG_SEGMENT* seg41 = tri4->get_segment1();
338     MG_SEGMENT* seg42 = tri4->get_segment2();
339     MG_SEGMENT* seg43 = tri4->get_segment3();
340     int ori_tri4 = tri4->get_origine();
341     MG_TRIANGLE* tri44 = new MG_TRIANGLE(NULL,no41,no42,no43,seg41,seg42,seg43,ori_tri4);
342     mai3->ajouter_mg_triangle(tri11);
343     mai3->ajouter_mg_triangle(tri22);
344     mai3->ajouter_mg_triangle(tri33);
345     mai3->ajouter_mg_triangle(tri44);
346     int ori_tet = mgtet_j->get_origine();
347     MG_TETRA* mgtet = new MG_TETRA(NULL,noeud11,noeud22,noeud33,noeud44,tri11,tri22,tri33,tri44,ori_tet);
348     mai3->ajouter_mg_tetra(mgtet);
349     }
350     gest3.enregistrer(nomfichier);
351     }
352    
353 picher 248 void MG_LISSAGE::lisse(FEM_MAILLAGE* mai, MG_GESTIONNAIRE& gest2, double epsilon, double sigma, double sigmaf, double sigmag, double gamma, int iter_max, int reactiv, int bruit, int liss,int opti,int imp,int m_auto)
354 picher 233 {
355 francois 234 //Menu de la m�thode de lissage
356 francois 232 affiche("Extraction du maillage de surface");
357 picher 231 int coderesu = 0;
358 picher 233 double frac_min = 0.005;
359     int mai2_id;
360 francois 232 if (opti==1) conserve(OPTIMISE);
361     if (imp==1) conserve(IMPOSE);
362 picher 233 if (m_auto==1) conserve(MAILLEUR_AUTO);
363     copieorigine(mai,gest2);
364     MG_MAILLAGE* mg_mai = (MG_MAILLAGE*)mai->get_mg_maillage();
365     gain_poids(mg_mai,gest2);
366     if (reactiv == 1)
367     {
368     reactivation(mg_mai,gest2);
369     }
370 picher 231 do
371     {
372 picher 233 coderesu = extract_skin(mg_mai,gest2,frac_min,&mai2_id);
373 francois 234 gain_poids(mg_mai,gest2);
374 picher 231 }
375     while (coderesu == 0);
376 picher 233
377     MG_MAILLAGE* mg_mai2=gest2.get_mg_maillageid(mai2_id);
378 francois 232 if (bruit == 1)
379 picher 231 {
380 francois 232 affiche("Bruitage du maillage initial");
381 picher 233 bruitage(mg_mai2,gest2);
382 picher 231 }
383 picher 248 //Choix de la méthode de lissage (1 = chen(2005), 2 = Jones(2007), 3 = chen(2008))
384 francois 232 if (liss == 1)
385 picher 231 {
386 francois 234 affiche("Procedure de lissage");
387 picher 233 lissage(mg_mai2,gest2,epsilon,sigma,iter_max);
388 picher 231 }
389 picher 248 if (liss == 2)
390     {
391     affiche("Procedure de lissage");
392     lissage2(mg_mai2,gest2,sigmaf,sigmag,iter_max);
393     }
394     if (liss == 3)
395     {
396     affiche("Procedure de lissage");
397     lissage3(mg_mai2,gest2,sigma,gamma,epsilon,iter_max);
398     }
399 picher 231 }
400    
401 picher 230 void MG_LISSAGE::bruitage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
402     {
403 francois 234 //M�thode test pour bruiter un maillage initial et tester la m�thode de lissage
404 picher 230 srand(time(NULL));
405     LISTE_MG_NOEUD::iterator it_no;
406     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
407     {
408 francois 234 //On modifie la position des noeuds l�g�rement (ajout d'un bruit au maillage initial)
409 picher 230 int signe;
410     int delta;
411     double delta_x,delta_y,delta_z;
412     signe = rand()%2 + 1;
413     delta = rand()%2 + 1;
414     if(signe == 1)
415     {
416 picher 248 delta_x = delta*0.5;
417 picher 230 }
418     else //(signe == 2)
419     {
420 picher 248 delta_x = -delta*0.5;
421 picher 230 }
422     signe = rand()%2 + 1;
423     delta = rand()%2 + 1;
424     if(signe == 1)
425     {
426 picher 248 delta_y = delta*0.5;
427 picher 230 }
428     else //(signe == 2)
429     {
430 picher 248 delta_y = -delta*0.5;
431 picher 230 }
432     signe = rand()%2 + 1;
433     delta = rand()%2 + 1;
434     if(signe == 1)
435     {
436 picher 248 delta_z = delta*0.5;
437 picher 230 }
438     else //(signe == 2)
439     {
440 picher 248 delta_z = -delta*0.5;
441 picher 230 }
442     double x_new = noeud_i->get_x() + delta_x;
443     double y_new = noeud_i->get_y() + delta_y;
444     double z_new = noeud_i->get_z() + delta_z;
445     noeud_i->change_x(x_new);
446     noeud_i->change_y(y_new);
447     noeud_i->change_z(z_new);
448     }
449     }
450 francois 224
451 picher 230 void MG_LISSAGE::lissage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max)
452     {
453 francois 234 //M�thode de lissage inspir�e de Chen(2005)
454 picher 230 double un_sur_pi = 1./M_PI;
455     int compteur = 0;
456     int fin = 0;
457 picher 231 FILE *out = fopen("convergence.txt","wt");
458 francois 224
459 picher 230 do
460     {
461     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
462 picher 233 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
463 picher 230 TPL_LISTE_ENTITE<double> liste_wij;
464     LISTE_MG_TRIANGLE::iterator it_tri;
465     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
466     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
467     {
468     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
469     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
470 picher 233 normal_f_i.norme();
471     liste_normales2.ajouter(normal_f_i);
472 picher 230 //Remplissage de la liste des voisins du triangle i
473     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
474     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
475     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
476     for (int j = 0;j<nb_voisins1;j++)
477     {
478     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
479     liste_voisins.ajouter(mgtri_1);
480     }
481     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
482     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
483     for (int j = 0;j<nb_voisins2;j++)
484     {
485     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
486     liste_voisins.ajouter(mgtri_2);
487     }
488     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
489     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
490     for (int j = 0;j<nb_voisins3;j++)
491     {
492     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
493     liste_voisins.ajouter(mgtri_3);
494     }
495     liste_voisins.supprimer(mgtri_i);
496     int nb_voisins = liste_voisins.get_nb();
497     double w_ij = 1./nb_voisins;
498     double phi_i_min = 10.;
499     double s_i = 0.0;
500     double phi_im = 0.0;
501 picher 231 double *phi_ij = new double[nb_voisins];
502 picher 230 OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
503     OT_VECTEUR_3D eta_i(0.,0.,0.);
504     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
505     int j = 0;
506     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
507     {
508     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
509     //1-Calculer la normale moyenne pour chaque triangle
510     normal_f_i_mean = normal_f_i_mean + w_ij*normal_f_j;
511 francois 234 //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 � Nb_voisins
512 picher 230 double prod_scalaire = normal_f_i*normal_f_j;
513 picher 233 if (prod_scalaire > 1.)
514     {
515     prod_scalaire = 1.;
516     }
517     if (prod_scalaire < -1.)
518     {
519     prod_scalaire = -1.;
520     }
521 picher 230 phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
522     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
523     if (phi_ij[j] < phi_i_min)
524     {
525     phi_i_min = phi_ij[j];
526     eta_i = normal_f_j;
527     }
528     //3.1-On calcul l'angle moyen phi_im
529     phi_im = phi_im + w_ij*phi_ij[j];
530     j++;
531     }
532     normal_f_i_mean.norme();
533     j = 0;
534     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
535     {
536     //3.2-Calcul de s_i selon la variance
537     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
538     j++;
539     }
540 picher 231 delete[] phi_ij;
541 picher 230 //4-On calcule une nouvelle normale pour chaque triangle
542     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
543     normal_f_i_new.norme();
544     liste_normales.ajouter(normal_f_i_new);
545     liste_wij.ajouter(w_ij);
546     mgtri->change_nouveau_numero(k);
547     k++;
548     }
549    
550     LISTE_MG_NOEUD::iterator it_no;
551     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
552     {
553     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
554     double w_ij_prime = 0.0;
555     OT_VECTEUR_3D v_temp(0.,0.,0.);
556     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
557     for(int j=0;j<nb_voisins_j;j++)
558     {
559     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
560 francois 234 //On calcule le centro�de cj du triangle mgtri_j
561 picher 230 MG_NOEUD* n1 = mgtri_j->get_noeud1();
562     MG_NOEUD* n2 = mgtri_j->get_noeud2();
563     MG_NOEUD* n3 = mgtri_j->get_noeud3();
564     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
565     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
566     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
567     //On forme le vecteur vi_cj
568     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
569     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
570 francois 234 // w_ij_prime correspond � la somme des aires des triangles voisins du noeuds
571 picher 233 OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
572     OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
573     OT_VECTEUR_3D prodvect = AB&AC;
574     double w_ij = 0.5*prodvect.get_longueur();
575 picher 230 w_ij_prime = w_ij_prime + w_ij;
576     v_temp = v_temp + w_ij*(vi_cj*normal_f_i_new)*normal_f_i_new;
577     }
578 francois 234 //5-On met � jour la position des noeuds
579 picher 230 v_i = v_i + v_temp/w_ij_prime;
580 picher 233 int origine = noeud_i->get_origine();
581     if (origine != IMPOSE)
582     {
583     noeud_i->change_x(v_i.get_x());
584     noeud_i->change_y(v_i.get_y());
585     noeud_i->change_z(v_i.get_z());
586     }
587 picher 230 }
588    
589 francois 234 //Crit�re d'arr�t de l'algorithme
590 picher 230 int l=0;
591     int nb_tri = mg_mai->get_nb_mg_triangle();
592     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
593     {
594     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
595 picher 233 OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
596 picher 230 OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
597 picher 248 double critere = 1. - normal_f_i*normal_f_i_new;
598 picher 231 fprintf(out,"%lf\n",critere);
599 picher 230 if (critere <= epsilon) l++;
600 picher 233 if (critere >= 1.)
601     {
602     int no_id = mgtri->get_noeud1()->get_id();
603     char str[50];
604     sprintf (str, "test%d.magic", no_id);
605     }
606 picher 230 }
607 picher 248 double tolerance = 0.01*nb_tri;
608     if (nb_tri - l <= 0) fin = 1;
609 picher 230 compteur++;
610 picher 248 cout << "Fin de l'iteration #" << compteur << endl;
611     fprintf(out,"\nFin de l'iteration #%d\n\n",compteur);
612 picher 230 }
613     while ((fin == 0) && (compteur < iter_max));
614    
615     if (fin == 1)
616     {
617 picher 248 cout << "Convergence de la methode en " << compteur << " iterations" << endl;
618 picher 230 }
619     else
620     {
621 picher 248 cout << "Arret de la procedure apres " << compteur << " iterations" << endl;
622 picher 230 }
623    
624 picher 231 fclose(out);
625 picher 230 } //Fin de lissage
626    
627 picher 248 void MG_LISSAGE::lissage2(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double sigmaf, double sigmag, int iter_max)
628     {
629     //Non-Iterative Smoothing (NIS) (Jones(2007))
630     LISTE_MG_NOEUD::iterator it_no;
631     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_noeuds;
632     int l = 0;
633     //1ere passe de NIS
634     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
635     {
636     double somme_num_x = 0.;
637     double somme_num_y = 0.;
638     double somme_num_z = 0.;
639     double somme_denom = 0.;
640     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
641     int nb_voisins1 = noeud_i->get_lien_triangle()->get_nb();
642     for (int j = 0;j<nb_voisins1;j++)
643     {
644     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
645     liste_voisins.ajouter(mgtri_1);
646     MG_NOEUD* no1 = mgtri_1->get_noeud1();
647     int nb_vois1 = no1->get_lien_triangle()->get_nb();
648     for (int k = 0;k<nb_vois1;k++)
649     {
650     MG_TRIANGLE_PEAU* mgtri1 = (MG_TRIANGLE_PEAU*) no1->get_lien_triangle()->get(k);
651     liste_voisins.ajouter(mgtri1);
652     }
653     MG_NOEUD* no2 = mgtri_1->get_noeud2();
654     int nb_vois2 = no2->get_lien_triangle()->get_nb();
655     for (int k = 0;k<nb_vois2;k++)
656     {
657     MG_TRIANGLE_PEAU* mgtri2 = (MG_TRIANGLE_PEAU*) no2->get_lien_triangle()->get(k);
658     liste_voisins.ajouter(mgtri2);
659     }
660     MG_NOEUD* no3 = mgtri_1->get_noeud3();
661     int nb_vois3 = no3->get_lien_triangle()->get_nb();
662     for (int k = 0;k<nb_vois3;k++)
663     {
664     MG_TRIANGLE_PEAU* mgtri3 = (MG_TRIANGLE_PEAU*) no3->get_lien_triangle()->get(k);
665     liste_voisins.ajouter(mgtri3);
666     }
667     }
668     int nb_voisins = liste_voisins.get_nb();
669     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
670     for(MG_TRIANGLE_PEAU* mgtri_peau=liste_voisins.get_premier(it);mgtri_peau!=NULL;mgtri_peau=liste_voisins.get_suivant(it))
671     {
672     OT_VECTEUR_3D normale_j = mgtri_peau->calcul_normal();
673     //Centroide du triangle
674     MG_NOEUD* n1 = mgtri_peau->get_noeud1();
675     MG_NOEUD* n2 = mgtri_peau->get_noeud2();
676     MG_NOEUD* n3 = mgtri_peau->get_noeud3();
677     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
678     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
679     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
680     //spatial weight
681     double px = noeud_i->get_x();
682     double py = noeud_i->get_y();
683     double pz = noeud_i->get_z();
684     double spatial = sqrt(pow(cj_x - px,2) + pow(cj_y - py,2) +pow(cj_z - pz,2));
685     double f = ponderation_gaussian(spatial,sigmaf/2.);
686     //calcul de l'aire du triangle
687     OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
688     OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
689     OT_VECTEUR_3D prodvect = AB&AC;
690     double aire_q = 0.5*prodvect.get_longueur();
691     somme_num_x = somme_num_x + cj_x*aire_q*f;
692     somme_num_y = somme_num_y + cj_y*aire_q*f;
693     somme_num_z = somme_num_z + cj_z*aire_q*f;
694     somme_denom = somme_denom + aire_q*f;
695     }
696     //On ne change pas réellement la position des noeuds lors de la première passe, mais on conserve les normales mollifiées !!
697     OT_VECTEUR_3D p_prime(somme_num_x/somme_denom,somme_num_y/somme_denom,somme_num_z/somme_denom);
698     liste_noeuds.ajouter(p_prime);
699     noeud_i->change_nouveau_numero(l);
700     l++;
701     }
702     //Calcul des nouvelles normales
703     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
704     LISTE_MG_TRIANGLE::iterator it_tri;
705     l = 0;
706     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
707     {
708     MG_NOEUD* n1 = mgtri->get_noeud1();
709     OT_VECTEUR_3D p1 = liste_noeuds.get(n1->get_nouveau_numero());
710     MG_NOEUD* n2 = mgtri->get_noeud2();
711     OT_VECTEUR_3D p2 = liste_noeuds.get(n2->get_nouveau_numero());
712     MG_NOEUD* n3 = mgtri->get_noeud3();
713     OT_VECTEUR_3D p3 = liste_noeuds.get(n3->get_nouveau_numero());
714     OT_VECTEUR_3D AB(p2.get_x() - p1.get_x(),p2.get_y() - p1.get_y(),p2.get_z() - p1.get_z());
715     OT_VECTEUR_3D AC(p3.get_x() - p1.get_x(),p3.get_y() - p1.get_y(),p3.get_z() - p1.get_z());
716     OT_VECTEUR_3D normale = AB&AC;
717     normale.norme();
718     liste_normales.ajouter(normale);
719     mgtri->change_nouveau_numero(l);
720     l++;
721     }
722     //2e passe de NIS
723     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_pprime;
724     l = 0;
725     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
726     {
727     double somme_num_x = 0.;
728     double somme_num_y = 0.;
729     double somme_num_z = 0.;
730     double somme_denom = 0.;
731     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
732     int nb_voisins1 = noeud_i->get_lien_triangle()->get_nb();
733     for (int j = 0;j<nb_voisins1;j++)
734     {
735     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
736     liste_voisins.ajouter(mgtri_1);
737     MG_NOEUD* no1 = mgtri_1->get_noeud1();
738     int nb_vois1 = no1->get_lien_triangle()->get_nb();
739     for (int k = 0;k<nb_vois1;k++)
740     {
741     MG_TRIANGLE_PEAU* mgtri1 = (MG_TRIANGLE_PEAU*) no1->get_lien_triangle()->get(k);
742     liste_voisins.ajouter(mgtri1);
743     }
744     MG_NOEUD* no2 = mgtri_1->get_noeud2();
745     int nb_vois2 = no2->get_lien_triangle()->get_nb();
746     for (int k = 0;k<nb_vois2;k++)
747     {
748     MG_TRIANGLE_PEAU* mgtri2 = (MG_TRIANGLE_PEAU*) no2->get_lien_triangle()->get(k);
749     liste_voisins.ajouter(mgtri2);
750     }
751     MG_NOEUD* no3 = mgtri_1->get_noeud3();
752     int nb_vois3 = no3->get_lien_triangle()->get_nb();
753     for (int k = 0;k<nb_vois3;k++)
754     {
755     MG_TRIANGLE_PEAU* mgtri3 = (MG_TRIANGLE_PEAU*) no3->get_lien_triangle()->get(k);
756     liste_voisins.ajouter(mgtri3);
757     }
758     }
759     int nb_voisins = liste_voisins.get_nb();
760     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
761     for(MG_TRIANGLE_PEAU* mgtri_peau=liste_voisins.get_premier(it);mgtri_peau!=NULL;mgtri_peau=liste_voisins.get_suivant(it))
762     {
763     OT_VECTEUR_3D normale_j = liste_normales.get(mgtri_peau->get_nouveau_numero());
764     //Centroide du triangle
765     MG_NOEUD* n1 = mgtri_peau->get_noeud1();
766     MG_NOEUD* n2 = mgtri_peau->get_noeud2();
767     MG_NOEUD* n3 = mgtri_peau->get_noeud3();
768     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
769     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
770     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
771     //spatial weight
772     double px = noeud_i->get_x();
773     double py = noeud_i->get_y();
774     double pz = noeud_i->get_z();
775     double spatial = sqrt(pow(cj_x - px,2) + pow(cj_y - py,2) +pow(cj_z - pz,2));
776     double f = ponderation_gaussian(spatial,sigmaf);
777     //calcul de l'aire du triangle
778     OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
779     OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
780     OT_VECTEUR_3D prodvect = AB&AC;
781     double aire_q = 0.5*prodvect.get_longueur();
782     //Calcul de la projection du noeuds sur le plan tangent au triangle
783     double a = normale_j.get_x();
784     double b = normale_j.get_y();
785     double c = normale_j.get_z();
786     double d = -(a*n1->get_x() + b*n1->get_y() + c*n1->get_z());
787     double k = -(a*px + b*py + c*pz + d)/(a*a + b*b + c*c);
788     double proj_x = px + k*a;
789     double proj_y = py + k*b;
790     double proj_z = pz + k*c;
791     //influence weight
792     double influence = sqrt(pow(proj_x - px,2) + pow(proj_y - py,2) + pow(proj_z - pz,2));
793     double g = ponderation_gaussian(influence,sigmag);
794     somme_num_x = somme_num_x + proj_x*aire_q*f*g;
795     somme_num_y = somme_num_y + proj_y*aire_q*f*g;
796     somme_num_z = somme_num_z + proj_z*aire_q*f*g;
797     somme_denom = somme_denom + aire_q*f*g;
798     }
799     OT_VECTEUR_3D pprime(somme_num_x/somme_denom,somme_num_y/somme_denom,somme_num_z/somme_denom);
800     liste_pprime.ajouter(pprime);
801     noeud_i->change_nouveau_numero(l);
802     l++;
803     }
804     //Mise à jour en bloc de la position des noeuds
805     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
806     {
807     int origine = noeud_i->get_origine();
808     if (origine != IMPOSE)
809     {
810     OT_VECTEUR_3D pprime = liste_pprime.get(noeud_i->get_nouveau_numero());
811     noeud_i->change_x(pprime.get_x());
812     noeud_i->change_y(pprime.get_y());
813     noeud_i->change_z(pprime.get_z());
814     }
815     }
816    
817     } //Fin de lissage2
818    
819     void MG_LISSAGE::lissage3(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double sigma, double gamma, double epsilon, int iter_max)
820     {
821     //M�thode de lissage inspir�e de Chen(2008)
822     double un_sur_pi = 1./M_PI;
823     int compteur = 0;
824     int fin = 0;
825     FILE *out = fopen("convergence.txt","wt");
826    
827     do
828     {
829     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
830     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
831     TPL_LISTE_ENTITE<double> liste_wij;
832     LISTE_MG_TRIANGLE::iterator it_tri;
833     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
834     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
835     {
836     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
837     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
838     normal_f_i.norme();
839     liste_normales2.ajouter(normal_f_i);
840     //Remplissage de la liste des voisins du triangle i
841     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
842     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
843     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
844     for (int j = 0;j<nb_voisins1;j++)
845     {
846     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
847     liste_voisins.ajouter(mgtri_1);
848     }
849     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
850     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
851     for (int j = 0;j<nb_voisins2;j++)
852     {
853     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
854     liste_voisins.ajouter(mgtri_2);
855     }
856     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
857     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
858     for (int j = 0;j<nb_voisins3;j++)
859     {
860     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
861     liste_voisins.ajouter(mgtri_3);
862     }
863     liste_voisins.supprimer(mgtri_i);
864     int nb_voisins = liste_voisins.get_nb();
865     double w_ij = 1./nb_voisins;
866     double phi_i_min = 10.;
867     double s_i = 0.0;
868     double phi_im = 0.0;
869     double *phi_ij = new double[nb_voisins];
870     OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
871     normal_f_i_mean = normal_f_i;
872     OT_VECTEUR_3D eta_i(0.,0.,0.);
873     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
874     int j = 0;
875     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
876     {
877     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
878     //1-Calculer la normale moyenne pour chaque triangle
879     normal_f_i_mean = normal_f_i_mean + normal_f_j;
880     //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 � Nb_voisins
881     double prod_scalaire = normal_f_i*normal_f_j;
882     if (prod_scalaire > 1.)
883     {
884     prod_scalaire = 1.;
885     }
886     if (prod_scalaire < -1.)
887     {
888     prod_scalaire = -1.;
889     }
890     phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
891     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
892     if (phi_ij[j] < phi_i_min)
893     {
894     phi_i_min = phi_ij[j];
895     eta_i = normal_f_j;
896     }
897     //3.1-On calcul l'angle moyen phi_im
898     phi_im = phi_im + w_ij*phi_ij[j];
899     j++;
900     }
901     normal_f_i_mean.norme();
902     j = 0;
903     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
904     {
905     //3.2-Calcul de s_i selon la variance
906     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
907     j++;
908     }
909     delete[] phi_ij;
910     //4-On calcule une nouvelle normale pour chaque triangle
911     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
912     normal_f_i_new.norme();
913     liste_normales.ajouter(normal_f_i_new);
914     liste_wij.ajouter(w_ij);
915     mgtri->change_nouveau_numero(k);
916     k++;
917     }
918    
919     LISTE_MG_NOEUD::iterator it_no;
920     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
921     {
922     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
923     double w_ij_prime = 0.0;
924     OT_VECTEUR_3D v_temp(0.,0.,0.);
925     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
926     for(int j=0;j<nb_voisins_j;j++)
927     {
928     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
929     //On calcule le centro�de cj du triangle mgtri_j
930     MG_NOEUD* n1 = mgtri_j->get_noeud1();
931     MG_NOEUD* n2 = mgtri_j->get_noeud2();
932     MG_NOEUD* n3 = mgtri_j->get_noeud3();
933     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
934     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
935     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
936     //On forme le vecteur vi_cj
937     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
938     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
939     v_temp = v_temp + (vi_cj*normal_f_i_new)*normal_f_i_new;
940     }
941     //5-On met � jour la position des noeuds
942     v_i = v_i + (gamma/(2*nb_voisins_j))*v_temp;
943     int origine = noeud_i->get_origine();
944     if (origine != IMPOSE)
945     {
946     noeud_i->change_x(v_i.get_x());
947     noeud_i->change_y(v_i.get_y());
948     noeud_i->change_z(v_i.get_z());
949     }
950     }
951    
952     //Crit�re d'arr�t de l'algorithme
953     int l=0;
954     int nb_tri = mg_mai->get_nb_mg_triangle();
955     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
956     {
957     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
958     OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
959     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
960     double critere = 1. - normal_f_i*normal_f_i_new;
961     fprintf(out,"%lf\n",critere);
962     if (critere <= epsilon) l++;
963     if (critere >= 1.)
964     {
965     int no_id = mgtri->get_noeud1()->get_id();
966     char str[50];
967     sprintf (str, "test%d.magic", no_id);
968     }
969     }
970     double tolerance = 0.01*nb_tri;
971     if (nb_tri - l <= 0) fin = 1;
972     compteur++;
973     cout << "Fin de l'iteration #" << compteur << endl;
974     fprintf(out,"\nFin de l'iteration #%d\n\n",compteur);
975     }
976     while ((fin == 0) && (compteur < iter_max));
977    
978     if (fin == 1)
979     {
980     cout << "Convergence de la methode en " << compteur << " iterations" << endl;
981     }
982     else
983     {
984     cout << "Arret de la procedure apres " << compteur << " iterations" << endl;
985     }
986    
987     fclose(out);
988     }
989    
990 picher 233 int MG_LISSAGE::extract_skin(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2,double frac_min, int *mai2_id)
991 francois 222 {
992    
993     //�tape 0 - On commence par mettre � z�ro tous les nouveau_numero des triangles et des noeuds du maillage
994     LISTE_MG_TRIANGLE::iterator it_tri;
995     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
996     {
997     mgtri->change_nouveau_numero(0);
998 francois 224 mgtri->change_origine(MAILLEUR_AUTO);
999 francois 222 }
1000     LISTE_MG_NOEUD::iterator it_noeud;
1001     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
1002     {
1003     mgnoeud->change_nouveau_numero(0);
1004     }
1005    
1006     //�tape 1 - On boucle ensuite tous les t�tra�dres pour ajouter un compteur du nombre de fois qu'un triangle appartient � 1 t�tra
1007 picher 233 LISTE_MG_TETRA::iterator it_tetra;
1008     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
1009 francois 222 {
1010     int origine = mgtet->get_origine();
1011 francois 224 if (origine==IMPOSE)
1012 francois 222 {
1013 francois 224 mgtet->get_triangle1()->change_origine(IMPOSE);
1014     mgtet->get_triangle2()->change_origine(IMPOSE);
1015     mgtet->get_triangle3()->change_origine(IMPOSE);
1016     mgtet->get_triangle4()->change_origine(IMPOSE);
1017     }
1018     if (origine==OPTIMISE)
1019     {
1020     if (mgtet->get_triangle1()->get_origine()!=IMPOSE) mgtet->get_triangle1()->change_origine(OPTIMISE);
1021     if (mgtet->get_triangle2()->get_origine()!=IMPOSE) mgtet->get_triangle2()->change_origine(OPTIMISE);
1022     if (mgtet->get_triangle3()->get_origine()!=IMPOSE) mgtet->get_triangle3()->change_origine(OPTIMISE);
1023     if (mgtet->get_triangle4()->get_origine()!=IMPOSE) mgtet->get_triangle4()->change_origine(OPTIMISE);
1024     }
1025    
1026     if (((origine == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
1027     ((origine == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
1028     ((origine == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
1029     ((origine == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
1030     ((origine == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
1031     ((origine == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
1032    
1033     {
1034 francois 222 int num1 = mgtet->get_triangle1()->get_nouveau_numero();
1035     int num2 = mgtet->get_triangle2()->get_nouveau_numero();
1036     int num3 = mgtet->get_triangle3()->get_nouveau_numero();
1037     int num4 = mgtet->get_triangle4()->get_nouveau_numero();
1038     num1++;
1039     num2++;
1040     num3++;
1041     num4++;
1042     mgtet->get_triangle1()->change_nouveau_numero(num1);
1043     mgtet->get_triangle2()->change_nouveau_numero(num2);
1044     mgtet->get_triangle3()->change_nouveau_numero(num3);
1045     mgtet->get_triangle4()->change_nouveau_numero(num4);
1046     }
1047     }
1048    
1049     //�tape 2 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour identifier les noeuds leur appartenant
1050     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1051     {
1052     int num = mgtri->get_nouveau_numero();
1053     if (num == 1)
1054     {
1055     MG_NOEUD* noeud1 = mgtri->get_noeud1();
1056     MG_NOEUD* noeud2 = mgtri->get_noeud2();
1057     MG_NOEUD* noeud3 = mgtri->get_noeud3();
1058     noeud1->change_nouveau_numero(1);
1059     noeud2->change_nouveau_numero(1);
1060     noeud3->change_nouveau_numero(1);
1061 picher 248 if (mgtri->get_origine()==IMPOSE)
1062     {
1063     noeud1->change_origine(IMPOSE);
1064     noeud2->change_origine(IMPOSE);
1065     noeud3->change_origine(IMPOSE);
1066     }
1067 francois 222 }
1068     }
1069    
1070 francois 234 //�tape 3 - On cr�e un nouveau maillage pour la peau
1071 francois 222 MG_MAILLAGE* mai2 = new MG_MAILLAGE(NULL);
1072     gest2.ajouter_mg_maillage(mai2);
1073    
1074     // //�tape 4 - On boucle l'ensemble des noeuds identifi�s � l'�tape 2 pour les recr�er dans le second maillage
1075     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
1076     {
1077     int num = mgnoeud->get_nouveau_numero();
1078     if (num == 1)
1079     {
1080     double x = mgnoeud->get_x();
1081     double y = mgnoeud->get_y();
1082     double z = mgnoeud->get_z();
1083 francois 224 int origine=TRIANGULATION;
1084     if (mgnoeud->get_origine()==IMPOSE) origine=IMPOSE;
1085     MG_NOEUD* noeud1 = new MG_NOEUD(NULL,x,y,z,origine);
1086 francois 222 mai2->ajouter_mg_noeud(noeud1);
1087     mgnoeud->change_nouveau_numero(noeud1->get_id());
1088 francois 234 noeud1->change_nouveau_numero(mgnoeud->get_id());
1089 francois 222 }
1090     }
1091    
1092     // //�tape 5 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour les recr�er dans le maillage 2
1093 picher 233 for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
1094 francois 222 {
1095 francois 224 int originetet=mgtet->get_origine();
1096     if (((originetet == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
1097     ((originetet == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
1098     ((originetet == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
1099     ((originetet == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
1100     ((originetet == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
1101     ((originetet == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
1102 francois 222 {
1103 francois 224 MG_NOEUD* noeud1 = mgtet->get_noeud1();
1104     MG_NOEUD* noeud2 = mgtet->get_noeud2();
1105     MG_NOEUD* noeud3 = mgtet->get_noeud3();
1106     MG_NOEUD* noeud4 = mgtet->get_noeud4();
1107 francois 222 MG_NOEUD* node1 = mai2->get_mg_noeudid(noeud1->get_nouveau_numero());
1108     MG_NOEUD* node2 = mai2->get_mg_noeudid(noeud2->get_nouveau_numero());
1109     MG_NOEUD* node3 = mai2->get_mg_noeudid(noeud3->get_nouveau_numero());
1110 francois 224 MG_NOEUD* node4 = mai2->get_mg_noeudid(noeud4->get_nouveau_numero());
1111     MG_TRIANGLE* tri1=mgtet->get_triangle1();
1112     MG_TRIANGLE* tri2=mgtet->get_triangle2();
1113     MG_TRIANGLE* tri3=mgtet->get_triangle3();
1114     MG_TRIANGLE* tri4=mgtet->get_triangle4();
1115     if (tri1->get_nouveau_numero()==1)
1116     {
1117     int origine=TRIANGULATION;
1118     if (tri1->get_origine()==IMPOSE) origine=IMPOSE;
1119     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node2,node3,mai2,origine);
1120     tripeau->change_nouveau_numero(0);
1121     tri1->change_nouveau_numero(0);
1122     }
1123     if (tri2->get_nouveau_numero()==1)
1124     {
1125     int origine=TRIANGULATION;
1126     if (tri2->get_origine()==IMPOSE) origine=IMPOSE;
1127     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node4,node2,mai2,origine);
1128     tripeau->change_nouveau_numero(0);
1129     tri2->change_nouveau_numero(0);
1130     }
1131     if (tri3->get_nouveau_numero()==1)
1132     {
1133     int origine=TRIANGULATION;
1134     if (tri3->get_origine()==IMPOSE) origine=IMPOSE;
1135     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node2,node4,node3,mai2,origine);
1136     tripeau->change_nouveau_numero(0);
1137     tri3->change_nouveau_numero(0);
1138     }
1139     if (tri4->get_nouveau_numero()==1)
1140     {
1141     int origine=TRIANGULATION;
1142     if (tri4->get_origine()==IMPOSE) origine=IMPOSE;
1143     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node3,node4,mai2,origine);
1144     tripeau->change_nouveau_numero(0);
1145     tri4->change_nouveau_numero(0);
1146     }
1147 francois 222 }
1148 francois 224 }
1149     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
1150     // etape 6 recherche des vosins
1151     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
1152     {
1153     MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
1154     MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
1155     MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
1156     MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
1157     tripeau->change_voisin1(voisin1);
1158     tripeau->change_voisin2(voisin2);
1159     tripeau->change_voisin3(voisin3);
1160     }
1161     // //�tape 7 - On recherche les peaux
1162     int fin;
1163     do
1164     {
1165     fin=1;
1166     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
1167     {
1168     MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
1169     if (tripeau->get_nouveau_numero()==0)
1170     {
1171     fin=0;
1172     std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
1173     lst_peau.push_back(peau);
1174     tripeau->change_nouveau_numero(1);
1175     peau->push_back(tripeau);
1176     determine_peau(peau);
1177     }
1178     }
1179     }
1180     while (fin==0);
1181 picher 233
1182 francois 224 std::cout << lst_peau.size() << " peau" << endl;
1183     for (int i=0;i<lst_peau.size();i++)
1184     std::cout << " peau " << i << " " << lst_peau[i]->size() << " triangles" << endl;
1185     //test de manifold
1186     LISTE_MG_SEGMENT::iterator itseg;
1187     for (MG_SEGMENT* seg=mai2->get_premier_segment(itseg);seg!=NULL;seg=mai2->get_suivant_segment(itseg))
1188     seg->change_nouveau_numero(0);
1189     TPL_MAP_ENTITE<MG_SEGMENT*> nonmanifoldarete;
1190     TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeud;
1191 francois 234 TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeuddepuisarete;
1192 francois 224 int nbpeau=lst_peau.size();
1193     for (int i=0;i<nbpeau;i++)
1194     {
1195     int nbtri=lst_peau[i]->size();
1196     for (int j=0;j<nbtri;j++)
1197     {
1198     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
1199     tripeau->get_segment1()->change_nouveau_numero(tripeau->get_segment1()->get_nouveau_numero()+1);
1200     tripeau->get_segment2()->change_nouveau_numero(tripeau->get_segment2()->get_nouveau_numero()+1);
1201     tripeau->get_segment3()->change_nouveau_numero(tripeau->get_segment3()->get_nouveau_numero()+1);
1202     if (tripeau->get_segment1()->get_nouveau_numero()>2)
1203     nonmanifoldarete.ajouter(tripeau->get_segment1());
1204     if (tripeau->get_segment2()->get_nouveau_numero()>2)
1205     nonmanifoldarete.ajouter(tripeau->get_segment2());
1206     if (tripeau->get_segment3()->get_nouveau_numero()>2)
1207     nonmanifoldarete.ajouter(tripeau->get_segment3());
1208     }
1209     }
1210 picher 233
1211 francois 224 int nbnonmanifoldarete=nonmanifoldarete.get_nb();
1212     for (int i=0;i<nbnonmanifoldarete;i++)
1213     {
1214     MG_SEGMENT* seg=nonmanifoldarete.get(i);
1215     MG_NOEUD* n1=seg->get_noeud1();
1216     MG_NOEUD* n2=seg->get_noeud2();
1217 francois 234 nonmanifoldnoeuddepuisarete.ajouter(n1);
1218     nonmanifoldnoeuddepuisarete.ajouter(n2);
1219 francois 224 }
1220 francois 234 LISTE_MG_NOEUD::iterator itnoeud;
1221     for (MG_NOEUD* no=mai2->get_premier_noeud(itnoeud);no!=NULL;no=mai2->get_suivant_noeud(itnoeud))
1222     {
1223     if (nonmanifoldnoeuddepuisarete.existe(no)) continue;
1224     if (est_non_manifold(no)) nonmanifoldnoeud.ajouter(no);
1225     }
1226 francois 224
1227     cout << "Non manifold par arete " << nonmanifoldarete.get_nb() << endl;
1228 picher 233 cout << "Non manifold par noeud " << nonmanifoldnoeud.get_nb() << endl;
1229 picher 230
1230 francois 234 //�tape 8 - Traitement des cas de non manifold
1231     //non manifold par ar�te
1232 picher 233 for (int i=0;i<nbnonmanifoldarete;i++)
1233     {
1234     MG_SEGMENT* segment=nonmanifoldarete.get(i);
1235 francois 234 MG_NOEUD* noeud1 = mg_mai->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
1236     MG_NOEUD* noeud2 = mg_mai->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
1237     int nb_tetra = noeud1->get_lien_tetra()->get_nb();
1238     for (int j=0;j<nb_tetra;j++)
1239 picher 233 {
1240 francois 234 MG_TETRA* mgtet =noeud1->get_lien_tetra()->get(j);
1241     int originetet=mgtet->get_origine();
1242     if (originetet == MAILLEUR_AUTO)
1243 picher 233 {
1244 francois 234 //On r�active le tetra si l'autre noeud du segment lui appartient aussi
1245     MG_NOEUD* no1 = mgtet->get_noeud1();
1246     MG_NOEUD* no2 = mgtet->get_noeud2();
1247     MG_NOEUD* no3 = mgtet->get_noeud3();
1248     MG_NOEUD* no4 = mgtet->get_noeud4();
1249     if((no1 == noeud2) ||(no2 == noeud2) || (no3 == noeud2) || (no4 == noeud2))
1250     mgtet->change_origine(OPTIMISE);
1251 picher 233 }
1252     }
1253     }
1254    
1255     //non manifold par noeud
1256 picher 231 int nbnonmanifoldnoeud=nonmanifoldnoeud.get_nb();
1257     for (int i=0;i<nbnonmanifoldnoeud;i++)
1258     {
1259 francois 234 MG_NOEUD* noeud=mg_mai->get_mg_noeudid(nonmanifoldnoeud.get(i)->get_nouveau_numero());
1260     int nb_tetra = noeud->get_lien_tetra()->get_nb();
1261     for (int j = 0;j<nb_tetra;j++)
1262 picher 231 {
1263 francois 234 MG_TETRA* mgtet =noeud->get_lien_tetra()->get(j);
1264     int originetet=mgtet->get_origine();
1265     if (originetet == MAILLEUR_AUTO)
1266 picher 231 {
1267 francois 234 mgtet->change_origine(OPTIMISE);
1268 picher 231 }
1269     }
1270     }
1271 francois 234 *mai2_id = mai2->get_id();
1272     if ((nbnonmanifoldarete != 0) || (nbnonmanifoldnoeud != 0))
1273     {
1274     for (int i=0;i<lst_peau.size();i++)
1275     {
1276     delete lst_peau[i];
1277     }
1278     lst_peau.clear();
1279     gest2.supprimer_mg_maillageid(*mai2_id);
1280     return 0;
1281     }
1282 picher 230
1283 francois 234 //�tape 9 - Suppression des peaux de tr�s faible taille (mati�re flottante)
1284 picher 233 cout << "Suppression des peaux de faible taille" << endl;
1285     int nbtri_total = mai2->get_nb_mg_triangle();
1286 picher 230 for (int i=0;i<nbpeau;i++)
1287     {
1288     int nbtri=lst_peau[i]->size();
1289 picher 233 if (nbtri < frac_min*nbtri_total)
1290 picher 230 {
1291     for (int j=0;j<nbtri;j++)
1292     {
1293     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
1294 picher 233 mai2->supprimer_mg_triangleid(tripeau->get_id());
1295 picher 230 }
1296     }
1297 picher 233 }
1298     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
1299    
1300 francois 234 return 1;
1301 francois 222 }
1302    
1303 francois 224 void MG_LISSAGE::determine_peau(std::vector<MG_TRIANGLE_PEAU*> * peau)
1304 francois 222 {
1305 francois 224 int num=0;
1306     while (num!=peau->size())
1307     {
1308     MG_TRIANGLE_PEAU* p=(*peau)[num];
1309     if (p->get_voisin1()->get_nouveau_numero()==0)
1310     {
1311     peau->push_back(p->get_voisin1());
1312     p->get_voisin1()->change_nouveau_numero(1);
1313     }
1314     if (p->get_voisin2()->get_nouveau_numero()==0)
1315     {
1316     peau->push_back(p->get_voisin2());
1317     p->get_voisin2()->change_nouveau_numero(1);
1318     }
1319     if (p->get_voisin3()->get_nouveau_numero()==0)
1320     {
1321     peau->push_back(p->get_voisin3());
1322     p->get_voisin3()->change_nouveau_numero(1);
1323     }
1324     num++;
1325     }
1326     }
1327    
1328    
1329     MG_TRIANGLE_PEAU* MG_LISSAGE::insere_triangle(MG_ELEMENT_TOPOLOGIQUE* topo,class MG_NOEUD *mgnoeud1,class MG_NOEUD *mgnoeud2,class MG_NOEUD *mgnoeud3,MG_MAILLAGE* mg_maillage,int origine)
1330     {
1331 francois 222 MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
1332     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
1333     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
1334     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
1335     if (mgsegment1==NULL)
1336 francois 224 mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,origine);
1337 francois 222 if (mgsegment2==NULL)
1338 francois 224 mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,origine);
1339 francois 222 if (mgsegment3==NULL)
1340 francois 224 mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,origine);
1341     MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
1342 francois 222 mg_maillage->ajouter_mg_triangle(mtriangle);
1343     return mtriangle;
1344     }
1345    
1346 francois 224 MG_TRIANGLE_PEAU* MG_LISSAGE::recherche_voisin(MG_NOEUD* mg_noeud1,MG_NOEUD* mg_noeud2,MG_TRIANGLE_PEAU* triref)
1347 francois 222 {
1348 francois 224 MG_TRIANGLE_PEAU* trisol=NULL;
1349     double angleref=4.*M_PI;
1350 francois 222 int nb1=mg_noeud1->get_lien_triangle()->get_nb();
1351     int nb2=mg_noeud2->get_lien_triangle()->get_nb();
1352     for (int i=0;i<nb1;i++)
1353     for (int j=0;j<nb2;j++)
1354 francois 224 if (mg_noeud1->get_lien_triangle()->get(i)==mg_noeud2->get_lien_triangle()->get(j))
1355     if ((MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i)!=triref)
1356     {
1357     double angle=calcul_angle(triref,(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i));
1358     if (angle<angleref)
1359     {
1360     angleref=angle;
1361     trisol=(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i);
1362     }
1363     }
1364     return trisol;
1365 francois 222 }
1366 francois 223
1367 francois 234 int MG_LISSAGE::est_non_manifold(MG_NOEUD* no)
1368     {
1369     static int compteur=0;
1370     compteur++;
1371     int nb_tri=no->get_lien_triangle()->get_nb();
1372     TPL_MAP_ENTITE<MG_TRIANGLE*> lst;
1373     for (int i=0;i<nb_tri;i++)
1374     lst.ajouter(no->get_lien_triangle()->get(i));
1375     MG_TRIANGLE_PEAU* premier_triangle=(MG_TRIANGLE_PEAU*)lst.get(0);
1376     lst.supprimer(premier_triangle);
1377     MG_TRIANGLE_PEAU* tricour=(MG_TRIANGLE_PEAU*)premier_triangle;
1378     int i=0;
1379     do
1380     {
1381     if (lst.existe(tricour->get_voisin1()) || ((tricour->get_voisin1()==premier_triangle)&& (i>1)))
1382     {
1383     tricour=tricour->get_voisin1();
1384     lst.supprimer(tricour);
1385     }
1386     else if (lst.existe(tricour->get_voisin2()) || ((tricour->get_voisin2()==premier_triangle)&& (i>1)))
1387     {
1388     tricour=tricour->get_voisin2();
1389     lst.supprimer(tricour);
1390     }
1391     else if (lst.existe(tricour->get_voisin3()) || ((tricour->get_voisin3()==premier_triangle)&& (i>1)))
1392     {
1393     tricour=tricour->get_voisin3();
1394     lst.supprimer(tricour);
1395     }
1396     i++;
1397     }
1398     while (tricour!=premier_triangle);
1399     if (lst.get_nb()>0)
1400     return 1;
1401     return 0;
1402     }
1403 francois 224
1404     double MG_LISSAGE::calcul_angle(MG_TRIANGLE_PEAU* ft1,MG_TRIANGLE_PEAU* ft2)
1405 francois 223 {
1406 francois 224 MG_NOEUD* noeud1=ft1->get_noeud1();
1407     MG_NOEUD* noeud2=ft1->get_noeud2();
1408     MG_NOEUD* noeud3=ft1->get_noeud3();
1409     MG_NOEUD* noeud4=ft2->get_noeud1();
1410     MG_NOEUD* noeud5=ft2->get_noeud2();
1411     MG_NOEUD* noeud6=ft2->get_noeud3();
1412    
1413     MG_NOEUD* noeuda=NULL;
1414     MG_NOEUD* noeudb=NULL;
1415     MG_NOEUD* noeudc=NULL;
1416     MG_NOEUD* noeudd=NULL;
1417    
1418    
1419     if (!(noeud1==noeud4))
1420     if (!(noeud1==noeud5))
1421     if (!(noeud1==noeud6))
1422     {
1423     noeuda=noeud1;
1424     noeudb=noeud2;
1425     noeudc=noeud3;
1426     }
1427     if (!(noeud2==noeud4))
1428     if (!(noeud2==noeud5))
1429     if (!(noeud2==noeud6))
1430     {
1431     noeuda=noeud2;
1432     noeudb=noeud1;
1433     noeudc=noeud3;
1434     }
1435     if (!(noeud3==noeud4))
1436     if (!(noeud3==noeud5))
1437     if (!(noeud3==noeud6))
1438     {
1439     noeuda=noeud3;
1440     noeudb=noeud1;
1441     noeudc=noeud2;
1442     }
1443     if (noeuda==NULL) return 2*M_PI;
1444     if (!(noeud4==noeudb))
1445     if (!(noeud4==noeudc)) noeudd=noeud4;
1446     if (!(noeud5==noeudb))
1447     if (!(noeud5==noeudc)) noeudd=noeud5;
1448     if (!(noeud6==noeudb))
1449     if (!(noeud6==noeudc)) noeudd=noeud6;
1450    
1451     OT_VECTEUR_3D milieu(0.5*(noeudb->get_x()+noeudc->get_x()),0.5*(noeudb->get_y()+noeudc->get_y()),0.5*(noeudb->get_z()+noeudc->get_z()));
1452     OT_VECTEUR_3D bc(noeudc->get_x()-noeudb->get_x(),noeudc->get_y()-noeudb->get_y(),noeudc->get_z()-noeudb->get_z());
1453     bc.norme();
1454     double xyz[3];
1455     xyz[0]=noeuda->get_x();
1456     xyz[1]=noeuda->get_y();
1457     xyz[2]=noeuda->get_z();
1458     OT_VECTEUR_3D n1(xyz,milieu.get_xyz());
1459     n1.norme();
1460     OT_VECTEUR_3D n=bc&n1;
1461     n1=n&bc;
1462     n1.norme();
1463     xyz[0]=noeudd->get_x();
1464     xyz[1]=noeudd->get_y();
1465     xyz[2]=noeudd->get_z();
1466     OT_VECTEUR_3D n2(xyz,milieu.get_xyz());
1467     n2.norme();
1468     n=bc&n2;
1469     n2=n&bc;
1470     n2.norme();
1471    
1472     double ps=n1*n2;
1473     if (ps>1.) ps=1.;
1474     if (ps<-1.) ps=-1.;
1475     double angle=acos(ps);
1476     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
1477     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
1478     OT_VECTEUR_3D nbnd(noeudd->get_x()-noeudb->get_x(),noeudd->get_y()-noeudb->get_y(),noeudd->get_z()-noeudb->get_z());
1479     n=n1n2&n1n3;
1480     n.norme();
1481     nbnd.norme();
1482    
1483     ps=n*nbnd;
1484     if (ps<-0.0001) angle=2*M_PI-angle;
1485     return angle;
1486    
1487    
1488 francois 223 }
1489 picher 230
1490     double MG_LISSAGE::ponderation_gaussian(double s,double sigma)
1491     {
1492     double w_s;
1493     w_s = exp(-(s*s)/(2.*sigma*sigma));
1494     return w_s;
1495     }
1496    
1497     double MG_LISSAGE::ponderation_laplacian(double s,double sigma)
1498     {
1499     double w_s;
1500     w_s = exp(-(sqrt(2.)*fabs(s))/sigma);
1501     return w_s;
1502     }
1503    
1504     double MG_LISSAGE::ponderation_elfallahford(double s,double sigma)
1505     {
1506     double w_s;
1507     w_s = 1./sqrt(1+pow((s/sigma),2));
1508     return w_s;
1509     }