ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mg_lissage.cpp
Revision: 335
Committed: Tue May 8 16:29:46 2012 UTC (13 years ago) by francois
File size: 53442 byte(s)
Log Message:
Generalisation du calcul de l'angle entre 2 triangles orientés

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