Au cours des dernières semaines, j'ai travaillé sur l'implémentation de l'algorithme de Fortune en C ++. Cet algorithme prend beaucoup de points 2D et construit à partir d'eux un diagramme de Voronoi . Si vous ne savez pas ce qu'est un diagramme de Voronoi, regardez la figure:


Pour chaque point d'entrée, appelé «site», nous devons trouver de nombreux points plus proches de cet endroit que de tout le monde. Ces ensembles de points forment les cellules montrées dans l'image ci-dessus.

Il est remarquable dans l'algorithme de Fortune qu'il construit de tels diagrammes dans le temps $$$O (n \ log n)$(ce qui est optimal pour un algorithme de comparaison), où $$$n$Est le nombre de places.

J'écris cet article car je considère que la mise en œuvre de cet algorithme est une tâche très difficile. Pour le moment, c'est le plus compliqué des algorithmes que j'ai dû implémenter. Par conséquent, je veux partager les problèmes que j'ai rencontrés et comment les résoudre.

Comme d'habitude, le code est affiché sur github , et tous les documents de référence que j'ai utilisés sont répertoriés à la fin de l'article.

Description de l'algorithme Fortune

Je ne vais pas expliquer comment fonctionne l'algorithme, car d'autres personnes l'ont déjà bien fait. Je peux recommander d'étudier ces deux articles: ici et ici . La seconde est très intéressante - l'auteur a écrit une démo interactive en Javascript, utile pour comprendre le fonctionnement de l'algorithme. Si vous avez besoin d'une approche plus formelle avec toutes les preuves, je vous recommande de lire le chapitre 7 de Computational Geometry, 3e édition .

De plus, je préfère traiter des détails d'implémentation qui ne sont pas bien documentés. Et ce sont eux qui rendent la mise en œuvre de l'algorithme si complexe. Je me concentrerai en particulier sur les structures de données utilisées.

Je viens d'écrire un pseudo-code de l'algorithme pour vous faire une idée de sa structure globale:

  ajouter un événement de lieu à la file d'attente d'événements pour chaque lieu
 jusqu'à ce que la file d'attente d'événements soit vide
     récupérer le meilleur événement
     si l'événement est un événement de lieu
         insérer un nouvel arc dans le littoral
         vérifier les nouveaux événements du cercle
     sinon
         créer un sommet dans le diagramme
         on retire du littoral un arc resserré
         supprimer les événements invalides
         vérifier les nouveaux événements du cercle 

Structure des données du graphique

Le premier problème que j'ai rencontré a été de choisir la façon de stocker le diagramme de Voronoi.

J'ai décidé d'utiliser une structure de données largement utilisée en géométrie numérique appelée la liste des bords doublement connectés (DCEL).

Ma classe VoronoiDiagram utilise quatre conteneurs comme champs:

 class VoronoiDiagram { public: // ... private: std::vector<Site> mSites; std::vector<Face> mFaces; std::list<Vertex> mVertices; std::list<HalfEdge> mHalfEdges; } 

Je vais parler en détail de chacun d'eux.

La classe Site décrit le point d'entrée. Chaque lieu a un index, qui est utile pour le mettre dans la file d'attente, des coordonnées et un pointeur vers la cellule ( face ):

 struct Site { std::size_t index; Vector2 point; Face* face; }; 

Les sommets de la cellule sont représentés par la classe Vertex , ils n'ont qu'un champ de coordonnées:

 struct Vertex { Vector2 point; }; 

Voici l'implémentation des demi-bords:

 struct HalfEdge { Vertex* origin; Vertex* destination; HalfEdge* twin; Face* incidentFace; HalfEdge* prev; HalfEdge* next; }; 

Vous vous demandez peut-être ce qu'est une demi-côte? Un bord dans le diagramme de Voronoi est commun à deux cellules voisines. Dans la structure de données DCEL, nous divisons ces bords en deux demi-bords, un pour chaque cellule, et ils sont liés par un pointeur twin . De plus, le demi-bord a un point de départ et un point d'arrivée. Le champ incidentFace indique la face à laquelle appartient le demi-bord. Les cellules dans DCEL sont implémentées en tant que liste cyclique de demi-bords liés, dans laquelle les demi-bords adjacents sont connectés ensemble. Par conséquent, les champs précédent et suivant indiquent les demi-bords précédent et suivant dans la cellule.

La figure ci-dessous montre tous ces champs pour le demi-bord rouge:


Enfin, la classe Face définit la cellule. Il contient simplement un pointeur sur sa place et un autre sur l'une de ses demi-côtes. Peu importe lequel des demi-bords est sélectionné, car la cellule est un polygone fermé. Ainsi, nous avons accès à toutes les demi-arêtes tout en parcourant une liste chaînée cyclique.

 struct Face { Site* site; HalfEdge* outerComponent; }; 

File d'attente des événements

La manière standard d'implémenter une file d'attente d'événements est avec une file d'attente prioritaire. Lors du traitement des événements de lieu et de cercle, il se peut que nous devions supprimer les événements de cercle de la file d'attente car ils ne sont plus valides. Mais la plupart des implémentations de file d'attente prioritaire standard ne vous permettent pas de supprimer un élément qui n'est pas en haut. Cela s'applique en particulier à std::priority_queue .

Il existe deux façons de résoudre ce problème. La première, la plus simple, consiste à ajouter un indicateur valid aux événements. valid est initialement défini sur true . Ensuite, au lieu de supprimer l'événement Circle de la file d'attente, nous pouvons simplement définir son indicateur sur false . Enfin, lors du traitement de tous les événements de la boucle principale, nous vérifions si la valeur de drapeau valid de l'événement est false , et si c'est le cas, sautez-la simplement et traitez la suivante.

La deuxième méthode que j'ai appliquée n'était pas d'utiliser std::priority_queue . Au lieu de cela, j'ai implémenté ma propre file d'attente prioritaire, qui prend en charge la suppression de tout élément qu'elle contient. La mise en place d'une telle file d'attente est assez simple. J'ai choisi cette méthode car elle rend le code de l'algorithme plus propre.

Littoral

La structure des données du littoral est une partie complexe de l'algorithme. En cas d'implémentation incorrecte, rien ne garantit que l'algorithme sera exécuté en $$$O (n \ log n)$. La clé pour obtenir cette complexité temporelle est d'utiliser un arbre à équilibrage automatique. Mais c'est plus facile à dire qu'à faire!

La plupart des ressources que j'ai étudiées (les deux articles ci-dessus et le livre sur la géométrie computationnelle ) vous conseillent d'implémenter le littoral comme un arbre dans lequel les nœuds internes indiquent les points de rupture et les feuilles indiquent les arcs. Mais ils ne disent rien sur la façon d'équilibrer un arbre. Je pense qu'un tel modèle n'est pas le meilleur, et voici pourquoi:

• il y a des informations redondantes: nous savons qu'il y a un point de rupture entre deux arcs adjacents, il n'est donc pas nécessaire de représenter ces points comme des nœuds
• il est inadéquat pour l'auto-équilibrage: seul le sous-arbre formé par les points de rupture peut être équilibré. Nous ne pouvons vraiment pas équilibrer tout l'arbre, sinon les arcs peuvent devenir des nœuds internes et des feuilles de points de rupture. Écrire un algorithme pour équilibrer uniquement le sous-arbre formé par les nœuds internes me semble un cauchemar.

J'ai donc décidé de présenter le littoral différemment. Dans mon implémentation, le littoral est également un arbre, mais tous les nœuds sont des arcs. Un tel modèle ne présente aucun des inconvénients énumérés.

Voici la définition de l'arc Arc dans mon implémentation:

 struct Arc { enum class Color{RED, BLACK}; // Hierarchy Arc* parent; Arc* left; Arc* right; // Diagram VoronoiDiagram::Site* site; VoronoiDiagram::HalfEdge* leftHalfEdge; VoronoiDiagram::HalfEdge* rightHalfEdge; Event* event; // Optimizations Arc* prev; Arc* next; // Only for balancing Color color; }; 

Les trois premiers champs sont utilisés pour structurer l'arbre. Le champ leftHalfEdge indique la demi-arête dessinée par le point le plus à gauche de l'arc. Et rightHalfEdge est sur le demi-bord dessiné par le point d'extrême droite. Deux pointeurs, précédent et suivant next sont utilisés pour accéder directement à l'arc précédent et suivant du littoral. En outre, ils vous permettent également de contourner le littoral en tant que liste doublement liée. Enfin, chaque arc a une couleur qui sert à équilibrer le littoral.

Pour équilibrer le littoral, j'ai décidé d'utiliser un schéma rouge-noir . Lors de l'écriture de code, je me suis inspiré du livre Introduction to Algorithms . Le chapitre 13 décrit deux algorithmes intéressants, insertFixup et deleteFixup , qui équilibrent l'arbre après l'insertion ou la suppression.

Cependant, je ne peux pas utiliser la méthode d' insert indiquée dans le livre, car les clés sont utilisées pour trouver le point d'insertion du nœud. Il n'y a pas de clés dans l'algorithme de Fortune, nous savons seulement que nous devons insérer un arc avant ou après un autre dans le littoral. Pour implémenter cela, j'ai créé les insertAfter insertBefore et insertAfter :

 void Beachline::insertBefore(Arc* x, Arc* y) { // Find the right place if (isNil(x->left)) { x->left = y; y->parent = x; } else { x->prev->right = y; y->parent = x->prev; } // Set the pointers y->prev = x->prev; if (!isNil(y->prev)) y->prev->next = y; y->next = x; x->prev = y; // Balance the tree insertFixup(y); } 

L'insertion de y avant x se fait en trois étapes:

1. Trouvez un endroit pour insérer un nouveau nœud. Pour ce faire, j'ai utilisé l'observation suivante: l'enfant gauche x ou l'enfant droit x->prev est Nil , et celui qui est Nil est avant x et après x->prev .
2. À l'intérieur du littoral, nous conservons la structure d'une liste doublement liée, nous devons donc mettre à jour en conséquence les pointeurs prev et next des éléments x->prev , y et x .
3. Enfin, nous appelons simplement la méthode insertFixup décrite dans le livre pour équilibrer l'arbre.

insertAfter est implémenté de manière similaire.

La méthode de suppression extraite du livre peut être implémentée sans modifications.

Limite du graphique

Voici la sortie de l'algorithme de Fortune décrit ci-dessus:


Il y a un petit problème avec certains bords de cellules sur le bord de l'image: ils ne sont pas dessinés car ils sont infinis.

Pire, une cellule peut ne pas être un fragment unique. Par exemple, si nous prenons trois points qui sont sur la même ligne, le point médian aura deux demi-bords infinis qui ne sont pas connectés ensemble. Cela ne nous convient pas beaucoup, car nous ne pourrons pas accéder à l'un des demi-bords, car la cellule est une liste chaînée d'arêtes.

Pour résoudre ces problèmes, nous limiterons le diagramme. J'entends par là que nous limiterons chaque cellule du diagramme afin qu'elles n'aient plus de bords infinis et que chaque cellule soit un polygone fermé.

Heureusement, l'algorithme de Fortune nous permet de trouver rapidement des bords sans fin: ils correspondent à des demi-bords toujours dans le littoral à la fin de l'algorithme.

Mon algorithme de restriction reçoit une boîte en entrée et se compose de trois étapes:

1. Il fournit le placement de chaque sommet du diagramme à l'intérieur du rectangle.
2. Coupez chaque bord infini.
3. Ferme les cellules.

L'étape 1 est triviale - nous avons juste besoin d'agrandir le rectangle s'il ne contient pas de sommet.

L'étape 2 est également assez simple - elle consiste à calculer les intersections entre les rayons et le rectangle.

L'étape 3 n'est pas non plus très compliquée, ne nécessite qu'une attention particulière. Je l'exécute en deux temps. Tout d'abord, j'ajoute les points d'angle du rectangle aux cellules aux sommets dont ils devraient être. Ensuite, je m'assure que tous les sommets de la cellule sont reliés par des demi-bords.

Je vous recommande d'étudier le code et de poser des questions si vous avez besoin de détails sur cette partie.

Voici le diagramme de sortie de l'algorithme de délimitation:


Nous voyons maintenant que toutes les arêtes sont dessinées. Et si vous effectuez un zoom arrière, vous pouvez vous assurer que toutes les cellules sont fermées:

Intersection avec rectangle

Super! Mais la première image du début de l'article est meilleure, non?

Dans de nombreuses applications, il est utile d'avoir l'intersection entre le diagramme de Voronoi et le rectangle, comme indiqué dans la première image.

La bonne chose est qu'après avoir restreint le graphique, c'est beaucoup plus facile à faire. La mauvaise nouvelle est que bien que l'algorithme ne soit pas très compliqué, nous devons être prudents.

L'idée est la suivante: nous contournons son demi-bord de chaque cellule et vérifions l'intersection entre le demi-bord et le rectangle. Cinq cas sont possibles:

1. La demi-côte est complètement à l'intérieur du rectangle: on enregistre une telle demi-côte
2. La demi-côte est complètement en dehors du rectangle: on jette une telle demi-côte
3. La demi-nervure sort du rectangle: on tronque la demi-nervure et la sauvegarde comme dernière demi-nervure qui sort .
4. La demi-nervure va à l'intérieur du rectangle: on tronque la demi-nervure pour la relier à la dernière demi-nervure qui est sortie (on la sauvegarde dans le cas 3 ou 5)
5. La demi-nervure traverse le rectangle deux fois: nous tronquons la demi-nervure et ajoutons une demi-nervure pour la relier à la dernière demi-nervure qui est sortie , puis l'enregistrons en tant que nouvelle dernière demi-nervure qui est sortie .

Oui, il y a eu de nombreux cas. J'ai créé une image pour les montrer tous:


Le polygone orange est la cellule d'origine et le rouge est la cellule tronquée. Les demi-côtes tronquées sont marquées en rouge. Des nervures vertes ont été ajoutées pour relier les demi-nervures entrant dans le rectangle avec les demi-nervures qui sortent.

En appliquant cet algorithme à un diagramme borné, nous obtenons le résultat attendu:

Conclusion

L'article s'est avéré assez long. Et je suis sûr que de nombreux aspects ne sont toujours pas clairs pour vous. J'espère néanmoins que cela vous sera utile. Examinez le code pour plus de détails.

Pour résumer et m'assurer que nous n'avons pas perdu de temps en vain, j'ai mesuré sur mon portable (pas cher) le temps de calculer le diagramme de Voronoi pour un nombre de places différent:

• $$$n = 1000$: 2 ms
• $$$n = 10 000$: 33 ms
• $$ : 450 ms
• $$ : 6600 ms

Je n'ai rien avec quoi comparer ces indicateurs, mais il semble que ce soit incroyablement rapide!

Les références

• Article original de Stephen Fortune
• Computational Geometry, 3e édition par Mark de Berg, Otfried Cheong, Marc van Kreveld et Mark Overmars
• Algorithme Fortunes: une explication intuitive sur jacquesheunis.com
• Algorithme et implémentation de Fortune sur blog.ivank.net
• Introduction to Algorithms, 3e édition par Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest et Clifford Stein