triangle

La librairie "triangle" est une librairie python (également disponible dans d'autres langages) permettant de calculer des triangulations de Delaunay en imposant diverses contraintes. Cette librairie a été développée par J.R. Shewchuk de l'université de Californie, Berkeley. La documentation complète peut être trouvée sur le site de l'auteur. Plusieurs papiers donnent également des spécifications sur les algorithmes utilisés. Une liste non exaustive est donnée en fin de page.

Sommaire

1 - Installation

L'installation via l'invité de commande se fait comme pour toute autre librairie python. Lancez l'invité de commande windows et tapez l'instruction suivante. De plus amples instructions sont disponibles sur la page sur PyPi pour d'autres méthodes d'installation.


pip install triangle
						

2 - Utilisation simple

Comme pour tout autre librairie python, il est nécessaire de l'importer dans vos scripts.


import triangle as tr
						

2.1 - Structure des données d'entrée

La structure de données d'entrée minimale est la suivante :


init=dict (vertices= [(XA,YA),(XB,YB),(XC,YC),(XD,YD),...])
						

Avec :

  • (XA, YA), (XB, YB), ..., les coordonnées x et y des noeuds du maillage

2.2 - Triangulation d'un ensemble de points

Un maillage se génère grâce à la commande :


mesh = tr.triangulate(init)
						

La présente triangulation est une triangulation de Delaunay. Une telle triangulation (Delaunay, 1934) satisfait les conditions suivantes. Pour tout triangle, il existe un cercle circonscrit à ce triangle tel qu’aucun nœud ne se trouve à l’intérieur de ce cercle. Pour un ensemble quelconque de nœuds, cette triangulation existe et est unique à l’exception des cas où l’on a un ensemble de quatre nœuds cocycliques (par exemple formant un carré, lequel peut être triangulé selon l’une ou l’autre de ses diagonales).

2.3 - Structure des données de sortie

Les données de sortie sont stockées sous la forme d'un dictionnaire contenant deux array (structure de données de la librairie numpy.


mesh["vertices"]  #Noeuds de la triangulation
mesh["triangles"] #Description des triangles
						

Les noeuds de la triangulation sont identifiés par leur coordonnées (stockés sous forme d'array). Les triangles sont décrit par les identifiants des noeuds les composants. Ces identifiants correspondent à la position des noeuds dans mesh["vertices"]
Attention : ces identifiants commencent à 0, comme il est d'usage dans python.

2.4 - Exemple commenté


#Importation de la librairie triangle
import triangle as tr 

#Importation de la librairie random (pour générer des coordonnées aléatoires)
import random as rd 

#Importation de la librairie matplotlib pour l'affichage des résultats
import matplotlib.pyplot as plt 

#Génération aléatoire de 10 points aux coordonnées entières
#comprises entre 0 et 25
vert=[[rd.randint(0,25),rd.randint(0,25)] for i in range(10)]

#Triangulation de Delaunay
init=dict (vertices= vert)
mesh=tr.triangulate(init)

#Affichage des résultats
Figure = plt.figure()
Axes = Figure.add_subplot(111)
for tri in mesh["triangles"] : 
	for j in range(len(tri)):
		XA=mesh["vertices"][tri[j-1]][0] #Coordonnée X de la première extrémité
		XB=mesh["vertices"][tri[j]][0]   #Coordonnée X de la seconde extrémité
		YA=mesh["vertices"][tri[j-1]][1] #Coordonnée Y de la première extrémité
		YB=mesh["vertices"][tri[j]][1]   #Coordonnée Y de la seconde extrémité
		Axes.plot([XA,XB],[YA,YB],c="b")
plt.show()
#NB : ici certains arcs sont représentés deux fois s'ils appartiennent à deux triangles. 
#Cette double représentation à été laissée afin de garder ce premier exemple simple.
						

3 - Triangulation sous contraintes et paramètres de triangulation

3.1 - Qu'est ce qu'une triangulation sous contraintes ?

Une triangulation sous contraintes peut être de différentes natures. Le terme recouvre le plus souvent les triangulations pour lesquels on impose que certains segments reliant deux noeuds du domaines soient correspondent après triangulation à un côté d'un triangle du maillage. Par exemple, on peut souhaiter imposer des bordures au domaine.
On peut aussi vouloir imposer des contraintes d'angles minimaux dans les triangles ou des contraintes d'aires minimales. La librairie triangle permet de gérer cela en spécifiant le type de triangulation.


mesh = tr.triangulate(init, 'type')
						

'type', spécifie le type de triangulation et les contraintes sur celle-ci. Certaines de ces modalités sont détaillées dans les paragraphes suivants. La documentation officielle les détaillent toutes.

3.2 - Spécifier une bordure du domaine de triangulation

Pour spécifier les bordures du domaine, il suffit de renseigner, en plus de la liste des noeuds, une liste de segments dans la géométrie initiale. Cette liste se structure sous la forme d'un ensemble de couples de noeuds identifiés par leur position dans la liste des noeuds.


vert=[(XA,YA),(XB,YB),(XC,YC),(XD,YD),...]
borderExt = [(Indice_1,Indice_2),(Indice_2, Indice_3),...]
init=dict (vertices=vert , segments= [borderExt])
mesh = tr.triangulate(init)
						

La prise en compte de la bordure dans la triangulation se fait automatiquement dès lors que celle-ci forme un polygone convexe. Dans le cas contraire, il est nécéssaire d'utiliser le 'p' Switch (voir ci-dessous).

3.3 - Spécifier un trou et 'p' Switch

Le paramètre 'p', permet de prendre en compte les segments, les concavités et les espaces vides lors de la triangulation. Par défaut l'algorithme ne considère que la bordure extérieur.


vert=[(XA,YA),(XB,YB),(XC,YC),(XD,YD),...]
border = [(Indice_1,Indice_2),(Indice_2, Indice_3),...] #Bordure extérieur
border.extend((Indice_1,Indice_2),(Indice_2, Indice_3),...) #Ajout des bordures ou segments internes
trous = [[XH,YH]]
init=dict (vertices=vert , segments= border, holes=trous)
mesh = tr.triangulate(init, 'p')
						
  • border, les différents segments composant la bordure extérieure du polygone. Sont ajoutés dans un second temps, les différents segments composant une bordure intérieure du polygone. Chaque segment est identifié par les indices des noeuds constituant ses extrémités.
  • trous, la liste des points à l'intérieur des polygones à évider. Ces points doivent se trouver à l'intérieur d'un polygone décrit par une bordure intérieure. La triangulation ne sera ainsi pas appliqué à l'intérieur de ce polygone.

Attention : une telle triangulation ne suit généralement plus les propriétés de Delaunay. Celles-ci peuvent être vérifiées localement (pour un ou plusieurs triangles particulier), mais il est rare que ce soit le cas de manière globale. Pour que ce soit le cas, il est nécessaire d'utiliser en complèment le 'D' switch.

3.4 - Contraindre les angles et 'q' Switch

Le paramètre 'q', permet de contraindre les triangles à avoir des angles infèrieurs à 20°.
Un angle alternatif minimum peut être spécifié de la manière suivante :


mesh = tr.triangulate(poly, 'q14.3') #imposera un angle minimum de 14,3°
						

Si l'angle est inférieur à 20,7°, l'algorithme de triangulation doit normalement se terminer si on suppose une précision arithmetique infinie. En réalité, si on dépasse la précision maximum, cela peut causer une interruption de l'algorithme et la triangulation peut échouer. En pratique, l'algorithme fonctionne bien pour des angles supérieurs minimaux inférieur à 33°. Le taux d'échec est beaucoup plus important au delà.
Attention : python ne retournera pas nécessairement d'erreur. Si l'algorithme s'arrête de manière innopinée, cet aspect d'angle minimum peut en être la cause.

Dans les exemples ci-dessus, les contraintes d'angles et de bordures sont cumulées. Pour spécifier plusieurs options de triangulation, il suffit de noter les différents switchs les uns derrière les autres.


mesh = tr.triangulate(poly, 'pq') #Options de tringulation utilisées pour la figure 10
mesh = tr.triangulate(poly, 'pq30') #Options de tringulation utilisées pour la figure 11
mesh = tr.triangulate(poly, 'pq33') #Options de tringulation utilisées pour la figure 12
						

3.5 - Contraindre les aires et 'a' Switch

Le paramètre 'aX', permet de définir l'aire maximale des triangles : aucun triangle de généré n'aura une aire dépassant la valeur de X.


mesh = tr.triangulate(poly, 'a0.1')
						

3.6 - Triangulation Delaunay-contrainte et 'D' Switch

Le paramètre 'D' permet de s'assurer que la triangulation est bien Delaunay-conforme, c'est à dire qu'elle possède toutes les propriétés d'une triangulation de Delaunay. L'utilisation de ce paramètre autorise l'ajout de noeuds supplémentaires à la triangulation (là où l'algorithme de base réalise la triangulation uniquement des noeuds passés en entrée). Il est intéressant d'utiliser ce paramètre en complément du paramètre 'p'.


mesh = tr.triangulate(poly, 'D')
						

4 - Autres paramètres

4.1 - 'V' Switch

'V', permet de retourner des messages lors du calcul du maillage. Il est possible d'ajouter plusieurs 'V' afin de retourner des messages plus précis dans la console. Ces messages donnent des retours sur l'avancement des grandes étapes du calcul du maillage et sur les éventuelles erreur. Cette option est très utile pour le debugguage.


mesh = tr.triangulate(poly, 'VV')  #Permet d'afficher 2 niveaux de messages.
						

5 - Quelques spécifications théoriques

Cette section est conscrée à quelques explications techniques sur les algorithmes à l'oeuvre dans la librairie. L'ensemble de ces spécifications sont issues des articles écrits par l'auteur de la librairie ou des textes cités dans la partie références.

5.1 - Implémentation des algorithmes

Les algorithmes de triangulation sont implémentés de sorte à éliminer les noeuds dupliqués passés en entrée. Ils sont également en mesure de détecter les intersections de segments (bordures extérieures ou interieures). Dans ce cas, un nouveau noeud est créé à l'intersection de ces segments et les segments sont divisés en plusieurs parties dont le nouveau est une extrémité.

5.2 - Prise en compte des contraintes de segments

La prise en compte des contraintes fonctionne de deux manières différentes. suivante.

Première méthode : d'abord une triangulation de Delaunay de l'ensemble des points passés en entrée est réalisée. Ensuite, on identifie par mis les segments contraints ceux ne faisant pas partie d’un côté d’un triangle. Pour ces segments un point au milieu est inséré. L'algorithme d'insertion incrémental de Lawson est utilisé pour cela. Cet algorithme permet de maintenir les propriétés d'une triangulation de Delaunay. De cette manière, le segment est divisé par la moitié et les deux sous segments peuvent apparaître dans le maillage chacun comme un côté d'un triangle. Si ce n’est pas le cas la procédure est itérée sur les sous-segments jusqu’à tant que le segment original soit représenté par une séquence de segments appartenant au maillage.

La deuxième manière de faire est de réaliser d’emblée une triangulation de Delaunay sous contrainte avec un algorithme approprié. Chaque segment est inséré en supprimant les triangles qu’il recouvre et l’ensemble de la région aux alentours du segment est re-triangulée. La librairie Triangle utilise cette méthode par défaut.

La prise en compte des trous, se fait après coup. L'ensemble de la zone est donc triangulée en prenant en compte les contraintes de bordures interieurs et extérieurs. Dans un second temps, les triangles se situant dans les zones identifiées commes des trous sont supprimés.

5.3 - Prise en compte des contraintes d'angle

Les contraintes d'angles sont assurées une fois les étapes précédentes réalisées par insertion de noeuds intermédiaires grâce à l'algorithme de Lawson. L’insertion de nœuds est gouvernée par deux règles.

  • Propriété du cercle empiétant : un segment entre deux noeuds est dit empiétant si un troisième point au moins se trouve dans le cercle de diamètre ce segment. Chaque segment empiétant est immédiatement divisé par l’insertion d’un nouveau nœud en son milieu. Les deux segments résultant ont un cercle dont le segment-diamètre est plus petit est potentiellement n’est plus empiétant.
  • Propriété du cercle circonscrit : le cercle circonscrit au triangle est l’unique cercle qui passe par les trois nœuds du triangle. Un triangle est dit mauvais s’il a un angle trop petit ou une surface trop importante pour satisfaire les contraintes renseignées par l’utilisateur. Dans ce cas, le triangle est divisé par l’insertion d’un nouveau nœud au centre du cercle.

Il y a priorité de la propriété d’empiétement sur celle des mauvais triangles. La procédure de raffinement des mauvais triangles (selon leur angle) fonctionne pour une contrainte d’angle allant jusqu’à 20,7°. En pratique, selon Shewchuk, l’algorithme réussi l’opération la plupart du temps pour des contraintes allant jusqu’à 33,8°, mais échoue à partir de 33.9°. Les triangles situés dans les trous et les concavités sont supprimés avant l’étape de raffinement.

6 - Notes et références

6.1 - Voir aussi

6.2 - Références

  • Shewchuk, J. (2002a). Constrained Delaunay Tetrahedralizations and Provably Good Boundary Recovery. Proceedings of the 11th International Meshing Roundtable.
  • Shewchuk, J. (2002b). What Is a Good Linear Finite Element ? - Interpolation, Conditioning, Anisotropy, and Quality Measures. Proceedings of the 11th International Meshing Roundtable, 73.
  • Shewchuk, J. R. (1996). Triangle : Engineering a 2D Quality Mesh Generator and Delaunay Triangulator. In M. C. Lin & D. Manocha (Éds.), Applied Computational Geometry : Towards Geometric Engineering (Vol. 1148, p. 203‑222). Springer-Verlag.
  • Shewchuk, J. R. (1998). A condition guaranteeing the existence of higher-dimensional constrained Delaunay triangulations. Proceedings of the Fourteenth Annual Symposium on Computational Geometry  - SCG ’98, 76‑85. https://doi.org/10.1145/276884.276893