ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mg_lissage.cpp
Revision: 233
Committed: Fri Feb 26 14:35:04 2010 UTC (15 years, 2 months ago) by picher
Original Path: magic/lib/optimisation/optimisation/src/mg_lissage.cpp
File size: 42203 byte(s)
Log Message:
Mise a jour de la classe mg_lissage et ajout de la possibilite d'importer la repartition de densite dans la classe mg_import

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     }
129     if (compteur == 0) mgtet->change_origine(OPTIMISE);
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     //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     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     void MG_LISSAGE::lisse(FEM_MAILLAGE* mai, MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max, int reactiv, int bruit, int liss,int opti,int imp,int m_auto)
354     {
355 picher 231 //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 //visualisation(mg_mai,gest2,"nm_noeud.magic");
373     coderesu = extract_skin(mg_mai,gest2,frac_min,&mai2_id);
374     //visualisation(mg_mai,gest2,"nm_noeud_resolu.magic");
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 232 affiche("Bruitage du maillage initial");
382 picher 233 bruitage(mg_mai2,gest2);
383 picher 231 }
384 francois 232 if (liss == 1)
385 picher 231 {
386 francois 232 affiche("Procédure de lissage");
387 picher 233 //visualisation2(mg_mai2,gest2,"test_lissage.magic",18615);
388     lissage(mg_mai2,gest2,epsilon,sigma,iter_max);
389 picher 231 }
390     }
391    
392 picher 230 void MG_LISSAGE::bruitage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
393     {
394     //Méthode test pour bruiter un maillage initial et tester la méthode de lissage
395     srand(time(NULL));
396     LISTE_MG_NOEUD::iterator it_no;
397     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
398     {
399     //On modifie la position des noeuds légèrement (ajout d'un bruit au maillage initial)
400     int signe;
401     int delta;
402     double delta_x,delta_y,delta_z;
403     signe = rand()%2 + 1;
404     delta = rand()%2 + 1;
405     if(signe == 1)
406     {
407 picher 233 delta_x = delta*0.6;
408 picher 230 }
409     else //(signe == 2)
410     {
411 picher 233 delta_x = -delta*0.6;
412 picher 230 }
413     signe = rand()%2 + 1;
414     delta = rand()%2 + 1;
415     if(signe == 1)
416     {
417 picher 233 delta_y = delta*0.6;
418 picher 230 }
419     else //(signe == 2)
420     {
421 picher 233 delta_y = -delta*0.6;
422 picher 230 }
423     signe = rand()%2 + 1;
424     delta = rand()%2 + 1;
425     if(signe == 1)
426     {
427 picher 233 delta_z = delta*0.6;
428 picher 230 }
429     else //(signe == 2)
430     {
431 picher 233 delta_z = -delta*0.6;
432 picher 230 }
433     double x_new = noeud_i->get_x() + delta_x;
434     double y_new = noeud_i->get_y() + delta_y;
435     double z_new = noeud_i->get_z() + delta_z;
436     noeud_i->change_x(x_new);
437     noeud_i->change_y(y_new);
438     noeud_i->change_z(z_new);
439     }
440     }
441 francois 224
442 picher 230 void MG_LISSAGE::lissage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max)
443     {
444     //Méthode de lissage inspirée de Chen(2005)
445     double un_sur_pi = 1./M_PI;
446     int compteur = 0;
447     int fin = 0;
448 picher 231 FILE *out = fopen("convergence.txt","wt");
449 francois 224
450 picher 231
451 picher 230 do
452     {
453     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
454 picher 233 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
455 picher 230 TPL_LISTE_ENTITE<double> liste_wij;
456     LISTE_MG_TRIANGLE::iterator it_tri;
457     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
458     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
459     {
460     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
461     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
462 picher 233 normal_f_i.norme();
463     liste_normales2.ajouter(normal_f_i);
464 picher 230 //Remplissage de la liste des voisins du triangle i
465     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
466     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
467     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
468     for (int j = 0;j<nb_voisins1;j++)
469     {
470     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
471     liste_voisins.ajouter(mgtri_1);
472     }
473     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
474     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
475     for (int j = 0;j<nb_voisins2;j++)
476     {
477     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
478     liste_voisins.ajouter(mgtri_2);
479     }
480     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
481     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
482     for (int j = 0;j<nb_voisins3;j++)
483     {
484     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
485     liste_voisins.ajouter(mgtri_3);
486     }
487     liste_voisins.supprimer(mgtri_i);
488     int nb_voisins = liste_voisins.get_nb();
489     double w_ij = 1./nb_voisins;
490     double phi_i_min = 10.;
491     double s_i = 0.0;
492     double phi_im = 0.0;
493 picher 231 double *phi_ij = new double[nb_voisins];
494 picher 230 OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
495     OT_VECTEUR_3D eta_i(0.,0.,0.);
496     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
497     int j = 0;
498     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
499     {
500     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
501     //1-Calculer la normale moyenne pour chaque triangle
502     normal_f_i_mean = normal_f_i_mean + w_ij*normal_f_j;
503     //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 à Nb_voisins
504     double prod_scalaire = normal_f_i*normal_f_j;
505 picher 233 if (prod_scalaire > 1.)
506     {
507     prod_scalaire = 1.;
508     }
509     if (prod_scalaire < -1.)
510     {
511     prod_scalaire = -1.;
512     }
513 picher 230 phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
514     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
515     if (phi_ij[j] < phi_i_min)
516     {
517     phi_i_min = phi_ij[j];
518     eta_i = normal_f_j;
519     }
520     //3.1-On calcul l'angle moyen phi_im
521     phi_im = phi_im + w_ij*phi_ij[j];
522     j++;
523     }
524     normal_f_i_mean.norme();
525     j = 0;
526     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
527     {
528     //3.2-Calcul de s_i selon la variance
529     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
530     j++;
531     }
532 picher 231 delete[] phi_ij;
533 picher 230 //4-On calcule une nouvelle normale pour chaque triangle
534     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
535     normal_f_i_new.norme();
536     liste_normales.ajouter(normal_f_i_new);
537     liste_wij.ajouter(w_ij);
538     mgtri->change_nouveau_numero(k);
539     k++;
540     }
541    
542     LISTE_MG_NOEUD::iterator it_no;
543     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
544     {
545     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
546     double w_ij_prime = 0.0;
547     OT_VECTEUR_3D v_temp(0.,0.,0.);
548     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
549     for(int j=0;j<nb_voisins_j;j++)
550     {
551     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
552     //On calcule le centroïde cj du triangle mgtri_j
553     MG_NOEUD* n1 = mgtri_j->get_noeud1();
554     MG_NOEUD* n2 = mgtri_j->get_noeud2();
555     MG_NOEUD* n3 = mgtri_j->get_noeud3();
556     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
557     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
558     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
559     //On forme le vecteur vi_cj
560     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
561     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
562 picher 233 // w_ij_prime correspond à la somme des aires des triangles voisins du noeuds
563     OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
564     OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
565     OT_VECTEUR_3D prodvect = AB&AC;
566     double w_ij = 0.5*prodvect.get_longueur();
567 picher 230 w_ij_prime = w_ij_prime + w_ij;
568     v_temp = v_temp + w_ij*(vi_cj*normal_f_i_new)*normal_f_i_new;
569     }
570     //5-On met à jour la position des noeuds
571     v_i = v_i + v_temp/w_ij_prime;
572 picher 233 int origine = noeud_i->get_origine();
573     if (origine != IMPOSE)
574     {
575     noeud_i->change_x(v_i.get_x());
576     noeud_i->change_y(v_i.get_y());
577     noeud_i->change_z(v_i.get_z());
578     }
579 picher 230 }
580    
581     //Critère d'arrêt de l'algorithme
582     int l=0;
583     int nb_tri = mg_mai->get_nb_mg_triangle();
584     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
585     {
586     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
587 picher 233 //OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
588     OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
589 picher 230 OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
590 picher 233 double critere = 1. - fabs(normal_f_i*normal_f_i_new);
591 picher 231 fprintf(out,"%lf\n",critere);
592 picher 230 if (critere <= epsilon) l++;
593 picher 233 if (critere >= 1.)
594     {
595     int no_id = mgtri->get_noeud1()->get_id();
596     cout << "no_id : " << no_id << endl;
597     char str[50];
598     sprintf (str, "test%d.magic", no_id);
599     //if (no_id == 18615) visualisation2(mg_mai,gest2,str,no_id);
600     }
601 picher 230 }
602 picher 233 if (nb_tri == l) fin = 1;
603 picher 230 compteur++;
604     cout << "Fin de l'itération #" << compteur << endl;
605 picher 231 fprintf(out,"\nFin de l'itération #%d\n\n",compteur);
606 picher 230 }
607     while ((fin == 0) && (compteur < iter_max));
608    
609     if (fin == 1)
610     {
611     cout << "Convergence de la méthode en " << compteur << " itérations" << endl;
612     }
613     else
614     {
615     cout << "Arrêt de la procédure après " << compteur << " itérations" << endl;
616     }
617    
618 picher 231 fclose(out);
619 picher 230 } //Fin de lissage
620    
621 picher 233 int MG_LISSAGE::extract_skin(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2,double frac_min, int *mai2_id)
622 francois 222 {
623    
624     //�tape 0 - On commence par mettre � z�ro tous les nouveau_numero des triangles et des noeuds du maillage
625     LISTE_MG_TRIANGLE::iterator it_tri;
626     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
627     {
628     mgtri->change_nouveau_numero(0);
629 francois 224 mgtri->change_origine(MAILLEUR_AUTO);
630 francois 222 }
631     LISTE_MG_NOEUD::iterator it_noeud;
632     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
633     {
634     mgnoeud->change_nouveau_numero(0);
635     }
636    
637     //�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
638 picher 233 LISTE_MG_TETRA::iterator it_tetra;
639     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
640 francois 222 {
641     int origine = mgtet->get_origine();
642 francois 224 if (origine==IMPOSE)
643 francois 222 {
644 francois 224 mgtet->get_triangle1()->change_origine(IMPOSE);
645     mgtet->get_triangle2()->change_origine(IMPOSE);
646     mgtet->get_triangle3()->change_origine(IMPOSE);
647     mgtet->get_triangle4()->change_origine(IMPOSE);
648     }
649     if (origine==OPTIMISE)
650     {
651     if (mgtet->get_triangle1()->get_origine()!=IMPOSE) mgtet->get_triangle1()->change_origine(OPTIMISE);
652     if (mgtet->get_triangle2()->get_origine()!=IMPOSE) mgtet->get_triangle2()->change_origine(OPTIMISE);
653     if (mgtet->get_triangle3()->get_origine()!=IMPOSE) mgtet->get_triangle3()->change_origine(OPTIMISE);
654     if (mgtet->get_triangle4()->get_origine()!=IMPOSE) mgtet->get_triangle4()->change_origine(OPTIMISE);
655     }
656    
657     if (((origine == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
658     ((origine == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
659     ((origine == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
660     ((origine == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
661     ((origine == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
662     ((origine == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
663    
664     {
665 francois 222 int num1 = mgtet->get_triangle1()->get_nouveau_numero();
666     int num2 = mgtet->get_triangle2()->get_nouveau_numero();
667     int num3 = mgtet->get_triangle3()->get_nouveau_numero();
668     int num4 = mgtet->get_triangle4()->get_nouveau_numero();
669     num1++;
670     num2++;
671     num3++;
672     num4++;
673     mgtet->get_triangle1()->change_nouveau_numero(num1);
674     mgtet->get_triangle2()->change_nouveau_numero(num2);
675     mgtet->get_triangle3()->change_nouveau_numero(num3);
676     mgtet->get_triangle4()->change_nouveau_numero(num4);
677     }
678     }
679    
680     //�tape 2 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour identifier les noeuds leur appartenant
681     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
682     {
683     int num = mgtri->get_nouveau_numero();
684     if (num == 1)
685     {
686     MG_NOEUD* noeud1 = mgtri->get_noeud1();
687     MG_NOEUD* noeud2 = mgtri->get_noeud2();
688     MG_NOEUD* noeud3 = mgtri->get_noeud3();
689     noeud1->change_nouveau_numero(1);
690     noeud2->change_nouveau_numero(1);
691     noeud3->change_nouveau_numero(1);
692     }
693     }
694    
695 picher 233 //�tape 3 - On crée un nouveau maillage pour la peau
696 francois 222 MG_MAILLAGE* mai2 = new MG_MAILLAGE(NULL);
697     gest2.ajouter_mg_maillage(mai2);
698    
699     // //�tape 4 - On boucle l'ensemble des noeuds identifi�s � l'�tape 2 pour les recr�er dans le second maillage
700     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
701     {
702     int num = mgnoeud->get_nouveau_numero();
703     if (num == 1)
704     {
705     double x = mgnoeud->get_x();
706     double y = mgnoeud->get_y();
707     double z = mgnoeud->get_z();
708 francois 224 int origine=TRIANGULATION;
709     if (mgnoeud->get_origine()==IMPOSE) origine=IMPOSE;
710     MG_NOEUD* noeud1 = new MG_NOEUD(NULL,x,y,z,origine);
711 francois 222 mai2->ajouter_mg_noeud(noeud1);
712     mgnoeud->change_nouveau_numero(noeud1->get_id());
713     }
714     }
715    
716     // //�tape 5 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour les recr�er dans le maillage 2
717 picher 233 for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
718 francois 222 {
719 francois 224 int originetet=mgtet->get_origine();
720     if (((originetet == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
721     ((originetet == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
722     ((originetet == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
723     ((originetet == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
724     ((originetet == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
725     ((originetet == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
726 francois 222 {
727 francois 224 MG_NOEUD* noeud1 = mgtet->get_noeud1();
728     MG_NOEUD* noeud2 = mgtet->get_noeud2();
729     MG_NOEUD* noeud3 = mgtet->get_noeud3();
730     MG_NOEUD* noeud4 = mgtet->get_noeud4();
731 francois 222 MG_NOEUD* node1 = mai2->get_mg_noeudid(noeud1->get_nouveau_numero());
732     MG_NOEUD* node2 = mai2->get_mg_noeudid(noeud2->get_nouveau_numero());
733     MG_NOEUD* node3 = mai2->get_mg_noeudid(noeud3->get_nouveau_numero());
734 francois 224 MG_NOEUD* node4 = mai2->get_mg_noeudid(noeud4->get_nouveau_numero());
735     MG_TRIANGLE* tri1=mgtet->get_triangle1();
736     MG_TRIANGLE* tri2=mgtet->get_triangle2();
737     MG_TRIANGLE* tri3=mgtet->get_triangle3();
738     MG_TRIANGLE* tri4=mgtet->get_triangle4();
739     if (tri1->get_nouveau_numero()==1)
740     {
741     int origine=TRIANGULATION;
742     if (tri1->get_origine()==IMPOSE) origine=IMPOSE;
743     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node2,node3,mai2,origine);
744     tripeau->change_nouveau_numero(0);
745     tri1->change_nouveau_numero(0);
746     }
747     if (tri2->get_nouveau_numero()==1)
748     {
749     int origine=TRIANGULATION;
750     if (tri2->get_origine()==IMPOSE) origine=IMPOSE;
751     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node4,node2,mai2,origine);
752     tripeau->change_nouveau_numero(0);
753     tri2->change_nouveau_numero(0);
754     }
755     if (tri3->get_nouveau_numero()==1)
756     {
757     int origine=TRIANGULATION;
758     if (tri3->get_origine()==IMPOSE) origine=IMPOSE;
759     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node2,node4,node3,mai2,origine);
760     tripeau->change_nouveau_numero(0);
761     tri3->change_nouveau_numero(0);
762     }
763     if (tri4->get_nouveau_numero()==1)
764     {
765     int origine=TRIANGULATION;
766     if (tri4->get_origine()==IMPOSE) origine=IMPOSE;
767     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node3,node4,mai2,origine);
768     tripeau->change_nouveau_numero(0);
769     tri4->change_nouveau_numero(0);
770     }
771 francois 222 }
772 francois 224 }
773     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
774     // etape 6 recherche des vosins
775     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
776     {
777     MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
778     MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
779     MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
780     MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
781     tripeau->change_voisin1(voisin1);
782     tripeau->change_voisin2(voisin2);
783     tripeau->change_voisin3(voisin3);
784     }
785     // //�tape 7 - On recherche les peaux
786     int fin;
787     do
788     {
789     fin=1;
790     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
791     {
792     MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
793     if (tripeau->get_nouveau_numero()==0)
794     {
795     fin=0;
796     std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
797     lst_peau.push_back(peau);
798     tripeau->change_nouveau_numero(1);
799     peau->push_back(tripeau);
800     determine_peau(peau);
801     }
802     }
803     }
804     while (fin==0);
805 picher 233
806 francois 224 std::cout << lst_peau.size() << " peau" << endl;
807     for (int i=0;i<lst_peau.size();i++)
808     std::cout << " peau " << i << " " << lst_peau[i]->size() << " triangles" << endl;
809     //test de manifold
810     for (MG_NOEUD* mgnoeud=mai2->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mai2->get_suivant_noeud(it_noeud))
811     mgnoeud->change_nouveau_numero(0);
812     LISTE_MG_SEGMENT::iterator itseg;
813     for (MG_SEGMENT* seg=mai2->get_premier_segment(itseg);seg!=NULL;seg=mai2->get_suivant_segment(itseg))
814     seg->change_nouveau_numero(0);
815     TPL_MAP_ENTITE<MG_SEGMENT*> nonmanifoldarete;
816     TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeud;
817     int nbpeau=lst_peau.size();
818     for (int i=0;i<nbpeau;i++)
819     {
820     int nbtri=lst_peau[i]->size();
821     for (int j=0;j<nbtri;j++)
822     {
823     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
824     tripeau->get_segment1()->change_nouveau_numero(tripeau->get_segment1()->get_nouveau_numero()+1);
825     tripeau->get_segment2()->change_nouveau_numero(tripeau->get_segment2()->get_nouveau_numero()+1);
826     tripeau->get_segment3()->change_nouveau_numero(tripeau->get_segment3()->get_nouveau_numero()+1);
827     if (tripeau->get_segment1()->get_nouveau_numero()>2)
828     nonmanifoldarete.ajouter(tripeau->get_segment1());
829     if (tripeau->get_segment2()->get_nouveau_numero()>2)
830     nonmanifoldarete.ajouter(tripeau->get_segment2());
831     if (tripeau->get_segment3()->get_nouveau_numero()>2)
832     nonmanifoldarete.ajouter(tripeau->get_segment3());
833     if ((tripeau->get_noeud1()->get_nouveau_numero()!=0) && (tripeau->get_noeud1()->get_nouveau_numero()!=i+1))
834     nonmanifoldnoeud.ajouter(tripeau->get_noeud1());
835     if ((tripeau->get_noeud2()->get_nouveau_numero()!=0) && (tripeau->get_noeud2()->get_nouveau_numero()!=i+1))
836     nonmanifoldnoeud.ajouter(tripeau->get_noeud2());
837     if ((tripeau->get_noeud3()->get_nouveau_numero()!=0) && (tripeau->get_noeud3()->get_nouveau_numero()!=i+1))
838     nonmanifoldnoeud.ajouter(tripeau->get_noeud3());
839     tripeau->get_noeud1()->change_nouveau_numero(i+1);
840     tripeau->get_noeud2()->change_nouveau_numero(i+1);
841     tripeau->get_noeud3()->change_nouveau_numero(i+1);
842     }
843     }
844 picher 233
845 francois 224 int nbnonmanifoldarete=nonmanifoldarete.get_nb();
846     for (int i=0;i<nbnonmanifoldarete;i++)
847     {
848     MG_SEGMENT* seg=nonmanifoldarete.get(i);
849     MG_NOEUD* n1=seg->get_noeud1();
850     MG_NOEUD* n2=seg->get_noeud2();
851     nonmanifoldnoeud.supprimer(n1);
852     nonmanifoldnoeud.supprimer(n2);
853     }
854 francois 222
855 francois 224
856     cout << "Non manifold par arete " << nonmanifoldarete.get_nb() << endl;
857 picher 233 cout << "Non manifold par noeud " << nonmanifoldnoeud.get_nb() << endl;
858 picher 230
859     //Étape 8 - Traitement des cas de non manifold
860 picher 233 //non manifold par arête
861     for (int i=0;i<nbnonmanifoldarete;i++)
862     {
863     MG_SEGMENT* segment=nonmanifoldarete.get(i);
864     MG_NOEUD* noeud1 = segment->get_noeud1();
865     int id = noeud1->get_id();
866     MG_NOEUD* noeud2 = segment->get_noeud2();
867     int id2 = noeud2->get_id();
868     int id_1;
869     int id_2;
870     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
871     {
872     int num = mgnoeud->get_nouveau_numero();
873     if (num == id)
874     {
875     id_1 = mgnoeud->get_id();
876     int nb_tetra = mgnoeud->get_lien_tetra()->get_nb();
877     for (int j=0;j<nb_tetra;j++)
878     {
879     MG_TETRA* mgtet = (MG_TETRA*) mgnoeud->get_lien_tetra()->get(j);
880     int originetet=mgtet->get_origine();
881     if (originetet == MAILLEUR_AUTO)
882     {
883     //On réactive le tetra si l'autre noeud du segment lui appartient aussi
884     for (MG_NOEUD* mgnoeud2=mg_mai->get_premier_noeud(it_noeud);mgnoeud2!=NULL;mgnoeud2=mg_mai->get_suivant_noeud(it_noeud))
885     {
886     int num2 = mgnoeud2->get_nouveau_numero();
887     if (num2 == id2)
888     {
889     id_2 = mgnoeud2->get_id();
890     MG_NOEUD* no1 = mgtet->get_noeud1();
891     MG_NOEUD* no2 = mgtet->get_noeud2();
892     MG_NOEUD* no3 = mgtet->get_noeud3();
893     MG_NOEUD* no4 = mgtet->get_noeud4();
894     if((no1 == mgnoeud2) ||(no2 == mgnoeud2) || (no3 == mgnoeud2) || (no4 == mgnoeud2))
895     {
896     mgtet->change_origine(OPTIMISE);
897     }
898     }
899     }
900     }
901     }
902     }
903     }
904     //cout << "no1_id : " << id_1 << " " << "no2_id : " << id_2 << endl;
905     }
906    
907     //non manifold par noeud
908 picher 231 int nbnonmanifoldnoeud=nonmanifoldnoeud.get_nb();
909 picher 233 //if (nbnonmanifoldnoeud > 1) visualisation2(mai2,gest2,"nm_noeud.magic",39199);
910 picher 231 for (int i=0;i<nbnonmanifoldnoeud;i++)
911     {
912 picher 233 MG_NOEUD* noeud=nonmanifoldnoeud.get(i);
913     int id = noeud->get_id();
914     int id_1;
915     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
916 picher 231 {
917 picher 233 int num = mgnoeud->get_nouveau_numero();
918     if (num == id)
919 picher 231 {
920 picher 233 id_1 = mgnoeud->get_id();
921     int nb_tetra = mgnoeud->get_lien_tetra()->get_nb();
922     for (int j = 0;j<nb_tetra;j++)
923 picher 231 {
924 picher 233 MG_TETRA* mgtet = (MG_TETRA*) mgnoeud->get_lien_tetra()->get(j);
925     int originetet=mgtet->get_origine();
926     if (originetet == MAILLEUR_AUTO)
927     {
928     //On réactive le tetra en question
929     int id2 = mgtet->get_id();
930     mgtet->change_origine(OPTIMISE);
931     }
932 picher 231 }
933     }
934 picher 233 //if (num == 20467) cout << "test : " << mgnoeud->get_id() << endl;
935 picher 231 }
936 picher 233 cout << "id1 : " << id << endl;
937     cout << "id_1 : " << id_1 << endl;
938 picher 231 }
939 picher 233 //if (nbnonmanifoldnoeud > 1) visualisation2(mai2,gest2,"nm_noeud_resolu.magic",39199);
940     // for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
941     // {
942     // int id = mgnoeud->get_id();
943     // if (id == 37985)
944     // {
945     // int id1 = mgnoeud->get_nouveau_numero();
946     // cout << "nouv_id : " << id1 << endl;
947     // }
948     // }
949 picher 230
950 picher 233 //if ((nbnonmanifoldnoeud == 0) && (nbnonmanifoldarete == 0)) visualisation2(mai2,gest2,"nm_noeud_resolu.magic",39199);
951 picher 231
952 picher 233 //Étape 9 - Suppression des peaux de très faible taille (matière flottante)
953     cout << "Suppression des peaux de faible taille" << endl;
954     int nbtri_total = mai2->get_nb_mg_triangle();
955 picher 230 for (int i=0;i<nbpeau;i++)
956     {
957     int nbtri=lst_peau[i]->size();
958 picher 233 if (nbtri < frac_min*nbtri_total)
959 picher 230 {
960     for (int j=0;j<nbtri;j++)
961     {
962     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
963 picher 233 mai2->supprimer_mg_triangleid(tripeau->get_id());
964 picher 230 }
965     }
966 picher 233 }
967     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
968    
969     *mai2_id = mai2->get_id();
970     if ((nbnonmanifoldarete == 0) && (nbnonmanifoldnoeud == 0))
971     {
972     for (int i=0;i<lst_peau.size();i++)
973     {
974     delete lst_peau[i];
975     }
976     lst_peau.clear();
977     return 1;
978     }
979     else
980     {
981     for (int i=0;i<lst_peau.size();i++)
982     {
983     delete lst_peau[i];
984     }
985     lst_peau.clear();
986     gest2.supprimer_mg_maillageid(*mai2_id);
987     return 0;
988     }
989 francois 222 }
990    
991 francois 224 void MG_LISSAGE::determine_peau(std::vector<MG_TRIANGLE_PEAU*> * peau)
992 francois 222 {
993 francois 224 int num=0;
994     while (num!=peau->size())
995     {
996     MG_TRIANGLE_PEAU* p=(*peau)[num];
997     if (p->get_voisin1()->get_nouveau_numero()==0)
998     {
999     peau->push_back(p->get_voisin1());
1000     p->get_voisin1()->change_nouveau_numero(1);
1001     }
1002     if (p->get_voisin2()->get_nouveau_numero()==0)
1003     {
1004     peau->push_back(p->get_voisin2());
1005     p->get_voisin2()->change_nouveau_numero(1);
1006     }
1007     if (p->get_voisin3()->get_nouveau_numero()==0)
1008     {
1009     peau->push_back(p->get_voisin3());
1010     p->get_voisin3()->change_nouveau_numero(1);
1011     }
1012     num++;
1013     }
1014     }
1015    
1016    
1017     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)
1018     {
1019 francois 222 MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
1020     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
1021     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
1022     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
1023     if (mgsegment1==NULL)
1024 francois 224 mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,origine);
1025 francois 222 if (mgsegment2==NULL)
1026 francois 224 mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,origine);
1027 francois 222 if (mgsegment3==NULL)
1028 francois 224 mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,origine);
1029     MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
1030 francois 222 mg_maillage->ajouter_mg_triangle(mtriangle);
1031     return mtriangle;
1032     }
1033    
1034 francois 224 MG_TRIANGLE_PEAU* MG_LISSAGE::recherche_voisin(MG_NOEUD* mg_noeud1,MG_NOEUD* mg_noeud2,MG_TRIANGLE_PEAU* triref)
1035 francois 222 {
1036 francois 224 MG_TRIANGLE_PEAU* trisol=NULL;
1037     double angleref=4.*M_PI;
1038 francois 222 int nb1=mg_noeud1->get_lien_triangle()->get_nb();
1039     int nb2=mg_noeud2->get_lien_triangle()->get_nb();
1040     for (int i=0;i<nb1;i++)
1041     for (int j=0;j<nb2;j++)
1042 francois 224 if (mg_noeud1->get_lien_triangle()->get(i)==mg_noeud2->get_lien_triangle()->get(j))
1043     if ((MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i)!=triref)
1044     {
1045     double angle=calcul_angle(triref,(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i));
1046     if (angle<angleref)
1047     {
1048     angleref=angle;
1049     trisol=(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i);
1050     }
1051     }
1052     return trisol;
1053 francois 222 }
1054 francois 223
1055 francois 224
1056    
1057     double MG_LISSAGE::calcul_angle(MG_TRIANGLE_PEAU* ft1,MG_TRIANGLE_PEAU* ft2)
1058 francois 223 {
1059 francois 224 MG_NOEUD* noeud1=ft1->get_noeud1();
1060     MG_NOEUD* noeud2=ft1->get_noeud2();
1061     MG_NOEUD* noeud3=ft1->get_noeud3();
1062     MG_NOEUD* noeud4=ft2->get_noeud1();
1063     MG_NOEUD* noeud5=ft2->get_noeud2();
1064     MG_NOEUD* noeud6=ft2->get_noeud3();
1065    
1066     MG_NOEUD* noeuda=NULL;
1067     MG_NOEUD* noeudb=NULL;
1068     MG_NOEUD* noeudc=NULL;
1069     MG_NOEUD* noeudd=NULL;
1070    
1071    
1072     if (!(noeud1==noeud4))
1073     if (!(noeud1==noeud5))
1074     if (!(noeud1==noeud6))
1075     {
1076     noeuda=noeud1;
1077     noeudb=noeud2;
1078     noeudc=noeud3;
1079     }
1080     if (!(noeud2==noeud4))
1081     if (!(noeud2==noeud5))
1082     if (!(noeud2==noeud6))
1083     {
1084     noeuda=noeud2;
1085     noeudb=noeud1;
1086     noeudc=noeud3;
1087     }
1088     if (!(noeud3==noeud4))
1089     if (!(noeud3==noeud5))
1090     if (!(noeud3==noeud6))
1091     {
1092     noeuda=noeud3;
1093     noeudb=noeud1;
1094     noeudc=noeud2;
1095     }
1096     if (noeuda==NULL) return 2*M_PI;
1097     if (!(noeud4==noeudb))
1098     if (!(noeud4==noeudc)) noeudd=noeud4;
1099     if (!(noeud5==noeudb))
1100     if (!(noeud5==noeudc)) noeudd=noeud5;
1101     if (!(noeud6==noeudb))
1102     if (!(noeud6==noeudc)) noeudd=noeud6;
1103    
1104     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()));
1105     OT_VECTEUR_3D bc(noeudc->get_x()-noeudb->get_x(),noeudc->get_y()-noeudb->get_y(),noeudc->get_z()-noeudb->get_z());
1106     bc.norme();
1107     double xyz[3];
1108     xyz[0]=noeuda->get_x();
1109     xyz[1]=noeuda->get_y();
1110     xyz[2]=noeuda->get_z();
1111     OT_VECTEUR_3D n1(xyz,milieu.get_xyz());
1112     n1.norme();
1113     OT_VECTEUR_3D n=bc&n1;
1114     n1=n&bc;
1115     n1.norme();
1116     xyz[0]=noeudd->get_x();
1117     xyz[1]=noeudd->get_y();
1118     xyz[2]=noeudd->get_z();
1119     OT_VECTEUR_3D n2(xyz,milieu.get_xyz());
1120     n2.norme();
1121     n=bc&n2;
1122     n2=n&bc;
1123     n2.norme();
1124    
1125     double ps=n1*n2;
1126     if (ps>1.) ps=1.;
1127     if (ps<-1.) ps=-1.;
1128     double angle=acos(ps);
1129     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
1130     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
1131     OT_VECTEUR_3D nbnd(noeudd->get_x()-noeudb->get_x(),noeudd->get_y()-noeudb->get_y(),noeudd->get_z()-noeudb->get_z());
1132     n=n1n2&n1n3;
1133     n.norme();
1134     nbnd.norme();
1135    
1136     ps=n*nbnd;
1137     if (ps<-0.0001) angle=2*M_PI-angle;
1138     return angle;
1139    
1140    
1141 francois 223 }
1142 picher 230
1143     double MG_LISSAGE::ponderation_gaussian(double s,double sigma)
1144     {
1145     double w_s;
1146     w_s = exp(-(s*s)/(2.*sigma*sigma));
1147     return w_s;
1148     }
1149    
1150     double MG_LISSAGE::ponderation_laplacian(double s,double sigma)
1151     {
1152     double w_s;
1153     w_s = exp(-(sqrt(2.)*fabs(s))/sigma);
1154     return w_s;
1155     }
1156    
1157     double MG_LISSAGE::ponderation_elfallahford(double s,double sigma)
1158     {
1159     double w_s;
1160     w_s = 1./sqrt(1+pow((s/sigma),2));
1161     return w_s;
1162     }