Il m'arrive souvent de voir passer de jolies animations avec des diagrammes de Voronoï, mais je n'avais jamais trouvé l'occasion d'en faire une par moi même,. c'est bien dommage car la manipulation des diagramme de Voronoï est excellent un exercice de modélisation, qui peut se faire avec des structures de données très simples et permet ensuite d'obtenir des sorties graphiques amusantes. Donc aujourd'hui je vous propose d'explorer ce sujet avec l'algorithme de Lloyd.
Un diagramme de Voronoï est un découpage d'une partie $\Omega$ du plan (disons un rectangle) en polygones $C_i$ ($i=1,\dots n$) où chaque polygone est associé à un point $m_i$ (de coordonnées $(x_i,y_i)$) et correspond aux points de $\Omega$ qui sont plus proche de $m_i$ que des autres points $m_j$:
$$ \Omega=\cup_{i=1}^n C_i,~~C_i=\left\{(x,y) \in\Omega \vert \,\forall i\neq j,\,\Vert (x,y)-(x_i,y_i)\Vert\leq \Vert (x,y)-(x_j,y_j)\Vert \right\}$$
Un diagramme de Voronoï peut donc être juste défini par une liste de points (où les listes de leurs coordonnées), mais dans ce cas il faut disposer d'une fonction simple pour identifier les points du plan correspondants à chaque polygone $C_i$ et ensuite tracer le diagramme. Tout cela peut se faire en utilisant des tableaux à une ou deux dimensions, je vais l'illustrer avec scilab mais peut se faire dans n'importe quel langage.
Calcul et tracé d'un diagramme de Voronoï
L'idée la plus simple pour le faire est de discrétiser le domaine $\Omega$ par un quadrillage régulier de points et d'associer à chaque point le numéro du polygone $C_i$ auquel il appartient. Cela peut se faire facilement avec plusieurs tableaux à deux dimensions de même taille :
- deux tableaux X,Y contenant les coordonnées des points quadrillant $\Omega$
- un tableau Vdiag (de même taille que X et Y) contenant le numéro $i$ du polygone $C_i$ correspondant à chaque point
la difficulté est de calculer ce tableau Vdiag à partir de la seule liste des points $m_i$. L'algorithme le plus simple consiste à calculer les distances de chaque pixel aux n points $m_i$ pour trouver celui qui est le plus proche. C'est loin d'être optimal mais quitte à prendre une grille pas trop fine cela ne prendra pas trop de temps à calculer :
function [Vdiag, X, Y]=ComputeVoronoiDiag(x, y, tx, ty) // tx,ty= for the meshgrid of the domain // X,Y= the meshgrid itself // x,y= n centers coordinnates (in X,Y) // Vdiag=matrix (same size as X and Y) // = k for points in (x(k),y(k)) polygon // create matrices n=length(x) [X,Y]=meshgrid(tx,ty);// the mesh mx=length(tx) my=length(ty) Vdiag=0*X; // compute the distances to each of the n centers D=matrix(zeros(1,n*mx*my),mx,my,n); for k=1:n D(:,:,k)=((X-x(k)).^2+(Y-y(k)).^2); end // select the closest center for each pixel for ix=1:mx for iy=1:my [val,k]=min(D(ix,iy,:)) Vdiag(ix,iy)=k; end end endfunction
Il est bien pratique ici de pouvoir utiliser une hypermatrice pour stocker les distances de chaque points du quadrillage aux n points $m_i$ et trouver celui qui est le plus proche. c'est la commande matrix en scilab qui permet de créer cette structure de tableau à 3 dimensions (qu'on pour voir comme l'empilement de n matrices de tailles mx*my ).
Une fois la Matrice Vdiag construite le tracé se fait en scilab avec la commande graphique comme plot3d(X,Y,Vdiag) ou surf(X,Y,Vdiag) à condition d'avoir une table des couleurs ayant plus de couleurs que de polygones à différencier pour que tous les polygones soient bien différents.
function PlotVoronoiDiag(x, y, X, Y, Vdiag, titre) // X,Y= meshgrid of the domain // x,y= n centers coordinnates (in X,Y) // Vdiag=matrix (same size as X and Y) // = k for points in (x(k),y(k)) polygon // centers position tx=X(1,:) ty=Y(:,1) n=length(x)// number of polygons for k=1:n [val,indx]=min(abs(tx-x(k))) [val,indy]=min(abs(ty-y(k))) Vdiag(indy,indx)=0;// black pixel at center end //plotting diagram drawlater() clf F=gcf(); F.color_map=[0 0 0; jetcolormap(n)];// black + n colors A=gca(); surf(X,Y,Vdiag) A.children.color_mode=-1;// no grid A.rotation_angles=[0,0];// flat display colorbar() drawnow() endfunction
L'algorithme de Lloyd
si on choisit les points $m_i$ au hasard on va obtenir des polygones de tailles diverses formant un ensemble pas forcément très équilibrées. L'algorithme de Lloyd a pour objectif en partant d'une liste de points quelconque $(m_i)_{i=1,\dots,n}$ de la faire évoluer vers une liste de points associé à un diagramme de Voronoï mieux équilibrée. La formulation de L'algorithme est très simple :
tant que un des $m_i$ est modifié
pour chaque polygone $C_i$
remplacer le point $m_i$ par le barycentre de $C_i$
$x_i={1\over \vert C_i\vert}\int_{C_i} x\,dxdy$
$y_i={1\over \vert C_i\vert}\int_{C_i} y\,dydy$
numériquement les calculs d'intégrales reviennent à une simple somme sur les points de polygone concerné:
$${1\over \vert C_i\vert}\int_{C_i} x\,dxdy\rightarrow {1\over \#\{(k,l)\vert Vdiag(k,l)=i\} } \sum_{Vdiag(k,l)=i} X(k,l) $$
Une fois l'lagorithme implémenté il suffit ensuite d'afficher le polygone à chaque étape pour obtenir une animation de l'algorithme de Lloyd.
function [x, y, Vdiag]=Lloyd_algo(x, y, X, Y, Vdiag) // X,Y= the meshgrid itself // x,y= n centers coordinnates (in X,Y) // Vdiag=matrix (same size as X and Y) // = k for points in (x(k),y(k)) polygon // compute a well spaced Voronoi Diagram n=length(x)// nomber of polygons erreur=1 q=0 while erreur>0 q=q+1 Vdiag0=Vdiag // compute new centers for k=1:n ind=find(Vdiag==k) C=length(ind) x(k)= sum(X(ind))/C y(k)=sum(Y(ind))/C end //compute new polygones tx=X(1,:),ty=Y(:,1) [Vdiag,X,Y]=ComputeVoronoiDiag(x,y,tx,ty) PlotVoronoiDiag(x,y,X,Y,Vdiag) A=gca() A.title.text="$$ \text{Lloyd''s Algorithm step}="+msprintf("%3.d",q)+"\\\;$$" A.title.font_size=6 erreur=length(find(Vdiag<>Vdiag0))// number of modifications end endfunction
Il ne reste plus qu'à appeler les fonctions précédentes avec une liste de points $m_i$ pour obtenir l'animation :
n=13// number of polygons //create lists of m_i coordinnates x=0.2+0.8*grand(1,n,'unf',0,1); y=0.2+0.8*grand(1,n,'unf',0,1); // datas for the mesh dt=0.02; // mesh step tx=[0:dt:1];ty=tx;
// start the algorythm
[Vdiag,X,Y]=ComputeVoronoiDiag(x,y,tx,ty); PlotVoronoiDiag(x,y,X,Y,Vdiag) [x,y,Vdiag]=Lloyd_algo(x,y,X,Y,Vdiag); //redraw diagramm with smooth mesh dt=0.002;tx=[0:dt:1];ty=tx; [Vdiag,X,Y]=ComputeVoronoiDiag(x,y,tx,ty); PlotVoronoiDiag(x,y,X,Y,Vdiag)
Si $\Omega$ est un carré et que le nombre de polygones est un carré $n=p^2$ l'algorithme devrait conduire à un découpage de $\Omega$ en n carrés sur p lignes et p colonnes. C'est ce qu'on observe approximativement pour n pas trop grand, le résultat final étant limité par la taille du pas de discrétisation de $\Omega$.
Très intéressant!!!!!!! Merci.
RépondreSupprimer