ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_posttraitement.cpp
Revision: 230
Committed: Fri Dec 18 19:24:29 2009 UTC (15 years, 4 months ago) by picher
Original Path: magic/lib/optimisation/optimisation/src/mg_lissage.cpp
File size: 23987 byte(s)
Log Message:
Ajout d'une methode de lissage a mg_lissage et ajout de l'operation de lissage a mg_operation

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 222 MG_LISSAGE::MG_LISSAGE()
18     {
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    
29 francois 224 void MG_LISSAGE::conserve(int origine)
30     {
31     etat[(origine-1000)/10]=1;
32     }
33    
34 picher 230 void MG_LISSAGE::bruitage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
35     {
36     //Méthode test pour bruiter un maillage initial et tester la méthode de lissage
37     srand(time(NULL));
38     LISTE_MG_NOEUD::iterator it_no;
39     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
40     {
41     //On modifie la position des noeuds légèrement (ajout d'un bruit au maillage initial)
42     int signe;
43     int delta;
44     double delta_x,delta_y,delta_z;
45     signe = rand()%2 + 1;
46     delta = rand()%2 + 1;
47     if(signe == 1)
48     {
49     delta_x = delta*0.5;
50     }
51     else //(signe == 2)
52     {
53     delta_x = -delta*0.5;
54     }
55     signe = rand()%2 + 1;
56     delta = rand()%2 + 1;
57     if(signe == 1)
58     {
59     delta_y = delta*0.5;
60     }
61     else //(signe == 2)
62     {
63     delta_y = -delta*0.5;
64     }
65     signe = rand()%2 + 1;
66     delta = rand()%2 + 1;
67     if(signe == 1)
68     {
69     delta_z = delta*0.5;
70     }
71     else //(signe == 2)
72     {
73     delta_z = -delta*0.5;
74     }
75     double x_new = noeud_i->get_x() + delta_x;
76     double y_new = noeud_i->get_y() + delta_y;
77     double z_new = noeud_i->get_z() + delta_z;
78     noeud_i->change_x(x_new);
79     noeud_i->change_y(y_new);
80     noeud_i->change_z(z_new);
81     }
82     }
83 francois 224
84 picher 230 void MG_LISSAGE::lissage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max)
85     {
86     //Méthode de lissage inspirée de Chen(2005)
87     double un_sur_pi = 1./M_PI;
88     int compteur = 0;
89     int fin = 0;
90 francois 224
91 picher 230 do
92     {
93     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
94     TPL_LISTE_ENTITE<double> liste_wij;
95     LISTE_MG_TRIANGLE::iterator it_tri;
96     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
97     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
98     {
99     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
100     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
101     //Remplissage de la liste des voisins du triangle i
102     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
103     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
104     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
105     for (int j = 0;j<nb_voisins1;j++)
106     {
107     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
108     liste_voisins.ajouter(mgtri_1);
109     }
110     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
111     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
112     for (int j = 0;j<nb_voisins2;j++)
113     {
114     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
115     liste_voisins.ajouter(mgtri_2);
116     }
117     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
118     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
119     for (int j = 0;j<nb_voisins3;j++)
120     {
121     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
122     liste_voisins.ajouter(mgtri_3);
123     }
124     liste_voisins.supprimer(mgtri_i);
125     int nb_voisins = liste_voisins.get_nb();
126     double w_ij = 1./nb_voisins;
127     double phi_i_min = 10.;
128     double s_i = 0.0;
129     double phi_im = 0.0;
130     double *phi_ij = (double*) malloc (sizeof(double)*(nb_voisins));
131     OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
132     OT_VECTEUR_3D eta_i(0.,0.,0.);
133     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
134     int j = 0;
135     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
136     {
137     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
138     //1-Calculer la normale moyenne pour chaque triangle
139     normal_f_i_mean = normal_f_i_mean + w_ij*normal_f_j;
140     //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 à Nb_voisins
141     double prod_scalaire = normal_f_i*normal_f_j;
142     if (fabs(prod_scalaire) > 1.) prod_scalaire = 1.;
143     phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
144     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
145     if (phi_ij[j] < phi_i_min)
146     {
147     phi_i_min = phi_ij[j];
148     eta_i = normal_f_j;
149     }
150     //3.1-On calcul l'angle moyen phi_im
151     phi_im = phi_im + w_ij*phi_ij[j];
152     j++;
153     }
154     normal_f_i_mean.norme();
155     j = 0;
156     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
157     {
158     //3.2-Calcul de s_i selon la variance
159     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
160     j++;
161     }
162     //4-On calcule une nouvelle normale pour chaque triangle
163     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
164     normal_f_i_new.norme();
165     liste_normales.ajouter(normal_f_i_new);
166     liste_wij.ajouter(w_ij);
167     mgtri->change_nouveau_numero(k);
168     k++;
169     }
170    
171     LISTE_MG_NOEUD::iterator it_no;
172     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
173     {
174     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
175     double w_ij_prime = 0.0;
176     OT_VECTEUR_3D v_temp(0.,0.,0.);
177     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
178     for(int j=0;j<nb_voisins_j;j++)
179     {
180     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
181     //On calcule le centroïde cj du triangle mgtri_j
182     MG_NOEUD* n1 = mgtri_j->get_noeud1();
183     MG_NOEUD* n2 = mgtri_j->get_noeud2();
184     MG_NOEUD* n3 = mgtri_j->get_noeud3();
185     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
186     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
187     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
188     //On forme le vecteur vi_cj
189     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
190     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
191     double w_ij = liste_wij.get(mgtri_j->get_nouveau_numero());
192     w_ij_prime = w_ij_prime + w_ij;
193     v_temp = v_temp + w_ij*(vi_cj*normal_f_i_new)*normal_f_i_new;
194     }
195     //5-On met à jour la position des noeuds
196     v_i = v_i + v_temp/w_ij_prime;
197     noeud_i->change_x(v_i.get_x());
198     noeud_i->change_y(v_i.get_y());
199     noeud_i->change_z(v_i.get_z());
200     }
201    
202     //Critère d'arrêt de l'algorithme
203     int l=0;
204     int nb_tri = mg_mai->get_nb_mg_triangle();
205     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
206     {
207     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
208     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
209     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
210     double critere = 1. - normal_f_i*normal_f_i_new;
211     if (critere <= epsilon) l++;
212     }
213     if (l == nb_tri) fin = 1;
214     compteur++;
215     cout << "Fin de l'itération #" << compteur << endl;
216     }
217     while ((fin == 0) && (compteur < iter_max));
218    
219     if (fin == 1)
220     {
221     cout << "Convergence de la méthode en " << compteur << " itérations" << endl;
222     }
223     else
224     {
225     cout << "Arrêt de la procédure après " << compteur << " itérations" << endl;
226     }
227    
228     } //Fin de lissage
229    
230 francois 222 void MG_LISSAGE::extract_skin(FEM_MAILLAGE* mai,MG_GESTIONNAIRE& gest2)
231     {
232    
233     //�tape 0 - On commence par mettre � z�ro tous les nouveau_numero des triangles et des noeuds du maillage
234     MG_MAILLAGE* mg_mai = (MG_MAILLAGE*)mai->get_mg_maillage();
235     LISTE_MG_TRIANGLE::iterator it_tri;
236     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
237     {
238     mgtri->change_nouveau_numero(0);
239 francois 224 mgtri->change_origine(MAILLEUR_AUTO);
240 francois 222 }
241     LISTE_MG_NOEUD::iterator it_noeud;
242     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
243     {
244     mgnoeud->change_nouveau_numero(0);
245     }
246    
247     //�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
248     LISTE_FEM_TETRA::iterator it_tetra;
249     for (FEM_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
250     {
251     MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
252     int origine = mgtet->get_origine();
253 francois 224 if (origine==IMPOSE)
254 francois 222 {
255 francois 224 mgtet->get_triangle1()->change_origine(IMPOSE);
256     mgtet->get_triangle2()->change_origine(IMPOSE);
257     mgtet->get_triangle3()->change_origine(IMPOSE);
258     mgtet->get_triangle4()->change_origine(IMPOSE);
259     }
260     if (origine==OPTIMISE)
261     {
262     if (mgtet->get_triangle1()->get_origine()!=IMPOSE) mgtet->get_triangle1()->change_origine(OPTIMISE);
263     if (mgtet->get_triangle2()->get_origine()!=IMPOSE) mgtet->get_triangle2()->change_origine(OPTIMISE);
264     if (mgtet->get_triangle3()->get_origine()!=IMPOSE) mgtet->get_triangle3()->change_origine(OPTIMISE);
265     if (mgtet->get_triangle4()->get_origine()!=IMPOSE) mgtet->get_triangle4()->change_origine(OPTIMISE);
266     }
267    
268     if (((origine == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
269     ((origine == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
270     ((origine == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
271     ((origine == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
272     ((origine == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
273     ((origine == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
274    
275     {
276 francois 222 int num1 = mgtet->get_triangle1()->get_nouveau_numero();
277     int num2 = mgtet->get_triangle2()->get_nouveau_numero();
278     int num3 = mgtet->get_triangle3()->get_nouveau_numero();
279     int num4 = mgtet->get_triangle4()->get_nouveau_numero();
280     num1++;
281     num2++;
282     num3++;
283     num4++;
284     mgtet->get_triangle1()->change_nouveau_numero(num1);
285     mgtet->get_triangle2()->change_nouveau_numero(num2);
286     mgtet->get_triangle3()->change_nouveau_numero(num3);
287     mgtet->get_triangle4()->change_nouveau_numero(num4);
288     }
289     }
290    
291     //�tape 2 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour identifier les noeuds leur appartenant
292     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
293     {
294     int num = mgtri->get_nouveau_numero();
295     if (num == 1)
296     {
297     MG_NOEUD* noeud1 = mgtri->get_noeud1();
298     MG_NOEUD* noeud2 = mgtri->get_noeud2();
299     MG_NOEUD* noeud3 = mgtri->get_noeud3();
300     noeud1->change_nouveau_numero(1);
301     noeud2->change_nouveau_numero(1);
302     noeud3->change_nouveau_numero(1);
303     }
304     }
305    
306     //�tape 3 - On cr�e un nouveau gestionnaire et un nouveau maillage
307    
308    
309     MG_MAILLAGE* mai2 = new MG_MAILLAGE(NULL);
310     gest2.ajouter_mg_maillage(mai2);
311    
312     // //�tape 4 - On boucle l'ensemble des noeuds identifi�s � l'�tape 2 pour les recr�er dans le second maillage
313     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
314     {
315     int num = mgnoeud->get_nouveau_numero();
316     if (num == 1)
317     {
318     double x = mgnoeud->get_x();
319     double y = mgnoeud->get_y();
320     double z = mgnoeud->get_z();
321 francois 224 int origine=TRIANGULATION;
322     if (mgnoeud->get_origine()==IMPOSE) origine=IMPOSE;
323     MG_NOEUD* noeud1 = new MG_NOEUD(NULL,x,y,z,origine);
324 francois 222 mai2->ajouter_mg_noeud(noeud1);
325     mgnoeud->change_nouveau_numero(noeud1->get_id());
326     }
327     }
328    
329     // //�tape 5 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour les recr�er dans le maillage 2
330 francois 224 for (FEM_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
331 francois 222 {
332 francois 224 MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
333     int originetet=mgtet->get_origine();
334     if (((originetet == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
335     ((originetet == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
336     ((originetet == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
337     ((originetet == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
338     ((originetet == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
339     ((originetet == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
340 francois 222 {
341 francois 224 MG_NOEUD* noeud1 = mgtet->get_noeud1();
342     MG_NOEUD* noeud2 = mgtet->get_noeud2();
343     MG_NOEUD* noeud3 = mgtet->get_noeud3();
344     MG_NOEUD* noeud4 = mgtet->get_noeud4();
345 francois 222 MG_NOEUD* node1 = mai2->get_mg_noeudid(noeud1->get_nouveau_numero());
346     MG_NOEUD* node2 = mai2->get_mg_noeudid(noeud2->get_nouveau_numero());
347     MG_NOEUD* node3 = mai2->get_mg_noeudid(noeud3->get_nouveau_numero());
348 francois 224 MG_NOEUD* node4 = mai2->get_mg_noeudid(noeud4->get_nouveau_numero());
349     MG_TRIANGLE* tri1=mgtet->get_triangle1();
350     MG_TRIANGLE* tri2=mgtet->get_triangle2();
351     MG_TRIANGLE* tri3=mgtet->get_triangle3();
352     MG_TRIANGLE* tri4=mgtet->get_triangle4();
353     if (tri1->get_nouveau_numero()==1)
354     {
355     int origine=TRIANGULATION;
356     if (tri1->get_origine()==IMPOSE) origine=IMPOSE;
357     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node2,node3,mai2,origine);
358     tripeau->change_nouveau_numero(0);
359     tri1->change_nouveau_numero(0);
360     }
361     if (tri2->get_nouveau_numero()==1)
362     {
363     int origine=TRIANGULATION;
364     if (tri2->get_origine()==IMPOSE) origine=IMPOSE;
365     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node4,node2,mai2,origine);
366     tripeau->change_nouveau_numero(0);
367     tri2->change_nouveau_numero(0);
368     }
369     if (tri3->get_nouveau_numero()==1)
370     {
371     int origine=TRIANGULATION;
372     if (tri3->get_origine()==IMPOSE) origine=IMPOSE;
373     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node2,node4,node3,mai2,origine);
374     tripeau->change_nouveau_numero(0);
375     tri3->change_nouveau_numero(0);
376     }
377     if (tri4->get_nouveau_numero()==1)
378     {
379     int origine=TRIANGULATION;
380     if (tri4->get_origine()==IMPOSE) origine=IMPOSE;
381     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node3,node4,mai2,origine);
382     tripeau->change_nouveau_numero(0);
383     tri4->change_nouveau_numero(0);
384     }
385 francois 222 }
386 francois 224 }
387     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
388     // etape 6 recherche des vosins
389     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
390     {
391     MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
392     MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
393     MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
394     MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
395     tripeau->change_voisin1(voisin1);
396     tripeau->change_voisin2(voisin2);
397     tripeau->change_voisin3(voisin3);
398     }
399     // //�tape 7 - On recherche les peaux
400     int fin;
401     do
402     {
403     fin=1;
404     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
405     {
406     MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
407     if (tripeau->get_nouveau_numero()==0)
408     {
409     fin=0;
410     std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
411     lst_peau.push_back(peau);
412     tripeau->change_nouveau_numero(1);
413     peau->push_back(tripeau);
414     determine_peau(peau);
415     }
416     }
417     }
418     while (fin==0);
419     std::cout << lst_peau.size() << " peau" << endl;
420     for (int i=0;i<lst_peau.size();i++)
421     std::cout << " peau " << i << " " << lst_peau[i]->size() << " triangles" << endl;
422     //test de manifold
423     for (MG_NOEUD* mgnoeud=mai2->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mai2->get_suivant_noeud(it_noeud))
424     mgnoeud->change_nouveau_numero(0);
425     LISTE_MG_SEGMENT::iterator itseg;
426     for (MG_SEGMENT* seg=mai2->get_premier_segment(itseg);seg!=NULL;seg=mai2->get_suivant_segment(itseg))
427     seg->change_nouveau_numero(0);
428     TPL_MAP_ENTITE<MG_SEGMENT*> nonmanifoldarete;
429     TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeud;
430     int nbpeau=lst_peau.size();
431     for (int i=0;i<nbpeau;i++)
432     {
433     int nbtri=lst_peau[i]->size();
434     for (int j=0;j<nbtri;j++)
435     {
436     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
437     tripeau->get_segment1()->change_nouveau_numero(tripeau->get_segment1()->get_nouveau_numero()+1);
438     tripeau->get_segment2()->change_nouveau_numero(tripeau->get_segment2()->get_nouveau_numero()+1);
439     tripeau->get_segment3()->change_nouveau_numero(tripeau->get_segment3()->get_nouveau_numero()+1);
440     if (tripeau->get_segment1()->get_nouveau_numero()>2)
441     nonmanifoldarete.ajouter(tripeau->get_segment1());
442     if (tripeau->get_segment2()->get_nouveau_numero()>2)
443     nonmanifoldarete.ajouter(tripeau->get_segment2());
444     if (tripeau->get_segment3()->get_nouveau_numero()>2)
445     nonmanifoldarete.ajouter(tripeau->get_segment3());
446     if ((tripeau->get_noeud1()->get_nouveau_numero()!=0) && (tripeau->get_noeud1()->get_nouveau_numero()!=i+1))
447     nonmanifoldnoeud.ajouter(tripeau->get_noeud1());
448     if ((tripeau->get_noeud2()->get_nouveau_numero()!=0) && (tripeau->get_noeud2()->get_nouveau_numero()!=i+1))
449     nonmanifoldnoeud.ajouter(tripeau->get_noeud2());
450     if ((tripeau->get_noeud3()->get_nouveau_numero()!=0) && (tripeau->get_noeud3()->get_nouveau_numero()!=i+1))
451     nonmanifoldnoeud.ajouter(tripeau->get_noeud3());
452     tripeau->get_noeud1()->change_nouveau_numero(i+1);
453     tripeau->get_noeud2()->change_nouveau_numero(i+1);
454     tripeau->get_noeud3()->change_nouveau_numero(i+1);
455     }
456     }
457     int nbnonmanifoldarete=nonmanifoldarete.get_nb();
458     for (int i=0;i<nbnonmanifoldarete;i++)
459     {
460     MG_SEGMENT* seg=nonmanifoldarete.get(i);
461     MG_NOEUD* n1=seg->get_noeud1();
462     MG_NOEUD* n2=seg->get_noeud2();
463     nonmanifoldnoeud.supprimer(n1);
464     nonmanifoldnoeud.supprimer(n2);
465     }
466 francois 222
467 francois 224
468     cout << "Non manifold par arete " << nonmanifoldarete.get_nb() << endl;
469     cout << "Non manifold par noeud " << nonmanifoldnoeud.get_nb() << endl;;
470 picher 230
471     //Étape 8 - Traitement des cas de non manifold
472     // À faire !!!!
473    
474     //Étape 9 - Suppression des peaux de très faible taille
475     int nbtri_total = mai2->get_nb_mg_triangle();
476     for (int i=0;i<nbpeau;i++)
477     {
478     int nbtri=lst_peau[i]->size();
479     if (nbtri < 0.005*nbtri_total)
480     {
481     for (int j=0;j<nbtri;j++)
482     {
483     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
484     int id = tripeau->get_id();
485     cout << "tripeau id:" << id <<endl;
486     mai2->supprimer_mg_triangleid(id);
487     }
488     }
489     }
490 francois 222 }
491    
492 francois 224
493     void MG_LISSAGE::determine_peau(std::vector<MG_TRIANGLE_PEAU*> * peau)
494 francois 222 {
495 francois 224 int num=0;
496     while (num!=peau->size())
497     {
498     MG_TRIANGLE_PEAU* p=(*peau)[num];
499     if (p->get_voisin1()->get_nouveau_numero()==0)
500     {
501     peau->push_back(p->get_voisin1());
502     p->get_voisin1()->change_nouveau_numero(1);
503     }
504     if (p->get_voisin2()->get_nouveau_numero()==0)
505     {
506     peau->push_back(p->get_voisin2());
507     p->get_voisin2()->change_nouveau_numero(1);
508     }
509     if (p->get_voisin3()->get_nouveau_numero()==0)
510     {
511     peau->push_back(p->get_voisin3());
512     p->get_voisin3()->change_nouveau_numero(1);
513     }
514     num++;
515     }
516     }
517    
518    
519     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)
520     {
521 francois 222 MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
522     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
523     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
524     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
525     if (mgsegment1==NULL)
526 francois 224 mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,origine);
527 francois 222 if (mgsegment2==NULL)
528 francois 224 mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,origine);
529 francois 222 if (mgsegment3==NULL)
530 francois 224 mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,origine);
531     MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
532 francois 222 mg_maillage->ajouter_mg_triangle(mtriangle);
533     return mtriangle;
534     }
535    
536 francois 224 MG_TRIANGLE_PEAU* MG_LISSAGE::recherche_voisin(MG_NOEUD* mg_noeud1,MG_NOEUD* mg_noeud2,MG_TRIANGLE_PEAU* triref)
537 francois 222 {
538 francois 224 MG_TRIANGLE_PEAU* trisol=NULL;
539     double angleref=4.*M_PI;
540 francois 222 int nb1=mg_noeud1->get_lien_triangle()->get_nb();
541     int nb2=mg_noeud2->get_lien_triangle()->get_nb();
542     for (int i=0;i<nb1;i++)
543     for (int j=0;j<nb2;j++)
544 francois 224 if (mg_noeud1->get_lien_triangle()->get(i)==mg_noeud2->get_lien_triangle()->get(j))
545     if ((MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i)!=triref)
546     {
547     double angle=calcul_angle(triref,(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i));
548     if (angle<angleref)
549     {
550     angleref=angle;
551     trisol=(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i);
552     }
553     }
554     return trisol;
555 francois 222 }
556 francois 223
557 francois 224
558    
559     double MG_LISSAGE::calcul_angle(MG_TRIANGLE_PEAU* ft1,MG_TRIANGLE_PEAU* ft2)
560 francois 223 {
561 francois 224 MG_NOEUD* noeud1=ft1->get_noeud1();
562     MG_NOEUD* noeud2=ft1->get_noeud2();
563     MG_NOEUD* noeud3=ft1->get_noeud3();
564     MG_NOEUD* noeud4=ft2->get_noeud1();
565     MG_NOEUD* noeud5=ft2->get_noeud2();
566     MG_NOEUD* noeud6=ft2->get_noeud3();
567    
568     MG_NOEUD* noeuda=NULL;
569     MG_NOEUD* noeudb=NULL;
570     MG_NOEUD* noeudc=NULL;
571     MG_NOEUD* noeudd=NULL;
572    
573    
574     if (!(noeud1==noeud4))
575     if (!(noeud1==noeud5))
576     if (!(noeud1==noeud6))
577     {
578     noeuda=noeud1;
579     noeudb=noeud2;
580     noeudc=noeud3;
581     }
582     if (!(noeud2==noeud4))
583     if (!(noeud2==noeud5))
584     if (!(noeud2==noeud6))
585     {
586     noeuda=noeud2;
587     noeudb=noeud1;
588     noeudc=noeud3;
589     }
590     if (!(noeud3==noeud4))
591     if (!(noeud3==noeud5))
592     if (!(noeud3==noeud6))
593     {
594     noeuda=noeud3;
595     noeudb=noeud1;
596     noeudc=noeud2;
597     }
598     if (noeuda==NULL) return 2*M_PI;
599     if (!(noeud4==noeudb))
600     if (!(noeud4==noeudc)) noeudd=noeud4;
601     if (!(noeud5==noeudb))
602     if (!(noeud5==noeudc)) noeudd=noeud5;
603     if (!(noeud6==noeudb))
604     if (!(noeud6==noeudc)) noeudd=noeud6;
605    
606     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()));
607     OT_VECTEUR_3D bc(noeudc->get_x()-noeudb->get_x(),noeudc->get_y()-noeudb->get_y(),noeudc->get_z()-noeudb->get_z());
608     bc.norme();
609     double xyz[3];
610     xyz[0]=noeuda->get_x();
611     xyz[1]=noeuda->get_y();
612     xyz[2]=noeuda->get_z();
613     OT_VECTEUR_3D n1(xyz,milieu.get_xyz());
614     n1.norme();
615     OT_VECTEUR_3D n=bc&n1;
616     n1=n&bc;
617     n1.norme();
618     xyz[0]=noeudd->get_x();
619     xyz[1]=noeudd->get_y();
620     xyz[2]=noeudd->get_z();
621     OT_VECTEUR_3D n2(xyz,milieu.get_xyz());
622     n2.norme();
623     n=bc&n2;
624     n2=n&bc;
625     n2.norme();
626    
627     double ps=n1*n2;
628     if (ps>1.) ps=1.;
629     if (ps<-1.) ps=-1.;
630     double angle=acos(ps);
631     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
632     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
633     OT_VECTEUR_3D nbnd(noeudd->get_x()-noeudb->get_x(),noeudd->get_y()-noeudb->get_y(),noeudd->get_z()-noeudb->get_z());
634     n=n1n2&n1n3;
635     n.norme();
636     nbnd.norme();
637    
638     ps=n*nbnd;
639     if (ps<-0.0001) angle=2*M_PI-angle;
640     return angle;
641    
642    
643 francois 223 }
644 picher 230
645     double MG_LISSAGE::ponderation_gaussian(double s,double sigma)
646     {
647     double w_s;
648     w_s = exp(-(s*s)/(2.*sigma*sigma));
649     return w_s;
650     }
651    
652     double MG_LISSAGE::ponderation_laplacian(double s,double sigma)
653     {
654     double w_s;
655     w_s = exp(-(sqrt(2.)*fabs(s))/sigma);
656     return w_s;
657     }
658    
659     double MG_LISSAGE::ponderation_elfallahford(double s,double sigma)
660     {
661     double w_s;
662     w_s = 1./sqrt(1+pow((s/sigma),2));
663     return w_s;
664     }