1 // MEFISTO2: a library to compute 2D triangulation from segmented boundaries
4 // Copyright (C) 2006 Laboratoire J.-L. Lions UPMC Paris
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
22 // File : aptrte.cxx le C++ de l'appel du trianguleur plan
24 // Author : Alain PERRONNET
25 // Date : 13 novembre 2006
29 #include "utilities.h"
41 areteideale()//( R3 xyz, R3 direction )
46 //calcul de la longueur ideale de l'arete au sommet xyz (z ici inactif)
47 //dans la direction donnee
48 //a ajuster pour chaque surface plane et selon l'entier notysu (voir plus bas)
51 static double cpunew, cpuold=0;
53 void tempscpu_( double & tempsec )
54 //Retourne le temps CPU utilise en secondes
56 tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
57 //MESSAGE( "temps cpu=" << tempsec );
61 void deltacpu_( R & dtcpu )
62 //Retourne le temps CPU utilise en secondes depuis le precedent appel
65 dtcpu = R( cpunew - cpuold );
67 //MESSAGE( "delta temps cpu=" << dtcpu );
72 void aptrte( Z nutysu, R aretmx,
73 Z nblf, Z * nudslf, R2 * uvslf,
75 Z & nbst, R2 * & uvst,
78 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
79 // but : appel de la triangulation par un arbre-4 recouvrant
80 // ----- de triangles equilateraux
81 // le contour du domaine plan est defini par des lignes fermees
82 // la premiere ligne etant l'enveloppe de toutes les autres
83 // la fonction areteideale(s,d) donne la taille d'arete
84 // au point s dans la direction (actuellement inactive) d
85 // des lors toute arete issue d'un sommet s devrait avoir une longueur
86 // comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
89 // Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
90 // De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
94 // nutysu : numero de traitement de areteideale_(s,d) selon le type de surface
95 // 0 pas d'emploi de la fonction areteideale_() et aretmx est active
96 // 1 il existe une fonction areteideale_(s,d)
97 // dont seules les 2 premieres composantes de uv sont actives
98 // ... autres options a definir ...
99 // aretmx : longueur maximale des aretes de la future triangulation
100 // nblf : nombre de lignes fermees de la surface
101 // nudslf : numero du dernier sommet de chacune des nblf lignes fermees
102 // nudslf(0)=0 pour permettre la difference sans test
103 // Attention le dernier sommet de chaque ligne est raccorde au premier
104 // tous les sommets et les points internes ont des coordonnees
105 // UV differentes <=> Pas de point double!
106 // uvslf : uv des nudslf(nblf) sommets des lignes fermees
107 // nbpti : nombre de points internes futurs sommets de la triangulation
108 // uvpti : uv des points internes futurs sommets de la triangulation
112 // nbst : nombre de sommets de la triangulation finale
113 // uvst : coordonnees uv des nbst sommets de la triangulation
114 // nbt : nombre de triangles de la triangulation finale
115 // nust : 4 numeros dans uvst des sommets des nbt triangles
116 // s1, s2, s3, 0: no dans uvst des 3 sommets et 0 car quadrangle!
117 // ierr : 0 si pas d'erreur
119 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
120 // auteur : Alain Perronnet Laboratoire J.-L. LIONS Paris UPMC mars 2006
121 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
123 Z nbsttria=4; //Attention: 4 sommets stockes par triangle
124 //no st1, st2, st3, 0 (non quadrangle)
127 R3 direction=R3(0,0,0); //direction pour areteideale() inactive ici!
128 Z nbarfr=nudslf[nblf]; //nombre total d'aretes des lignes fermees
129 Z mxtrou = Max( 1024, nblf ); //nombre maximal de trous dans la surface
132 Z *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
133 Z *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
134 Z *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE
135 Z *mnqueu=NULL, mxqueu;
137 Z *mnarcf=NULL, mxarcf;
146 R3 comxmi[2]; //coordonnees UV Min et Maximales
147 R aremin, aremax; //longueur minimale et maximale des aretes
148 R airemx; //aire maximale souhaitee d'un triangle
152 Z i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
153 Z mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
157 // initialisation du temps cpu
161 // quelques reservations de tableaux pour faire les calculs
162 // ========================================================
163 // declaration du tableau des coordonnees des sommets de la frontiere
164 // puis des sommets internes ajoutes
165 // majoration empirique du nombre de sommets de la triangulation
167 mxsomm = Max( 20000, 64*nbpti+i*i );
168 MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
169 MESSAGE( "nutysu=" << nutysu << " aretmx=" << aretmx
170 << " mxsomm=" << mxsomm );
171 MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
174 //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
175 if( mnpxyd!=NULL ) delete [] mnpxyd;
176 mnpxyd = new R3[mxsomm];
177 if( mnpxyd==NULL ) goto ERREUR;
179 // le tableau mnsoar des aretes des triangles
180 // 1: sommet 1 dans pxyd,
181 // 2: sommet 2 dans pxyd,
182 // 3: numero de 1 a nblf de la ligne qui supporte l'arete
183 // 4: numero dans mnartr du triangle 1 partageant cette arete,
184 // 5: numero dans mnartr du triangle 2 partageant cette arete,
185 // 6: chainage des aretes frontalieres ou internes ou
186 // des aretes simples des etoiles de triangles,
187 // 7: chainage du hachage des aretes
188 // nombre d'aretes = 3 ( nombre de sommets - 1 + nombre de trous )
189 // pour le hachage des aretes mxsoar doit etre > 3*mxsomm!
190 // h(ns1,ns2) = min( ns1, ns2 )
191 if( mnsoar!=NULL ) delete [] mnsoar;
192 mxsoar = 3 * ( mxsomm + mxtrou );
193 mnsoar = new Z[mosoar*mxsoar];
194 if( mnsoar==NULL ) goto ERREUR;
195 //initialiser le tableau mnsoar pour le hachage des aretes
196 insoar( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
198 // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
199 if( mnarst!=NULL ) delete [] mnarst;
200 mnarst = new Z[1+mxsomm];
201 if( mnarst==NULL ) goto ERREUR;
205 // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
206 // ou no du point si interne forc'e par l'utilisateur
207 // ou 0 si interne cree par le module
208 if( mnslig!=NULL ) delete [] mnslig;
209 mnslig = new Z[mxsomm];
210 if( mnslig==NULL ) goto ERREUR;
211 azeroi( mxsomm, mnslig );
213 // initialisation des aretes frontalieres de la triangulation future
214 // renumerotation des sommets des aretes des lignes pour la triangulation
215 // mise a l'echelle des coordonnees des sommets pour obtenir une
216 // meilleure precision lors des calculs + quelques verifications
217 // boucle sur les lignes fermees qui forment la frontiere
218 // ======================================================================
223 for (n=1; n<=nblf; n++)
225 //l'initialisation de la premiere arete de la ligne n dans la triangulation
226 //-------------------------------------------------------------------------
227 //le sommet ns0 est le numero de l'origine de la ligne
229 mnpxyd[ns0].x = uvslf[ns0].x;
230 mnpxyd[ns0].y = uvslf[ns0].y;
231 mnpxyd[ns0].z = areteideale();//( mnpxyd[ns0], direction );
232 // MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
233 // << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
235 //carre de la longueur de l'arete 1 de la ligne fermee n
236 d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 )
237 + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
238 aremin = Min( aremin, d );
239 aremax = Max( aremax, d );
241 //le numero des 2 sommets (ns1,ns2) de la premiere arete de la ligne
242 //initialisation de la 1-ere arete ns1-ns1+1 de cette ligne fermee n
243 //le numero des 2 sommets ns1 ns2 de la 1-ere arete
244 //Attention: les numeros ns debutent a 1 (ils ont >0)
245 // les tableaux c++ demarrent a zero!
246 // les tableaux fortran demarrent ou l'on veut!
251 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
252 mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
253 fasoar( ns1, ns2, moins1, moins1, n,
254 mosoar, mxsoar, n1soar, mnsoar, mnarst,
256 //pas de test sur ierr car pas de saturation possible a ce niveau
258 //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
261 //la nouvelle arete est la suivante de l'arete definie juste avant
263 mnsoar[mosoar * noar - mosoar + 5] = noar0;
265 //l'initialisation des aretes suivantes de la ligne dans la triangulation
266 //-----------------------------------------------------------------------
267 nbarli = nudslf[n] - nudslf[n-1]; //nombre d'aretes=sommets de la ligne n
268 for (i=2; i<=nbarli; i++)
270 ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
272 //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
275 //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
278 //l'arete precedente est dotee de sa suivante:celle cree ensuite
279 //les 2 coordonnees du sommet ns2 de la ligne
281 //debut ajout 5/10/2006 ................................................
282 nuds = Max( nuds, ns ); //le numero du dernier sommet traite
283 //fin ajout 5/10/2006 ................................................
284 mnpxyd[ns].x = uvslf[ns].x;
285 mnpxyd[ns].y = uvslf[ns].y;
286 mnpxyd[ns].z = areteideale();//( mnpxyd[ns], direction );
287 // MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
288 // << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
290 //carre de la longueur de l'arete
291 d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2)
292 + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
293 aremin = Min( aremin, d );
294 aremax = Max( aremax, d );
296 //debut ajout du 5/10/2006 .............................................
297 //la longueur de l'arete ns1-ns2
299 //longueur arete = Min ( aretmx, aretes incidentes )
300 mnpxyd[ns ].z = Min( mnpxyd[ns ].z, d );
301 mnpxyd[ns2-1].z = Min( mnpxyd[ns2-1].z, d );
302 //fin ajout du 5/10/2006 ...............................................
304 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
305 mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
307 //ajout de l'arete dans la liste
308 fasoar( ns1, ns2, moins1, moins1, n,
309 mosoar, mxsoar, n1soar, mnsoar,
310 mnarst, noar, ierr );
311 //pas de test sur ierr car pas de saturation possible a ce niveau
313 //chainage des aretes frontalieres en position 6 du tableau mnsoar
314 //la nouvelle arete est la suivante de l'arete definie juste avant
315 mnsoar[ mosoar * noar0 - mosoar + 5 ] = noar;
318 //attention: la derniere arete de la ligne fermee enveloppe
319 // devient en fait la premiere arete de cette ligne
320 // dans le chainage des aretes de la frontiere!
322 if( ierr != 0 ) goto ERREUR;
324 aremin = sqrt( aremin ); //longueur minimale d'une arete des lignes fermees
325 aremax = sqrt( aremax ); //longueur maximale d'une arete
327 //debut ajout 9/11/2006 ................................................
328 // devenu un commentaire aretmx = Min( aretmx, aremax ); //pour homogeneiser
330 // protection contre une arete max desiree trop grande ou trop petite
331 if( aretmx > aremax*2.05 ) aretmx = aremax;
333 // protection contre une arete max desiree trop petite
334 if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 )
335 aretmx =(aremin+aremax*2)/3.0;
337 if( aretmx < aremin && aremin > 0 )
340 //sauvegarde pour la fonction areteideale_
341 aretemaxface_ = aretmx;
343 //aire maximale souhaitee des triangles
344 airemx = aretmx * aretmx * sqrt(3.0) / 2.0; //Aire triangle equilateral
346 for(i=0; i<=nuds; i++ )
347 mnpxyd[i].z = Min( mnpxyd[i].z, aretmx );
348 //MESSAGE("Numero du dernier sommet frontalier=" << nuds+1);
349 //fin ajout 9/11/2006 .................................................
352 MESSAGE("Sur le bord: arete min=" << aremin << " arete max=" << aremax );
353 MESSAGE("Triangulation: arete mx=" << aretmx
354 << " triangle aire mx=" << airemx );
356 //chainage des aretes frontalieres : la derniere arete frontaliere
357 mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
359 //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
360 //reservation du tableau des numeros des 3 aretes de chaque triangle
361 //mnartr( moartr, mxartr )
362 //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
363 // 3Triangles = 2 Aretes internes + Aretes frontalieres
364 // d'ou 3T/2 < AI + AF => T < 3T/2 - Sommets + 1-Trous
365 //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
366 if( mnartr!=NULL ) delete [] mnartr;
367 mxartr = 2 * ( mxsomm + mxtrou );
368 mnartr = new Z[moartr*mxartr];
369 if( mnartr==NULL ) goto ERREUR;
371 //Ajout des points internes
372 ns1 = nudslf[ nblf ];
373 for (i=0; i<nbpti; i++)
375 //les 2 coordonnees du point i de sommet nbs
376 mnpxyd[ns1].x = uvpti[i].x;
377 mnpxyd[ns1].y = uvpti[i].y;
378 mnpxyd[ns1].z = areteideale();//( mnpxyd[ns1], direction );
379 //le numero i du point interne
384 //nombre de sommets de la frontiere et internes
387 // creation de l'arbre-4 des te (tableau letree)
388 // ajout dans les te des sommets des lignes et des points internes imposes
389 // =======================================================================
390 // premiere estimation de mxtree
393 NEWTREE: //en cas de saturation de l'un des tableaux, on boucle
394 MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
395 if( mntree != NULL ) delete [] mntree;
397 mntree = new Z[motree*(1+mxtree)];
398 if( mntree==NULL ) goto ERREUR;
400 //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
401 teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
407 //saturation de letree => sa taille est augmentee et relance
410 MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
416 MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
417 if( ierr != 0 ) goto ERREUR;
418 //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
420 // homogeneisation de l'arbre des te a un saut de taille au plus
421 // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
422 // ===========================================================================
423 // reservation de la queue pour parcourir les te de l'arbre
424 if( mnqueu != NULL ) delete [] mnqueu;
426 mnqueu = new Z[mxqueu];
427 if( mnqueu==NULL) goto ERREUR;
429 tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
431 mntree, mxqueu, mnqueu,
436 MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
437 << d << " secondes");
440 //destruction du tableau auxiliaire et de l'arbre
445 MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
453 // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
454 // et des points de la frontiere, des points internes imposes interieurs
455 // ==========================================================================
456 tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
457 mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
458 moartr, mxartr, n1artr, mnartr, mnarst,
461 // destruction de la queue et de l'arbre devenus inutiles
462 delete [] mnqueu; mnqueu=NULL;
463 delete [] mntree; mntree=NULL;
468 MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
470 // ierr =0 si pas d'erreur
471 // =1 si le tableau mnsoar est sature
472 // =2 si le tableau mnartr est sature
473 // =3 si aucun des triangles ne contient l'un des points internes
474 // =5 si saturation de la queue de parcours de l'arbre des te
475 if( ierr != 0 ) goto ERREUR;
477 //qualites de la triangulation actuelle
478 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
479 nbt, quamoy, quamin );
481 // boucle sur les aretes internes (non sur une ligne de la frontiere)
482 // avec echange des 2 diagonales afin de rendre la triangulation delaunay
483 // ======================================================================
484 // formation du chainage 6 des aretes internes a echanger eventuellement
485 aisoar( mosoar, mxsoar, mnsoar, na );
486 tedela( mnpxyd, mnarst,
487 mosoar, mxsoar, n1soar, mnsoar, na,
488 moartr, mxartr, n1artr, mnartr, n );
490 MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
493 MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
494 << d << " secondes");
496 //qualites de la triangulation actuelle
497 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
498 nbt, quamoy, quamin );
500 // detection des aretes frontalieres initiales perdues
501 // triangulation frontale pour les restaurer
502 // ===================================================
504 if( mn1arcf != NULL ) delete [] mn1arcf;
505 if( mnarcf != NULL ) delete [] mnarcf;
506 if( mnarcf1 != NULL ) delete [] mnarcf1;
507 if( mnarcf2 != NULL ) delete [] mnarcf2;
508 mn1arcf = new Z[1+mxarcf];
509 if( mn1arcf == NULL ) goto ERREUR;
510 mnarcf = new Z[3*mxarcf];
511 if( mnarcf == NULL ) goto ERREUR;
512 mnarcf1 = new Z[mxarcf];
513 if( mnarcf1 == NULL ) goto ERREUR;
514 mnarcf2 = new Z[mxarcf];
515 if( mnarcf2 == NULL ) goto ERREUR;
517 terefr( nbarpi, mnpxyd,
518 mosoar, mxsoar, n1soar, mnsoar,
519 moartr, mxartr, n1artr, mnartr, mnarst,
520 mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
523 MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere ierr=" << ierr );
526 MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
527 << d << " secondes");
529 if( ierr != 0 ) goto ERREUR;
531 //qualites de la triangulation actuelle
532 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
533 nbt, quamoy, quamin );
535 // fin de la triangulation avec respect des aretes initiales frontalieres
537 // suppression des triangles externes a la surface
538 // ===============================================
539 // recherche du dernier triangle utilise
540 mn = mxartr * moartr;
541 for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
544 if( mnartr[mn] != 0 ) break;
547 if( mntrsu != NULL ) delete [] mntrsu;
548 mntrsu = new Z[ndtri0];
549 if( mntrsu == NULL ) goto ERREUR;
551 if( mnlftr != NULL ) delete [] mnlftr;
552 mnlftr = new Z[nblf];
553 if( mnlftr == NULL ) goto ERREUR;
555 for (n=0; n<nblf; n++) //numero de la ligne fermee de 1 a nblf
558 tesuex( nblf, mnlftr,
559 ndtri0, nbsomm, mnpxyd, mnslig,
560 mosoar, mxsoar, mnsoar,
561 moartr, mxartr, n1artr, mnartr, mnarst,
564 delete [] mnlftr; mnlftr=NULL;
565 delete [] mntrsu; mntrsu=NULL;
569 MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
570 if( ierr != 0 ) goto ERREUR;
572 //qualites de la triangulation actuelle
573 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
574 nbt, quamoy, quamin );
576 // amelioration de la qualite de la triangulation par
577 // barycentrage des sommets internes a la triangulation
578 // suppression des aretes trop longues ou trop courtes
579 // modification de la topologie des groupes de triangles
580 // mise en delaunay de la triangulation
581 // =====================================================
582 mnarcf3 = new Z[mxarcf];
583 if( mnarcf3 == NULL )
585 MESSAGE ( "aptrte: MC saturee mnarcf3=" << mnarcf3 );
588 teamqt( nutysu, aretmx, airemx,
589 mnarst, mosoar, mxsoar, n1soar, mnsoar,
590 moartr, mxartr, n1artr, mnartr,
591 mxarcf, mnarcf2, mnarcf3,
592 mn1arcf, mnarcf, mnarcf1,
593 nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
595 if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
596 if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
597 if( mnarcf != NULL ) {delete [] mnarcf; mnarcf =NULL;}
598 if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
599 if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
603 MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
604 if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
605 if( ierr != 0 ) goto ERREUR;
607 //qualites de la triangulation finale
608 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
609 nbt, quamoy, quamin );
611 // renumerotation des sommets internes: mnarst(i)=numero final du sommet
612 // ===================================
613 for (i=0; i<=nbsomm; i++)
616 for (nt=1; nt<=mxartr; nt++)
618 if( mnartr[nt*moartr-moartr] != 0 )
620 //le numero des 3 sommets du triangle nt
621 nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
622 //les 3 sommets du triangle sont actifs
623 mnarst[ nosotr[0] ] = 1;
624 mnarst[ nosotr[1] ] = 1;
625 mnarst[ nosotr[2] ] = 1;
629 for (i=1; i<=nbsomm; i++)
635 // generation du tableau uvst de la surface triangulee
636 // ---------------------------------------------------
637 if( uvst != NULL ) delete [] uvst;
639 if( uvst == NULL ) goto ERREUR;
642 for (i=0; i<nbsomm; i++ )
647 uvst[nbst].x = mnpxyd[i].x;
648 uvst[nbst].y = mnpxyd[i].y;
650 //si le sommet est un point ou appartient a une ligne
651 //ses coordonnees initiales sont restaurees
658 //retour aux coordonnees initiales dans uvslf
660 n = n - 1000000 * l + nudslf[l-1] - 1;
661 uvst[nbst].x = uvslf[n].x;
662 uvst[nbst].y = uvslf[n].y;
666 //point utilisateur n interne impose
667 //retour aux coordonnees initiales dans uvpti
668 uvst[nbst].x = uvpti[n-1].x;
669 uvst[nbst].y = uvpti[n-1].y;
676 // generation du tableau 'nsef' de la surface triangulee
677 // -----------------------------------------------------
678 // boucle sur les triangles occupes (internes et externes)
679 if( nust != NULL ) delete [] nust;
680 nust = new Z[nbsttria*nbt];
681 if( nust == NULL ) goto ERREUR;
683 for (i=1; i<=mxartr; i++)
685 //le triangle i de mnartr
686 if( mnartr[i*moartr-moartr] != 0 )
688 //le triangle i est interne => nosotr numero de ses 3 sommets
689 nusotr( i, mosoar, mnsoar, moartr, mnartr, nosotr );
690 nust[nbt++] = mnarst[ nosotr[0] ];
691 nust[nbt++] = mnarst[ nosotr[1] ];
692 nust[nbt++] = mnarst[ nosotr[2] ];
696 nbt /= nbsttria; //le nombre final de triangles de la surface
697 MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
698 << nbt << " triangles" );
701 MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
703 // destruction des tableaux auxiliaires
704 // ------------------------------------
706 if( mnarst != NULL ) delete [] mnarst;
707 if( mnartr != NULL ) delete [] mnartr;
708 if( mnslig != NULL ) delete [] mnslig;
709 if( mnsoar != NULL ) delete [] mnsoar;
710 if( mnpxyd != NULL ) delete [] mnpxyd;
711 if( mntree != NULL ) delete [] mntree;
712 if( mnqueu != NULL ) delete [] mnqueu;
713 if( mntrsu != NULL ) delete [] mntrsu;
714 if( mnlftr != NULL ) delete [] mnlftr;
715 if( mn1arcf != NULL ) delete [] mn1arcf;
716 if( mnarcf != NULL ) delete [] mnarcf;
717 if( mnarcf1 != NULL ) delete [] mnarcf1;
718 if( mnarcf2 != NULL ) delete [] mnarcf2;
719 if( mnarcf3 != NULL ) delete [] mnarcf3;
723 if( ierr == 51 || ierr == 52 )
725 //saturation des sommets => redepart avec 2 fois plus de sommets
732 MESSAGE( "APTRTE: Triangulation NON REALISEE avec erreur=" << ierr );
733 if( ierr == 0 ) ierr=1;
741 qualitetrte( R3 *mnpxyd,
742 Z & mosoar, Z & mxsoar, Z *mnsoar,
743 Z & moartr, Z & mxartr, Z *mnartr,
744 Z & nbtria, R & quamoy, R & quamin )
745 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
746 // but : calculer la qualite moyenne et minimale de la triangulation
747 // ----- actuelle definie par les tableaux mnsoar et mnartr
750 // mnpxyd : tableau des coordonnees 2d des points
751 // par point : x y distance_souhaitee
752 // mosoar : nombre maximal d'entiers par arete et
753 // indice dans mnsoar de l'arete suivante dans le hachage
754 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
755 // attention: mxsoar>3*mxsomm obligatoire!
756 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
757 // chainage des aretes frontalieres, chainage du hachage des aretes
758 // hachage des aretes = mnsoar(1)+mnsoar(2)*2
759 // avec mxsoar>=3*mxsomm
760 // une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
761 // mnsoar(2,arete vide)=l'arete vide qui precede
762 // mnsoar(3,arete vide)=l'arete vide qui suit
763 // moartr : nombre maximal d'entiers par arete du tableau mnartr
764 // mxartr : nombre maximal de triangles declarables
765 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
766 // arete1 = 0 si triangle vide => arete2 = triangle vide suivant
769 // nbtria : nombre de triangles internes au domaine
770 // quamoy : qualite moyenne des triangles actuels
771 // quamin : qualite minimale des triangles actuels
772 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
775 Z nosotr[3], mn, nbtrianeg, nt, ntqmin;
785 for ( nt=1; nt<=mxartr; nt++ )
790 //un triangle occupe de plus
793 //le numero des 3 sommets du triangle nt
794 nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
796 //la qualite du triangle ns1 ns2 ns3
797 qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
803 //la qualite minimale
804 if( qualite < quamin )
810 //aire signee du triangle nt
811 d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
814 //un triangle d'aire negative de plus
816 MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
817 << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
818 << " a une aire " << d <<"<=0");
821 //aire des triangles actuels
828 MESSAGE("Qualite moyenne=" << quamoy
829 << " Qualite minimale=" << quamin
830 << " des " << nbtria << " triangles de surface plane totale="
835 //le numero des 3 sommets du triangle ntqmin de qualite minimale
836 nusotr(ntqmin, mosoar, mnsoar, moartr, mnartr, nosotr );
837 MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
838 <<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
839 for (int i=0;i<3;i++)
840 MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
841 <<" y="<< mnpxyd[nosotr[i]-1].y);
845 MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );