Algorithmes géométriques

Les algorithmes géométriques sont utiles pour les applications en graphisme par ordinateur ainsi que pour les applications de CAO qui travaillent sur des structures géométriques en 2 ou 3D. Les principaux algorithmes utilisés sont discutés ici, à l'exception des algorithmes classiques de dessin, découpage et surfaces cachées (qui sont couverts spécifiquement par les cours d'infographie).

Les structures géométriques de base

Les structures suivantes sont à la base des algorithmes présentés par la suite:

-

Le point (coordonnée x, y, z dans l'espace cartésien). Ces valeurs peuvent être entières si l'espace adressé s'y prête (écran avec pixels) ou réelles.

-

La ligne définie par les deux points aux extrémités. Elle peut aussi s'exprimer sous forme paramétrique: $x = a_xs + b_x$ et $y = a_ys + b_y$, s allant de 0 à 1.

-

Un polygone est une surface définie par une liste de points qui englobent en sens horaire (ou anti-horaire selon la convention utilisée) sa surface interne.

-

Une surface peut aussi être définie par une suite de courbes paramétriques du troisième degré (souvent suivant la formulation de Bézier). On peut bien sûr avoir aussi des courbes du second degré (les lignes elles sont du 1er degré) ou des courbes de degré plus élevé mais le troisième degré représente un excellent compromis puisqu'elles se manipulent assez bien et permettent de garantir la continuité de position et de pente.

-

On peut aussi avoir des solides formés avec des surfaces paramétriques.

Retour sur les courbes de Bézier

-

Une courbe paramétrisée du troisième degré est représentée par 4 points: P1, P2, P3, P4.

-

En effet cela permet de fixer les quatres coefficients de l'équation paramétrique: $x = a_x s^3 + b_x s^2 + c_x s + d_x$, même chose en y, avec s entre 0 et 1.

-

P1 est le début de la courbe $x(0) = P1 = d_x$.

-

P4 est la fin de la courbe $x(1) = P4 = a_x + b_x + c_x + d_x$.

-

P2 indique la direction de la courbe en sortant de P1: $x'(0) = 3 (P2 - P1) = c_x$.

-

P3 indique la même chose en allant vers P4: $x'(1) = 3 (P4 - P3) = 3 a_x + 2 b_x + c_x$.

-

Nous avons donc les équations matricielles suivantes:

$$\left[ \begin{array}{r} P_1 P_2 P_3 P_4 \end{array... ... \begin{array}{r} a_x b_x c_x d_x \end{array}\right]$$

$$\left[ \begin{array}{r} a_x b_x c_x d_x \end{array... ... \begin{array}{r} P_1 P_2 P_3 P_4 \end{array}\right]$$

Intersection de deux lignes

Afin de déterminer si deux lignes s'intersectent, on peut calculer leur équation paramétrique et vérifier pour quelle valeur de t les deux droites on leur x et y identiques. Si cela arrive pour t entre 0 et 1 dans les deux droites, elles s'intersectent.

On peut aussi au préalable vérifier que les intervalles en x et y des deux segments de droite s'intersectent. Si on s'attend à ce qu'il y ait peu d'intersections, ce simple test peut rapidement éliminer plusieurs cas où les droites sont éloignées, ne demandant le calcul plus complexe seulement lorsque les droites ont une forte probabilité de se recouper.

Il existe aussi un autre algorithme intéressant:

int ccw(point p0, point p1, point p2)

{
  int dx1, dx2, dy1, dy2;

  dx1 = p1.x - p0.x;
  dx2 = p2.x - p1.x;
  dy1 = p1.y - p0.y;
  dy2 = p2.y - p1.y;

  if(dx1 * dy2 > dy1 * dx2) return(1);
  if(dx1 * dy2 < dy1 * dx2) return(-1);

  if(dx1 * dx2 < 0 || dy1 * dy2 < 0) return(-1);
  if(dx1 * dx1 + dy1 * dy1 >= dx2 * dx2 + dy2 * dy2) return(0);
  return(1);
}

int intersect(ligne l1, ligne l2)

{
  return(ccw(l1.p1, l1.p2, l2.p1) * ccw(l1.p1,l1.p2,l2.p2) <= 0 &&
         ccw(l2.p1, l2.p2, l1.p1) * ccw(l2.p1,l2.p2,l1.p2) <= 0)
}

La première fonction compare la pente des deux lignes formées par les trois points et détermine donc si elles forment un tournant en sens horaire ou non. La seconde fonction vérifie si on tourne dans des direction différentes d'une ligne vers les deux points de l'autre (ce qui indique que la droite associée à la première ligne coupe la seconde ligne) et vice-versa. Si c'est le cas pour les deux, on a automatiquement intersection.

Cet algorithme géométrique est intéressant et pas nécessairement intuitif. Il est robuste et simple à programmer mais, sous sa forme présente n'est pas nécessairement le plus efficace.

-

Si d'une ligne vers chaque extrémité de l'autre on tourne d'un côté différent, le prolongement de la première ligne coupe la seconde ligne.

-

Si de la seconde ligne on tourne encore de côtés différents en allant vers les extrémités de la première, les deux lignes intersectent.

Intersection de courbes de Bézier

La plupart des applications qui calculent l'intersection entre deux courbes de Bézier transforment ces courbes en une suite de lignes ou retardent cette étape pour travailler au niveau pixel. Ceci s'explique par le fait que l'intersection entre deux courbes peut produire 9 points et son calcul est équivalent à trouver la racine d'un polynome du troisième degré.

Heureusement il est aussi possible de découper une telle courbe en quelques segments ayant des propriétés particulières de monotonocité se sorte qu'on puisse utiliser un algorithme itératif entre chaque paire de segments avec une convergence assurée. Une telle technique est très peu utilisée cependant.

-

On peut remplacer la courbe par une approximation formée de plusieurs petits segments droits et vérifier ces segments deux à deux d'une courbe à l'autre.

-

On peut résoudre le système formé des équations des deux courbes, ce qui équivaut à chercher les racines d'un polynome du neuvième degré.

-

On peut découper les courbes en sections concaves et convexes et utiliser un algorithme itératif à convergence assurée.

Le problème du chemin simple fermé

Il arrive qu'on ait un ensemble de points par lesquels on veut passer sans jamais croiser son chemin. Pour trouver un tel chemin (et non pas pour trouver le meilleur chemin de ce genre) il existe un algorithme simple et relativement efficace. Il s'agit de prendre un point de départ puis de tracer une ligne de ce point vers chaque autre point pour calculer l'angle associé et ensuite visiter les points par ordre croissant de cet angle.

Dans la mesure où seulement le tri nous intéresse, il est possible d'utiliser une fonction autre que arctangente(dy/dx) à la condition que cette fonction ait des propriétés similaires.

-

choisir un point de départ.

-

tracer une ligne de ce point vers chaque autre.

-

calculer l'angle de chaque ligne avec l'horizontale.

-

visiter les points par ordre croissant d'angle.

-

complexité du tri.

float theta(point p1, point p2)

{
  float dx, dy, ax, ay;
 
  dx = p2.x - p1.x;
  ax = abs(dx);
  dy = p2.y - p1.y;
  ay = abs(dy);
  if(dx == 0 && dy == 0) t = 0;
  else t = dy / (ax + ay);
  if(dx < 0) t = 2 - t;
  else if(dy < 0) t = 4 + t;
  theta = t * 90;
}

Ceci évite de calculer la relation non linéaire entre l'angle et sa tangeante. Il suffit de placer le tout dans le bon quadrant et de multiplier par 90 si on veut, pour des raisons cosmétiques.

Inclusion dans un polygone

Pour savoir si un point est à l'intérieur d'un polygone, la méthode la plus utilisée est de tracer une droite de ce point vers l'infini et de compter les intersections entre cette droite et les segments du polygone qui la croisent. Il faut bien sûr faire attention aux coins qui tombent juste sur la demi-droite.

-

tracer une droite du point vers l'infini.

-

calculer le nombre d'intersections entre cette droite et les segments formant le polygone.

-

un nombre impair indique qu'on est à l'intérieur.

int inside(point t, polygon p, int N)

{
  count = 0;
  j = 0;
  p[0] = p[N];
  p[N + 1] = p[1];
  lt.p1 = t;
  lt.p2.x = maxint;
  lt.p2.y = t.y;

  for(i = 1; i <= N ; i++)
   { lp.p1 = p[i]; lp.p2 = p[i];
     if(intersect(lt,lp)) continue;
     lp.p2 = p[j];
     j = i;
     if(intersect(lp,lt)) count++;
   }
  return(count % 1);
}

Trouver le polygone convexe circonscrit

Lorsqu'on a un nuage de points, un sous-ensemble de ceux-ci forment le polygone convexe qui les englobe tous. Une droite, qui vient de l'infini et se rapproche des points, tombe automatiquement sur un des points de ce polygone convexe. On peut dont prendre les points avec x min, x max, y min ou y max pour obtenir immédiatement jusqu'à 4 de ces points. Pour le reste, les algorithmes ressemblent passablement au tri (en terme de complexité également qui est n log n).

La première méthode est celle de l' emballage. On prend un point garanti d'être sur ce polygone et on calcule l'angle entre ce point et chaque autre point pour enfin choisir le minimum (utilisant la fonction theta présentée plus tot). Le point choisi devient le nouveau point de départ. Le polygone est terminé lorsque le point choisi est notre point de départ initial. La complexité de cet algorithme est N * M ou N est le nomrbe de points dans le nuage et M (M $\leq$ N) est le nombre de points sur le polygone.

-

choisir un point qui est sur le polygone circonscrit (xmin, xmax, ymin ou ymax).

-

calculer l'angle entre l'horizontale et la droite reliant ce point à chaque autre.

-

retenir l'autre point pour lequel l'angle est le plus petit.

-

recommencer avec ce nouveau point.

-

arrêter lorsqu'on revient au point de départ.

-

complexité $O(M \times N)$ M nombre de points sur le polygone et N le nombre total de points.

Un second algorithme est appelé le balayage de graham. On commence par faire un polygone fermé avec tous les points, tel que vu précédemment. D'un point à l'autre on vérifie que l'on tourne toujours vers la gauche. Sinon, on obtient une portion concave à éliminer en enlevant le point correspondant. Il arrive parfois que plusieurs points soient enlevés lorsqu'un nouveau point est ajouté et que l'on vérifie la convexité.

La complexité de l'algorithme est donc n log n pour le tri initial puis n pour le reste puisque chaque point peut être visité et ensuite enlevé une seule fois. Cette procédure peut s'exprimer comme suit:

-

trouver le chemin fermé passant par tous les points en commençant par un point qui est sur le polygone circonscrit.

-

complexité $O(N \log (N))$

-

repasser par chaque point sur le chemin fermé et enlever les points qui forment des creux concaves (tournant anti-horaire).

-

cette deuxième étape est linéaire.

min = 1;

for(i = 2 ; i <= N ; i++)
 { if(p[i].y =< p[min].y) 
   if(p[i].y != p[min].y || p[i].x > p[min].x) min = i;
 }

t = p[1];
p[1] = p[min];
p[min] = t;

sort(p,N);

p[0] = p[N];
M = 3;

for(i = 4 ; i <= N ; i++)
 { while( ccw(p[M],p[M - 1],p[i]) >= 0 ) M--;
   M++;
   t = p[M];
   p[M] = p[i];
   p[i] = t;
 }

On peut remarquer que le premier point rencontré est nécessairement sur le polygone circonscrit. Les points candidats sont accumulés au début du vecteur et sont enlevés par la suite et mis entre M et i lorsqu'on découvre qu'ils sont dans une baie concave.

Une dernière technique mérite d'être examinée puisqu'elle ne peut garantir une meilleure complexité mais permet en général de réduire considérablement le nombre de points à examiner. L'idée est fort simple, il s'agit de choisir 4 points et d'enlever tous les points qui tombent à l'intérieur de ce quadrilatère.

-

choisir quatre points.

-

enlever les points qui tombent à l'intérieur.

-

appliquer une des deux méthodes précédentes sur les points restants.

En résumé, la méthode de l'emballage donne en moyenne d'aussi bons résultats que le balayage lorsque les points sont distribués uniformément et que M est d'environ log N. Par ailleurs, la méthode du quadrilatère suivie du balayage donne en moyenne une complexité linéaire est est la plus avantageuse si les points sont bien répartis sur le plan.

-

Lorsque les points sont distribués uniformément, M est comparable à $\log(N)$.

-

Dans ce cas, l'emballage et le balayage de Graham donnent des résultats similaires.

-

En prenant quatre points extrêmes (xmin, xmax, ymin, ymax), une majorité de points sont enlevés et on obtient une complexité souvent linéaire

La recherche dans l'espace 2D ou 3D

Dans l'espace unidimensionnel, la recherche des objets qui sont dans un intervalle donné revient à une recherche dans un arbre binaire avec un intervalle de clés plutôt qu'une clé unique. La complexité est donc n log n pour bâtir l'arbre et R log n pour une recherche d'un intervalle qui retourne R éléments.

Regardons maintenant le problème en 2 dimensions, avec les axes x et y. On veut trouver les objets qui tombent dans un rectangle donné. Il est possible de faire le tri sur x et une recherche linéaire pour vérifier ceux qui satisfont la condition sur y. Il est cependant possible de faire mieux.

On peut dessiner une grille et avoir une matrice qui donne pour chaque carreau la liste des objets contenus à l'intérieur. Le problème est de choisir la granularité de la grille, un peu comme pour les hash table.

Plus flexible est la méthode de diviser l'espace récursivement de l'arbre binaire/ quaternaire qui divise récursivement l'espace en deux (quatre) selon x et y en alternance.

{
  divide = X;
  val = p.x;
  for(;;)
   { old_val = val;
     if(divide == X) 
      { val = p.x;
        divide = Y;
      }
     else 
      { val = p.y;
        divide = X;
      }
     
     if(old_val < node.boundary) 
      { if(node->left == NULL) 
         { node->left = new(node);
           node->left.boundary = val;
           return;
         }
        node = node->left;
      }
     else
      { if(node->right == NULL)
         { node->right = new(node);
           node->right.boundary = val;
           return;
         }
        node = node->right;
      }
   }
}

-

1D: arbre de tri par grandeur croissante de position. Pour trouver les points qui intersectent un intervalle, il suffit de faire une recherche d'intervalle.

-

Complexité de $n \log(n)$ pour construire l'arbre et de $R \log(n)$ pour une recherche d'intervalle qui retourne $R$ objets.

-

2D: arbre de tri selon x et recherche linéaire selon y.

-

arbre quaternaire ou alternant selon x et y. Complexité $\log(n)$ si l'arbre est balancé.

-

grille ou table 2D avec fonction de dispersement; complexité linéaire si la distribution est uniforme.

Le problème devient plus complexe lorsque l'on veut trouver tous les rectangles qui sont dans (intersectent) une certaine région. En effet, lorsqu'on applique les subdivisions selon x et y, il arrive souvent qu'on tombe sur un rectangle qui chevauche la frontière. Dans un tel cas, on peut:

-

conserver le rectangle dans une liste associée au noeud frontière correspondant.

-

diviser le rectangle le long de la frontière en espérant que le nombre résultant de rectangles une fois que tout est entré dans l'arbre ne soit pas trop grand comparé au nombre original de rectangles.

-

envoyer le rectangle (un pointeur) de chaque coté de la frontière, ce qui revient pratiquement au même.

-

si les rectangles ne s'intersectent pas, une structure différente est possible. En effet, on peut mettre des rectangles de vide entre les autres et ainsi obtenir une région continue de rectangles qui se pointent les uns aux autres. La recherche a alors une complexité de $\sqrt n$ pour une région de surface $n$.

Une telle application se retrouve dans l'affichage de formes géométriques sur un écran. En effet, il faut trouver tous les rectangles dans le plan dont une portion intersecte la fenêtre virtuelle associée à l'écran.

Une autre méthode efficace pour cette application est celle du R-tree qui permet à plusieurs régions de se recouper plutôt que d'avoir les rectangles qui recoupent plus d'une région. Le principe d'opération est le suivant pour bâtir l'arbre ou lorsque des rectangles sont ajoutés:

-

Une région rectangulaire est associée à chaque noeud dans l'arbre.

-

Les feuilles contiennent une liste de rectangles, K au plus.

-

Lorsqu'on commence avec une liste ou lorsque la liste contient trop d'éléments, on choisit les deux rectangles les plus éloignés pour former deux nouvelles sous-régions. Ensuite, chaque autre rectangle est ajouté à la sous région qui s'expand le moins par son addition.

-

Pour une recherche, on fouille toutes les branches dont la région rectangulaire intersecte l'intervalle désiré. Rendu dans les feuilles, on continue avec une recherche linéaire pour les environ 10 rectangles contenus.

L'intersection de polygones

Dans plusieurs applications comme la vérification des règles de dessin des masques de fabrication en micro-électronique, il faut vérifier si des intersections (court-circuits) sont présentes parmi un très grand nombre de polygones. On peut simplement vérifier tous les objets par paires mais avec des millions de rectangles c'est plutôt impratiquable ($n^2$).

La première méthode vue est celle du livre et concerne des segments de droites horizontaux et verticaux. La procédure est comme suit:

-

Tous les segments sont triés par ordre croissant de y minimum.

-

On avance la ligne horizontale de balayage d'une valeur de y à la prochaine rencontrée dans la liste triée.

-

Le segment avec cette nouvelle valeur de y est ajouté dans un arbre qui trie les segments en balayage selon x.

-

Lorsqu'un nouveau segment vertical est ajouté, il faut vérifier s'il en existe un autre en balayage avec la même valeur de x.

-

Lorsqu'un nouveau segment horizontal est rencontré, on vérifie l'intervalle correspondant dans l'arbre pour un segment horizontal ou vertical qui l'intersecterait.

-

Les segments horizontaux ne restent qu'un instant dans l'arbre alors que les segments verticaux sont enlevés lorsque leur valeur y max est rencontrée.

Lorsque l'on veut trouver l'intersection entre des polygones, il faut utiliser un arbre plus complexe puisque les rectangles peuvent intersecter. On peut aussi utiliser une simple liste qui ne fournit pas la performance de n log n + Intersections de l'algorithme précédent.

-

chaque polygone est défini par une suite de segments qui entourent une région dans le sens horaire.

-

la ligne de balayage se déplace de $y_{min}$ à $y_{max}$.

-

les segments sont triés selon leur $y_{min}$.

-

le balayage se déplace jusqu'au prochain élément dans cette liste triée ou jusqu'au prochain élément de la liste de balayage qui doit être retiré.

-

un segment de la liste triée dont le $y_{min}$ est rejoint est ajouté à la liste de balayage (segments coupés par la ligne de balayage).

-

ou un segment de la liste de balayage dont le $y_{max}$ est dépassé est retiré.

-

la liste de balayage doit être maintenue triée en $x$.

-

on part de $x$ très petit avec un indice de recouvrement égal à 0.

-

on augmente de 1 à chaque segment rencontré qui va vers le haut et on diminue de 1 à chaque segment qui va vers le bas.

-

un indice de recouvrement de 2 ou plus indique intersection; la section correspondante du segment de polygone est marquée comme couverte et un lien est établi avec le segment qui l'intersecte.

-

on peut maintenir un indice séparé pour deux types de segments: objet à dessiner et fenêtre d'affichage. Une portion d'objet est alors affichée seulement lorsque les deux indices de recouvrement sont à 1.

-

complexité entre $n \log(n)$ et $n^2$ selon le genre de polygones.

-

en moyenne on a $\sqrt n$ étapes de balayage avec $\sqrt n$ polygones à la fois dans la liste de balayage.

-

on peut aussi étendre des polygones de $d/2$ et vérifier leur intersection pour savoir si une distance minimale de $d$ est maintenue entre eux.

Recherche des proches voisins

Il arrive souvent que l'on veuille trouver les deux plus proches voisins dans un nuage. Pour ce faire, on divise récursivement en 2 l'espace selon x et on maintient les y triés à mesure des fusions en x. Ceci donne la procédure:

-

trier selon x

-

trier selon y en divisant récursivement en 2 l'espace déja trié en x. A chaque division on sépare l'espace selon x et on en retire deux listes triées en y qui seront jointes par fusion.

-

Dans chaque moitié, on sait que les deux points les plus proches sont à une distance min. On retient le meilleur des deux.

-

Il faut ensuite vérifier les couples formés par des points de part et d'autre de la séparation en x, pas plus éloignés que de min.

-

On traverse la liste fusionnée en y pour examiner les points consécutifs qui sont dans l'intervalle min de la séparation en x. Ceci permet de trouver tous les couples de points de distance inférieure à min qui sont de part et d'autre de la séparation.

-

Comme cette procédure suit le merge sort, elle est de complexité n log n.

On peut généraliser la procédure précédente pour trouver toutes les paires de voisins les plus proches. Ceci est un problème similaire à celui de trouver le diagramme de Voronoi.

Il suffit de décomposer récursivement l'espace en 2 selon x puis de relier les deux points laissés ensemble dans une partition. Ensuite, de partition en partition, les points en périphérie des polygones circonscrits formés, qui sont du côté de la frontière où se fait la fusion sont joints en suivant les y triés.

Conclusion

En conclusion, les algorithmes géométriques forment un domaine encore jeune et très actif qui s'est développé avec la CAO et l'infographie. La plupart des algorithmes requièrent un temps de l'ordre de $n$, $n log n$ ou $n^2$, ce qui n'est pas trop mal.

Ce domaine se situe entre l'infographie qui s'intéresse aux aspects plus visuels de la chose et la théorie des graphes qui fait abstraction de la position géométrique des objets et se concentre sur les liens qui les unissent.