Triangulation de Delaunay

Triangulation d’un robot
selon la méthode de Delaunay

Prince Mickaël 2017

GLOSSAIRE

Ce projet contient des notions de mathématiques et d’informatique qui nécessite d’inclure certaines définitions en amont afin de facilité la compréhension générale.
Algorithme itératif : une méthode itérative est un procédé algorithmique utilisé pour résoudre un problème, par exemple la recherche d’une solution d’un système d'équations ou d’un problème d’optimisation.
Algorithme récursif : Un algorithme récursif est un algorithme qui résout un problème en calculant des solutions d'instances plus petites du même problème
Cercle circonscrit à un triangle : étant donné un triangle non plat (i.e dont les sommets ne sont pas alignés), le cercle circonscrit à ce triangle est l'unique cercle passant par les sommets de ce triangle.
Complexité algorithmique : L'analyse de la complexité d'un algorithme consiste en l'étude formelle de la quantité de ressources (par exemple de temps ou d'espace) nécessaire à l'exécution de cet algorithme.
Convexité : En mathématiques, le mot « convexe » est utilisé dans la désignation de deux notions bien distinctes quoique apparentées :

  • lorsqu'il se rapporte à une forme géométrique, un ensemble de points, il renvoie au concept d'ensemble convexe ;
  •  lorsqu'il se rapporte à une fonction, il renvoie au concept de fonction convexe.

Dichotomie : procédé de recherche itératif par lequel on découpe en deux parties l'ensemble de recherche pour le restreindre à chaque étape.
diviser pour régner : Une technique algorithmique consistant à :

  • Diviser : découper un problème initial en sous-problèmes ;
  • Régner : résoudre les sous-problèmes (récursivement ou directement s'ils sont assez petits) ;
  • Combiner : calculer une solution au problème initial à partir des solutions des sousproblèmes.

Maillage : Un maillage est la discrétisation spatiale d'un milieu continu, ou aussi, une modélisation
géométrique d’un domaine par des éléments proportionnés finis et bien définis.
Points cocycliques : En géométrie, des points du plan sont dits cocycliques s'ils appartiennent à un
même cercle. Trois points non alignés du plan sont cocycliques. En effet, tout triangle possède
un cercle circonscrit.
Templates : en langage C, il s'agit d'un patron définissant un modèle de fonction ou de classe dont
certaines parties sont des paramètres (type traité, taille maximale).

INTRODUCTION

Ce projet est une réécriture de la Triangulation selon la méthode Delaunay, elle n’a pas pour but d’affirmer ou d’infirmer les dires de mes pairs mais de renouveler l’explication des méthodes et techniques utilisées pour parvenir au rendu final.
La triangulation selon Delaunay est un ensemble d’arêtes et points d’intersection – sous la forme de triangles - formant un maillage fini. Les points, placés aléatoirement sur le schéma, dans le but d’augmenter la densité des points dans la zone d’intérêt et ainsi d’améliorer la qualité de la
solution à apporter.
Il est à noter que la forme du maillage fini importe peu, il est donc nécessaire d’en prendre compte lors du traitement de la solution.
L’objectif de ce compte rendu est de développer comment un logiciel portable permettant de générer un maillage peut être utilisé pour un système de triangulation dans le but trouver l’emplacement d’un robot.
La première partie sera dirigé vers la présentation du projet, de ses aspect théoriques et de la méthode de travail pour arriver au rendu final ; puis dans un second temps seront détaillés les méthodes d’implémentations utilisé. Enfin, un exemple d’application sera développé pour visualisé le résultat et constituer une analyse précise du rendu.

I - PRÉSENTATION

I1 - La triangulation de Delaunay

Il existe de nombreuses méthodes pour appliquer une triangulation à partir d’un ensemble de points. Cependant, toutes ne sont pas exploitables. Pour cet exercice, il est nécessaire que l’ensemble des triangle dans le maillage soit relativement uniformes. Si ce n’est pas le cas, les calculs fournis seraient beaucoup moins précis.
LA grande contrainte de la méthode de Delaunay est : pour chaque triangle du maillage, son cercle circonscrit ne doit contenir aucun élément de l'ensemble des sommets comme le montre la figure suivante :

Représentation de la propriété de Delaunay

Cette triangulation possède des propriétés intéressantes pour construire un algorithme de génération de maillage [4]. On se restreint ici aux propriétés de la triangulation dans un espace à deux dimensions.
Soit n le nombre de points :

  • La triangulation est unique s'il n'y a pas trois points alignés ou quatre points cocycliques.
  • La triangulation contient au plus n simplexes.
  • L'union des simplexes de la triangulation forme l'enveloppe convexe de l'ensemble des points.

I2 - Algorithmes

Afin de générer une triangulation il est à la fois possible d’utiliser des algorithmesitératifs et récursifs [voir glossaire]. Dans les deux cas, on utilise un ensemble de points donnant unensemble de triangles (la triangulation).
L’algorithme itératif démarre d’un grand triangle à l’intérieur duquel va s’ajouter de nombreux triangles plus petits afin d’obtenir au final, le rendu désiré.
L’algorithme récursif quand à lui, adopte la stratégie du « diviser pour régner ». Une division par dichotomie de l'ensemble de points est effectuée avant fusion des ensembles.

I3 – Algorithme itératifs

L'ensemble de points étant fini dans le plan. On peut ainsi définir deux triangles initiaux englobant l'ensemble des points, Puis, on choisit un point dans l'ensemble initial (le
point rouge en image2). On peut localiser le triangle dans lequel il se trouve et créer les nouveaux triangles formés entre le point et les sommets de ce triangle.On continue de la sorte jusqu'à ne plus
avoir de point à choisir. On obtient la triangulation de Delaunay en rendu.


