ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_posttraitement.cpp
Revision: 234
Committed: Fri Feb 26 22:20:15 2010 UTC (15 years, 2 months ago) by francois
Original Path: magic/lib/optimisation/optimisation/src/mg_lissage.cpp
File size: 41656 byte(s)
Log Message:
Correction du non manifold par noeud dans le lissage plus resolution de bug divers

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 francois 234 //MG_NOEUD* noeud_nm2 = mg_mai->get_mg_noeudid(12437);//placer ici le ID d'un second noeud non manifold (pour une ar�te)
189 picher 233 TPL_MAP_ENTITE<MG_TETRA*> liste_voisins;
190     int nb_voisins1 = noeud_nm1->get_lien_tetra()->get_nb();
191     for (int j = 0;j<nb_voisins1;j++)
192     {
193     MG_TETRA* mgtet_1 = (MG_TETRA*) noeud_nm1->get_lien_tetra()->get(j);
194     int origine = mgtet_1->get_origine();
195     if ((origine == OPTIMISE) || (origine == IMPOSE))
196     {
197     liste_voisins.ajouter(mgtet_1);
198     MG_NOEUD* no1 = mgtet_1->get_noeud1();
199     int nb_vois1 = no1->get_lien_tetra()->get_nb();
200     for (int k = 0;k<nb_vois1;k++)
201     {
202     MG_TETRA* mgtet1 = (MG_TETRA*) no1->get_lien_tetra()->get(k);
203     int origine1 = mgtet1->get_origine();
204     if ((origine1 == OPTIMISE) || (origine1 == IMPOSE)) liste_voisins.ajouter(mgtet1);
205     }
206     MG_NOEUD* no2 = mgtet_1->get_noeud2();
207     int nb_vois2 = no2->get_lien_tetra()->get_nb();
208     for (int k = 0;k<nb_vois2;k++)
209     {
210     MG_TETRA* mgtet2 = (MG_TETRA*) no2->get_lien_tetra()->get(k);
211     int origine2 = mgtet2->get_origine();
212     if ((origine2 == OPTIMISE) || (origine2 == IMPOSE)) liste_voisins.ajouter(mgtet2);
213     }
214     MG_NOEUD* no3 = mgtet_1->get_noeud3();
215     int nb_vois3 = no3->get_lien_tetra()->get_nb();
216     for (int k = 0;k<nb_vois3;k++)
217     {
218     MG_TETRA* mgtet3 = (MG_TETRA*) no3->get_lien_tetra()->get(k);
219     int origine3 = mgtet3->get_origine();
220     if ((origine3 == OPTIMISE) || (origine3 == IMPOSE)) liste_voisins.ajouter(mgtet3);
221     }
222     MG_NOEUD* no4 = mgtet_1->get_noeud4();
223     int nb_vois4 = no4->get_lien_tetra()->get_nb();
224     for (int k = 0;k<nb_vois4;k++)
225     {
226     MG_TETRA* mgtet4 = (MG_TETRA*) no4->get_lien_tetra()->get(k);
227     int origine4 = mgtet4->get_origine();
228     if ((origine4 == OPTIMISE) || (origine4 == IMPOSE)) liste_voisins.ajouter(mgtet4);
229     }
230     }
231     }
232     // int nb_voisins2 = noeud_nm2->get_lien_tetra()->get_nb();
233     // for (int j = 0;j<nb_voisins2;j++)
234     // {
235     // MG_TETRA* mgtet_1 = (MG_TETRA*) noeud_nm2->get_lien_tetra()->get(j);
236     // int origine = mgtet_1->get_origine();
237     // if ((origine == OPTIMISE) || (origine == IMPOSE))
238     // {
239     // liste_voisins.ajouter(mgtet_1);
240     // // MG_NOEUD* no1 = mgtet_1->get_noeud1();
241     // // int nb_vois1 = no1->get_lien_tetra()->get_nb();
242     // // for (int k = 0;k<nb_vois1;k++)
243     // // {
244     // // MG_TETRA* mgtet1 = (MG_TETRA*) no1->get_lien_tetra()->get(k);
245     // // int origine1 = mgtet1->get_origine();
246     // // if ((origine1 == OPTIMISE) || (origine1 == IMPOSE)) liste_voisins.ajouter(mgtet1);
247     // // }
248     // // MG_NOEUD* no2 = mgtet_1->get_noeud2();
249     // // int nb_vois2 = no2->get_lien_tetra()->get_nb();
250     // // for (int k = 0;k<nb_vois2;k++)
251     // // {
252     // // MG_TETRA* mgtet2 = (MG_TETRA*) no2->get_lien_tetra()->get(k);
253     // // int origine2 = mgtet2->get_origine();
254     // // if ((origine2 == OPTIMISE) || (origine2 == IMPOSE)) liste_voisins.ajouter(mgtet2);
255     // // }
256     // // MG_NOEUD* no3 = mgtet_1->get_noeud3();
257     // // int nb_vois3 = no3->get_lien_tetra()->get_nb();
258     // // for (int k = 0;k<nb_vois3;k++)
259     // // {
260     // // MG_TETRA* mgtet3 = (MG_TETRA*) no3->get_lien_tetra()->get(k);
261     // // int origine3 = mgtet3->get_origine();
262     // // if ((origine3 == OPTIMISE) || (origine3 == IMPOSE)) liste_voisins.ajouter(mgtet3);
263     // // }
264     // // MG_NOEUD* no4 = mgtet_1->get_noeud4();
265     // // int nb_vois4 = no4->get_lien_tetra()->get_nb();
266     // // for (int k = 0;k<nb_vois4;k++)
267     // // {
268     // // MG_TETRA* mgtet4 = (MG_TETRA*) no4->get_lien_tetra()->get(k);
269     // // int origine4 = mgtet4->get_origine();
270     // // if ((origine4 == OPTIMISE) || (origine4 == IMPOSE)) liste_voisins.ajouter(mgtet4);
271     // // }
272     // }
273     // }
274     int nb_voisins = liste_voisins.get_nb();
275     TPL_MAP_ENTITE<MG_TETRA*>::ITERATEUR it;
276     for(MG_TETRA* mgtet_j=liste_voisins.get_premier(it);mgtet_j!=NULL;mgtet_j=liste_voisins.get_suivant(it))
277     {
278     MG_NOEUD* noeud1 = mgtet_j->get_noeud1();
279     double xx1 = noeud1->get_x();
280     double yy1 = noeud1->get_y();
281     double zz1 = noeud1->get_z();
282     int ori_no1 = noeud1->get_origine();
283     MG_NOEUD* noeud11 = new MG_NOEUD(NULL,xx1,yy1,zz1,ori_no1);
284     MG_NOEUD* noeud2 = mgtet_j->get_noeud2();
285     double xx2 = noeud2->get_x();
286     double yy2 = noeud2->get_y();
287     double zz2 = noeud2->get_z();
288     int ori_no2 = noeud2->get_origine();
289     MG_NOEUD* noeud22 = new MG_NOEUD(NULL,xx2,yy2,zz2,ori_no2);
290     MG_NOEUD* noeud3 = mgtet_j->get_noeud3();
291     double xx3 = noeud3->get_x();
292     double yy3 = noeud3->get_y();
293     double zz3 = noeud3->get_z();
294     int ori_no3 = noeud3->get_origine();
295     MG_NOEUD* noeud33 = new MG_NOEUD(NULL,xx3,yy3,zz3,ori_no3);
296     MG_NOEUD* noeud4 = mgtet_j->get_noeud4();
297     double xx4 = noeud4->get_x();
298     double yy4 = noeud4->get_y();
299     double zz4 = noeud4->get_z();
300     int ori_no4 = noeud4->get_origine();
301     MG_NOEUD* noeud44 = new MG_NOEUD(NULL,xx4,yy4,zz4,ori_no4);
302     mai3->ajouter_mg_noeud(noeud11);
303     mai3->ajouter_mg_noeud(noeud22);
304     mai3->ajouter_mg_noeud(noeud33);
305     mai3->ajouter_mg_noeud(noeud44);
306     MG_TRIANGLE* tri1 = mgtet_j->get_triangle1();
307     MG_NOEUD* no11 = tri1->get_noeud1();
308     MG_NOEUD* no12 = tri1->get_noeud2();
309     MG_NOEUD* no13 = tri1->get_noeud3();
310     MG_SEGMENT* seg11 = tri1->get_segment1();
311     MG_SEGMENT* seg12 = tri1->get_segment2();
312     MG_SEGMENT* seg13 = tri1->get_segment3();
313     int ori_tri1 = tri1->get_origine();
314     MG_TRIANGLE* tri11 = new MG_TRIANGLE(NULL,no11,no12,no13,seg11,seg12,seg13,ori_tri1);
315     MG_TRIANGLE* tri2 = mgtet_j->get_triangle2();
316     MG_NOEUD* no21 = tri2->get_noeud1();
317     MG_NOEUD* no22 = tri2->get_noeud2();
318     MG_NOEUD* no23 = tri2->get_noeud3();
319     MG_SEGMENT* seg21 = tri2->get_segment1();
320     MG_SEGMENT* seg22 = tri2->get_segment2();
321     MG_SEGMENT* seg23 = tri2->get_segment3();
322     int ori_tri2 = tri2->get_origine();
323     MG_TRIANGLE* tri22 = new MG_TRIANGLE(NULL,no21,no22,no23,seg21,seg22,seg23,ori_tri2);
324     MG_TRIANGLE* tri3 = mgtet_j->get_triangle3();
325     MG_NOEUD* no31 = tri3->get_noeud1();
326     MG_NOEUD* no32 = tri3->get_noeud2();
327     MG_NOEUD* no33 = tri3->get_noeud3();
328     MG_SEGMENT* seg31 = tri3->get_segment1();
329     MG_SEGMENT* seg32 = tri3->get_segment2();
330     MG_SEGMENT* seg33 = tri3->get_segment3();
331     int ori_tri3 = tri3->get_origine();
332     MG_TRIANGLE* tri33 = new MG_TRIANGLE(NULL,no31,no32,no33,seg31,seg32,seg33,ori_tri3);
333     MG_TRIANGLE* tri4 = mgtet_j->get_triangle4();
334     MG_NOEUD* no41 = tri4->get_noeud1();
335     MG_NOEUD* no42 = tri4->get_noeud2();
336     MG_NOEUD* no43 = tri4->get_noeud3();
337     MG_SEGMENT* seg41 = tri4->get_segment1();
338     MG_SEGMENT* seg42 = tri4->get_segment2();
339     MG_SEGMENT* seg43 = tri4->get_segment3();
340     int ori_tri4 = tri4->get_origine();
341     MG_TRIANGLE* tri44 = new MG_TRIANGLE(NULL,no41,no42,no43,seg41,seg42,seg43,ori_tri4);
342     mai3->ajouter_mg_triangle(tri11);
343     mai3->ajouter_mg_triangle(tri22);
344     mai3->ajouter_mg_triangle(tri33);
345     mai3->ajouter_mg_triangle(tri44);
346     int ori_tet = mgtet_j->get_origine();
347     MG_TETRA* mgtet = new MG_TETRA(NULL,noeud11,noeud22,noeud33,noeud44,tri11,tri22,tri33,tri44,ori_tet);
348     mai3->ajouter_mg_tetra(mgtet);
349     }
350     gest3.enregistrer(nomfichier);
351     }
352    
353     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 francois 234 //Menu de la m�thode de lissage
356 francois 232 affiche("Extraction du maillage de surface");
357 picher 231 int coderesu = 0;
358 picher 233 double frac_min = 0.005;
359     int mai2_id;
360 francois 232 if (opti==1) conserve(OPTIMISE);
361     if (imp==1) conserve(IMPOSE);
362 picher 233 if (m_auto==1) conserve(MAILLEUR_AUTO);
363     copieorigine(mai,gest2);
364     MG_MAILLAGE* mg_mai = (MG_MAILLAGE*)mai->get_mg_maillage();
365     gain_poids(mg_mai,gest2);
366     if (reactiv == 1)
367     {
368     reactivation(mg_mai,gest2);
369     }
370 picher 231 do
371     {
372 picher 233 //visualisation(mg_mai,gest2,"nm_noeud.magic");
373     coderesu = extract_skin(mg_mai,gest2,frac_min,&mai2_id);
374 francois 234 gain_poids(mg_mai,gest2);
375 picher 233 //visualisation(mg_mai,gest2,"nm_noeud_resolu.magic");
376 picher 231 }
377     while (coderesu == 0);
378 picher 233
379     MG_MAILLAGE* mg_mai2=gest2.get_mg_maillageid(mai2_id);
380 francois 232 if (bruit == 1)
381 picher 231 {
382 francois 232 affiche("Bruitage du maillage initial");
383 picher 233 bruitage(mg_mai2,gest2);
384 picher 231 }
385 francois 232 if (liss == 1)
386 picher 231 {
387 francois 234 affiche("Procedure de lissage");
388 picher 233 //visualisation2(mg_mai2,gest2,"test_lissage.magic",18615);
389     lissage(mg_mai2,gest2,epsilon,sigma,iter_max);
390 picher 231 }
391     }
392    
393 picher 230 void MG_LISSAGE::bruitage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
394     {
395 francois 234 //M�thode test pour bruiter un maillage initial et tester la m�thode de lissage
396 picher 230 srand(time(NULL));
397     LISTE_MG_NOEUD::iterator it_no;
398     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
399     {
400 francois 234 //On modifie la position des noeuds l�g�rement (ajout d'un bruit au maillage initial)
401 picher 230 int signe;
402     int delta;
403     double delta_x,delta_y,delta_z;
404     signe = rand()%2 + 1;
405     delta = rand()%2 + 1;
406     if(signe == 1)
407     {
408 picher 233 delta_x = delta*0.6;
409 picher 230 }
410     else //(signe == 2)
411     {
412 picher 233 delta_x = -delta*0.6;
413 picher 230 }
414     signe = rand()%2 + 1;
415     delta = rand()%2 + 1;
416     if(signe == 1)
417     {
418 picher 233 delta_y = delta*0.6;
419 picher 230 }
420     else //(signe == 2)
421     {
422 picher 233 delta_y = -delta*0.6;
423 picher 230 }
424     signe = rand()%2 + 1;
425     delta = rand()%2 + 1;
426     if(signe == 1)
427     {
428 picher 233 delta_z = delta*0.6;
429 picher 230 }
430     else //(signe == 2)
431     {
432 picher 233 delta_z = -delta*0.6;
433 picher 230 }
434     double x_new = noeud_i->get_x() + delta_x;
435     double y_new = noeud_i->get_y() + delta_y;
436     double z_new = noeud_i->get_z() + delta_z;
437     noeud_i->change_x(x_new);
438     noeud_i->change_y(y_new);
439     noeud_i->change_z(z_new);
440     }
441     }
442 francois 224
443 picher 230 void MG_LISSAGE::lissage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max)
444     {
445 francois 234 //M�thode de lissage inspir�e de Chen(2005)
446 picher 230 double un_sur_pi = 1./M_PI;
447     int compteur = 0;
448     int fin = 0;
449 picher 231 FILE *out = fopen("convergence.txt","wt");
450 francois 224
451 picher 231
452 picher 230 do
453     {
454     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
455 picher 233 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
456 picher 230 TPL_LISTE_ENTITE<double> liste_wij;
457     LISTE_MG_TRIANGLE::iterator it_tri;
458     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
459     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
460     {
461     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
462     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
463 picher 233 normal_f_i.norme();
464     liste_normales2.ajouter(normal_f_i);
465 picher 230 //Remplissage de la liste des voisins du triangle i
466     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
467     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
468     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
469     for (int j = 0;j<nb_voisins1;j++)
470     {
471     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
472     liste_voisins.ajouter(mgtri_1);
473     }
474     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
475     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
476     for (int j = 0;j<nb_voisins2;j++)
477     {
478     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
479     liste_voisins.ajouter(mgtri_2);
480     }
481     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
482     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
483     for (int j = 0;j<nb_voisins3;j++)
484     {
485     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
486     liste_voisins.ajouter(mgtri_3);
487     }
488     liste_voisins.supprimer(mgtri_i);
489     int nb_voisins = liste_voisins.get_nb();
490     double w_ij = 1./nb_voisins;
491     double phi_i_min = 10.;
492     double s_i = 0.0;
493     double phi_im = 0.0;
494 picher 231 double *phi_ij = new double[nb_voisins];
495 picher 230 OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
496     OT_VECTEUR_3D eta_i(0.,0.,0.);
497     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
498     int j = 0;
499     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
500     {
501     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
502     //1-Calculer la normale moyenne pour chaque triangle
503     normal_f_i_mean = normal_f_i_mean + w_ij*normal_f_j;
504 francois 234 //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 � Nb_voisins
505 picher 230 double prod_scalaire = normal_f_i*normal_f_j;
506 picher 233 if (prod_scalaire > 1.)
507     {
508     prod_scalaire = 1.;
509     }
510     if (prod_scalaire < -1.)
511     {
512     prod_scalaire = -1.;
513     }
514 picher 230 phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
515     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
516     if (phi_ij[j] < phi_i_min)
517     {
518     phi_i_min = phi_ij[j];
519     eta_i = normal_f_j;
520     }
521     //3.1-On calcul l'angle moyen phi_im
522     phi_im = phi_im + w_ij*phi_ij[j];
523     j++;
524     }
525     normal_f_i_mean.norme();
526     j = 0;
527     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
528     {
529     //3.2-Calcul de s_i selon la variance
530     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
531     j++;
532     }
533 picher 231 delete[] phi_ij;
534 picher 230 //4-On calcule une nouvelle normale pour chaque triangle
535     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
536     normal_f_i_new.norme();
537     liste_normales.ajouter(normal_f_i_new);
538     liste_wij.ajouter(w_ij);
539     mgtri->change_nouveau_numero(k);
540     k++;
541     }
542    
543     LISTE_MG_NOEUD::iterator it_no;
544     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
545     {
546     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
547     double w_ij_prime = 0.0;
548     OT_VECTEUR_3D v_temp(0.,0.,0.);
549     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
550     for(int j=0;j<nb_voisins_j;j++)
551     {
552     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
553 francois 234 //On calcule le centro�de cj du triangle mgtri_j
554 picher 230 MG_NOEUD* n1 = mgtri_j->get_noeud1();
555     MG_NOEUD* n2 = mgtri_j->get_noeud2();
556     MG_NOEUD* n3 = mgtri_j->get_noeud3();
557     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
558     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
559     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
560     //On forme le vecteur vi_cj
561     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
562     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
563 francois 234 // w_ij_prime correspond � la somme des aires des triangles voisins du noeuds
564 picher 233 OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
565     OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
566     OT_VECTEUR_3D prodvect = AB&AC;
567     double w_ij = 0.5*prodvect.get_longueur();
568 picher 230 w_ij_prime = w_ij_prime + w_ij;
569     v_temp = v_temp + w_ij*(vi_cj*normal_f_i_new)*normal_f_i_new;
570     }
571 francois 234 //5-On met � jour la position des noeuds
572 picher 230 v_i = v_i + v_temp/w_ij_prime;
573 picher 233 int origine = noeud_i->get_origine();
574     if (origine != IMPOSE)
575     {
576     noeud_i->change_x(v_i.get_x());
577     noeud_i->change_y(v_i.get_y());
578     noeud_i->change_z(v_i.get_z());
579     }
580 picher 230 }
581    
582 francois 234 //Crit�re d'arr�t de l'algorithme
583 picher 230 int l=0;
584     int nb_tri = mg_mai->get_nb_mg_triangle();
585     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
586     {
587     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
588 picher 233 //OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
589     OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
590 picher 230 OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
591 picher 233 double critere = 1. - fabs(normal_f_i*normal_f_i_new);
592 picher 231 fprintf(out,"%lf\n",critere);
593 picher 230 if (critere <= epsilon) l++;
594 picher 233 if (critere >= 1.)
595     {
596     int no_id = mgtri->get_noeud1()->get_id();
597     cout << "no_id : " << no_id << endl;
598     char str[50];
599     sprintf (str, "test%d.magic", no_id);
600     //if (no_id == 18615) visualisation2(mg_mai,gest2,str,no_id);
601     }
602 picher 230 }
603 picher 233 if (nb_tri == l) fin = 1;
604 picher 230 compteur++;
605 francois 234 cout << "Fin de l'it�ration #" << compteur << endl;
606     fprintf(out,"\nFin de l'it�ration #%d\n\n",compteur);
607 picher 230 }
608     while ((fin == 0) && (compteur < iter_max));
609    
610     if (fin == 1)
611     {
612 francois 234 cout << "Convergence de la m�thode en " << compteur << " it�rations" << endl;
613 picher 230 }
614     else
615     {
616 francois 234 cout << "Arr�t de la proc�dure apr�s " << compteur << " it�rations" << endl;
617 picher 230 }
618    
619 picher 231 fclose(out);
620 picher 230 } //Fin de lissage
621    
622 picher 233 int MG_LISSAGE::extract_skin(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2,double frac_min, int *mai2_id)
623 francois 222 {
624    
625     //�tape 0 - On commence par mettre � z�ro tous les nouveau_numero des triangles et des noeuds du maillage
626     LISTE_MG_TRIANGLE::iterator it_tri;
627     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
628     {
629     mgtri->change_nouveau_numero(0);
630 francois 224 mgtri->change_origine(MAILLEUR_AUTO);
631 francois 222 }
632     LISTE_MG_NOEUD::iterator it_noeud;
633     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
634     {
635     mgnoeud->change_nouveau_numero(0);
636     }
637    
638     //�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
639 picher 233 LISTE_MG_TETRA::iterator it_tetra;
640     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
641 francois 222 {
642     int origine = mgtet->get_origine();
643 francois 224 if (origine==IMPOSE)
644 francois 222 {
645 francois 224 mgtet->get_triangle1()->change_origine(IMPOSE);
646     mgtet->get_triangle2()->change_origine(IMPOSE);
647     mgtet->get_triangle3()->change_origine(IMPOSE);
648     mgtet->get_triangle4()->change_origine(IMPOSE);
649     }
650     if (origine==OPTIMISE)
651     {
652     if (mgtet->get_triangle1()->get_origine()!=IMPOSE) mgtet->get_triangle1()->change_origine(OPTIMISE);
653     if (mgtet->get_triangle2()->get_origine()!=IMPOSE) mgtet->get_triangle2()->change_origine(OPTIMISE);
654     if (mgtet->get_triangle3()->get_origine()!=IMPOSE) mgtet->get_triangle3()->change_origine(OPTIMISE);
655     if (mgtet->get_triangle4()->get_origine()!=IMPOSE) mgtet->get_triangle4()->change_origine(OPTIMISE);
656     }
657    
658     if (((origine == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
659     ((origine == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
660     ((origine == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
661     ((origine == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
662     ((origine == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
663     ((origine == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
664    
665     {
666 francois 222 int num1 = mgtet->get_triangle1()->get_nouveau_numero();
667     int num2 = mgtet->get_triangle2()->get_nouveau_numero();
668     int num3 = mgtet->get_triangle3()->get_nouveau_numero();
669     int num4 = mgtet->get_triangle4()->get_nouveau_numero();
670     num1++;
671     num2++;
672     num3++;
673     num4++;
674     mgtet->get_triangle1()->change_nouveau_numero(num1);
675     mgtet->get_triangle2()->change_nouveau_numero(num2);
676     mgtet->get_triangle3()->change_nouveau_numero(num3);
677     mgtet->get_triangle4()->change_nouveau_numero(num4);
678     }
679     }
680    
681     //�tape 2 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour identifier les noeuds leur appartenant
682     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
683     {
684     int num = mgtri->get_nouveau_numero();
685     if (num == 1)
686     {
687     MG_NOEUD* noeud1 = mgtri->get_noeud1();
688     MG_NOEUD* noeud2 = mgtri->get_noeud2();
689     MG_NOEUD* noeud3 = mgtri->get_noeud3();
690     noeud1->change_nouveau_numero(1);
691     noeud2->change_nouveau_numero(1);
692     noeud3->change_nouveau_numero(1);
693     }
694     }
695    
696 francois 234 //�tape 3 - On cr�e un nouveau maillage pour la peau
697 francois 222 MG_MAILLAGE* mai2 = new MG_MAILLAGE(NULL);
698     gest2.ajouter_mg_maillage(mai2);
699    
700     // //�tape 4 - On boucle l'ensemble des noeuds identifi�s � l'�tape 2 pour les recr�er dans le second maillage
701     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
702     {
703     int num = mgnoeud->get_nouveau_numero();
704     if (num == 1)
705     {
706     double x = mgnoeud->get_x();
707     double y = mgnoeud->get_y();
708     double z = mgnoeud->get_z();
709 francois 224 int origine=TRIANGULATION;
710     if (mgnoeud->get_origine()==IMPOSE) origine=IMPOSE;
711     MG_NOEUD* noeud1 = new MG_NOEUD(NULL,x,y,z,origine);
712 francois 222 mai2->ajouter_mg_noeud(noeud1);
713     mgnoeud->change_nouveau_numero(noeud1->get_id());
714 francois 234 noeud1->change_nouveau_numero(mgnoeud->get_id());
715 francois 222 }
716     }
717    
718     // //�tape 5 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour les recr�er dans le maillage 2
719 picher 233 for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
720 francois 222 {
721 francois 224 int originetet=mgtet->get_origine();
722     if (((originetet == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
723     ((originetet == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
724     ((originetet == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
725     ((originetet == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
726     ((originetet == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
727     ((originetet == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
728 francois 222 {
729 francois 224 MG_NOEUD* noeud1 = mgtet->get_noeud1();
730     MG_NOEUD* noeud2 = mgtet->get_noeud2();
731     MG_NOEUD* noeud3 = mgtet->get_noeud3();
732     MG_NOEUD* noeud4 = mgtet->get_noeud4();
733 francois 222 MG_NOEUD* node1 = mai2->get_mg_noeudid(noeud1->get_nouveau_numero());
734     MG_NOEUD* node2 = mai2->get_mg_noeudid(noeud2->get_nouveau_numero());
735     MG_NOEUD* node3 = mai2->get_mg_noeudid(noeud3->get_nouveau_numero());
736 francois 224 MG_NOEUD* node4 = mai2->get_mg_noeudid(noeud4->get_nouveau_numero());
737     MG_TRIANGLE* tri1=mgtet->get_triangle1();
738     MG_TRIANGLE* tri2=mgtet->get_triangle2();
739     MG_TRIANGLE* tri3=mgtet->get_triangle3();
740     MG_TRIANGLE* tri4=mgtet->get_triangle4();
741     if (tri1->get_nouveau_numero()==1)
742     {
743     int origine=TRIANGULATION;
744     if (tri1->get_origine()==IMPOSE) origine=IMPOSE;
745     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node2,node3,mai2,origine);
746     tripeau->change_nouveau_numero(0);
747     tri1->change_nouveau_numero(0);
748     }
749     if (tri2->get_nouveau_numero()==1)
750     {
751     int origine=TRIANGULATION;
752     if (tri2->get_origine()==IMPOSE) origine=IMPOSE;
753     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node4,node2,mai2,origine);
754     tripeau->change_nouveau_numero(0);
755     tri2->change_nouveau_numero(0);
756     }
757     if (tri3->get_nouveau_numero()==1)
758     {
759     int origine=TRIANGULATION;
760     if (tri3->get_origine()==IMPOSE) origine=IMPOSE;
761     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node2,node4,node3,mai2,origine);
762     tripeau->change_nouveau_numero(0);
763     tri3->change_nouveau_numero(0);
764     }
765     if (tri4->get_nouveau_numero()==1)
766     {
767     int origine=TRIANGULATION;
768     if (tri4->get_origine()==IMPOSE) origine=IMPOSE;
769     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node3,node4,mai2,origine);
770     tripeau->change_nouveau_numero(0);
771     tri4->change_nouveau_numero(0);
772     }
773 francois 222 }
774 francois 224 }
775     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
776     // etape 6 recherche des vosins
777     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
778     {
779     MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
780     MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
781     MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
782     MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
783     tripeau->change_voisin1(voisin1);
784     tripeau->change_voisin2(voisin2);
785     tripeau->change_voisin3(voisin3);
786     }
787     // //�tape 7 - On recherche les peaux
788     int fin;
789     do
790     {
791     fin=1;
792     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
793     {
794     MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
795     if (tripeau->get_nouveau_numero()==0)
796     {
797     fin=0;
798     std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
799     lst_peau.push_back(peau);
800     tripeau->change_nouveau_numero(1);
801     peau->push_back(tripeau);
802     determine_peau(peau);
803     }
804     }
805     }
806     while (fin==0);
807 picher 233
808 francois 224 std::cout << lst_peau.size() << " peau" << endl;
809     for (int i=0;i<lst_peau.size();i++)
810     std::cout << " peau " << i << " " << lst_peau[i]->size() << " triangles" << endl;
811     //test de manifold
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 francois 234 TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeuddepuisarete;
818 francois 224 int nbpeau=lst_peau.size();
819     for (int i=0;i<nbpeau;i++)
820     {
821     int nbtri=lst_peau[i]->size();
822     for (int j=0;j<nbtri;j++)
823     {
824     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
825     tripeau->get_segment1()->change_nouveau_numero(tripeau->get_segment1()->get_nouveau_numero()+1);
826     tripeau->get_segment2()->change_nouveau_numero(tripeau->get_segment2()->get_nouveau_numero()+1);
827     tripeau->get_segment3()->change_nouveau_numero(tripeau->get_segment3()->get_nouveau_numero()+1);
828     if (tripeau->get_segment1()->get_nouveau_numero()>2)
829     nonmanifoldarete.ajouter(tripeau->get_segment1());
830     if (tripeau->get_segment2()->get_nouveau_numero()>2)
831     nonmanifoldarete.ajouter(tripeau->get_segment2());
832     if (tripeau->get_segment3()->get_nouveau_numero()>2)
833     nonmanifoldarete.ajouter(tripeau->get_segment3());
834     }
835     }
836 picher 233
837 francois 224 int nbnonmanifoldarete=nonmanifoldarete.get_nb();
838     for (int i=0;i<nbnonmanifoldarete;i++)
839     {
840     MG_SEGMENT* seg=nonmanifoldarete.get(i);
841     MG_NOEUD* n1=seg->get_noeud1();
842     MG_NOEUD* n2=seg->get_noeud2();
843 francois 234 nonmanifoldnoeuddepuisarete.ajouter(n1);
844     nonmanifoldnoeuddepuisarete.ajouter(n2);
845 francois 224 }
846 francois 234 LISTE_MG_NOEUD::iterator itnoeud;
847     for (MG_NOEUD* no=mai2->get_premier_noeud(itnoeud);no!=NULL;no=mai2->get_suivant_noeud(itnoeud))
848     {
849     if (nonmanifoldnoeuddepuisarete.existe(no)) continue;
850     if (est_non_manifold(no)) nonmanifoldnoeud.ajouter(no);
851     }
852 francois 224
853     cout << "Non manifold par arete " << nonmanifoldarete.get_nb() << endl;
854 picher 233 cout << "Non manifold par noeud " << nonmanifoldnoeud.get_nb() << endl;
855 picher 230
856 francois 234 //�tape 8 - Traitement des cas de non manifold
857     //non manifold par ar�te
858 picher 233 for (int i=0;i<nbnonmanifoldarete;i++)
859     {
860     MG_SEGMENT* segment=nonmanifoldarete.get(i);
861 francois 234 MG_NOEUD* noeud1 = mg_mai->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
862     MG_NOEUD* noeud2 = mg_mai->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
863     int nb_tetra = noeud1->get_lien_tetra()->get_nb();
864     for (int j=0;j<nb_tetra;j++)
865 picher 233 {
866 francois 234 MG_TETRA* mgtet =noeud1->get_lien_tetra()->get(j);
867     int originetet=mgtet->get_origine();
868     if (originetet == MAILLEUR_AUTO)
869 picher 233 {
870 francois 234 //On r�active le tetra si l'autre noeud du segment lui appartient aussi
871     MG_NOEUD* no1 = mgtet->get_noeud1();
872     MG_NOEUD* no2 = mgtet->get_noeud2();
873     MG_NOEUD* no3 = mgtet->get_noeud3();
874     MG_NOEUD* no4 = mgtet->get_noeud4();
875     if((no1 == noeud2) ||(no2 == noeud2) || (no3 == noeud2) || (no4 == noeud2))
876     mgtet->change_origine(OPTIMISE);
877 picher 233 }
878     }
879     }
880    
881     //non manifold par noeud
882 picher 231 int nbnonmanifoldnoeud=nonmanifoldnoeud.get_nb();
883 picher 233 //if (nbnonmanifoldnoeud > 1) visualisation2(mai2,gest2,"nm_noeud.magic",39199);
884 picher 231 for (int i=0;i<nbnonmanifoldnoeud;i++)
885     {
886 francois 234 MG_NOEUD* noeud=mg_mai->get_mg_noeudid(nonmanifoldnoeud.get(i)->get_nouveau_numero());
887     int nb_tetra = noeud->get_lien_tetra()->get_nb();
888     for (int j = 0;j<nb_tetra;j++)
889 picher 231 {
890 francois 234 MG_TETRA* mgtet =noeud->get_lien_tetra()->get(j);
891     int originetet=mgtet->get_origine();
892     if (originetet == MAILLEUR_AUTO)
893 picher 231 {
894 francois 234 mgtet->change_origine(OPTIMISE);
895 picher 231 }
896     }
897     }
898 francois 234 *mai2_id = mai2->get_id();
899     if ((nbnonmanifoldarete != 0) || (nbnonmanifoldnoeud != 0))
900     {
901     for (int i=0;i<lst_peau.size();i++)
902     {
903     delete lst_peau[i];
904     }
905     lst_peau.clear();
906     gest2.supprimer_mg_maillageid(*mai2_id);
907     return 0;
908     }
909 picher 233 //if (nbnonmanifoldnoeud > 1) visualisation2(mai2,gest2,"nm_noeud_resolu.magic",39199);
910     // for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
911     // {
912     // int id = mgnoeud->get_id();
913     // if (id == 37985)
914     // {
915     // int id1 = mgnoeud->get_nouveau_numero();
916     // cout << "nouv_id : " << id1 << endl;
917     // }
918     // }
919 picher 230
920 picher 233 //if ((nbnonmanifoldnoeud == 0) && (nbnonmanifoldarete == 0)) visualisation2(mai2,gest2,"nm_noeud_resolu.magic",39199);
921 picher 231
922 francois 234 //�tape 9 - Suppression des peaux de tr�s faible taille (mati�re flottante)
923 picher 233 cout << "Suppression des peaux de faible taille" << endl;
924     int nbtri_total = mai2->get_nb_mg_triangle();
925 picher 230 for (int i=0;i<nbpeau;i++)
926     {
927     int nbtri=lst_peau[i]->size();
928 picher 233 if (nbtri < frac_min*nbtri_total)
929 picher 230 {
930     for (int j=0;j<nbtri;j++)
931     {
932     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
933 picher 233 mai2->supprimer_mg_triangleid(tripeau->get_id());
934 picher 230 }
935     }
936 picher 233 }
937     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
938    
939 francois 234 return 1;
940 francois 222 }
941    
942 francois 224 void MG_LISSAGE::determine_peau(std::vector<MG_TRIANGLE_PEAU*> * peau)
943 francois 222 {
944 francois 224 int num=0;
945     while (num!=peau->size())
946     {
947     MG_TRIANGLE_PEAU* p=(*peau)[num];
948     if (p->get_voisin1()->get_nouveau_numero()==0)
949     {
950     peau->push_back(p->get_voisin1());
951     p->get_voisin1()->change_nouveau_numero(1);
952     }
953     if (p->get_voisin2()->get_nouveau_numero()==0)
954     {
955     peau->push_back(p->get_voisin2());
956     p->get_voisin2()->change_nouveau_numero(1);
957     }
958     if (p->get_voisin3()->get_nouveau_numero()==0)
959     {
960     peau->push_back(p->get_voisin3());
961     p->get_voisin3()->change_nouveau_numero(1);
962     }
963     num++;
964     }
965     }
966    
967    
968     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)
969     {
970 francois 222 MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
971     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
972     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
973     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
974     if (mgsegment1==NULL)
975 francois 224 mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,origine);
976 francois 222 if (mgsegment2==NULL)
977 francois 224 mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,origine);
978 francois 222 if (mgsegment3==NULL)
979 francois 224 mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,origine);
980     MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
981 francois 222 mg_maillage->ajouter_mg_triangle(mtriangle);
982     return mtriangle;
983     }
984    
985 francois 224 MG_TRIANGLE_PEAU* MG_LISSAGE::recherche_voisin(MG_NOEUD* mg_noeud1,MG_NOEUD* mg_noeud2,MG_TRIANGLE_PEAU* triref)
986 francois 222 {
987 francois 224 MG_TRIANGLE_PEAU* trisol=NULL;
988     double angleref=4.*M_PI;
989 francois 222 int nb1=mg_noeud1->get_lien_triangle()->get_nb();
990     int nb2=mg_noeud2->get_lien_triangle()->get_nb();
991     for (int i=0;i<nb1;i++)
992     for (int j=0;j<nb2;j++)
993 francois 224 if (mg_noeud1->get_lien_triangle()->get(i)==mg_noeud2->get_lien_triangle()->get(j))
994     if ((MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i)!=triref)
995     {
996     double angle=calcul_angle(triref,(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i));
997     if (angle<angleref)
998     {
999     angleref=angle;
1000     trisol=(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i);
1001     }
1002     }
1003     return trisol;
1004 francois 222 }
1005 francois 223
1006 francois 234 int MG_LISSAGE::est_non_manifold(MG_NOEUD* no)
1007     {
1008     static int compteur=0;
1009     compteur++;
1010     int nb_tri=no->get_lien_triangle()->get_nb();
1011     TPL_MAP_ENTITE<MG_TRIANGLE*> lst;
1012     for (int i=0;i<nb_tri;i++)
1013     lst.ajouter(no->get_lien_triangle()->get(i));
1014     MG_TRIANGLE_PEAU* premier_triangle=(MG_TRIANGLE_PEAU*)lst.get(0);
1015     lst.supprimer(premier_triangle);
1016     MG_TRIANGLE_PEAU* tricour=(MG_TRIANGLE_PEAU*)premier_triangle;
1017     int i=0;
1018     do
1019     {
1020     if (lst.existe(tricour->get_voisin1()) || ((tricour->get_voisin1()==premier_triangle)&& (i>1)))
1021     {
1022     tricour=tricour->get_voisin1();
1023     lst.supprimer(tricour);
1024     }
1025     else if (lst.existe(tricour->get_voisin2()) || ((tricour->get_voisin2()==premier_triangle)&& (i>1)))
1026     {
1027     tricour=tricour->get_voisin2();
1028     lst.supprimer(tricour);
1029     }
1030     else if (lst.existe(tricour->get_voisin3()) || ((tricour->get_voisin3()==premier_triangle)&& (i>1)))
1031     {
1032     tricour=tricour->get_voisin3();
1033     lst.supprimer(tricour);
1034     }
1035     i++;
1036     }
1037     while (tricour!=premier_triangle);
1038     if (lst.get_nb()>0)
1039     return 1;
1040     return 0;
1041     }
1042 francois 224
1043     double MG_LISSAGE::calcul_angle(MG_TRIANGLE_PEAU* ft1,MG_TRIANGLE_PEAU* ft2)
1044 francois 223 {
1045 francois 224 MG_NOEUD* noeud1=ft1->get_noeud1();
1046     MG_NOEUD* noeud2=ft1->get_noeud2();
1047     MG_NOEUD* noeud3=ft1->get_noeud3();
1048     MG_NOEUD* noeud4=ft2->get_noeud1();
1049     MG_NOEUD* noeud5=ft2->get_noeud2();
1050     MG_NOEUD* noeud6=ft2->get_noeud3();
1051    
1052     MG_NOEUD* noeuda=NULL;
1053     MG_NOEUD* noeudb=NULL;
1054     MG_NOEUD* noeudc=NULL;
1055     MG_NOEUD* noeudd=NULL;
1056    
1057    
1058     if (!(noeud1==noeud4))
1059     if (!(noeud1==noeud5))
1060     if (!(noeud1==noeud6))
1061     {
1062     noeuda=noeud1;
1063     noeudb=noeud2;
1064     noeudc=noeud3;
1065     }
1066     if (!(noeud2==noeud4))
1067     if (!(noeud2==noeud5))
1068     if (!(noeud2==noeud6))
1069     {
1070     noeuda=noeud2;
1071     noeudb=noeud1;
1072     noeudc=noeud3;
1073     }
1074     if (!(noeud3==noeud4))
1075     if (!(noeud3==noeud5))
1076     if (!(noeud3==noeud6))
1077     {
1078     noeuda=noeud3;
1079     noeudb=noeud1;
1080     noeudc=noeud2;
1081     }
1082     if (noeuda==NULL) return 2*M_PI;
1083     if (!(noeud4==noeudb))
1084     if (!(noeud4==noeudc)) noeudd=noeud4;
1085     if (!(noeud5==noeudb))
1086     if (!(noeud5==noeudc)) noeudd=noeud5;
1087     if (!(noeud6==noeudb))
1088     if (!(noeud6==noeudc)) noeudd=noeud6;
1089    
1090     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()));
1091     OT_VECTEUR_3D bc(noeudc->get_x()-noeudb->get_x(),noeudc->get_y()-noeudb->get_y(),noeudc->get_z()-noeudb->get_z());
1092     bc.norme();
1093     double xyz[3];
1094     xyz[0]=noeuda->get_x();
1095     xyz[1]=noeuda->get_y();
1096     xyz[2]=noeuda->get_z();
1097     OT_VECTEUR_3D n1(xyz,milieu.get_xyz());
1098     n1.norme();
1099     OT_VECTEUR_3D n=bc&n1;
1100     n1=n&bc;
1101     n1.norme();
1102     xyz[0]=noeudd->get_x();
1103     xyz[1]=noeudd->get_y();
1104     xyz[2]=noeudd->get_z();
1105     OT_VECTEUR_3D n2(xyz,milieu.get_xyz());
1106     n2.norme();
1107     n=bc&n2;
1108     n2=n&bc;
1109     n2.norme();
1110    
1111     double ps=n1*n2;
1112     if (ps>1.) ps=1.;
1113     if (ps<-1.) ps=-1.;
1114     double angle=acos(ps);
1115     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
1116     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
1117     OT_VECTEUR_3D nbnd(noeudd->get_x()-noeudb->get_x(),noeudd->get_y()-noeudb->get_y(),noeudd->get_z()-noeudb->get_z());
1118     n=n1n2&n1n3;
1119     n.norme();
1120     nbnd.norme();
1121    
1122     ps=n*nbnd;
1123     if (ps<-0.0001) angle=2*M_PI-angle;
1124     return angle;
1125    
1126    
1127 francois 223 }
1128 picher 230
1129     double MG_LISSAGE::ponderation_gaussian(double s,double sigma)
1130     {
1131     double w_s;
1132     w_s = exp(-(s*s)/(2.*sigma*sigma));
1133     return w_s;
1134     }
1135    
1136     double MG_LISSAGE::ponderation_laplacian(double s,double sigma)
1137     {
1138     double w_s;
1139     w_s = exp(-(sqrt(2.)*fabs(s))/sigma);
1140     return w_s;
1141     }
1142    
1143     double MG_LISSAGE::ponderation_elfallahford(double s,double sigma)
1144     {
1145     double w_s;
1146     w_s = 1./sqrt(1+pow((s/sigma),2));
1147     return w_s;
1148     }