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"
44 areteideale()//( R3 xyz, R3 direction )
49 //calcul de la longueur ideale de l'arete au sommet xyz (z ici inactif)
50 //dans la direction donnee
51 //a ajuster pour chaque surface plane et selon l'entier notysu (voir plus bas)
54 static double cpunew, cpuold=0;
63 tempscpu_( double & tempsec )
64 //Retourne le temps CPU utilise en secondes
66 tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
67 //MESSAGE( "temps cpu=" << tempsec );
78 deltacpu_( R & dtcpu )
79 //Retourne le temps CPU utilise en secondes depuis le precedent appel
82 dtcpu = R( cpunew - cpuold );
84 //MESSAGE( "delta temps cpu=" << dtcpu );
89 void aptrte( Z nutysu, R aretmx,
90 Z nblf, Z * nudslf, R2 * uvslf,
92 Z & nbst, R2 * & uvst,
95 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
96 // but : appel de la triangulation par un arbre-4 recouvrant
97 // ----- de triangles equilateraux
98 // le contour du domaine plan est defini par des lignes fermees
99 // la premiere ligne etant l'enveloppe de toutes les autres
100 // la fonction areteideale(s,d) donne la taille d'arete
101 // au point s dans la direction (actuellement inactive) d
102 // des lors toute arete issue d'un sommet s devrait avoir une longueur
103 // comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
106 // Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
107 // De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
111 // nutysu : numero de traitement de areteideale_(s,d) selon le type de surface
112 // 0 pas d'emploi de la fonction areteideale_() et aretmx est active
113 // 1 il existe une fonction areteideale_(s,d)
114 // dont seules les 2 premieres composantes de uv sont actives
115 // ... autres options a definir ...
116 // aretmx : longueur maximale des aretes de la future triangulation
117 // nblf : nombre de lignes fermees de la surface
118 // nudslf : numero du dernier sommet de chacune des nblf lignes fermees
119 // nudslf(0)=0 pour permettre la difference sans test
120 // Attention le dernier sommet de chaque ligne est raccorde au premier
121 // tous les sommets et les points internes ont des coordonnees
122 // UV differentes <=> Pas de point double!
123 // uvslf : uv des nudslf(nblf) sommets des lignes fermees
124 // nbpti : nombre de points internes futurs sommets de la triangulation
125 // uvpti : uv des points internes futurs sommets de la triangulation
129 // nbst : nombre de sommets de la triangulation finale
130 // uvst : coordonnees uv des nbst sommets de la triangulation
131 // nbt : nombre de triangles de la triangulation finale
132 // nust : 4 numeros dans uvst des sommets des nbt triangles
133 // s1, s2, s3, 0: no dans uvst des 3 sommets et 0 car quadrangle!
134 // ierr : 0 si pas d'erreur
136 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
137 // auteur : Alain Perronnet Laboratoire J.-L. LIONS Paris UPMC mars 2006
138 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
140 Z nbsttria=4; //Attention: 4 sommets stockes par triangle
141 //no st1, st2, st3, 0 (non quadrangle)
144 R3 direction=R3(0,0,0); //direction pour areteideale() inactive ici!
145 Z nbarfr=nudslf[nblf]; //nombre total d'aretes des lignes fermees
146 Z mxtrou = Max( 1024, nblf ); //nombre maximal de trous dans la surface
149 Z *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
150 Z *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
151 Z *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE
152 Z *mnqueu=NULL, mxqueu;
154 Z *mnarcf=NULL, mxarcf;
163 R3 comxmi[2]; //coordonnees UV Min et Maximales
164 R aremin, aremax; //longueur minimale et maximale des aretes
165 R airemx; //aire maximale souhaitee d'un triangle
169 Z i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
170 Z mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
174 // initialisation du temps cpu
178 // quelques reservations de tableaux pour faire les calculs
179 // ========================================================
180 // declaration du tableau des coordonnees des sommets de la frontiere
181 // puis des sommets internes ajoutes
182 // majoration empirique du nombre de sommets de la triangulation
184 mxsomm = Max( 20000, 64*nbpti+i*i );
185 MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
186 MESSAGE( "nutysu=" << nutysu << " aretmx=" << aretmx
187 << " mxsomm=" << mxsomm );
188 MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
191 //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
192 if( mnpxyd!=NULL ) delete [] mnpxyd;
193 mnpxyd = new R3[mxsomm];
194 if( mnpxyd==NULL ) goto ERREUR;
196 // le tableau mnsoar des aretes des triangles
197 // 1: sommet 1 dans pxyd,
198 // 2: sommet 2 dans pxyd,
199 // 3: numero de 1 a nblf de la ligne qui supporte l'arete
200 // 4: numero dans mnartr du triangle 1 partageant cette arete,
201 // 5: numero dans mnartr du triangle 2 partageant cette arete,
202 // 6: chainage des aretes frontalieres ou internes ou
203 // des aretes simples des etoiles de triangles,
204 // 7: chainage du hachage des aretes
205 // nombre d'aretes = 3 ( nombre de sommets - 1 + nombre de trous )
206 // pour le hachage des aretes mxsoar doit etre > 3*mxsomm!
207 // h(ns1,ns2) = min( ns1, ns2 )
208 if( mnsoar!=NULL ) delete [] mnsoar;
209 mxsoar = 3 * ( mxsomm + mxtrou );
210 mnsoar = new Z[mosoar*mxsoar];
211 if( mnsoar==NULL ) goto ERREUR;
212 //initialiser le tableau mnsoar pour le hachage des aretes
213 insoar( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
215 // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
216 if( mnarst!=NULL ) delete [] mnarst;
217 mnarst = new Z[1+mxsomm];
218 if( mnarst==NULL ) goto ERREUR;
222 // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
223 // ou no du point si interne forc'e par l'utilisateur
224 // ou 0 si interne cree par le module
225 if( mnslig!=NULL ) delete [] mnslig;
226 mnslig = new Z[mxsomm];
227 if( mnslig==NULL ) goto ERREUR;
228 azeroi( mxsomm, mnslig );
230 // initialisation des aretes frontalieres de la triangulation future
231 // renumerotation des sommets des aretes des lignes pour la triangulation
232 // mise a l'echelle des coordonnees des sommets pour obtenir une
233 // meilleure precision lors des calculs + quelques verifications
234 // boucle sur les lignes fermees qui forment la frontiere
235 // ======================================================================
240 for (n=1; n<=nblf; n++)
242 //l'initialisation de la premiere arete de la ligne n dans la triangulation
243 //-------------------------------------------------------------------------
244 //le sommet ns0 est le numero de l'origine de la ligne
246 mnpxyd[ns0].x = uvslf[ns0].x;
247 mnpxyd[ns0].y = uvslf[ns0].y;
248 mnpxyd[ns0].z = areteideale();//( mnpxyd[ns0], direction );
249 // MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
250 // << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
252 //carre de la longueur de l'arete 1 de la ligne fermee n
253 d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 )
254 + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
255 aremin = Min( aremin, d );
256 aremax = Max( aremax, d );
258 //le numero des 2 sommets (ns1,ns2) de la premiere arete de la ligne
259 //initialisation de la 1-ere arete ns1-ns1+1 de cette ligne fermee n
260 //le numero des 2 sommets ns1 ns2 de la 1-ere arete
261 //Attention: les numeros ns debutent a 1 (ils ont >0)
262 // les tableaux c++ demarrent a zero!
263 // les tableaux fortran demarrent ou l'on veut!
268 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
269 mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
270 fasoar( ns1, ns2, moins1, moins1, n,
271 mosoar, mxsoar, n1soar, mnsoar, mnarst,
273 //pas de test sur ierr car pas de saturation possible a ce niveau
275 //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
278 //la nouvelle arete est la suivante de l'arete definie juste avant
280 mnsoar[mosoar * noar - mosoar + 5] = noar0;
282 //l'initialisation des aretes suivantes de la ligne dans la triangulation
283 //-----------------------------------------------------------------------
284 nbarli = nudslf[n] - nudslf[n-1]; //nombre d'aretes=sommets de la ligne n
285 for (i=2; i<=nbarli; i++)
287 ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
289 //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
292 //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
295 //l'arete precedente est dotee de sa suivante:celle cree ensuite
296 //les 2 coordonnees du sommet ns2 de la ligne
298 //debut ajout 5/10/2006 ................................................
299 nuds = Max( nuds, ns ); //le numero du dernier sommet traite
300 //fin ajout 5/10/2006 ................................................
301 mnpxyd[ns].x = uvslf[ns].x;
302 mnpxyd[ns].y = uvslf[ns].y;
303 mnpxyd[ns].z = areteideale();//( mnpxyd[ns], direction );
304 // MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
305 // << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
307 //carre de la longueur de l'arete
308 d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2)
309 + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
310 aremin = Min( aremin, d );
311 aremax = Max( aremax, d );
313 //debut ajout du 5/10/2006 .............................................
314 //la longueur de l'arete ns1-ns2
316 //longueur arete = Min ( aretmx, aretes incidentes )
317 mnpxyd[ns ].z = Min( mnpxyd[ns ].z, d );
318 mnpxyd[ns2-1].z = Min( mnpxyd[ns2-1].z, d );
319 //fin ajout du 5/10/2006 ...............................................
321 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
322 mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
324 //ajout de l'arete dans la liste
325 fasoar( ns1, ns2, moins1, moins1, n,
326 mosoar, mxsoar, n1soar, mnsoar,
327 mnarst, noar, ierr );
328 //pas de test sur ierr car pas de saturation possible a ce niveau
330 //chainage des aretes frontalieres en position 6 du tableau mnsoar
331 //la nouvelle arete est la suivante de l'arete definie juste avant
332 mnsoar[ mosoar * noar0 - mosoar + 5 ] = noar;
335 //attention: la derniere arete de la ligne fermee enveloppe
336 // devient en fait la premiere arete de cette ligne
337 // dans le chainage des aretes de la frontiere!
339 if( ierr != 0 ) goto ERREUR;
341 aremin = sqrt( aremin ); //longueur minimale d'une arete des lignes fermees
342 aremax = sqrt( aremax ); //longueur maximale d'une arete
344 //debut ajout 9/11/2006 ................................................
345 // devenu un commentaire aretmx = Min( aretmx, aremax ); //pour homogeneiser
347 // protection contre une arete max desiree trop grande ou trop petite
348 if( aretmx > aremax*2.05 ) aretmx = aremax;
350 // protection contre une arete max desiree trop petite
351 if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 )
352 aretmx =(aremin+aremax*2)/3.0;
354 if( aretmx < aremin && aremin > 0 )
357 //sauvegarde pour la fonction areteideale_
358 aretemaxface_ = aretmx;
360 //aire maximale souhaitee des triangles
361 airemx = aretmx * aretmx * sqrt(3.0) / 2.0; //Aire triangle equilateral
363 for(i=0; i<=nuds; i++ )
364 mnpxyd[i].z = Min( mnpxyd[i].z, aretmx );
365 //MESSAGE("Numero du dernier sommet frontalier=" << nuds+1);
366 //fin ajout 9/11/2006 .................................................
369 MESSAGE("Sur le bord: arete min=" << aremin << " arete max=" << aremax );
370 MESSAGE("Triangulation: arete mx=" << aretmx
371 << " triangle aire mx=" << airemx );
373 //chainage des aretes frontalieres : la derniere arete frontaliere
374 mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
376 //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
377 //reservation du tableau des numeros des 3 aretes de chaque triangle
378 //mnartr( moartr, mxartr )
379 //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
380 // 3Triangles = 2 Aretes internes + Aretes frontalieres
381 // d'ou 3T/2 < AI + AF => T < 3T/2 - Sommets + 1-Trous
382 //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
383 if( mnartr!=NULL ) delete [] mnartr;
384 mxartr = 2 * ( mxsomm + mxtrou );
385 mnartr = new Z[moartr*mxartr];
386 if( mnartr==NULL ) goto ERREUR;
388 //Ajout des points internes
389 ns1 = nudslf[ nblf ];
390 for (i=0; i<nbpti; i++)
392 //les 2 coordonnees du point i de sommet nbs
393 mnpxyd[ns1].x = uvpti[i].x;
394 mnpxyd[ns1].y = uvpti[i].y;
395 mnpxyd[ns1].z = areteideale();//( mnpxyd[ns1], direction );
396 //le numero i du point interne
401 //nombre de sommets de la frontiere et internes
404 // creation de l'arbre-4 des te (tableau letree)
405 // ajout dans les te des sommets des lignes et des points internes imposes
406 // =======================================================================
407 // premiere estimation de mxtree
410 NEWTREE: //en cas de saturation de l'un des tableaux, on boucle
411 MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
412 if( mntree != NULL ) delete [] mntree;
414 mntree = new Z[motree*(1+mxtree)];
415 if( mntree==NULL ) goto ERREUR;
417 //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
418 teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
424 //saturation de letree => sa taille est augmentee et relance
427 MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
433 MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
434 if( ierr != 0 ) goto ERREUR;
435 //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
437 // homogeneisation de l'arbre des te a un saut de taille au plus
438 // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
439 // ===========================================================================
440 // reservation de la queue pour parcourir les te de l'arbre
441 if( mnqueu != NULL ) delete [] mnqueu;
443 mnqueu = new Z[mxqueu];
444 if( mnqueu==NULL) goto ERREUR;
446 tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
448 mntree, mxqueu, mnqueu,
453 MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
454 << d << " secondes");
457 //destruction du tableau auxiliaire et de l'arbre
462 MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
470 // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
471 // et des points de la frontiere, des points internes imposes interieurs
472 // ==========================================================================
473 tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
474 mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
475 moartr, mxartr, n1artr, mnartr, mnarst,
478 // destruction de la queue et de l'arbre devenus inutiles
479 delete [] mnqueu; mnqueu=NULL;
480 delete [] mntree; mntree=NULL;
485 MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
487 // ierr =0 si pas d'erreur
488 // =1 si le tableau mnsoar est sature
489 // =2 si le tableau mnartr est sature
490 // =3 si aucun des triangles ne contient l'un des points internes
491 // =5 si saturation de la queue de parcours de l'arbre des te
492 if( ierr != 0 ) goto ERREUR;
494 //qualites de la triangulation actuelle
495 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
496 nbt, quamoy, quamin );
498 // boucle sur les aretes internes (non sur une ligne de la frontiere)
499 // avec echange des 2 diagonales afin de rendre la triangulation delaunay
500 // ======================================================================
501 // formation du chainage 6 des aretes internes a echanger eventuellement
502 aisoar( mosoar, mxsoar, mnsoar, na );
503 tedela( mnpxyd, mnarst,
504 mosoar, mxsoar, n1soar, mnsoar, na,
505 moartr, mxartr, n1artr, mnartr, n );
507 MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
510 MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
511 << d << " secondes");
513 //qualites de la triangulation actuelle
514 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
515 nbt, quamoy, quamin );
517 // detection des aretes frontalieres initiales perdues
518 // triangulation frontale pour les restaurer
519 // ===================================================
521 if( mn1arcf != NULL ) delete [] mn1arcf;
522 if( mnarcf != NULL ) delete [] mnarcf;
523 if( mnarcf1 != NULL ) delete [] mnarcf1;
524 if( mnarcf2 != NULL ) delete [] mnarcf2;
525 mn1arcf = new Z[1+mxarcf];
526 if( mn1arcf == NULL ) goto ERREUR;
527 mnarcf = new Z[3*mxarcf];
528 if( mnarcf == NULL ) goto ERREUR;
529 mnarcf1 = new Z[mxarcf];
530 if( mnarcf1 == NULL ) goto ERREUR;
531 mnarcf2 = new Z[mxarcf];
532 if( mnarcf2 == NULL ) goto ERREUR;
534 terefr( nbarpi, mnpxyd,
535 mosoar, mxsoar, n1soar, mnsoar,
536 moartr, mxartr, n1artr, mnartr, mnarst,
537 mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
540 MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere ierr=" << ierr );
543 MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
544 << d << " secondes");
546 if( ierr != 0 ) goto ERREUR;
548 //qualites de la triangulation actuelle
549 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
550 nbt, quamoy, quamin );
552 // fin de la triangulation avec respect des aretes initiales frontalieres
554 // suppression des triangles externes a la surface
555 // ===============================================
556 // recherche du dernier triangle utilise
557 mn = mxartr * moartr;
558 for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
561 if( mnartr[mn] != 0 ) break;
564 if( mntrsu != NULL ) delete [] mntrsu;
565 mntrsu = new Z[ndtri0];
566 if( mntrsu == NULL ) goto ERREUR;
568 if( mnlftr != NULL ) delete [] mnlftr;
569 mnlftr = new Z[nblf];
570 if( mnlftr == NULL ) goto ERREUR;
572 for (n=0; n<nblf; n++) //numero de la ligne fermee de 1 a nblf
575 tesuex( nblf, mnlftr,
576 ndtri0, nbsomm, mnpxyd, mnslig,
577 mosoar, mxsoar, mnsoar,
578 moartr, mxartr, n1artr, mnartr, mnarst,
581 delete [] mnlftr; mnlftr=NULL;
582 delete [] mntrsu; mntrsu=NULL;
586 MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
587 if( ierr != 0 ) goto ERREUR;
589 //qualites de la triangulation actuelle
590 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
591 nbt, quamoy, quamin );
593 // amelioration de la qualite de la triangulation par
594 // barycentrage des sommets internes a la triangulation
595 // suppression des aretes trop longues ou trop courtes
596 // modification de la topologie des groupes de triangles
597 // mise en delaunay de la triangulation
598 // =====================================================
599 mnarcf3 = new Z[mxarcf];
600 if( mnarcf3 == NULL )
602 MESSAGE ( "aptrte: MC saturee mnarcf3=" << mnarcf3 );
605 teamqt( nutysu, aretmx, airemx,
606 mnarst, mosoar, mxsoar, n1soar, mnsoar,
607 moartr, mxartr, n1artr, mnartr,
608 mxarcf, mnarcf2, mnarcf3,
609 mn1arcf, mnarcf, mnarcf1,
610 nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
612 if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
613 if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
614 if( mnarcf != NULL ) {delete [] mnarcf; mnarcf =NULL;}
615 if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
616 if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
620 MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
621 if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
622 if( ierr != 0 ) goto ERREUR;
624 //qualites de la triangulation finale
625 qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
626 nbt, quamoy, quamin );
628 // renumerotation des sommets internes: mnarst(i)=numero final du sommet
629 // ===================================
630 for (i=0; i<=nbsomm; i++)
633 for (nt=1; nt<=mxartr; nt++)
635 if( mnartr[nt*moartr-moartr] != 0 )
637 //le numero des 3 sommets du triangle nt
638 nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
639 //les 3 sommets du triangle sont actifs
640 mnarst[ nosotr[0] ] = 1;
641 mnarst[ nosotr[1] ] = 1;
642 mnarst[ nosotr[2] ] = 1;
646 for (i=1; i<=nbsomm; i++)
652 // generation du tableau uvst de la surface triangulee
653 // ---------------------------------------------------
654 if( uvst != NULL ) delete [] uvst;
656 if( uvst == NULL ) goto ERREUR;
659 for (i=0; i<nbsomm; i++ )
664 uvst[nbst].x = mnpxyd[i].x;
665 uvst[nbst].y = mnpxyd[i].y;
667 //si le sommet est un point ou appartient a une ligne
668 //ses coordonnees initiales sont restaurees
675 //retour aux coordonnees initiales dans uvslf
677 n = n - 1000000 * l + nudslf[l-1] - 1;
678 uvst[nbst].x = uvslf[n].x;
679 uvst[nbst].y = uvslf[n].y;
683 //point utilisateur n interne impose
684 //retour aux coordonnees initiales dans uvpti
685 uvst[nbst].x = uvpti[n-1].x;
686 uvst[nbst].y = uvpti[n-1].y;
693 // generation du tableau 'nsef' de la surface triangulee
694 // -----------------------------------------------------
695 // boucle sur les triangles occupes (internes et externes)
696 if( nust != NULL ) delete [] nust;
697 nust = new Z[nbsttria*nbt];
698 if( nust == NULL ) goto ERREUR;
700 for (i=1; i<=mxartr; i++)
702 //le triangle i de mnartr
703 if( mnartr[i*moartr-moartr] != 0 )
705 //le triangle i est interne => nosotr numero de ses 3 sommets
706 nusotr( i, mosoar, mnsoar, moartr, mnartr, nosotr );
707 nust[nbt++] = mnarst[ nosotr[0] ];
708 nust[nbt++] = mnarst[ nosotr[1] ];
709 nust[nbt++] = mnarst[ nosotr[2] ];
713 nbt /= nbsttria; //le nombre final de triangles de la surface
714 MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
715 << nbt << " triangles" );
718 MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
720 // destruction des tableaux auxiliaires
721 // ------------------------------------
723 if( mnarst != NULL ) delete [] mnarst;
724 if( mnartr != NULL ) delete [] mnartr;
725 if( mnslig != NULL ) delete [] mnslig;
726 if( mnsoar != NULL ) delete [] mnsoar;
727 if( mnpxyd != NULL ) delete [] mnpxyd;
728 if( mntree != NULL ) delete [] mntree;
729 if( mnqueu != NULL ) delete [] mnqueu;
730 if( mntrsu != NULL ) delete [] mntrsu;
731 if( mnlftr != NULL ) delete [] mnlftr;
732 if( mn1arcf != NULL ) delete [] mn1arcf;
733 if( mnarcf != NULL ) delete [] mnarcf;
734 if( mnarcf1 != NULL ) delete [] mnarcf1;
735 if( mnarcf2 != NULL ) delete [] mnarcf2;
736 if( mnarcf3 != NULL ) delete [] mnarcf3;
740 if( ierr == 51 || ierr == 52 )
742 //saturation des sommets => redepart avec 2 fois plus de sommets
749 MESSAGE( "APTRTE: Triangulation NON REALISEE avec erreur=" << ierr );
750 if( ierr == 0 ) ierr=1;
761 qualitetrte( R3 *mnpxyd,
762 Z & mosoar, Z & mxsoar, Z *mnsoar,
763 Z & moartr, Z & mxartr, Z *mnartr,
764 Z & nbtria, R & quamoy, R & quamin )
765 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
766 // but : calculer la qualite moyenne et minimale de la triangulation
767 // ----- actuelle definie par les tableaux mnsoar et mnartr
770 // mnpxyd : tableau des coordonnees 2d des points
771 // par point : x y distance_souhaitee
772 // mosoar : nombre maximal d'entiers par arete et
773 // indice dans mnsoar de l'arete suivante dans le hachage
774 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
775 // attention: mxsoar>3*mxsomm obligatoire!
776 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
777 // chainage des aretes frontalieres, chainage du hachage des aretes
778 // hachage des aretes = mnsoar(1)+mnsoar(2)*2
779 // avec mxsoar>=3*mxsomm
780 // une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
781 // mnsoar(2,arete vide)=l'arete vide qui precede
782 // mnsoar(3,arete vide)=l'arete vide qui suit
783 // moartr : nombre maximal d'entiers par arete du tableau mnartr
784 // mxartr : nombre maximal de triangles declarables
785 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
786 // arete1 = 0 si triangle vide => arete2 = triangle vide suivant
789 // nbtria : nombre de triangles internes au domaine
790 // quamoy : qualite moyenne des triangles actuels
791 // quamin : qualite minimale des triangles actuels
792 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
795 Z nosotr[3], mn, nbtrianeg, nt, ntqmin;
805 for ( nt=1; nt<=mxartr; nt++ )
810 //un triangle occupe de plus
813 //le numero des 3 sommets du triangle nt
814 nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
816 //la qualite du triangle ns1 ns2 ns3
817 qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
823 //la qualite minimale
824 if( qualite < quamin )
830 //aire signee du triangle nt
831 d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
834 //un triangle d'aire negative de plus
836 MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
837 << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
838 << " a une aire " << d <<"<=0");
841 //aire des triangles actuels
848 MESSAGE("Qualite moyenne=" << quamoy
849 << " Qualite minimale=" << quamin
850 << " des " << nbtria << " triangles de surface plane totale="
855 //le numero des 3 sommets du triangle ntqmin de qualite minimale
856 nusotr(ntqmin, mosoar, mnsoar, moartr, mnartr, nosotr );
857 MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
858 <<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
859 for (int i=0;i<3;i++)
860 MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
861 <<" y="<< mnpxyd[nosotr[i]-1].y);
865 MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );