Pour plus de détails...

Considérons un nuage de points du plan. On souhaite délimiter des frontières. Ces frontières délimitent des zones. Sur la figure suivante, la zone de A est l'ensemble des points du plan qui sont plus proches de A que de n'importe quel autre point du nuage. La zone de B est l'ensemble des points les plus proches de B.

La frontière entre A et B est formée des points à égale distance de A et B. Ce segment suit donc la médiatrice de [AB].

Le point marqué en noir n'est pas un point du nuage : c'est le point à égale distance de A, B et C. C'est donc le centre du cercle-circonscrit à ABC. Ce point est à la rencontre des frontières AB, AC et BC.

Wikipedia

Il existe d'autres algorithmes mais celui-ci est le plus utilisé. Il permet une résolution en $n\log(n)$ mais atteindre cette limite théorique nécessite d'être très soigneux : par exemple nous aurons des listes au contenu changeant et dans lesquelles il faudra chercher. Selon la méthode de recherche, on peut perdre du temps.

Toutefois, l'algorithme n'est pas simple à bien comprendre aussi ne chercherons nous pas à trop compliqué. Retenez simplement que c'est un bon algorithme.

L'algorithme de Fortune repose sur un balayage. Nous ferons un balayage horizontal, de gauche à droite. Nous aurons donc une ligne d'équation $x = k$ avec $k$ qui augmentera au cours de l'exécution.

La droite de balayage étant placée. Appelons là $\Delta$. Considérons un point A du nuage. On s'intéresse à la ligne des points à égale distance de A et de $\Delta$.

Cette ligne est une parabole de foyer $A$ et de directrice $\Delta$. Connaissant les coordonnées $(x_A ; y_A)$ et l'équation de $\Delta$, on peut calculer l'équation de cette parabole.

Je ne vous demande pas de le faire. Je vais fournir toutes les fonctions utiles.

Lors du balayage, plusieurs points forment ce qu'on appellera la plage.

Dans l'exemple à gauche, il suffira de préciser, de bas en haut, ABCA. Cela veut dire ici que le bas de la plage est constitué par un arc de parabole de foyer A (en bleu), il est suivi d'un arc de parabole de foyer B (rouge), puis un arc de parabole de foyer C (vert) et enfin de nouveau un arc de parabole de foyer A (bleu). On peut calculer l'équation de chaque parabole et les coordonnées des intersections a, b, c.

Les positions des points a, b, c sera importante car ce sont ces points qui tracent les segments du diagramme de Voronoï.

Donc, en même temps que l'on gardera la progression de la plage ABCA, il faudra conserver la progression des points d'intersections. Ici on pourrait dire A-a-B-b-C-c-A.

Notez bien sur le graphique que la plage est constituée des arcs de paraboles les plus à droite.

Là encore, par ces explications je vous donne le principe de l'algorithme. Vous n'aurez pas à tracer de paraboles et je vous fournirai le code qui prend en charge tout l'aspect mathématique.

Supposons maintenant qu'avec la plage telle qu'elle était précédemment, la droite de balayage balaie un nouveau point D. Il faut prendre en compte ce nouveau point en ajoutant la parabole correspondante dans la plage.

Puisque D est pour l'instant sur $\Delta$, la parabole est dégénérée (c'est un terme utilisé en mathématiques pour les cas particuliers comme un triangle aplati).

Il nous faut chercher, sur la plage en cours, où tombe la parabole rose de foyer D. Autrement dit, on cherche le point sur la plage qui a l'ordonnée $y = y_D$.

Il faut donc parcourir la plage A-a-B-b-C-c-A. Puisque $y_a \leqslant y_D \leqslant y_b$, on en conclut que le point recherché doit être sur l'arc de foyer B (rouge). Maintenant que l'on sait cela, on peut calculer les coordonnées du point d. On sait déjà que $y_d = y_D$. Il s'agit de calculer, pour la parabole de foyer B, le $x$ correspondant.

Maintenant que l'on sait cela, on doit mettre à jour la plage. L'apparition de D coupe B en deux. On aura donc : A-a-B-d-D-d-B-b-C-c-A.

On sait que le point d est le début du polygone qui va faire la frontière de la zone de D dans le diagramme de Voronoï. Pour l'instant, sur la plage, D a sur sa gauche et sur sa droite le même point d. Cela peut surprendre. Mais pensez que $\Delta$ va continuer de progresser vers la droite. Quand $\Delta$ avance, les paraboles s'élargissent. La parabole en rose qui est pour l'instant toute aplatie va donc s'élargir. Les deux frontières de part et d'autre de D vont s'écarter.

Sur cette figure, la droite de balayage a légèrement progressé vers la droite. Cela a pour effet de mettre à jour toutes les paraboles constituant la plage et donc les points d'intersections.

J'ai laissé les points d'intersections antérieurs. Ainsi on voit le parcours du pour $a$ depuis sa position précédente : il a suivi le segment $[a_0a]$. Ce segment formera la frontière entre A et B.

Les points $e$ et $f$ sont les intersections de part et d'autre de D. Quand D est apparu (quand $\Delta$ est passé sur D), les points $e$ et $f$ étaient confondus en $d$. Avec la progression de $\Delta$, ils se sont séparés. On remarque toutefois que $e$, $d$ et $f$ restent alignés.

Mais ici, le point le plus important est que certains arcs de parabole constituant la plage vont en diminuant. Par exemple le premier arc rouge (en bas) rétrécit. Les points $a$ et $e$ se sont rapprochés et finiront par se rejoindre. Il en va de même pour le second arc rouge : les points $f$ et $b$ vont se rejoindre.

En revanche, l'arc vert a grandi. Les points $b$ et $c$ ne se rapprochent pas, au contraire.

Certains arcs tendent à disparaître, d'autres non. Nous avons besoin de savoir :

• quels arcs disparaissent,
• quand ils disparaissent.

Sur la figure précédente, la plage était A-a-B-e-D-f-B-b-C-c-A. Pour identifier les cercles, nous allons parcourir les triplets de foyers consécutifs et voir comment ils sont parcourus. Cela nous permettra de répondre.

• premier triplet, ABD. Le parcours ABD se fait en sens horaire (regardez la figure). Dans ce cas, on sait que dans A-a-B-e-D, les points $a$ et $e$ vont se rejoindre et que l'arc de foyer B va disparaître.
• deuxième triplet, BDB. Ces points sont alignés (il n'y en a que 2 différents). On peut passer.
• 3e, c'est DBC. Là encore, le parcours est en sens horaire. Donc dans D-f-B-b-C, on sait que $f$ et $b$ vont se rejoindre et que l'arc de foyer B va disparaître.
• 4e et dernier, BCA, le parcours est en sens anti-horaire. On passe.

Je n'ai pas cherché à prouver qu'il en serait forcément ainsi. Le critère de sens horaire / anti-horaire fonctionne néanmoins.

Nous savons donc que 2 arcs vont disparaître. Mais que faire de cela ?

Prenons l'exemple du premier triplet, ABD. Nous devons déterminer le cercle passant par ces trois points. Sur la figure, ce cercle a le centre $g$. Nous savons que $a$ et $e$ vont se rejoindre en $g$.

Autrement dit, le segment qui fait la frontière entre A et B et qui suit le point $a$, va se prolonger jusque $g$ et s'arrêter là. De même la frontière entre B et D et qui suit le point $d$, va se prolonger jusque $g$.

Quand cela arrivera, on passera de A-a-B-e-D-… à A-g-D-… L'arc B aura disparu.

On sait aussi que cela se produira quand $\Delta$ aura atteint le bord droit du cercle circonscrit.

La structure retenue pour le diagramme de Voronoï peut être plus ou moins compliquée, selon les besoins. On peut par exemple chercher pour chaque point du nuage à former des polygones fermés.

Nous nous contenterons ici de plus simple : nous enregistrerons simplement des segments.

Sur la plage, entre deux paraboles, il y a toujours un segment qui fait la frontière.

Il n'est pas utile d'enregistrer tous les points du segment, seulement son point de départ et son point d'arrivée. Quand un nouveau point est pris en compte et ajouté dans la plage, on démarre deux nouveau segment de part et d'autre de ce point. Quand un cercle est pris en compte, on ferme deux segment et on en ouvre un nouveau.

Je vais donner plus loin un exemple détaillé, vous comprendrez mieux.

Il ne sera pas utile de faite balayer $\Delta$ pixel à pixel. Au lieu de de cela nous allons créer une file avec priorité, une queue. L'ordre de priorité sera d'abord la coordonnées $x$, puis la coordonnées $y$. Dans cette queue commencerons par placer les points du nuage. Dans l'exemple : A, B, C, D. Puis d'éventuels cercles s'ajouteront.

Considérons les points A, B, C, D comme précédemment.

Nous aurons simultanément une queue. Celle si commence avec [A, B, C, D] dans cet ordre (x croissant). Et la plage, qui est vide pour l'instant.

La ligne de balayage se place sur A qui doit être ajouté à la plage. La plage est vide. Il suffit donc de l'initialiser avec A.

la ligne de balayage se place sur B. On cherche, dans la plage, quel arc est intersecté par B. C'est forcément A. On cherche l'intersection. C'est le point a.

La plage est donc ABA, on démarre deux segments de chaque côté de B. On pourrait le noter ainsi : A - [a,?] - B [a, ?] - A.

Autour de B que l'on vient d'insérer, on regarde s'il n'y aurait pas un cercle valable. ABA n'est pas valable.

On dépile C. On cherche dans la plage où intersecte C. Pour cela, comme la plage est ABA, on cherche l'intersection pour A-B, ce qui nous donne b. L'ordonnée de b est inférieure à celle de C on peut continuer. L'intersection pour B-A est le point c. Son ordonnée est encore inférieure à celle de C. On conclut que C intersecte le 2e arc A de la plage. L'intersection est e. On obtient donc la nouvelle plage :

A - [a,?] - B - [a,?] - A - [e,?] - C - [e, ?] - A

Autour de C on cherche un cercle qui serait valable : BAC, ACA (inutile de regarder ceux qui ne contiennent pas C). BAC est valable et est ajouté dans la queue qui contient maintenant [BAC, D] car BAC devra être pris en compte pour l'abscisse correspondant à son bord droit, donc avant D.

On dépile un cercle. On sait que ce cercle, BAC, a le centre f et est lié à l'arc A. Voici ce que l'on doit faire :

• plage : A - [a,?] - B - [a,?] - A - [e,?] - C - [e, ?] - A,
• recherche de la séquence BAC,
• fermer autour de A : A - [a,?] - B - [a,f] - A - [e,f] - C - [e, ?] - A
• supprimer A et les segments (qu'il faut garder quelque part),
• boucher le trou avec un nouveau segment commençant en f : A - [a,?] - B - [f,?] - C - [e, ?] - A

Notez que dans une situation compliquée, il pourrait arriver que l'on doive traiter un événement cercle comme BAC mais que la séquence BAC ait disparu (la plage n'arrête pas de changer). Dans ce cas pas de problème, si la séquence est introuvable, le cercle peut être ignoré.

Il y aurait moyen de ne pas avoir à chercher BAC et pointer directement A dans la plage mais je préfères ralentir un peu notre algo pour le rendre plus simple.

Notez également que chaque fois que l'on retire ou que l'on ajoute quelque chose à la plage, il faut vérifier, autour de la modification, la présence d'éventuels cercles à traiter. En faisant cela, il sera peut-être compliqué d'éviter de considérer plusieurs fois le même cercle. On fera donc en sorte que, si on entre deux fois un même cercle dans la queue, la 2e sera ignorée.

Ici, ABC et BCA ne sont pas valides.

• La plage est A - B - C - A.
• A - B donne le point b, en dessous de D,
• B - C donne h, au-dessus de D. Donc D intersecte B.
• l'intersection B - D est en i.

La plage devient :

A - [a,?] - B - [i,?] - D - [i,?] - B - [f,?] - C - [e, ?] - A

On pense à supprimer d'éventuels événements cercle auraient été lié à l'arc B qui vient d'être coupé en deux (il n'y en avait pas). Puis on cherche les cercles parmi ABD, BDB et DBC. ABD et DBC sont valables et sont placés dans la queue. ABD est prioritaire (bord droit plus à gauche). La queue est donc [ABD, DBC]

On trouve ABD dans la plage. C'est le point au centre qui est susceptible d'être enlevé : B. Le centre du cercle est j.

A - [a,?] - B - [i,?] - D - [i,?] - B - [f,?] - C - [e, ?] - A

devient A - [a,j] - B - [i,j] - D - [i,?] - B - [f,?] - C - [e, ?] - A.

Les segments [a,j] et [i,j] sont gardés quelque part. On peut remplacer dans la plage :

A - [j,?] - D - [i,?] - B - [f,?] - C - [e, ?] - A.

Pas de cercles à effacer ni de cercle à ajouter.

Idem, on trouve DBC. C'est B qui devra être enlevé et le centre du cercle est k.

A - [j,?] - D - [i,?] - B - [f,?] - C - [e, ?] - A

devient A - [j,?] - D - [i,k] - B - [f,k] - C - [e, ?] - A.

Les segments [i,k] et [f,k] sont gardés quelque part. On peut remplacer dans la plage :

A - [j,?] - D - [k,?] - C - [e, ?] - A.

Pas de cercles à effacer ni de cercle à ajouter.

La queue est terminée. Tous les segments intérieurs sont terminés. Nous avons ainsi [a,f], [e,f], [a,j], [i,j], [i,k], [f,k]. Mais les segments [j,?], [k,?], [e,?] restent ouverts.

Une possibilité est de :

• demander à $\Delta$ de se décaler encore un peu
• clôturer tous les segments avec les intersections correspondantes. Par exemple pour A - [j,?] - D, on cherche l'intersection A - D, par exemple p (en dehors du graphique) et on clôture [j,p].
• De la même façon, D - C donne n, donc on clôture [k,n] ; C - A donne q (en dehors de la figure), donc on clôture [e,q].

Ces segments supplémentaires vont théoriquement jusqu'à l'infini. Si on veut aller jusqu'au bout, il faudrait déterminer une zone de cadrage et interrompre les segments convenablement sur le bord du cadre. Certains des segments calculés peuvent d'ailleurs être entièrement hors cadre. Je ne cherche pas à couvrir ce problème.

Suite à cela, nos segments sont [a,f], [e,f], [a,j], [i,j], [i,k], [f,k], [j,p], [k,n] et [e,q]. Certains comme [i,j] et [i,k] pourraient être réduits en un segment [j,k] (on est certain que les 3 points i, j, k sont alignés) mais je laisse aussi de côté ce problème.

Comme on connaît tous les segments, on peut faire un tracé.

Je vous propose d'utiliser les classes Points, Circle, PriorityQueue et Segment ainsi que du module parabole disponibles dans ce fichier voronoi.zip.

Ci-dessous quelques explications.

Ces classes contiennent beaucoup de méthodes privées (celles dont le nom est en __). Vous pouvez bien sûr chercher à comprendre mais ce n'est pas indispensable. L'idée de la différence privée / public est que je me contente de vous garantir que les fonctions privées font ce qu'elles doivent, sans avoir à détailler ce qui reste secret dans le privé.

Définition d'un point avec ses deux coordonnées x et y. Le point est immuable : une fois le point créé, les coordonnées ne peuvent être changées.

On peut tester l'égalité de deux points. Deux points sont considérés égaux s'ils sont même coordonnées.

La classe contient d'autres fonctions pour faire des calculs. Je m'en suis servi dans les autres classes. Il s'agit surtout de calculs vectoriels. Mais vous ne devriez pas en avoir besoin.

>>> from point import Point
>>> A = Point(4, 2)
>>> B = Point(2, 6)
>>> C = Point(4, 2)
>>> A == B
False
>>> A == C
True

Définition d'un segment selon ses extrémités qui sont des points. À la création, on doit indiquer le point de départ du segment. On peut aussi préciser le point de fin mais pas forcément. Plus tard, on peut indiquer le point final du segment. Notez que les extrémités ne peuvent être changées une fois écrites.

J'ai ajouté une méthode boxed qui permet de découper un segment pour qu'il rentre dans une boîte de taille fixé en précisant les valeurs : gauche, droite, haut, bas.

>>> from point import Point
>>> from segment import Segment
>>> A = Point(4,2)
>>> B = Point(6,12)
>>> s = Segment(A,B) # fourni complet tout de suite
>>> s2 = Segment(A)  # seulement le début est fourni
>>> s2.finish(B)     # achevé plus tard

Une classe de cercle. On lui fournit trois points qui doivent être valables dans le cadre du problème : pas alignés et en sens horaire. La classe se charge de calculer centre et rayon du point.

Il ne faut pas chercher à créer un cercle avec des points non valables, donc il faut tester avant.

>>> from point import Point
>>> from circle import Circle
>>> A = Point(0,1)
>>> B = Point(10,1)
>>> C = Point(9,4)
>>> Circle.points_are_valid(A,B,C) # mauvais sens
False
>>> Circle.points_are_valid(B,A,C)
True
>>> c = Circle(B, A, C)
>>> c.center
Point[5.0, 1.0]
>>> c.rayon
5.0

Gestion de la queue pour les événement. Une fois créée, vous pouvez y insérer des objets avec push. uniquement des Point et des Circle !. Pour extraire le prochain événement, utilisez pop. Vous pouvez vérifier que la queue est vide avec empty.

>>> from point import Point
>>> from priorityqueue import PriorityQueue
>>> A = Point(10,2)
>>> B = Point(5,1)
>>> q = PriorityQueue()
>>> q.push(A)
>>> q.push(B)
>>> q.pop()
Point[5,1]
>>> q.empty()
False
>>> q.pop()
Point[10,2]
>>> q.empty()
True

Ce n'est pas une classe mais un module contenant toutes les fonctions utiles liées au parabole selon l'algorithme de Fortune. Il n'y en a que 2 : get_x et intersection. On définit les paraboles par leur foyer (de type Point) et par la position de leur directrice, c'est à dire l'abscisse de la ligne de balayage, xd, un float.

Fonction get_x

Elle reçoit le foyer f, une ordonnée y et l'abscisse de la directrice xd. Considérant la parabole de foyer f et de directrice en xd, on veut connaître la position du point d'ordonnée y. Puisque l'on connaît déjà l'ordonnée du point, on ne demande que son abscisse. La fonction renvoie l'abscisse.

Fonction intersection

On fournit deux foyers f0 et f1 et la directrice. L'ordre des foyers est important comme on l'a expliqué plus haut. Il faudra donner f0 et f1 dans l'ordre, de bas en haut, où ils viennent sur la plage. La fonction renvoie le Point d'intersection.

Pour l'utiliser vous pouvez faire quelque chose comme :

>>> import parabole
>>> from point import Point
>>> f = Point(4,8)
>>> parabole.get_x(f, 9, 6) # y = 9, directrice en 6
4.75

Autrement dit, le point de coordonnées (4.75,9) est à égale distance du point f et de la ligne verticale en xd (je vous rappelle qu'on parle de paraboles parce que ce sont justement des lignes à égale distance d'un point et d'une droite)

Je vous recommande de faire 2 choses : une classe ou un module que l'on pourra appeler voronoi et une classe que l'on appellera beach.

• voronoi prendra en charge la structure globale de l'algorithme.
• beach devra gérer la plage.

Ce module doit dérouler le fonctionnement suivant :

• À la création, on fournit un nuage de points, c'est à dire une liste de paires (x,y), ainsi, éventuellement, qu'un cadre.
• Création d'une queue q.
• Dès la création, les paires (x,y) sont convertis en Point et ajoutés à la queue. Pas besoin de se soucier d'ordre, c'est la queue qui s'en charge.
• Création d'une plage b initialement vide.

Ensuite, lancement du processus suivant :

TANT QUE q n'est pas vide:
    soit e le prochain événement dans q
    SI e est un Point:
        insérer e dans la plage b
        détecter tous les cercles valables autour de cette insertion
        insérer les cercles dans la queue
    SINON:
        c'est un cercle...
        charger la plage de prendre en charge ce cercle
        détecter tous les cercles valables autour de cette modification
        insérer les cercles dans la queue
    FIN
FIN
Avancer la ligne de balayage encore un peu pour clôturer les segments encore ouverts

Ce module doit conserver la liste des points constituant la plage et aussi la liste des segments faisant frontière entre les points.

Quand certains segments sont terminés (cela arrive quand on gère un cercle), ils sont retirés de la plage elle-même. Il faut les conserver quelque part. Par exemple on prévoit une list de segments terminés que l'on pourra récupéré quant tout sera terminé.

création de la plage

On se contente d'initialiser 3 conteneurs, par exemple des list. Un pour les points, un autre pour les segments en cours de construction, un autre pour les segments terminés.

insertion d'un point

Si la plage est vide, s'il suffit d'ajouter le point, c'est tout.

Sinon, il faut chercher le point correspondant à la parabole face à laquelle tombe le point inséré.

Une fois trouvé, il faut insérer le point comme l'a expliqué plus haut et créer deux segments correspondant.

détection de cercles

Idéalement il faudrait ne considérer que les points autour des modifications faites dans la plage. Mais comme c'est déjà assez difficile, on peut se contenter d'une solution moins optimisée et prévoir une fonction qui parcours tous les triplets à la recherche de cercles.

Cela fait perdre un peu de temps à l'exécution et pourrait occasionner qu'un même cercle serait ajouté plusieurs fois. Mais la queue est prévue pour ignorer les ajouts répétés, donc pas de problème de ce côté.

prise en compte d'un cercle

Si un cercle est par exemple ABC, il faut chercher la séquence ABC dans les points. Si on ne la trouve pas, ne rien faire. Si on la trouve, il faut : supprimer B de la plage ; fermer les deux segments autour de B par le centre du cercle ; placer entre A et C un nouveau segment commençant au centre du cercle.

dernier pas

pour clôturer à la fin, avancer un peu la ligne de balayage puis par exemple si la plage est ABCD :

• calculer l'intersection A, B et terminer le segment entre A et B par ce point,
• idem pour l'intersection entre B et C et le segment entre B et C
• ainsi de suite

Il sera difficile de juger du résultat si on ne peut pas le voir graphiquement.

Prévoyez donc une démonstration où il s'agira de créer un nuage de points, de représenter ce nuage, de faire fonctionner l'algorithme de fortune, de tracer enfin tous les segments.

Pensez que le principe du diagramme de Voronoï est basé sur des distances. Si on veut que les distances apparaissent bien sur le graphique, celui-ci doit utiliser des axes orthonormés.

import matplotlib.pyplot as plt

...

# en supposant que nuage est la liste de Point du diagramme
xs_nuage = [pt.x for pt in nuage]
ys_nuage = [pt.y for pt in nuage]
plt.scatter(xs_nuage, ys_nuage)

# en supposant que segs est la liste des segments
for s in segs:
    p1, p2 = s.points
    plt.plot([p1.x, p2.x], [p1.y, p2.y])

# zoom sur la zone utile, en supposant que
# left, right, top, bottom définissent le cadre
plt.xlim(left, right)
plt.ylim(bottom, top)
    
# orthonormé :
ax = plt.gca()
ax.set_aspect(1)

# affichage
plt.show()

Si vous avez réussi, vous aurez un beau diagramme de Voronoï. Peut-être certains segments ne vont-ils pas au bord ou bien dépassent du bord. Vous pouvez aller plus loin en essayant d'exploiter la fonction boxed de la classe Segment pour extraire juste ce qu'il faut pour que les segments touchent les bords et pas plus loin.