Première itération de l'algorithme incrémental

La complexité théorique de ce genre d'algorithme est quadratique.

I4 – Algorithme récursif

On démarre par un ensemble de points auquel on impose une ordonnance selon l’axe des abysses puis des ordonnées avant de séparer en deux sous-ensembles les points de manière strictement égale numériquement. On répète l’opération jusqu’à ce qu’il n’y ai plus que des ensemble de 3 points maximum.
Puis, on relie chaque groupe de points pour obtenir la triangulation de Delaunay relative à ces deux ensembles en supprimant certaines arêtes en cas de nécessité. On obtient ainsi un nouvel ensemble constitué de triangles. On réitère la fusion entre deux ensembles et cela jusqu'à ce qu'il n'y ait plus d'ensemble à fusionner.
Par construction, on obtient la triangulation de Delaunay de l'ensemble de points initial.


Division des ensembles puis fusion

II – DÉVELOPPEMENT

II1 – Structure du programme

Pour développer une triangulation nous devons disposer d’un ensemble de points repérés par leurs coordonnées sur le plans.
En entrée, notre programme prend :

  • Quatre limites de réels (Xmin,Ymin,Xmax,Ymax) qui représente l'intervalle I de |R2 dans lequel on veut générer la triangulation.
  • Une fonction densité sur cet intervalle. Il s'agit d'une fonction d : I → [0,1] qui représente la probabilité qu'un point de I existe en fonction de ses coordonnées. Cette densité aura une valeur proche de 1 sur les zones où l'on veut créer beaucoup de points (zone de grand intérêt) et une valeur proche de 0 dans une zone où le nombre de points est désiré faible.
  • Une fonction distance définie sur I. Il s'agit d'une fonction d' : I → |R qui représente la distance signée d'un point à un ensemble. Plus concrètement, cette fonction aide à définir une forme quelconque sur l'intervalle comme l'objet à mailler. Étant donné un point P de I, d'(P) sera positif si le P se situe en dehors du domaine défini par l'ensemble, négatif s'il se situe à l'intérieur et nul à sa frontière.

Le programme va alors quadriller l’intervalle et créer une grille sur laquelle seront généré l’ensemble des points qui servira de base pour la triangulation.

À cet ensemble de points sont retirés ceux qui ne respectent pas la densité. Ainsi, on obtient un nuage de points qui représente globalement la densité souhaitée.

Puis, l'algorithme de triangulation de Delaunay est appliqué sur ce nuage de points triés. On obtient ainsi un maillage sur tout l'intervalle I.

Enfin, le programme adapte ce maillage en utilisant la fonction distance d'. Pour chaque triangle du maillage, on calcule les distances des ses sommets avec la fonction d'. Ainsi, on sait quels sont ses sommets qui se trouvent à l'extérieur de l'ensemble. Si tous ses sommets sont à l'extérieur, cela implique que le triangle tout entier est à l'extérieur et on le supprime du maillage final. Dans le cas particulier où un triangle chevauche la frontière de l'ensemble, il faut projeter les sommets situés à l'extérieur sur la frontière (Illustration 4).

II2 Triangulation de Delaunay

Ci-dessous, le schéma général des classes dont la plupart hérite de classes appartenant à la libraire de la view.

La classe Mesh permet de stocker et de manipuler des maillages d'éléments finis. zzMesh est une classe qui hérite de Mesh que nous avons implémenté. Son rôle est de gérer toutes les actions relatives au maillage comme l'ajout et la suppression d'un triangle.
Grid est une classe qui permet de gérer une grille. zzCloud en hérite  et peut ainsi gérer l'ensemble des points d'une grille.
Element a pour rôle de gérer et de stocker élément fini. zzTriangle, classe fille de Element, est la classe qui représente un élément fini te type triangle.
Point2D définit un point dans le plan. zzVertex en hérite et comprend toutes les méthodes relatives à la gestion d'un point du plan.
Le fichier « delaunay.cpp » contient les fonctions relatives à la triangulation.
Comme nous l'avons vu précédemment, l'algorithme réalise d'abord une division récursive de l'ensemble des points avant de fusionner ces ensembles deux à deux [3].
Pour fusionner deux ensembles, l'algorithme n'a besoin de connaître que leur enveloppe convexe et l'ensemble des triangles. Le principe est de partir du segment le plus bas reliant les ensembles et de le faire remonter en construisant au fur et à mesure des triangles respectant la condition de Delaunay.
L'algorithme de fusion est détaillé dans le schéma ci-après (Schéma 2).

Entrées : E1, E2, les enveloppes convexes des ensembles à fusionner.
Mesh, l'ensemble des triangles existants.
Sortie : E, l'enveloppe convexe après fusion de E1 et E2.
Tant Que non Stop
Soit (BG,BD) la base inférieure de la fusion
Soit (HG,HD) la base supérieure de la fusion
Rechercher le candidat CG dans E1
Rechercher le candidat CD dans E2
Supprimer les triangles qui croisent (CG,BG,BD) et ceux qui croisent (CD,BG,BD)
Si CG = HG et CD = HD : Stop
Sinon:
Si CD est dans le cercle circonscrit à (CG,BG,BD) :
Ajouter (CD,BG,BD) à Mesh
BD ← CD
Sinon
Ajouter (CG,BG,BD) à Mesh
BG ← CG
Fin Tant Que
Schéma 2: Algorithme Fusion

Remarques : La base inférieure de la fusion est le couple de points (BG,BD) (pour Bas Gauche et Bas Droit) tels que tout point de E1 et tout point de E2 se trouvent au dessus de la droite (BG,BD) . Symétriquement, la base supérieure est le couple de points (HG,HD) (pour Haut Gauche et Haut Droit) tels que tout point de E1 et tout point de E2 se trouvent au dessous de la droite (HG,HD)

 Bases inférieure et supérieure de la fusion

Pour « monter », la base doit déterminer un candidat parmi les voisins de BG et ceux de BD. Cette détermination du candidat dans l'ensemble de gauche se fait symétriquement de la même manière que dans l'ensemble de droite. Le schéma suivant (Schéma 3) décrit l'algorithme de recherche de candidat pour l'ensemble de gauche, donc à partir de BG :

Entrées : BD, autre point de la base
Sortie : CG, candidat de l'ensemble de gauche
Initialisation : Angle_min ← 180°
Pour chaque voisin V de BG :
Si  BD , BG , V < Angle_min :
Angle_min ← BD , BG , V
CG = V
Fin Pour
Angle_min ← 180°
Pour chaque voisin V de BG :
Si  CG , BG ,V < Angle_min :
Angle_min ← CG , BG ,V
C2 ← V
Fin Pour
Si C2 appartient au cercle circonscrit à (CG,BG,BD), alors CG ← C2
Schéma 3: Algorithme de détermination du candidat dans l'ensemble de gauche

La base (BG,BD) évolue ainsi en remontant vers (HG,HD). Les nouveaux triangles créés sont ajoutés à la triangulation et les anciens triangles chevauchant les nouveaux sont supprimés. Il faut aussi noter que les points situés sur les faces en regard des deux ensembles à fusionner ne seront pas ajoutés à l'enveloppe convexe de l'ensemble final.


Fusion des deux ensembles

Il reste à adapter le maillage à la forme définie par la fonction distance.
L'algorithme d'adaptation choisi est le suivant (Schéma 4):

Pour chaque triangle (A,B,C) de T
Si les 3 sommets sont à l'extérieur : Supprimer (A,B,C) de T.
Si 2 sommets A et B sont à l'extérieur : Soit PA et PB les projections respectives sur la frontière de A et B selon ->AC et -> BC . Supprimer (A,B,C) à T. Ajouter (PA,PB,C) à T.
Si le sommet A seulement est à l'extérieur : Soit PA1 et PA2 les projections sur la frontière de A selon respectivement -> AB et ->AC . Supprimer (A,B,C) à T. Ajouter (PA1,B,C) et (PA2,B,C) à T.
Sinon : Ne rien faire.
Fin Pour

On voit ici que la difficulté repose sur le chevauchement des triangles avec la frontière.
Dans ce cas, les deux situations à envisager (deux point à l'extérieur ou un unique à l'extérieur) utilisent une projection selon une direction (Illustration 13, Illustration 14). Comme on veut que le point extérieur soit projeté sur la frontière, en notant A le point à l'extérieur et B le point à l'intérieur, la direction de projection est ->AB . Voici le schéma de principe de la projection.

Entrée : Point A externe au domaine. Point B interne au domaine.
Sortie : Point A projeté sur la frontière
Faire
Soit I le milieu de [A, B].
Si d(I)>0 : A ← I
Sinon : B ← I
Tant que |d(I)| > ε

On remarque qu'il s'agit d'une recherche d'un point sur un segment. Ce type de recherche est rapide et précise. L'epsilon choisi dans le programme est de l'ordre de 10 −1 .
Il est important de noter que les triangles modifiés sur la frontière peuvent ne plus respecter la condition de Delaunay. Cependant, cette restructuration ne peut pas altérer la qualité de la triangulation car cela concentre les points sur la frontière qui est une zone de grand intérêt.


III RÉSULTATS

III1 Triangulation d’un robot

Nous allons maintenant utiliser le programme dans une étude concrète. L'objet considéré est une pièce robotique circulaire et on désire étudier sa déformation sous l'effet d'une force appliquée sur une zone de son contour.

On se restreint à une étude de la déformation dans le plan de la pièce robotique. Il nous faut tout d'abord définir un intervalle du plan. Prenons le pavé I=[0,100]×[0,100] (Illustration 16). Le quadruplet (Xmin,Ymin,Xmax,Ymax) est donc (0,0,100,100).
Dans le repère créé, la pièce robotique a pour centre le point (50,50) et pour rayon 40.
Dans cet intervalle, compte tenu de la géométrie de la pièce robotique, il est aisé de définir la fonction distance qui la représente. Ici :
Pour tout (x,y) de I, d'(x,y) = (x50)2 + (y50)2 – 402 = (x50)2 + (y50)2 – 1600Il faut aussi choisir le nombre de points que la grille initiale comportera, par exemple : 50 points sur l'axe des abscisses et 50 points sur l'axe des ordonnées.Il nous faut définir une fonction densité pour la génération des points du maillage. Comme la zone d'intérêt se situe autour du point d'application de la force (90,50), la densité doit être plus importante dans cette région. Prenons la fonction suivante :
Pour tout (x,y) de I, d(x,y) =  x −10 2 y−50 2/10000 Autour du point d'intérêt (90,50), la valeur de d est très proche de 1 et plus on s'éloigne de ce point, plus la densité décroit. Une division par une grande valeur est nécessaire pour que les valeurs soient comprises entre 0 et 1. Ici, cette valeur a été déterminée empiriquement pour avoir une bonne concentration de points autour de la zone d'intérêt et peu de points ailleurs (à gauche de la pièce robotique).Les paramètres étant déterminés, on peut les entrer dans le programme avant de le lancer. On peut observer le résultat obtenu.

Affichage du résultat

À partir de la liste des triangles générée, l'utilisateur peut appliquer une méthode de résolution par éléments finis pour calculer la déformation causée par la force F sur la pièce robotique.

III2 Analyses

Il est intéressant d'analyser la complexité du programme pour évaluer son efficacité.
Le programme mêle des fonctions itératives et récursives et il devient très vite compliqué d'évaluer avec précision sa complexité théorique. Cependant une rapide analyse de la structure du programme permet d'en donner une estimation.
La génération de la grille est fonction de la densité et du nombre de points initiaux voulus. Avec m le nombre de points sur l'axe des abscisses et n le nombre de points sur l'axe des ordonnées, la génération de la grille est de l'ordre de mxn calculs.
Soit N le nombre de points restants après génération de la grille. La division de Delaunay se fait de manière récursive et admet une complexité théorique en N×log(N).La fusion de deux ensembles convexes, appelée à chaque étape de la division, requiert :

  • Une recherche des bases inférieure et supérieure entre les ensembles : avec p et q les nombres respectifs de points du premier et du second ensemble convexe, la complexité est de l'ordre de p2 + q2.
  • Pour chaque point considéré, il faut rechercher parmi ses voisins le candidat suivant pour déterminer la nouvelle base. Cela et linéaire en fonction du nombre de voisins. Cette recherche est effectuée 4 fois par étape.
  • Lorsqu'un candidat est trouvé et qu'un triangle est ajouté, il faut chercher parmi tous les triangles déjà construits s'il n'y a pas chevauchement de triangle. Cela est linéaire en fonction du nombre de triangles existants.
  • La création de l'ensemble convexe pour l'étape suivante se fait par réunion des deux ensembles convexes traités. Avec p et q les nombres respectifs de points du premier et du
    p/2+q/2

second ensemble convexe, le nombre de calculs peut être estimé de l'ordre de L'adaptation de la triangulation à une forme parcourt tous les triangles une fois. Pour chaque triangle,
dans les cas de chevauchement avec la frontière de la forme, une recherche dichotomique du(des) point(s) d'intersection s'effectue entre le point externe A et le point interne B. Avec k le rapport entre la longueur du segment [A,B] et l'epsilon choisi pour la précision de la distance à la frontière, cette recherche est de l'ordre de k×log(k).
Finalement, trop de paramètres indéterminés sont à considérer pour le calcul exact de la complexité comme le nombre moyen de voisins d'un point, le nombre de points moyens
contenus dans l'enveloppe convexe d'un ensemble, etc.
Cependant, une mesure du temps d'exécution en fonction de nombre de point entrés initialement peut indiquer comment se comporte le programme quand le nombre de points augmente.

Cette mesure a été faite pour un maillage sur l'intervalle [0,100]×[0,100], avec une densité
d = x−50 2 y−50 2/10000
Les temps d'exécution variant considérablement entre deux mêmes lancements du programme et sur deux ordinateurs différents avec les mêmes paramètres, les valeurs obtenues ne sont pas significatives. Seule l'allure de la courbe est représentative du fonctionnement du programme.
On observe ici que la complexité réelle n'atteint pas la complexité théorique en n×log(n) prévue par l'algorithme de Delaunay et le comportement du programme avec un nombre de point supérieur à 100000 devient légèrement longue. Cependant, le résultat n'est pas insatisfaisant, comptetenu de la difficulté à adapter les structures de données offertes par un langage de programmation tel que le C++ à un algorithme tel que celuici.
Certaines parties de notre programme peuvent être optimisées en adoptant une autre stratégie.Par exemple, la fonction checkIntersect(A,B) qui parcourt tous les triangles existant pour supprimer ceux qui possèdent un côté qui croise le segment [A,B], est très gourmande en calculs. Cette fonction aurait pu être évitée en mettant à jour un voisinage de triangles pour chaque triangle. Seul le parcourt de la liste des voisins aurait été utile. Le fait que les points soient générés sur une grille induit l'apparition de triangles plats. Ces triangles particuliers sont autant de sous cas à traiter indépendamment des autres. Par exemple, le test d'appartenance au cercle circonscrit ne peut pas être appliqué comme habituellement puisque le cercle n'est pas défini. C'est une des raisons pour lesquelles la fonction de recherche de candidat lors de la fusion arbore une structure si particulière. En effet, beaucoup de tests sont à effectuer pour tester l'alignement de sommets, et de calculs de produits scalaires dans ces cas précis. Cela a pour conséquence une complexification considérable du code et cela le rend difficile à interpréter. Une clarification de cette fonction est envisageable pour augmenter sa lisibilité et son efficacité.La première implémentation du programme n'utilisait pas la librairie OFELI. L'incorporation de nos classes dans les classes préexistantes de la librairie nous a posé quelques problèmes, notamment lors de la généralisation des types utilisés. En effet, les points étaient repérés par des coordonnées de type 'double'. Or la librairie OFELI utilise des templates pour une plus grande généralisation. Nous avons donc restructuré l'ensemble de notre code.

REMERCIEMENTS

Je tiens à remercier mes formateurs et responsables pour leur soutient tout au long de l’élaboration de ce rendu, leurs conseils et leurs flexibilité qui m’auront permis d’étendre mes connaissances et compétences.

SOURCES

https://fr.wikipedia.org/wiki/M%C3%A9thode_it%C3%A9rative
https://fr.wikipedia.org/wiki/Algorithme_r%C3%A9cursif
https://fr.wikipedia.org/wiki/Analyse_de_la_complexit%C3%A9_des_algorithmes
https://fr.wikipedia.org/wiki/Convexit%C3%A9
https://fr.wikipedia.org/wiki/Diviser_pour_r%C3%A9gner_(informatique)
https://fr.wikipedia.org/wiki/Maillage
https://fr.wikipedia.org/wiki/Cocyclique
https://fr.ISMA.fr

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *